diff --git a/exercises/exercise-coupling-ff-pm/CMakeLists.txt b/exercises/exercise-coupling-ff-pm/CMakeLists.txt
index 8c868cd47010970d5aae114671048ae0bc2c17f1..bae4c263a6bf7ae6a06d90af66a2e45d224f9b2e 100644
--- a/exercises/exercise-coupling-ff-pm/CMakeLists.txt
+++ b/exercises/exercise-coupling-ff-pm/CMakeLists.txt
@@ -1,2 +1,3 @@
add_subdirectory(interface)
add_subdirectory(models)
+add_subdirectory(turbulence)
diff --git a/exercises/exercise-coupling-ff-pm/README.md b/exercises/exercise-coupling-ff-pm/README.md
index a1ea80399dc4fe426defda0d70d2f3ae2df9d3ef..b325d29033122ebe5010cac7d2e96b5fc3d20683 100644
--- a/exercises/exercise-coupling-ff-pm/README.md
+++ b/exercises/exercise-coupling-ff-pm/README.md
@@ -206,7 +206,7 @@ In the first task, the porous-medium model will be changed from a 1p2c system to
Although a 2p2c system is plugged in, we still want to simulate the same situation as before.
The following changes have to be made in the porous-medium model (`ex_models_pmproblem.hh`):
* Include the 2pnc model: include the respective headers and inherit from the new model `TwoPNC`
-* Exchange the spatial parameters for the 1-phase system by those for a 2-phase system (hint: two occurrences).
+* Exchange the spatial parameters for the 1-phase system by those for a 2-phase system (hint: two occurrences).
One big difference between the 1p and 2p model is, that the primary variables for which
the problem is solved, are not fixed.
@@ -285,7 +285,7 @@ The regularization has a strong influence on how long liquid water is present at
The porous-medium model can now reach a liquid saturation of zero. As already explained above,
for this case the liquid saturation cannot serve as primary variable anymore. However,
manually adapting the primary variable states and values is inconvenient.
-[Class et al. (2007)](http://dx.doi.org/10.1016/S0309-1708(02)00014-3)
+[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.
- Replace the current implementation of the Newton solver with the version which can handle
primary variable switches: `dumux/multidomain/privarswitchnewtonsolver.hh`.
@@ -295,4 +295,95 @@ describe an algorithm to switch the primary variables, if phases should appear o
Now you are able to simulate a complete drying of the porous medium.
-### 4. ...
+### 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, and take the inertial forces into account.
+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 last entry in the `NEW_TYPE_TAG` 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.
+
+Additionally, you have to update the static and dynamic wall properties with
+```c++
+stokesProblem->updateStaticWallProperties();
+stokesProblem->updateDynamicWallProperties(stokesSol);
+```
+after applying the initial solution to the Stokes problem in the main file.
+
+In order to update these properties, the Stokes problem file needs to be extended to provide an isOnWall() method:
+```c++
+bool isOnWall(const GlobalPosition& globalPos) const
+{
+ return (onLowerBoundary_(globalPos) || onUpperBoundary_(globalPos));
+}
+```
+
+Since the dynamic wall properties can change during the simulation and depend on the solution, they have to be updated in each time step.
+Include the following lines into your main file (after `// update dynamic wall properties`):
+```c++
+stokesProblem->updateDynamicWallProperties(stokesSol);
+```
+
+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 choose turbulence model specific quantities (e.g. the turbulent viscosity `nu_t`) for the free flow domain.
+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 `isOnWall(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/exercises/exercise-coupling-ff-pm/turbulence/CMakeLists.txt b/exercises/exercise-coupling-ff-pm/turbulence/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..3fb93b14e46940ca5110f1b2ab3f69143b089593
--- /dev/null
+++ b/exercises/exercise-coupling-ff-pm/turbulence/CMakeLists.txt
@@ -0,0 +1,9 @@
+add_input_file_links()
+
+# executables for ex_interface_coupling_ff-pm
+dune_add_test(NAME ex_turbulence_coupling_ff-pm
+ SOURCES ex_turbulence_coupling_ff-pm.cc
+ CMD_ARGS ex_turbulence_coupling_ff-pm.input)
+
+# add tutorial to the common target
+add_dependencies(test_exercises ex_turbulence_coupling_ff-pm)
diff --git a/exercises/exercise-coupling-ff-pm/turbulence/ex_turbulence_coupling_ff-pm.cc b/exercises/exercise-coupling-ff-pm/turbulence/ex_turbulence_coupling_ff-pm.cc
new file mode 100644
index 0000000000000000000000000000000000000000..80a5b9be520e8573a020b083c213fa95c2dfee9e
--- /dev/null
+++ b/exercises/exercise-coupling-ff-pm/turbulence/ex_turbulence_coupling_ff-pm.cc
@@ -0,0 +1,300 @@
+// -*- 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/geometry/diameter.hh>
+#include <dumux/linear/seqsolverbackend.hh>
+#include <dumux/assembly/fvassembler.hh>
+#include <dumux/assembly/diffmethod.hh>
+#include <dumux/discretization/methods.hh>
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/staggeredvtkoutputmodule.hh>
+#include <dumux/io/grid/gridmanager.hh>
+
+#include <dumux/multidomain/staggeredtraits.hh>
+#include <dumux/multidomain/fvassembler.hh>
+#include <dumux/multidomain/privarswitchnewtonsolver.hh>
+
+#include <dumux/multidomain/boundary/stokesdarcy/couplingmanager.hh>
+
+#include "ex_turbulence_pmproblem.hh"
+#include "ex_turbulence_ffproblem.hh"
+
+namespace Dumux {
+namespace Properties {
+
+SET_PROP(ZeroEqTypeTag, CouplingManager)
+{
+ using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, TTAG(DarcyTwoPTwoCTypeTag)>;
+ using type = Dumux::StokesDarcyCouplingManager<Traits>;
+};
+
+SET_PROP(DarcyTwoPTwoCTypeTag, CouplingManager)
+{
+ using Traits = StaggeredMultiDomainTraits<TTAG(ZeroEqTypeTag), TTAG(ZeroEqTypeTag), TypeTag>;
+ using type = Dumux::StokesDarcyCouplingManager<Traits>;
+};
+
+} // end namespace Properties
+} // end namespace Dumux
+
+int main(int argc, char** argv) try
+{
+ using namespace Dumux;
+
+ // initialize MPI, finalize is done automatically on exit
+ const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+ // print dumux start message
+ if (mpiHelper.rank() == 0)
+ DumuxMessage::print(/*firstCall=*/true);
+
+ // parse command line arguments and input file
+ Parameters::init(argc, argv);
+
+ // Define the sub problem type tags
+ using StokesTypeTag = TTAG(ZeroEqTypeTag);
+ using DarcyTypeTag = TTAG(DarcyTwoPTwoCTypeTag);
+
+ // try to create a grid (from the given grid file or the input file)
+ // for both sub-domains
+ using DarcyGridManager = Dumux::GridManager<typename GET_PROP_TYPE(DarcyTypeTag, Grid)>;
+ DarcyGridManager darcyGridManager;
+ darcyGridManager.init("Darcy"); // pass parameter group
+
+ using StokesGridManager = Dumux::GridManager<typename GET_PROP_TYPE(StokesTypeTag, Grid)>;
+ StokesGridManager stokesGridManager;
+ stokesGridManager.init("Stokes"); // pass parameter group
+
+ // we compute on the leaf grid view
+ const auto& darcyGridView = darcyGridManager.grid().leafGridView();
+ const auto& stokesGridView = stokesGridManager.grid().leafGridView();
+
+ // create the finite volume grid geometry
+ using StokesFVGridGeometry = typename GET_PROP_TYPE(StokesTypeTag, FVGridGeometry);
+ auto stokesFvGridGeometry = std::make_shared<StokesFVGridGeometry>(stokesGridView);
+ stokesFvGridGeometry->update();
+ using DarcyFVGridGeometry = typename GET_PROP_TYPE(DarcyTypeTag, FVGridGeometry);
+ auto darcyFvGridGeometry = std::make_shared<DarcyFVGridGeometry>(darcyGridView);
+ darcyFvGridGeometry->update();
+
+ using Traits = StaggeredMultiDomainTraits<StokesTypeTag, StokesTypeTag, DarcyTypeTag>;
+
+ // the coupling manager
+ using CouplingManager = StokesDarcyCouplingManager<Traits>;
+ auto couplingManager = std::make_shared<CouplingManager>(stokesFvGridGeometry, darcyFvGridGeometry);
+
+ // the indices
+ constexpr auto stokesCellCenterIdx = CouplingManager::stokesCellCenterIdx;
+ constexpr auto stokesFaceIdx = CouplingManager::stokesFaceIdx;
+ constexpr auto darcyIdx = CouplingManager::darcyIdx;
+
+ // the problem (initial and boundary conditions)
+ using StokesProblem = typename GET_PROP_TYPE(StokesTypeTag, Problem);
+ auto stokesProblem = std::make_shared<StokesProblem>(stokesFvGridGeometry, couplingManager);
+ using DarcyProblem = typename GET_PROP_TYPE(DarcyTypeTag, Problem);
+ auto darcyProblem = std::make_shared<DarcyProblem>(darcyFvGridGeometry, couplingManager);
+
+ // initialize the fluidsystem (tabulation)
+ GET_PROP_TYPE(StokesTypeTag, FluidSystem)::init();
+
+ // get some time loop parameters
+ using Scalar = typename GET_PROP_TYPE(StokesTypeTag, Scalar);
+ const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+ const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+ auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+
+ // check if we are about to restart a previously interrupted simulation
+ Scalar restartTime = 0;
+ if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart"))
+ restartTime = getParam<Scalar>("TimeLoop.Restart");
+
+ // instantiate time loop
+ auto timeLoop = std::make_shared<TimeLoop<Scalar>>(restartTime, 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());
+
+ // apply initial solution for instationary problems
+ // auxiliary free flow solution vector
+ typename GET_PROP_TYPE(StokesTypeTag, SolutionVector) stokesSol;
+ stokesSol[stokesCellCenterIdx].resize(sol[stokesCellCenterIdx].size());
+ stokesSol[stokesFaceIdx].resize(sol[stokesFaceIdx].size());
+ stokesProblem->applyInitialSolution(stokesSol);
+ auto solStokesOld = stokesSol;
+ sol[stokesCellCenterIdx] = stokesSol[stokesCellCenterIdx];
+ sol[stokesFaceIdx] = stokesSol[stokesFaceIdx];
+
+ // TODO: update static wall properties
+ // TODO: update dynamic wall properties
+
+ darcyProblem->applyInitialSolution(sol[darcyIdx]);
+ auto solDarcyOld = sol[darcyIdx];
+
+ auto solOld = sol;
+
+ couplingManager->init(stokesProblem, darcyProblem, sol);
+
+ // the grid variables
+ using StokesGridVariables = typename GET_PROP_TYPE(StokesTypeTag, GridVariables);
+ auto stokesGridVariables = std::make_shared<StokesGridVariables>(stokesProblem, stokesFvGridGeometry);
+ stokesGridVariables->init(stokesSol, solStokesOld);
+ using DarcyGridVariables = typename GET_PROP_TYPE(DarcyTypeTag, GridVariables);
+ auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyFvGridGeometry);
+ darcyGridVariables->init(sol[darcyIdx], solDarcyOld);
+
+ // intialize the vtk output module
+ const auto stokesName = getParam<std::string>("Problem.Name") + "_" + stokesProblem->name();
+ const auto darcyName = getParam<std::string>("Problem.Name") + "_" + darcyProblem->name();
+
+ StaggeredVtkOutputModule<StokesTypeTag, GET_PROP_VALUE(StokesTypeTag, PhaseIdx)> stokesVtkWriter(*stokesProblem, *stokesFvGridGeometry, *stokesGridVariables, stokesSol, stokesName);
+ GET_PROP_TYPE(StokesTypeTag, VtkOutputFields)::init(stokesVtkWriter);
+ stokesVtkWriter.write(0.0);
+
+ VtkOutputModule<DarcyTypeTag> darcyVtkWriter(*darcyProblem, *darcyFvGridGeometry, *darcyGridVariables, sol[darcyIdx], darcyName);
+ GET_PROP_TYPE(DarcyTypeTag, VtkOutputFields)::init(darcyVtkWriter);
+ darcyVtkWriter.write(0.0);
+
+ // the assembler with time loop for instationary problem
+ using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
+ auto assembler = std::make_shared<Assembler>(std::make_tuple(stokesProblem, stokesProblem, darcyProblem),
+ std::make_tuple(stokesFvGridGeometry->cellCenterFVGridGeometryPtr(),
+ stokesFvGridGeometry->faceFVGridGeometryPtr(),
+ darcyFvGridGeometry),
+ std::make_tuple(stokesGridVariables->cellCenterGridVariablesPtr(),
+ stokesGridVariables->faceGridVariablesPtr(),
+ darcyGridVariables),
+ couplingManager,
+ timeLoop);
+
+ // the linear solver
+ using LinearSolver = UMFPackBackend;
+ auto linearSolver = std::make_shared<LinearSolver>();
+
+ // the primary variable switches used by the sub models
+ using PriVarSwitchTuple = std::tuple<NoPrimaryVariableSwitch, NoPrimaryVariableSwitch, typename GET_PROP_TYPE(DarcyTypeTag, PrimaryVariableSwitch)>;
+
+ // the non-linear solver
+ using NewtonSolver = MultiDomainPriVarSwitchNewtonSolver<Assembler, LinearSolver, CouplingManager, PriVarSwitchTuple>;
+ 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;
+
+ // update the auxiliary free flow solution vector
+ stokesSol[stokesCellCenterIdx] = sol[stokesCellCenterIdx];
+ stokesSol[stokesFaceIdx] = sol[stokesFaceIdx];
+
+ // 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/exercises/exercise-coupling-ff-pm/turbulence/ex_turbulence_coupling_ff-pm.input b/exercises/exercise-coupling-ff-pm/turbulence/ex_turbulence_coupling_ff-pm.input
new file mode 100644
index 0000000000000000000000000000000000000000..514cf156d55211777c2ec0980f37d1a95124867e
--- /dev/null
+++ b/exercises/exercise-coupling-ff-pm/turbulence/ex_turbulence_coupling_ff-pm.input
@@ -0,0 +1,62 @@
+[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]
+
+[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
diff --git a/exercises/exercise-coupling-ff-pm/turbulence/ex_turbulence_ffproblem.hh b/exercises/exercise-coupling-ff-pm/turbulence/ex_turbulence_ffproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..117f518e8ac6e7f87c772f79752f971d3bf5dc8e
--- /dev/null
+++ b/exercises/exercise-coupling-ff-pm/turbulence/ex_turbulence_ffproblem.hh
@@ -0,0 +1,378 @@
+// -*- 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/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
+{
+NEW_TYPE_TAG(ZeroEqTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, NavierStokesNCNI));
+
+// Set the grid type
+SET_TYPE_PROP(ZeroEqTypeTag, Grid, Dune::YaspGrid<2, Dune::TensorProductCoordinates<typename GET_PROP_TYPE(TypeTag, Scalar), 2> >);
+
+// set the fluid system
+SET_TYPE_PROP(ZeroEqTypeTag, FluidSystem, FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>);
+
+// set phase index (air)
+SET_INT_PROP(ZeroEqTypeTag, PhaseIdx, GET_PROP_TYPE(TypeTag, FluidSystem)::gasPhaseIdx);
+
+SET_INT_PROP(ZeroEqTypeTag, ReplaceCompEqIdx, 3);
+
+// Use formulation based on mass fractions
+SET_BOOL_PROP(ZeroEqTypeTag, UseMoles, true);
+
+// Set the problem property
+SET_TYPE_PROP(ZeroEqTypeTag, Problem, Dumux::FreeFlowSubProblem<TypeTag> );
+
+SET_BOOL_PROP(ZeroEqTypeTag, EnableFVGridGeometryCache, true);
+SET_BOOL_PROP(ZeroEqTypeTag, EnableGridFluxVariablesCache, true);
+SET_BOOL_PROP(ZeroEqTypeTag, EnableGridVolumeVariablesCache, true);
+
+SET_BOOL_PROP(ZeroEqTypeTag, EnableInertiaTerms, true);
+}
+
+/*!
+ * \brief The free-flow sub problem
+ */
+template <class TypeTag>
+class FreeFlowSubProblem : public NavierStokesProblem<TypeTag>
+{
+ using ParentType = NavierStokesProblem<TypeTag>;
+
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
+ using ElementFaceVariables = typename GET_PROP_TYPE(TypeTag, GridFaceVariables)::LocalView;
+ using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
+
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+
+ using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager);
+ using TimeLoopPtr = std::shared_ptr<TimeLoop<Scalar>>;
+
+ using DiffusionCoefficientAveragingType = typename StokesDarcyCouplingOptions::DiffusionCoefficientAveragingType;
+
+ static constexpr bool useMoles = GET_PROP_TYPE(TypeTag, ModelTraits)::useMoles();
+
+ static constexpr auto phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
+ static constexpr auto transportCompIdx = 1 - phaseIdx;
+
+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);
+ values.setDirichlet(Indices::energyBalanceIdx);
+ }
+
+ if (onLowerBoundary_(globalPos))
+ {
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ values.setNeumann(Indices::conti0EqIdx);
+ values.setNeumann(Indices::conti0EqIdx + 1);
+ values.setNeumann(Indices::energyBalanceIdx);
+ }
+
+ if (onUpperBoundary_(globalPos))
+ {
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ values.setNeumann(Indices::conti0EqIdx);
+ values.setNeumann(Indices::conti0EqIdx + 1);
+ values.setNeumann(Indices::energyBalanceIdx);
+ }
+
+ if (onRightBoundary_(globalPos))
+ {
+ values.setDirichlet(Indices::pressureIdx);
+ values.setOutflow(Indices::conti0EqIdx);
+ values.setOutflow(Indices::energyBalanceIdx);
+ }
+
+ if (couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
+ {
+ values.setCouplingNeumann(Indices::conti0EqIdx);
+ values.setCouplingNeumann(Indices::conti0EqIdx + 1);
+ values.setCouplingNeumann(Indices::energyBalanceIdx);
+ values.setCouplingNeumann(Indices::momentumYBalanceIdx);
+ values.setBJS(Indices::momentumXBalanceIdx);
+ }
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Dirichlet control volume.
+ *
+ * \param element The element
+ * \param scvf The subcontrolvolume face
+ */
+ PrimaryVariables dirichletAtPos(const GlobalPosition& pos) const
+ {
+ PrimaryVariables values(0.0);
+ values = initialAtPos(pos);
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Neumann control volume.
+ *
+ * \param element The element for which the Neumann boundary condition is set
+ * \param fvGeomentry The fvGeometry
+ * \param elemVolVars The element volume variables
+ * \param elemFaceVars The element face variables
+ * \param scvf The boundary sub control volume face
+ */
+ NumEqVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFaceVariables& elemFaceVars,
+ const SubControlVolumeFace& scvf) const
+ {
+ PrimaryVariables values(0.0);
+ if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
+ {
+ values[Indices::momentumYBalanceIdx] = couplingManager().couplingData().momentumCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf);
+
+ const auto massFlux = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf, diffCoeffAvgType_);
+ values[Indices::conti0EqIdx] = massFlux[0];
+ values[Indices::conti0EqIdx + 1] = massFlux[1];
+ values[Indices::energyBalanceIdx] = couplingManager().couplingData().energyCouplingCondition(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, phaseIdx);
+
+ PrimaryVariables values(0.0);
+ values[Indices::pressureIdx] = refPressure() + density*this->gravity()[1]*(globalPos[1] - this->fvGridGeometry().bBoxMin()[1]);
+ values[Indices::conti0EqIdx] = refMoleFrac();
+ values[Indices::velocityXIdx] = refVelocity();
+ values[Indices::temperatureIdx] = refTemperature();
+
+ if(onUpperBoundary_(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 SubControlVolumeFace& scvf) const
+ {
+ return couplingManager().couplingData().darcyPermeability(scvf);
+ }
+
+ /*!
+ * \brief Returns the alpha value required as input parameter for the Beavers-Joseph-Saffman boundary condition
+ */
+ Scalar alphaBJ(const SubControlVolumeFace& scvf) const
+ {
+ return couplingManager().problem(CouplingManager::darcyIdx).spatialParams().beaversJosephCoeffAtPos(scvf.center());
+ }
+
+ // \}
+
+private:
+ bool onLeftBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_; }
+
+ bool onRightBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_; }
+
+ bool onLowerBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] < this->fvGridGeometry().bBoxMin()[1] + eps_; }
+
+ bool onUpperBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_; }
+
+ //! \brief updates the fluid state to obtain required quantities for IC/BC
+ void updateFluidStateForBC_(FluidState& fluidState, const Scalar temperature,
+ const Scalar pressure, const Scalar moleFraction) const
+ {
+ fluidState.setTemperature(temperature);
+ fluidState.setPressure(phaseIdx, pressure);
+ fluidState.setSaturation(phaseIdx, 1.0);
+ fluidState.setMoleFraction(phaseIdx, transportCompIdx, moleFraction);
+ fluidState.setMoleFraction(phaseIdx, phaseIdx, 1.0 - moleFraction);
+
+ typename FluidSystem::ParameterCache paramCache;
+ paramCache.updatePhase(fluidState, phaseIdx);
+
+ const Scalar density = FluidSystem::density(fluidState, paramCache, phaseIdx);
+ fluidState.setDensity(phaseIdx, density);
+
+ const Scalar molarDensity = FluidSystem::molarDensity(fluidState, paramCache, phaseIdx);
+ fluidState.setMolarDensity(phaseIdx, molarDensity);
+
+ const Scalar enthalpy = FluidSystem::enthalpy(fluidState, paramCache, phaseIdx);
+ fluidState.setEnthalpy(phaseIdx, enthalpy);
+ }
+
+ // the height of the free-flow domain
+ const Scalar height_() const
+ { return this->fvGridGeometry().bBoxMax()[1] - this->fvGridGeometry().bBoxMin()[1]; }
+
+ Scalar eps_;
+
+ Scalar 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/exercises/exercise-coupling-ff-pm/turbulence/ex_turbulence_pmproblem.hh b/exercises/exercise-coupling-ff-pm/turbulence/ex_turbulence_pmproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..2d3a84d25c4884441720cc2ad10b0ccb8ab7a4ab
--- /dev/null
+++ b/exercises/exercise-coupling-ff-pm/turbulence/ex_turbulence_pmproblem.hh
@@ -0,0 +1,324 @@
+// -*- 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/cellcentered/tpfa/properties.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
+{
+NEW_TYPE_TAG(DarcyTwoPTwoCTypeTag, INHERITS_FROM(CCTpfaModel, TwoPTwoCNI));
+
+// Set the problem property
+SET_TYPE_PROP(DarcyTwoPTwoCTypeTag, Problem, Dumux::DarcySubProblem<TypeTag>);
+
+// the fluid system
+SET_TYPE_PROP(DarcyTwoPTwoCTypeTag, FluidSystem, FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>);
+
+//! Set the default formulation to pw-Sn: This can be over written in the problem.
+SET_PROP(DarcyTwoPTwoCTypeTag, Formulation)
+{ static constexpr auto value = TwoPFormulation::p1s0; };
+
+// The gas component balance (air) is replaced by the total mass balance
+SET_INT_PROP(DarcyTwoPTwoCTypeTag, ReplaceCompEqIdx, 3);
+
+// Set the grid type
+SET_TYPE_PROP(DarcyTwoPTwoCTypeTag, Grid, Dune::YaspGrid<2, Dune::TensorProductCoordinates<typename GET_PROP_TYPE(TypeTag, Scalar), 2> >);
+
+SET_BOOL_PROP(DarcyTwoPTwoCTypeTag, UseMoles, true);
+
+SET_TYPE_PROP(DarcyTwoPTwoCTypeTag, SpatialParams, TwoPSpatialParams<TypeTag>);
+}
+
+/*!
+ * \brief The porous medium sub problem
+ */
+template <class TypeTag>
+class DarcySubProblem : public PorousMediumFlowProblem<TypeTag>
+{
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
+
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+
+ // copy some indices for convenience
+ using Indices = typename GET_PROP_TYPE(TypeTag, 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 = typename GET_PROP_TYPE(TypeTag, 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(fvGeometry, elemVolVars, scvf, diffCoeffAvgType_);
+
+ for(int i = 0; i< massFlux.size(); ++i)
+ values[i] = massFlux[i];
+
+ values[Indices::energyEqIdx] = couplingManager().couplingData().energyCouplingCondition(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/exercises/extradoc/ex_ff-pm-turb_viscosity.png b/exercises/extradoc/ex_ff-pm-turb_viscosity.png
new file mode 100644
index 0000000000000000000000000000000000000000..1e7ddab16db7d8a718490ffbf5c2a57ddb3aeacf
Binary files /dev/null and b/exercises/extradoc/ex_ff-pm-turb_viscosity.png differ
diff --git a/exercises/solution/exercise-coupling-ff-pm/CMakeLists.txt b/exercises/solution/exercise-coupling-ff-pm/CMakeLists.txt
index 8c868cd47010970d5aae114671048ae0bc2c17f1..bae4c263a6bf7ae6a06d90af66a2e45d224f9b2e 100644
--- a/exercises/solution/exercise-coupling-ff-pm/CMakeLists.txt
+++ b/exercises/solution/exercise-coupling-ff-pm/CMakeLists.txt
@@ -1,2 +1,3 @@
add_subdirectory(interface)
add_subdirectory(models)
+add_subdirectory(turbulence)
diff --git a/exercises/solution/exercise-coupling-ff-pm/turbulence/CMakeLists.txt b/exercises/solution/exercise-coupling-ff-pm/turbulence/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..e19d81c519756532c8f6fd2b7f28926723f56952
--- /dev/null
+++ b/exercises/solution/exercise-coupling-ff-pm/turbulence/CMakeLists.txt
@@ -0,0 +1,19 @@
+add_input_file_links()
+
+dune_add_test(NAME orig_turbulence_coupling_ff-pm
+ SOURCES ex_turbulence_coupling_ff-pm.cc
+ COMPILE_DEFINITIONS EXNUMBER=0
+ CMD_ARGS ex_turbulence_coupling_ff-pm.input)
+
+dune_add_test(NAME sol_a_ex_turbulence_coupling_ff-pm
+ SOURCES ex_turbulence_coupling_ff-pm.cc
+ COMPILE_DEFINITIONS EXNUMBER=1
+ CMD_ARGS ex_turbulence_coupling_ff-pm.input)
+
+dune_add_test(NAME sol_b_ex_turbulence_coupling_ff-pm
+ SOURCES ex_turbulence_coupling_ff-pm.cc
+ COMPILE_DEFINITIONS EXNUMBER=2
+ CMD_ARGS ex_turbulence_coupling_ff-pm.input)
+
+# 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)
diff --git a/exercises/solution/exercise-coupling-ff-pm/turbulence/ex_turbulence_coupling_ff-pm.cc b/exercises/solution/exercise-coupling-ff-pm/turbulence/ex_turbulence_coupling_ff-pm.cc
new file mode 100644
index 0000000000000000000000000000000000000000..9b02336e65a9ac918d99a2cd34c6353d66e258d7
--- /dev/null
+++ b/exercises/solution/exercise-coupling-ff-pm/turbulence/ex_turbulence_coupling_ff-pm.cc
@@ -0,0 +1,306 @@
+// -*- 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/geometry/diameter.hh>
+#include <dumux/linear/seqsolverbackend.hh>
+#include <dumux/assembly/fvassembler.hh>
+#include <dumux/assembly/diffmethod.hh>
+#include <dumux/discretization/methods.hh>
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/staggeredvtkoutputmodule.hh>
+#include <dumux/io/grid/gridmanager.hh>
+
+#include <dumux/multidomain/staggeredtraits.hh>
+#include <dumux/multidomain/fvassembler.hh>
+#include <dumux/multidomain/privarswitchnewtonsolver.hh>
+
+#include <dumux/multidomain/boundary/stokesdarcy/couplingmanager.hh>
+
+#include "ex_turbulence_pmproblem.hh"
+#include "ex_turbulence_ffproblem.hh"
+
+namespace Dumux {
+namespace Properties {
+
+SET_PROP(ZeroEqTypeTag, CouplingManager)
+{
+ using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, TTAG(DarcyTwoPTwoCTypeTag)>;
+ using type = Dumux::StokesDarcyCouplingManager<Traits>;
+};
+
+SET_PROP(DarcyTwoPTwoCTypeTag, CouplingManager)
+{
+ using Traits = StaggeredMultiDomainTraits<TTAG(ZeroEqTypeTag), TTAG(ZeroEqTypeTag), TypeTag>;
+ using type = Dumux::StokesDarcyCouplingManager<Traits>;
+};
+
+} // end namespace Properties
+} // end namespace Dumux
+
+int main(int argc, char** argv) try
+{
+ using namespace Dumux;
+
+ // initialize MPI, finalize is done automatically on exit
+ const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+ // print dumux start message
+ if (mpiHelper.rank() == 0)
+ DumuxMessage::print(/*firstCall=*/true);
+
+ // parse command line arguments and input file
+ Parameters::init(argc, argv);
+
+ // Define the sub problem type tags
+ using StokesTypeTag = TTAG(ZeroEqTypeTag);
+ using DarcyTypeTag = TTAG(DarcyTwoPTwoCTypeTag);
+
+ // try to create a grid (from the given grid file or the input file)
+ // for both sub-domains
+ using DarcyGridManager = Dumux::GridManager<typename GET_PROP_TYPE(DarcyTypeTag, Grid)>;
+ DarcyGridManager darcyGridManager;
+ darcyGridManager.init("Darcy"); // pass parameter group
+
+ using StokesGridManager = Dumux::GridManager<typename GET_PROP_TYPE(StokesTypeTag, Grid)>;
+ StokesGridManager stokesGridManager;
+ stokesGridManager.init("Stokes"); // pass parameter group
+
+ // we compute on the leaf grid view
+ const auto& darcyGridView = darcyGridManager.grid().leafGridView();
+ const auto& stokesGridView = stokesGridManager.grid().leafGridView();
+
+ // create the finite volume grid geometry
+ using StokesFVGridGeometry = typename GET_PROP_TYPE(StokesTypeTag, FVGridGeometry);
+ auto stokesFvGridGeometry = std::make_shared<StokesFVGridGeometry>(stokesGridView);
+ stokesFvGridGeometry->update();
+ using DarcyFVGridGeometry = typename GET_PROP_TYPE(DarcyTypeTag, FVGridGeometry);
+ auto darcyFvGridGeometry = std::make_shared<DarcyFVGridGeometry>(darcyGridView);
+ darcyFvGridGeometry->update();
+
+ using Traits = StaggeredMultiDomainTraits<StokesTypeTag, StokesTypeTag, DarcyTypeTag>;
+
+ // the coupling manager
+ using CouplingManager = StokesDarcyCouplingManager<Traits>;
+ auto couplingManager = std::make_shared<CouplingManager>(stokesFvGridGeometry, darcyFvGridGeometry);
+
+ // the indices
+ constexpr auto stokesCellCenterIdx = CouplingManager::stokesCellCenterIdx;
+ constexpr auto stokesFaceIdx = CouplingManager::stokesFaceIdx;
+ constexpr auto darcyIdx = CouplingManager::darcyIdx;
+
+ // the problem (initial and boundary conditions)
+ using StokesProblem = typename GET_PROP_TYPE(StokesTypeTag, Problem);
+ auto stokesProblem = std::make_shared<StokesProblem>(stokesFvGridGeometry, couplingManager);
+ using DarcyProblem = typename GET_PROP_TYPE(DarcyTypeTag, Problem);
+ auto darcyProblem = std::make_shared<DarcyProblem>(darcyFvGridGeometry, couplingManager);
+
+ // initialize the fluidsystem (tabulation)
+ GET_PROP_TYPE(StokesTypeTag, FluidSystem)::init();
+
+ // get some time loop parameters
+ using Scalar = typename GET_PROP_TYPE(StokesTypeTag, Scalar);
+ const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+ const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+ auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+
+ // check if we are about to restart a previously interrupted simulation
+ Scalar restartTime = 0;
+ if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart"))
+ restartTime = getParam<Scalar>("TimeLoop.Restart");
+
+ // instantiate time loop
+ auto timeLoop = std::make_shared<TimeLoop<Scalar>>(restartTime, 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());
+
+ // apply initial solution for instationary problems
+ // auxiliary free flow solution vector
+ typename GET_PROP_TYPE(StokesTypeTag, SolutionVector) stokesSol;
+ stokesSol[stokesCellCenterIdx].resize(sol[stokesCellCenterIdx].size());
+ stokesSol[stokesFaceIdx].resize(sol[stokesFaceIdx].size());
+ stokesProblem->applyInitialSolution(stokesSol);
+ auto solStokesOld = stokesSol;
+ sol[stokesCellCenterIdx] = stokesSol[stokesCellCenterIdx];
+ sol[stokesFaceIdx] = stokesSol[stokesFaceIdx];
+ // TODO: update static wall properties
+ // TODO: update dynamic wall properties
+#if EXNUMBER >= 1
+ stokesProblem->updateStaticWallProperties();
+ stokesProblem->updateDynamicWallProperties(stokesSol);
+#endif
+
+ darcyProblem->applyInitialSolution(sol[darcyIdx]);
+ auto solDarcyOld = sol[darcyIdx];
+
+ auto solOld = sol;
+
+ couplingManager->init(stokesProblem, darcyProblem, sol);
+
+ // the grid variables
+ using StokesGridVariables = typename GET_PROP_TYPE(StokesTypeTag, GridVariables);
+ auto stokesGridVariables = std::make_shared<StokesGridVariables>(stokesProblem, stokesFvGridGeometry);
+ stokesGridVariables->init(stokesSol, solStokesOld);
+ using DarcyGridVariables = typename GET_PROP_TYPE(DarcyTypeTag, GridVariables);
+ auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyFvGridGeometry);
+ darcyGridVariables->init(sol[darcyIdx], solDarcyOld);
+
+ // intialize the vtk output module
+ const auto stokesName = getParam<std::string>("Problem.Name") + "_" + stokesProblem->name();
+ const auto darcyName = getParam<std::string>("Problem.Name") + "_" + darcyProblem->name();
+
+ StaggeredVtkOutputModule<StokesTypeTag, GET_PROP_VALUE(StokesTypeTag, PhaseIdx)> stokesVtkWriter(*stokesProblem, *stokesFvGridGeometry, *stokesGridVariables, stokesSol, stokesName);
+ GET_PROP_TYPE(StokesTypeTag, VtkOutputFields)::init(stokesVtkWriter);
+ stokesVtkWriter.write(0.0);
+
+ VtkOutputModule<DarcyTypeTag> darcyVtkWriter(*darcyProblem, *darcyFvGridGeometry, *darcyGridVariables, sol[darcyIdx], darcyName);
+ GET_PROP_TYPE(DarcyTypeTag, VtkOutputFields)::init(darcyVtkWriter);
+ darcyVtkWriter.write(0.0);
+
+ // the assembler with time loop for instationary problem
+ using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
+ auto assembler = std::make_shared<Assembler>(std::make_tuple(stokesProblem, stokesProblem, darcyProblem),
+ std::make_tuple(stokesFvGridGeometry->cellCenterFVGridGeometryPtr(),
+ stokesFvGridGeometry->faceFVGridGeometryPtr(),
+ darcyFvGridGeometry),
+ std::make_tuple(stokesGridVariables->cellCenterGridVariablesPtr(),
+ stokesGridVariables->faceGridVariablesPtr(),
+ darcyGridVariables),
+ couplingManager,
+ timeLoop);
+
+ // the linear solver
+ using LinearSolver = UMFPackBackend;
+ auto linearSolver = std::make_shared<LinearSolver>();
+
+ // the primary variable switches used by the sub models
+ using PriVarSwitchTuple = std::tuple<NoPrimaryVariableSwitch, NoPrimaryVariableSwitch, typename GET_PROP_TYPE(DarcyTypeTag, PrimaryVariableSwitch)>;
+
+ // the non-linear solver
+ using NewtonSolver = MultiDomainPriVarSwitchNewtonSolver<Assembler, LinearSolver, CouplingManager, PriVarSwitchTuple>;
+ 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;
+
+ // update the auxiliary free flow solution vector
+ stokesSol[stokesCellCenterIdx] = sol[stokesCellCenterIdx];
+ stokesSol[stokesFaceIdx] = sol[stokesFaceIdx];
+
+#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/exercises/solution/exercise-coupling-ff-pm/turbulence/ex_turbulence_coupling_ff-pm.input b/exercises/solution/exercise-coupling-ff-pm/turbulence/ex_turbulence_coupling_ff-pm.input
new file mode 100644
index 0000000000000000000000000000000000000000..514cf156d55211777c2ec0980f37d1a95124867e
--- /dev/null
+++ b/exercises/solution/exercise-coupling-ff-pm/turbulence/ex_turbulence_coupling_ff-pm.input
@@ -0,0 +1,62 @@
+[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]
+
+[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
diff --git a/exercises/solution/exercise-coupling-ff-pm/turbulence/ex_turbulence_ffproblem.hh b/exercises/solution/exercise-coupling-ff-pm/turbulence/ex_turbulence_ffproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..604cd024a422fefca6c1244981527de306be968d
--- /dev/null
+++ b/exercises/solution/exercise-coupling-ff-pm/turbulence/ex_turbulence_ffproblem.hh
@@ -0,0 +1,409 @@
+// -*- 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/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
+{
+#if EXNUMBER >= 1
+NEW_TYPE_TAG(ZeroEqTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, ZeroEqNCNI));
+#else
+NEW_TYPE_TAG(ZeroEqTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, NavierStokesNCNI));
+#endif
+
+// Set the grid type
+SET_TYPE_PROP(ZeroEqTypeTag, Grid, Dune::YaspGrid<2, Dune::TensorProductCoordinates<typename GET_PROP_TYPE(TypeTag, Scalar), 2> >);
+
+// set the fluid system
+SET_TYPE_PROP(ZeroEqTypeTag, FluidSystem, FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>);
+
+// set phase index (air)
+SET_INT_PROP(ZeroEqTypeTag, PhaseIdx, GET_PROP_TYPE(TypeTag, FluidSystem)::gasPhaseIdx);
+
+SET_INT_PROP(ZeroEqTypeTag, ReplaceCompEqIdx, 3);
+
+// Use formulation based on mass fractions
+SET_BOOL_PROP(ZeroEqTypeTag, UseMoles, true);
+
+// Set the problem property
+SET_TYPE_PROP(ZeroEqTypeTag, Problem, Dumux::FreeFlowSubProblem<TypeTag> );
+
+SET_BOOL_PROP(ZeroEqTypeTag, EnableFVGridGeometryCache, true);
+SET_BOOL_PROP(ZeroEqTypeTag, EnableGridFluxVariablesCache, true);
+SET_BOOL_PROP(ZeroEqTypeTag, EnableGridVolumeVariablesCache, true);
+
+SET_BOOL_PROP(ZeroEqTypeTag, EnableInertiaTerms, 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 = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
+ using ElementFaceVariables = typename GET_PROP_TYPE(TypeTag, GridFaceVariables)::LocalView;
+ using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
+
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+
+ using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager);
+ using TimeLoopPtr = std::shared_ptr<TimeLoop<Scalar>>;
+
+ using DiffusionCoefficientAveragingType = typename StokesDarcyCouplingOptions::DiffusionCoefficientAveragingType;
+
+ static constexpr bool useMoles = GET_PROP_TYPE(TypeTag, ModelTraits)::useMoles();
+
+ static constexpr auto phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIdx);
+ static constexpr auto transportCompIdx = 1 - phaseIdx;
+
+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);
+ values.setDirichlet(Indices::energyBalanceIdx);
+ }
+
+ if (onLowerBoundary_(globalPos))
+ {
+ values.setDirichlet(Indices::velocityXIdx);
+ values.setDirichlet(Indices::velocityYIdx);
+ values.setNeumann(Indices::conti0EqIdx);
+ values.setNeumann(Indices::conti0EqIdx + 1);
+ values.setNeumann(Indices::energyBalanceIdx);
+ }
+
+ 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::energyBalanceIdx);
+#endif
+ }
+
+ if (onRightBoundary_(globalPos))
+ {
+ values.setDirichlet(Indices::pressureIdx);
+ values.setOutflow(Indices::conti0EqIdx);
+ values.setOutflow(Indices::energyBalanceIdx);
+ }
+
+ if (couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
+ {
+ values.setCouplingNeumann(Indices::conti0EqIdx);
+ values.setCouplingNeumann(Indices::conti0EqIdx + 1);
+ values.setCouplingNeumann(Indices::energyBalanceIdx);
+ values.setCouplingNeumann(Indices::momentumYBalanceIdx);
+ values.setBJS(Indices::momentumXBalanceIdx);
+ }
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Dirichlet control volume.
+ *
+ * \param element The element
+ * \param scvf The subcontrolvolume face
+ */
+ PrimaryVariables dirichletAtPos(const GlobalPosition& pos) const
+ {
+ PrimaryVariables values(0.0);
+ values = initialAtPos(pos);
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a Neumann control volume.
+ *
+ * \param element The element for which the Neumann boundary condition is set
+ * \param fvGeomentry The fvGeometry
+ * \param elemVolVars The element volume variables
+ * \param elemFaceVars The element face variables
+ * \param scvf The boundary sub control volume face
+ */
+ NumEqVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFaceVariables& elemFaceVars,
+ const SubControlVolumeFace& scvf) const
+ {
+ PrimaryVariables values(0.0);
+ if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
+ {
+ values[Indices::momentumYBalanceIdx] = couplingManager().couplingData().momentumCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf);
+
+ const auto massFlux = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf, diffCoeffAvgType_);
+ values[Indices::conti0EqIdx] = massFlux[0];
+ values[Indices::conti0EqIdx + 1] = massFlux[1];
+ values[Indices::energyBalanceIdx] = couplingManager().couplingData().energyCouplingCondition(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 isOnWall(const GlobalPosition& globalPos) const
+ {
+ return (onLowerBoundary_(globalPos));
+ }
+#elif EXNUMBER >= 1
+ bool isOnWall(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, phaseIdx);
+
+ PrimaryVariables values(0.0);
+ values[Indices::pressureIdx] = refPressure() + density*this->gravity()[1]*(globalPos[1] - this->fvGridGeometry().bBoxMin()[1]);
+ values[Indices::conti0EqIdx] = refMoleFrac();
+ values[Indices::velocityXIdx] = refVelocity();
+ values[Indices::temperatureIdx] = refTemperature();
+
+ if(onUpperBoundary_(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 SubControlVolumeFace& scvf) const
+ {
+ return couplingManager().couplingData().darcyPermeability(scvf);
+ }
+
+ /*!
+ * \brief Returns the alpha value required as input parameter for the Beavers-Joseph-Saffman boundary condition
+ */
+ Scalar alphaBJ(const SubControlVolumeFace& scvf) const
+ {
+ return couplingManager().problem(CouplingManager::darcyIdx).spatialParams().beaversJosephCoeffAtPos(scvf.center());
+ }
+
+ // \}
+
+private:
+ bool onLeftBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_; }
+
+ bool onRightBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_; }
+
+ bool onLowerBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] < this->fvGridGeometry().bBoxMin()[1] + eps_; }
+
+ bool onUpperBoundary_(const GlobalPosition &globalPos) const
+ { return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_; }
+
+ //! \brief updates the fluid state to obtain required quantities for IC/BC
+ void updateFluidStateForBC_(FluidState& fluidState, const Scalar temperature,
+ const Scalar pressure, const Scalar moleFraction) const
+ {
+ fluidState.setTemperature(temperature);
+ fluidState.setPressure(phaseIdx, pressure);
+ fluidState.setSaturation(phaseIdx, 1.0);
+ fluidState.setMoleFraction(phaseIdx, transportCompIdx, moleFraction);
+ fluidState.setMoleFraction(phaseIdx, phaseIdx, 1.0 - moleFraction);
+
+ typename FluidSystem::ParameterCache paramCache;
+ paramCache.updatePhase(fluidState, phaseIdx);
+
+ const Scalar density = FluidSystem::density(fluidState, paramCache, phaseIdx);
+ fluidState.setDensity(phaseIdx, density);
+
+ const Scalar molarDensity = FluidSystem::molarDensity(fluidState, paramCache, phaseIdx);
+ fluidState.setMolarDensity(phaseIdx, molarDensity);
+
+ const Scalar enthalpy = FluidSystem::enthalpy(fluidState, paramCache, phaseIdx);
+ fluidState.setEnthalpy(phaseIdx, enthalpy);
+ }
+
+ // the height of the free-flow domain
+ const Scalar height_() const
+ { return this->fvGridGeometry().bBoxMax()[1] - this->fvGridGeometry().bBoxMin()[1]; }
+
+ Scalar eps_;
+
+ Scalar 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/exercises/solution/exercise-coupling-ff-pm/turbulence/ex_turbulence_pmproblem.hh b/exercises/solution/exercise-coupling-ff-pm/turbulence/ex_turbulence_pmproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..bd589d6303f95b7c9ba96b94102174af24c8a129
--- /dev/null
+++ b/exercises/solution/exercise-coupling-ff-pm/turbulence/ex_turbulence_pmproblem.hh
@@ -0,0 +1,324 @@
+// -*- 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/cellcentered/tpfa/properties.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
+{
+NEW_TYPE_TAG(DarcyTwoPTwoCTypeTag, INHERITS_FROM(CCTpfaModel, TwoPTwoCNI));
+
+// Set the problem property
+SET_TYPE_PROP(DarcyTwoPTwoCTypeTag, Problem, Dumux::DarcySubProblem<TypeTag>);
+
+// the fluid system
+SET_TYPE_PROP(DarcyTwoPTwoCTypeTag, FluidSystem, FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>);
+
+//! Set the default formulation to pw-Sn: This can be over written in the problem.
+SET_PROP(DarcyTwoPTwoCTypeTag, Formulation)
+{ static constexpr auto value = TwoPFormulation::p1s0; };
+
+// The gas component balance (air) is replaced by the total mass balance
+SET_INT_PROP(DarcyTwoPTwoCTypeTag, ReplaceCompEqIdx, 3);
+
+// Set the grid type
+SET_TYPE_PROP(DarcyTwoPTwoCTypeTag, Grid, Dune::YaspGrid<2, Dune::TensorProductCoordinates<typename GET_PROP_TYPE(TypeTag, Scalar), 2> >);
+
+SET_BOOL_PROP(DarcyTwoPTwoCTypeTag, UseMoles, true);
+
+SET_TYPE_PROP(DarcyTwoPTwoCTypeTag, SpatialParams, TwoPSpatialParams<TypeTag>);
+}
+
+/*!
+ * \brief The porous medium sub problem
+ */
+template <class TypeTag>
+class DarcySubProblem : public PorousMediumFlowProblem<TypeTag>
+{
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
+
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+
+ // copy some indices for convenience
+ using Indices = typename GET_PROP_TYPE(TypeTag, 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 = typename GET_PROP_TYPE(TypeTag, 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(fvGeometry, elemVolVars, scvf, diffCoeffAvgType_);
+
+ for(int i = 0; i< massFlux.size(); ++i)
+ values[i] = massFlux[i];
+
+ values[Indices::energyEqIdx] = couplingManager().couplingData().energyCouplingCondition(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