diff --git a/appl/coupling-ff-pm/CMakeLists.txt b/appl/coupling-ff-pm/CMakeLists.txt index 46693285f4b230c23d125b2c5c8591d4d1627f1d..fe2d97cd2d64bd439644eaa5b5616e72350f390c 100644 --- a/appl/coupling-ff-pm/CMakeLists.txt +++ b/appl/coupling-ff-pm/CMakeLists.txt @@ -2,4 +2,4 @@ add_custom_target(test_exercises) add_subdirectory(monolithic) -add_subdirectory(solution) +add_subdirectory(iterative) diff --git a/appl/coupling-ff-pm/iterative/1pspatialparams.hh b/appl/coupling-ff-pm/iterative/1pspatialparams.hh new file mode 100644 index 0000000000000000000000000000000000000000..a042388f1eac51bb38ef980cb0f183756aaf6362 --- /dev/null +++ b/appl/coupling-ff-pm/iterative/1pspatialparams.hh @@ -0,0 +1,91 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see <http://www.gnu.org/licenses/>. * + *****************************************************************************/ +/*! + * \file + * \ingroup OnePTests + * \brief The spatial parameters class for the test problem using the 1p cc model + */ +#ifndef DUMUX_1P_TEST_SPATIALPARAMS_HH +#define DUMUX_1P_TEST_SPATIALPARAMS_HH + +#include <dumux/material/spatialparams/fv1p.hh> + +namespace Dumux { + +/*! + * \ingroup OnePModel + * + * \brief The spatial parameters class for the test problem using the + * 1p cc model + */ +template<class FVGridGeometry, class Scalar> +class OnePSpatialParams +: public FVSpatialParamsOneP<FVGridGeometry, Scalar, OnePSpatialParams<FVGridGeometry, Scalar>> +{ + using GridView = typename FVGridGeometry::GridView; + using ParentType = FVSpatialParamsOneP<FVGridGeometry, Scalar, OnePSpatialParams<FVGridGeometry, Scalar>>; + + using Element = typename GridView::template Codim<0>::Entity; + using GlobalPosition = typename Element::Geometry::GlobalCoordinate; + +public: + // export permeability type + using PermeabilityType = Scalar; + + OnePSpatialParams(std::shared_ptr<const FVGridGeometry> fvGridGeometry) + : ParentType(fvGridGeometry) + { + permeability_ = getParam<Scalar>("Darcy.SpatialParams.Permeability"); + porosity_ = getParam<Scalar>("Darcy.SpatialParams.Porosity"); + alphaBJ_ = getParam<Scalar>("Darcy.SpatialParams.AlphaBeaversJoseph"); + } + + /*! + * \brief Function for defining the (intrinsic) permeability \f$[m^2]\f$. + * + * \param globalPos The global position + * \return the intrinsic permeability + */ + PermeabilityType permeabilityAtPos(const GlobalPosition& globalPos) const + { return permeability_; } + + /*! \brief Define the porosity in [-]. + * + * \param globalPos The global position + */ + Scalar porosityAtPos(const GlobalPosition& globalPos) const + { return porosity_; } + + /*! \brief Define the Beavers-Joseph coefficient in [-]. + * + * \param globalPos The global position + */ + Scalar beaversJosephCoeffAtPos(const GlobalPosition& globalPos) const + { return alphaBJ_; } + + +private: + Scalar permeability_; + Scalar porosity_; + Scalar alphaBJ_; +}; + +} // end namespace + +#endif diff --git a/appl/coupling-ff-pm/iterative/CMakeLists.txt b/appl/coupling-ff-pm/iterative/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..ff70d2b3274017964ebc540a07e253f41d66dffe --- /dev/null +++ b/appl/coupling-ff-pm/iterative/CMakeLists.txt @@ -0,0 +1,11 @@ +# executables for the iterative case +dune_add_test(NAME test_ff + SOURCES main_ff.cc + COMPILE_DEFINITIONS ENABLEMONOLITHIC=0) + +dune_add_test(NAME test_pm + SOURCES main_pm.cc + COMPILE_DEFINITIONS ENABLEMONOLITHIC=0) + +# add a symlink for each input file +add_input_file_links() diff --git a/appl/coupling-ff-pm/iterative/main_ff.cc b/appl/coupling-ff-pm/iterative/main_ff.cc new file mode 100644 index 0000000000000000000000000000000000000000..f15f6295fdd6a167822d60311178fc17c15b7ec4 --- /dev/null +++ b/appl/coupling-ff-pm/iterative/main_ff.cc @@ -0,0 +1,158 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see <http://www.gnu.org/licenses/>. * + *****************************************************************************/ +/*! + * \file + * + * \brief A test problem for the coupled Stokes/Darcy problem (1p) + */ + #include <config.h> + + #include <ctime> + #include <iostream> + + #include <dune/common/parallel/mpihelper.hh> + #include <dune/common/timer.hh> + #include <dune/istl/io.hh> + + #include <dumux/common/properties.hh> + #include <dumux/common/parameters.hh> + #include <dumux/common/partial.hh> + #include <dumux/common/dumuxmessage.hh> + #include <dumux/linear/seqsolverbackend.hh> + #include <dumux/assembly/fvassembler.hh> + #include <dumux/assembly/diffmethod.hh> + #include <dumux/discretization/method.hh> + #include <dumux/io/vtkoutputmodule.hh> + #include <dumux/io/staggeredvtkoutputmodule.hh> + #include <dumux/io/grid/gridmanager.hh> + + #include <dumux/assembly/staggeredfvassembler.hh> + #include <dumux/nonlinear/newtonsolver.hh> + +#include "../monolithic/ffproblem.hh" + +int main(int argc, char** argv) try +{ + using namespace Dumux; + + // initialize MPI, finalize is done automatically on exit + const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv); + + // print dumux start message + if (mpiHelper.rank() == 0) + DumuxMessage::print(/*firstCall=*/true); + + // parse command line arguments and input file + Parameters::init(argc, argv); + + // Define the sub problem type tags + using FreeFlowTypeTag = Properties::TTag::FreeFlowModel; + + // try to create a grid (from the given grid file or the input file) + using FreeFlowGridManager = Dumux::GridManager<GetPropType<FreeFlowTypeTag, Properties::Grid>>; + FreeFlowGridManager freeFlowGridManager; + freeFlowGridManager.init("FreeFlow"); // pass parameter group + + // we compute on the leaf grid view + const auto& freeFlowGridView = freeFlowGridManager.grid().leafGridView(); + + // create the finite volume grid geometry + using FreeFlowFVGridGeometry = GetPropType<FreeFlowTypeTag, Properties::FVGridGeometry>; + auto freeFlowFvGridGeometry = std::make_shared<FreeFlowFVGridGeometry>(freeFlowGridView); + freeFlowFvGridGeometry->update(); + + // the problem (initial and boundary conditions) + using FreeFlowProblem = GetPropType<FreeFlowTypeTag, Properties::Problem>; + auto freeFlowProblem = std::make_shared<FreeFlowProblem>(freeFlowFvGridGeometry); + + // the solution vector + GetPropType<FreeFlowTypeTag, Properties::SolutionVector> sol; + sol[FreeFlowFVGridGeometry::cellCenterIdx()].resize(freeFlowFvGridGeometry->numCellCenterDofs()); + sol[FreeFlowFVGridGeometry::faceIdx()].resize(freeFlowFvGridGeometry->numFaceDofs()); + + + // apply initial solution for instationary problems + freeFlowProblem->applyInitialSolution(sol); + + // the grid variables + using FreeFlowGridVariables = GetPropType<FreeFlowTypeTag, Properties::GridVariables>; + auto freeFlowGridVariables = std::make_shared<FreeFlowGridVariables>(freeFlowProblem, freeFlowFvGridGeometry); + freeFlowGridVariables->init(sol); + + // intialize the vtk output module + StaggeredVtkOutputModule<FreeFlowGridVariables, decltype(sol)> freeFlowVtkWriter(*freeFlowGridVariables, sol, freeFlowProblem->name()); + GetPropType<FreeFlowTypeTag, Properties::IOFields>::initOutputModule(freeFlowVtkWriter); + freeFlowVtkWriter.addField(freeFlowProblem->getAnalyticalVelocityX(), "analyticalV_x"); + freeFlowVtkWriter.write(0.0); + + // the assembler for a stationary problem + using Assembler = StaggeredFVAssembler<FreeFlowTypeTag, DiffMethod::numeric>; + auto assembler = std::make_shared<Assembler>(freeFlowProblem, freeFlowFvGridGeometry, freeFlowGridVariables); + + // the linear solver + using LinearSolver = UMFPackBackend; + auto linearSolver = std::make_shared<LinearSolver>(); + + // the non-linear solver + using NewtonSolver = NewtonSolver<Assembler, LinearSolver>; + NewtonSolver nonLinearSolver(assembler, linearSolver); + + // solve the non-linear system + nonLinearSolver.solve(sol); + + // write vtk output + freeFlowVtkWriter.write(1.0); + + //////////////////////////////////////////////////////////// + // finalize, print dumux message to say goodbye + //////////////////////////////////////////////////////////// + + // print dumux end message + if (mpiHelper.rank() == 0) + { + Parameters::print(); + DumuxMessage::print(/*firstCall=*/false); + } + + return 0; +} // end main +catch (Dumux::ParameterException &e) +{ + std::cerr << std::endl << e << " ---> Abort!" << std::endl; + return 1; +} +catch (Dune::DGFException & e) +{ + std::cerr << "DGF exception thrown (" << e << + "). Most likely, the DGF file name is wrong " + "or the DGF file is corrupted, " + "e.g. missing hash at end of file or wrong number (dimensions) of entries." + << " ---> Abort!" << std::endl; + return 2; +} +catch (Dune::Exception &e) +{ + std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl; + return 3; +} +catch (...) +{ + std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl; + return 4; +} diff --git a/appl/coupling-ff-pm/iterative/main_pm.cc b/appl/coupling-ff-pm/iterative/main_pm.cc new file mode 100644 index 0000000000000000000000000000000000000000..fbdbb2c20df0fcc43ce873687f60a1311d35bf5d --- /dev/null +++ b/appl/coupling-ff-pm/iterative/main_pm.cc @@ -0,0 +1,192 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 2 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see <http://www.gnu.org/licenses/>. * + *****************************************************************************/ +/*! + * \file + * + * \brief A test problem for the coupled Stokes/Darcy problem (1p) + */ + #include <config.h> + + #include <ctime> + #include <iostream> + + #include <dune/common/parallel/mpihelper.hh> + #include <dune/common/timer.hh> + #include <dune/grid/io/file/dgfparser/dgfexception.hh> + #include <dune/grid/io/file/vtk.hh> + #include <dune/istl/io.hh> + + #include <dumux/common/properties.hh> + #include <dumux/common/parameters.hh> + #include <dumux/common/valgrind.hh> + #include <dumux/common/dumuxmessage.hh> + #include <dumux/common/defaultusagemessage.hh> + + #include <dumux/linear/amgbackend.hh> + #include <dumux/nonlinear/newtonsolver.hh> + + #include <dumux/assembly/fvassembler.hh> + #include <dumux/assembly/diffmethod.hh> + + #include <dumux/discretization/method.hh> + + #include <dumux/io/vtkoutputmodule.hh> + #include <dumux/io/grid/gridmanager.hh> + + #include "../monolithic/pmproblem.hh" + + /*! + * \brief Returns the pressure at the interface using Darcy's law for reconstruction + */ + template<class Problem, class Element, class SubControlVolumeFace, class VolumeVariables, class Scalar> + Scalar pressureAtInterface(const Problem& problem, + const Element& element, + const SubControlVolumeFace& scvf, + const VolumeVariables& volVars, + const Scalar boundaryFlux) + { + const Scalar cellCenterPressure = volVars.pressure(); + + // v = -K/mu * (gradP + rho*g) + auto velocity = scvf.unitOuterNormal(); + velocity *= boundaryFlux; // TODO check sign + + const Scalar ccPressure = volVars.pressure(); + const Scalar mobility = volVars.mobility(); + const Scalar density = volVars.density(); + const auto K = volVars.permeability(); + + // v = -kr/mu*K * (gradP + rho*g) = -mobility*K * (gradP + rho*g) + + const auto alpha = vtmv(scvf.unitOuterNormal(), K, problem().spatialParams.gravity(scvf.center())); + + auto distanceVector = scvf.center() - element.geometry().center(); + distanceVector /= distanceVector.two_norm2(); + const Scalar ti = vtmv(distanceVector, K, scvf.unitOuterNormal()); + + return (1/mobility * (scvf.unitOuterNormal() * velocity) + density * alpha)/ti + + ccPressure; + } + +int main(int argc, char** argv) try +{ + using namespace Dumux; + + // initialize MPI, finalize is done automatically on exit + const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv); + + // print dumux start message + if (mpiHelper.rank() == 0) + DumuxMessage::print(/*firstCall=*/true); + + // parse command line arguments and input file + Parameters::init(argc, argv); + + using DarcyTypeTag = Properties::TTag::DarcyOneP; + + using DarcyGridManager = Dumux::GridManager<GetPropType<DarcyTypeTag, Properties::Grid>>; + DarcyGridManager darcyGridManager; + darcyGridManager.init("Darcy"); // pass parameter group + + // we compute on the leaf grid view + const auto& darcyGridView = darcyGridManager.grid().leafGridView(); + + // create the finite volume grid geometry + using DarcyFVGridGeometry = GetPropType<DarcyTypeTag, Properties::FVGridGeometry>; + auto darcyFvGridGeometry = std::make_shared<DarcyFVGridGeometry>(darcyGridView); + darcyFvGridGeometry->update(); + + using DarcyProblem = GetPropType<DarcyTypeTag, Properties::Problem>; + auto darcyProblem = std::make_shared<DarcyProblem>(darcyFvGridGeometry); + + // the solution vector + GetPropType<DarcyTypeTag, Properties::SolutionVector> sol; + sol.resize(darcyFvGridGeometry->numDofs()); + + darcyProblem->applyInitialSolution(sol); + + // the grid variables + using DarcyGridVariables = GetPropType<DarcyTypeTag, Properties::GridVariables>; + auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyFvGridGeometry); + darcyGridVariables->init(sol); + + // intialize the vtk output module + const auto darcyName = getParam<std::string>("Problem.Name") + "_" + darcyProblem->name(); + + VtkOutputModule<DarcyGridVariables, GetPropType<DarcyTypeTag, Properties::SolutionVector>> darcyVtkWriter(*darcyGridVariables, sol, darcyName); + using DarcyVelocityOutput = GetPropType<DarcyTypeTag, Properties::VelocityOutput>; + darcyVtkWriter.addVelocityOutput(std::make_shared<DarcyVelocityOutput>(*darcyGridVariables)); + GetPropType<DarcyTypeTag, Properties::IOFields>::initOutputModule(darcyVtkWriter); + darcyVtkWriter.write(0.0); + + // the assembler for a stationary problem + using Assembler = FVAssembler<DarcyTypeTag, DiffMethod::numeric>; + auto assembler = std::make_shared<Assembler>(darcyProblem, darcyFvGridGeometry, darcyGridVariables); + + // the linear solver + using LinearSolver = UMFPackBackend; + auto linearSolver = std::make_shared<LinearSolver>(); + + // the non-linear solver + using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>; + NewtonSolver nonLinearSolver(assembler, linearSolver); + + // solve the non-linear system + nonLinearSolver.solve(sol); + + // write vtk output + darcyVtkWriter.write(1.0); + + //////////////////////////////////////////////////////////// + // finalize, print dumux message to say goodbye + //////////////////////////////////////////////////////////// + + // print dumux end message + if (mpiHelper.rank() == 0) + { + Parameters::print(); + DumuxMessage::print(/*firstCall=*/false); + } + + return 0; +} // end main +catch (Dumux::ParameterException &e) +{ + std::cerr << std::endl << e << " ---> Abort!" << std::endl; + return 1; +} +catch (Dune::DGFException & e) +{ + std::cerr << "DGF exception thrown (" << e << + "). Most likely, the DGF file name is wrong " + "or the DGF file is corrupted, " + "e.g. missing hash at end of file or wrong number (dimensions) of entries." + << " ---> Abort!" << std::endl; + return 2; +} +catch (Dune::Exception &e) +{ + std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl; + return 3; +} +catch (...) +{ + std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl; + return 4; +} diff --git a/appl/coupling-ff-pm/iterative/params.input b/appl/coupling-ff-pm/iterative/params.input new file mode 100644 index 0000000000000000000000000000000000000000..4117b7f17c501817b738dfdb6836cc2c3223e4a8 --- /dev/null +++ b/appl/coupling-ff-pm/iterative/params.input @@ -0,0 +1,46 @@ + # for dune-subgrid + [Grid] +Positions0 = 0 1 +Positions1 = 0 0.2 0.3 0.65 +Cells0 = 100 +Cells1 = 10 50 18 +Baseline = 0.25 # [m] +Amplitude = 0.04 # [m] +Offset = 0.5 # [m] +Scaling = 0.2 #[m] + +[Stokes.Grid] +Verbosity = true +Positions0 = 0.0 1.0 +Positions1 = 1.0 2.0 +Cells0 = 20 +Cells1 = 100 +Grading1 = 1 + +[Darcy.Grid] +Verbosity = true +Positions0 = 0.0 1.0 +Positions1 = 0.0 1.0 +Cells0 = 20 +Cells1 = 20 +Grading1 = 1 + +[Stokes.Problem] +Name = stokes +PressureDifference = 1e-9 +EnableInertiaTerms = false + +[Darcy.Problem] +Name = darcy + +[Darcy.SpatialParams] +Permeability = 1e-6 # m^2 +Porosity = 0.4 +AlphaBeaversJoseph = 1.0 + +[Problem] +Name = ex_ff-pm-interface +EnableGravity = false + +[Vtk] +AddVelocity = 1 diff --git a/appl/coupling-ff-pm/monolithic/CMakeLists.txt b/appl/coupling-ff-pm/monolithic/CMakeLists.txt index bae4c263a6bf7ae6a06d90af66a2e45d224f9b2e..4d6a867cb8db5d2ec8339c0bb18a7ccbc980a70e 100644 --- a/appl/coupling-ff-pm/monolithic/CMakeLists.txt +++ b/appl/coupling-ff-pm/monolithic/CMakeLists.txt @@ -1,3 +1,6 @@ -add_subdirectory(interface) -add_subdirectory(models) -add_subdirectory(turbulence) +# executables for the monolithic case +dune_add_test(NAME ff-pm_monolithic + SOURCES main.cc) + +# add a symlink for each input file +add_input_file_links() diff --git a/appl/coupling-ff-pm/monolithic/ffproblem.hh b/appl/coupling-ff-pm/monolithic/ffproblem.hh index cbe437eeb9ce1d90db2328932e8866d151560f61..35a2a43efdba19b928a5e88cd4985eb30acaa75d 100644 --- a/appl/coupling-ff-pm/monolithic/ffproblem.hh +++ b/appl/coupling-ff-pm/monolithic/ffproblem.hh @@ -23,10 +23,11 @@ #ifndef DUMUX_STOKES_SUBPROBLEM_HH #define DUMUX_STOKES_SUBPROBLEM_HH -#include <dune/grid/yaspgrid.hh> +#ifndef ENABLEMONOLITHIC +#define ENABLEMONOLITHIC 1 +#endif -//****** uncomment for the last exercise *****// -// #include <dumux/io/grid/subgridgridcreator.hh> +#include <dune/grid/yaspgrid.hh> #include <dumux/material/fluidsystems/1pliquid.hh> #include <dumux/material/components/simpleh2o.hh> @@ -44,12 +45,12 @@ namespace Properties { // Create new type tags namespace TTag { -struct StokesOneP { using InheritsFrom = std::tuple<NavierStokes, StaggeredFreeFlowModel>; }; +struct FreeFlowModel { using InheritsFrom = std::tuple<NavierStokes, StaggeredFreeFlowModel>; }; } // end namespace TTag // the fluid system template<class TypeTag> -struct FluidSystem<TypeTag, TTag::StokesOneP> +struct FluidSystem<TypeTag, TTag::FreeFlowModel> { using Scalar = GetPropType<TypeTag, Properties::Scalar>; using type = FluidSystems::OnePLiquid<Scalar, Dumux::Components::SimpleH2O<Scalar> > ; @@ -57,7 +58,7 @@ struct FluidSystem<TypeTag, TTag::StokesOneP> // Set the grid type template<class TypeTag> -struct Grid<TypeTag, TTag::StokesOneP> +struct Grid<TypeTag, TTag::FreeFlowModel> { static constexpr auto dim = 2; using Scalar = GetPropType<TypeTag, Properties::Scalar>; @@ -73,14 +74,14 @@ struct Grid<TypeTag, TTag::StokesOneP> // Set the problem property template<class TypeTag> -struct Problem<TypeTag, TTag::StokesOneP> { using type = Dumux::StokesSubProblem<TypeTag> ; }; +struct Problem<TypeTag, TTag::FreeFlowModel> { using type = Dumux::StokesSubProblem<TypeTag> ; }; template<class TypeTag> -struct EnableFVGridGeometryCache<TypeTag, TTag::StokesOneP> { static constexpr bool value = true; }; +struct EnableFVGridGeometryCache<TypeTag, TTag::FreeFlowModel> { static constexpr bool value = true; }; template<class TypeTag> -struct EnableGridFluxVariablesCache<TypeTag, TTag::StokesOneP> { static constexpr bool value = true; }; +struct EnableGridFluxVariablesCache<TypeTag, TTag::FreeFlowModel> { static constexpr bool value = true; }; template<class TypeTag> -struct EnableGridVolumeVariablesCache<TypeTag, TTag::StokesOneP> { static constexpr bool value = true; }; +struct EnableGridVolumeVariablesCache<TypeTag, TTag::FreeFlowModel> { static constexpr bool value = true; }; } /*! @@ -109,11 +110,18 @@ class StokesSubProblem : public NavierStokesProblem<TypeTag> using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>; using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>; +#if ENABLEMONOLITHIC using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>; +#endif public: +#if ENABLEMONOLITHIC StokesSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry, std::shared_ptr<CouplingManager> couplingManager) : ParentType(fvGridGeometry, "Stokes"), eps_(1e-6), couplingManager_(couplingManager) +#else + StokesSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry) + : ParentType(fvGridGeometry, "Stokes"), eps_(1e-6) +#endif { deltaP_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.PressureDifference"); } @@ -162,8 +170,6 @@ public: // left/right wall if (onRightBoundary_(globalPos) || (onLeftBoundary_(globalPos))) { -// values.setDirichlet(Indices::velocityXIdx); -// values.setDirichlet(Indices::velocityYIdx); values.setDirichlet(Indices::pressureIdx); } else @@ -173,13 +179,22 @@ public: } // coupling interface +#if ENABLEMONOLITHIC if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf)) { values.setCouplingNeumann(Indices::conti0EqIdx); values.setCouplingNeumann(Indices::momentumYBalanceIdx); -// values.setDirichlet(Indices::velocityXIdx); // assume no slip on interface - values.setBJS(Indices::momentumXBalanceIdx); // assume no slip on interface + values.setBJS(Indices::momentumXBalanceIdx); } +#else + // // TODO do preCICE stuff in analogy to heat transfer + // if(/*preCice*/) + // { + // values.setCouplingNeumann(Indices::conti0EqIdx); + // values.setCouplingNeumann(Indices::momentumYBalanceIdx); + // values.setBJS(Indices::momentumXBalanceIdx); + // } +#endif return values; } @@ -193,18 +208,6 @@ public: { PrimaryVariables values(0.0); values = initialAtPos(globalPos); - -// if(onUpperBoundary_(globalPos)) -// { -// values[Indices::velocityXIdx] = 0.0; -// values[Indices::velocityYIdx] = 0.0; -// } - -// if(onLeftBoundary_(globalPos)) -// values[Indices::pressureIdx] = deltaP_; -// if(onRightBoundary_(globalPos)) -// values[Indices::pressureIdx] = 0.0; - return values; } @@ -226,24 +229,31 @@ public: { NumEqVector values(0.0); +#if ENABLEMONOLITHIC if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf)) { values[Indices::conti0EqIdx] = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf); values[Indices::momentumYBalanceIdx] = couplingManager().couplingData().momentumCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf); } +#else + // if(/*preCICE*/) + // { + // const Scalar density = 1000; // TODO how to handle compressible fluids? + // values[Indices::conti0EqIdx] = elemFaceVars[scvf].velocitySelf() * scvf.directionSign() * density; + // values[Indices::momentumYBalanceIdx] = /*pressure from Darcy*/ + // } +#endif return values; } // \} - //! Set the coupling manager - void setCouplingManager(std::shared_ptr<CouplingManager> cm) - { couplingManager_ = cm; } - +#if ENABLEMONOLITHIC //! Get the coupling manager const CouplingManager& couplingManager() const { return *couplingManager_; } +#endif /*! * \name Volume terms @@ -272,7 +282,11 @@ public: */ Scalar permeability(const Element& element, const SubControlVolumeFace& scvf) const { +#if ENABLEMONOLITHIC return couplingManager().couplingData().darcyPermeability(element, scvf); +#else + return 1e-10; // TODO transfer information or just use constant value +#endif } /*! @@ -280,7 +294,11 @@ public: */ Scalar alphaBJ(const SubControlVolumeFace& scvf) const { +#if ENABLEMONOLITHIC return couplingManager().problem(CouplingManager::darcyIdx).spatialParams().beaversJosephCoeffAtPos(scvf.center()); +#else + return 1.0; // TODO transfer information or just use constant value +#endif } /*! @@ -338,7 +356,9 @@ private: Scalar eps_; Scalar deltaP_; +#if ENABLEMONOLITHIC std::shared_ptr<CouplingManager> couplingManager_; +#endif mutable std::vector<Scalar> analyticalVelocityX_; }; diff --git a/appl/coupling-ff-pm/monolithic/main.cc b/appl/coupling-ff-pm/monolithic/main.cc index 8c44bfbf7b259e0d99cdc3efb157ecf381eb1f81..6d228a9f0a940e01860d02904913ed0fd3a4149e 100644 --- a/appl/coupling-ff-pm/monolithic/main.cc +++ b/appl/coupling-ff-pm/monolithic/main.cc @@ -46,14 +46,14 @@ #include <dumux/multidomain/boundary/stokesdarcy/couplingmanager.hh> -#include "ex_interface_pmproblem.hh" -#include "ex_interface_ffproblem.hh" +#include "pmproblem.hh" +#include "ffproblem.hh" namespace Dumux { namespace Properties { template<class TypeTag> -struct CouplingManager<TypeTag, TTag::StokesOneP> +struct CouplingManager<TypeTag, TTag::FreeFlowModel> { using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, Properties::TTag::DarcyOneP>; using type = Dumux::StokesDarcyCouplingManager<Traits>; @@ -62,7 +62,7 @@ struct CouplingManager<TypeTag, TTag::StokesOneP> template<class TypeTag> struct CouplingManager<TypeTag, TTag::DarcyOneP> { - using Traits = StaggeredMultiDomainTraits<Properties::TTag::StokesOneP, Properties::TTag::StokesOneP, TypeTag>; + using Traits = StaggeredMultiDomainTraits<Properties::TTag::FreeFlowModel, Properties::TTag::FreeFlowModel, TypeTag>; using type = Dumux::StokesDarcyCouplingManager<Traits>; }; @@ -84,7 +84,7 @@ int main(int argc, char** argv) try Parameters::init(argc, argv); // Define the sub problem type tags - using StokesTypeTag = Properties::TTag::StokesOneP; + using StokesTypeTag = Properties::TTag::FreeFlowModel; using DarcyTypeTag = Properties::TTag::DarcyOneP; diff --git a/appl/coupling-ff-pm/monolithic/pmproblem.hh b/appl/coupling-ff-pm/monolithic/pmproblem.hh index 5b56c5aad1d6e4d0cf62b50239470be82ad5ce3b..49565ca35aba6ef9f177430a8b8400849f9e96fb 100644 --- a/appl/coupling-ff-pm/monolithic/pmproblem.hh +++ b/appl/coupling-ff-pm/monolithic/pmproblem.hh @@ -24,6 +24,10 @@ #ifndef DUMUX_DARCY_SUBPROBLEM_HH #define DUMUX_DARCY_SUBPROBLEM_HH +#ifndef ENABLEMONOLITHIC +#define ENABLEMONOLITHIC 1 +#endif + #include <dune/grid/yaspgrid.hh> //****** uncomment for the last exercise *****// @@ -34,7 +38,7 @@ #include <dumux/porousmediumflow/1p/model.hh> #include <dumux/porousmediumflow/problem.hh> -#include "../1pspatialparams.hh" +#include "1pspatialparams.hh" #include <dumux/material/components/simpleh2o.hh> #include <dumux/material/fluidsystems/1pliquid.hh> @@ -109,12 +113,19 @@ class DarcySubProblem : public PorousMediumFlowProblem<TypeTag> using Element = typename GridView::template Codim<0>::Entity; using GlobalPosition = typename Element::Geometry::GlobalCoordinate; +#if ENABLEMONOLITHIC using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>; +#endif public: +#if ENABLEMONOLITHIC DarcySubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry, std::shared_ptr<CouplingManager> couplingManager) : ParentType(fvGridGeometry, "Darcy"), eps_(1e-7), couplingManager_(couplingManager) +#else +DarcySubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry) + : ParentType(fvGridGeometry, "Darcy"), eps_(1e-7) +#endif {} /*! @@ -149,14 +160,11 @@ public: // set Neumann BCs to all boundaries first values.setAllNeumann(); +#if ENABLEMONOLITHIC // set the coupling boundary condition at the interface if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf)) values.setAllCouplingNeumann(); - - // set a Dirichlet boundary condition at the bottom - //if (onLowerBoundary_(scvf.center())) - // values.setAllDirichlet(); - +#endif return values; } @@ -196,10 +204,16 @@ public: // no-flow everywhere ... NumEqVector values(0.0); +#if ENABLEMONOLITHIC // ... except at the coupling interface if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf)) values[Indices::conti0EqIdx] = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, scvf); - +#else + // if (/*preCICE*/) + // { + // values[Indices::conti0EqIdx] = /*mass flux from stokes*/; + // } +#endif return values; } @@ -242,13 +256,11 @@ public: // \} - //! Set the coupling manager - void setCouplingManager(std::shared_ptr<CouplingManager> cm) - { couplingManager_ = cm; } - +#if ENABLEMONOLITHIC //! Get the coupling manager const CouplingManager& couplingManager() const { return *couplingManager_; } +#endif private: bool onLeftBoundary_(const GlobalPosition &globalPos) const @@ -264,7 +276,10 @@ private: { return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_; } Scalar eps_; + +#if ENABLEMONOLITHIC std::shared_ptr<CouplingManager> couplingManager_; +#endif }; } //end namespace