Skip to content
GitLab
Commit 5df03184 authored by Thomas Fetzer's avatar Thomas Fetzer Committed by Kilian Weishaupt
Browse files

[coupling-ff-pm] Add model exercise

parent 3bb596eb
Hide whitespace changes
Inline Side-by-side
exercises/exercise-coupling-ff-pm/interface/1pspatialparams.hh → exercises/exercise-coupling-ff-pm/1pspatialparams.hh
View file @ 5df03184
......@@ -59,6 +59,7 @@ public:
: ParentType(fvGridGeometry)
{
permeability_ = getParam<Scalar>("Darcy.SpatialParams.Permeability");
porosity_ = getParam<Scalar>("Darcy.SpatialParams.Porosity");
alphaBJ_ = getParam<Scalar>("Darcy.SpatialParams.AlphaBeaversJoseph");
}
......@@ -76,7 +77,7 @@ public:
* \param globalPos The global position
*/
Scalar porosityAtPos(const GlobalPosition& globalPos) const
{ return 0.4; }
{ return porosity_; }
/*! \brief Define the Beavers-Joseph coefficient in [-].
*
......@@ -88,6 +89,7 @@ public:
private:
Scalar permeability_;
Scalar porosity_;
Scalar alphaBJ_;
};
......
exercises/exercise-coupling-ff-pm/2pspatialparams.hh 0 → 100644
View file @ 5df03184
// -*- 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 TwoPTests
* \brief The spatial parameters class for the test problem using the 2p cc model
*/
#ifndef DUMUX_TWOPHASE_SPATIAL_PARAMS_HH
#define DUMUX_TWOPHASE_SPATIAL_PARAMS_HH
#include <dumux/material/spatialparams/fv.hh>
#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
#include <dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh>
#include <dumux/material/fluidmatrixinteractions/2p/thermalconductivitysomerton.hh>
namespace Dumux
{
/*!
* \ingroup TwoPModel
* \ingroup ImplicitTestProblems
*
* \brief The spatial parameters class for the test problem using the 2p cc model
*/
template<class TypeTag>
class TwoPSpatialParams
: public FVSpatialParams<typename GET_PROP_TYPE(TypeTag, FVGridGeometry),
typename GET_PROP_TYPE(TypeTag, Scalar),
TwoPSpatialParams<TypeTag>>
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Element = typename GridView::template Codim<0>::Entity;
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using FVElementGeometry = typename FVGridGeometry::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using ParentType = FVSpatialParams<FVGridGeometry, Scalar, TwoPSpatialParams<TypeTag>>;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
using EffectiveLaw = RegularizedVanGenuchten<Scalar>;
public:
using MaterialLaw = EffToAbsLaw<EffectiveLaw>;
using MaterialLawParams = typename MaterialLaw::Params;
using PermeabilityType = Scalar;
TwoPSpatialParams(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");
// residual saturations
params_.setSwr(getParam<Scalar>("Darcy.SpatialParams.Swr"));
params_.setSnr(getParam<Scalar>("Darcy.SpatialParams.Snr"));
// parameters for the vanGenuchten law
params_.setVgAlpha(getParam<Scalar>("Darcy.SpatialParams.VgAlpha"));
params_.setVgn(getParam<Scalar>("Darcy.SpatialParams.VgN"));
params_.setPcLowSw(params_.swr()*5.0);
params_.setPcHighSw(1.0-params_.snr()*5.0);
}
/*!
* \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_; }
/*!
* \brief Returns the parameter object for the Brooks-Corey material law.
* In this test, we use element-wise distributed material parameters.
*
* \param element The current element
* \param scv The sub-control volume inside the element.
* \param elemSol The solution at the dofs connected to the element.
* \return the material parameters object
*/
template<class ElementSolutionVector>
const MaterialLawParams& materialLawParams(const Element& element,
const SubControlVolume& scv,
const ElementSolutionVector& elemSol) const
{ return params_; }
/*!
* \brief Function for defining which phase is to be considered as the wetting phase.
*
* \return the wetting phase index
* \param globalPos The global position
*/
template<class FluidSystem>
int wettingPhaseAtPos(const GlobalPosition& globalPos) const
{ return FluidSystem::phase0Idx; }
private:
Scalar permeability_;
Scalar porosity_;
Scalar alphaBJ_;
MaterialLawParams params_;
static constexpr Scalar eps_ = 1.0e-7;
};
} // end namespace
#endif
exercises/exercise-coupling-ff-pm/CMakeLists.txt
View file @ 5df03184
add_subdirectory(interface)
add_subdirectory(models)
exercises/exercise-coupling-ff-pm/README.md
View file @ 5df03184
# Exercise #5 (DuMuX course)
The aim of this exercise is to get familiar with the set up of coupled free flow/porous medium flow problems.
The aim of this exercise is to get familiar with setting up coupled free flow/porous medium flow problems.
## Problem set-up
The model domain consists of two non-overlapping subdomains.
Free flow is modeled in the upper subdomain, while the lower subdomain models a flow in a porous medium.
Both single-phase flow and two-phase flow will be considered in the porous domain.
![](../extradoc/exercise5_setup.png)
Free flow is modeled in the upper subdomain, while the lower subdomain models a flow in a porous medium. Both single-phase flow and two-phase flow will be considered in the porous domain.
### 0. Getting familiar with the code
......@@ -20,21 +15,24 @@ Both single-phase flow and two-phase flow will be considered in the porous domai
There are three sub folders: `interface`, `models` and `turbulence`.
The problem-related files for this exercise are:
* Three __main files__ for the three sub-tasks :`ex_interface_coupling_ff-pm.cc`, `ex2_coupling_ff-pm.cc`, `ex2_coupling_ff-pm.cc`,
* Three __free flow problem files__: `ex_interface_ffproblem.hh`, `ex2_ffproblem.hh`, `ex3_ffproblem.hh`
* Three __porous medium flow problem files__: `ex_interface_pmproblem.hh`, `ex2_pmproblem.hh`, `ex3_pmproblem.hh`
* The __input files__: ...
* Three __main files__ for the three sub-tasks :`ex_interface_coupling_ff-pm.cc`, `ex_models_coupling_ff-pm.cc`, `ex_turbulence_coupling_ff-pm.cc`,
* Three __free flow problem files__: `ex_interface_ffproblem.hh`, `ex_models_ffproblem.hh`, `ex_turbulence__ffproblem.hh`
* Three __porous medium flow problem files__: `ex_interface_pmproblem.hh`, `ex_models_pmproblem.hh`, `ex_turbulence_pmproblem.hh`
* The __input files__: `ex_interface_coupling_ff-pm.input`, `ex_models_coupling_ff-pm.input`, `ex_turbulence_coupling_ff-pm.input`,
* The __spatial parameters files__: `1pspatialparams.hh`, `2pspatialparams.hh`
In the main file, `TypeTags` for both submodels are defined. The same applies for types such as `GridManager`, `FVGridGeometry`, `Problem`, etc.. Since we use a monolithic coupling scheme, there is only one `Assembler` and one
`NewtonSolver`.
In the main file, `TypeTags` for both submodels are defined.
The same applies for types such as `GridManager`, `FVGridGeometry`, `Problem`, etc...
Since we use a monolithic coupling scheme, there is only one `Assembler` and one `NewtonSolver`.
The problem files very much look like "regular", uncoupled ones with the exception that they hold a pointer to the `CouplingManager` which allows to evaluate the coupling conditions and to exchange information between the coupled models. The coupling conditions are realized technically in terms of boundary condition. For instance, in lines 178 and 179
The problem files very much look like "regular", uncoupled ones with the exception that they hold a pointer to the `CouplingManager` which allows to evaluate the coupling conditions and to exchange information between the coupled models.
The coupling conditions are realized technically in terms of boundary condition. For instance, in lines 171 and 172
in `ex_interface_ffproblem.hh`, `couplingNeumann` boundary conditions are set, which means that the free flow models evaluates the
mass and momentum fluxes coming from the porous domain and uses these values as boundary conditions at the interface.
Note the certain checks are performed when combining different models, e.g., the fluid system has to be the same for both domains.
Note the certain checks are performed when combining different models, e.g., the fluid system has to be the same for both domains
and the sub-control-volume faces at the interface have to match.
We will use a staggered grid for the free flow and a cell-centered finite volume method for the porous medium.
......@@ -42,35 +40,28 @@ Keep in mind that the staggered grid implementation distinguishes between face v
cell center variables (all other variables), therefore in some cases either the `stokesCellCenterIdx`
or the `stokesFaceIdx` is used respectively, while for the porous medium all variables can be accessed with `darcyIdx`.
The main part of the Stokes/Darcy coupling is implemented in the folder `dumux/multidomain/boundary/stokesdarcy`.
__Task__:
Take a closer look at the Stokes/Darcy coupling files and try to answer the following questions before moving to the next part of the exercise:
* Where are the elements at the coupling interface paired with their respective coupling partner elements of the other domain?
* Where are the coupling conditions implemented and evaluated?
* You will find four implementations of the `massCouplingCondition`. What are the differences between them?
* Which data is stored and updated regularly such that it is accessible from the respective other domain?
* How can we access volume variables like density and pressure from the respective other domain?
* ...
Take a closer look at the Stokes/Darcy coupling files before moving to the next part of the exercise:
### 1. Changing the interface
In this part of the exercise, a simple coupled system consisting of a one-phase (1p) free flow and a one-phase flow in a porous medium is set up. Both subproblems have no-flow boundaries at the sides.
A velocity profile is set on the upper free flow boundary, which leads to a vertical flow into the porous medium.
Currently, a velocity profile is set on the upper free flow boundary, which leads to a vertical flow into the porous medium:
![](../extradoc/ex_ff-pm-vertical-flow.png)
* We will first change the flow direction such that the free flow is parallel to the porous medium.
* Afterwards, the Beavers-Joseph-Saffman condition will be used as an interface condition for the tangential momentum transfer.
* Last, we change the flat interface between the two domains to a wave-shaped one.
__Tasks__:
__Task A: Change the flow direction__
Open the file `ex_interface_ffproblem.hh` and navigate to line 155, where the types of boundary condition are set.
Open the file `ex_interface_ffproblem.hh` and navigate to line 148, where the types of boundary condition are set.
Instead of applying a fixed velocity profile at the top of the domain, we want to use fixed pressure boundary conditions
at the left and right side of the free flow domain, while the top represents an impermeable wall.
__Change the flow direction__
Set a Dirichlet boundary condition for the pressure at the left and ride side of the domain:
Set a Dirichlet boundary condition for the pressure at the left and right side of the domain:
``` cpp
if(onLeftBoundary_(globalPos) || onRightBoundary_(globalPos))
values.setDirichlet(Indices::pressureIdx);
......@@ -120,18 +111,19 @@ if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
values.setAllCouplingNeumann();
```
__Include slip-condition__
This should make the flow go from left to right.
This should make the flow go from left to right. However, we are still missing one important feature:
the Beavers-Joseph-Saffman slip condition at the interface:
__Task B: Include slip-condition__
However, we are still missing one important feature:
at the moment, the velocity component tangential to the interface gets a no-slip condition.
In the next step we want to implement the Beavers-Joseph-Saffman slip condition at the interface:
$`\frac{\partial v_x}{\partial y} = \frac{\alpha}{\sqrt K} (v_x - q_{pm})\quad`$ at $`\quad y=0`$
with $`\quad q_{pm}=0`$.
The include this, we just set:
To include this, just set a Beavers-Joseph-Saffman (BJS) boundary condition for the respective momentum balance:
``` cpp
values.setBJS(Indices::momentumXBalanceIdx);
```
......@@ -139,15 +131,15 @@ values.setBJS(Indices::momentumXBalanceIdx);
at the position where the coupling boundary conditions are set in `ex_interface_ffproblem.hh`.
To check if the simulation behaves as expected, we can compare the velocity profile $`v_x(y)`$ with the analytical solution provided by [Beavers and Joseph (1967)](https://doi.org/10.1017/S0022112067001375).
For doing so, we uncomment line 212
For doing so, we uncomment line 212 in `ex_interface_coupling_ff-pm.cc`.
```cpp
stokesVtkWriter.addField(stokesProblem->getAnalyticalVelocityX(), "analyticalV_x");
```
in `ex_interface_coupling_ff-pm.cc`. After re-compiling and re-running the executable, we should be able to see also
After re-compiling and re-running the executable, we should be able to see also
the analytical solution of $`v_x`$ on the free flow domain. Play around with the grid resolution to see how that affects the velocity profile.
__Cange shape of interface__
__Task C: Cange shape of interface__
Now we want to include a non-flat interface between the two domains. We use `dune-subgrid` to construct
two grids for the two domains from one common host grid. Comment out lines 93-103 in `ex_interface_coupling_ff-pm.cc` and comment lines 112-147 in the same file. This will instantiate a host grid and define two helper lambda functions that are used to choose elements from to host grid for the respective sub grid. In the given case,
......@@ -167,14 +159,14 @@ auto elementSelectorDarcy = [&](const auto& element)
};
```
Make sure to comment in line 30 in both problem files
Make sure, that you have uncommented the lines including the gridcreators in both problem files
```cpp
#include <dumux/io/grid/subgridgridcreator.hh>
```
and do the changes in the respective lines for the `Grid` property.
The problem should no compile and run. However, an error occurs due to the coupling conditions.
The problem should now compile and run. However, an error occurs due to the coupling conditions.
So far, we assumed a flat interface, therefore the normal momentum coupling condition
$`[\sigma \cdot \mathbf{n}]^{FF} = p^{PM}`$
......@@ -193,15 +185,114 @@ values.setCouplingNeumann(1 - scvf.directionIndex());
```
The final result should look something like this:
![](../extradoc/ex_ff-pm-wave-interface.png)
### 3. Change the models
### 2. Changing the porous medium model
1p2c/1p2c -->1p2c/2p2c
In this part of the exercise, the coupled system will be extended such that the transport of components
in both parts is included and the presence of a second phase in the porous medium is considered.
This enables the simulation of the drying of a porous medium (however no energy balance is included yet).
__Task__:
We start with an example in which the transport of water vapor in the gas phase is considered.
The porous medium is filled with gas, initially the mole fraction of water vapor is $`x^w_g = 0.1`$.
Above the porous medium, a dry gas flows and by diffusion, the porous medium dries out.
__Task A: Change the model__:
In the first task, the porous-medium model will be changed from a 1p2c system to a 2p2c system.
Although a 2p2c system is plugged in, we still want to simulate the same situation as before.
The following changes have to be made in the porous-medium model (`ex_models_pmproblem.hh`):
* Include the 2pnc model: include the respective headers and inherit from the new model `TwoPNC`
* Exchange the spatial parameters for the 1-phase system by those for a 2-phase system (hint: two occurrences).
One big difference between the 1p and 2p model is, that the primary variables for which
the problem is solved, are not fixed.
It is possible to use different formulations, e.g. with the gas pressure and the
liquid saturation as primary variables (p_g-S_l -> `p1s0`) or vice versa.
* Set the property
```
SET_PROP(DarcyTypeTag, Formulation)
{ static constexpr auto value = TwoPFormulation::p1s0; };
```
in the Properties section in the problem file.
In contrast to the formulation, which stays the same during one simulation, the meaning of
the primary variables may vary during one simulation.
In the case we are investigating, we want to use the gas pressure and the liquid saturation as primary variables.
However, if only the gas phase is present, the liquid saturation is always zero.
In this case, the chosen formulation will set the given value as the mole fraction of water vapor in the gas phase.
* To tell to program which phases are present in which parts of the domain at the beginning of the simulation,
you have to call `values.setState(MY_PHASE_PRESENCE);` in `initialAtPos()`. Have a look at the `indices.hh`
in the `2p2c` model to figure out which is the correct value of `MY_PHASE_PRESENCE` for the presence of
a gas-phase only (hint: the numbering of phase indices begins with 0, the numbering of the phases begins
with 1. Take a look at your formulation to find out which phase index to use for the gas phase.)
* Since two phases are involved now, the index `fluidSystemPhaseIdx` does not make sense any more and needs to be removed
from the problem file. Do so, and set the `transportCompIdx` to `Indices::switchIdx`.
Compile and run the new simulation.
(Don't mind the compiler warning, we will deal with it in the next task.)
__Task B: Add output__:
In the next step, we want to add some output to the simulation.
First we want to know the water mass in the (porous-medium) system.
Therefore, we evaluate the storage term of the water component.
* Have a look at the function `evaluateWaterMassStorageTerm()` in the porous medium subproblem.
Then implement a calculation of the total water mass:
$`\sum_{\alpha \in \textrm{g,l}} \left( \phi S_\alpha X^\text{w}_\alpha \varrho_\alpha V_\textrm{scv} \right)`$.
Afterwards, call the method from `init()` and save the result in the respective private variable needed
by the `evaluateWaterMassStorageTerm` method. In addition, you need to set an initial value for
another variable used in the same method.
We also want to investigate the temporal evolution of the water mass.
The following steps need to be done to do so.
Check if all instructions are implemented accordingly:
* Calculate the initial water mass at the beginning of the simulation and add the water mass loss to the output.
Based on the water mass loss you can derive the evaporation rate. Because despite at the interface, all
boundaries have Neumann no-flow conditions and no sources are present, the water can only leave the system
via the porous-medium free-flow interface. The evaporation in [mm/d] can be calculated by:
$`e = \frac{\textrm{d}\, M^\text{w}}{\textrm{d}\, t} / A_\textrm{interface}`$.
Finally we want to know the distribution of the water mass fluxes across the interface.
* Have a look at the function `evaluateInterfaceFluxes()` in the porous medium problem.
Use the facilities therein to return the values of ...massCouplingCondition... from the `couplingManager`
for each coupling scvf. Then the fluxes are visualized with gnuplot, when setting `Problem.PlotFluxes = true`.
If the simulation is to fast, you can have a look at the flux*.png files after the simulation.
* You can use the property `Problem.PlotStorage = true` to see the temporal evolution of the evaporation rate
and the cumulative water mass loss.
Compile and run the simulation and take a look at the results.
__Task C: Let it dry out__:
In this exercise we want to completely dry an initial water-filled porous medium.
- Change the initial condition such that it is started with a two-phase system.
Set the initial saturation to $`S_\text{l} = 0.1`$.
If one wants to simulate the complete drying of a porous medium, the standard capillary pressure--saturation
relationships have to be regularized. If no regularization is used, then the capillary pressure would
approach infinity when the liquid saturation goes to zero. This means that water can flow to the
interface from any region of the porous medium and of course this poses numerical problems.
Therefore, the capillary pressure--saturation relationship has to be regularized for low saturations.
The regularization has a strong influence on how long liquid water is present at the interface, see
[Mosthaf (2014)](http://dx.doi.org/10.18419/opus-519).
* Take a look at how the regularization is set in the `2pspatialparams.hh` and see how adapating
the parameter for `pcLowSw` and `pcHighSw` affects the results.
The porous-medium model can now reach a liquid saturation of zero. As already explained above,
for this case the liquid saturation cannot serve as primary variable anymore. However,
manually adapting the primary variable states and values is inconvenient.
[Class et al. (2007)](http://dx.doi.org/10.1016/S0309-1708(02)00014-3)
describe an algorithm to switch the primary variables, if phases should appear or disappear during a simulation.
- Replace the current implementation of the Newton solver with the version which can handle
primary variable switches: `dumux/multidomain/privarswitchnewtonsolver.hh`.
You also have to uncomment the line containing the `PriVarSwitchTuple` and to overwrite
the last argument with the `PrimaryVariableSwitch` property from the Darcy model.
Now you are able to simulate a complete drying of the porous medium.
...
### 4. ...
exercises/exercise-coupling-ff-pm/interface/ex_interface_coupling_ff-pm.input
View file @ 5df03184
......@@ -34,6 +34,7 @@ Name = darcy
[Darcy.SpatialParams]
Permeability = 1e-6 # m^2
Porosity = 0.4
AlphaBeaversJoseph = 1.0
[Problem]
......
exercises/exercise-coupling-ff-pm/interface/ex_interface_pmproblem.hh
View file @ 5df03184
......@@ -34,7 +34,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>
......
exercises/exercise-coupling-ff-pm/models/CMakeLists.txt 0 → 100644
View file @ 5df03184
add_input_file_links()
# executables for ex_interface_coupling_ff-pm
dune_add_test(NAME ex_models_coupling_ff-pm
SOURCES ex_models_coupling_ff-pm.cc
CMD_ARGS ex_models_coupling_ff-pm.input)
# add tutorial to the common target
add_dependencies(test_exercises ex_models_coupling_ff-pm)
exercises/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.cc 0 → 100644
View file @ 5df03184
// -*- 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 <fstream>
#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/dumuxmessage.hh>
#include <dumux/common/geometry/diameter.hh>
#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/assembly/fvassembler.hh>
#include <dumux/assembly/diffmethod.hh>
#include <dumux/discretization/methods.hh>
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/io/staggeredvtkoutputmodule.hh>
#include <dumux/io/grid/gridmanager.hh>
#include <dumux/multidomain/staggeredtraits.hh>
#include <dumux/multidomain/fvassembler.hh>
#include <dumux/multidomain/newtonsolver.hh>
#include <dumux/multidomain/boundary/stokesdarcy/couplingmanager.hh>
#include "ex_models_pmproblem.hh"
#include "ex_models_ffproblem.hh"
namespace Dumux {
namespace Properties {
SET_PROP(StokesTypeTag, CouplingManager)
{
using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, TTAG(DarcyTypeTag)>;
using type = Dumux::StokesDarcyCouplingManager<Traits>;
};
SET_PROP(DarcyTypeTag, CouplingManager)
{
using Traits = StaggeredMultiDomainTraits<TTAG(StokesTypeTag), TTAG(StokesTypeTag), TypeTag>;
using type = Dumux::StokesDarcyCouplingManager<Traits>;
};
} // end namespace Properties
} // end namespace Dumux
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 StokesTypeTag = TTAG(StokesTypeTag);
using DarcyTypeTag = TTAG(DarcyTypeTag);
// try to create a grid (from the given grid file or the input file)
// for both sub-domains
using DarcyGridManager = Dumux::GridManager<typename GET_PROP_TYPE(DarcyTypeTag, Grid)>;
DarcyGridManager darcyGridManager;
darcyGridManager.init("Darcy"); // pass parameter group
using StokesGridManager = Dumux::GridManager<typename GET_PROP_TYPE(StokesTypeTag, Grid)>;
StokesGridManager stokesGridManager;
stokesGridManager.init("Stokes"); // pass parameter group
// we compute on the leaf grid view
const auto& darcyGridView = darcyGridManager.grid().leafGridView();
const auto& stokesGridView = stokesGridManager.grid().leafGridView();
// create the finite volume grid geometry
using StokesFVGridGeometry = typename GET_PROP_TYPE(StokesTypeTag, FVGridGeometry);
auto stokesFvGridGeometry = std::make_shared<StokesFVGridGeometry>(stokesGridView);
stokesFvGridGeometry->update();
using DarcyFVGridGeometry = typename GET_PROP_TYPE(DarcyTypeTag, FVGridGeometry);
auto darcyFvGridGeometry = std::make_shared<DarcyFVGridGeometry>(darcyGridView);
darcyFvGridGeometry->update();
using Traits = StaggeredMultiDomainTraits<StokesTypeTag, StokesTypeTag, DarcyTypeTag>;
// the coupling manager
using CouplingManager = StokesDarcyCouplingManager<Traits>;
auto couplingManager = std::make_shared<CouplingManager>(stokesFvGridGeometry, darcyFvGridGeometry);
// the indices
constexpr auto stokesCellCenterIdx = CouplingManager::stokesCellCenterIdx;
constexpr auto stokesFaceIdx = CouplingManager::stokesFaceIdx;
constexpr auto darcyIdx = CouplingManager::darcyIdx;
// the problem (initial and boundary conditions)
using StokesProblem = typename GET_PROP_TYPE(StokesTypeTag, Problem);
auto stokesProblem = std::make_shared<StokesProblem>(stokesFvGridGeometry, couplingManager);
using DarcyProblem = typename GET_PROP_TYPE(DarcyTypeTag, Problem);
auto darcyProblem = std::make_shared<DarcyProblem>(darcyFvGridGeometry, couplingManager);
// initialize the fluidsystem (tabulation