From f6a82454c75a9334de9e931ccd630f5216ce36f9 Mon Sep 17 00:00:00 2001
From: Alexander Jaust <alexander.jaust@ipvs.uni-stuttgart.de>
Date: Fri, 10 May 2019 13:24:09 +0200
Subject: [PATCH] added free flow + porous media flow coupling case from
tutorial; The case should be at the state 'Task B: Include slip-condition';
Can we verify the case and split it into two parts?
---
appl/CMakeLists.txt | 1 +
appl/conjugateheattransfer/CMakeLists.txt | 1 +
appl/coupling-ff-pm/CMakeLists.txt | 5 +
.../monolithic/1pspatialparams.hh | 91 ++++
.../monolithic/2pspatialparams.hh | 132 +++++
appl/coupling-ff-pm/monolithic/CMakeLists.txt | 3 +
appl/coupling-ff-pm/monolithic/README.md | 403 +++++++++++++++
.../monolithic/interface/CMakeLists.txt | 9 +
.../interface/ex_interface_coupling_ff-pm.cc | 280 +++++++++++
.../ex_interface_coupling_ff-pm.input | 46 ++
.../interface/ex_interface_ffproblem.hh | 347 +++++++++++++
.../interface/ex_interface_pmproblem.hh | 271 ++++++++++
.../monolithic/models/CMakeLists.txt | 9 +
.../models/ex_models_coupling_ff-pm.cc | 283 +++++++++++
.../models/ex_models_coupling_ff-pm.input | 54 ++
.../monolithic/models/ex_models_ffproblem.hh | 354 +++++++++++++
.../monolithic/models/ex_models_pmproblem.hh | 404 +++++++++++++++
.../monolithic/turbulence/CMakeLists.txt | 10 +
.../ex_turbulence_coupling_ff-pm.cc | 284 +++++++++++
.../ex_turbulence_coupling_ff-pm.input | 68 +++
.../turbulence/ex_turbulence_ffproblem.hh | 387 ++++++++++++++
.../turbulence/ex_turbulence_pmproblem.hh | 337 +++++++++++++
appl/coupling-ff-pm/solution/CMakeLists.txt | 1 +
.../solution/monolithic/1pspatialparams.hh | 91 ++++
.../solution/monolithic/2pspatialparams.hh | 132 +++++
.../solution/monolithic/CMakeLists.txt | 3 +
.../monolithic/interface/CMakeLists.txt | 22 +
.../interface/ex_interface_coupling_ff-pm.cc | 274 ++++++++++
.../interface/ex_interface_ffproblem.hh | 367 ++++++++++++++
.../interface/ex_interface_pmproblem.hh | 266 ++++++++++
.../monolithic/interface/params.input | 46 ++
.../solution/monolithic/models/CMakeLists.txt | 21 +
.../models/ex_models_coupling_ff-pm.cc | 289 +++++++++++
.../monolithic/models/ex_models_ffproblem.hh | 357 +++++++++++++
.../monolithic/models/ex_models_pmproblem.hh | 475 ++++++++++++++++++
.../solution/monolithic/models/params.input | 50 ++
.../monolithic/turbulence/CMakeLists.txt | 17 +
.../ex_turbulence_coupling_ff-pm.cc | 290 +++++++++++
.../turbulence/ex_turbulence_ffproblem.hh | 423 ++++++++++++++++
.../turbulence/ex_turbulence_pmproblem.hh | 337 +++++++++++++
.../monolithic/turbulence/params.input | 68 +++
41 files changed, 7308 insertions(+)
create mode 100644 appl/coupling-ff-pm/CMakeLists.txt
create mode 100644 appl/coupling-ff-pm/monolithic/1pspatialparams.hh
create mode 100644 appl/coupling-ff-pm/monolithic/2pspatialparams.hh
create mode 100644 appl/coupling-ff-pm/monolithic/CMakeLists.txt
create mode 100644 appl/coupling-ff-pm/monolithic/README.md
create mode 100644 appl/coupling-ff-pm/monolithic/interface/CMakeLists.txt
create mode 100644 appl/coupling-ff-pm/monolithic/interface/ex_interface_coupling_ff-pm.cc
create mode 100644 appl/coupling-ff-pm/monolithic/interface/ex_interface_coupling_ff-pm.input
create mode 100644 appl/coupling-ff-pm/monolithic/interface/ex_interface_ffproblem.hh
create mode 100644 appl/coupling-ff-pm/monolithic/interface/ex_interface_pmproblem.hh
create mode 100644 appl/coupling-ff-pm/monolithic/models/CMakeLists.txt
create mode 100644 appl/coupling-ff-pm/monolithic/models/ex_models_coupling_ff-pm.cc
create mode 100644 appl/coupling-ff-pm/monolithic/models/ex_models_coupling_ff-pm.input
create mode 100644 appl/coupling-ff-pm/monolithic/models/ex_models_ffproblem.hh
create mode 100644 appl/coupling-ff-pm/monolithic/models/ex_models_pmproblem.hh
create mode 100644 appl/coupling-ff-pm/monolithic/turbulence/CMakeLists.txt
create mode 100644 appl/coupling-ff-pm/monolithic/turbulence/ex_turbulence_coupling_ff-pm.cc
create mode 100644 appl/coupling-ff-pm/monolithic/turbulence/ex_turbulence_coupling_ff-pm.input
create mode 100644 appl/coupling-ff-pm/monolithic/turbulence/ex_turbulence_ffproblem.hh
create mode 100644 appl/coupling-ff-pm/monolithic/turbulence/ex_turbulence_pmproblem.hh
create mode 100644 appl/coupling-ff-pm/solution/CMakeLists.txt
create mode 100644 appl/coupling-ff-pm/solution/monolithic/1pspatialparams.hh
create mode 100644 appl/coupling-ff-pm/solution/monolithic/2pspatialparams.hh
create mode 100644 appl/coupling-ff-pm/solution/monolithic/CMakeLists.txt
create mode 100644 appl/coupling-ff-pm/solution/monolithic/interface/CMakeLists.txt
create mode 100644 appl/coupling-ff-pm/solution/monolithic/interface/ex_interface_coupling_ff-pm.cc
create mode 100644 appl/coupling-ff-pm/solution/monolithic/interface/ex_interface_ffproblem.hh
create mode 100644 appl/coupling-ff-pm/solution/monolithic/interface/ex_interface_pmproblem.hh
create mode 100644 appl/coupling-ff-pm/solution/monolithic/interface/params.input
create mode 100644 appl/coupling-ff-pm/solution/monolithic/models/CMakeLists.txt
create mode 100644 appl/coupling-ff-pm/solution/monolithic/models/ex_models_coupling_ff-pm.cc
create mode 100644 appl/coupling-ff-pm/solution/monolithic/models/ex_models_ffproblem.hh
create mode 100644 appl/coupling-ff-pm/solution/monolithic/models/ex_models_pmproblem.hh
create mode 100644 appl/coupling-ff-pm/solution/monolithic/models/params.input
create mode 100644 appl/coupling-ff-pm/solution/monolithic/turbulence/CMakeLists.txt
create mode 100644 appl/coupling-ff-pm/solution/monolithic/turbulence/ex_turbulence_coupling_ff-pm.cc
create mode 100644 appl/coupling-ff-pm/solution/monolithic/turbulence/ex_turbulence_ffproblem.hh
create mode 100644 appl/coupling-ff-pm/solution/monolithic/turbulence/ex_turbulence_pmproblem.hh
create mode 100644 appl/coupling-ff-pm/solution/monolithic/turbulence/params.input
diff --git a/appl/CMakeLists.txt b/appl/CMakeLists.txt
index c85eea6..933871d 100644
--- a/appl/CMakeLists.txt
+++ b/appl/CMakeLists.txt
@@ -1 +1,2 @@
add_subdirectory(conjugateheattransfer)
+add_subdirectory(coupling-ff-pm)
diff --git a/appl/conjugateheattransfer/CMakeLists.txt b/appl/conjugateheattransfer/CMakeLists.txt
index 25a2da4..84fbd12 100644
--- a/appl/conjugateheattransfer/CMakeLists.txt
+++ b/appl/conjugateheattransfer/CMakeLists.txt
@@ -1,2 +1,3 @@
add_subdirectory(monolithic)
+add_subdirectory(monolithic-const-dt)
add_subdirectory(iterative)
diff --git a/appl/coupling-ff-pm/CMakeLists.txt b/appl/coupling-ff-pm/CMakeLists.txt
new file mode 100644
index 0000000..4669328
--- /dev/null
+++ b/appl/coupling-ff-pm/CMakeLists.txt
@@ -0,0 +1,5 @@
+# add a target that builds all exercise solutions
+add_custom_target(test_exercises)
+
+add_subdirectory(monolithic)
+add_subdirectory(solution)
diff --git a/appl/coupling-ff-pm/monolithic/1pspatialparams.hh b/appl/coupling-ff-pm/monolithic/1pspatialparams.hh
new file mode 100644
index 0000000..a042388
--- /dev/null
+++ b/appl/coupling-ff-pm/monolithic/1pspatialparams.hh
@@ -0,0 +1,91 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup OnePTests
+ * \brief The spatial parameters class for the test problem using the 1p cc model
+ */
+#ifndef DUMUX_1P_TEST_SPATIALPARAMS_HH
+#define DUMUX_1P_TEST_SPATIALPARAMS_HH
+
+#include <dumux/material/spatialparams/fv1p.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup OnePModel
+ *
+ * \brief The spatial parameters class for the test problem using the
+ * 1p cc model
+ */
+template<class FVGridGeometry, class Scalar>
+class OnePSpatialParams
+: public FVSpatialParamsOneP<FVGridGeometry, Scalar, OnePSpatialParams<FVGridGeometry, Scalar>>
+{
+ using GridView = typename FVGridGeometry::GridView;
+ using ParentType = FVSpatialParamsOneP<FVGridGeometry, Scalar, OnePSpatialParams<FVGridGeometry, Scalar>>;
+
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+public:
+ // export permeability type
+ using PermeabilityType = Scalar;
+
+ OnePSpatialParams(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+ : ParentType(fvGridGeometry)
+ {
+ permeability_ = getParam<Scalar>("Darcy.SpatialParams.Permeability");
+ porosity_ = getParam<Scalar>("Darcy.SpatialParams.Porosity");
+ alphaBJ_ = getParam<Scalar>("Darcy.SpatialParams.AlphaBeaversJoseph");
+ }
+
+ /*!
+ * \brief Function for defining the (intrinsic) permeability \f$[m^2]\f$.
+ *
+ * \param globalPos The global position
+ * \return the intrinsic permeability
+ */
+ PermeabilityType permeabilityAtPos(const GlobalPosition& globalPos) const
+ { return permeability_; }
+
+ /*! \brief Define the porosity in [-].
+ *
+ * \param globalPos The global position
+ */
+ Scalar porosityAtPos(const GlobalPosition& globalPos) const
+ { return porosity_; }
+
+ /*! \brief Define the Beavers-Joseph coefficient in [-].
+ *
+ * \param globalPos The global position
+ */
+ Scalar beaversJosephCoeffAtPos(const GlobalPosition& globalPos) const
+ { return alphaBJ_; }
+
+
+private:
+ Scalar permeability_;
+ Scalar porosity_;
+ Scalar alphaBJ_;
+};
+
+} // end namespace
+
+#endif
diff --git a/appl/coupling-ff-pm/monolithic/2pspatialparams.hh b/appl/coupling-ff-pm/monolithic/2pspatialparams.hh
new file mode 100644
index 0000000..622e4f5
--- /dev/null
+++ b/appl/coupling-ff-pm/monolithic/2pspatialparams.hh
@@ -0,0 +1,132 @@
+// -*- 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
+ *
+ * \brief The spatial parameters class for the test problem using the 2p cc model
+ */
+template<class FVGridGeometry, class Scalar>
+class TwoPSpatialParams
+: public FVSpatialParams<FVGridGeometry, Scalar, TwoPSpatialParams<FVGridGeometry, Scalar>>
+{
+ using GridView = typename FVGridGeometry::GridView;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using ParentType = FVSpatialParams<FVGridGeometry, Scalar, TwoPSpatialParams<FVGridGeometry, Scalar>>;
+
+ 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 Dumux
+
+#endif
diff --git a/appl/coupling-ff-pm/monolithic/CMakeLists.txt b/appl/coupling-ff-pm/monolithic/CMakeLists.txt
new file mode 100644
index 0000000..bae4c26
--- /dev/null
+++ b/appl/coupling-ff-pm/monolithic/CMakeLists.txt
@@ -0,0 +1,3 @@
+add_subdirectory(interface)
+add_subdirectory(models)
+add_subdirectory(turbulence)
diff --git a/appl/coupling-ff-pm/monolithic/README.md b/appl/coupling-ff-pm/monolithic/README.md
new file mode 100644
index 0000000..cbd87c5
--- /dev/null
+++ b/appl/coupling-ff-pm/monolithic/README.md
@@ -0,0 +1,403 @@
+# Exercise Coupling free flow/porous medium flow
+
+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.
+
+
+### 0. Getting familiar with the code
+
+* Navigate to the directory `exercises/exercise-coupling-ff-pm`
+
+There are three sub folders: `interface` (Exercise 1), `models` (Exercise 2) and `turbulence` (Exercise 3).
+
+The problem-related files for this exercise are:
+* 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`.
+
+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
+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
+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.
+Keep in mind that the staggered grid implementation distinguishes between face variables (velocity components) and
+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`.
+
+__Task__:
+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.
+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.
+
+__Task A: Change the flow direction__
+
+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.
+
+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);
+
+```
+
+Set a Dirichlet boundary condition for the velocities at the top:
+``` cpp
+if(onUpperBoundary_(globalPos))
+{
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+}
+```
+
+Keep the coupling boundary condition:
+``` cpp
+if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
+{
+ values.setCouplingNeumann(Indices::conti0EqIdx);
+ values.setCouplingNeumann(Indices::momentumYBalanceIdx);
+ values.setDirichlet(Indices::velocityXIdx); // assume no slip on interface
+}
+```
+
+Having changed the types of boundary conditions, we must now assign the correct values for them.
+
+Set a no-slip, no-flow condition for the velocity at the top:
+``` cpp
+values[Indices::velocityXIdx] = 0.0;
+values[Indices::velocityYIdx] = 0.0;
+```
+Apply a fixed pressure difference between the inlet and outlet, e.g.:
+``` cpp
+if(onLeftBoundary_(globalPos))
+ values[Indices::pressureIdx] = deltaP_;
+if(onRightBoundary_(globalPos))
+ values[Indices::pressureIdx] = 0.0;
+```
+
+For changing the flow direction, the boundary conditions for the porous medium have to be changed as well.
+
+Use Neumann no-flow boundaries everywhere, keep the coupling conditions.
+``` cpp
+values.setAllNeumann();
+
+if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
+ values.setAllCouplingNeumann();
+```
+
+This should make the flow go from left to right.
+
+__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`$.
+
+To include this, just replace the no-slip condition at the interface
+``` cpp
+values.setDirichlet(Indices::velocityXIdx); // assume no slip on interface
+```
+with a Beavers-Joseph-Saffman (BJS) boundary condition for the respective momentum balance:
+``` cpp
+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 in `ex_interface_coupling_ff-pm.cc`.
+```cpp
+stokesVtkWriter.addField(stokesProblem->getAnalyticalVelocityX(), "analyticalV_x");
+```
+
+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.
+
+__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 96-106 in `ex_interface_coupling_ff-pm.cc` and comment lines 114-149 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,
+the domain is split in two haves, separated by a sinusoidal interface.
+
+```cpp
+auto elementSelectorStokes = [&](const auto& element)
+{
+ double interface = params.amplitude * std::sin(( element.geometry().center()[0] -params.offset) / params.scaling * 2.0 * M_PI) + params.baseline;
+ return element.geometry().center()[1] > interface;
+};
+
+auto elementSelectorDarcy = [&](const auto& element)
+{
+ double interface = params.amplitude * std::sin(( element.geometry().center()[0] - params.offset) / params.scaling * 2.0 * M_PI) + params.baseline;
+ return element.geometry().center()[1] < interface;
+};
+```
+
+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 now compile. 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}`$
+
+ was always set for a fixed $`\mathbf{n} = (0,1)^T`$. We need to account for the curvature of the interface and thus replace
+ ```cpp
+values.setCouplingNeumann(Indices::momentumYBalanceIdx);
+ ```
+ with
+ ```cpp
+values.setCouplingNeumann(scvf.directionIndex());
+ ```
+The same if true for the BJS condition, however, here we need to consider the tangential direction:
+```cpp
+values.setBJS(1 - scvf.directionIndex());
+```
+
+The final result should look something like this:
+
+
+
+
+### 2. Changing the porous medium model
+
+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).
+
+This part of the exercise consists of the following steps:
+* replacing the 1pnc porous-medium model by a 2pnc porous-medium model,
+* adding the output for fluxes and storage terms (having gnuplot is recommended),
+* enable a variable switch for a complete drying.
+
+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.
+(Don't mind the compiler warning, we will deal with it in task B.)
+
+
+__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, i.e., air with water vapor in both domains.
+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).
+* Since two phases are involved now, we do not need to use the OnePAdapter anymore. Change to property of the FluidSystem such that `H2OAir` is used directly.
+ Afterwards, set the `transportCompIdx` to `Indices::switchIdx`.
+
+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
+
+```
+template<class TypeTag>
+struct Formulation<TypeTag, TTag::DarcyOnePNC>
+{ 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 phase presence states begins
+ with 1. Take a look at your formulation to find out which phase index to use for the gas phase.)
+
+__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, adapt the method `init()` such that the variable `initialWaterContent_` is initialized correctly using the `evaluateWaterMassStorageTerm()` method and assign that value also to the variable `lastWaterMass_`.
+
+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 too fast, you can have a look at the flux*.png files after the simulation.
+* You can use the parameter `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. (2002)](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.
+
+Now you are able to simulate a complete drying of the porous medium.
+
+
+### 4. Use a turbulence model in the free flow domain
+
+Several RANS turbulence models are implemented in DuMuX.
+This part of the exercise consists of the following steps:
+* replacing the Navier-Stokes model by the zero equation turbulence model,
+* switching to a symmetry boundary condition,
+* applying a grid refinement towards the interface,
+* subsequently refining the grid (convergence study).
+
+We will work with a `1p2cni/2p2cni` coupled problem, where `ni` stands for non-isothermal.
+All the prepared files can be found in the subfolder `exercise-coupling-ff-pm/turbulence`.
+
+__Task A__:
+
+The file `ex_turbulence_ffproblem.hh` is your free flow problem file within this exercise.
+
+For using the compositional zero equation turbulence model, the following header files need to be included:
+```
+#include <dumux/freeflow/compositional/zeroeqncmodel.hh>
+#include <dumux/freeflow/rans/zeroeq/problem.hh>
+```
+The includes for the NavierStokesNC model and the NavierStokesProblem are no longer needed and can be removed.
+
+Make sure your free flow problem inherits from the correct parent type:
+* Change the entry in the `StokesZeroEq` definition accordingly (non-isothermal zero equation model)
+* Adapt the inheritance of the problem class (hint: two occurrences)
+
+Take a look into the two headers included above to see how the correct TypeTag and the Problem class the inherit from are called.
+
+Here, the turbulent free flow is wall bounded which means that the main reason for the development
+of turbulent flow is the presence of walls.
+The porous medium at the bottom and also the channel wall are such walls.
+Inside the problem you have to tell the turbulence model, where these walls are located
+by providing an `isOnWallAtPos()` function:
+```c++
+bool isOnWallAtPos(const GlobalPosition& globalPos) const
+{
+ return (onLowerBoundary_(globalPos) || onUpperBoundary_(globalPos));
+}
+```
+
+In addition, especially for the zero-equation models, any element in the free-flow domain interacts with the walls,
+e.g. this defines the wall distance which is needed to calculate the eddy viscosity.
+To get all these interactions, you have to call
+ ```cpp
+ stokesProblem->updateStaticWallProperties()
+ ```
+in `ex_turbulence_coupling_ff-pm.cc`.
+However, there is also a solution-dependent component of these interactions, e.g. for a correct
+damping of the eddy viscosity toward the wall, the velocity gradient at the wall and inside the
+cells is needed.
+These dynamic interactions are to be updated by calling
+```cpp
+stokesProblem->updateDynamicWallProperties(stokesSol)
+```
+in the time loop (after `// update dynamic wall properties`).
+
+Compile and run your new coupled problem and take a look at the results in Paraview.
+In addition to the standard variables and parameters, you can now analyze turbulence model specific quantities
+(e.g. the turbulent viscosity $`\nu_\textrm{t}`$ or the turbulent diffusivity $`D_\textrm{t}`$) for the free flow domain.
+In paraview you may compare the magnitude of $`D`$ and $`D_\textrm{t}`$ to see where the transport is affected by turbulence.
+The result for the turbulent viscosity should look like this:
+
+
+
+__Task B__:
+
+Instead of computing the whole cross-section of a channel, you can use symmetric boundary conditions at the top boundary of your free flow domain by replacing all previous boundary conditions with
+```c++
+values.setAllSymmetry();
+```
+
+In addition, you have to remove the condition `onUpperBoundary_(globalPos)` from the `isOnWallAtPos(globalPos)`
+and `initialAtPos(globalPos)` method.
+
+__Task C__:
+
+Choose `Surface With Edges` instead of `Surface` in the Paraview toolbar to see the discretization grid.
+We will refine the grid now in order to better resolve the processes at the coupling interface.
+Since not much is happening at the upper and lower boundaries of the whole domain, we want to keep the resolution low in these areas to save some computation time.
+
+A grid refinement towards the interface is called _grading_.
+Try different gradings by changing the values in the respective sections in the input file:
+```c++
+Grading0 = 1.0
+Grading1 = 1.0
+```
+
+__Task D__:
+
+For the grid convergence study, run various simulations with the following grading parameters:
+```c++
+* [Stokes.Grid] Grading1 = 1.2, [Darcy.Grid] Grading1 = -1.2
+* [Stokes.Grid] Grading1 = 1.3, [Darcy.Grid] Grading1 = -1.3
+* [Stokes.Grid] Grading1 = 1.4, [Darcy.Grid] Grading1 = -1.4
+* [Stokes.Grid] Grading1 = 1.5, [Darcy.Grid] Grading1 = -1.5
+* [Stokes.Grid] Grading1 = 1.6, [Darcy.Grid] Grading1 = -1.6
+```
+
+By changing the parameter `Problem.Name` for each grading factor, you avoid loosing the `.vtu` and `.pvd` files of the previous simulation runs.
+Check the first lines of the output to see how the grading factors change the height of your grid cells.
+Compare the velocity fields with Paraview.
diff --git a/appl/coupling-ff-pm/monolithic/interface/CMakeLists.txt b/appl/coupling-ff-pm/monolithic/interface/CMakeLists.txt
new file mode 100644
index 0000000..2865c75
--- /dev/null
+++ b/appl/coupling-ff-pm/monolithic/interface/CMakeLists.txt
@@ -0,0 +1,9 @@
+# executables for ex_interface_coupling_ff-pm
+dune_add_test(NAME ex_interface_coupling_ff-pm
+ SOURCES ex_interface_coupling_ff-pm.cc)
+
+# add tutorial to the common target
+add_dependencies(test_exercises ex_interface_coupling_ff-pm)
+
+# add a symlink for each input file
+add_input_file_links()
diff --git a/appl/coupling-ff-pm/monolithic/interface/ex_interface_coupling_ff-pm.cc b/appl/coupling-ff-pm/monolithic/interface/ex_interface_coupling_ff-pm.cc
new file mode 100644
index 0000000..8c44bfb
--- /dev/null
+++ b/appl/coupling-ff-pm/monolithic/interface/ex_interface_coupling_ff-pm.cc
@@ -0,0 +1,280 @@
+// -*- 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 <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/istl/io.hh>
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/dumuxmessage.hh>
+#include <dumux/common/partial.hh>
+#include <dumux/linear/seqsolverbackend.hh>
+#include <dumux/assembly/fvassembler.hh>
+#include <dumux/assembly/diffmethod.hh>
+#include <dumux/discretization/method.hh>
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/staggeredvtkoutputmodule.hh>
+#include <dumux/io/grid/gridmanager.hh>
+
+#include <dumux/multidomain/staggeredtraits.hh>
+#include <dumux/multidomain/fvassembler.hh>
+#include <dumux/multidomain/newtonsolver.hh>
+
+#include <dumux/multidomain/boundary/stokesdarcy/couplingmanager.hh>
+
+#include "ex_interface_pmproblem.hh"
+#include "ex_interface_ffproblem.hh"
+
+namespace Dumux {
+namespace Properties {
+
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::StokesOneP>
+{
+ using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, Properties::TTag::DarcyOneP>;
+ using type = Dumux::StokesDarcyCouplingManager<Traits>;
+};
+
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::DarcyOneP>
+{
+ using Traits = StaggeredMultiDomainTraits<Properties::TTag::StokesOneP, Properties::TTag::StokesOneP, 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 = Properties::TTag::StokesOneP;
+ using DarcyTypeTag = Properties::TTag::DarcyOneP;
+
+
+
+ // ******************** comment-out this section for the last exercise **************** //
+
+ // create two individual grids (from the given grid file or the input file)
+ // for both sub-domains
+ using DarcyGridManager = Dumux::GridManager<GetPropType<DarcyTypeTag, Properties::Grid>>;
+ DarcyGridManager darcyGridManager;
+ darcyGridManager.init("Darcy"); // pass parameter group
+
+ using StokesGridManager = Dumux::GridManager<GetPropType<StokesTypeTag, Properties::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();
+
+ // ************************************************************************************ //
+
+
+ // ******************** uncomment this section for the last exercise ****************** //
+
+ // // use dune-subgrid to create the individual grids
+ // static constexpr int dim = 2;
+ // using HostGrid = Dune::YaspGrid<2, Dune::TensorProductCoordinates<double, dim> >;
+ // using HostGridManager = Dumux::GridManager<HostGrid>;
+ // HostGridManager hostGridManager;
+ // hostGridManager.init();
+ // auto& hostGrid = hostGridManager.grid();
+ //
+ // struct Params
+ // {
+ // double amplitude = getParam<double>("Grid.Amplitude");
+ // double baseline = getParam<double>("Grid.Baseline");
+ // double offset = getParam<double>("Grid.Offset");
+ // double scaling = getParam<double>("Grid.Scaling");
+ // };
+ //
+ // Params params;
+ //
+ // auto elementSelectorStokes = [&](const auto& element)
+ // {
+ // double interface = params.amplitude * std::sin(( element.geometry().center()[0] -params.offset) / params.scaling * 2.0 * M_PI) + params.baseline;
+ // return element.geometry().center()[1] > interface;
+ // };
+ //
+ // auto elementSelectorDarcy = [&](const auto& element)
+ // {
+ // double interface = params.amplitude * std::sin(( element.geometry().center()[0] - params.offset) / params.scaling * 2.0 * M_PI) + params.baseline;
+ // return element.geometry().center()[1] < interface;
+ // };
+ //
+ // // subgrid Pointer
+ // auto stokesGridPtr = SubgridGridCreator<HostGrid>::makeGrid(hostGrid, elementSelectorStokes, "Stokes");
+ // auto darcyGridPtr = SubgridGridCreator<HostGrid>::makeGrid(hostGrid, elementSelectorDarcy, "Darcy");
+ //
+ // // we compute on the leaf grid view
+ // const auto& darcyGridView = darcyGridPtr->leafGridView();
+ // const auto& stokesGridView = stokesGridPtr->leafGridView();
+
+ // ************************************************************************************ //
+
+
+ // create the finite volume grid geometry
+ using StokesFVGridGeometry = GetPropType<StokesTypeTag, Properties::FVGridGeometry>;
+ auto stokesFvGridGeometry = std::make_shared<StokesFVGridGeometry>(stokesGridView);
+ stokesFvGridGeometry->update();
+ using DarcyFVGridGeometry = GetPropType<DarcyTypeTag, Properties::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 = GetPropType<StokesTypeTag, Properties::Problem>;
+ auto stokesProblem = std::make_shared<StokesProblem>(stokesFvGridGeometry, couplingManager);
+ using DarcyProblem = GetPropType<DarcyTypeTag, Properties::Problem>;
+ auto darcyProblem = std::make_shared<DarcyProblem>(darcyFvGridGeometry, couplingManager);
+
+ // the solution vector
+ Traits::SolutionVector sol;
+ sol[stokesCellCenterIdx].resize(stokesFvGridGeometry->numCellCenterDofs());
+ sol[stokesFaceIdx].resize(stokesFvGridGeometry->numFaceDofs());
+ sol[darcyIdx].resize(darcyFvGridGeometry->numDofs());
+
+ auto stokesSol = partial(sol, stokesCellCenterIdx, stokesFaceIdx);
+
+ stokesProblem->applyInitialSolution(stokesSol);
+ darcyProblem->applyInitialSolution(sol[darcyIdx]);
+
+ couplingManager->init(stokesProblem, darcyProblem, sol);
+
+ // the grid variables
+ using StokesGridVariables = GetPropType<StokesTypeTag, Properties::GridVariables>;
+ auto stokesGridVariables = std::make_shared<StokesGridVariables>(stokesProblem, stokesFvGridGeometry);
+ stokesGridVariables->init(stokesSol);
+ using DarcyGridVariables = GetPropType<DarcyTypeTag, Properties::GridVariables>;
+ auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyFvGridGeometry);
+ darcyGridVariables->init(sol[darcyIdx]);
+
+ // 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<StokesGridVariables, decltype(stokesSol)> stokesVtkWriter(*stokesGridVariables, stokesSol, stokesName);
+ GetPropType<StokesTypeTag, Properties::IOFields>::initOutputModule(stokesVtkWriter);
+
+ //****** uncomment the add analytical solution of v_x *****//
+ stokesVtkWriter.addField(stokesProblem->getAnalyticalVelocityX(), "analyticalV_x");
+
+ stokesVtkWriter.write(0.0);
+
+ VtkOutputModule<DarcyGridVariables, GetPropType<DarcyTypeTag, Properties::SolutionVector>> darcyVtkWriter(*darcyGridVariables, sol[darcyIdx], darcyName);
+ using DarcyVelocityOutput = GetPropType<DarcyTypeTag, Properties::VelocityOutput>;
+ darcyVtkWriter.addVelocityOutput(std::make_shared<DarcyVelocityOutput>(*darcyGridVariables));
+ GetPropType<DarcyTypeTag, Properties::IOFields>::initOutputModule(darcyVtkWriter);
+ darcyVtkWriter.write(0.0);
+
+ // the assembler for a stationary problem
+ using Assembler = 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);
+
+ // the linear solver
+ using LinearSolver = UMFPackBackend;
+ auto linearSolver = std::make_shared<LinearSolver>();
+
+ // the non-linear solver
+ using NewtonSolver = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
+ NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
+
+ // solve the non-linear system
+ nonLinearSolver.solve(sol);
+
+ // write vtk output
+ stokesVtkWriter.write(1.0);
+ darcyVtkWriter.write(1.0);
+
+ ////////////////////////////////////////////////////////////
+ // finalize, print dumux message to say goodbye
+ ////////////////////////////////////////////////////////////
+
+ // print dumux end message
+ if (mpiHelper.rank() == 0)
+ {
+ Parameters::print();
+ DumuxMessage::print(/*firstCall=*/false);
+ }
+
+ return 0;
+} // end main
+catch (Dumux::ParameterException &e)
+{
+ std::cerr << std::endl << e << " ---> Abort!" << std::endl;
+ return 1;
+}
+catch (Dune::DGFException & e)
+{
+ std::cerr << "DGF exception thrown (" << e <<
+ "). Most likely, the DGF file name is wrong "
+ "or the DGF file is corrupted, "
+ "e.g. missing hash at end of file or wrong number (dimensions) of entries."
+ << " ---> Abort!" << std::endl;
+ return 2;
+}
+catch (Dune::Exception &e)
+{
+ std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
+ return 3;
+}
+catch (...)
+{
+ std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
+ return 4;
+}
diff --git a/appl/coupling-ff-pm/monolithic/interface/ex_interface_coupling_ff-pm.input b/appl/coupling-ff-pm/monolithic/interface/ex_interface_coupling_ff-pm.input
new file mode 100644
index 0000000..4117b7f
--- /dev/null
+++ b/appl/coupling-ff-pm/monolithic/interface/ex_interface_coupling_ff-pm.input
@@ -0,0 +1,46 @@
+ # for dune-subgrid
+ [Grid]
+Positions0 = 0 1
+Positions1 = 0 0.2 0.3 0.65
+Cells0 = 100
+Cells1 = 10 50 18
+Baseline = 0.25 # [m]
+Amplitude = 0.04 # [m]
+Offset = 0.5 # [m]
+Scaling = 0.2 #[m]
+
+[Stokes.Grid]
+Verbosity = true
+Positions0 = 0.0 1.0
+Positions1 = 1.0 2.0
+Cells0 = 20
+Cells1 = 100
+Grading1 = 1
+
+[Darcy.Grid]
+Verbosity = true
+Positions0 = 0.0 1.0
+Positions1 = 0.0 1.0
+Cells0 = 20
+Cells1 = 20
+Grading1 = 1
+
+[Stokes.Problem]
+Name = stokes
+PressureDifference = 1e-9
+EnableInertiaTerms = false
+
+[Darcy.Problem]
+Name = darcy
+
+[Darcy.SpatialParams]
+Permeability = 1e-6 # m^2
+Porosity = 0.4
+AlphaBeaversJoseph = 1.0
+
+[Problem]
+Name = ex_ff-pm-interface
+EnableGravity = false
+
+[Vtk]
+AddVelocity = 1
diff --git a/appl/coupling-ff-pm/monolithic/interface/ex_interface_ffproblem.hh b/appl/coupling-ff-pm/monolithic/interface/ex_interface_ffproblem.hh
new file mode 100644
index 0000000..cbe437e
--- /dev/null
+++ b/appl/coupling-ff-pm/monolithic/interface/ex_interface_ffproblem.hh
@@ -0,0 +1,347 @@
+// -*- 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 The free flow sub problem
+ */
+#ifndef DUMUX_STOKES_SUBPROBLEM_HH
+#define DUMUX_STOKES_SUBPROBLEM_HH
+
+#include <dune/grid/yaspgrid.hh>
+
+//****** uncomment for the last exercise *****//
+// #include <dumux/io/grid/subgridgridcreator.hh>
+
+#include <dumux/material/fluidsystems/1pliquid.hh>
+#include <dumux/material/components/simpleh2o.hh>
+
+#include <dumux/freeflow/navierstokes/problem.hh>
+#include <dumux/discretization/staggered/freeflow/properties.hh>
+#include <dumux/freeflow/navierstokes/model.hh>
+
+namespace Dumux
+{
+template <class TypeTag>
+class StokesSubProblem;
+
+namespace Properties
+{
+// Create new type tags
+namespace TTag {
+struct StokesOneP { using InheritsFrom = std::tuple<NavierStokes, StaggeredFreeFlowModel>; };
+} // end namespace TTag
+
+// the fluid system
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::StokesOneP>
+{
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using type = FluidSystems::OnePLiquid<Scalar, Dumux::Components::SimpleH2O<Scalar> > ;
+};
+
+// Set the grid type
+template<class TypeTag>
+struct Grid<TypeTag, TTag::StokesOneP>
+{
+ static constexpr auto dim = 2;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using TensorGrid = Dune::YaspGrid<2, Dune::TensorProductCoordinates<Scalar, dim> >;
+
+//****** comment out for the last exercise *****//
+ using type = TensorGrid;
+
+//****** uncomment for the last exercise *****//
+ // using HostGrid = TensorGrid;
+ // using type = Dune::SubGrid<dim, HostGrid>;
+};
+
+// Set the problem property
+template<class TypeTag>
+struct Problem<TypeTag, TTag::StokesOneP> { using type = Dumux::StokesSubProblem<TypeTag> ; };
+
+template<class TypeTag>
+struct EnableFVGridGeometryCache<TypeTag, TTag::StokesOneP> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridFluxVariablesCache<TypeTag, TTag::StokesOneP> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridVolumeVariablesCache<TypeTag, TTag::StokesOneP> { static constexpr bool value = true; };
+}
+
+/*!
+ * \brief The free flow sub problem
+ */
+template <class TypeTag>
+class StokesSubProblem : public NavierStokesProblem<TypeTag>
+{
+ using ParentType = NavierStokesProblem<TypeTag>;
+
+ using GridView = GetPropType<TypeTag, Properties::GridView>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+
+ using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+
+ using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
+
+ using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using Element = typename GridView::template Codim<0>::Entity;
+
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
+ using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+
+ using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
+
+public:
+ StokesSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry, std::shared_ptr<CouplingManager> couplingManager)
+ : ParentType(fvGridGeometry, "Stokes"), eps_(1e-6), couplingManager_(couplingManager)
+ {
+ deltaP_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.PressureDifference");
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief Return the temperature within the domain in [K].
+ *
+ * This problem assumes a temperature of 10 degrees Celsius.
+ */
+ Scalar temperature() const
+ { return 273.15 + 10; } // 10°C
+
+ /*!
+ * \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.dofPosition();
+
+ // left/right wall
+ if (onRightBoundary_(globalPos) || (onLeftBoundary_(globalPos)))
+ {
+// values.setDirichlet(Indices::velocityXIdx);
+// values.setDirichlet(Indices::velocityYIdx);
+ values.setDirichlet(Indices::pressureIdx);
+ }
+ else
+ {
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ }
+
+ // coupling interface
+ if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
+ {
+ values.setCouplingNeumann(Indices::conti0EqIdx);
+ values.setCouplingNeumann(Indices::momentumYBalanceIdx);
+// values.setDirichlet(Indices::velocityXIdx); // assume no slip on interface
+ values.setBJS(Indices::momentumXBalanceIdx); // assume no slip on interface
+ }
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Dirichlet control volume.
+ *
+ * \param globalPos The global position
+ */
+ PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
+ {
+ PrimaryVariables values(0.0);
+ values = initialAtPos(globalPos);
+
+// if(onUpperBoundary_(globalPos))
+// {
+// values[Indices::velocityXIdx] = 0.0;
+// values[Indices::velocityYIdx] = 0.0;
+// }
+
+// if(onLeftBoundary_(globalPos))
+// values[Indices::pressureIdx] = deltaP_;
+// if(onRightBoundary_(globalPos))
+// values[Indices::pressureIdx] = 0.0;
+
+ return values;
+ }
+
+ /*!
+ * \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
+ */
+ template<class ElementVolumeVariables, class ElementFaceVariables>
+ NumEqVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFaceVariables& elemFaceVars,
+ const SubControlVolumeFace& scvf) const
+ {
+ NumEqVector values(0.0);
+
+ if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
+ {
+ values[Indices::conti0EqIdx] = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf);
+ values[Indices::momentumYBalanceIdx] = couplingManager().couplingData().momentumCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf);
+ }
+
+ 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_; }
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * \param globalPos The global position
+ */
+ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+ {
+ PrimaryVariables values(0.0);
+ //values[Indices::velocityYIdx] = -1e-6 * globalPos[0] * (this->fvGridGeometry().bBoxMax()[0] - globalPos[0]);
+ if(onLeftBoundary_(globalPos))
+ values[Indices::pressureIdx] = deltaP_;
+ if(onRightBoundary_(globalPos))
+ values[Indices::pressureIdx] = 0.0;
+
+ return values;
+ }
+
+ /*!
+ * \brief Returns the intrinsic permeability of required as input parameter for the Beavers-Joseph-Saffman boundary condition
+ */
+ Scalar permeability(const Element& element, const SubControlVolumeFace& scvf) const
+ {
+ return couplingManager().couplingData().darcyPermeability(element, 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());
+ }
+
+ /*!
+ * \brief calculate the analytical velocity in x direction based on Beavers & Joseph (1967)
+ */
+ void calculateAnalyticalVelocityX() const
+ {
+ analyticalVelocityX_.resize(this->fvGridGeometry().gridView().size(0));
+
+ using std::sqrt;
+ const Scalar dPdX = -deltaP_ / (this->fvGridGeometry().bBoxMax()[0] - this->fvGridGeometry().bBoxMin()[0]);
+ static const Scalar mu = FluidSystem::viscosity(temperature(), 1e5);
+ static const Scalar alpha = getParam<Scalar>("Darcy.SpatialParams.AlphaBeaversJoseph");
+ static const Scalar K = getParam<Scalar>("Darcy.SpatialParams.Permeability");
+ static const Scalar sqrtK = sqrt(K);
+ const Scalar sigma = (this->fvGridGeometry().bBoxMax()[1] - this->fvGridGeometry().bBoxMin()[1])/sqrtK;
+
+ const Scalar uB = -K/(2.0*mu) * ((sigma*sigma + 2.0*alpha*sigma) / (1.0 + alpha*sigma)) * dPdX;
+
+ for (const auto& element : elements(this->fvGridGeometry().gridView()))
+ {
+ const auto eIdx = this->fvGridGeometry().gridView().indexSet().index(element);
+ const Scalar y = element.geometry().center()[1] - this->fvGridGeometry().bBoxMin()[1];
+
+ const Scalar u = uB*(1.0 + alpha/sqrtK*y) + 1.0/(2.0*mu) * (y*y + 2*alpha*y*sqrtK) * dPdX;
+ analyticalVelocityX_[eIdx] = u;
+ }
+ }
+
+ /*!
+ * \brief Get the analytical velocity in x direction
+ */
+ const std::vector<Scalar>& getAnalyticalVelocityX() const
+ {
+ if(analyticalVelocityX_.empty())
+ calculateAnalyticalVelocityX();
+ return analyticalVelocityX_;
+ }
+
+ // \}
+
+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 deltaP_;
+
+ std::shared_ptr<CouplingManager> couplingManager_;
+
+ mutable std::vector<Scalar> analyticalVelocityX_;
+};
+} //end namespace
+
+#endif // DUMUX_STOKES_SUBPROBLEM_HH
diff --git a/appl/coupling-ff-pm/monolithic/interface/ex_interface_pmproblem.hh b/appl/coupling-ff-pm/monolithic/interface/ex_interface_pmproblem.hh
new file mode 100644
index 0000000..5b56c5a
--- /dev/null
+++ b/appl/coupling-ff-pm/monolithic/interface/ex_interface_pmproblem.hh
@@ -0,0 +1,271 @@
+// -*- 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 The porous medium flow sub problem
+ */
+#ifndef DUMUX_DARCY_SUBPROBLEM_HH
+#define DUMUX_DARCY_SUBPROBLEM_HH
+
+#include <dune/grid/yaspgrid.hh>
+
+//****** uncomment for the last exercise *****//
+// #include <dumux/io/grid/subgridgridcreator.hh>
+
+#include <dumux/discretization/cctpfa.hh>
+
+#include <dumux/porousmediumflow/1p/model.hh>
+#include <dumux/porousmediumflow/problem.hh>
+
+#include "../1pspatialparams.hh"
+
+#include <dumux/material/components/simpleh2o.hh>
+#include <dumux/material/fluidsystems/1pliquid.hh>
+
+namespace Dumux
+{
+template <class TypeTag>
+class DarcySubProblem;
+
+namespace Properties
+{
+// Create new type tags
+namespace TTag {
+struct DarcyOneP { using InheritsFrom = std::tuple<OneP, CCTpfaModel>; };
+} // end namespace TTag
+
+// Set the problem property
+template<class TypeTag>
+struct Problem<TypeTag, TTag::DarcyOneP> { using type = Dumux::DarcySubProblem<TypeTag>; };
+
+// the fluid system
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::DarcyOneP>
+{
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using type = FluidSystems::OnePLiquid<Scalar, Dumux::Components::SimpleH2O<Scalar> > ;
+};
+
+// Set the grid type
+template<class TypeTag>
+struct Grid<TypeTag, TTag::DarcyOneP>
+{
+ static constexpr auto dim = 2;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using TensorGrid = Dune::YaspGrid<2, Dune::TensorProductCoordinates<Scalar, dim> >;
+
+//****** comment out for the last exercise *****//
+ using type = TensorGrid;
+
+//****** uncomment for the last exercise *****//
+ // using HostGrid = TensorGrid;
+ // using type = Dune::SubGrid<dim, HostGrid>;
+};
+
+template<class TypeTag>
+struct SpatialParams<TypeTag, TTag::DarcyOneP> {
+ using type = OnePSpatialParams<GetPropType<TypeTag, FVGridGeometry>, GetPropType<TypeTag, Scalar>>;
+};
+
+} // end namespace Properties
+
+/*!
+ * \brief The porous medium flow sub problem
+ */
+template <class TypeTag>
+class DarcySubProblem : public PorousMediumFlowProblem<TypeTag>
+{
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
+ using GridView = GetPropType<TypeTag, Properties::GridView>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
+ using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
+ using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
+ using FVElementGeometry = typename GetPropType<TypeTag, Properties::FVGridGeometry>::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
+
+ using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
+
+public:
+ DarcySubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry,
+ std::shared_ptr<CouplingManager> couplingManager)
+ : ParentType(fvGridGeometry, "Darcy"), eps_(1e-7), couplingManager_(couplingManager)
+ {}
+
+ /*!
+ * \name Simulation steering
+ */
+ // \{
+
+ /*!
+ * \brief Return the temperature within the domain in [K].
+ *
+ */
+ Scalar temperature() const
+ { return 273.15 + 10; } // 10°C
+ // \}
+
+ /*!
+ * \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;
+
+ // set Neumann BCs to all boundaries first
+ values.setAllNeumann();
+
+ // set the coupling boundary condition at the interface
+ if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
+ values.setAllCouplingNeumann();
+
+ // set a Dirichlet boundary condition at the bottom
+ //if (onLowerBoundary_(scvf.center()))
+ // values.setAllDirichlet();
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Dirichlet control volume.
+ *
+ * \param element The element for which the Dirichlet boundary condition is set
+ * \param scvf The boundary subcontrolvolumeface
+ *
+ * For this method, the \a values parameter stores primary variables.
+ */
+ PrimaryVariables dirichlet(const Element &element, const SubControlVolumeFace &scvf) const
+ {
+ // set p = 0 at the bottom
+ PrimaryVariables values(0.0);
+ values = initial(element);
+
+ 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
+ {
+ // no-flow everywhere ...
+ NumEqVector values(0.0);
+
+ // ... except at the coupling interface
+ if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
+ values[Indices::conti0EqIdx] = couplingManager().couplingData().massCouplingCondition(element, 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 initial(const Element &element) const
+ {
+ return PrimaryVariables(0.0);
+ }
+
+ // \}
+
+ //! Set the coupling manager
+ void setCouplingManager(std::shared_ptr<CouplingManager> cm)
+ { couplingManager_ = cm; }
+
+ //! Get the coupling manager
+ const CouplingManager& couplingManager() const
+ { return *couplingManager_; }
+
+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_;
+ std::shared_ptr<CouplingManager> couplingManager_;
+};
+} //end namespace
+
+#endif //DUMUX_DARCY_SUBPROBLEM_HH
diff --git a/appl/coupling-ff-pm/monolithic/models/CMakeLists.txt b/appl/coupling-ff-pm/monolithic/models/CMakeLists.txt
new file mode 100644
index 0000000..c0f7fde
--- /dev/null
+++ b/appl/coupling-ff-pm/monolithic/models/CMakeLists.txt
@@ -0,0 +1,9 @@
+# executables for ex_interface_coupling_ff-pm
+dune_add_test(NAME ex_models_coupling_ff-pm
+ SOURCES ex_models_coupling_ff-pm.cc)
+
+# add tutorial to the common target
+add_dependencies(test_exercises ex_models_coupling_ff-pm)
+
+# add a symlink for each input file
+add_input_file_links()
diff --git a/appl/coupling-ff-pm/monolithic/models/ex_models_coupling_ff-pm.cc b/appl/coupling-ff-pm/monolithic/models/ex_models_coupling_ff-pm.cc
new file mode 100644
index 0000000..6f48a6d
--- /dev/null
+++ b/appl/coupling-ff-pm/monolithic/models/ex_models_coupling_ff-pm.cc
@@ -0,0 +1,283 @@
+// -*- 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/partial.hh>
+#include <dumux/linear/seqsolverbackend.hh>
+#include <dumux/assembly/fvassembler.hh>
+#include <dumux/assembly/diffmethod.hh>
+#include <dumux/discretization/method.hh>
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/staggeredvtkoutputmodule.hh>
+#include <dumux/io/grid/gridmanager.hh>
+
+#include <dumux/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 {
+
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::StokesNC>
+{
+ using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, Properties::TTag::DarcyOnePNC>;
+ using type = Dumux::StokesDarcyCouplingManager<Traits>;
+};
+
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::DarcyOnePNC>
+{
+ using Traits = StaggeredMultiDomainTraits<Properties::TTag::StokesNC, Properties::TTag::StokesNC, 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 = Properties::TTag::StokesNC;
+ using DarcyTypeTag = Properties::TTag::DarcyOnePNC;
+
+ // try to create a grid (from the given grid file or the input file)
+ // for both sub-domains
+ using DarcyGridManager = Dumux::GridManager<GetPropType<DarcyTypeTag, Properties::Grid>>;
+ DarcyGridManager darcyGridManager;
+ darcyGridManager.init("Darcy"); // pass parameter group
+
+ using StokesGridManager = Dumux::GridManager<GetPropType<StokesTypeTag, Properties::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 = GetPropType<StokesTypeTag, Properties::FVGridGeometry>;
+ auto stokesFvGridGeometry = std::make_shared<StokesFVGridGeometry>(stokesGridView);
+ stokesFvGridGeometry->update();
+ using DarcyFVGridGeometry = GetPropType<DarcyTypeTag, Properties::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 = GetPropType<StokesTypeTag, Properties::Problem>;
+ auto stokesProblem = std::make_shared<StokesProblem>(stokesFvGridGeometry, couplingManager);
+ using DarcyProblem = GetPropType<DarcyTypeTag, Properties::Problem>;
+ auto darcyProblem = std::make_shared<DarcyProblem>(darcyFvGridGeometry, couplingManager);
+
+ // initialize the fluidsystem (tabulation)
+ GetPropType<StokesTypeTag, Properties::FluidSystem>::init();
+
+ // get some time loop parameters
+ using Scalar = GetPropType<StokesTypeTag, Properties::Scalar>;
+ const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+ const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+ auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+
+ // instantiate time loop
+ auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(0, 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());
+
+ auto stokesSol = partial(sol, stokesCellCenterIdx, stokesFaceIdx);
+
+ stokesProblem->applyInitialSolution(stokesSol);
+ darcyProblem->applyInitialSolution(sol[darcyIdx]);
+
+ auto solOld = sol;
+
+ couplingManager->init(stokesProblem, darcyProblem, sol);
+
+ // the grid variables
+ using StokesGridVariables = GetPropType<StokesTypeTag, Properties::GridVariables>;
+ auto stokesGridVariables = std::make_shared<StokesGridVariables>(stokesProblem, stokesFvGridGeometry);
+ stokesGridVariables->init(stokesSol);
+ using DarcyGridVariables = GetPropType<DarcyTypeTag, Properties::GridVariables>;
+ auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyFvGridGeometry);
+ darcyGridVariables->init(sol[darcyIdx]);
+
+ // 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<StokesGridVariables, decltype(stokesSol)> stokesVtkWriter(*stokesGridVariables, stokesSol, stokesName);
+ GetPropType<StokesTypeTag, Properties::IOFields>::initOutputModule(stokesVtkWriter);
+ stokesVtkWriter.write(0.0);
+
+ VtkOutputModule<DarcyGridVariables, GetPropType<DarcyTypeTag, Properties::SolutionVector>> darcyVtkWriter(*darcyGridVariables, sol[darcyIdx], darcyName);
+ using DarcyVelocityOutput = GetPropType<DarcyTypeTag, Properties::VelocityOutput>;
+ darcyVtkWriter.addVelocityOutput(std::make_shared<DarcyVelocityOutput>(*darcyGridVariables));
+ GetPropType<DarcyTypeTag, Properties::IOFields>::initOutputModule(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>();
+
+ 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/appl/coupling-ff-pm/monolithic/models/ex_models_coupling_ff-pm.input b/appl/coupling-ff-pm/monolithic/models/ex_models_coupling_ff-pm.input
new file mode 100644
index 0000000..07bb4dd
--- /dev/null
+++ b/appl/coupling-ff-pm/monolithic/models/ex_models_coupling_ff-pm.input
@@ -0,0 +1,54 @@
+[TimeLoop]
+DtInitial = 1 # 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 = true
+EnableInertiaTerms = true
+MoleFraction = 0.0 # -
+Pressure = 1e5 # Pa
+Velocity = 1 # 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-7
+TargetSteps = 5
+
+[Vtk]
+AddVelocity = 1
+
+[Assembly]
+NumericDifference.BaseEpsilon = 1e-8
diff --git a/appl/coupling-ff-pm/monolithic/models/ex_models_ffproblem.hh b/appl/coupling-ff-pm/monolithic/models/ex_models_ffproblem.hh
new file mode 100644
index 0000000..853bc86
--- /dev/null
+++ b/appl/coupling-ff-pm/monolithic/models/ex_models_ffproblem.hh
@@ -0,0 +1,354 @@
+// -*- 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/1padapter.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
+{
+// Create new type tags
+namespace TTag {
+struct StokesNC { using InheritsFrom = std::tuple<NavierStokesNC, StaggeredFreeFlowModel>; };
+} // end namespace TTag
+
+// Set the grid type
+template<class TypeTag>
+struct Grid<TypeTag, TTag::StokesNC> { using type = Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
+
+// The fluid system
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::StokesNC>
+{
+ using H2OAir = FluidSystems::H2OAir<GetPropType<TypeTag, Properties::Scalar>>;
+ using type = FluidSystems::OnePAdapter<H2OAir, H2OAir::gasPhaseIdx>;
+};
+
+// Do not replace one equation with a total mass balance
+template<class TypeTag>
+struct ReplaceCompEqIdx<TypeTag, TTag::StokesNC> { static constexpr int value = 3; };
+
+// Use formulation based on mass fractions
+template<class TypeTag>
+struct UseMoles<TypeTag, TTag::StokesNC> { static constexpr bool value = true; };
+
+// Set the problem property
+template<class TypeTag>
+struct Problem<TypeTag, TTag::StokesNC> { using type = Dumux::StokesSubProblem<TypeTag> ; };
+
+template<class TypeTag>
+struct EnableFVGridGeometryCache<TypeTag, TTag::StokesNC> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridFluxVariablesCache<TypeTag, TTag::StokesNC> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridVolumeVariablesCache<TypeTag, TTag::StokesNC> { static constexpr bool value = true; };
+}
+
+/*!
+ * \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 = GetPropType<TypeTag, Properties::GridView>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+ using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+ using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
+
+ using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
+ using ElementFaceVariables = typename GetPropType<TypeTag, Properties::GridFaceVariables>::LocalView;
+ using FluidState = GetPropType<TypeTag, Properties::FluidState>;
+
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
+
+ using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
+ using TimeLoopPtr = std::shared_ptr<TimeLoop<Scalar>>;
+
+ static constexpr bool useMoles = GetPropType<TypeTag, Properties::ModelTraits>::useMoles();
+
+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
+ */
+ // \{
+
+ /*!
+ * \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(Indices::conti0EqIdx + 1);
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ }
+ else if(onRightBoundary_(globalPos))
+ {
+ values.setDirichlet(Indices::pressureIdx);
+ values.setOutflow(Indices::conti0EqIdx + 1);
+ }
+ else
+ {
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ values.setNeumann(Indices::conti0EqIdx);
+ values.setNeumann(Indices::conti0EqIdx + 1);
+ }
+
+ if (couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
+ {
+ values.setCouplingNeumann(Indices::conti0EqIdx);
+ values.setCouplingNeumann(Indices::conti0EqIdx + 1);
+ 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(element, fvGeometry, elemVolVars, elemFaceVars, scvf);
+
+ const auto massFlux = couplingManager().couplingData().massCouplingCondition(element, 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
+ {
+ PrimaryVariables values(0.0);
+
+ // This is only an approximation of the actual hydorostatic pressure gradient.
+ // Air is compressible (the density depends on pressure), thus the actual
+ // vertical pressure profile is non-linear.
+ // This discrepancy can lead to spurious flows at the outlet if the inlet
+ // velocity is chosen too small.
+ FluidState fluidState;
+ updateFluidStateForBC_(fluidState, pressure_);
+ const Scalar density = FluidSystem::density(fluidState, 0);
+ values[Indices::pressureIdx] = pressure_ - density*this->gravity()[1]*(this->fvGridGeometry().bBoxMax()[1] - globalPos[1]);
+
+
+ // gravity has negative sign
+ values[Indices::conti0EqIdx + 1] = 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 Element& element, const SubControlVolumeFace& scvf) const
+ {
+ return couplingManager().couplingData().darcyPermeability(element, 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(0, pressure);
+ fluidState.setSaturation(0, 1.0);
+ fluidState.setMoleFraction(0, 1, moleFraction_);
+ fluidState.setMoleFraction(0, 0, 1.0 - moleFraction_);
+
+ typename FluidSystem::ParameterCache paramCache;
+ paramCache.updatePhase(fluidState, 0);
+
+ const Scalar density = FluidSystem::density(fluidState, paramCache, 0);
+ fluidState.setDensity(0, density);
+
+ const Scalar molarDensity = FluidSystem::molarDensity(fluidState, paramCache, 0);
+ fluidState.setMolarDensity(0, molarDensity);
+
+ const Scalar enthalpy = FluidSystem::enthalpy(fluidState, paramCache, 0);
+ fluidState.setEnthalpy(0, 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/appl/coupling-ff-pm/monolithic/models/ex_models_pmproblem.hh b/appl/coupling-ff-pm/monolithic/models/ex_models_pmproblem.hh
new file mode 100644
index 0000000..f25401b
--- /dev/null
+++ b/appl/coupling-ff-pm/monolithic/models/ex_models_pmproblem.hh
@@ -0,0 +1,404 @@
+// -*- 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/cctpfa.hh>
+#include <dumux/io/gnuplotinterface.hh>
+#include <dumux/material/fluidsystems/1padapter.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
+{
+// Create new type tags
+namespace TTag {
+struct DarcyOnePNC { using InheritsFrom = std::tuple<OnePNC, CCTpfaModel>; };
+} // end namespace TTag
+
+// Set the problem property
+template<class TypeTag>
+struct Problem<TypeTag, TTag::DarcyOnePNC> { using type = Dumux::DarcySubProblem<TypeTag>; };
+
+// The fluid system
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::DarcyOnePNC>
+{
+ using H2OAir = FluidSystems::H2OAir<GetPropType<TypeTag, Properties::Scalar>>;
+ using type = FluidSystems::OnePAdapter<H2OAir, H2OAir::gasPhaseIdx>;
+};
+
+// Use moles
+template<class TypeTag>
+struct UseMoles<TypeTag, TTag::DarcyOnePNC> { static constexpr bool value = true; };
+
+// Do not replace one equation with a total mass balance
+template<class TypeTag>
+struct ReplaceCompEqIdx<TypeTag, TTag::DarcyOnePNC> { static constexpr int value = 3; };
+
+//! Use a model with constant tortuosity for the effective diffusivity
+template<class TypeTag>
+struct EffectiveDiffusivityModel<TypeTag, TTag::DarcyOnePNC>
+{ using type = DiffusivityConstantTortuosity<GetPropType<TypeTag, Properties::Scalar>>; };
+// Set the grid type
+template<class TypeTag>
+struct Grid<TypeTag, TTag::DarcyOnePNC> { using type = Dune::YaspGrid<2>; };
+
+// Set the spatial paramaters type
+template<class TypeTag>
+struct SpatialParams<TypeTag, TTag::DarcyOnePNC> {
+ using type = OnePSpatialParams<GetPropType<TypeTag, FVGridGeometry>, GetPropType<TypeTag, Scalar>>;
+};
+
+} // end namespace Properties
+
+template <class TypeTag>
+class DarcySubProblem : public PorousMediumFlowProblem<TypeTag>
+{
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
+ using GridView = GetPropType<TypeTag, Properties::GridView>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+ using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
+ using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
+ using FVElementGeometry = typename GetPropType<TypeTag, Properties::FVGridGeometry>::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
+
+ // copy some indices for convenience
+ using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+ enum {
+ // grid and world dimension
+ dim = GridView::dimension,
+ dimworld = GridView::dimensionworld,
+
+ // primary variable indices
+ conti0EqIdx = Indices::conti0EqIdx,
+ pressureIdx = Indices::pressureIdx,
+ phaseIdx = 0,
+ transportCompIdx = 1
+ };
+
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
+
+ using CouplingManager = GetPropType<TypeTag, Properties::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");
+
+ // 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];
+ // 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;
+
+ 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()));
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \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(element, 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
+ {
+ static const Scalar stokesPressure = getParamFromGroup<Scalar>("Stokes", "Problem.Pressure");
+
+ PrimaryVariables values(0.0);
+ values[Indices::pressureIdx] = stokesPressure;
+ values[transportCompIdx] = moleFraction_;
+ 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 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 Dumux
+
+#endif //DUMUX_DARCY_SUBPROBLEM_HH
diff --git a/appl/coupling-ff-pm/monolithic/turbulence/CMakeLists.txt b/appl/coupling-ff-pm/monolithic/turbulence/CMakeLists.txt
new file mode 100644
index 0000000..3315f6a
--- /dev/null
+++ b/appl/coupling-ff-pm/monolithic/turbulence/CMakeLists.txt
@@ -0,0 +1,10 @@
+# executables for ex_interface_coupling_ff-pm
+dune_add_test(NAME ex_turbulence_coupling_ff-pm
+ SOURCES ex_turbulence_coupling_ff-pm.cc)
+
+
+# add a symlink for each input file
+add_input_file_links()
+
+# add tutorial to the common target
+add_dependencies(test_exercises ex_turbulence_coupling_ff-pm)
diff --git a/appl/coupling-ff-pm/monolithic/turbulence/ex_turbulence_coupling_ff-pm.cc b/appl/coupling-ff-pm/monolithic/turbulence/ex_turbulence_coupling_ff-pm.cc
new file mode 100644
index 0000000..69a4063
--- /dev/null
+++ b/appl/coupling-ff-pm/monolithic/turbulence/ex_turbulence_coupling_ff-pm.cc
@@ -0,0 +1,284 @@
+// -*- 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 isothermal coupled Stokes/Darcy problem (1p2c/2p2c)
+ */
+#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/partial.hh>
+#include <dumux/linear/seqsolverbackend.hh>
+#include <dumux/assembly/fvassembler.hh>
+#include <dumux/assembly/diffmethod.hh>
+#include <dumux/discretization/method.hh>
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/staggeredvtkoutputmodule.hh>
+#include <dumux/io/grid/gridmanager.hh>
+
+#include <dumux/multidomain/staggeredtraits.hh>
+#include <dumux/multidomain/fvassembler.hh>
+#include <dumux/multidomain/newtonsolver.hh>
+
+#include <dumux/multidomain/boundary/stokesdarcy/couplingmanager.hh>
+
+#include "ex_turbulence_pmproblem.hh"
+#include "ex_turbulence_ffproblem.hh"
+
+namespace Dumux {
+namespace Properties {
+
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::StokesZeroEq>
+{
+ using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, Properties::TTag::DarcyTwoPTwoCNI>;
+ using type = Dumux::StokesDarcyCouplingManager<Traits>;
+};
+
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::DarcyTwoPTwoCNI>
+{
+ using Traits = StaggeredMultiDomainTraits<Properties::TTag::StokesZeroEq, Properties::TTag::StokesZeroEq, 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 = Properties::TTag::StokesZeroEq;
+ using DarcyTypeTag = Properties::TTag::DarcyTwoPTwoCNI;
+
+ // try to create a grid (from the given grid file or the input file)
+ // for both sub-domains
+ using DarcyGridManager = Dumux::GridManager<GetPropType<DarcyTypeTag, Properties::Grid>>;
+ DarcyGridManager darcyGridManager;
+ darcyGridManager.init("Darcy"); // pass parameter group
+
+ using StokesGridManager = Dumux::GridManager<GetPropType<StokesTypeTag, Properties::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 = GetPropType<StokesTypeTag, Properties::FVGridGeometry>;
+ auto stokesFvGridGeometry = std::make_shared<StokesFVGridGeometry>(stokesGridView);
+ stokesFvGridGeometry->update();
+ using DarcyFVGridGeometry = GetPropType<DarcyTypeTag, Properties::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 = GetPropType<StokesTypeTag, Properties::Problem>;
+ auto stokesProblem = std::make_shared<StokesProblem>(stokesFvGridGeometry, couplingManager);
+ using DarcyProblem = GetPropType<DarcyTypeTag, Properties::Problem>;
+ auto darcyProblem = std::make_shared<DarcyProblem>(darcyFvGridGeometry, couplingManager);
+
+ // initialize the fluidsystem (tabulation)
+ GetPropType<StokesTypeTag, Properties::FluidSystem>::init();
+
+ // get some time loop parameters
+ using Scalar = GetPropType<StokesTypeTag, Properties::Scalar>;
+ const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+ const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+ auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+
+ // instantiate time loop
+ auto timeLoop = std::make_shared<TimeLoop<Scalar>>(0, dt, tEnd);
+ timeLoop->setMaxTimeStepSize(maxDt);
+
+ // set timeloop for the subproblems, needed for boundary value variations
+ 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());
+
+ auto stokesSol = partial(sol, stokesCellCenterIdx, stokesFaceIdx);
+
+ stokesProblem->applyInitialSolution(stokesSol);
+ darcyProblem->applyInitialSolution(sol[darcyIdx]);
+
+ auto solOld = sol;
+
+ couplingManager->init(stokesProblem, darcyProblem, sol);
+
+ // TODO: update static wall properties
+ // TODO: update dynamic wall properties
+
+ // the grid variables
+ using StokesGridVariables = GetPropType<StokesTypeTag, Properties::GridVariables>;
+ auto stokesGridVariables = std::make_shared<StokesGridVariables>(stokesProblem, stokesFvGridGeometry);
+ stokesGridVariables->init(stokesSol);
+ using DarcyGridVariables = GetPropType<DarcyTypeTag, Properties::GridVariables>;
+ auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyFvGridGeometry);
+ darcyGridVariables->init(sol[darcyIdx]);
+
+ // 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<StokesGridVariables, decltype(stokesSol)> stokesVtkWriter(*stokesGridVariables, stokesSol, stokesName);
+ GetPropType<StokesTypeTag, Properties::IOFields>::initOutputModule(stokesVtkWriter);
+ stokesVtkWriter.write(0.0);
+
+ VtkOutputModule<DarcyGridVariables, GetPropType<DarcyTypeTag, Properties::SolutionVector>> darcyVtkWriter(*darcyGridVariables, sol[darcyIdx], darcyName);
+ using DarcyVelocityOutput = GetPropType<DarcyTypeTag, Properties::VelocityOutput>;
+ darcyVtkWriter.addVelocityOutput(std::make_shared<DarcyVelocityOutput>(*darcyGridVariables));
+ GetPropType<DarcyTypeTag, Properties::IOFields>::initOutputModule(darcyVtkWriter);
+ darcyVtkWriter.write(0.0);
+
+ // the assembler 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 non-linear solver
+ using NewtonSolver = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
+ NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
+
+ // time loop
+ 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;
+
+ // TODO: update dynamic wall properties
+
+ // post time step treatment of Darcy problem
+ darcyProblem->postTimeStep(sol[darcyIdx], *darcyGridVariables, timeLoop->timeStepSize());
+
+ // advance grid variables to the next time step
+ stokesGridVariables->advanceTimeStep();
+ darcyGridVariables->advanceTimeStep();
+
+ // advance to the time loop to the next step
+ timeLoop->advanceTimeStep();
+
+ // write vtk output
+ 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/appl/coupling-ff-pm/monolithic/turbulence/ex_turbulence_coupling_ff-pm.input b/appl/coupling-ff-pm/monolithic/turbulence/ex_turbulence_coupling_ff-pm.input
new file mode 100644
index 0000000..00e0246
--- /dev/null
+++ b/appl/coupling-ff-pm/monolithic/turbulence/ex_turbulence_coupling_ff-pm.input
@@ -0,0 +1,68 @@
+[TimeLoop]
+DtInitial = 1e-1 # [s]
+MaxTimeStepSize = 43200 # [s] (12 hours)
+TEnd = 864000 # [s] (6 days)
+
+[Stokes.Grid]
+Positions0 = 0.0 0.25
+Positions1 = 0.25 0.5
+Grading0 = 1.0
+Grading1 = 1.0
+Cells0 = 15
+Cells1 = 20
+Verbosity = true
+
+[Darcy.Grid]
+Positions0 = 0.0 0.25
+Positions1 = 0.0 0.25
+Cells0 = 15
+Cells1 = 10
+Grading0 = 1.0
+Grading1 = 1.0
+Verbosity = true
+
+[Stokes.Problem]
+Name = stokes
+RefVelocity = 3.5 # [m/s]
+RefPressure = 1e5 # [Pa]
+refMoleFrac = 0 # [-]
+RefTemperature = 298.15 # [K]
+EnableInertiaTerms = true
+
+[Darcy.Problem]
+Name = darcy
+Pressure = 1e5
+Saturation = 0.5 # initial Sw
+Temperature = 298.15 # [K]
+InitPhasePresence = 3 # bothPhases
+
+[Darcy.SpatialParams]
+Porosity = 0.41
+Permeability = 2.65e-10
+AlphaBeaversJoseph = 1.0
+Swr = 0.005
+Snr = 0.01
+VgAlpha = 6.371e-4
+VgN = 6.9
+
+[Problem]
+Name = ex_coupling_turbulence_ff-pm
+EnableGravity = true
+InterfaceDiffusionCoefficientAvg = Harmonic
+
+[Vtk]
+AddVelocity = true
+WriteFaceData = false
+
+[Newton]
+MaxSteps = 12
+MaxRelativeShift = 1e-5
+
+[Assembly]
+NumericDifferenceMethod = 0
+NumericDifference.BaseEpsilon = 1e-8
+
+[Component]
+SolidDensity = 2700
+SolidThermalConductivity = 2.8
+SolidHeatCapacity = 790
diff --git a/appl/coupling-ff-pm/monolithic/turbulence/ex_turbulence_ffproblem.hh b/appl/coupling-ff-pm/monolithic/turbulence/ex_turbulence_ffproblem.hh
new file mode 100644
index 0000000..e6abb8e
--- /dev/null
+++ b/appl/coupling-ff-pm/monolithic/turbulence/ex_turbulence_ffproblem.hh
@@ -0,0 +1,387 @@
+// -*- 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 The free-flow sub problem
+ */
+#ifndef DUMUX_FREEFLOW1P2C_SUBPROBLEM_HH
+#define DUMUX_FREEFLOW1P2C_SUBPROBLEM_HH
+
+#include <dune/grid/yaspgrid.hh>
+
+#include <dumux/material/fluidsystems/1padapter.hh>
+#include <dumux/material/fluidsystems/h2oair.hh>
+#include <dumux/discretization/staggered/freeflow/properties.hh>
+
+#include <dumux/freeflow/compositional/navierstokesncmodel.hh>
+#include <dumux/freeflow/navierstokes/problem.hh>
+
+namespace Dumux {
+
+template <class TypeTag>
+class FreeFlowSubProblem;
+
+namespace Properties {
+
+// Create new type tags
+namespace TTag {
+struct StokesZeroEq { using InheritsFrom = std::tuple<NavierStokesNCNI, StaggeredFreeFlowModel>; };
+} // end namespace TTag
+
+// Set the grid type
+template<class TypeTag>
+struct Grid<TypeTag, TTag::StokesZeroEq> { using type = Dune::YaspGrid<2, Dune::TensorProductCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
+
+// The fluid system
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::StokesZeroEq>
+{
+ using H2OAir = FluidSystems::H2OAir<GetPropType<TypeTag, Properties::Scalar>>;
+ static constexpr auto phaseIdx = H2OAir::gasPhaseIdx; // simulate the air phase
+ using type = FluidSystems::OnePAdapter<H2OAir, phaseIdx>;
+};
+
+template<class TypeTag>
+struct ReplaceCompEqIdx<TypeTag, TTag::StokesZeroEq> { static constexpr int value = 3; };
+
+// Use formulation based on mass fractions
+template<class TypeTag>
+struct UseMoles<TypeTag, TTag::StokesZeroEq> { static constexpr bool value = true; };
+
+// Set the problem property
+template<class TypeTag>
+struct Problem<TypeTag, TTag::StokesZeroEq> { using type = Dumux::FreeFlowSubProblem<TypeTag> ; };
+
+template<class TypeTag>
+struct EnableFVGridGeometryCache<TypeTag, TTag::StokesZeroEq> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridFluxVariablesCache<TypeTag, TTag::StokesZeroEq> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridVolumeVariablesCache<TypeTag, TTag::StokesZeroEq> { static constexpr bool value = true; };
+}
+
+/*!
+ * \brief The free-flow sub problem
+ */
+template <class TypeTag>
+class FreeFlowSubProblem : public NavierStokesProblem<TypeTag>
+{
+ using ParentType = NavierStokesProblem<TypeTag>;
+
+ using GridView = GetPropType<TypeTag, Properties::GridView>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+ using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+ using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
+
+ using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
+ using ElementFaceVariables = typename GetPropType<TypeTag, Properties::GridFaceVariables>::LocalView;
+ using FluidState = GetPropType<TypeTag, Properties::FluidState>;
+
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
+
+ using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
+ using TimeLoopPtr = std::shared_ptr<TimeLoop<Scalar>>;
+
+ using DiffusionCoefficientAveragingType = typename StokesDarcyCouplingOptions::DiffusionCoefficientAveragingType;
+
+ static constexpr bool useMoles = GetPropType<TypeTag, Properties::ModelTraits>::useMoles();
+
+public:
+ FreeFlowSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry, std::shared_ptr<CouplingManager> couplingManager)
+ : ParentType(fvGridGeometry, "Stokes"), eps_(1e-6), couplingManager_(couplingManager)
+ {
+ refVelocity_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.RefVelocity");
+ refPressure_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.RefPressure");
+ refMoleFrac_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.refMoleFrac");
+ refTemperature_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.RefTemperature");
+
+ diffCoeffAvgType_ = StokesDarcyCouplingOptions::stringToEnum(DiffusionCoefficientAveragingType{},
+ getParamFromGroup<std::string>(this->paramGroup(), "Problem.InterfaceDiffusionCoefficientAvg"));
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief Return the temperature within the domain in [K].
+ */
+ Scalar temperature() const
+ { return refTemperature_; }
+
+ /*!
+ * \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(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ values.setDirichlet(Indices::conti0EqIdx + 1);
+ values.setDirichlet(Indices::energyEqIdx);
+ }
+
+ if (onLowerBoundary_(globalPos))
+ {
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ values.setNeumann(Indices::conti0EqIdx);
+ values.setNeumann(Indices::conti0EqIdx + 1);
+ values.setNeumann(Indices::energyEqIdx);
+ }
+
+ if (onUpperBoundary_(globalPos))
+ {
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ values.setNeumann(Indices::conti0EqIdx);
+ values.setNeumann(Indices::conti0EqIdx + 1);
+ values.setNeumann(Indices::energyEqIdx);
+ }
+
+ if (onRightBoundary_(globalPos))
+ {
+ values.setDirichlet(Indices::pressureIdx);
+ values.setOutflow(Indices::conti0EqIdx + 1);
+ values.setOutflow(Indices::energyEqIdx);
+ }
+
+ if (couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
+ {
+ values.setCouplingNeumann(Indices::conti0EqIdx);
+ values.setCouplingNeumann(Indices::conti0EqIdx + 1);
+ values.setCouplingNeumann(Indices::energyEqIdx);
+ 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(element, fvGeometry, elemVolVars, elemFaceVars, scvf);
+
+ const auto massFlux = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf, diffCoeffAvgType_);
+ values[Indices::conti0EqIdx] = massFlux[0];
+ values[Indices::conti0EqIdx + 1] = massFlux[1];
+ values[Indices::energyEqIdx] = couplingManager().couplingData().energyCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf, diffCoeffAvgType_);
+ }
+ 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, refTemperature(), refPressure(), refMoleFrac());
+
+ const Scalar density = FluidSystem::density(fluidState, 0);
+
+ PrimaryVariables values(0.0);
+ values[Indices::pressureIdx] = refPressure() + density*this->gravity()[1]*(globalPos[1] - this->fvGridGeometry().bBoxMin()[1]);
+ values[Indices::conti0EqIdx + 1] = refMoleFrac();
+ values[Indices::velocityXIdx] = refVelocity();
+ values[Indices::temperatureIdx] = refTemperature();
+
+ if(onUpperBoundary_(globalPos) || onLowerBoundary_(globalPos))
+ values[Indices::velocityXIdx] = 0.0;
+
+ return values;
+ }
+
+ //! \brief Returns the reference velocity.
+ const Scalar refVelocity() const
+ { return refVelocity_ ;}
+
+ //! \brief Returns the reference pressure.
+ const Scalar refPressure() const
+ { return refPressure_; }
+
+ //! \brief Returns the reference mass fraction.
+ const Scalar refMoleFrac() const
+ { return refMoleFrac_; }
+
+ //! \brief Returns the reference temperature.
+ const Scalar refTemperature() const
+ { return refTemperature_; }
+
+
+ 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 Element& element, const SubControlVolumeFace& scvf) const
+ {
+ return couplingManager().problem(CouplingManager::darcyIdx).spatialParams().permeabilityAtPos(scvf.center());
+ }
+
+ /*!
+ * \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 temperature,
+ const Scalar pressure, const Scalar moleFraction) const
+ {
+ fluidState.setTemperature(temperature);
+ fluidState.setPressure(0, pressure);
+ fluidState.setSaturation(0, 1.0);
+ fluidState.setMoleFraction(0, 1, moleFraction);
+ fluidState.setMoleFraction(0, 0, 1.0 - moleFraction);
+
+ typename FluidSystem::ParameterCache paramCache;
+ paramCache.updatePhase(fluidState, 0);
+
+ const Scalar density = FluidSystem::density(fluidState, paramCache, 0);
+ fluidState.setDensity(0, density);
+
+ const Scalar molarDensity = FluidSystem::molarDensity(fluidState, paramCache, 0);
+ fluidState.setMolarDensity(0, molarDensity);
+
+ const Scalar enthalpy = FluidSystem::enthalpy(fluidState, paramCache, 0);
+ fluidState.setEnthalpy(0, enthalpy);
+ }
+
+ // the height of the free-flow domain
+ const Scalar height_() const
+ { return this->fvGridGeometry().bBoxMax()[1] - this->fvGridGeometry().bBoxMin()[1]; }
+
+ Scalar eps_;
+
+ Scalar refVelocity_;
+ Scalar refPressure_;
+ Scalar refMoleFrac_;
+ Scalar refTemperature_;
+
+ TimeLoopPtr timeLoop_;
+
+ std::shared_ptr<CouplingManager> couplingManager_;
+
+ DiffusionCoefficientAveragingType diffCoeffAvgType_;
+};
+} //end namespace
+
+#endif // DUMUX_STOKES1P2C_SUBPROBLEM_HH
diff --git a/appl/coupling-ff-pm/monolithic/turbulence/ex_turbulence_pmproblem.hh b/appl/coupling-ff-pm/monolithic/turbulence/ex_turbulence_pmproblem.hh
new file mode 100644
index 0000000..0ded6f6
--- /dev/null
+++ b/appl/coupling-ff-pm/monolithic/turbulence/ex_turbulence_pmproblem.hh
@@ -0,0 +1,337 @@
+// -*- 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 The porous medium sub problem
+ */
+#ifndef DUMUX_DARCY2P2C_SUBPROBLEM_HH
+#define DUMUX_DARCY2P2C_SUBPROBLEM_HH
+
+#include <dune/grid/yaspgrid.hh>
+
+#include <dumux/discretization/cctpfa.hh>
+
+#include <dumux/porousmediumflow/2p2c/model.hh>
+#include <dumux/porousmediumflow/problem.hh>
+
+#include <dumux/material/fluidsystems/h2oair.hh>
+
+#include "../2pspatialparams.hh"
+
+namespace Dumux
+{
+template <class TypeTag>
+class DarcySubProblem;
+
+namespace Properties
+{
+// Create new type tags
+namespace TTag {
+struct DarcyTwoPTwoCNI { using InheritsFrom = std::tuple<TwoPTwoCNI, CCTpfaModel>; };
+} // end namespace TTag
+
+// Set the problem property
+template<class TypeTag>
+struct Problem<TypeTag, TTag::DarcyTwoPTwoCNI> { using type = Dumux::DarcySubProblem<TypeTag>; };
+
+// the fluid system
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::DarcyTwoPTwoCNI> { using type = FluidSystems::H2OAir<GetPropType<TypeTag, Properties::Scalar>>; };
+
+//! Set the default formulation to pw-Sn: This can be over written in the problem.
+template<class TypeTag>
+struct Formulation<TypeTag, TTag::DarcyTwoPTwoCNI>
+{ static constexpr auto value = TwoPFormulation::p1s0; };
+
+// The gas component balance (air) is replaced by the total mass balance
+template<class TypeTag>
+struct ReplaceCompEqIdx<TypeTag, TTag::DarcyTwoPTwoCNI> { static constexpr int value = 3; };
+
+// Set the grid type
+template<class TypeTag>
+struct Grid<TypeTag, TTag::DarcyTwoPTwoCNI> { using type = Dune::YaspGrid<2, Dune::TensorProductCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
+
+template<class TypeTag>
+struct UseMoles<TypeTag, TTag::DarcyTwoPTwoCNI> { static constexpr bool value = true; };
+
+template<class TypeTag>
+struct SpatialParams<TypeTag, TTag::DarcyTwoPTwoCNI> {
+ using type = TwoPSpatialParams<GetPropType<TypeTag, FVGridGeometry>, GetPropType<TypeTag, Scalar>>;
+};
+
+} // end namespace Properties
+
+/*!
+ * \brief The porous medium sub problem
+ */
+template <class TypeTag>
+class DarcySubProblem : public PorousMediumFlowProblem<TypeTag>
+{
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
+ using GridView = GetPropType<TypeTag, Properties::GridView>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
+ using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
+ using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
+ using FVElementGeometry = typename GetPropType<TypeTag, Properties::FVGridGeometry>::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
+ using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
+
+ using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+
+ // copy some indices for convenience
+ using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+ enum {
+ // primary variable indices
+ conti0EqIdx = Indices::conti0EqIdx,
+ contiWEqIdx = Indices::conti0EqIdx + FluidSystem::H2OIdx,
+ contiNEqIdx = Indices::conti0EqIdx + FluidSystem::AirIdx,
+ pressureIdx = Indices::pressureIdx,
+ switchIdx = Indices::switchIdx
+ };
+
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
+ using TimeLoopPtr = std::shared_ptr<TimeLoop<Scalar>>;
+
+ using DiffusionCoefficientAveragingType = typename StokesDarcyCouplingOptions::DiffusionCoefficientAveragingType;
+
+public:
+ DarcySubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry,
+ std::shared_ptr<CouplingManager> couplingManager)
+ : ParentType(fvGridGeometry, "Darcy"), eps_(1e-7), couplingManager_(couplingManager)
+ {
+ pressure_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Pressure");
+ initialSw_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Saturation");
+ temperature_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Temperature");
+ initialPhasePresence_ = getParamFromGroup<int>(this->paramGroup(), "Problem.InitPhasePresence");
+
+ diffCoeffAvgType_ = StokesDarcyCouplingOptions::stringToEnum(DiffusionCoefficientAveragingType{},
+ getParamFromGroup<std::string>(this->paramGroup(), "Problem.InterfaceDiffusionCoefficientAvg"));
+ }
+
+ /*!
+ * \name Simulation steering
+ */
+ // \{
+
+ template<class SolutionVector, class GridVariables>
+ void postTimeStep(const SolutionVector& curSol,
+ const GridVariables& gridVariables,
+ const Scalar timeStepSize)
+
+ {
+ // compute the mass in the entire domain
+ Scalar massWater = 0.0;
+
+ // bulk elements
+ 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)
+ {
+ massWater += volVars.massFraction(phaseIdx, FluidSystem::H2OIdx)*volVars.density(phaseIdx)
+ * scv.volume() * volVars.saturation(phaseIdx) * volVars.porosity() * volVars.extrusionFactor();
+ }
+ }
+ }
+
+ std::cout << "mass of water is: " << massWater << std::endl;
+ }
+
+ /*!
+ * \brief Return the temperature within the domain in [K].
+ */
+ Scalar temperature() const
+ { return temperature_; }
+ // \}
+
+ /*!
+ * \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 Dirichlet control volume.
+ *
+ * \param element The element for which the Dirichlet boundary condition is set
+ * \param scvf The boundary subcontrolvolumeface
+ *
+ * For this method, the \a values parameter stores primary variables.
+ */
+ PrimaryVariables dirichlet(const Element &element, const SubControlVolumeFace &scvf) const
+ {
+ PrimaryVariables values(0.0);
+ values = initialAtPos(scvf.center());
+
+ 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
+ *
+ */
+ 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))
+ {
+ const auto massFlux = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, scvf, diffCoeffAvgType_);
+
+ for(int i = 0; i< massFlux.size(); ++i)
+ values[i] = massFlux[i];
+
+ values[Indices::energyEqIdx] = couplingManager().couplingData().energyCouplingCondition(element, fvGeometry, elemVolVars, scvf, diffCoeffAvgType_);
+ }
+
+ 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
+ *
+ * For this method, the \a values variable stores the rate mass
+ * of a component is generated or annihilate per volume
+ * unit. Positive values mean that mass is created, negative ones
+ * mean that it vanishes.
+ */
+ 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.
+ *
+ * For this method, the \a priVars parameter stores primary
+ * variables.
+ */
+ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+ {
+ PrimaryVariables values(0.0);
+ values.setState(initialPhasePresence_);
+
+ values[pressureIdx] = pressure_ + 1. * this->gravity()[1] * (globalPos[1] - this->fvGridGeometry().bBoxMax()[1]);
+ values[switchIdx] = initialSw_;
+ values[Indices::temperatureIdx] = temperature_;
+
+ 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_; }
+
+ 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 pressure_;
+ Scalar initialSw_;
+ Scalar temperature_;
+ int initialPhasePresence_;
+
+ TimeLoopPtr timeLoop_;
+
+ Scalar eps_;
+
+ std::shared_ptr<CouplingManager> couplingManager_;
+ DiffusionCoefficientAveragingType diffCoeffAvgType_;
+};
+} //end namespace
+
+#endif //DUMUX_DARCY2P2C_SUBPROBLEM_HH
diff --git a/appl/coupling-ff-pm/solution/CMakeLists.txt b/appl/coupling-ff-pm/solution/CMakeLists.txt
new file mode 100644
index 0000000..2ffa314
--- /dev/null
+++ b/appl/coupling-ff-pm/solution/CMakeLists.txt
@@ -0,0 +1 @@
+add_subdirectory(monolithic)
diff --git a/appl/coupling-ff-pm/solution/monolithic/1pspatialparams.hh b/appl/coupling-ff-pm/solution/monolithic/1pspatialparams.hh
new file mode 100644
index 0000000..a042388
--- /dev/null
+++ b/appl/coupling-ff-pm/solution/monolithic/1pspatialparams.hh
@@ -0,0 +1,91 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup OnePTests
+ * \brief The spatial parameters class for the test problem using the 1p cc model
+ */
+#ifndef DUMUX_1P_TEST_SPATIALPARAMS_HH
+#define DUMUX_1P_TEST_SPATIALPARAMS_HH
+
+#include <dumux/material/spatialparams/fv1p.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup OnePModel
+ *
+ * \brief The spatial parameters class for the test problem using the
+ * 1p cc model
+ */
+template<class FVGridGeometry, class Scalar>
+class OnePSpatialParams
+: public FVSpatialParamsOneP<FVGridGeometry, Scalar, OnePSpatialParams<FVGridGeometry, Scalar>>
+{
+ using GridView = typename FVGridGeometry::GridView;
+ using ParentType = FVSpatialParamsOneP<FVGridGeometry, Scalar, OnePSpatialParams<FVGridGeometry, Scalar>>;
+
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+public:
+ // export permeability type
+ using PermeabilityType = Scalar;
+
+ OnePSpatialParams(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+ : ParentType(fvGridGeometry)
+ {
+ permeability_ = getParam<Scalar>("Darcy.SpatialParams.Permeability");
+ porosity_ = getParam<Scalar>("Darcy.SpatialParams.Porosity");
+ alphaBJ_ = getParam<Scalar>("Darcy.SpatialParams.AlphaBeaversJoseph");
+ }
+
+ /*!
+ * \brief Function for defining the (intrinsic) permeability \f$[m^2]\f$.
+ *
+ * \param globalPos The global position
+ * \return the intrinsic permeability
+ */
+ PermeabilityType permeabilityAtPos(const GlobalPosition& globalPos) const
+ { return permeability_; }
+
+ /*! \brief Define the porosity in [-].
+ *
+ * \param globalPos The global position
+ */
+ Scalar porosityAtPos(const GlobalPosition& globalPos) const
+ { return porosity_; }
+
+ /*! \brief Define the Beavers-Joseph coefficient in [-].
+ *
+ * \param globalPos The global position
+ */
+ Scalar beaversJosephCoeffAtPos(const GlobalPosition& globalPos) const
+ { return alphaBJ_; }
+
+
+private:
+ Scalar permeability_;
+ Scalar porosity_;
+ Scalar alphaBJ_;
+};
+
+} // end namespace
+
+#endif
diff --git a/appl/coupling-ff-pm/solution/monolithic/2pspatialparams.hh b/appl/coupling-ff-pm/solution/monolithic/2pspatialparams.hh
new file mode 100644
index 0000000..622e4f5
--- /dev/null
+++ b/appl/coupling-ff-pm/solution/monolithic/2pspatialparams.hh
@@ -0,0 +1,132 @@
+// -*- 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
+ *
+ * \brief The spatial parameters class for the test problem using the 2p cc model
+ */
+template<class FVGridGeometry, class Scalar>
+class TwoPSpatialParams
+: public FVSpatialParams<FVGridGeometry, Scalar, TwoPSpatialParams<FVGridGeometry, Scalar>>
+{
+ using GridView = typename FVGridGeometry::GridView;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using ParentType = FVSpatialParams<FVGridGeometry, Scalar, TwoPSpatialParams<FVGridGeometry, Scalar>>;
+
+ 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 Dumux
+
+#endif
diff --git a/appl/coupling-ff-pm/solution/monolithic/CMakeLists.txt b/appl/coupling-ff-pm/solution/monolithic/CMakeLists.txt
new file mode 100644
index 0000000..bae4c26
--- /dev/null
+++ b/appl/coupling-ff-pm/solution/monolithic/CMakeLists.txt
@@ -0,0 +1,3 @@
+add_subdirectory(interface)
+add_subdirectory(models)
+add_subdirectory(turbulence)
diff --git a/appl/coupling-ff-pm/solution/monolithic/interface/CMakeLists.txt b/appl/coupling-ff-pm/solution/monolithic/interface/CMakeLists.txt
new file mode 100644
index 0000000..b0eb179
--- /dev/null
+++ b/appl/coupling-ff-pm/solution/monolithic/interface/CMakeLists.txt
@@ -0,0 +1,22 @@
+dune_add_test(NAME orig_ex_interface_coupling_ff-pm
+ SOURCES ex_interface_coupling_ff-pm.cc
+ COMPILE_DEFINITIONS EXNUMBER=0)
+
+dune_add_test(NAME sol_a_ex_interface_coupling_ff-pm
+ SOURCES ex_interface_coupling_ff-pm.cc
+ COMPILE_DEFINITIONS EXNUMBER=1)
+
+dune_add_test(NAME sol_b_ex_interface_coupling_ff-pm
+ SOURCES ex_interface_coupling_ff-pm.cc
+ COMPILE_DEFINITIONS EXNUMBER=2)
+
+dune_add_test(NAME sol_c_ex_interface_coupling_ff-pm
+ SOURCES ex_interface_coupling_ff-pm.cc
+ COMPILE_DEFINITIONS EXNUMBER=3)
+
+
+# add exercise to the common target
+add_dependencies(test_exercises sol_a_ex_interface_coupling_ff-pm sol_b_ex_interface_coupling_ff-pm sol_b_ex_interface_coupling_ff-pm)
+
+# add a symlink for each input file
+add_input_file_links()
diff --git a/appl/coupling-ff-pm/solution/monolithic/interface/ex_interface_coupling_ff-pm.cc b/appl/coupling-ff-pm/solution/monolithic/interface/ex_interface_coupling_ff-pm.cc
new file mode 100644
index 0000000..868ccc9
--- /dev/null
+++ b/appl/coupling-ff-pm/solution/monolithic/interface/ex_interface_coupling_ff-pm.cc
@@ -0,0 +1,274 @@
+// -*- 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 <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/istl/io.hh>
+#include <dune/grid/yaspgrid.hh>
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/dumuxmessage.hh>
+#include <dumux/common/partial.hh>
+#include <dumux/linear/seqsolverbackend.hh>
+#include <dumux/assembly/fvassembler.hh>
+#include <dumux/assembly/diffmethod.hh>
+#include <dumux/discretization/method.hh>
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/staggeredvtkoutputmodule.hh>
+#include <dumux/io/grid/gridmanager.hh>
+
+#include <dumux/multidomain/staggeredtraits.hh>
+#include <dumux/multidomain/fvassembler.hh>
+#include <dumux/multidomain/newtonsolver.hh>
+
+#include <dumux/multidomain/boundary/stokesdarcy/couplingmanager.hh>
+
+#include "ex_interface_pmproblem.hh"
+#include "ex_interface_ffproblem.hh"
+
+namespace Dumux {
+namespace Properties {
+
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::StokesOneP>
+{
+ using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, Properties::TTag::DarcyOneP>;
+ using type = Dumux::StokesDarcyCouplingManager<Traits>;
+};
+
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::DarcyOneP>
+{
+ using Traits = StaggeredMultiDomainTraits<Properties::TTag::StokesOneP, Properties::TTag::StokesOneP, 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 = Properties::TTag::StokesOneP;
+ using DarcyTypeTag = Properties::TTag::DarcyOneP;
+
+#if EXNUMBER < 3
+ // try to create a grid (from the given grid file or the input file)
+ // for both sub-domains
+ using DarcyGridManager = Dumux::GridManager<GetPropType<DarcyTypeTag, Properties::Grid>>;
+ DarcyGridManager darcyGridManager;
+ darcyGridManager.init("Darcy"); // pass parameter group
+
+ using StokesGridManager = Dumux::GridManager<GetPropType<StokesTypeTag, Properties::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();
+#else
+ // use dune-subgrid to create the individual grids
+ static constexpr int dim = 2;
+ using HostGrid = Dune::YaspGrid<2, Dune::TensorProductCoordinates<double, dim> >;
+ using HostGridManager = Dumux::GridManager<HostGrid>;
+ HostGridManager hostGridManager;
+ hostGridManager.init();
+ auto& hostGrid = hostGridManager.grid();
+
+ struct Params
+ {
+ double amplitude = getParam<double>("Grid.Amplitude");
+ double baseline = getParam<double>("Grid.Baseline");
+ double offset = getParam<double>("Grid.Offset");
+ double scaling = getParam<double>("Grid.Scaling");
+ };
+
+ Params params;
+
+ auto elementSelectorStokes = [&](const auto& element)
+ {
+ double interface = params.amplitude * std::sin(( element.geometry().center()[0] -params.offset) / params.scaling * 2.0 * M_PI) + params.baseline;
+ return element.geometry().center()[1] > interface;
+ };
+
+ auto elementSelectorDarcy = [&](const auto& element)
+ {
+ double interface = params.amplitude * std::sin(( element.geometry().center()[0] - params.offset) / params.scaling * 2.0 * M_PI) + params.baseline;
+ return element.geometry().center()[1] < interface;
+ };
+
+ // subgrid Pointer
+ auto stokesGridPtr = SubgridGridCreator<HostGrid>::makeGrid(hostGrid, elementSelectorStokes, "Stokes");
+ auto darcyGridPtr = SubgridGridCreator<HostGrid>::makeGrid(hostGrid, elementSelectorDarcy, "Darcy");
+
+ // we compute on the leaf grid view
+ const auto& darcyGridView = darcyGridPtr->leafGridView();
+ const auto& stokesGridView = stokesGridPtr->leafGridView();
+
+#endif
+
+
+ // create the finite volume grid geometry
+ using StokesFVGridGeometry = GetPropType<StokesTypeTag, Properties::FVGridGeometry>;
+ auto stokesFvGridGeometry = std::make_shared<StokesFVGridGeometry>(stokesGridView);
+ stokesFvGridGeometry->update();
+ using DarcyFVGridGeometry = GetPropType<DarcyTypeTag, Properties::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 = GetPropType<StokesTypeTag, Properties::Problem>;
+ auto stokesProblem = std::make_shared<StokesProblem>(stokesFvGridGeometry, couplingManager);
+ using DarcyProblem = GetPropType<DarcyTypeTag, Properties::Problem>;
+ auto darcyProblem = std::make_shared<DarcyProblem>(darcyFvGridGeometry, couplingManager);
+
+ // the solution vector
+ Traits::SolutionVector sol;
+ sol[stokesCellCenterIdx].resize(stokesFvGridGeometry->numCellCenterDofs());
+ sol[stokesFaceIdx].resize(stokesFvGridGeometry->numFaceDofs());
+ sol[darcyIdx].resize(darcyFvGridGeometry->numDofs());
+
+ auto stokesSol = partial(sol, stokesCellCenterIdx, stokesFaceIdx);
+
+ stokesProblem->applyInitialSolution(stokesSol);
+ darcyProblem->applyInitialSolution(sol[darcyIdx]);
+
+ couplingManager->init(stokesProblem, darcyProblem, sol);
+
+ // the grid variables
+ using StokesGridVariables = GetPropType<StokesTypeTag, Properties::GridVariables>;
+ auto stokesGridVariables = std::make_shared<StokesGridVariables>(stokesProblem, stokesFvGridGeometry);
+ stokesGridVariables->init(stokesSol);
+ using DarcyGridVariables = GetPropType<DarcyTypeTag, Properties::GridVariables>;
+ auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyFvGridGeometry);
+ darcyGridVariables->init(sol[darcyIdx]);
+
+ // 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<StokesGridVariables, decltype(stokesSol)> stokesVtkWriter(*stokesGridVariables, stokesSol, stokesName);
+ GetPropType<StokesTypeTag, Properties::IOFields>::initOutputModule(stokesVtkWriter);
+
+#if EXNUMBER >= 2
+ stokesVtkWriter.addField(stokesProblem->getAnalyticalVelocityX(), "analyticalV_x");
+#endif
+
+ stokesVtkWriter.write(0.0);
+
+ VtkOutputModule<DarcyGridVariables, GetPropType<DarcyTypeTag, Properties::SolutionVector>> darcyVtkWriter(*darcyGridVariables, sol[darcyIdx], darcyName);
+ using DarcyVelocityOutput = GetPropType<DarcyTypeTag, Properties::VelocityOutput>;
+ darcyVtkWriter.addVelocityOutput(std::make_shared<DarcyVelocityOutput>(*darcyGridVariables));
+ GetPropType<DarcyTypeTag, Properties::IOFields>::initOutputModule(darcyVtkWriter);
+ darcyVtkWriter.write(0.0);
+
+ // the assembler for a stationary problem
+ using Assembler = 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);
+
+ // the linear solver
+ using LinearSolver = UMFPackBackend;
+ auto linearSolver = std::make_shared<LinearSolver>();
+
+ // the non-linear solver
+ using NewtonSolver = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
+ NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
+
+ // solve the non-linear system
+ nonLinearSolver.solve(sol);
+
+ // write vtk output
+ stokesVtkWriter.write(1.0);
+ darcyVtkWriter.write(1.0);
+
+ ////////////////////////////////////////////////////////////
+ // finalize, print dumux message to say goodbye
+ ////////////////////////////////////////////////////////////
+
+ // print dumux end message
+ if (mpiHelper.rank() == 0)
+ {
+ Parameters::print();
+ DumuxMessage::print(/*firstCall=*/false);
+ }
+
+ return 0;
+} // end main
+catch (Dumux::ParameterException &e)
+{
+ std::cerr << std::endl << e << " ---> Abort!" << std::endl;
+ return 1;
+}
+catch (Dune::DGFException & e)
+{
+ std::cerr << "DGF exception thrown (" << e <<
+ "). Most likely, the DGF file name is wrong "
+ "or the DGF file is corrupted, "
+ "e.g. missing hash at end of file or wrong number (dimensions) of entries."
+ << " ---> Abort!" << std::endl;
+ return 2;
+}
+catch (Dune::Exception &e)
+{
+ std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
+ return 3;
+}
+catch (...)
+{
+ std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
+ return 4;
+}
diff --git a/appl/coupling-ff-pm/solution/monolithic/interface/ex_interface_ffproblem.hh b/appl/coupling-ff-pm/solution/monolithic/interface/ex_interface_ffproblem.hh
new file mode 100644
index 0000000..79e7c8e
--- /dev/null
+++ b/appl/coupling-ff-pm/solution/monolithic/interface/ex_interface_ffproblem.hh
@@ -0,0 +1,367 @@
+// -*- 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 The free flow sub problem
+ */
+#ifndef DUMUX_STOKES_SUBPROBLEM_HH
+#define DUMUX_STOKES_SUBPROBLEM_HH
+
+#include <dune/grid/yaspgrid.hh>
+
+#if EXNUMBER >= 3
+#include <dumux/io/grid/subgridgridcreator.hh>
+#endif
+
+#include <dumux/material/fluidsystems/1pliquid.hh>
+#include <dumux/material/components/simpleh2o.hh>
+
+#include <dumux/freeflow/navierstokes/problem.hh>
+#include <dumux/discretization/staggered/freeflow/properties.hh>
+#include <dumux/freeflow/navierstokes/model.hh>
+
+namespace Dumux
+{
+template <class TypeTag>
+class StokesSubProblem;
+
+namespace Properties
+{
+// Create new type tags
+namespace TTag {
+struct StokesOneP { using InheritsFrom = std::tuple<NavierStokes, StaggeredFreeFlowModel>; };
+} // end namespace TTag
+
+// the fluid system
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::StokesOneP>
+{
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using type = FluidSystems::OnePLiquid<Scalar, Dumux::Components::SimpleH2O<Scalar> > ;
+};
+
+// Set the grid type
+template<class TypeTag>
+struct Grid<TypeTag, TTag::StokesOneP>
+{
+ static constexpr auto dim = 2;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using TensorGrid = Dune::YaspGrid<2, Dune::TensorProductCoordinates<Scalar, dim> >;
+
+#if EXNUMBER < 3 // use "normal" grid
+ using type = TensorGrid;
+#else // use dune-subgrid
+ using HostGrid = TensorGrid;
+ using type = Dune::SubGrid<dim, HostGrid>;
+#endif
+};
+
+// Set the problem property
+template<class TypeTag>
+struct Problem<TypeTag, TTag::StokesOneP> { using type = Dumux::StokesSubProblem<TypeTag> ; };
+
+template<class TypeTag>
+struct EnableFVGridGeometryCache<TypeTag, TTag::StokesOneP> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridFluxVariablesCache<TypeTag, TTag::StokesOneP> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridVolumeVariablesCache<TypeTag, TTag::StokesOneP> { static constexpr bool value = true; };
+}
+
+/*!
+ * \brief The free flow sub problem
+ */
+template <class TypeTag>
+class StokesSubProblem : public NavierStokesProblem<TypeTag>
+{
+ using ParentType = NavierStokesProblem<TypeTag>;
+
+ using GridView = GetPropType<TypeTag, Properties::GridView>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+
+ using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+
+ using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
+
+ using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using Element = typename GridView::template Codim<0>::Entity;
+
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
+ using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+
+ using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
+
+public:
+ StokesSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry, std::shared_ptr<CouplingManager> couplingManager)
+ : ParentType(fvGridGeometry, "Stokes"), eps_(1e-6), couplingManager_(couplingManager)
+ {
+ deltaP_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.PressureDifference");
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief Return the temperature within the domain in [K].
+ *
+ * This problem assumes a temperature of 10 degrees Celsius.
+ */
+ Scalar temperature() const
+ { return 273.15 + 10; } // 10°C
+
+ /*!
+ * \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.dofPosition();
+
+#if EXNUMBER == 0 // flow from top to bottom
+ if(onUpperBoundary_(globalPos))
+ {
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ }
+
+ if (onRightBoundary_(globalPos) || (onLeftBoundary_(globalPos)))
+ {
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ }
+#else // flow flom left to right
+ if(onLeftBoundary_(globalPos) || onRightBoundary_(globalPos))
+ values.setDirichlet(Indices::pressureIdx);
+ else
+ {
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ }
+#endif
+
+ if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
+ {
+ values.setCouplingNeumann(Indices::conti0EqIdx);
+#if EXNUMBER < 3
+ values.setCouplingNeumann(Indices::momentumYBalanceIdx);
+#else
+ //consider orientation of face automatically
+ values.setCouplingNeumann(scvf.directionIndex());
+#endif
+
+#if EXNUMBER < 2
+ values.setDirichlet(Indices::velocityXIdx); // assume no slip on interface
+#elif EXNUMBER == 2
+ // set the Beaver-Joseph-Saffman slip condition for the tangential momentum balance equation
+ values.setBJS(Indices::momentumXBalanceIdx);
+#else
+ // set the Beaver-Joseph-Saffman slip condition for the tangential momentum balance equation,
+ // consider orientation of face automatically
+ values.setBJS(1 - scvf.directionIndex());
+#endif
+ }
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Dirichlet control volume.
+ *
+ * \param globalPos The global position
+ */
+ PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
+ {
+ PrimaryVariables values(0.0);
+ values = initialAtPos(globalPos);
+
+ 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
+ */
+ template<class ElementVolumeVariables, class ElementFaceVariables>
+ NumEqVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFaceVariables& elemFaceVars,
+ const SubControlVolumeFace& scvf) const
+ {
+ NumEqVector values(0.0);
+
+ if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
+ {
+ values[Indices::conti0EqIdx] = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf);
+#if EXNUMBER < 3
+ values[Indices::momentumYBalanceIdx] = couplingManager().couplingData().momentumCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf);
+#else
+ values[scvf.directionIndex()] = couplingManager().couplingData().momentumCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf);
+#endif
+
+ }
+ 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_; }
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * \param globalPos The global position
+ */
+ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+ {
+ PrimaryVariables values(0.0);
+#if EXNUMBER == 0
+ values[Indices::velocityYIdx] = -1e-6 * globalPos[0] * (this->fvGridGeometry().bBoxMax()[0] - globalPos[0]);
+#else
+ // set fixed pressures on the left and right boundary
+ if(onLeftBoundary_(globalPos))
+ values[Indices::pressureIdx] = deltaP_;
+ if(onRightBoundary_(globalPos))
+ values[Indices::pressureIdx] = 0.0;
+#endif
+
+ return values;
+ }
+
+ /*!
+ * \brief Returns the intrinsic permeability of required as input parameter for the Beavers-Joseph-Saffman boundary condition
+ */
+ Scalar permeability(const Element& element, const SubControlVolumeFace& scvf) const
+ {
+ return couplingManager().couplingData().darcyPermeability(element, 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());
+ }
+
+ /*!
+ * \brief calculate the analytical velocity in x direction based on Beavers & Joseph (1967)
+ */
+ void calculateAnalyticalVelocityX() const
+ {
+ analyticalVelocityX_.resize(this->fvGridGeometry().gridView().size(0));
+
+ using std::sqrt;
+ const Scalar dPdX = -deltaP_ / (this->fvGridGeometry().bBoxMax()[0] - this->fvGridGeometry().bBoxMin()[0]);
+ static const Scalar mu = FluidSystem::viscosity(temperature(), 1e5);
+ static const Scalar alpha = getParam<Scalar>("Darcy.SpatialParams.AlphaBeaversJoseph");
+ static const Scalar K = getParam<Scalar>("Darcy.SpatialParams.Permeability");
+ static const Scalar sqrtK = sqrt(K);
+ const Scalar sigma = (this->fvGridGeometry().bBoxMax()[1] - this->fvGridGeometry().bBoxMin()[1])/sqrtK;
+
+ const Scalar uB = -K/(2.0*mu) * ((sigma*sigma + 2.0*alpha*sigma) / (1.0 + alpha*sigma)) * dPdX;
+
+ for (const auto& element : elements(this->fvGridGeometry().gridView()))
+ {
+ const auto eIdx = this->fvGridGeometry().gridView().indexSet().index(element);
+ const Scalar y = element.geometry().center()[1] - this->fvGridGeometry().bBoxMin()[1];
+
+ const Scalar u = uB*(1.0 + alpha/sqrtK*y) + 1.0/(2.0*mu) * (y*y + 2*alpha*y*sqrtK) * dPdX;
+ analyticalVelocityX_[eIdx] = u;
+ }
+ }
+
+ /*!
+ * \brief Get the analytical velocity in x direction
+ */
+ const std::vector<Scalar>& getAnalyticalVelocityX() const
+ {
+ if(analyticalVelocityX_.empty())
+ calculateAnalyticalVelocityX();
+ return analyticalVelocityX_;
+ }
+
+ // \}
+
+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 deltaP_;
+
+ std::shared_ptr<CouplingManager> couplingManager_;
+
+ mutable std::vector<Scalar> analyticalVelocityX_;
+};
+} //end namespace
+
+#endif // DUMUX_STOKES_SUBPROBLEM_HH
diff --git a/appl/coupling-ff-pm/solution/monolithic/interface/ex_interface_pmproblem.hh b/appl/coupling-ff-pm/solution/monolithic/interface/ex_interface_pmproblem.hh
new file mode 100644
index 0000000..9649b0e
--- /dev/null
+++ b/appl/coupling-ff-pm/solution/monolithic/interface/ex_interface_pmproblem.hh
@@ -0,0 +1,266 @@
+// -*- 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 The porous medium flow sub problem
+ */
+#ifndef DUMUX_DARCY_SUBPROBLEM_HH
+#define DUMUX_DARCY_SUBPROBLEM_HH
+
+#include <dune/grid/yaspgrid.hh>
+
+#if EXNUMBER >= 3
+#include <dumux/io/grid/subgridgridcreator.hh>
+#endif
+
+#include <dumux/discretization/cctpfa.hh>
+
+#include <dumux/porousmediumflow/1p/model.hh>
+#include <dumux/porousmediumflow/problem.hh>
+
+#include "../1pspatialparams.hh"
+
+#include <dumux/material/components/simpleh2o.hh>
+#include <dumux/material/fluidsystems/1pliquid.hh>
+
+namespace Dumux
+{
+template <class TypeTag>
+class DarcySubProblem;
+
+namespace Properties
+{
+// Create new type tags
+namespace TTag {
+struct DarcyOneP { using InheritsFrom = std::tuple<OneP, CCTpfaModel>; };
+} // end namespace TTag
+
+// Set the problem property
+template<class TypeTag>
+struct Problem<TypeTag, TTag::DarcyOneP> { using type = Dumux::DarcySubProblem<TypeTag>; };
+
+// the fluid system
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::DarcyOneP>
+{
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using type = FluidSystems::OnePLiquid<Scalar, Dumux::Components::SimpleH2O<Scalar> > ;
+};
+
+// Set the grid type
+template<class TypeTag>
+struct Grid<TypeTag, TTag::DarcyOneP>
+{
+ static constexpr auto dim = 2;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using TensorGrid = Dune::YaspGrid<2, Dune::TensorProductCoordinates<Scalar, dim> >;
+
+#if EXNUMBER < 3 // use "normal" grid
+ using type = TensorGrid;
+#else // use dune-subgrid
+ using HostGrid = TensorGrid;
+ using type = Dune::SubGrid<dim, HostGrid>;
+#endif
+};
+
+template<class TypeTag>
+struct SpatialParams<TypeTag, TTag::DarcyOneP> {
+ using type = OnePSpatialParams<GetPropType<TypeTag, FVGridGeometry>, GetPropType<TypeTag, Scalar>>;
+};
+} // end namespace Properties
+
+/*!
+ * \brief The porous medium flow sub problem
+ */
+template <class TypeTag>
+class DarcySubProblem : public PorousMediumFlowProblem<TypeTag>
+{
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
+ using GridView = GetPropType<TypeTag, Properties::GridView>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
+ using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
+ using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
+ using FVElementGeometry = typename GetPropType<TypeTag, Properties::FVGridGeometry>::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
+
+ using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
+
+public:
+ DarcySubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry,
+ std::shared_ptr<CouplingManager> couplingManager)
+ : ParentType(fvGridGeometry, "Darcy"), eps_(1e-7), couplingManager_(couplingManager)
+ {}
+
+ /*!
+ * \name Simulation steering
+ */
+ // \{
+
+ /*!
+ * \brief Return the temperature within the domain in [K].
+ *
+ */
+ Scalar temperature() const
+ { return 273.15 + 10; } // 10°C
+ // \}
+
+ /*!
+ * \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 EXNUMBER == 0 // flow from top to bottom
+ if (onLowerBoundary_(scvf.center()))
+ values.setAllDirichlet();
+#endif
+
+ if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
+ values.setAllCouplingNeumann();
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Dirichlet control volume.
+ *
+ * \param element The element for which the Dirichlet boundary condition is set
+ * \param scvf The boundary subcontrolvolumeface
+ *
+ * For this method, the \a values parameter stores primary variables.
+ */
+ PrimaryVariables dirichlet(const Element &element, const SubControlVolumeFace &scvf) const
+ {
+ PrimaryVariables values(0.0);
+ values = initial(element);
+
+ 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[Indices::conti0EqIdx] = couplingManager().couplingData().massCouplingCondition(element, 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 initial(const Element &element) const
+ {
+ return PrimaryVariables(0.0);
+ }
+
+ // \}
+
+ //! Set the coupling manager
+ void setCouplingManager(std::shared_ptr<CouplingManager> cm)
+ { couplingManager_ = cm; }
+
+ //! Get the coupling manager
+ const CouplingManager& couplingManager() const
+ { return *couplingManager_; }
+
+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_;
+ std::shared_ptr<CouplingManager> couplingManager_;
+};
+} //end namespace
+
+#endif //DUMUX_DARCY_SUBPROBLEM_HH
diff --git a/appl/coupling-ff-pm/solution/monolithic/interface/params.input b/appl/coupling-ff-pm/solution/monolithic/interface/params.input
new file mode 100644
index 0000000..4117b7f
--- /dev/null
+++ b/appl/coupling-ff-pm/solution/monolithic/interface/params.input
@@ -0,0 +1,46 @@
+ # for dune-subgrid
+ [Grid]
+Positions0 = 0 1
+Positions1 = 0 0.2 0.3 0.65
+Cells0 = 100
+Cells1 = 10 50 18
+Baseline = 0.25 # [m]
+Amplitude = 0.04 # [m]
+Offset = 0.5 # [m]
+Scaling = 0.2 #[m]
+
+[Stokes.Grid]
+Verbosity = true
+Positions0 = 0.0 1.0
+Positions1 = 1.0 2.0
+Cells0 = 20
+Cells1 = 100
+Grading1 = 1
+
+[Darcy.Grid]
+Verbosity = true
+Positions0 = 0.0 1.0
+Positions1 = 0.0 1.0
+Cells0 = 20
+Cells1 = 20
+Grading1 = 1
+
+[Stokes.Problem]
+Name = stokes
+PressureDifference = 1e-9
+EnableInertiaTerms = false
+
+[Darcy.Problem]
+Name = darcy
+
+[Darcy.SpatialParams]
+Permeability = 1e-6 # m^2
+Porosity = 0.4
+AlphaBeaversJoseph = 1.0
+
+[Problem]
+Name = ex_ff-pm-interface
+EnableGravity = false
+
+[Vtk]
+AddVelocity = 1
diff --git a/appl/coupling-ff-pm/solution/monolithic/models/CMakeLists.txt b/appl/coupling-ff-pm/solution/monolithic/models/CMakeLists.txt
new file mode 100644
index 0000000..dc982cb
--- /dev/null
+++ b/appl/coupling-ff-pm/solution/monolithic/models/CMakeLists.txt
@@ -0,0 +1,21 @@
+dune_add_test(NAME orig_ex_models_coupling_ff-pm
+ SOURCES ex_models_coupling_ff-pm.cc
+ COMPILE_DEFINITIONS EXNUMBER=0)
+
+dune_add_test(NAME sol_a_ex_models_coupling_ff-pm
+ SOURCES ex_models_coupling_ff-pm.cc
+ COMPILE_DEFINITIONS EXNUMBER=1)
+
+dune_add_test(NAME sol_b_ex_models_coupling_ff-pm
+ SOURCES ex_models_coupling_ff-pm.cc
+ COMPILE_DEFINITIONS EXNUMBER=2)
+
+dune_add_test(NAME sol_c_ex_models_coupling_ff-pm
+ SOURCES ex_models_coupling_ff-pm.cc
+ COMPILE_DEFINITIONS EXNUMBER=3)
+
+# 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)
+
+# add a symlink for each input file
+add_input_file_links()
diff --git a/appl/coupling-ff-pm/solution/monolithic/models/ex_models_coupling_ff-pm.cc b/appl/coupling-ff-pm/solution/monolithic/models/ex_models_coupling_ff-pm.cc
new file mode 100644
index 0000000..11b7076
--- /dev/null
+++ b/appl/coupling-ff-pm/solution/monolithic/models/ex_models_coupling_ff-pm.cc
@@ -0,0 +1,289 @@
+// -*- 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/partial.hh>
+#include <dumux/linear/seqsolverbackend.hh>
+#include <dumux/assembly/fvassembler.hh>
+#include <dumux/assembly/diffmethod.hh>
+#include <dumux/discretization/method.hh>
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/staggeredvtkoutputmodule.hh>
+#include <dumux/io/grid/gridmanager.hh>
+
+#include <dumux/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 {
+
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::StokesNC>
+{
+ using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, Properties::TTag::DarcyOnePNC>;
+ using type = Dumux::StokesDarcyCouplingManager<Traits>;
+};
+
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::DarcyOnePNC>
+{
+ using Traits = StaggeredMultiDomainTraits<Properties::TTag::StokesNC, Properties::TTag::StokesNC, 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 = Properties::TTag::StokesNC;
+ using DarcyTypeTag = Properties::TTag::DarcyOnePNC;
+
+ // try to create a grid (from the given grid file or the input file)
+ // for both sub-domains
+ using DarcyGridManager = Dumux::GridManager<GetPropType<DarcyTypeTag, Properties::Grid>>;
+ DarcyGridManager darcyGridManager;
+ darcyGridManager.init("Darcy"); // pass parameter group
+
+ using StokesGridManager = Dumux::GridManager<GetPropType<StokesTypeTag, Properties::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 = GetPropType<StokesTypeTag, Properties::FVGridGeometry>;
+ auto stokesFvGridGeometry = std::make_shared<StokesFVGridGeometry>(stokesGridView);
+ stokesFvGridGeometry->update();
+ using DarcyFVGridGeometry = GetPropType<DarcyTypeTag, Properties::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 = GetPropType<StokesTypeTag, Properties::Problem>;
+ auto stokesProblem = std::make_shared<StokesProblem>(stokesFvGridGeometry, couplingManager);
+ using DarcyProblem = GetPropType<DarcyTypeTag, Properties::Problem>;
+ auto darcyProblem = std::make_shared<DarcyProblem>(darcyFvGridGeometry, couplingManager);
+
+ // initialize the fluidsystem (tabulation)
+ GetPropType<StokesTypeTag, Properties::FluidSystem>::init();
+
+ // get some time loop parameters
+ using Scalar = GetPropType<StokesTypeTag, Properties::Scalar>;
+ const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+ const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+ auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+
+ // instantiate time loop
+ auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(0, 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());
+
+ auto stokesSol = partial(sol, stokesCellCenterIdx, stokesFaceIdx);
+
+ stokesProblem->applyInitialSolution(stokesSol);
+ darcyProblem->applyInitialSolution(sol[darcyIdx]);
+
+ auto solOld = sol;
+
+ couplingManager->init(stokesProblem, darcyProblem, sol);
+
+ // the grid variables
+ using StokesGridVariables = GetPropType<StokesTypeTag, Properties::GridVariables>;
+ auto stokesGridVariables = std::make_shared<StokesGridVariables>(stokesProblem, stokesFvGridGeometry);
+ stokesGridVariables->init(stokesSol);
+ using DarcyGridVariables = GetPropType<DarcyTypeTag, Properties::GridVariables>;
+ auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyFvGridGeometry);
+ darcyGridVariables->init(sol[darcyIdx]);
+
+ // 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<StokesGridVariables, decltype(stokesSol)> stokesVtkWriter(*stokesGridVariables, stokesSol, stokesName);
+ GetPropType<StokesTypeTag, Properties::IOFields>::initOutputModule(stokesVtkWriter);
+ stokesVtkWriter.write(0);
+
+ VtkOutputModule<DarcyGridVariables, GetPropType<DarcyTypeTag, Properties::SolutionVector>> darcyVtkWriter(*darcyGridVariables, sol[darcyIdx], darcyName);
+ using DarcyVelocityOutput = GetPropType<DarcyTypeTag, Properties::VelocityOutput>;
+ darcyVtkWriter.addVelocityOutput(std::make_shared<DarcyVelocityOutput>(*darcyGridVariables));
+ GetPropType<DarcyTypeTag, Properties::IOFields>::initOutputModule(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>();
+
+ 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/appl/coupling-ff-pm/solution/monolithic/models/ex_models_ffproblem.hh b/appl/coupling-ff-pm/solution/monolithic/models/ex_models_ffproblem.hh
new file mode 100644
index 0000000..bc978d9
--- /dev/null
+++ b/appl/coupling-ff-pm/solution/monolithic/models/ex_models_ffproblem.hh
@@ -0,0 +1,357 @@
+// -*- 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/1padapter.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
+{
+// Create new type tags
+namespace TTag {
+struct StokesNC { using InheritsFrom = std::tuple<NavierStokesNC, StaggeredFreeFlowModel>; };
+} // end namespace TTag
+
+// Set the grid type
+template<class TypeTag>
+struct Grid<TypeTag, TTag::StokesNC> { using type = Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
+
+// The fluid system
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::StokesNC>
+{
+ using H2OAir = FluidSystems::H2OAir<GetPropType<TypeTag, Properties::Scalar>>;
+ using type = FluidSystems::OnePAdapter<H2OAir, H2OAir::gasPhaseIdx>;
+};
+
+// Do not replace one equation with a total mass balance
+template<class TypeTag>
+struct ReplaceCompEqIdx<TypeTag, TTag::StokesNC> { static constexpr int value = 3; };
+
+// Use formulation based on mass fractions
+template<class TypeTag>
+struct UseMoles<TypeTag, TTag::StokesNC> { static constexpr bool value = true; };
+
+// Set the problem property
+template<class TypeTag>
+struct Problem<TypeTag, TTag::StokesNC> { using type = Dumux::StokesSubProblem<TypeTag> ; };
+
+template<class TypeTag>
+struct EnableFVGridGeometryCache<TypeTag, TTag::StokesNC> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridFluxVariablesCache<TypeTag, TTag::StokesNC> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridVolumeVariablesCache<TypeTag, TTag::StokesNC> { static constexpr bool value = true; };
+}
+
+/*!
+ * \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 = GetPropType<TypeTag, Properties::GridView>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+ using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+ using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
+
+ using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
+ using ElementFaceVariables = typename GetPropType<TypeTag, Properties::GridFaceVariables>::LocalView;
+ using FluidState = GetPropType<TypeTag, Properties::FluidState>;
+
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
+
+ using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
+ using TimeLoopPtr = std::shared_ptr<TimeLoop<Scalar>>;
+
+ static constexpr bool useMoles = GetPropType<TypeTag, Properties::ModelTraits>::useMoles();
+
+ static constexpr auto dim = GetPropType<TypeTag, Properties::ModelTraits>::dim();
+ static constexpr auto transportCompIdx = 1;
+
+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
+ */
+ // \{
+
+ /*!
+ * \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(Indices::conti0EqIdx + 1);
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ }
+ else if(onRightBoundary_(globalPos))
+ {
+ values.setDirichlet(Indices::pressureIdx);
+ values.setOutflow(Indices::conti0EqIdx + 1);
+ }
+ else
+ {
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ values.setNeumann(Indices::conti0EqIdx);
+ values.setNeumann(Indices::conti0EqIdx + 1);
+ }
+
+ if (couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
+ {
+ values.setCouplingNeumann(Indices::conti0EqIdx);
+ values.setCouplingNeumann(Indices::conti0EqIdx + 1);
+ 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(element, fvGeometry, elemVolVars, elemFaceVars, scvf);
+
+ const auto massFlux = couplingManager().couplingData().massCouplingCondition(element, 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
+ {
+ PrimaryVariables values(0.0);
+
+ // This is only an approximation of the actual hydorostatic pressure gradient.
+ // Air is compressible (the density depends on pressure), thus the actual
+ // vertical pressure profile is non-linear.
+ // This discrepancy can lead to spurious flows at the outlet if the inlet
+ // velocity is chosen too small.
+ FluidState fluidState;
+ updateFluidStateForBC_(fluidState, pressure_);
+ const Scalar density = FluidSystem::density(fluidState, 0);
+ values[Indices::pressureIdx] = pressure_ - density*this->gravity()[1]*(this->fvGridGeometry().bBoxMax()[1] - globalPos[1]);
+
+
+ // gravity has negative sign
+ values[Indices::conti0EqIdx + 1] = 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 Element& element, const SubControlVolumeFace& scvf) const
+ {
+ return couplingManager().couplingData().darcyPermeability(element, 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(0, pressure);
+ fluidState.setSaturation(0, 1.0);
+ fluidState.setMoleFraction(0, 1, moleFraction_);
+ fluidState.setMoleFraction(0, 0, 1.0 - moleFraction_);
+
+ typename FluidSystem::ParameterCache paramCache;
+ paramCache.updatePhase(fluidState, 0);
+
+ const Scalar density = FluidSystem::density(fluidState, paramCache, 0);
+ fluidState.setDensity(0, density);
+
+ const Scalar molarDensity = FluidSystem::molarDensity(fluidState, paramCache, 0);
+ fluidState.setMolarDensity(0, molarDensity);
+
+ const Scalar enthalpy = FluidSystem::enthalpy(fluidState, paramCache, 0);
+ fluidState.setEnthalpy(0, 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/appl/coupling-ff-pm/solution/monolithic/models/ex_models_pmproblem.hh b/appl/coupling-ff-pm/solution/monolithic/models/ex_models_pmproblem.hh
new file mode 100644
index 0000000..7377192
--- /dev/null
+++ b/appl/coupling-ff-pm/solution/monolithic/models/ex_models_pmproblem.hh
@@ -0,0 +1,475 @@
+// -*- 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/cctpfa.hh>
+#include <dumux/io/gnuplotinterface.hh>
+#include <dumux/material/fluidsystems/1padapter.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
+{
+// Create new type tags
+namespace TTag {
+#if EXNUMBER >= 1
+struct DarcyOnePNC { using InheritsFrom = std::tuple<TwoPNC, CCTpfaModel>; };
+#else
+struct DarcyOnePNC { using InheritsFrom = std::tuple<OnePNC, CCTpfaModel>; };
+#endif
+} // end namespace TTag
+
+// Set the problem property
+template<class TypeTag>
+struct Problem<TypeTag, TTag::DarcyOnePNC> { using type = Dumux::DarcySubProblem<TypeTag>; };
+
+// The fluid system
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::DarcyOnePNC>
+{
+ using H2OAir = FluidSystems::H2OAir<GetPropType<TypeTag, Properties::Scalar>>;
+#if EXNUMBER == 0
+ using type = FluidSystems::OnePAdapter<H2OAir, H2OAir::gasPhaseIdx>;
+#else
+ using type = H2OAir;
+#endif
+};
+
+// Use moles
+template<class TypeTag>
+struct UseMoles<TypeTag, TTag::DarcyOnePNC> { static constexpr bool value = true; };
+
+// Do not replace one equation with a total mass balance
+template<class TypeTag>
+struct ReplaceCompEqIdx<TypeTag, TTag::DarcyOnePNC> { static constexpr int value = 3; };
+
+//! Use a model with constant tortuosity for the effective diffusivity
+template<class TypeTag>
+struct EffectiveDiffusivityModel<TypeTag, TTag::DarcyOnePNC>
+{ using type = DiffusivityConstantTortuosity<GetPropType<TypeTag, Properties::Scalar>>; };
+// Set the grid type
+template<class TypeTag>
+struct Grid<TypeTag, TTag::DarcyOnePNC> { using type = Dune::YaspGrid<2>; };
+
+#if EXNUMBER >= 1
+//! Set the default formulation to pw-Sn: This can be over written in the problem.
+template<class TypeTag>
+struct Formulation<TypeTag, TTag::DarcyOnePNC>
+{ static constexpr auto value = TwoPFormulation::p1s0; };
+#endif
+
+// Set the spatial paramaters type
+#if EXNUMBER >= 1
+template<class TypeTag>
+struct SpatialParams<TypeTag, TTag::DarcyOnePNC> {
+ using type = TwoPSpatialParams<GetPropType<TypeTag, FVGridGeometry>, GetPropType<TypeTag, Scalar>>;
+};
+#else
+template<class TypeTag>
+struct SpatialParams<TypeTag, TTag::DarcyOnePNC> {
+ using type = OnePSpatialParams<GetPropType<TypeTag, FVGridGeometry>, GetPropType<TypeTag, Scalar>>;
+};
+#endif
+
+} // end namespace Properties
+
+template <class TypeTag>
+class DarcySubProblem : public PorousMediumFlowProblem<TypeTag>
+{
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
+ using GridView = GetPropType<TypeTag, Properties::GridView>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+ using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
+ using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
+ using FVElementGeometry = typename GetPropType<TypeTag, Properties::FVGridGeometry>::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
+
+ // copy some indices for convenience
+ using Indices = typename GetPropType<TypeTag, Properties::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 = 0,
+ transportCompIdx = 1
+#endif
+ };
+
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
+
+ using CouplingManager = GetPropType<TypeTag, Properties::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
+
+ // 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))
+ {
+ for(int phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx)
+ {
+ // insert calculation of the water mass here
+#if EXNUMBER >= 2
+ const auto& volVars = elemVolVars[scv];
+ 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;
+
+#if EXNUMBER >= 2
+ NumEqVector flux = couplingManager().couplingData().massCouplingCondition(element, 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()));
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \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(element, 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
+ {
+ static const Scalar stokesPressure = getParamFromGroup<Scalar>("Stokes", "Problem.Pressure");
+
+ 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] = stokesPressure;
+ 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 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 Dumux
+
+#endif //DUMUX_DARCY_SUBPROBLEM_HH
diff --git a/appl/coupling-ff-pm/solution/monolithic/models/params.input b/appl/coupling-ff-pm/solution/monolithic/models/params.input
new file mode 100644
index 0000000..7dda6b1
--- /dev/null
+++ b/appl/coupling-ff-pm/solution/monolithic/models/params.input
@@ -0,0 +1,50 @@
+[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 = true
+EnableInertiaTerms = true
+MoleFraction = 0.0 # -
+Pressure = 1e5 # Pa
+Velocity = 1 # 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 = true
+PlotStorage = true
+
+[Newton]
+MaxSteps = 12
+MaxRelativeShift = 1e-5
+
+[Vtk]
+AddVelocity = 1
diff --git a/appl/coupling-ff-pm/solution/monolithic/turbulence/CMakeLists.txt b/appl/coupling-ff-pm/solution/monolithic/turbulence/CMakeLists.txt
new file mode 100644
index 0000000..f39c13f
--- /dev/null
+++ b/appl/coupling-ff-pm/solution/monolithic/turbulence/CMakeLists.txt
@@ -0,0 +1,17 @@
+dune_add_test(NAME orig_turbulence_coupling_ff-pm
+ SOURCES ex_turbulence_coupling_ff-pm.cc
+ COMPILE_DEFINITIONS EXNUMBER=0)
+
+dune_add_test(NAME sol_a_ex_turbulence_coupling_ff-pm
+ SOURCES ex_turbulence_coupling_ff-pm.cc
+ COMPILE_DEFINITIONS EXNUMBER=1)
+
+dune_add_test(NAME sol_b_ex_turbulence_coupling_ff-pm
+ SOURCES ex_turbulence_coupling_ff-pm.cc
+ COMPILE_DEFINITIONS EXNUMBER=2)
+
+# add exercise to the common target
+add_dependencies(test_exercises sol_a_ex_turbulence_coupling_ff-pm sol_b_ex_turbulence_coupling_ff-pm sol_b_ex_turbulence_coupling_ff-pm)
+
+# add a symlink for each input file
+add_input_file_links()
diff --git a/appl/coupling-ff-pm/solution/monolithic/turbulence/ex_turbulence_coupling_ff-pm.cc b/appl/coupling-ff-pm/solution/monolithic/turbulence/ex_turbulence_coupling_ff-pm.cc
new file mode 100644
index 0000000..995a760
--- /dev/null
+++ b/appl/coupling-ff-pm/solution/monolithic/turbulence/ex_turbulence_coupling_ff-pm.cc
@@ -0,0 +1,290 @@
+// -*- 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 isothermal coupled Stokes/Darcy problem (1p2c/2p2c)
+ */
+#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/partial.hh>
+#include <dumux/linear/seqsolverbackend.hh>
+#include <dumux/assembly/fvassembler.hh>
+#include <dumux/assembly/diffmethod.hh>
+#include <dumux/discretization/method.hh>
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/staggeredvtkoutputmodule.hh>
+#include <dumux/io/grid/gridmanager.hh>
+
+#include <dumux/multidomain/staggeredtraits.hh>
+#include <dumux/multidomain/fvassembler.hh>
+#include <dumux/multidomain/newtonsolver.hh>
+
+#include <dumux/multidomain/boundary/stokesdarcy/couplingmanager.hh>
+
+#include "ex_turbulence_pmproblem.hh"
+#include "ex_turbulence_ffproblem.hh"
+
+namespace Dumux {
+namespace Properties {
+
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::StokesZeroEq>
+{
+ using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, Properties::TTag::DarcyTwoPTwoCNI>;
+ using type = Dumux::StokesDarcyCouplingManager<Traits>;
+};
+
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::DarcyTwoPTwoCNI>
+{
+ using Traits = StaggeredMultiDomainTraits<Properties::TTag::StokesZeroEq, Properties::TTag::StokesZeroEq, 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 = Properties::TTag::StokesZeroEq;
+ using DarcyTypeTag = Properties::TTag::DarcyTwoPTwoCNI;
+
+ // try to create a grid (from the given grid file or the input file)
+ // for both sub-domains
+ using DarcyGridManager = Dumux::GridManager<GetPropType<DarcyTypeTag, Properties::Grid>>;
+ DarcyGridManager darcyGridManager;
+ darcyGridManager.init("Darcy"); // pass parameter group
+
+ using StokesGridManager = Dumux::GridManager<GetPropType<StokesTypeTag, Properties::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 = GetPropType<StokesTypeTag, Properties::FVGridGeometry>;
+ auto stokesFvGridGeometry = std::make_shared<StokesFVGridGeometry>(stokesGridView);
+ stokesFvGridGeometry->update();
+ using DarcyFVGridGeometry = GetPropType<DarcyTypeTag, Properties::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 = GetPropType<StokesTypeTag, Properties::Problem>;
+ auto stokesProblem = std::make_shared<StokesProblem>(stokesFvGridGeometry, couplingManager);
+ using DarcyProblem = GetPropType<DarcyTypeTag, Properties::Problem>;
+ auto darcyProblem = std::make_shared<DarcyProblem>(darcyFvGridGeometry, couplingManager);
+
+ // initialize the fluidsystem (tabulation)
+ GetPropType<StokesTypeTag, Properties::FluidSystem>::init();
+
+ // get some time loop parameters
+ using Scalar = GetPropType<StokesTypeTag, Properties::Scalar>;
+ const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+ const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+ auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+
+ // instantiate time loop
+ auto timeLoop = std::make_shared<TimeLoop<Scalar>>(0, dt, tEnd);
+ timeLoop->setMaxTimeStepSize(maxDt);
+
+ // set timeloop for the subproblems, needed for boundary value variations
+ 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());
+
+ auto stokesSol = partial(sol, stokesCellCenterIdx, stokesFaceIdx);
+
+ stokesProblem->applyInitialSolution(stokesSol);
+ darcyProblem->applyInitialSolution(sol[darcyIdx]);
+
+ auto solOld = sol;
+
+ couplingManager->init(stokesProblem, darcyProblem, sol);
+ // TODO: update static wall properties
+ // TODO: update dynamic wall properties
+#if EXNUMBER >= 1
+ stokesProblem->updateStaticWallProperties();
+ stokesProblem->updateDynamicWallProperties(stokesSol);
+#endif
+
+ // the grid variables
+ using StokesGridVariables = GetPropType<StokesTypeTag, Properties::GridVariables>;
+ auto stokesGridVariables = std::make_shared<StokesGridVariables>(stokesProblem, stokesFvGridGeometry);
+ stokesGridVariables->init(stokesSol);
+ using DarcyGridVariables = GetPropType<DarcyTypeTag, Properties::GridVariables>;
+ auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyFvGridGeometry);
+ darcyGridVariables->init(sol[darcyIdx]);
+
+ // 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<StokesGridVariables, decltype(stokesSol)> stokesVtkWriter(*stokesGridVariables, stokesSol, stokesName);
+ GetPropType<StokesTypeTag, Properties::IOFields>::initOutputModule(stokesVtkWriter);
+ stokesVtkWriter.write(0.0);
+
+ VtkOutputModule<DarcyGridVariables, GetPropType<DarcyTypeTag, Properties::SolutionVector>> darcyVtkWriter(*darcyGridVariables, sol[darcyIdx], darcyName);
+ using DarcyVelocityOutput = GetPropType<DarcyTypeTag, Properties::VelocityOutput>;
+ darcyVtkWriter.addVelocityOutput(std::make_shared<DarcyVelocityOutput>(*darcyGridVariables));
+ GetPropType<DarcyTypeTag, Properties::IOFields>::initOutputModule(darcyVtkWriter);
+ darcyVtkWriter.write(0.0);
+
+ // the assembler 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 non-linear solver
+ using NewtonSolver = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
+ NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
+
+ // time loop
+ 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;
+
+#if EXNUMBER >= 1
+ // TODO: update dynamic wall properties
+ stokesProblem->updateDynamicWallProperties(stokesSol);
+#endif
+
+ // post time step treatment of Darcy problem
+ darcyProblem->postTimeStep(sol[darcyIdx], *darcyGridVariables, timeLoop->timeStepSize());
+
+ // advance grid variables to the next time step
+ stokesGridVariables->advanceTimeStep();
+ darcyGridVariables->advanceTimeStep();
+
+ // advance to the time loop to the next step
+ timeLoop->advanceTimeStep();
+
+ // write vtk output
+ 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/appl/coupling-ff-pm/solution/monolithic/turbulence/ex_turbulence_ffproblem.hh b/appl/coupling-ff-pm/solution/monolithic/turbulence/ex_turbulence_ffproblem.hh
new file mode 100644
index 0000000..07213ac
--- /dev/null
+++ b/appl/coupling-ff-pm/solution/monolithic/turbulence/ex_turbulence_ffproblem.hh
@@ -0,0 +1,423 @@
+// -*- 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 The free-flow sub problem
+ */
+#ifndef DUMUX_FREEFLOW1P2C_SUBPROBLEM_HH
+#define DUMUX_FREEFLOW1P2C_SUBPROBLEM_HH
+
+#include <dune/grid/yaspgrid.hh>
+
+#include <dumux/material/fluidsystems/1padapter.hh>
+#include <dumux/material/fluidsystems/h2oair.hh>
+#include <dumux/discretization/staggered/freeflow/properties.hh>
+
+#if EXNUMBER >= 1
+#include <dumux/freeflow/compositional/zeroeqncmodel.hh>
+#include <dumux/freeflow/rans/zeroeq/problem.hh>
+#else
+#include <dumux/freeflow/compositional/navierstokesncmodel.hh>
+#include <dumux/freeflow/navierstokes/problem.hh>
+#endif
+
+namespace Dumux {
+
+template <class TypeTag>
+class FreeFlowSubProblem;
+
+namespace Properties {
+
+// Create new type tags
+namespace TTag {
+#if EXNUMBER >= 1
+struct StokesZeroEq { using InheritsFrom = std::tuple<ZeroEqNCNI, StaggeredFreeFlowModel>; };
+#else
+struct StokesZeroEq { using InheritsFrom = std::tuple<NavierStokesNCNI, StaggeredFreeFlowModel>; };
+#endif
+} // end namespace TTag
+
+// Set the grid type
+template<class TypeTag>
+struct Grid<TypeTag, TTag::StokesZeroEq> { using type = Dune::YaspGrid<2, Dune::TensorProductCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
+
+// The fluid system
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::StokesZeroEq>
+{
+ using H2OAir = FluidSystems::H2OAir<GetPropType<TypeTag, Properties::Scalar>>;
+ static constexpr auto phaseIdx = H2OAir::gasPhaseIdx; // simulate the air phase
+ using type = FluidSystems::OnePAdapter<H2OAir, phaseIdx>;
+};
+
+template<class TypeTag>
+struct ReplaceCompEqIdx<TypeTag, TTag::StokesZeroEq> { static constexpr int value = 3; };
+
+// Use formulation based on mass fractions
+template<class TypeTag>
+struct UseMoles<TypeTag, TTag::StokesZeroEq> { static constexpr bool value = true; };
+
+// Set the problem property
+template<class TypeTag>
+struct Problem<TypeTag, TTag::StokesZeroEq> { using type = Dumux::FreeFlowSubProblem<TypeTag> ; };
+
+template<class TypeTag>
+struct EnableFVGridGeometryCache<TypeTag, TTag::StokesZeroEq> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridFluxVariablesCache<TypeTag, TTag::StokesZeroEq> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridVolumeVariablesCache<TypeTag, TTag::StokesZeroEq> { static constexpr bool value = true; };
+}
+
+/*!
+ * \brief The free-flow sub problem
+ */
+template <class TypeTag>
+#if EXNUMBER >= 1
+class FreeFlowSubProblem : public ZeroEqProblem<TypeTag>
+{
+ using ParentType = ZeroEqProblem<TypeTag>;
+#else
+class FreeFlowSubProblem : public NavierStokesProblem<TypeTag>
+{
+ using ParentType = NavierStokesProblem<TypeTag>;
+#endif
+
+ using GridView = GetPropType<TypeTag, Properties::GridView>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+ using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+ using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
+
+ using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
+ using ElementFaceVariables = typename GetPropType<TypeTag, Properties::GridFaceVariables>::LocalView;
+ using FluidState = GetPropType<TypeTag, Properties::FluidState>;
+
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
+
+ using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
+ using TimeLoopPtr = std::shared_ptr<TimeLoop<Scalar>>;
+
+ using DiffusionCoefficientAveragingType = typename StokesDarcyCouplingOptions::DiffusionCoefficientAveragingType;
+
+ static constexpr bool useMoles = GetPropType<TypeTag, Properties::ModelTraits>::useMoles();
+
+public:
+ FreeFlowSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry, std::shared_ptr<CouplingManager> couplingManager)
+ : ParentType(fvGridGeometry, "Stokes"), eps_(1e-6), couplingManager_(couplingManager)
+ {
+ refVelocity_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.RefVelocity");
+ refPressure_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.RefPressure");
+ refMoleFrac_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.refMoleFrac");
+ refTemperature_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.RefTemperature");
+
+ diffCoeffAvgType_ = StokesDarcyCouplingOptions::stringToEnum(DiffusionCoefficientAveragingType{},
+ getParamFromGroup<std::string>(this->paramGroup(), "Problem.InterfaceDiffusionCoefficientAvg"));
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief Return the temperature within the domain in [K].
+ */
+ Scalar temperature() const
+ { return refTemperature_; }
+
+ /*!
+ * \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(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ values.setDirichlet(Indices::conti0EqIdx + 1);
+ values.setDirichlet(Indices::energyEqIdx);
+ }
+
+ if (onLowerBoundary_(globalPos))
+ {
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ values.setNeumann(Indices::conti0EqIdx);
+ values.setNeumann(Indices::conti0EqIdx + 1);
+ values.setNeumann(Indices::energyEqIdx);
+ }
+
+ if (onUpperBoundary_(globalPos))
+ {
+#if EXNUMBER >=2
+ values.setAllSymmetry();
+#else
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ values.setNeumann(Indices::conti0EqIdx);
+ values.setNeumann(Indices::conti0EqIdx + 1);
+ values.setNeumann(Indices::energyEqIdx);
+#endif
+ }
+
+ if (onRightBoundary_(globalPos))
+ {
+ values.setDirichlet(Indices::pressureIdx);
+ values.setOutflow(Indices::conti0EqIdx + 1);
+ values.setOutflow(Indices::energyEqIdx);
+ }
+
+ if (couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
+ {
+ values.setCouplingNeumann(Indices::conti0EqIdx);
+ values.setCouplingNeumann(Indices::conti0EqIdx + 1);
+ values.setCouplingNeumann(Indices::energyEqIdx);
+ 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(element, fvGeometry, elemVolVars, elemFaceVars, scvf);
+
+ const auto massFlux = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf, diffCoeffAvgType_);
+ values[Indices::conti0EqIdx] = massFlux[0];
+ values[Indices::conti0EqIdx + 1] = massFlux[1];
+ values[Indices::energyEqIdx] = couplingManager().couplingData().energyCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf, diffCoeffAvgType_);
+ }
+ 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_; }
+
+#if EXNUMBER >= 2
+ bool isOnWallAtPos(const GlobalPosition& globalPos) const
+ {
+ return (onLowerBoundary_(globalPos));
+ }
+#elif EXNUMBER >= 1
+ bool isOnWallAtPos(const GlobalPosition& globalPos) const
+ {
+ return (onLowerBoundary_(globalPos) || onUpperBoundary_(globalPos));
+ }
+#endif
+
+ /*!
+ * \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, refTemperature(), refPressure(), refMoleFrac());
+
+ const Scalar density = FluidSystem::density(fluidState, 0);
+
+ PrimaryVariables values(0.0);
+ values[Indices::pressureIdx] = refPressure() + density*this->gravity()[1]*(globalPos[1] - this->fvGridGeometry().bBoxMin()[1]);
+ values[Indices::conti0EqIdx + 1] = refMoleFrac();
+ values[Indices::velocityXIdx] = refVelocity();
+ values[Indices::temperatureIdx] = refTemperature();
+
+#if EXNUMBER >= 2
+ if(onLowerBoundary_(globalPos))
+ values[Indices::velocityXIdx] = 0.0;
+#else
+ if(onUpperBoundary_(globalPos) || onLowerBoundary_(globalPos))
+ values[Indices::velocityXIdx] = 0.0;
+#endif
+
+ return values;
+ }
+
+ //! \brief Returns the reference velocity.
+ const Scalar refVelocity() const
+ { return refVelocity_ ;}
+
+ //! \brief Returns the reference pressure.
+ const Scalar refPressure() const
+ { return refPressure_; }
+
+ //! \brief Returns the reference mass fraction.
+ const Scalar refMoleFrac() const
+ { return refMoleFrac_; }
+
+ //! \brief Returns the reference temperature.
+ const Scalar refTemperature() const
+ { return refTemperature_; }
+
+
+ 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 Element& element, const SubControlVolumeFace& scvf) const
+ {
+ return couplingManager().problem(CouplingManager::darcyIdx).spatialParams().permeabilityAtPos(scvf.center());
+ }
+
+ /*!
+ * \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 temperature,
+ const Scalar pressure, const Scalar moleFraction) const
+ {
+ fluidState.setTemperature(temperature);
+ fluidState.setPressure(0, pressure);
+ fluidState.setSaturation(0, 1.0);
+ fluidState.setMoleFraction(0, 1, moleFraction);
+ fluidState.setMoleFraction(0, 0, 1.0 - moleFraction);
+
+ typename FluidSystem::ParameterCache paramCache;
+ paramCache.updatePhase(fluidState, 0);
+
+ const Scalar density = FluidSystem::density(fluidState, paramCache, 0);
+ fluidState.setDensity(0, density);
+
+ const Scalar molarDensity = FluidSystem::molarDensity(fluidState, paramCache, 0);
+ fluidState.setMolarDensity(0, molarDensity);
+
+ const Scalar enthalpy = FluidSystem::enthalpy(fluidState, paramCache, 0);
+ fluidState.setEnthalpy(0, enthalpy);
+ }
+
+ // the height of the free-flow domain
+ const Scalar height_() const
+ { return this->fvGridGeometry().bBoxMax()[1] - this->fvGridGeometry().bBoxMin()[1]; }
+
+ Scalar eps_;
+
+ Scalar refVelocity_;
+ Scalar refPressure_;
+ Scalar refMoleFrac_;
+ Scalar refTemperature_;
+
+ TimeLoopPtr timeLoop_;
+
+ std::shared_ptr<CouplingManager> couplingManager_;
+
+ DiffusionCoefficientAveragingType diffCoeffAvgType_;
+};
+} //end namespace
+
+#endif // DUMUX_STOKES1P2C_SUBPROBLEM_HH
diff --git a/appl/coupling-ff-pm/solution/monolithic/turbulence/ex_turbulence_pmproblem.hh b/appl/coupling-ff-pm/solution/monolithic/turbulence/ex_turbulence_pmproblem.hh
new file mode 100644
index 0000000..0ded6f6
--- /dev/null
+++ b/appl/coupling-ff-pm/solution/monolithic/turbulence/ex_turbulence_pmproblem.hh
@@ -0,0 +1,337 @@
+// -*- 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 The porous medium sub problem
+ */
+#ifndef DUMUX_DARCY2P2C_SUBPROBLEM_HH
+#define DUMUX_DARCY2P2C_SUBPROBLEM_HH
+
+#include <dune/grid/yaspgrid.hh>
+
+#include <dumux/discretization/cctpfa.hh>
+
+#include <dumux/porousmediumflow/2p2c/model.hh>
+#include <dumux/porousmediumflow/problem.hh>
+
+#include <dumux/material/fluidsystems/h2oair.hh>
+
+#include "../2pspatialparams.hh"
+
+namespace Dumux
+{
+template <class TypeTag>
+class DarcySubProblem;
+
+namespace Properties
+{
+// Create new type tags
+namespace TTag {
+struct DarcyTwoPTwoCNI { using InheritsFrom = std::tuple<TwoPTwoCNI, CCTpfaModel>; };
+} // end namespace TTag
+
+// Set the problem property
+template<class TypeTag>
+struct Problem<TypeTag, TTag::DarcyTwoPTwoCNI> { using type = Dumux::DarcySubProblem<TypeTag>; };
+
+// the fluid system
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::DarcyTwoPTwoCNI> { using type = FluidSystems::H2OAir<GetPropType<TypeTag, Properties::Scalar>>; };
+
+//! Set the default formulation to pw-Sn: This can be over written in the problem.
+template<class TypeTag>
+struct Formulation<TypeTag, TTag::DarcyTwoPTwoCNI>
+{ static constexpr auto value = TwoPFormulation::p1s0; };
+
+// The gas component balance (air) is replaced by the total mass balance
+template<class TypeTag>
+struct ReplaceCompEqIdx<TypeTag, TTag::DarcyTwoPTwoCNI> { static constexpr int value = 3; };
+
+// Set the grid type
+template<class TypeTag>
+struct Grid<TypeTag, TTag::DarcyTwoPTwoCNI> { using type = Dune::YaspGrid<2, Dune::TensorProductCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
+
+template<class TypeTag>
+struct UseMoles<TypeTag, TTag::DarcyTwoPTwoCNI> { static constexpr bool value = true; };
+
+template<class TypeTag>
+struct SpatialParams<TypeTag, TTag::DarcyTwoPTwoCNI> {
+ using type = TwoPSpatialParams<GetPropType<TypeTag, FVGridGeometry>, GetPropType<TypeTag, Scalar>>;
+};
+
+} // end namespace Properties
+
+/*!
+ * \brief The porous medium sub problem
+ */
+template <class TypeTag>
+class DarcySubProblem : public PorousMediumFlowProblem<TypeTag>
+{
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
+ using GridView = GetPropType<TypeTag, Properties::GridView>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
+ using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
+ using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
+ using FVElementGeometry = typename GetPropType<TypeTag, Properties::FVGridGeometry>::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
+ using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
+
+ using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+
+ // copy some indices for convenience
+ using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+ enum {
+ // primary variable indices
+ conti0EqIdx = Indices::conti0EqIdx,
+ contiWEqIdx = Indices::conti0EqIdx + FluidSystem::H2OIdx,
+ contiNEqIdx = Indices::conti0EqIdx + FluidSystem::AirIdx,
+ pressureIdx = Indices::pressureIdx,
+ switchIdx = Indices::switchIdx
+ };
+
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
+ using TimeLoopPtr = std::shared_ptr<TimeLoop<Scalar>>;
+
+ using DiffusionCoefficientAveragingType = typename StokesDarcyCouplingOptions::DiffusionCoefficientAveragingType;
+
+public:
+ DarcySubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry,
+ std::shared_ptr<CouplingManager> couplingManager)
+ : ParentType(fvGridGeometry, "Darcy"), eps_(1e-7), couplingManager_(couplingManager)
+ {
+ pressure_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Pressure");
+ initialSw_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Saturation");
+ temperature_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Temperature");
+ initialPhasePresence_ = getParamFromGroup<int>(this->paramGroup(), "Problem.InitPhasePresence");
+
+ diffCoeffAvgType_ = StokesDarcyCouplingOptions::stringToEnum(DiffusionCoefficientAveragingType{},
+ getParamFromGroup<std::string>(this->paramGroup(), "Problem.InterfaceDiffusionCoefficientAvg"));
+ }
+
+ /*!
+ * \name Simulation steering
+ */
+ // \{
+
+ template<class SolutionVector, class GridVariables>
+ void postTimeStep(const SolutionVector& curSol,
+ const GridVariables& gridVariables,
+ const Scalar timeStepSize)
+
+ {
+ // compute the mass in the entire domain
+ Scalar massWater = 0.0;
+
+ // bulk elements
+ 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)
+ {
+ massWater += volVars.massFraction(phaseIdx, FluidSystem::H2OIdx)*volVars.density(phaseIdx)
+ * scv.volume() * volVars.saturation(phaseIdx) * volVars.porosity() * volVars.extrusionFactor();
+ }
+ }
+ }
+
+ std::cout << "mass of water is: " << massWater << std::endl;
+ }
+
+ /*!
+ * \brief Return the temperature within the domain in [K].
+ */
+ Scalar temperature() const
+ { return temperature_; }
+ // \}
+
+ /*!
+ * \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 Dirichlet control volume.
+ *
+ * \param element The element for which the Dirichlet boundary condition is set
+ * \param scvf The boundary subcontrolvolumeface
+ *
+ * For this method, the \a values parameter stores primary variables.
+ */
+ PrimaryVariables dirichlet(const Element &element, const SubControlVolumeFace &scvf) const
+ {
+ PrimaryVariables values(0.0);
+ values = initialAtPos(scvf.center());
+
+ 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
+ *
+ */
+ 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))
+ {
+ const auto massFlux = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, scvf, diffCoeffAvgType_);
+
+ for(int i = 0; i< massFlux.size(); ++i)
+ values[i] = massFlux[i];
+
+ values[Indices::energyEqIdx] = couplingManager().couplingData().energyCouplingCondition(element, fvGeometry, elemVolVars, scvf, diffCoeffAvgType_);
+ }
+
+ 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
+ *
+ * For this method, the \a values variable stores the rate mass
+ * of a component is generated or annihilate per volume
+ * unit. Positive values mean that mass is created, negative ones
+ * mean that it vanishes.
+ */
+ 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.
+ *
+ * For this method, the \a priVars parameter stores primary
+ * variables.
+ */
+ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+ {
+ PrimaryVariables values(0.0);
+ values.setState(initialPhasePresence_);
+
+ values[pressureIdx] = pressure_ + 1. * this->gravity()[1] * (globalPos[1] - this->fvGridGeometry().bBoxMax()[1]);
+ values[switchIdx] = initialSw_;
+ values[Indices::temperatureIdx] = temperature_;
+
+ 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_; }
+
+ 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 pressure_;
+ Scalar initialSw_;
+ Scalar temperature_;
+ int initialPhasePresence_;
+
+ TimeLoopPtr timeLoop_;
+
+ Scalar eps_;
+
+ std::shared_ptr<CouplingManager> couplingManager_;
+ DiffusionCoefficientAveragingType diffCoeffAvgType_;
+};
+} //end namespace
+
+#endif //DUMUX_DARCY2P2C_SUBPROBLEM_HH
diff --git a/appl/coupling-ff-pm/solution/monolithic/turbulence/params.input b/appl/coupling-ff-pm/solution/monolithic/turbulence/params.input
new file mode 100644
index 0000000..00e0246
--- /dev/null
+++ b/appl/coupling-ff-pm/solution/monolithic/turbulence/params.input
@@ -0,0 +1,68 @@
+[TimeLoop]
+DtInitial = 1e-1 # [s]
+MaxTimeStepSize = 43200 # [s] (12 hours)
+TEnd = 864000 # [s] (6 days)
+
+[Stokes.Grid]
+Positions0 = 0.0 0.25
+Positions1 = 0.25 0.5
+Grading0 = 1.0
+Grading1 = 1.0
+Cells0 = 15
+Cells1 = 20
+Verbosity = true
+
+[Darcy.Grid]
+Positions0 = 0.0 0.25
+Positions1 = 0.0 0.25
+Cells0 = 15
+Cells1 = 10
+Grading0 = 1.0
+Grading1 = 1.0
+Verbosity = true
+
+[Stokes.Problem]
+Name = stokes
+RefVelocity = 3.5 # [m/s]
+RefPressure = 1e5 # [Pa]
+refMoleFrac = 0 # [-]
+RefTemperature = 298.15 # [K]
+EnableInertiaTerms = true
+
+[Darcy.Problem]
+Name = darcy
+Pressure = 1e5
+Saturation = 0.5 # initial Sw
+Temperature = 298.15 # [K]
+InitPhasePresence = 3 # bothPhases
+
+[Darcy.SpatialParams]
+Porosity = 0.41
+Permeability = 2.65e-10
+AlphaBeaversJoseph = 1.0
+Swr = 0.005
+Snr = 0.01
+VgAlpha = 6.371e-4
+VgN = 6.9
+
+[Problem]
+Name = ex_coupling_turbulence_ff-pm
+EnableGravity = true
+InterfaceDiffusionCoefficientAvg = Harmonic
+
+[Vtk]
+AddVelocity = true
+WriteFaceData = false
+
+[Newton]
+MaxSteps = 12
+MaxRelativeShift = 1e-5
+
+[Assembly]
+NumericDifferenceMethod = 0
+NumericDifference.BaseEpsilon = 1e-8
+
+[Component]
+SolidDensity = 2700
+SolidThermalConductivity = 2.8
+SolidHeatCapacity = 790
--
GitLab