From 5df0318497a5aa3dabe9d194ae587051f75c5689 Mon Sep 17 00:00:00 2001
From: Thomas Fetzer <thomas.fetzer@iws.uni-stuttgart.de>
Date: Mon, 16 Jul 2018 07:31:14 +0200
Subject: [PATCH] [coupling-ff-pm] Add model exercise
---
.../{interface => }/1pspatialparams.hh | 4 +-
.../2pspatialparams.hh | 139 ++++++
.../exercise-coupling-ff-pm/CMakeLists.txt | 1 +
exercises/exercise-coupling-ff-pm/README.md | 179 +++++--
.../ex_interface_coupling_ff-pm.input | 1 +
.../interface/ex_interface_pmproblem.hh | 2 +-
.../models/CMakeLists.txt | 9 +
.../models/ex_models_coupling_ff-pm.cc | 296 +++++++++++
.../models/ex_models_coupling_ff-pm.input | 49 ++
.../models/ex_models_ffproblem.hh | 343 +++++++++++++
.../models/ex_models_pmproblem.hh | 405 +++++++++++++++
.../{interface => }/1pspatialparams.hh | 4 +-
.../2pspatialparams.hh | 139 ++++++
.../exercise-coupling-ff-pm/CMakeLists.txt | 1 +
.../ex_interface_coupling_ff-pm.input | 1 +
.../interface/ex_interface_pmproblem.hh | 2 +-
.../models/CMakeLists.txt | 24 +
.../models/ex_models_coupling_ff-pm.cc | 311 ++++++++++++
.../models/ex_models_coupling_ff-pm.input | 49 ++
.../models/ex_models_ffproblem.hh | 343 +++++++++++++
.../models/ex_models_pmproblem.hh | 467 ++++++++++++++++++
21 files changed, 2721 insertions(+), 48 deletions(-)
rename exercises/exercise-coupling-ff-pm/{interface => }/1pspatialparams.hh (96%)
create mode 100644 exercises/exercise-coupling-ff-pm/2pspatialparams.hh
create mode 100644 exercises/exercise-coupling-ff-pm/models/CMakeLists.txt
create mode 100644 exercises/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.cc
create mode 100644 exercises/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.input
create mode 100644 exercises/exercise-coupling-ff-pm/models/ex_models_ffproblem.hh
create mode 100644 exercises/exercise-coupling-ff-pm/models/ex_models_pmproblem.hh
rename exercises/solution/exercise-coupling-ff-pm/{interface => }/1pspatialparams.hh (96%)
create mode 100644 exercises/solution/exercise-coupling-ff-pm/2pspatialparams.hh
create mode 100644 exercises/solution/exercise-coupling-ff-pm/models/CMakeLists.txt
create mode 100644 exercises/solution/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.cc
create mode 100644 exercises/solution/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.input
create mode 100644 exercises/solution/exercise-coupling-ff-pm/models/ex_models_ffproblem.hh
create mode 100644 exercises/solution/exercise-coupling-ff-pm/models/ex_models_pmproblem.hh
diff --git a/exercises/exercise-coupling-ff-pm/interface/1pspatialparams.hh b/exercises/exercise-coupling-ff-pm/1pspatialparams.hh
similarity index 96%
rename from exercises/exercise-coupling-ff-pm/interface/1pspatialparams.hh
rename to exercises/exercise-coupling-ff-pm/1pspatialparams.hh
index 9708c05b..64887e33 100644
--- a/exercises/exercise-coupling-ff-pm/interface/1pspatialparams.hh
+++ b/exercises/exercise-coupling-ff-pm/1pspatialparams.hh
@@ -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_;
};
diff --git a/exercises/exercise-coupling-ff-pm/2pspatialparams.hh b/exercises/exercise-coupling-ff-pm/2pspatialparams.hh
new file mode 100644
index 00000000..1830f3cb
--- /dev/null
+++ b/exercises/exercise-coupling-ff-pm/2pspatialparams.hh
@@ -0,0 +1,139 @@
+// -*- 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
diff --git a/exercises/exercise-coupling-ff-pm/CMakeLists.txt b/exercises/exercise-coupling-ff-pm/CMakeLists.txt
index cc41f095..8c868cd4 100644
--- a/exercises/exercise-coupling-ff-pm/CMakeLists.txt
+++ b/exercises/exercise-coupling-ff-pm/CMakeLists.txt
@@ -1 +1,2 @@
add_subdirectory(interface)
+add_subdirectory(models)
diff --git a/exercises/exercise-coupling-ff-pm/README.md b/exercises/exercise-coupling-ff-pm/README.md
index d1c2d498..a1ea8039 100644
--- a/exercises/exercise-coupling-ff-pm/README.md
+++ b/exercises/exercise-coupling-ff-pm/README.md
@@ -1,16 +1,11 @@
# 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.
-
-
-
-
+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:
+
+
* 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:
+
-### 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. ...
diff --git a/exercises/exercise-coupling-ff-pm/interface/ex_interface_coupling_ff-pm.input b/exercises/exercise-coupling-ff-pm/interface/ex_interface_coupling_ff-pm.input
index bd16b9e8..e3813950 100644
--- a/exercises/exercise-coupling-ff-pm/interface/ex_interface_coupling_ff-pm.input
+++ b/exercises/exercise-coupling-ff-pm/interface/ex_interface_coupling_ff-pm.input
@@ -34,6 +34,7 @@ Name = darcy
[Darcy.SpatialParams]
Permeability = 1e-6 # m^2
+Porosity = 0.4
AlphaBeaversJoseph = 1.0
[Problem]
diff --git a/exercises/exercise-coupling-ff-pm/interface/ex_interface_pmproblem.hh b/exercises/exercise-coupling-ff-pm/interface/ex_interface_pmproblem.hh
index 8013d5aa..7dd08bc3 100644
--- a/exercises/exercise-coupling-ff-pm/interface/ex_interface_pmproblem.hh
+++ b/exercises/exercise-coupling-ff-pm/interface/ex_interface_pmproblem.hh
@@ -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>
diff --git a/exercises/exercise-coupling-ff-pm/models/CMakeLists.txt b/exercises/exercise-coupling-ff-pm/models/CMakeLists.txt
new file mode 100644
index 00000000..1776059b
--- /dev/null
+++ b/exercises/exercise-coupling-ff-pm/models/CMakeLists.txt
@@ -0,0 +1,9 @@
+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)
diff --git a/exercises/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.cc b/exercises/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.cc
new file mode 100644
index 00000000..b54ae19e
--- /dev/null
+++ b/exercises/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.cc
@@ -0,0 +1,296 @@
+// -*- 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)
+ GET_PROP_TYPE(StokesTypeTag, FluidSystem)::init();
+
+ // get some time loop parameters
+ using Scalar = typename GET_PROP_TYPE(StokesTypeTag, Scalar);
+ const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+ const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+ auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+
+ // check if we are about to restart a previously interrupted simulation
+ Scalar restartTime = 0;
+ if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart"))
+ restartTime = getParam<Scalar>("TimeLoop.Restart");
+
+ // instantiate time loop
+ auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(restartTime, dt, tEnd);
+ timeLoop->setMaxTimeStepSize(maxDt);
+ stokesProblem->setTimeLoop(timeLoop);
+ darcyProblem->setTimeLoop(timeLoop);
+
+ // the solution vector
+ Traits::SolutionVector sol;
+ sol[stokesCellCenterIdx].resize(stokesFvGridGeometry->numCellCenterDofs());
+ sol[stokesFaceIdx].resize(stokesFvGridGeometry->numFaceDofs());
+ sol[darcyIdx].resize(darcyFvGridGeometry->numDofs());
+
+ const auto& cellCenterSol = sol[stokesCellCenterIdx];
+ const auto& faceSol = sol[stokesFaceIdx];
+
+ // apply initial solution for instationary problems
+ typename GET_PROP_TYPE(StokesTypeTag, SolutionVector) stokesSol;
+ std::get<0>(stokesSol) = cellCenterSol;
+ std::get<1>(stokesSol) = faceSol;
+ stokesProblem->applyInitialSolution(stokesSol);
+ auto solStokesOld = stokesSol;
+ sol[stokesCellCenterIdx] = stokesSol[stokesCellCenterIdx];
+ sol[stokesFaceIdx] = stokesSol[stokesFaceIdx];
+
+ darcyProblem->applyInitialSolution(sol[darcyIdx]);
+ auto solDarcyOld = sol[darcyIdx];
+
+ auto solOld = sol;
+
+ couplingManager->init(stokesProblem, darcyProblem, sol);
+
+ // the grid variables
+ using StokesGridVariables = typename GET_PROP_TYPE(StokesTypeTag, GridVariables);
+ auto stokesGridVariables = std::make_shared<StokesGridVariables>(stokesProblem, stokesFvGridGeometry);
+ stokesGridVariables->init(stokesSol, solStokesOld);
+ using DarcyGridVariables = typename GET_PROP_TYPE(DarcyTypeTag, GridVariables);
+ auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyFvGridGeometry);
+ darcyGridVariables->init(sol[darcyIdx], solDarcyOld);
+
+ // intialize the vtk output module
+ const auto stokesName = getParam<std::string>("Problem.Name") + "_" + stokesProblem->name();
+ const auto darcyName = getParam<std::string>("Problem.Name") + "_" + darcyProblem->name();
+
+ StaggeredVtkOutputModule<StokesTypeTag> stokesVtkWriter(*stokesProblem, *stokesFvGridGeometry, *stokesGridVariables, stokesSol, stokesName);
+ GET_PROP_TYPE(StokesTypeTag, VtkOutputFields)::init(stokesVtkWriter);
+ stokesVtkWriter.write(0.0);
+
+ VtkOutputModule<DarcyTypeTag> darcyVtkWriter(*darcyProblem, *darcyFvGridGeometry, *darcyGridVariables, sol[darcyIdx], darcyName);
+ GET_PROP_TYPE(DarcyTypeTag, VtkOutputFields)::init(darcyVtkWriter);
+ darcyVtkWriter.write(0.0);
+
+ // intialize the subproblems
+ darcyProblem->init(sol[darcyIdx], *darcyGridVariables);
+
+ // the assembler with time loop for instationary problem
+ using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
+ auto assembler = std::make_shared<Assembler>(std::make_tuple(stokesProblem, stokesProblem, darcyProblem),
+ std::make_tuple(stokesFvGridGeometry->cellCenterFVGridGeometryPtr(),
+ stokesFvGridGeometry->faceFVGridGeometryPtr(),
+ darcyFvGridGeometry),
+ std::make_tuple(stokesGridVariables->cellCenterGridVariablesPtr(),
+ stokesGridVariables->faceGridVariablesPtr(),
+ darcyGridVariables),
+ couplingManager,
+ timeLoop);
+
+ // the linear solver
+ using LinearSolver = UMFPackBackend;
+ auto linearSolver = std::make_shared<LinearSolver>();
+
+ // the primary variable switches used by the sub models and the non-linear solver
+// using PriVarSwitchTuple = std::tuple<NoPrimaryVariableSwitch, NoPrimaryVariableSwitch, NoPrimaryVariableSwitch>;
+ using NewtonSolver = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
+ NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
+
+ // time loop
+ const auto episodeLength = getParam<Scalar>("TimeLoop.EpisodeLength");
+ if (episodeLength > 0.0)
+ timeLoop->setPeriodicCheckPoint(episodeLength);
+ timeLoop->start(); do
+ {
+ // set previous solution for storage evaluations
+ assembler->setPreviousSolution(solOld);
+
+ // solve the non-linear system with time step control
+ nonLinearSolver.solve(sol, *timeLoop);
+
+ // make the new solution the old solution
+ solOld = sol;
+ stokesGridVariables->advanceTimeStep();
+ darcyGridVariables->advanceTimeStep();
+
+ // advance to the time loop to the next step
+ timeLoop->advanceTimeStep();
+
+ // call the postTimeStep routine for output
+ darcyProblem->postTimeStep(sol[darcyIdx], *darcyGridVariables);
+
+ // write vtk output
+ if (timeLoop->isCheckPoint() || timeLoop->finished() || episodeLength < 0.0)
+ {
+ stokesVtkWriter.write(timeLoop->time());
+ darcyVtkWriter.write(timeLoop->time());
+ }
+
+ // report statistics of this time step
+ timeLoop->reportTimeStep();
+
+ // set new dt as suggested by newton solver
+ timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
+
+ } while (!timeLoop->finished());
+
+ timeLoop->finalize(stokesGridView.comm());
+ timeLoop->finalize(darcyGridView.comm());
+
+ ////////////////////////////////////////////////////////////
+ // 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/exercises/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.input b/exercises/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.input
new file mode 100644
index 00000000..d3f68518
--- /dev/null
+++ b/exercises/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.input
@@ -0,0 +1,49 @@
+[TimeLoop]
+DtInitial = 100 # s
+EpisodeLength = -360 # s # 0.25 days
+TEnd = 256000 # s # 2 days
+
+[Stokes.Grid]
+LowerLeft = 0 1
+UpperRight = 1 2
+Cells = 16 16
+
+[Darcy.Grid]
+UpperRight = 1 1
+Cells = 16 16
+
+[Stokes.Problem]
+Name = stokes
+EnableGravity = false
+MoleFraction = 0.0 # -
+Pressure = 1e5 # Pa
+Velocity = 1e-3 # m/s
+
+[Darcy.Problem]
+Name = darcy
+EnableGravity = true
+Saturation = 0.1 # -
+MoleFraction = 0.1 # -
+Pressure = 1e5 # Pa
+
+[Darcy.SpatialParams]
+Permeability = 2.65e-10 # m^2
+Porosity = 0.4 # -
+AlphaBeaversJoseph = 1.0 # -
+# EXNUMBER >= 1
+Swr = 0.005
+Snr = 0.01
+VgAlpha = 6.5e-4
+VgN = 8.0
+
+[Problem]
+Name = models_coupling
+PlotFluxes = false
+PlotStorage = false
+
+[Newton]
+MaxSteps = 12
+MaxRelativeShift = 1e-5
+
+[Vtk]
+AddVelocity = 1
diff --git a/exercises/exercise-coupling-ff-pm/models/ex_models_ffproblem.hh b/exercises/exercise-coupling-ff-pm/models/ex_models_ffproblem.hh
new file mode 100644
index 00000000..f4efefb0
--- /dev/null
+++ b/exercises/exercise-coupling-ff-pm/models/ex_models_ffproblem.hh
@@ -0,0 +1,343 @@
+// -*- 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 NavierStokesTests
+ * \brief A simple Stokes test problem for the staggered grid (Navier-)Stokes model.
+ */
+#ifndef DUMUX_STOKES1P2C_SUBPROBLEM_HH
+#define DUMUX_STOKES1P2C_SUBPROBLEM_HH
+
+#include <dune/grid/yaspgrid.hh>
+
+#include <dumux/material/fluidsystems/h2oair.hh>
+
+#include <dumux/freeflow/navierstokes/problem.hh>
+#include <dumux/discretization/staggered/freeflow/properties.hh>
+#include <dumux/freeflow/compositional/navierstokesncmodel.hh>
+
+namespace Dumux
+{
+template <class TypeTag>
+class StokesSubProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(StokesTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, NavierStokesNC));
+
+// Set the grid type
+SET_TYPE_PROP(StokesTypeTag, Grid, Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<typename GET_PROP_TYPE(TypeTag, Scalar), 2> >);
+
+// Set the fluid system
+SET_TYPE_PROP(StokesTypeTag, FluidSystem, FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>);
+
+// Set phase index (gas)
+SET_INT_PROP(StokesTypeTag, PhaseIdx, GET_PROP_TYPE(TypeTag, FluidSystem)::gasPhaseIdx);
+
+// Do not replace one equation with a total mass balance
+SET_INT_PROP(StokesTypeTag, ReplaceCompEqIdx, 3);
+
+// Use formulation based on mass fractions
+SET_BOOL_PROP(StokesTypeTag, UseMoles, true);
+
+// Set the problem property
+SET_TYPE_PROP(StokesTypeTag, Problem, Dumux::StokesSubProblem<TypeTag> );
+
+SET_BOOL_PROP(StokesTypeTag, EnableFVGridGeometryCache, true);
+SET_BOOL_PROP(StokesTypeTag, EnableGridFluxVariablesCache, true);
+SET_BOOL_PROP(StokesTypeTag, EnableGridVolumeVariablesCache, true);
+
+SET_BOOL_PROP(StokesTypeTag, EnableInertiaTerms, false);
+}
+
+/*!
+ * \ingroup NavierStokesTests
+ * \brief Test problem for the one-phase compositional (Navier-) Stokes problem.
+ *
+ * Horizontal flow from left to right with a parabolic velocity profile.
+ */
+template <class TypeTag>
+class StokesSubProblem : public NavierStokesProblem<TypeTag>
+{
+ using ParentType = NavierStokesProblem<TypeTag>;
+
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
+ using ElementFaceVariables = typename GET_PROP_TYPE(TypeTag, GridFaceVariables)::LocalView;
+ using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
+
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+
+ using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager);
+ using TimeLoopPtr = std::shared_ptr<TimeLoop<Scalar>>;
+
+ static constexpr bool useMoles = GET_PROP_TYPE(TypeTag, ModelTraits)::useMoles();
+
+ static constexpr auto dim = GET_PROP_TYPE(TypeTag, ModelTraits)::dim();
+ static constexpr auto phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
+ static constexpr auto transportCompIdx = 1 - phaseIdx;
+
+public:
+ StokesSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry, std::shared_ptr<CouplingManager> couplingManager)
+ : ParentType(fvGridGeometry, "Stokes"), eps_(1e-6), couplingManager_(couplingManager)
+ {
+ velocity_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Velocity");
+ pressure_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Pressure");
+ moleFraction_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.MoleFraction");
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+
+ bool shouldWriteRestartFile() const
+ { return false; }
+
+ /*!
+ * \brief Return the temperature within the domain in [K].
+ */
+ Scalar temperature() const
+ { return 293.15; }
+
+ /*!
+ * \brief Return the sources within the domain.
+ *
+ * \param globalPos The global position
+ */
+ NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
+ { return NumEqVector(0.0); }
+
+ // \}
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
+
+ /*!
+ * \brief Specifies which kind of boundary condition should be
+ * used for which equation on a given boundary segment.
+ *
+ * \param element The finite element
+ * \param scvf The sub control volume face
+ */
+ BoundaryTypes boundaryTypes(const Element& element,
+ const SubControlVolumeFace& scvf) const
+ {
+ BoundaryTypes values;
+
+ const auto& globalPos = scvf.center();
+
+ if(onLeftBoundary_(globalPos))
+ {
+ values.setDirichlet(transportCompIdx + dim);
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ }
+ else if(onRightBoundary_(globalPos))
+ {
+ values.setDirichlet(Indices::pressureIdx);
+ values.setOutflow(transportCompIdx + dim);
+ }
+ else
+ {
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ values.setNeumann(phaseIdx + dim);
+ values.setNeumann(transportCompIdx + dim);
+ }
+
+ if (couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
+ {
+ values.setCouplingNeumann(phaseIdx + dim);
+ values.setCouplingNeumann(transportCompIdx + dim);
+ values.setCouplingNeumann(Indices::momentumYBalanceIdx);
+ values.setBJS(Indices::momentumXBalanceIdx);
+ }
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Dirichlet control volume.
+ *
+ * \param element The element
+ * \param scvf The subcontrolvolume face
+ */
+ PrimaryVariables dirichletAtPos(const GlobalPosition& pos) const
+ {
+ PrimaryVariables values(0.0);
+ values = initialAtPos(pos);
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Neumann control volume.
+ *
+ * \param element The element for which the Neumann boundary condition is set
+ * \param fvGeomentry The fvGeometry
+ * \param elemVolVars The element volume variables
+ * \param elemFaceVars The element face variables
+ * \param scvf The boundary sub control volume face
+ */
+ NumEqVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFaceVariables& elemFaceVars,
+ const SubControlVolumeFace& scvf) const
+ {
+ PrimaryVariables values(0.0);
+
+ if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
+ {
+ values[Indices::momentumYBalanceIdx] = couplingManager().couplingData().momentumCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf);
+
+ const auto massFlux = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf);
+ values[Indices::conti0EqIdx] = massFlux[0];
+ values[Indices::conti0EqIdx + 1] = massFlux[1];
+ }
+ return values;
+ }
+
+ // \}
+
+ /*!
+ * \brief Set the coupling manager
+ */
+ void setCouplingManager(std::shared_ptr<CouplingManager> cm)
+ { couplingManager_ = cm; }
+
+ /*!
+ * \brief Get the coupling manager
+ */
+ const CouplingManager& couplingManager() const
+ { return *couplingManager_; }
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * \param globalPos The global position
+ */
+ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+ {
+ FluidState fluidState;
+ updateFluidStateForBC_(fluidState, pressure_);
+ const Scalar density = FluidSystem::density(fluidState, phaseIdx);
+
+ PrimaryVariables values(0.0);
+ values[Indices::pressureIdx] = pressure_ + density*this->gravity()[1]*(globalPos[1] - this->fvGridGeometry().bBoxMin()[1]);
+ values[transportCompIdx + dim] = moleFraction_;
+ values[Indices::velocityXIdx] = 4.0 * velocity_ * (globalPos[1] - this->fvGridGeometry().bBoxMin()[1])
+ * (this->fvGridGeometry().bBoxMax()[1] - globalPos[1])
+ / (height_() * height_());
+
+ return values;
+ }
+
+ void setTimeLoop(TimeLoopPtr timeLoop)
+ { timeLoop_ = timeLoop; }
+
+ /*!
+ * \brief Returns the intrinsic permeability of required as input parameter for the Beavers-Joseph-Saffman boundary condition
+ */
+ Scalar permeability(const SubControlVolumeFace& scvf) const
+ {
+ return couplingManager().couplingData().darcyPermeability(scvf);
+ }
+
+ /*!
+ * \brief Returns the alpha value required as input parameter for the Beavers-Joseph-Saffman boundary condition
+ */
+ Scalar alphaBJ(const SubControlVolumeFace& scvf) const
+ {
+ return couplingManager().problem(CouplingManager::darcyIdx).spatialParams().beaversJosephCoeffAtPos(scvf.center());
+ }
+
+ // \}
+
+private:
+ bool onLeftBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_; }
+
+ bool onRightBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_; }
+
+ bool onLowerBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] < this->fvGridGeometry().bBoxMin()[1] + eps_; }
+
+ bool onUpperBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_; }
+
+ //! \brief updates the fluid state to obtain required quantities for IC/BC
+ void updateFluidStateForBC_(FluidState& fluidState, const Scalar pressure) const
+ {
+ fluidState.setTemperature(temperature());
+ fluidState.setPressure(phaseIdx, pressure);
+ fluidState.setSaturation(phaseIdx, 1.0);
+ fluidState.setMoleFraction(phaseIdx, transportCompIdx, moleFraction_);
+ fluidState.setMoleFraction(phaseIdx, phaseIdx, 1.0 - moleFraction_);
+
+ typename FluidSystem::ParameterCache paramCache;
+ paramCache.updatePhase(fluidState, phaseIdx);
+
+ const Scalar density = FluidSystem::density(fluidState, paramCache, phaseIdx);
+ fluidState.setDensity(phaseIdx, density);
+
+ const Scalar molarDensity = FluidSystem::molarDensity(fluidState, paramCache, phaseIdx);
+ fluidState.setMolarDensity(phaseIdx, molarDensity);
+
+ const Scalar enthalpy = FluidSystem::enthalpy(fluidState, paramCache, phaseIdx);
+ fluidState.setEnthalpy(phaseIdx, enthalpy);
+ }
+ // the height of the free-flow domain
+ const Scalar height_() const
+ { return this->fvGridGeometry().bBoxMax()[1] - this->fvGridGeometry().bBoxMin()[1]; }
+
+ Scalar eps_;
+
+ Scalar velocity_;
+ Scalar pressure_;
+ Scalar moleFraction_;
+
+ TimeLoopPtr timeLoop_;
+
+ std::shared_ptr<CouplingManager> couplingManager_;
+};
+} //end namespace
+
+#endif // DUMUX_STOKES1P2C_SUBPROBLEM_HH
diff --git a/exercises/exercise-coupling-ff-pm/models/ex_models_pmproblem.hh b/exercises/exercise-coupling-ff-pm/models/ex_models_pmproblem.hh
new file mode 100644
index 00000000..d95c094e
--- /dev/null
+++ b/exercises/exercise-coupling-ff-pm/models/ex_models_pmproblem.hh
@@ -0,0 +1,405 @@
+// -*- 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 simple Darcy test problem (cell-centered finite volume method).
+ */
+#ifndef DUMUX_DARCY_SUBPROBLEM_HH
+#define DUMUX_DARCY_SUBPROBLEM_HH
+
+#include <dune/grid/yaspgrid.hh>
+
+#include <dumux/discretization/cellcentered/tpfa/properties.hh>
+#include <dumux/io/gnuplotinterface.hh>
+#include <dumux/material/fluidsystems/h2oair.hh>
+#include <dumux/material/fluidmatrixinteractions/diffusivityconstanttortuosity.hh>
+
+#include <dumux/porousmediumflow/problem.hh>
+
+#include <dumux/porousmediumflow/1pnc/model.hh>
+#include "../1pspatialparams.hh"
+
+namespace Dumux
+{
+template <class TypeTag>
+class DarcySubProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(DarcyTypeTag, INHERITS_FROM(CCTpfaModel, OnePNC));
+
+// Set the problem property
+SET_TYPE_PROP(DarcyTypeTag, Problem, Dumux::DarcySubProblem<TypeTag>);
+
+// The fluid system
+SET_TYPE_PROP(DarcyTypeTag, FluidSystem, FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>);
+
+// Use gas as phase
+SET_INT_PROP(DarcyTypeTag, PhaseIdx, GET_PROP_TYPE(TypeTag, FluidSystem)::gasPhaseIdx);
+
+// Use moles
+SET_BOOL_PROP(DarcyTypeTag, UseMoles, true);
+
+// Do not replace one equation with a total mass balance
+SET_INT_PROP(DarcyTypeTag, ReplaceCompEqIdx, 3);
+
+//! Use a model with constant tortuosity for the effective diffusivity
+SET_TYPE_PROP(DarcyTypeTag, EffectiveDiffusivityModel,
+ DiffusivityConstantTortuosity<typename GET_PROP_TYPE(TypeTag, Scalar)>);
+
+// Set the grid type
+SET_TYPE_PROP(DarcyTypeTag, Grid, Dune::YaspGrid<2>);
+
+// Set the spatial paramaters type
+SET_TYPE_PROP(DarcyTypeTag, SpatialParams, OnePSpatialParams<TypeTag>);
+}
+
+template <class TypeTag>
+class DarcySubProblem : public PorousMediumFlowProblem<TypeTag>
+{
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+
+ // copy some indices for convenience
+ using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+ enum {
+ // grid and world dimension
+ dim = GridView::dimension,
+ dimworld = GridView::dimensionworld,
+
+ // primary variable indices
+ conti0EqIdx = Indices::conti0EqIdx,
+ pressureIdx = Indices::pressureIdx,
+ phaseIdx = Indices::fluidSystemPhaseIdx,
+ transportCompIdx = 1-phaseIdx
+ };
+
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
+
+ using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager);
+ using TimeLoopPtr = std::shared_ptr<TimeLoop<Scalar>>;
+
+public:
+ DarcySubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry,
+ std::shared_ptr<CouplingManager> couplingManager)
+ : ParentType(fvGridGeometry, "Darcy"), eps_(1e-7), couplingManager_(couplingManager)
+ {
+ moleFraction_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.MoleFraction");
+ pressure_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Pressure");
+
+ // initialize output file
+ plotFluxes_ = getParamFromGroup<bool>(this->paramGroup(), "Problem.PlotFluxes", false);
+ plotStorage_ = getParamFromGroup<bool>(this->paramGroup(), "Problem.PlotStorage", false);
+ storageFileName_ = "storage_" + getParam<std::string>("Problem.Name") + "_" + this->name() + ".csv";
+ storageFile_.open(storageFileName_);
+ storageFile_ << "#Time[s]" << ";"
+ << "WaterMass[kg]" << ";"
+ << "WaterMassLoss[kg]" << ";"
+ << "EvaporationRate[mm/d]"
+ << std::endl;
+ }
+
+ /*!
+ * \name Simulation steering
+ */
+ // \{
+
+ /*!
+ * \brief Initialize the problem.
+ */
+ template<class SolutionVector, class GridVariables>
+ void init(const SolutionVector& curSol,
+ const GridVariables& gridVariables)
+ { }
+
+ template<class SolutionVector, class GridVariables>
+ void postTimeStep(const SolutionVector& curSol,
+ const GridVariables& gridVariables)
+
+ {
+ evaluateWaterMassStorageTerm(curSol, gridVariables);
+ evaluateInterfaceFluxes(curSol, gridVariables);
+
+ gnuplotStorage_.resetPlot();
+ gnuplotStorage_.setDatafileSeparator(';');
+ gnuplotStorage_.setXlabel("time [d]");
+ gnuplotStorage_.setXRange(0.0, getParam<Scalar>("TimeLoop.TEnd"));
+ gnuplotStorage_.setYlabel("evaporation rate [mm/d]");
+ gnuplotStorage_.setOption("set yrange [0.0:]");
+ gnuplotStorage_.setOption("set y2label 'cumulative mass loss'");
+ gnuplotStorage_.setOption("set y2range [0.0:0.5]");
+ gnuplotStorage_.setOption("set y2range [0.0:0.5]");
+ gnuplotStorage_.addFileToPlot(storageFileName_, "using 1:4 with lines title 'evaporation rate'");
+ gnuplotStorage_.addFileToPlot(storageFileName_, "using 1:3 axes x1y2 with lines title 'cumulative mass loss'");
+ if (plotStorage_)
+ gnuplotStorage_.plot("temp");
+ }
+
+ template<class SolutionVector, class GridVariables>
+ Scalar evaluateWaterMassStorageTerm(const SolutionVector& curSol,
+ const GridVariables& gridVariables)
+
+ {
+ // compute the mass in the entire domain
+ Scalar waterMass = 0.0;
+
+ for (const auto& element : elements(this->fvGridGeometry().gridView()))
+ {
+ auto fvGeometry = localView(this->fvGridGeometry());
+ fvGeometry.bindElement(element);
+
+ auto elemVolVars = localView(gridVariables.curGridVolVars());
+ elemVolVars.bindElement(element, fvGeometry, curSol);
+
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ const auto& volVars = elemVolVars[scv];
+ for(int phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
+ {
+ // insert calculation of the water mass here
+ waterMass += 0.0;
+ }
+ }
+ }
+
+ Scalar cumMassLoss = initialWaterContent_ - waterMass;
+ Scalar evaporationRate = (lastWaterMass_ - waterMass) * 86400
+ / (this->fvGridGeometry().bBoxMax()[0] - this->fvGridGeometry().bBoxMin()[0])
+ / timeLoop_->timeStepSize();
+ lastWaterMass_ = waterMass;
+
+ storageFile_ << timeLoop_->time() << ";"
+ << waterMass << ";"
+ << cumMassLoss << ";"
+ << evaporationRate
+ << std::endl;
+
+ return waterMass;
+ }
+
+ template<class SolutionVector, class GridVariables>
+ void evaluateInterfaceFluxes(const SolutionVector& curSol,
+ const GridVariables& gridVariables)
+
+ {
+ std::vector<Scalar> x;
+ std::vector<Scalar> y;
+
+ for (const auto& element : elements(this->fvGridGeometry().gridView()))
+ {
+ auto fvGeometry = localView(this->fvGridGeometry());
+ fvGeometry.bindElement(element);
+
+ auto elemVolVars = localView(gridVariables.curGridVolVars());
+ elemVolVars.bindElement(element, fvGeometry, curSol);
+
+ for (auto&& scvf : scvfs(fvGeometry))
+ {
+ if (!couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
+ continue;
+
+ // NOTE: binding the coupling context is necessary
+ couplingManager_->bindCouplingContext(CouplingManager::darcyIdx, element);
+ NumEqVector flux(0.0); // use "massCouplingCondition" from the couplingManager here
+
+ x.push_back(scvf.center()[0]);
+ y.push_back(flux[transportCompIdx]);
+ }
+ }
+
+ gnuplotInterfaceFluxes_.resetPlot();
+ gnuplotInterfaceFluxes_.setXlabel("x-position [m]");
+ gnuplotInterfaceFluxes_.setXRange(this->fvGridGeometry().bBoxMin()[0], this->fvGridGeometry().bBoxMax()[0]);
+ gnuplotInterfaceFluxes_.setYlabel("flux [kg/(m^2 s)]");
+ gnuplotInterfaceFluxes_.setYRange(-5e-4, 0.0);
+ gnuplotInterfaceFluxes_.setOption("set label 'time: " + std::to_string(timeLoop_->time()/86400.) + "d' at graph 0.8,0.8 ");
+ std::string fluxFileName = "flux_" + std::to_string(timeLoop_->timeStepIndex()) +
+ "_" + getParam<std::string>("Problem.Name") + "_" + this->name() + ".csv";
+ gnuplotInterfaceFluxes_.addDataSetToPlot(x, y, fluxFileName, "with lines title 'water mass flux'");
+ if (plotFluxes_)
+ gnuplotInterfaceFluxes_.plot("flux_" + std::to_string(timeLoop_->timeStepIndex()));
+ }
+
+ /*!
+ * \brief Returns true if a restart file should be written to
+ * disk.
+ */
+ bool shouldWriteRestartFile() const
+ { return false; }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ bool shouldWriteOutput() const // define output
+ { return true; }
+
+ /*!
+ * \brief Return the temperature within the domain in [K].
+ *
+ */
+ Scalar temperature() const
+ { return 293.15; }
+ // \}
+
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
+
+ /*!
+ * \brief Specifies which kind of boundary condition should be
+ * used for which equation on a given boundary control volume.
+ *
+ * \param element The element
+ * \param scvf The boundary sub control volume face
+ */
+ BoundaryTypes boundaryTypes(const Element& element, const SubControlVolumeFace& scvf) const
+ {
+ BoundaryTypes values;
+ values.setAllNeumann();
+
+ if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
+ values.setAllCouplingNeumann();
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Neumann control volume.
+ *
+ * \param element The element for which the Neumann boundary condition is set
+ * \param fvGeomentry The fvGeometry
+ * \param elemVolVars The element volume variables
+ * \param scvf The boundary sub control volume face
+ *
+ * For this method, the \a values variable stores primary variables.
+ */
+ template<class ElementVolumeVariables>
+ NumEqVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf) const
+ {
+ NumEqVector values(0.0);
+
+ if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
+ values = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, scvf);
+
+ return values;
+ }
+
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+ /*!
+ * \brief Evaluate the source term for all phases within a given
+ * sub-control-volume.
+ *
+ * \param element The element for which the source term is set
+ * \param fvGeomentry The fvGeometry
+ * \param elemVolVars The element volume variables
+ * \param scv The subcontrolvolume
+ */
+ template<class ElementVolumeVariables>
+ NumEqVector source(const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolume &scv) const
+ { return NumEqVector(0.0); }
+
+ // \}
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * \param element The element
+ *
+ * For this method, the \a priVars parameter stores primary
+ * variables.
+ */
+ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+ {
+ PrimaryVariables values(0.0);
+ values[transportCompIdx] = moleFraction_;
+ values[pressureIdx] = pressure_;
+ return values;
+ }
+
+ // \}
+
+ //! Set the coupling manager
+ void setCouplingManager(std::shared_ptr<CouplingManager> cm)
+ { couplingManager_ = cm; }
+
+ //! Get the coupling manager
+ const CouplingManager& couplingManager() const
+ { return *couplingManager_; }
+
+ void setTimeLoop(TimeLoopPtr timeLoop)
+ { timeLoop_ = timeLoop; }
+
+private:
+ bool onLeftBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_; }
+
+ bool onRightBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_; }
+
+ bool onLowerBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] < this->fvGridGeometry().bBoxMin()[1] + eps_; }
+
+ bool onUpperBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_; }
+
+ Scalar eps_;
+ Scalar moleFraction_;
+ Scalar pressure_;
+
+ Scalar initialWaterContent_ = 0.0;
+ Scalar lastWaterMass_ = 0.0;
+
+ TimeLoopPtr timeLoop_;
+ std::shared_ptr<CouplingManager> couplingManager_;
+
+ std::string storageFileName_;
+ std::ofstream storageFile_;
+ bool plotFluxes_;
+ bool plotStorage_;
+ Dumux::GnuplotInterface<Scalar> gnuplotInterfaceFluxes_;
+ Dumux::GnuplotInterface<Scalar> gnuplotStorage_;
+};
+} //end namespace
+
+#endif //DUMUX_DARCY_SUBPROBLEM_HH
diff --git a/exercises/solution/exercise-coupling-ff-pm/interface/1pspatialparams.hh b/exercises/solution/exercise-coupling-ff-pm/1pspatialparams.hh
similarity index 96%
rename from exercises/solution/exercise-coupling-ff-pm/interface/1pspatialparams.hh
rename to exercises/solution/exercise-coupling-ff-pm/1pspatialparams.hh
index 9708c05b..64887e33 100644
--- a/exercises/solution/exercise-coupling-ff-pm/interface/1pspatialparams.hh
+++ b/exercises/solution/exercise-coupling-ff-pm/1pspatialparams.hh
@@ -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_;
};
diff --git a/exercises/solution/exercise-coupling-ff-pm/2pspatialparams.hh b/exercises/solution/exercise-coupling-ff-pm/2pspatialparams.hh
new file mode 100644
index 00000000..1830f3cb
--- /dev/null
+++ b/exercises/solution/exercise-coupling-ff-pm/2pspatialparams.hh
@@ -0,0 +1,139 @@
+// -*- 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
diff --git a/exercises/solution/exercise-coupling-ff-pm/CMakeLists.txt b/exercises/solution/exercise-coupling-ff-pm/CMakeLists.txt
index cc41f095..8c868cd4 100644
--- a/exercises/solution/exercise-coupling-ff-pm/CMakeLists.txt
+++ b/exercises/solution/exercise-coupling-ff-pm/CMakeLists.txt
@@ -1 +1,2 @@
add_subdirectory(interface)
+add_subdirectory(models)
diff --git a/exercises/solution/exercise-coupling-ff-pm/interface/ex_interface_coupling_ff-pm.input b/exercises/solution/exercise-coupling-ff-pm/interface/ex_interface_coupling_ff-pm.input
index bd16b9e8..e3813950 100644
--- a/exercises/solution/exercise-coupling-ff-pm/interface/ex_interface_coupling_ff-pm.input
+++ b/exercises/solution/exercise-coupling-ff-pm/interface/ex_interface_coupling_ff-pm.input
@@ -34,6 +34,7 @@ Name = darcy
[Darcy.SpatialParams]
Permeability = 1e-6 # m^2
+Porosity = 0.4
AlphaBeaversJoseph = 1.0
[Problem]
diff --git a/exercises/solution/exercise-coupling-ff-pm/interface/ex_interface_pmproblem.hh b/exercises/solution/exercise-coupling-ff-pm/interface/ex_interface_pmproblem.hh
index eb1c5a20..7a5c4baa 100644
--- a/exercises/solution/exercise-coupling-ff-pm/interface/ex_interface_pmproblem.hh
+++ b/exercises/solution/exercise-coupling-ff-pm/interface/ex_interface_pmproblem.hh
@@ -35,7 +35,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>
diff --git a/exercises/solution/exercise-coupling-ff-pm/models/CMakeLists.txt b/exercises/solution/exercise-coupling-ff-pm/models/CMakeLists.txt
new file mode 100644
index 00000000..f7e7af5f
--- /dev/null
+++ b/exercises/solution/exercise-coupling-ff-pm/models/CMakeLists.txt
@@ -0,0 +1,24 @@
+add_input_file_links()
+
+dune_add_test(NAME orig_ex_models_coupling_ff-pm
+ SOURCES ex_models_coupling_ff-pm.cc
+ COMPILE_DEFINITIONS EXNUMBER=0
+ CMD_ARGS ex_models_coupling_ff-pm.input)
+
+dune_add_test(NAME sol_a_ex_models_coupling_ff-pm
+ SOURCES ex_models_coupling_ff-pm.cc
+ COMPILE_DEFINITIONS EXNUMBER=1
+ CMD_ARGS ex_models_coupling_ff-pm.input)
+
+dune_add_test(NAME sol_b_ex_models_coupling_ff-pm
+ SOURCES ex_models_coupling_ff-pm.cc
+ COMPILE_DEFINITIONS EXNUMBER=2
+ CMD_ARGS ex_models_coupling_ff-pm.input)
+
+dune_add_test(NAME sol_c_ex_models_coupling_ff-pm
+ SOURCES ex_models_coupling_ff-pm.cc
+ COMPILE_DEFINITIONS EXNUMBER=3
+ CMD_ARGS ex_models_coupling_ff-pm.input)
+
+# add exercise to the common target
+add_dependencies(test_exercises sol_a_ex_models_coupling_ff-pm sol_b_ex_models_coupling_ff-pm sol_b_ex_models_coupling_ff-pm)
diff --git a/exercises/solution/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.cc b/exercises/solution/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.cc
new file mode 100644
index 00000000..cdb7abdd
--- /dev/null
+++ b/exercises/solution/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.cc
@@ -0,0 +1,311 @@
+// -*- 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>
+#if EXNUMBER >= 3
+#include <dumux/multidomain/privarswitchnewtonsolver.hh>
+#else
+#include <dumux/multidomain/newtonsolver.hh>
+#endif
+
+#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)
+ GET_PROP_TYPE(StokesTypeTag, FluidSystem)::init();
+
+ // get some time loop parameters
+ using Scalar = typename GET_PROP_TYPE(StokesTypeTag, Scalar);
+ const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+ const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+ auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+
+ // check if we are about to restart a previously interrupted simulation
+ Scalar restartTime = 0;
+ if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart"))
+ restartTime = getParam<Scalar>("TimeLoop.Restart");
+
+ // instantiate time loop
+ auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(restartTime, dt, tEnd);
+ timeLoop->setMaxTimeStepSize(maxDt);
+ stokesProblem->setTimeLoop(timeLoop);
+ darcyProblem->setTimeLoop(timeLoop);
+
+ // the solution vector
+ Traits::SolutionVector sol;
+ sol[stokesCellCenterIdx].resize(stokesFvGridGeometry->numCellCenterDofs());
+ sol[stokesFaceIdx].resize(stokesFvGridGeometry->numFaceDofs());
+ sol[darcyIdx].resize(darcyFvGridGeometry->numDofs());
+
+ const auto& cellCenterSol = sol[stokesCellCenterIdx];
+ const auto& faceSol = sol[stokesFaceIdx];
+
+ // apply initial solution for instationary problems
+ typename GET_PROP_TYPE(StokesTypeTag, SolutionVector) stokesSol;
+ std::get<0>(stokesSol) = cellCenterSol;
+ std::get<1>(stokesSol) = faceSol;
+ stokesProblem->applyInitialSolution(stokesSol);
+ auto solStokesOld = stokesSol;
+ sol[stokesCellCenterIdx] = stokesSol[stokesCellCenterIdx];
+ sol[stokesFaceIdx] = stokesSol[stokesFaceIdx];
+
+ darcyProblem->applyInitialSolution(sol[darcyIdx]);
+ auto solDarcyOld = sol[darcyIdx];
+
+ auto solOld = sol;
+
+ couplingManager->init(stokesProblem, darcyProblem, sol);
+
+ // the grid variables
+ using StokesGridVariables = typename GET_PROP_TYPE(StokesTypeTag, GridVariables);
+ auto stokesGridVariables = std::make_shared<StokesGridVariables>(stokesProblem, stokesFvGridGeometry);
+ stokesGridVariables->init(stokesSol, solStokesOld);
+ using DarcyGridVariables = typename GET_PROP_TYPE(DarcyTypeTag, GridVariables);
+ auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyFvGridGeometry);
+ darcyGridVariables->init(sol[darcyIdx], solDarcyOld);
+
+ // intialize the vtk output module
+#if EXNUMBER >= 1
+ const std::array<std::string, 3> part = {"a", "b", "c"};
+ const auto stokesName = "sol_" + part[EXNUMBER-1] + "_" + getParam<std::string>("Problem.Name") + "_" + stokesProblem->name();
+ const auto darcyName = "sol_" + part[EXNUMBER-1] + "_" + getParam<std::string>("Problem.Name") + "_" + darcyProblem->name();
+#else
+ const auto stokesName = "orig_" + getParam<std::string>("Problem.Name") + "_" + stokesProblem->name();
+ const auto darcyName = "orig_" + getParam<std::string>("Problem.Name") + "_" + darcyProblem->name();
+#endif
+
+ StaggeredVtkOutputModule<StokesTypeTag> stokesVtkWriter(*stokesProblem, *stokesFvGridGeometry, *stokesGridVariables, stokesSol, stokesName);
+ GET_PROP_TYPE(StokesTypeTag, VtkOutputFields)::init(stokesVtkWriter);
+ stokesVtkWriter.write(0.0);
+
+ VtkOutputModule<DarcyTypeTag> darcyVtkWriter(*darcyProblem, *darcyFvGridGeometry, *darcyGridVariables, sol[darcyIdx], darcyName);
+ GET_PROP_TYPE(DarcyTypeTag, VtkOutputFields)::init(darcyVtkWriter);
+ darcyVtkWriter.write(0.0);
+
+ // intialize the subproblems
+ darcyProblem->init(sol[darcyIdx], *darcyGridVariables);
+
+ // the assembler with time loop for instationary problem
+ using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
+ auto assembler = std::make_shared<Assembler>(std::make_tuple(stokesProblem, stokesProblem, darcyProblem),
+ std::make_tuple(stokesFvGridGeometry->cellCenterFVGridGeometryPtr(),
+ stokesFvGridGeometry->faceFVGridGeometryPtr(),
+ darcyFvGridGeometry),
+ std::make_tuple(stokesGridVariables->cellCenterGridVariablesPtr(),
+ stokesGridVariables->faceGridVariablesPtr(),
+ darcyGridVariables),
+ couplingManager,
+ timeLoop);
+
+ // the linear solver
+ using LinearSolver = UMFPackBackend;
+ auto linearSolver = std::make_shared<LinearSolver>();
+
+ // the primary variable switches used by the sub models and the non-linear solver
+#if EXNUMBER >= 3
+ using PriVarSwitchTuple = std::tuple<NoPrimaryVariableSwitch, NoPrimaryVariableSwitch, typename GET_PROP_TYPE(DarcyTypeTag, PrimaryVariableSwitch)>;
+ using NewtonSolver = MultiDomainPriVarSwitchNewtonSolver<Assembler, LinearSolver, CouplingManager, PriVarSwitchTuple>;
+#else
+// using PriVarSwitchTuple = std::tuple<NoPrimaryVariableSwitch, NoPrimaryVariableSwitch, NoPrimaryVariableSwitch>;
+ using NewtonSolver = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
+#endif
+ NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
+
+ // time loop
+ const auto episodeLength = getParam<Scalar>("TimeLoop.EpisodeLength");
+ if (episodeLength > 0.0)
+ timeLoop->setPeriodicCheckPoint(episodeLength);
+ timeLoop->start(); do
+ {
+ // set previous solution for storage evaluations
+ assembler->setPreviousSolution(solOld);
+
+ // solve the non-linear system with time step control
+ nonLinearSolver.solve(sol, *timeLoop);
+
+ // make the new solution the old solution
+ solOld = sol;
+ stokesGridVariables->advanceTimeStep();
+ darcyGridVariables->advanceTimeStep();
+
+ // advance to the time loop to the next step
+ timeLoop->advanceTimeStep();
+
+ // call the postTimeStep routine for output
+ darcyProblem->postTimeStep(sol[darcyIdx], *darcyGridVariables);
+
+ // write vtk output
+ if (timeLoop->isCheckPoint() || timeLoop->finished() || episodeLength < 0.0)
+ {
+ stokesVtkWriter.write(timeLoop->time());
+ darcyVtkWriter.write(timeLoop->time());
+ }
+
+ // report statistics of this time step
+ timeLoop->reportTimeStep();
+
+ // set new dt as suggested by newton solver
+ timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
+
+ } while (!timeLoop->finished());
+
+ timeLoop->finalize(stokesGridView.comm());
+ timeLoop->finalize(darcyGridView.comm());
+
+ ////////////////////////////////////////////////////////////
+ // 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/exercises/solution/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.input b/exercises/solution/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.input
new file mode 100644
index 00000000..8d08fe83
--- /dev/null
+++ b/exercises/solution/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.input
@@ -0,0 +1,49 @@
+[TimeLoop]
+DtInitial = 100 # s
+EpisodeLength = -360 # s # 0.25 days
+TEnd = 256000 # s # 2 days
+
+[Stokes.Grid]
+LowerLeft = 0 1
+UpperRight = 1 2
+Cells = 16 16
+
+[Darcy.Grid]
+UpperRight = 1 1
+Cells = 16 16
+
+[Stokes.Problem]
+Name = stokes
+EnableGravity = false
+MoleFraction = 0.0 # -
+Pressure = 1e5 # Pa
+Velocity = 1e-3 # m/s
+
+[Darcy.Problem]
+Name = darcy
+EnableGravity = true
+Saturation = 0.1 # -
+MoleFraction = 0.1 # -
+Pressure = 1e5 # Pa
+
+[Darcy.SpatialParams]
+Permeability = 2.65e-10 # m^2
+Porosity = 0.4 # -
+AlphaBeaversJoseph = 1.0 # -
+# EXNUMBER >= 1
+Swr = 0.005
+Snr = 0.01
+VgAlpha = 6.5e-4
+VgN = 8.0
+
+[Problem]
+Name = ex_models_coupling
+PlotFluxes = false
+PlotStorage = false
+
+[Newton]
+MaxSteps = 12
+MaxRelativeShift = 1e-5
+
+[Vtk]
+AddVelocity = 1
diff --git a/exercises/solution/exercise-coupling-ff-pm/models/ex_models_ffproblem.hh b/exercises/solution/exercise-coupling-ff-pm/models/ex_models_ffproblem.hh
new file mode 100644
index 00000000..f4efefb0
--- /dev/null
+++ b/exercises/solution/exercise-coupling-ff-pm/models/ex_models_ffproblem.hh
@@ -0,0 +1,343 @@
+// -*- 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 NavierStokesTests
+ * \brief A simple Stokes test problem for the staggered grid (Navier-)Stokes model.
+ */
+#ifndef DUMUX_STOKES1P2C_SUBPROBLEM_HH
+#define DUMUX_STOKES1P2C_SUBPROBLEM_HH
+
+#include <dune/grid/yaspgrid.hh>
+
+#include <dumux/material/fluidsystems/h2oair.hh>
+
+#include <dumux/freeflow/navierstokes/problem.hh>
+#include <dumux/discretization/staggered/freeflow/properties.hh>
+#include <dumux/freeflow/compositional/navierstokesncmodel.hh>
+
+namespace Dumux
+{
+template <class TypeTag>
+class StokesSubProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(StokesTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, NavierStokesNC));
+
+// Set the grid type
+SET_TYPE_PROP(StokesTypeTag, Grid, Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<typename GET_PROP_TYPE(TypeTag, Scalar), 2> >);
+
+// Set the fluid system
+SET_TYPE_PROP(StokesTypeTag, FluidSystem, FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>);
+
+// Set phase index (gas)
+SET_INT_PROP(StokesTypeTag, PhaseIdx, GET_PROP_TYPE(TypeTag, FluidSystem)::gasPhaseIdx);
+
+// Do not replace one equation with a total mass balance
+SET_INT_PROP(StokesTypeTag, ReplaceCompEqIdx, 3);
+
+// Use formulation based on mass fractions
+SET_BOOL_PROP(StokesTypeTag, UseMoles, true);
+
+// Set the problem property
+SET_TYPE_PROP(StokesTypeTag, Problem, Dumux::StokesSubProblem<TypeTag> );
+
+SET_BOOL_PROP(StokesTypeTag, EnableFVGridGeometryCache, true);
+SET_BOOL_PROP(StokesTypeTag, EnableGridFluxVariablesCache, true);
+SET_BOOL_PROP(StokesTypeTag, EnableGridVolumeVariablesCache, true);
+
+SET_BOOL_PROP(StokesTypeTag, EnableInertiaTerms, false);
+}
+
+/*!
+ * \ingroup NavierStokesTests
+ * \brief Test problem for the one-phase compositional (Navier-) Stokes problem.
+ *
+ * Horizontal flow from left to right with a parabolic velocity profile.
+ */
+template <class TypeTag>
+class StokesSubProblem : public NavierStokesProblem<TypeTag>
+{
+ using ParentType = NavierStokesProblem<TypeTag>;
+
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
+ using ElementFaceVariables = typename GET_PROP_TYPE(TypeTag, GridFaceVariables)::LocalView;
+ using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
+
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+
+ using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager);
+ using TimeLoopPtr = std::shared_ptr<TimeLoop<Scalar>>;
+
+ static constexpr bool useMoles = GET_PROP_TYPE(TypeTag, ModelTraits)::useMoles();
+
+ static constexpr auto dim = GET_PROP_TYPE(TypeTag, ModelTraits)::dim();
+ static constexpr auto phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
+ static constexpr auto transportCompIdx = 1 - phaseIdx;
+
+public:
+ StokesSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry, std::shared_ptr<CouplingManager> couplingManager)
+ : ParentType(fvGridGeometry, "Stokes"), eps_(1e-6), couplingManager_(couplingManager)
+ {
+ velocity_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Velocity");
+ pressure_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Pressure");
+ moleFraction_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.MoleFraction");
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+
+ bool shouldWriteRestartFile() const
+ { return false; }
+
+ /*!
+ * \brief Return the temperature within the domain in [K].
+ */
+ Scalar temperature() const
+ { return 293.15; }
+
+ /*!
+ * \brief Return the sources within the domain.
+ *
+ * \param globalPos The global position
+ */
+ NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
+ { return NumEqVector(0.0); }
+
+ // \}
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
+
+ /*!
+ * \brief Specifies which kind of boundary condition should be
+ * used for which equation on a given boundary segment.
+ *
+ * \param element The finite element
+ * \param scvf The sub control volume face
+ */
+ BoundaryTypes boundaryTypes(const Element& element,
+ const SubControlVolumeFace& scvf) const
+ {
+ BoundaryTypes values;
+
+ const auto& globalPos = scvf.center();
+
+ if(onLeftBoundary_(globalPos))
+ {
+ values.setDirichlet(transportCompIdx + dim);
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ }
+ else if(onRightBoundary_(globalPos))
+ {
+ values.setDirichlet(Indices::pressureIdx);
+ values.setOutflow(transportCompIdx + dim);
+ }
+ else
+ {
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ values.setNeumann(phaseIdx + dim);
+ values.setNeumann(transportCompIdx + dim);
+ }
+
+ if (couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
+ {
+ values.setCouplingNeumann(phaseIdx + dim);
+ values.setCouplingNeumann(transportCompIdx + dim);
+ values.setCouplingNeumann(Indices::momentumYBalanceIdx);
+ values.setBJS(Indices::momentumXBalanceIdx);
+ }
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Dirichlet control volume.
+ *
+ * \param element The element
+ * \param scvf The subcontrolvolume face
+ */
+ PrimaryVariables dirichletAtPos(const GlobalPosition& pos) const
+ {
+ PrimaryVariables values(0.0);
+ values = initialAtPos(pos);
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Neumann control volume.
+ *
+ * \param element The element for which the Neumann boundary condition is set
+ * \param fvGeomentry The fvGeometry
+ * \param elemVolVars The element volume variables
+ * \param elemFaceVars The element face variables
+ * \param scvf The boundary sub control volume face
+ */
+ NumEqVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFaceVariables& elemFaceVars,
+ const SubControlVolumeFace& scvf) const
+ {
+ PrimaryVariables values(0.0);
+
+ if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
+ {
+ values[Indices::momentumYBalanceIdx] = couplingManager().couplingData().momentumCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf);
+
+ const auto massFlux = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf);
+ values[Indices::conti0EqIdx] = massFlux[0];
+ values[Indices::conti0EqIdx + 1] = massFlux[1];
+ }
+ return values;
+ }
+
+ // \}
+
+ /*!
+ * \brief Set the coupling manager
+ */
+ void setCouplingManager(std::shared_ptr<CouplingManager> cm)
+ { couplingManager_ = cm; }
+
+ /*!
+ * \brief Get the coupling manager
+ */
+ const CouplingManager& couplingManager() const
+ { return *couplingManager_; }
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * \param globalPos The global position
+ */
+ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+ {
+ FluidState fluidState;
+ updateFluidStateForBC_(fluidState, pressure_);
+ const Scalar density = FluidSystem::density(fluidState, phaseIdx);
+
+ PrimaryVariables values(0.0);
+ values[Indices::pressureIdx] = pressure_ + density*this->gravity()[1]*(globalPos[1] - this->fvGridGeometry().bBoxMin()[1]);
+ values[transportCompIdx + dim] = moleFraction_;
+ values[Indices::velocityXIdx] = 4.0 * velocity_ * (globalPos[1] - this->fvGridGeometry().bBoxMin()[1])
+ * (this->fvGridGeometry().bBoxMax()[1] - globalPos[1])
+ / (height_() * height_());
+
+ return values;
+ }
+
+ void setTimeLoop(TimeLoopPtr timeLoop)
+ { timeLoop_ = timeLoop; }
+
+ /*!
+ * \brief Returns the intrinsic permeability of required as input parameter for the Beavers-Joseph-Saffman boundary condition
+ */
+ Scalar permeability(const SubControlVolumeFace& scvf) const
+ {
+ return couplingManager().couplingData().darcyPermeability(scvf);
+ }
+
+ /*!
+ * \brief Returns the alpha value required as input parameter for the Beavers-Joseph-Saffman boundary condition
+ */
+ Scalar alphaBJ(const SubControlVolumeFace& scvf) const
+ {
+ return couplingManager().problem(CouplingManager::darcyIdx).spatialParams().beaversJosephCoeffAtPos(scvf.center());
+ }
+
+ // \}
+
+private:
+ bool onLeftBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_; }
+
+ bool onRightBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_; }
+
+ bool onLowerBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] < this->fvGridGeometry().bBoxMin()[1] + eps_; }
+
+ bool onUpperBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_; }
+
+ //! \brief updates the fluid state to obtain required quantities for IC/BC
+ void updateFluidStateForBC_(FluidState& fluidState, const Scalar pressure) const
+ {
+ fluidState.setTemperature(temperature());
+ fluidState.setPressure(phaseIdx, pressure);
+ fluidState.setSaturation(phaseIdx, 1.0);
+ fluidState.setMoleFraction(phaseIdx, transportCompIdx, moleFraction_);
+ fluidState.setMoleFraction(phaseIdx, phaseIdx, 1.0 - moleFraction_);
+
+ typename FluidSystem::ParameterCache paramCache;
+ paramCache.updatePhase(fluidState, phaseIdx);
+
+ const Scalar density = FluidSystem::density(fluidState, paramCache, phaseIdx);
+ fluidState.setDensity(phaseIdx, density);
+
+ const Scalar molarDensity = FluidSystem::molarDensity(fluidState, paramCache, phaseIdx);
+ fluidState.setMolarDensity(phaseIdx, molarDensity);
+
+ const Scalar enthalpy = FluidSystem::enthalpy(fluidState, paramCache, phaseIdx);
+ fluidState.setEnthalpy(phaseIdx, enthalpy);
+ }
+ // the height of the free-flow domain
+ const Scalar height_() const
+ { return this->fvGridGeometry().bBoxMax()[1] - this->fvGridGeometry().bBoxMin()[1]; }
+
+ Scalar eps_;
+
+ Scalar velocity_;
+ Scalar pressure_;
+ Scalar moleFraction_;
+
+ TimeLoopPtr timeLoop_;
+
+ std::shared_ptr<CouplingManager> couplingManager_;
+};
+} //end namespace
+
+#endif // DUMUX_STOKES1P2C_SUBPROBLEM_HH
diff --git a/exercises/solution/exercise-coupling-ff-pm/models/ex_models_pmproblem.hh b/exercises/solution/exercise-coupling-ff-pm/models/ex_models_pmproblem.hh
new file mode 100644
index 00000000..110997ef
--- /dev/null
+++ b/exercises/solution/exercise-coupling-ff-pm/models/ex_models_pmproblem.hh
@@ -0,0 +1,467 @@
+// -*- 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 simple Darcy test problem (cell-centered finite volume method).
+ */
+#ifndef DUMUX_DARCY_SUBPROBLEM_HH
+#define DUMUX_DARCY_SUBPROBLEM_HH
+
+#include <dune/grid/yaspgrid.hh>
+
+#include <dumux/discretization/cellcentered/tpfa/properties.hh>
+#include <dumux/io/gnuplotinterface.hh>
+#include <dumux/material/fluidsystems/h2oair.hh>
+#include <dumux/material/fluidmatrixinteractions/diffusivityconstanttortuosity.hh>
+
+#include <dumux/porousmediumflow/problem.hh>
+
+#if EXNUMBER >= 1
+#include <dumux/porousmediumflow/2pnc/model.hh>
+#include "../2pspatialparams.hh"
+#else
+#include <dumux/porousmediumflow/1pnc/model.hh>
+#include "../1pspatialparams.hh"
+#endif
+
+namespace Dumux
+{
+template <class TypeTag>
+class DarcySubProblem;
+
+namespace Properties
+{
+#if EXNUMBER >= 1
+NEW_TYPE_TAG(DarcyTypeTag, INHERITS_FROM(CCTpfaModel, TwoPNC));
+#else
+NEW_TYPE_TAG(DarcyTypeTag, INHERITS_FROM(CCTpfaModel, OnePNC));
+#endif
+
+// Set the problem property
+SET_TYPE_PROP(DarcyTypeTag, Problem, Dumux::DarcySubProblem<TypeTag>);
+
+// The fluid system
+SET_TYPE_PROP(DarcyTypeTag, FluidSystem, FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>);
+
+// Use gas as phase
+#if EXNUMBER == 0
+SET_INT_PROP(DarcyTypeTag, PhaseIdx, GET_PROP_TYPE(TypeTag, FluidSystem)::gasPhaseIdx);
+#endif
+
+// Use moles
+SET_BOOL_PROP(DarcyTypeTag, UseMoles, true);
+
+// Do not replace one equation with a total mass balance
+SET_INT_PROP(DarcyTypeTag, ReplaceCompEqIdx, 3);
+
+//! Use a model with constant tortuosity for the effective diffusivity
+SET_TYPE_PROP(DarcyTypeTag, EffectiveDiffusivityModel,
+ DiffusivityConstantTortuosity<typename GET_PROP_TYPE(TypeTag, Scalar)>);
+
+// Set the grid type
+SET_TYPE_PROP(DarcyTypeTag, Grid, Dune::YaspGrid<2>);
+
+#if EXNUMBER >= 1
+//! Set the default formulation to pw-Sn: This can be over written in the problem.
+SET_PROP(DarcyTypeTag, Formulation)
+{ static constexpr auto value = TwoPFormulation::p1s0; };
+#endif
+
+// Set the spatial paramaters type
+#if EXNUMBER >= 1
+SET_TYPE_PROP(DarcyTypeTag, SpatialParams, TwoPSpatialParams<TypeTag>);
+#else
+SET_TYPE_PROP(DarcyTypeTag, SpatialParams, OnePSpatialParams<TypeTag>);
+#endif
+}
+
+template <class TypeTag>
+class DarcySubProblem : public PorousMediumFlowProblem<TypeTag>
+{
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+
+ // copy some indices for convenience
+ using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+ enum {
+ // grid and world dimension
+ dim = GridView::dimension,
+ dimworld = GridView::dimensionworld,
+
+ // primary variable indices
+ conti0EqIdx = Indices::conti0EqIdx,
+ pressureIdx = Indices::pressureIdx,
+#if EXNUMBER >= 3
+ saturationIdx = Indices::switchIdx,
+ transportCompIdx = Indices::switchIdx
+#elif EXNUMBER >= 1
+ transportCompIdx = Indices::switchIdx
+#else
+ phaseIdx = Indices::fluidSystemPhaseIdx,
+ transportCompIdx = 1-phaseIdx
+#endif
+ };
+
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
+
+ using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager);
+ using TimeLoopPtr = std::shared_ptr<TimeLoop<Scalar>>;
+
+public:
+ DarcySubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry,
+ std::shared_ptr<CouplingManager> couplingManager)
+ : ParentType(fvGridGeometry, "Darcy"), eps_(1e-7), couplingManager_(couplingManager)
+ {
+#if EXNUMBER >= 3
+ saturation_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Saturation");
+#else
+ moleFraction_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.MoleFraction");
+#endif
+ pressure_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Pressure");
+
+ // initialize output file
+ plotFluxes_ = getParamFromGroup<bool>(this->paramGroup(), "Problem.PlotFluxes", false);
+ plotStorage_ = getParamFromGroup<bool>(this->paramGroup(), "Problem.PlotStorage", false);
+ storageFileName_ = "storage_" + getParam<std::string>("Problem.Name") + "_" + this->name() + ".csv";
+ storageFile_.open(storageFileName_);
+ storageFile_ << "#Time[s]" << ";"
+ << "WaterMass[kg]" << ";"
+ << "WaterMassLoss[kg]" << ";"
+ << "EvaporationRate[mm/d]"
+ << std::endl;
+ }
+
+ /*!
+ * \name Simulation steering
+ */
+ // \{
+
+ /*!
+ * \brief Initialize the problem.
+ */
+ template<class SolutionVector, class GridVariables>
+ void init(const SolutionVector& curSol,
+ const GridVariables& gridVariables)
+ {
+#if EXNUMBER >= 2
+ initialWaterContent_ = evaluateWaterMassStorageTerm(curSol, gridVariables);
+ lastWaterMass_ = initialWaterContent_;
+#endif
+ }
+
+ template<class SolutionVector, class GridVariables>
+ void postTimeStep(const SolutionVector& curSol,
+ const GridVariables& gridVariables)
+
+ {
+ evaluateWaterMassStorageTerm(curSol, gridVariables);
+ evaluateInterfaceFluxes(curSol, gridVariables);
+
+ gnuplotStorage_.resetPlot();
+ gnuplotStorage_.setDatafileSeparator(';');
+ gnuplotStorage_.setXlabel("time [d]");
+ gnuplotStorage_.setXRange(0.0, getParam<Scalar>("TimeLoop.TEnd"));
+ gnuplotStorage_.setYlabel("evaporation rate [mm/d]");
+ gnuplotStorage_.setOption("set yrange [0.0:]");
+ gnuplotStorage_.setOption("set y2label 'cumulative mass loss'");
+ gnuplotStorage_.setOption("set y2range [0.0:0.5]");
+ gnuplotStorage_.setOption("set y2range [0.0:0.5]");
+ gnuplotStorage_.addFileToPlot(storageFileName_, "using 1:4 with lines title 'evaporation rate'");
+ gnuplotStorage_.addFileToPlot(storageFileName_, "using 1:3 axes x1y2 with lines title 'cumulative mass loss'");
+ if (plotStorage_)
+ gnuplotStorage_.plot("temp");
+ }
+
+ template<class SolutionVector, class GridVariables>
+ Scalar evaluateWaterMassStorageTerm(const SolutionVector& curSol,
+ const GridVariables& gridVariables)
+
+ {
+ // compute the mass in the entire domain
+ Scalar waterMass = 0.0;
+
+ for (const auto& element : elements(this->fvGridGeometry().gridView()))
+ {
+ auto fvGeometry = localView(this->fvGridGeometry());
+ fvGeometry.bindElement(element);
+
+ auto elemVolVars = localView(gridVariables.curGridVolVars());
+ elemVolVars.bindElement(element, fvGeometry, curSol);
+
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ const auto& volVars = elemVolVars[scv];
+ for(int phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
+ {
+ // insert calculation of the water mass here
+#if EXNUMBER >= 2
+ waterMass += volVars.massFraction(phaseIdx, FluidSystem::H2OIdx) * volVars.density(phaseIdx)
+ * volVars.saturation(phaseIdx) * volVars.porosity()
+ * scv.volume() * volVars.extrusionFactor();
+#else
+ waterMass += 0.0;
+#endif
+ }
+ }
+ }
+#if EXNUMBER >= 2
+ std::cout << "Mass of water is: " << waterMass << std::endl;
+#endif
+
+ Scalar cumMassLoss = initialWaterContent_ - waterMass;
+ Scalar evaporationRate = (lastWaterMass_ - waterMass) * 86400
+ / (this->fvGridGeometry().bBoxMax()[0] - this->fvGridGeometry().bBoxMin()[0])
+ / timeLoop_->timeStepSize();
+ lastWaterMass_ = waterMass;
+
+ storageFile_ << timeLoop_->time() << ";"
+ << waterMass << ";"
+ << cumMassLoss << ";"
+ << evaporationRate
+ << std::endl;
+
+ return waterMass;
+ }
+
+ template<class SolutionVector, class GridVariables>
+ void evaluateInterfaceFluxes(const SolutionVector& curSol,
+ const GridVariables& gridVariables)
+
+ {
+ std::vector<Scalar> x;
+ std::vector<Scalar> y;
+
+ for (const auto& element : elements(this->fvGridGeometry().gridView()))
+ {
+ auto fvGeometry = localView(this->fvGridGeometry());
+ fvGeometry.bindElement(element);
+
+ auto elemVolVars = localView(gridVariables.curGridVolVars());
+ elemVolVars.bindElement(element, fvGeometry, curSol);
+
+ for (auto&& scvf : scvfs(fvGeometry))
+ {
+ if (!couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
+ continue;
+
+ // NOTE: binding the coupling context is necessary
+ couplingManager_->bindCouplingContext(CouplingManager::darcyIdx, element);
+#if EXNUMBER >= 2
+ NumEqVector flux = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, scvf);
+#else
+ NumEqVector flux(0.0); // add "massCouplingCondition" from the couplingManager here
+#endif
+
+ x.push_back(scvf.center()[0]);
+ y.push_back(flux[transportCompIdx]);
+ }
+ }
+
+ gnuplotInterfaceFluxes_.resetPlot();
+ gnuplotInterfaceFluxes_.setXlabel("x-position [m]");
+ gnuplotInterfaceFluxes_.setXRange(this->fvGridGeometry().bBoxMin()[0], this->fvGridGeometry().bBoxMax()[0]);
+ gnuplotInterfaceFluxes_.setYlabel("flux [kg/(m^2 s)]");
+ gnuplotInterfaceFluxes_.setYRange(-5e-4, 0.0);
+ gnuplotInterfaceFluxes_.setOption("set label 'time: " + std::to_string(timeLoop_->time()/86400.) + "d' at graph 0.8,0.8 ");
+ std::string fluxFileName = "flux_" + std::to_string(timeLoop_->timeStepIndex()) +
+ "_" + getParam<std::string>("Problem.Name") + "_" + this->name() + ".csv";
+ gnuplotInterfaceFluxes_.addDataSetToPlot(x, y, fluxFileName, "with lines title 'water mass flux'");
+ if (plotFluxes_)
+ gnuplotInterfaceFluxes_.plot("flux_" + std::to_string(timeLoop_->timeStepIndex()));
+ }
+
+ /*!
+ * \brief Returns true if a restart file should be written to
+ * disk.
+ */
+ bool shouldWriteRestartFile() const
+ { return false; }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ bool shouldWriteOutput() const // define output
+ { return true; }
+
+ /*!
+ * \brief Return the temperature within the domain in [K].
+ *
+ */
+ Scalar temperature() const
+ { return 293.15; }
+ // \}
+
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
+
+ /*!
+ * \brief Specifies which kind of boundary condition should be
+ * used for which equation on a given boundary control volume.
+ *
+ * \param element The element
+ * \param scvf The boundary sub control volume face
+ */
+ BoundaryTypes boundaryTypes(const Element& element, const SubControlVolumeFace& scvf) const
+ {
+ BoundaryTypes values;
+ values.setAllNeumann();
+
+ if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
+ values.setAllCouplingNeumann();
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Neumann control volume.
+ *
+ * \param element The element for which the Neumann boundary condition is set
+ * \param fvGeomentry The fvGeometry
+ * \param elemVolVars The element volume variables
+ * \param scvf The boundary sub control volume face
+ *
+ * For this method, the \a values variable stores primary variables.
+ */
+ template<class ElementVolumeVariables>
+ NumEqVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf) const
+ {
+ NumEqVector values(0.0);
+
+ if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
+ values = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, scvf);
+
+ return values;
+ }
+
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+ /*!
+ * \brief Evaluate the source term for all phases within a given
+ * sub-control-volume.
+ *
+ * \param element The element for which the source term is set
+ * \param fvGeomentry The fvGeometry
+ * \param elemVolVars The element volume variables
+ * \param scv The subcontrolvolume
+ */
+ template<class ElementVolumeVariables>
+ NumEqVector source(const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolume &scv) const
+ { return NumEqVector(0.0); }
+
+ // \}
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * \param element The element
+ *
+ * For this method, the \a priVars parameter stores primary
+ * variables.
+ */
+ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+ {
+ PrimaryVariables values(0.0);
+#if EXNUMBER >= 3
+ values.setState(3/*bothPhases*/);
+ values[saturationIdx] = saturation_;
+#elif EXNUMBER >= 1
+ values.setState(2/*secondPhaseOnly*/);
+ values[transportCompIdx] = moleFraction_;
+#else
+ values[transportCompIdx] = moleFraction_;
+#endif
+ values[pressureIdx] = pressure_;
+ return values;
+ }
+
+ // \}
+
+ //! Set the coupling manager
+ void setCouplingManager(std::shared_ptr<CouplingManager> cm)
+ { couplingManager_ = cm; }
+
+ //! Get the coupling manager
+ const CouplingManager& couplingManager() const
+ { return *couplingManager_; }
+
+ void setTimeLoop(TimeLoopPtr timeLoop)
+ { timeLoop_ = timeLoop; }
+
+private:
+ bool onLeftBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_; }
+
+ bool onRightBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_; }
+
+ bool onLowerBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] < this->fvGridGeometry().bBoxMin()[1] + eps_; }
+
+ bool onUpperBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_; }
+
+ Scalar eps_;
+#if EXNUMBER >= 3
+ Scalar saturation_;
+#else
+ Scalar moleFraction_;
+#endif
+ Scalar pressure_;
+
+ Scalar initialWaterContent_ = 0.0;
+ Scalar lastWaterMass_ = 0.0;
+
+ TimeLoopPtr timeLoop_;
+ std::shared_ptr<CouplingManager> couplingManager_;
+
+ std::string storageFileName_;
+ std::ofstream storageFile_;
+ bool plotFluxes_;
+ bool plotStorage_;
+ Dumux::GnuplotInterface<Scalar> gnuplotInterfaceFluxes_;
+ Dumux::GnuplotInterface<Scalar> gnuplotStorage_;
+};
+} //end namespace
+
+#endif //DUMUX_DARCY_SUBPROBLEM_HH
--
GitLab