From 50ac07cc0a5a35967340d96e8347cf47c1d9e912 Mon Sep 17 00:00:00 2001
From: Anna Mareike Kostelecky <anmako96@web.de>
Date: Mon, 15 Jul 2024 09:23:57 +0000
Subject: [PATCH] [exercises] update coupling ffpm exercise
---
.../exercise-coupling-ff-pm.patch | 1454 ++++++++++++-----
exercises/exercise-coupling-ff-pm/README.md | 499 ++----
.../interface/freeflowsubproblem.hh | 18 +-
.../exercise-coupling-ff-pm/interface/main.cc | 7 +-
.../interface/porousmediumsubproblem.hh | 11 -
.../interface/properties.hh | 3 +
.../interface/readme.md | 169 ++
.../models/CMakeLists.txt | 2 +-
.../models/params.input | 10 +-
.../models/porousmediumsubproblem.hh | 2 +
.../exercise-coupling-ff-pm/models/readme.md | 118 ++
.../turbulence/params.input | 6 +-
.../turbulence/readme.md | 126 ++
.../interface/freeflowsubproblem.hh | 35 +-
.../exercise-coupling-ff-pm/interface/main.cc | 8 +-
.../interface/porousmediumsubproblem.hh | 11 -
.../interface/properties.hh | 3 +
.../models/CMakeLists.txt | 10 +-
.../models/{params.input => params_b_c.input} | 7 +-
.../models/params_orig_a.input | 56 +
.../turbulence/CMakeLists.txt | 19 +-
.../turbulence/params_b.input | 72 +
.../turbulence/params_c_d.input | 72 +
.../{params.input => params_orig_a.input} | 2 -
24 files changed, 1792 insertions(+), 928 deletions(-)
create mode 100644 exercises/exercise-coupling-ff-pm/interface/readme.md
create mode 100644 exercises/exercise-coupling-ff-pm/models/readme.md
create mode 100644 exercises/exercise-coupling-ff-pm/turbulence/readme.md
rename exercises/solution/exercise-coupling-ff-pm/models/{params.input => params_b_c.input} (92%)
create mode 100644 exercises/solution/exercise-coupling-ff-pm/models/params_orig_a.input
create mode 100644 exercises/solution/exercise-coupling-ff-pm/turbulence/params_b.input
create mode 100644 exercises/solution/exercise-coupling-ff-pm/turbulence/params_c_d.input
rename exercises/solution/exercise-coupling-ff-pm/turbulence/{params.input => params_orig_a.input} (96%)
diff --git a/.patches/exercise-coupling-ff-pm/exercise-coupling-ff-pm.patch b/.patches/exercise-coupling-ff-pm/exercise-coupling-ff-pm.patch
index 7dba7786..4dc2abd2 100644
--- a/.patches/exercise-coupling-ff-pm/exercise-coupling-ff-pm.patch
+++ b/.patches/exercise-coupling-ff-pm/exercise-coupling-ff-pm.patch
@@ -1,6 +1,6 @@
diff -ruN exercises/exercise-coupling-ff-pm/interface/CMakeLists.txt exercises/solution/exercise-coupling-ff-pm/interface/CMakeLists.txt
---- exercises/exercise-coupling-ff-pm/interface/CMakeLists.txt 2024-05-21 14:15:07.141554755 +0200
-+++ exercises/solution/exercise-coupling-ff-pm/interface/CMakeLists.txt 2024-05-21 14:15:07.209555087 +0200
+--- exercises/exercise-coupling-ff-pm/interface/CMakeLists.txt 2024-07-10 08:19:20.139465833 +0200
++++ exercises/solution/exercise-coupling-ff-pm/interface/CMakeLists.txt 2023-04-04 12:15:36.845203367 +0200
@@ -1,7 +1,27 @@
-# executables for ex_interface_coupling_ff-pm
-dumux_add_test(NAME exercise_interface_coupling_ff-pm
@@ -33,8 +33,8 @@ diff -ruN exercises/exercise-coupling-ff-pm/interface/CMakeLists.txt exercises/s
# add a symlink for each input file
add_input_file_links()
diff -ruN exercises/exercise-coupling-ff-pm/interface/freeflowsubproblem.hh exercises/solution/exercise-coupling-ff-pm/interface/freeflowsubproblem.hh
---- exercises/exercise-coupling-ff-pm/interface/freeflowsubproblem.hh 2024-05-29 14:31:49.241625242 +0200
-+++ exercises/solution/exercise-coupling-ff-pm/interface/freeflowsubproblem.hh 2024-05-21 14:15:07.209555087 +0200
+--- exercises/exercise-coupling-ff-pm/interface/freeflowsubproblem.hh 2024-07-15 10:49:12.269896347 +0200
++++ exercises/solution/exercise-coupling-ff-pm/interface/freeflowsubproblem.hh 2024-07-15 10:49:12.273896458 +0200
@@ -83,8 +83,7 @@
const auto& globalPos = scvf.dofPosition();
@@ -96,46 +96,44 @@ diff -ruN exercises/exercise-coupling-ff-pm/interface/freeflowsubproblem.hh exer
}
return values;
-@@ -150,9 +174,13 @@
- 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;
- }
-
-@@ -172,10 +200,18 @@
- PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+@@ -124,9 +148,18 @@
+ PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
{
PrimaryVariables values(0.0);
+#if EXNUMBER == 0
values[Indices::velocityYIdx] = -1e-6 * globalPos[0] * (this->gridGeometry().bBoxMax()[0] - globalPos[0]);
- // TODO: dumux-course-task 1.A
-- // Set fixed pressures on the left and right boundary
--
+#elif EXNUMBER == 4
+ values[Indices::velocityXIdx] = 1e-6 * (globalPos[1] - this->gridGeometry().bBoxMin()[1])
+ * (this->gridGeometry().bBoxMax()[1] - globalPos[1]);
+#else
-+ // set fixed pressures on the left and right boundary
+ // 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;
+ }
+@@ -152,9 +185,13 @@
+ 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;
}
diff -ruN exercises/exercise-coupling-ff-pm/interface/main.cc exercises/solution/exercise-coupling-ff-pm/interface/main.cc
---- exercises/exercise-coupling-ff-pm/interface/main.cc 2024-05-29 14:31:49.241625242 +0200
-+++ exercises/solution/exercise-coupling-ff-pm/interface/main.cc 2024-05-21 14:15:07.209555087 +0200
+--- exercises/exercise-coupling-ff-pm/interface/main.cc 2024-07-15 10:49:12.269896347 +0200
++++ exercises/solution/exercise-coupling-ff-pm/interface/main.cc 2024-07-15 10:49:12.273896458 +0200
@@ -64,10 +64,8 @@
using FreeflowTypeTag = Properties::TTag::FreeflowOneP;
using PorousMediumTypeTag = Properties::TTag::PorousMediumFlowOneP;
@@ -250,7 +248,7 @@ diff -ruN exercises/exercise-coupling-ff-pm/interface/main.cc exercises/solution
// create the finite volume grid geometry
using FreeflowFVGridGeometry = GetPropType<FreeflowTypeTag, Properties::GridGeometry>;
-@@ -183,9 +175,10 @@
+@@ -180,9 +172,9 @@
StaggeredVtkOutputModule<FreeflowGridVariables, decltype(freeflowSol)> freeflowVtkWriter(*freeflowGridVariables, freeflowSol, freeflowName);
GetPropType<FreeflowTypeTag, Properties::IOFields>::initOutputModule(freeflowVtkWriter);
@@ -260,13 +258,12 @@ diff -ruN exercises/exercise-coupling-ff-pm/interface/main.cc exercises/solution
+#if EXNUMBER >= 2
+ freeflowVtkWriter.addField(freeflowProblem->getAnalyticalVelocityX(), "analyticalV_x");
+#endif
-+
- freeflowVtkWriter.write(0.0);
using PorousMediumSolutionVector = GetPropType<PorousMediumTypeTag, Properties::SolutionVector>;
+ VtkOutputModule<PorousMediumGridVariables, PorousMediumSolutionVector> porousMediumVtkWriter(*porousMediumGridVariables,
diff -ruN exercises/exercise-coupling-ff-pm/interface/porousmediumsubproblem.hh exercises/solution/exercise-coupling-ff-pm/interface/porousmediumsubproblem.hh
---- exercises/exercise-coupling-ff-pm/interface/porousmediumsubproblem.hh 2024-05-29 14:31:49.241625242 +0200
-+++ exercises/solution/exercise-coupling-ff-pm/interface/porousmediumsubproblem.hh 2024-05-21 14:15:07.209555087 +0200
+--- exercises/exercise-coupling-ff-pm/interface/porousmediumsubproblem.hh 2024-07-15 10:49:12.269896347 +0200
++++ exercises/solution/exercise-coupling-ff-pm/interface/porousmediumsubproblem.hh 2024-07-15 10:49:12.273896458 +0200
@@ -80,13 +80,13 @@
// set Neumann BCs to all boundaries first
values.setAllNeumann();
@@ -287,8 +284,8 @@ diff -ruN exercises/exercise-coupling-ff-pm/interface/porousmediumsubproblem.hh
return values;
}
diff -ruN exercises/exercise-coupling-ff-pm/interface/properties.hh exercises/solution/exercise-coupling-ff-pm/interface/properties.hh
---- exercises/exercise-coupling-ff-pm/interface/properties.hh 2024-05-29 14:31:49.242625227 +0200
-+++ exercises/solution/exercise-coupling-ff-pm/interface/properties.hh 2024-05-21 14:15:07.209555087 +0200
+--- exercises/exercise-coupling-ff-pm/interface/properties.hh 2024-07-15 10:49:12.269896347 +0200
++++ exercises/solution/exercise-coupling-ff-pm/interface/properties.hh 2024-07-15 10:49:12.273896458 +0200
@@ -28,9 +28,9 @@
#include <dumux/multidomain/staggeredtraits.hh>
#include <dumux/multidomain/boundary/stokesdarcy/couplingmanager.hh>
@@ -340,11 +337,185 @@ diff -ruN exercises/exercise-coupling-ff-pm/interface/properties.hh exercises/so
};
template<class TypeTag>
+diff -ruN exercises/exercise-coupling-ff-pm/interface/readme.md exercises/solution/exercise-coupling-ff-pm/interface/readme.md
+--- exercises/exercise-coupling-ff-pm/interface/readme.md 2024-07-15 10:52:11.058859896 +0200
++++ exercises/solution/exercise-coupling-ff-pm/interface/readme.md 1970-01-01 01:00:00.000000000 +0100
+@@ -1,169 +0,0 @@
+-## 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:
+-
+-
+-
+-Note that we neglect the influence of gravity and only solve for the stationary problem in this exercise.
+-
+-You can also compile and run the test without changing anything yet and investigate the above mentioned velocity distribution with paraview.
+-
+-* 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 `interface/freeflowsubproblem.hh` and navigate to the part, 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.
+-You can also delete the vertical velocity from `dirichletAtPos`, as Dirichlet-type boundary conditions should be set for the pressure.
+-
+-Recompile and rerun your code. Then check your results with the applied changes by visualizing them e.g. with paraview.
+-
+-### 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.setBeaversJoseph(Indices::momentumXBalanceIdx);
+-```
+-
+-at the position where the coupling boundary conditions are set in `interface/freeflowsubproblem.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 the following line in `main.cc` in the subfolder `interface`.
+-```cpp
+-freeflowVtkWriter.addField(freeflowProblem->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 in the input file to see how that affects the resulting velocity profile.
+-
+-### Task C: Change the shape of interface
+-
+-Now we want to include a non-flat interface between the two domains.
+-We use [`dune-subgrid`](https://doi.org/10.1007/s00607-009-0067-2) to construct two grids for the two domains from one common host grid. Thus, for the following tasks subgrid needs to be correctly installed.
+-Our hostgrid will be a Dune-Yasp grid (`Dune::YaspGrid<2, Dune::TensorProductCoordinates<double, dim> >`)
+-and the bounds and resolution of the grid will be set in the `params.input`file under the group `[Grid]`.
+-This hostgrid, along with `elementSelector` functions defining some spatial cut of this domain, are passed to the grid manager to create each subdomain grid.
+-
+-To introduce this, open `main.cc` in the subfolder `interface` again and search for `TODO: dumux-course-task 1.C`.
+-Comment out the first code block and uncomment the second.
+-This will instantiate a host grid and define two helper lambda functions that are used to choose elements from the host grid for the respective sub grid.
+-In the given case, the domain is split in two halves, separated by a sinusoidal interface.
+-
+-```cpp
+-auto elementSelectorFreeflow = [&](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 elementSelectorPorousMedium = [&](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 line for including the grid manager in the properties file, i.e.
+-```cpp
+-#include <dumux/io/grid/gridmanager_sub.hh>
+-```
+-
+-and do the changes in the respective lines for the `Grid` property.
+-
+-The problem should now compile. However, a runtime 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());
+- ```
+-in `freeflowsubproblem.hh` in the subfolder `interface`.
+-
+-The same is true for the BJS condition, however, here we need to consider the tangential direction:
+-```cpp
+-values.setBeaversJoseph(1 - scvf.directionIndex());
+-```
+-
+-Recompile and rerun to check if the final result looks something like this:
+-
+-
+-
+-*Extra Points:*
+-Rather than enforcing a pressure difference across the domain, an inflow velocity profile could be set.
+-What changes to the left boundary conditions in the free-flow domain would you make to introduce this? What conditions can be enforced on the right boundary?
+-Hint: A relation between velocity and position is used for the vertical velocity component in the original form of the `dirichletAtPos` method.
+\ No newline at end of file
diff -ruN exercises/exercise-coupling-ff-pm/models/CMakeLists.txt exercises/solution/exercise-coupling-ff-pm/models/CMakeLists.txt
---- exercises/exercise-coupling-ff-pm/models/CMakeLists.txt 2024-05-21 14:15:07.141554755 +0200
-+++ exercises/solution/exercise-coupling-ff-pm/models/CMakeLists.txt 2024-05-21 14:15:07.209555087 +0200
+--- exercises/exercise-coupling-ff-pm/models/CMakeLists.txt 2024-07-15 10:49:12.269896347 +0200
++++ exercises/solution/exercise-coupling-ff-pm/models/CMakeLists.txt 2024-07-15 10:49:12.273896458 +0200
@@ -1,6 +1,29 @@
--# executables for ex_interface_coupling_ff-pm
+-# executables for exercise_models_coupling_ff-pm
-dumux_add_test(NAME exercise_models_coupling_ff-pm
+dumux_add_test(NAME exercise_models_coupling_ff-pm_original
SOURCES main.cc
@@ -352,33 +523,34 @@ diff -ruN exercises/exercise-coupling-ff-pm/models/CMakeLists.txt exercises/solu
+ LABELS ffpm
+ COMPILE_DEFINITIONS EXNUMBER=0
+ COMMAND ${CMAKE_CURRENT_BINARY_DIR}/exercise_models_coupling_ff-pm_original
-+ CMD_ARGS -Problem.ExportStorage false -Problem.PlotStorage 0 -Problem.ExportFluxes 0)
-+
++ CMD_ARGS params_orig_a.input)
+
+-dune_symlink_to_source_files(FILES "params.input" plotFluxes.py)
+dumux_add_test(NAME exercise_models_coupling_ff-pm_a_solution
+ SOURCES main.cc
+ LABELS ffpm
+ COMPILE_DEFINITIONS EXNUMBER=1
+ COMMAND ${CMAKE_CURRENT_BINARY_DIR}/exercise_models_coupling_ff-pm_a_solution
-+ CMD_ARGS -Problem.ExportStorage false -Problem.PlotStorage 0 -Problem.ExportFluxes 0)
++ CMD_ARGS params_orig_a.input)
+
+dumux_add_test(NAME exercise_models_coupling_ff-pm_b_solution
+ SOURCES main.cc
+ LABELS ffpm
+ COMPILE_DEFINITIONS EXNUMBER=2
+ COMMAND ${CMAKE_CURRENT_BINARY_DIR}/exercise_models_coupling_ff-pm_b_solution
-+ CMD_ARGS -Problem.ExportStorage false -Problem.PlotStorage 0 -Problem.ExportFluxes 0)
++ CMD_ARGS params_b_c.input)
+
+dumux_add_test(NAME exercise_models_coupling_ff-pm_c_solution
+ SOURCES main.cc
+ LABELS ffpm
+ COMPILE_DEFINITIONS EXNUMBER=3
+ COMMAND ${CMAKE_CURRENT_BINARY_DIR}/exercise_models_coupling_ff-pm_c_solution
-+ CMD_ARGS -Problem.ExportStorage false -Problem.PlotStorage 0 -Problem.ExportFluxes 0)
-
- dune_symlink_to_source_files(FILES "params.input" plotFluxes.py)
++ CMD_ARGS params_b_c.input)
++
++dune_symlink_to_source_files(FILES "params_orig_a.input" "params_b_c.input" plotFluxes.py)
diff -ruN exercises/exercise-coupling-ff-pm/models/main.cc exercises/solution/exercise-coupling-ff-pm/models/main.cc
---- exercises/exercise-coupling-ff-pm/models/main.cc 2024-05-21 14:15:07.141554755 +0200
-+++ exercises/solution/exercise-coupling-ff-pm/models/main.cc 2024-05-21 14:15:07.209555087 +0200
+--- exercises/exercise-coupling-ff-pm/models/main.cc 2023-04-04 12:15:36.845203367 +0200
++++ exercises/solution/exercise-coupling-ff-pm/models/main.cc 2023-07-26 11:31:46.014485560 +0200
@@ -139,9 +139,15 @@
auto porousMediumGridVariables = std::make_shared<PorousMediumGridVariables>(porousMediumProblem, porousMediumFvGridGeometry);
porousMediumGridVariables->init(sol[porousMediumIdx]);
@@ -398,19 +570,192 @@ diff -ruN exercises/exercise-coupling-ff-pm/models/main.cc exercises/solution/ex
StaggeredVtkOutputModule<FreeflowGridVariables, decltype(freeflowSol)> freeflowVtkWriter(*freeflowGridVariables, freeflowSol, freeflowName);
GetPropType<FreeflowTypeTag, Properties::IOFields>::initOutputModule(freeflowVtkWriter);
+diff -ruN exercises/exercise-coupling-ff-pm/models/params_b_c.input exercises/solution/exercise-coupling-ff-pm/models/params_b_c.input
+--- exercises/exercise-coupling-ff-pm/models/params_b_c.input 1970-01-01 01:00:00.000000000 +0100
++++ exercises/solution/exercise-coupling-ff-pm/models/params_b_c.input 2024-07-15 10:49:12.273896458 +0200
+@@ -0,0 +1,57 @@
++[TimeLoop]
++DtInitial = 100 # s
++EpisodeLength = -360 # s # 0.25 days
++TEnd = 256000 # s # 2 days
++
++[Freeflow.Grid]
++LowerLeft = 0 1
++UpperRight = 1 2
++Cells = 16 16
++
++[PorousMedium.Grid]
++UpperRight = 1 1
++Cells = 16 16
++
++[Freeflow.Problem]
++Name = freeflow
++EnableGravity = true
++EnableInertiaTerms = true
++MoleFraction = 0.0 # -
++Pressure = 1e5 # Pa
++Velocity = 1 # m/s
++
++[PorousMedium.Problem]
++Name = porousmedium
++EnableGravity = true
++MoleFraction = 0.1 # -
++Saturation = 0.1 # -
++
++[SpatialParams]
++Permeability = 2.65e-10 # m^2
++Porosity = 0.4 # -
++AlphaBeaversJoseph = 1.0 # -
++# EXNUMBER >= 1
++VanGenuchtenN = 8.0
++VanGenuchtenAlpha = 6.5e-4
++Swr = 0.005
++Snr = 0.01
++VanGenuchtenPcLowSweThreshold = 0.01
++VanGenuchtenPcHighSweThreshold = 0.99
++Temperature = 293.15
++
++[Problem]
++Name = models_coupling
++ExportStorage = true
++PlotStorage = true
++ExportFluxes = true
++
++[Newton]
++MaxSteps = 12
++MaxRelativeShift = 1e-7
++TargetSteps = 5
++
++[Vtk]
++AddVelocity = 1
++
++[Assembly]
++NumericDifference.BaseEpsilon = 1e-8
diff -ruN exercises/exercise-coupling-ff-pm/models/params.input exercises/solution/exercise-coupling-ff-pm/models/params.input
---- exercises/exercise-coupling-ff-pm/models/params.input 2024-05-21 14:15:07.141554755 +0200
-+++ exercises/solution/exercise-coupling-ff-pm/models/params.input 2024-05-21 14:15:07.209555087 +0200
-@@ -1,5 +1,5 @@
- [TimeLoop]
+--- exercises/exercise-coupling-ff-pm/models/params.input 2024-07-15 10:49:12.269896347 +0200
++++ exercises/solution/exercise-coupling-ff-pm/models/params.input 1970-01-01 01:00:00.000000000 +0100
+@@ -1,58 +0,0 @@
+-[TimeLoop]
-DtInitial = 1 # s
+-EpisodeLength = -360 # s # 0.25 days
+-TEnd = 256000 # s # 2 days
+-
+-[Freeflow.Grid]
+-LowerLeft = 0 1
+-UpperRight = 1 2
+-Cells = 16 16
+-
+-[PorousMedium.Grid]
+-UpperRight = 1 1
+-Cells = 16 16
+-
+-[Freeflow.Problem]
+-Name = freeflow
+-EnableGravity = true
+-EnableInertiaTerms = true
+-MoleFraction = 0.0 # -
+-Pressure = 1e5 # Pa
+-Velocity = 1 # m/s
+-
+-[PorousMedium.Problem]
+-Name = porousmedium
+-EnableGravity = true
+-Saturation = 0.1 # -
+-MoleFraction = 0.1 # -
+-Pressure = 1e5 # Pa
+-
+-[SpatialParams]
+-Permeability = 2.65e-10 # m^2
+-Porosity = 0.4 # -
+-AlphaBeaversJoseph = 1.0 # -
+-# EXNUMBER >= 1
+-VanGenuchtenN = 8.0
+-VanGenuchtenAlpha = 6.5e-4
+-Swr = 0.005
+-Snr = 0.01
+-VanGenuchtenPcLowSweThreshold = 0.01
+-VanGenuchtenPcHighSweThreshold = 0.99
+-Temperature = 293.15
+-
+-[Problem]
+-Name = models_coupling
+-ExportStorage = false
+-PlotStorage = false
+-ExportFluxes = false
+-
+-[Newton]
+-MaxSteps = 12
+-MaxRelativeShift = 1e-7
+-TargetSteps = 5
+-
+-[Vtk]
+-AddVelocity = 1
+-
+-[Assembly]
+-NumericDifference.BaseEpsilon = 1e-8
+diff -ruN exercises/exercise-coupling-ff-pm/models/params_orig_a.input exercises/solution/exercise-coupling-ff-pm/models/params_orig_a.input
+--- exercises/exercise-coupling-ff-pm/models/params_orig_a.input 1970-01-01 01:00:00.000000000 +0100
++++ exercises/solution/exercise-coupling-ff-pm/models/params_orig_a.input 2024-07-15 10:49:12.273896458 +0200
+@@ -0,0 +1,56 @@
++[TimeLoop]
+DtInitial = 100 # s
- EpisodeLength = -360 # s # 0.25 days
- TEnd = 256000 # s # 2 days
-
++EpisodeLength = -360 # s # 0.25 days
++TEnd = 256000 # s # 2 days
++
++[Freeflow.Grid]
++LowerLeft = 0 1
++UpperRight = 1 2
++Cells = 16 16
++
++[PorousMedium.Grid]
++UpperRight = 1 1
++Cells = 16 16
++
++[Freeflow.Problem]
++Name = freeflow
++EnableGravity = true
++EnableInertiaTerms = true
++MoleFraction = 0.0 # -
++Pressure = 1e5 # Pa
++Velocity = 1 # m/s
++
++[PorousMedium.Problem]
++Name = porousmedium
++EnableGravity = true
++MoleFraction = 0.1 # -
++
++[SpatialParams]
++Permeability = 2.65e-10 # m^2
++Porosity = 0.4 # -
++AlphaBeaversJoseph = 1.0 # -
++# EXNUMBER >= 1
++VanGenuchtenN = 8.0
++VanGenuchtenAlpha = 6.5e-4
++Swr = 0.005
++Snr = 0.01
++VanGenuchtenPcLowSweThreshold = 0.01
++VanGenuchtenPcHighSweThreshold = 0.99
++Temperature = 293.15
++
++[Problem]
++Name = models_coupling
++ExportStorage = false
++PlotStorage = false
++ExportFluxes = false
++
++[Newton]
++MaxSteps = 12
++MaxRelativeShift = 1e-7
++TargetSteps = 5
++
++[Vtk]
++AddVelocity = 1
++
++[Assembly]
++NumericDifference.BaseEpsilon = 1e-8
diff -ruN exercises/exercise-coupling-ff-pm/models/porousmediumsubproblem.hh exercises/solution/exercise-coupling-ff-pm/models/porousmediumsubproblem.hh
---- exercises/exercise-coupling-ff-pm/models/porousmediumsubproblem.hh 2024-05-29 14:31:49.243625213 +0200
-+++ exercises/solution/exercise-coupling-ff-pm/models/porousmediumsubproblem.hh 2024-05-21 14:15:07.209555087 +0200
+--- exercises/exercise-coupling-ff-pm/models/porousmediumsubproblem.hh 2024-07-15 10:49:12.269896347 +0200
++++ exercises/solution/exercise-coupling-ff-pm/models/porousmediumsubproblem.hh 2024-01-03 15:09:59.465499101 +0100
@@ -64,10 +64,15 @@
// primary variable indices
static constexpr int conti0EqIdx = Indices::conti0EqIdx;
@@ -499,13 +844,15 @@ diff -ruN exercises/exercise-coupling-ff-pm/models/porousmediumsubproblem.hh exe
faceEvaporation.push_back(flux[transportCompIdx]);
}
}
-@@ -345,12 +369,16 @@
+@@ -345,14 +369,16 @@
static const Scalar freeflowPressure = getParamFromGroup<Scalar>("Freeflow", "Problem.Pressure");
PrimaryVariables values(0.0);
-
- // TODO: dumux-course-task 2.A
- // Declare here which phases are present.
+- // TODO: dumux-course-task 2.C
+- // Change initial condition to 2p system with liquid saturation of 0.1
-
- values[Indices::pressureIdx] = freeflowPressure;
+#if EXNUMBER >= 3
@@ -521,7 +868,7 @@ diff -ruN exercises/exercise-coupling-ff-pm/models/porousmediumsubproblem.hh exe
return values;
}
-@@ -385,7 +413,12 @@
+@@ -387,7 +413,12 @@
{ return globalPos[1] > this->gridGeometry().bBoxMax()[1] - eps_; }
Scalar eps_;
@@ -535,8 +882,8 @@ diff -ruN exercises/exercise-coupling-ff-pm/models/porousmediumsubproblem.hh exe
std::shared_ptr<CouplingManager> couplingManager_;
diff -ruN exercises/exercise-coupling-ff-pm/models/properties.hh exercises/solution/exercise-coupling-ff-pm/models/properties.hh
---- exercises/exercise-coupling-ff-pm/models/properties.hh 2024-05-29 14:31:49.243625213 +0200
-+++ exercises/solution/exercise-coupling-ff-pm/models/properties.hh 2024-05-27 14:47:39.392451042 +0200
+--- exercises/exercise-coupling-ff-pm/models/properties.hh 2024-07-10 09:56:51.776442001 +0200
++++ exercises/solution/exercise-coupling-ff-pm/models/properties.hh 2023-07-26 11:31:46.014485560 +0200
@@ -34,15 +34,17 @@
#include <dumux/material/fluidsystems/h2oair.hh>
@@ -624,271 +971,64 @@ diff -ruN exercises/exercise-coupling-ff-pm/models/properties.hh exercises/solut
template<class TypeTag>
struct EnableGridGeometryCache<TypeTag, TTag::FreeflowNC> { static constexpr bool value = true; };
-diff -ruN exercises/exercise-coupling-ff-pm/README.md exercises/solution/exercise-coupling-ff-pm/README.md
---- exercises/exercise-coupling-ff-pm/README.md 2024-05-21 14:15:07.141554755 +0200
-+++ exercises/solution/exercise-coupling-ff-pm/README.md 1970-01-01 01:00:00.000000000 +0100
-@@ -1,438 +0,0 @@
--# Exercise Coupling free flow/porous medium flow (DuMuX Course)
+diff -ruN exercises/exercise-coupling-ff-pm/models/readme.md exercises/solution/exercise-coupling-ff-pm/models/readme.md
+--- exercises/exercise-coupling-ff-pm/models/readme.md 2024-07-15 10:52:11.058859896 +0200
++++ exercises/solution/exercise-coupling-ff-pm/models/readme.md 1970-01-01 01:00:00.000000000 +0100
+@@ -1,118 +0,0 @@
+-## 2. Changing the porous medium model
-
--The aim of this exercise is to get familiar with setting up coupled free flow/porous medium flow problems.
+-In this part of the exercise, the coupled system will be extended such that the transport of components
+-in both domains is included and the presence of a second phase in the porous medium domain is considered.
+-This enables the simulation of the drying of a porous medium (however no energy balance is included yet).
-
--## Problem set-up
+-Compared to exercise `1. Changing the interface` account now for the gravity effect as well as solving for the transient problem.
-
--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.
+-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.)
-
--### 0. Getting familiar with the code
-
--* Navigate to the directory `exercises/exercise-coupling-ff-pm`
+-### Task A: Change the model
-
--There are three sub folders: `interface` (Exercise 1), `models` (Exercise 2) and `turbulence` (Exercise 3).
+-In the first task, the porous-medium model will be changed from a 1p2c system to a 2p2c system.
+-Although a 2p2c system is used, 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 (`models/properties.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.
+-* Since two phases are involved now, we do not need to use the `OnePAdapter` anymore. Change the property of the `FluidSystem` such that `H2OAir` is used directly.
+- Afterwards, set the `transportCompIdx` to `Indices::switchIdx` in `porousmediumsubproblem.hh`.
-
--The task-related files for the simulation exercises are as follows:
--* Three __main files__ for the three sub-tasks :`interface/main.cc`, `models/main.cc`, `turbulence/main.cc`,
--* Three __free flow problem files__: `interface/freeflowsubproblem.hh`, `models/freeflowsubproblem.hh`, `turbulence/freeflowsubproblem.hh`
--* Three __porous medium flow problem files__: `interface/porousmediumsubproblem.hh`, `models/porousmediumsubproblem.hh`, `turbulence/porousmediumsubproblem.hh`
--* Three __properties files__: `interface/properties.hh`, `models/properties.hh`, `turbulence/properties.hh`
--* The __input files__: `interface/params.input`, `models/params.input`, `turbulence/params.input`,
--* The __spatial parameters files__: `1pspatialparams.hh`, `2pspatialparams.hh`
+-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
-
--In the main file, `TypeTags` for both submodels are defined, `FreeflowTypeTag` and `PorousMediumTypeTag`. These `TypeTags` collect all of the properties associated with each subproblem.
--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`, which help to assemble and solve the full coupled problem.
+-```
+-template<class TypeTag>
+-struct Formulation<TypeTag, TTag::PorousMediumOnePNC>
+-{ static constexpr auto value = TwoPFormulation::p1s0; };
+-```
+- in the properties file.
-
--The problem files very much look like the "regular", uncoupled problem files seen in previous exercises, with the exception that they hold a pointer to the `CouplingManager`.
--This allows them to evaluate the coupling conditions and to exchange information between the coupled models.
--The coupling conditions are realized technically in terms of boundary conditions. For instance,
--in `interface/freeflowsubproblem.hh`, `couplingNeumann` boundary conditions are set, which means that the free flow model evaluates the
--mass and momentum fluxes coming from the porous domain and uses these values as boundary conditions at the interface.
+-In contrast to the formulation, which remains the same during one simulation,
+-the primary variables may vary during one simulation.
+-In the case under investigation, 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 the 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()`. `MY_PHASE_PRESENCE` should be replaced with the correct value for the case where only a gas-phase (second phase of the fluid system) is present. For finding out if a gas or liquid phase is the first of second phase of a fluidsystem (here `h20air`), take a look at your fluidsystem (e.g. [h2oair](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/blob/master/dumux/material/fluidsystems/h2oair.hh?ref_type=heads#L75) ).
+- Moreover, have a look at the [indices.hh](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/blob/master/dumux/porousmediumflow/2pnc/indices.hh?ref_type=heads#L22) in the `2pnc` model (as the 2p2c model is a special case of the 2pnc model) in the subfolder `porousmediumflow` in your DuMuX directory to see which value to set for the gas-phase (here: second phase of our fluidsystem).
-
--Note that 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 to discretize the free-flow domain and a cell-centered finite volume method for the porous medium domain.
--Keep in mind that the staggered grid implementation distinguishes between face variables (velocity components) and
--cell center variables (all other variables).
--For this reason one distinguishes between `CouplingManager::stokesCellCenterIdx` and `CouplingManager::stokesFaceIdx` indices (see `main.cc`), while for the porous medium all variables can be accessed with `CouplingManager::darcyIdx`.
--
--__Task__:
--Take a closer look at the above listed 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 `interface/freeflowsubproblem.hh` and navigate to the part, 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.
--You can also delete the initial vertical velocity from `initialAtPos()`, to be consistent with the new boundary conditions.
--
--__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.setBeaversJoseph(Indices::momentumXBalanceIdx);
--```
--
--at the position where the coupling boundary conditions are set in `interface/freeflowsubproblem.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 the following line in `main.cc` in the subfolder `interface`.
--```cpp
--freeflowVtkWriter.addField(freeflowProblem->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: Change the shape of interface__
--
--Now we want to include a non-flat interface between the two domains.
--We use [`dune-subgrid`](https://doi.org/10.1007/s00607-009-0067-2) to construct two grids for the two domains from one common host grid. Thus, for the following tasks subgrid needs to be correctly installed.
--Our hostgrid will be a Dune-Yasp grid (`Dune::YaspGrid<2, Dune::TensorProductCoordinates<double, dim> >`)
--and the bounds and resolution of the domain will be set in the `params.input`file under the group `[Grid]`.
--This hostgrid, along with `elementSelector` functions defining some spatial cut of this domain, are passed to the grid manager to create each subdomain grid.
--
--To introduce this, open `main.cc` in the subfolder `interface` again and search for `TODO: dumux-course-task 1.C`.
--Comment out the first code block and uncomment the second.
--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 halves, separated by a sinusoidal interface.
--
--```cpp
--auto elementSelectorFreeflow = [&](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 your code is compiling and running. You can have a look at the results with paraview. How do the mass or mole fractions of water in air change over time?
-
--auto elementSelectorPorousMedium = [&](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 line for including the grid manager in the properties file, i.e.
--```cpp
--#include <dumux/io/grid/gridmanager_sub.hh>
--```
--
--and do the changes in the respective lines for the `Grid` property.
--
--The problem should now compile. However, a runtime 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());
-- ```
--in `freeflowsubproblem.hh` in the subfolder `interface`.
--
--The same is true for the BJS condition, however, here we need to consider the tangential direction:
--```cpp
--values.setBeaversJoseph(1 - scvf.directionIndex());
--```
--
--The final result should look something like this:
--
--
--
--*Extra Points:*
--Rather than enforcing a pressure difference across the domain, an inflow velocity profile could be set.
--What changes to the left boundary conditions in the free-flow domain would you make to introduce this? What conditions can be enforced on the right boundary?
--Hint: A relation between velocity and position is used for the vertical velocity component in the original form of the `initialAtPos` method.
--
--### 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 domains is included and the presence of a second phase in the porous medium domain 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 (`models/properties.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.
--* Since two phases are involved now, we do not need to use the `OnePAdapter` anymore. Change the property of the `FluidSystem` such that `H2OAir` is used directly.
-- Afterwards, set the `transportCompIdx` to `Indices::switchIdx` in `porousmediumsubproblem.hh`.
--
--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::PorousMediumOnePNC>
--{ static constexpr auto value = TwoPFormulation::p1s0; };
--```
-- in the properties 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 under investigation, 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()`. `MY_PHASE_PRESENCE` should be replaced with the correct value for the case where only a gas-phase (second phase of the fluid system) is present.
-- Have a look at the `indices.hh` in the `2pnc` model (as the 2p2c model is a special case of the 2pnc model) in the subfolder `porousmediumflow` in your DuMuX directory. (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__:
+-### Task B: Add output
-
-In the next step, we want to add some output to the simulation. The standard method for providing simulation output for visualization is via a `VtkWriter`. These tools take the grid geometries of each domain, as well as the solutions and write spatial output that one can view in visualization tools such as paraview.
-
@@ -905,6 +1045,7 @@ diff -ruN exercises/exercise-coupling-ff-pm/README.md exercises/solution/exercis
-* Have a look at the function `evaluateWaterMassStorageTerm()` in the `porousmediumsubproblem.hh`.
- Then implement a calculation of the total water mass:
- $`M^\text{w}:=\sum_{\alpha \in \textrm{g,l}} \left( \phi S_\alpha X^\text{w}_\alpha \varrho_\alpha V_\textrm{scv} \right)`$.
+- Hint: the volume is a property of a subcontrolvolume (`scv`), the other properties can be obtained through the volume variables (`volVars`). The functions for the `volVars` can be found in the header [dumux/porousmediumflow/2pnc/volumevariables.hh](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/blob/master/dumux/porousmediumflow/2pnc/volumevariables.hh?ref_type=heads)
-* Afterwards, adapt the method `startStorageEvaluation()` such that the variable `initialWaterContent_` is
- initialized correctly using the `evaluateWaterMassStorageTerm()` method and assign that value also to the variable `lastWaterMass_`.
-
@@ -918,15 +1059,19 @@ diff -ruN exercises/exercise-coupling-ff-pm/README.md exercises/solution/exercis
-Second, 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. Multiply this with the relevant face areas, extrusion factor, mass fraction, and seconds per day to get the [mm/d] units as seen in the storage evaluation.
+- for each coupling scvf. Multiply this with the relevant face areas, extrusion factor, molar Mass, and seconds per day to get the [mm/d] units as seen in the storage evaluation.
-* When the `[Problem.ExportFluxes] = true` parameter is enabled, simulation data will be exported for this simulation to a `.json` file.
- This file can flexibly handle data of different types, which in this case is helpful as we have both temporal and spatial data.
-* Use the python plotting script `plotFluxes.py` to visualize the flux distribution across the surface at different times.
-- This script uses `matplotlib`, a very popular python based visualization library.
+- This script uses `matplotlib`, a very popular python based visualization library and reads-in the data stored in the file `flux_models_coupling.json`.
+- You can execute the script by
+- ```
+- python3 ./plotFluxes.py
+- ```
-
-Compile and run the simulation and take a look at the results.
-
--__Task C: Let it dry out__:
+-### 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.
@@ -938,9 +1083,8 @@ diff -ruN exercises/exercise-coupling-ff-pm/README.md exercises/solution/exercis
-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 adapting
-- the parameter for `pcLowSw` and `pcHighSw` affects the results.
+-[Mosthaf (2014)](http://dx.doi.org/10.18419/opus-519). The regularization of the capillary pressure-saturation relationship with using the vanGenuchten model can be found in the file [dumux/material/fluidmatrixinteractions/2p/vangenuchten.hh](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/blob/master/dumux/material/fluidmatrixinteractions/2p/vangenuchten.hh?ref_type=heads#L304).
+-* See how adapting the parameter for `VanGenuchtenPcLowSweThreshold` and `VanGenuchtenPcHighSweThreshold` 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,
@@ -948,151 +1092,178 @@ diff -ruN exercises/exercise-coupling-ff-pm/README.md exercises/solution/exercis
-[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.
--
+-Now you are able to simulate a complete drying of the porous medium. Have a look the resulting liquid saturation distribution within the porous medium (with using paraview).
+\ No newline at end of file
+diff -ruN exercises/exercise-coupling-ff-pm/README.md exercises/solution/exercise-coupling-ff-pm/README.md
+--- exercises/exercise-coupling-ff-pm/README.md 2024-07-15 10:49:12.269896347 +0200
++++ exercises/solution/exercise-coupling-ff-pm/README.md 1970-01-01 01:00:00.000000000 +0100
+@@ -1,125 +0,0 @@
+-# Exercise Coupling free flow/porous medium flow (DuMuX Course)
-
--### 3. Use a turbulence model in the free flow domain
+-The aim of this exercise is to get familiar with setting up coupled free flow/porous medium flow problems.
-
--Several RANS turbulence models are implemented in DuMu<sup>x</sup>.
--This part of the exercise consists of the following steps:
--* replacing the Navier-Stokes model by the K-Omega SST turbulence model,
--* switching to a symmetry boundary condition,
--* applying a grid refinement towards the interface,
--* subsequently refining the grid (convergence study).
+-## Problem set-up
-
--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`.
+-> Note: The following problem setup holds for all sub-exercises 1-3.
+->
+-The model domain consists of two non-overlapping two-dimensional subdomains.
+-Free flow is modeled in the upper subdomain, while the lower subdomain models a flow within a porous medium.
+-Both single-phase flow and two-phase flow will be considered in the porous domain.
-
--__Task A__:
-
--The file `freeflowsubproblem.hh` is your free flow problem file and `properties.hh` is the properties file within this exercise.
+-## 0. Getting familiar with the code
-
--For using the compositional zero equation turbulence model, the following header files need to be included
-- in properties file:
--```
--#include <dumux/freeflow/compositional/sstncmodel.hh>
--```
--and in problem file:
--```
--#include <dumux/freeflow/turbulenceproperties.hh>
--#include <dumux/freeflow/rans/twoeq/sst/problem.hh>
--#include <dumux/freeflow/rans/boundarytypes.hh>
--```
+-* Navigate to the directory `exercises/exercise-coupling-ff-pm`
-
--The includes for the NavierStokesNC model and the NavierStokesProblem are no longer needed and can be removed.
+-There are three sub folders: `interface` (contains Exercise 1), `models` (contains Exercise 2) and `turbulence` (contains Exercise 3).
-
--Make sure your free flow problem inherits from the correct parent type:
--* Change the entry in the `FreeflowModel` definition accordingly (multi-component non-isothermal K-Omega SST model, `SSTNCNI`) in the properties file,
--* Adapt the inheritance of the problem class in problem file (Use `RANSProblem<TypeTag>` rather than `NavierStokesStaggeredProblem<TypeTag>`).
+-The folders of the three exercises contain the following files:
+-* A __main file__ (`interface/main.cc`, `models/main.cc`, `turbulence/main.cc`),
+-* one __problem file for the free-flow domain__ (`interface/freeflowsubproblem.hh`, `models/freeflowsubproblem.hh`, `turbulence/freeflowsubproblem.hh`),
+-* one __problem file for the porous medium domain__ (`interface/porousmediumsubproblem.hh`, `models/porousmediumsubproblem.hh`, `turbulence/porousmediumsubproblem.hh`),
+-* one __properties file__ (`interface/properties.hh`, `models/properties.hh`, `turbulence/properties.hh`),
+-* and one __input files__ (`interface/params.input`, `models/params.input`, `turbulence/params.input`).
+-Moreover all the exercises share
+-* the __spatial parameters files__ (`1pspatialparams.hh` and `2pspatialparams.hh`)
-
--Here, the turbulent free flow is wall-bounded, meaning shear stress and turbulence in the flow develop primarily due to the walls.
--The porous medium at the bottom and also the channel wall (upper boundary) act here as such walls.
--Within the problem the location of boundaries representing walls need to be defined.
--To do this, add the following function to the upper and lower boundaries within the `boundaryTypes` function in `freeflowsubproblem.hh`:
--```cpp
-- values.setWall();
+-In the [Exercise Mainfiles](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux-course/-/blob/feature/updateCouplingFFPMExercise/exercises/exercise-mainfile/README.md) the overall structure of a main-file was already introduced. For the coupled setup we need now to define properties related to one subproblem (free-flow or porous medium flow) for each of the subproblems.
+-E.g. in the main file, `TypeTags` for both submodels are defined, `FreeflowTypeTag` and `PorousMediumTypeTag`. These `TypeTags` collect all of the properties associated with each subproblem.
+-```c++
+- // Define the sub problem type tags
+- using FreeflowTypeTag = Properties::TTag::FreeflowNC;
+- using PorousMediumTypeTag = Properties::TTag::PorousMediumOnePNC;
-```
--
--With the locations of the wall designated, and the problem initialized, a series of constant spatial properties,
--in particular the distance to the nearest wall from each cell center, should be initialized.
--To do this, add the following to the `main.cc` file.
--```cpp
-- freeflowProblem->updateStaticWallProperties();
+-The same applies for types such as `GridManager`,
+-```c++
+- // try to create a grid (from the given grid file or the input file)
+- // for both sub-domains
+- using PorousMediumGridManager = Dumux::GridManager<GetPropType<PorousMediumTypeTag, Properties::Grid>>;
+- PorousMediumGridManager porousMediumGridManager;
+- porousMediumGridManager.init("PorousMedium"); // pass parameter group
+-
+- using FreeflowGridManager = Dumux::GridManager<GetPropType<FreeflowTypeTag, Properties::Grid>>;
+- FreeflowGridManager freeflowGridManager;
+- freeflowGridManager.init("Freeflow"); // pass parameter group
+-
+- // we compute on the leaf grid view
+- const auto& porousMediumGridView = porousMediumGridManager.grid().leafGridView();
+- const auto& freeflowGridView = freeflowGridManager.grid().leafGridView();
-```
--
--In addition, there is also a non-local solution-dependent aspect of the turbulence models which is to be updated at the end of each step.
--An example of this is a stencil extended velocity gradient, and other flow properties in relation to the wall.
--These dynamic interactions are to be initialized in the mainfile directly after the `updateStaticWallProperties()` call,
--as well as in the time loop after `// Update dynamic wall properties`:
--```cpp
--freeflowProblem->updateDynamicWallProperties(freeflowSol);
+-`FVGridGeometry`,
+-```c++
+- // create the finite volume grid geometry
+- using FreeflowFVGridGeometry = GetPropType<FreeflowTypeTag, Properties::GridGeometry>;
+- auto freeflowFvGridGeometry = std::make_shared<FreeflowFVGridGeometry>(freeflowGridView);
+- using PorousMediumFVGridGeometry = GetPropType<PorousMediumTypeTag, Properties::GridGeometry>;
+- auto porousMediumFvGridGeometry = std::make_shared<PorousMediumFVGridGeometry>(porousMediumGridView);
+-```
+-`Problem`, etc...
+-```c++
+- // the problem (initial and boundary conditions)
+- using FreeflowProblem = GetPropType<FreeflowTypeTag, Properties::Problem>;
+- auto freeflowProblem = std::make_shared<FreeflowProblem>(freeflowFvGridGeometry, couplingManager);
+- using PorousMediumProblem = GetPropType<PorousMediumTypeTag, Properties::Problem>;
+- auto porousMediumProblem = std::make_shared<PorousMediumProblem>(porousMediumFvGridGeometry, couplingManager);
-```
-
--In addition to designating the locations of walls,
--additional boundary conditions and initial conditions need to be set for the two new primary variables $k$ and $\omega$.
--In the `boundaryTypes` function, set both variables on all walls to be dirichlet, except for the right boundary, which should have outflow conditions.
--
--For the initial conditions, Reynolds number specific base conditions should be applied everywhere.
--In the problem constructor, uncomment the code calculating these terms,
--then apply the `turbulentKineticEnergy_`and `dissipation_` variables to their primary variables in all locations.
--Within the dirichlet function for cell faces (`dirichlet(element, scvf)`),
--we also need to specify that these variables should be fixed to 0 at the wall.
--
--In addition, dirichlet cell constraints for the dissipation or $\omega$ variable will be set for all wall adjacent cells.
--This is done in the `isDirichletCell` function, as well as the `dirichlet` function already, and requires no further changes.
--
--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__:
+-Since we use a monolithic coupling scheme, there is only one `Assembler` and one `NewtonSolver`, which help to assemble and solve the full coupled problem.
-
--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 at the top with
+-The problem files very much look like the "regular", uncoupled problem files seen in previous exercises, with the exception that they hold a pointer to the `CouplingManager`. E.g.
-```c++
--values.setAllSymmetry();
+- FreeFlowSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry,std::shared_ptr<CouplingManager> couplingManager)
+- : ParentType(fvGridGeometry, "Freeflow"),
+- eps_(1e-6),
+- couplingManager_(couplingManager)
+- {
+- //...
+- }
-```
+-This allows them to evaluate the coupling conditions and to exchange information between the coupled models.
+-The coupling conditions are realized technically in terms of boundary conditions. For instance,
+-in `interface/freeflowsubproblem.hh`, `couplingNeumann` boundary conditions are set, which means that the free flow model evaluates the
+-mass and momentum fluxes coming from the porous domain and uses these values as boundary conditions at the interface.
+-```c++
+- BoundaryTypes boundaryTypes(const Element& element,
+- const SubControlVolumeFace& scvf) const
+- {
+- BoundaryTypes values;
+- const auto& globalPos = scvf.center();
+-
+- //...
+- // 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
+- }
-
--In addition, you have to remove the condition `onUpperBoundary_(globalPos)` from the `initialAtPos(globalPos)` method and the `dirichlet(element, scvf)` method.
--
--__Task C__:
+- return values;
+- }
+-```
+-Note that 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.
-
--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.
+-We will use a staggered grid (also calles Marker-and-Cell method - MAC)to discretize the free-flow domain and a cell-centered finite volume (CCFV) method for the porous medium domain.
+-Keep in mind that the staggered grid implementation distinguishes between face variables (velocity components) and cell center variables (all other variables).
+-For this reason one distinguishes between `CouplingManager::stokesCellCenterIdx` and `CouplingManager::stokesFaceIdx` indices (see `main.cc`), while for the porous medium all variables can be accessed with `CouplingManager::darcyIdx`.
-
--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
+- // the indices
+- constexpr auto freeflowCellCenterIdx = CouplingManager::stokesCellCenterIdx;
+- constexpr auto freeflowFaceIdx = CouplingManager::stokesFaceIdx;
+- constexpr auto porousMediumIdx = CouplingManager::darcyIdx;
-```
-
--__Task D__:
+-__Task__:
+-Take a closer look at the above listed files before moving to the three exercises below.
-
--For the grid convergence study, run various simulations with the following grading parameters:
--```c++
--* [Freeflow.Grid] Grading1 = 1.2, [PorousMedium.Grid] Grading1 = -1.2
--* [Freeflow.Grid] Grading1 = 1.3, [PorousMedium.Grid] Grading1 = -1.3
--* [Freeflow.Grid] Grading1 = 1.4, [PorousMedium.Grid] Grading1 = -1.4
--* [Freeflow.Grid] Grading1 = 1.5, [PorousMedium.Grid] Grading1 = -1.5
--* [Freeflow.Grid] Grading1 = 1.6, [PorousMedium.Grid] Grading1 = -1.6
--```
+-## Sub-Exercises
-
--By changing the parameter `Problem.Name` for each grading factor, you avoid losing 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.
+-* [**Exercise 1:** Changing the interface between the free- and the porous medium domain](./interface/readme.md)
+-* [**Exercise 2:** Changing the porous medium model](./models/readme.md)
+-* [**Exercise 3:** Introducing a turbulence model in the free flow domain](./turbulence/readme.md)
+\ No newline at end of file
diff -ruN exercises/exercise-coupling-ff-pm/turbulence/CMakeLists.txt exercises/solution/exercise-coupling-ff-pm/turbulence/CMakeLists.txt
---- exercises/exercise-coupling-ff-pm/turbulence/CMakeLists.txt 2024-05-21 14:15:07.141554755 +0200
-+++ exercises/solution/exercise-coupling-ff-pm/turbulence/CMakeLists.txt 2024-05-21 14:15:07.209555087 +0200
-@@ -1,7 +1,17 @@
+--- exercises/exercise-coupling-ff-pm/turbulence/CMakeLists.txt 2023-04-04 12:15:36.845203367 +0200
++++ exercises/solution/exercise-coupling-ff-pm/turbulence/CMakeLists.txt 2024-07-15 10:49:12.273896458 +0200
+@@ -1,7 +1,30 @@
-# executables for ex_interface_coupling_ff-pm
-dumux_add_test(NAME exercise_turbulence_coupling_ff-pm
+dumux_add_test(NAME exercise_turbulence_coupling_ff-pm_original
SOURCES main.cc
- LABELS ffpm)
+ LABELS ffpm
-+ COMPILE_DEFINITIONS EXNUMBER=0)
++ COMPILE_DEFINITIONS EXNUMBER=0
++ COMMAND ./exercise_turbulence_coupling_ff-pm_original
++ CMD_ARGS params_orig_a.input)
+
+dumux_add_test(NAME exercise_turbulence_coupling_ff-pm_a_solution
+ SOURCES main.cc
+ LABELS ffpm
-+ COMPILE_DEFINITIONS EXNUMBER=1)
++ COMPILE_DEFINITIONS EXNUMBER=1
++ COMMAND ./exercise_turbulence_coupling_ff-pm_a_solution
++ CMD_ARGS params_orig_a.input)
+
+dumux_add_test(NAME exercise_turbulence_coupling_ff-pm_b_solution
+ SOURCES main.cc
+ LABELS ffpm
-+ COMPILE_DEFINITIONS EXNUMBER=2)
++ COMPILE_DEFINITIONS EXNUMBER=2
++ COMMAND ./exercise_turbulence_coupling_ff-pm_b_solution
++ CMD_ARGS params_b.input)
++
++dumux_add_test(NAME exercise_turbulence_coupling_ff-pm_c_d_solution
++ SOURCES main.cc
++ LABELS ffpm
++ COMPILE_DEFINITIONS EXNUMBER=3
++ COMMAND ./exercise_turbulence_coupling_ff-pm_c_d_solution
++ CMD_ARGS params_c_d.input)
# add a symlink for each input file
add_input_file_links()
diff -ruN exercises/exercise-coupling-ff-pm/turbulence/freeflowsubproblem.hh exercises/solution/exercise-coupling-ff-pm/turbulence/freeflowsubproblem.hh
---- exercises/exercise-coupling-ff-pm/turbulence/freeflowsubproblem.hh 2024-05-29 14:31:49.243625213 +0200
-+++ exercises/solution/exercise-coupling-ff-pm/turbulence/freeflowsubproblem.hh 2024-05-21 14:15:07.209555087 +0200
+--- exercises/exercise-coupling-ff-pm/turbulence/freeflowsubproblem.hh 2024-07-12 14:03:42.757939741 +0200
++++ exercises/solution/exercise-coupling-ff-pm/turbulence/freeflowsubproblem.hh 2024-01-03 15:09:59.469499202 +0100
@@ -20,8 +20,9 @@
* \file
* \brief The free-flow sub problem
@@ -1308,8 +1479,8 @@ diff -ruN exercises/exercise-coupling-ff-pm/turbulence/freeflowsubproblem.hh exe
return values;
}
diff -ruN exercises/exercise-coupling-ff-pm/turbulence/main.cc exercises/solution/exercise-coupling-ff-pm/turbulence/main.cc
---- exercises/exercise-coupling-ff-pm/turbulence/main.cc 2024-05-21 14:15:07.145554774 +0200
-+++ exercises/solution/exercise-coupling-ff-pm/turbulence/main.cc 2024-05-21 14:15:07.209555087 +0200
+--- exercises/exercise-coupling-ff-pm/turbulence/main.cc 2024-07-12 14:03:39.469855701 +0200
++++ exercises/solution/exercise-coupling-ff-pm/turbulence/main.cc 2023-07-26 11:31:46.014485560 +0200
@@ -133,12 +133,12 @@
auto solOld = sol;
@@ -1342,19 +1513,314 @@ diff -ruN exercises/exercise-coupling-ff-pm/turbulence/main.cc exercises/solutio
// post time step treatment of PorousMedium problem
porousMediumProblem->postTimeStep(sol[porousMediumIdx], *porousMediumGridVariables, timeLoop->timeStepSize());
+diff -ruN exercises/exercise-coupling-ff-pm/turbulence/params_b.input exercises/solution/exercise-coupling-ff-pm/turbulence/params_b.input
+--- exercises/exercise-coupling-ff-pm/turbulence/params_b.input 1970-01-01 01:00:00.000000000 +0100
++++ exercises/solution/exercise-coupling-ff-pm/turbulence/params_b.input 2024-07-15 10:49:12.273896458 +0200
+@@ -0,0 +1,72 @@
++[TimeLoop]
++DtInitial = 1e-1 # [s]
++MaxTimeStepSize = 43200 # [s] (12 hours)
++TEnd = 864000 # [s] (6 days)
++
++[Freeflow.Grid]
++Positions0 = 0.0 0.25
++Positions1 = 0.25 0.375
++Grading0 = 1.0
++Grading1 = 1.0
++Cells0 = 15
++Cells1 = 10
++Verbosity = true
++
++[PorousMedium.Grid]
++Positions0 = 0.0 0.25
++Positions1 = 0.0 0.25
++Cells0 = 15
++Cells1 = 10
++Grading0 = 1.0
++Grading1 = 1.0
++Verbosity = true
++
++[Freeflow.Problem]
++Name = freeflow
++RefVelocity = 3.5 # [m/s]
++RefPressure = 1e5 # [Pa]
++refMoleFrac = 0 # [-]
++RefTemperature = 298.15 # [K]
++EnableInertiaTerms = true
++
++[PorousMedium.Problem]
++Name = porousmedium
++Pressure = 1e5
++Saturation = 0.5 # initial Sw
++Temperature = 298.15 # [K]
++InitPhasePresence = 3 # bothPhases
++
++[SpatialParams]
++Porosity = 0.41
++Permeability = 2.65e-10
++AlphaBeaversJoseph = 1.0
++VanGenuchtenN = 6.9
++VanGenuchtenAlpha = 6.371e-4
++Swr = 0.005
++Snr = 0.01
++Temperature = 298.15 # [K]
++
++[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
++
++[RANS]
++IsFlatWallBounded = True
+diff -ruN exercises/exercise-coupling-ff-pm/turbulence/params_c_d.input exercises/solution/exercise-coupling-ff-pm/turbulence/params_c_d.input
+--- exercises/exercise-coupling-ff-pm/turbulence/params_c_d.input 1970-01-01 01:00:00.000000000 +0100
++++ exercises/solution/exercise-coupling-ff-pm/turbulence/params_c_d.input 2024-07-15 10:49:12.273896458 +0200
+@@ -0,0 +1,72 @@
++[TimeLoop]
++DtInitial = 1e-1 # [s]
++MaxTimeStepSize = 43200 # [s] (12 hours)
++TEnd = 864000 # [s] (6 days)
++
++[Freeflow.Grid]
++Positions0 = 0.0 0.25
++Positions1 = 0.25 0.375
++Grading0 = 1.0
++Grading1 = 1.2
++Cells0 = 15
++Cells1 = 10
++Verbosity = true
++
++[PorousMedium.Grid]
++Positions0 = 0.0 0.25
++Positions1 = 0.0 0.25
++Cells0 = 15
++Cells1 = 10
++Grading0 = 1.0
++Grading1 = -1.2
++Verbosity = true
++
++[Freeflow.Problem]
++Name = freeflow
++RefVelocity = 3.5 # [m/s]
++RefPressure = 1e5 # [Pa]
++refMoleFrac = 0 # [-]
++RefTemperature = 298.15 # [K]
++EnableInertiaTerms = true
++
++[PorousMedium.Problem]
++Name = porousmedium
++Pressure = 1e5
++Saturation = 0.5 # initial Sw
++Temperature = 298.15 # [K]
++InitPhasePresence = 3 # bothPhases
++
++[SpatialParams]
++Porosity = 0.41
++Permeability = 2.65e-10
++AlphaBeaversJoseph = 1.0
++VanGenuchtenN = 6.9
++VanGenuchtenAlpha = 6.371e-4
++Swr = 0.005
++Snr = 0.01
++Temperature = 298.15 # [K]
++
++[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
++
++[RANS]
++IsFlatWallBounded = True
diff -ruN exercises/exercise-coupling-ff-pm/turbulence/params.input exercises/solution/exercise-coupling-ff-pm/turbulence/params.input
---- exercises/exercise-coupling-ff-pm/turbulence/params.input 2024-05-21 14:15:07.145554774 +0200
-+++ exercises/solution/exercise-coupling-ff-pm/turbulence/params.input 2024-05-21 14:15:07.209555087 +0200
-@@ -69,3 +69,6 @@
- SolidDensity = 2700
- SolidThermalConductivity = 2.8
- SolidHeatCapacity = 790
+--- exercises/exercise-coupling-ff-pm/turbulence/params.input 2024-07-15 10:49:12.273896458 +0200
++++ exercises/solution/exercise-coupling-ff-pm/turbulence/params.input 1970-01-01 01:00:00.000000000 +0100
+@@ -1,73 +0,0 @@
+-[TimeLoop]
+-DtInitial = 1e-1 # [s]
+-MaxTimeStepSize = 43200 # [s] (12 hours)
+-TEnd = 864000 # [s] (6 days)
+-
+-[Freeflow.Grid]
+-Positions0 = 0.0 0.25
+-# TODO: dumux-course-task 3.B - use only half of FF domain height
+-Positions1 = 0.25 0.5
+-# TODO: dumux-course-task 3.C - refine towards interface
+-Grading0 = 1.0
+-Grading1 = 1.0
+-Cells0 = 15
+-# TODO: dumux-course-task 3.B - use only half of FF domain height and adapt cell number
+-Cells1 = 20
+-Verbosity = true
+-
+-[PorousMedium.Grid]
+-Positions0 = 0.0 0.25
+-Positions1 = 0.0 0.25
+-Cells0 = 15
+-Cells1 = 10
+-# TODO: dumux-course-task 3.C - refine towards interface
+-Grading0 = 1.0
+-Grading1 = 1.0
+-Verbosity = true
+-
+-[Freeflow.Problem]
+-Name = freeflow
+-RefVelocity = 3.5 # [m/s]
+-RefPressure = 1e5 # [Pa]
+-refMoleFrac = 0 # [-]
+-RefTemperature = 298.15 # [K]
+-EnableInertiaTerms = true
+-
+-[PorousMedium.Problem]
+-Name = porousmedium
+-Pressure = 1e5
+-Saturation = 0.5 # initial Sw
+-Temperature = 298.15 # [K]
+-InitPhasePresence = 3 # bothPhases
+-
+-[SpatialParams]
+-Porosity = 0.41
+-Permeability = 2.65e-10
+-AlphaBeaversJoseph = 1.0
+-VanGenuchtenN = 6.9
+-VanGenuchtenAlpha = 6.371e-4
+-Swr = 0.005
+-Snr = 0.01
+-Temperature = 298.15 # [K]
+-
+-[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 -ruN exercises/exercise-coupling-ff-pm/turbulence/params_orig_a.input exercises/solution/exercise-coupling-ff-pm/turbulence/params_orig_a.input
+--- exercises/exercise-coupling-ff-pm/turbulence/params_orig_a.input 1970-01-01 01:00:00.000000000 +0100
++++ exercises/solution/exercise-coupling-ff-pm/turbulence/params_orig_a.input 2024-07-15 10:49:12.273896458 +0200
+@@ -0,0 +1,72 @@
++[TimeLoop]
++DtInitial = 1e-1 # [s]
++MaxTimeStepSize = 43200 # [s] (12 hours)
++TEnd = 864000 # [s] (6 days)
++
++[Freeflow.Grid]
++Positions0 = 0.0 0.25
++Positions1 = 0.25 0.5
++Grading0 = 1.0
++Grading1 = 1.0
++Cells0 = 15
++Cells1 = 20
++Verbosity = true
++
++[PorousMedium.Grid]
++Positions0 = 0.0 0.25
++Positions1 = 0.0 0.25
++Cells0 = 15
++Cells1 = 10
++Grading0 = 1.0
++Grading1 = 1.0
++Verbosity = true
++
++[Freeflow.Problem]
++Name = freeflow
++RefVelocity = 3.5 # [m/s]
++RefPressure = 1e5 # [Pa]
++refMoleFrac = 0 # [-]
++RefTemperature = 298.15 # [K]
++EnableInertiaTerms = true
++
++[PorousMedium.Problem]
++Name = porousmedium
++Pressure = 1e5
++Saturation = 0.5 # initial Sw
++Temperature = 298.15 # [K]
++InitPhasePresence = 3 # bothPhases
++
++[SpatialParams]
++Porosity = 0.41
++Permeability = 2.65e-10
++AlphaBeaversJoseph = 1.0
++VanGenuchtenN = 6.9
++VanGenuchtenAlpha = 6.371e-4
++Swr = 0.005
++Snr = 0.01
++Temperature = 298.15 # [K]
++
++[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
+
+[RANS]
+IsFlatWallBounded = True
diff -ruN exercises/exercise-coupling-ff-pm/turbulence/properties.hh exercises/solution/exercise-coupling-ff-pm/turbulence/properties.hh
---- exercises/exercise-coupling-ff-pm/turbulence/properties.hh 2024-05-29 14:31:49.243625213 +0200
-+++ exercises/solution/exercise-coupling-ff-pm/turbulence/properties.hh 2024-05-21 14:15:07.209555087 +0200
+--- exercises/exercise-coupling-ff-pm/turbulence/properties.hh 2024-07-12 14:03:39.469855701 +0200
++++ exercises/solution/exercise-coupling-ff-pm/turbulence/properties.hh 2023-07-26 11:31:46.014485560 +0200
@@ -20,8 +20,8 @@
* \file
* \brief The coupled exercise properties file or the turbulent case.
@@ -1395,3 +1861,133 @@ diff -ruN exercises/exercise-coupling-ff-pm/turbulence/properties.hh exercises/s
} // end namespace TTag
// Set the coupling manager
+diff -ruN exercises/exercise-coupling-ff-pm/turbulence/readme.md exercises/solution/exercise-coupling-ff-pm/turbulence/readme.md
+--- exercises/exercise-coupling-ff-pm/turbulence/readme.md 2024-07-15 10:49:12.273896458 +0200
++++ exercises/solution/exercise-coupling-ff-pm/turbulence/readme.md 1970-01-01 01:00:00.000000000 +0100
+@@ -1,126 +0,0 @@
+-## 3. Use a turbulence model in the free flow domain
+-As in the previous exercise, we account for the gravity effect as well as solving for the transient problem hereafter.
+-
+-Several RANS (Reynolds Averaged Navier Stokes) turbulence models are implemented in DuMu<sup>x</sup>.
+-This part of the exercise consists of the following steps:
+-* replacing the Navier-Stokes model by the [K-Omega SST (shear stress transport) turbulence model](https://www.cfd-online.com/Wiki/SST_k-omega_model) of [Menter (1994)](https://doi.org/10.2514/3.12149), which is a two-equation RANS 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: Use a $k$-$\omega$ turbulence model
+-
+-For using the compositional $k$-$\omega$ turbulence model, the following header files need to be included in properties file (`properties.hh`):
+-```
+-#include <dumux/freeflow/compositional/sstncmodel.hh>
+-```
+-and in free-flow problem file (`freeflowsubproblem.hh`):
+-```
+-#include <dumux/freeflow/turbulenceproperties.hh>
+-#include <dumux/freeflow/rans/twoeq/sst/problem.hh>
+-#include <dumux/freeflow/rans/boundarytypes.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 `FreeflowModel` definition accordingly (multi-component non-isothermal K-Omega SST model, `SSTNCNI`) in the properties file,
+-* Adapt the inheritance of the problem class in problem file (Use `RANSProblem<TypeTag>` rather than `NavierStokesStaggeredProblem<TypeTag>`).
+-
+-Additionally make sure that the alias `BoundaryTypes` uses the new boundary types for the RANS model (as included by the respective header above) instead of `NavierStokesBoundaryTypes`.
+-
+-Here, the turbulent free flow is wall-bounded, meaning shear stress and turbulence in the flow develop primarily due to the walls.
+-The porous medium at the bottom and also the channel wall (upper boundary) act here as such walls.
+-Within the problem the location of boundaries representing walls need to be defined.
+-To do this, add the following function to the upper and lower boundaries within the `boundaryTypes` function in `freeflowsubproblem.hh`:
+-```cpp
+- values.setWall();
+-```
+-
+-With the locations of the wall designated, and the problem initialized, a series of constant spatial properties,
+-in particular the distance to the nearest wall from each cell center, should be initialized.
+-To do this, add the following to the `main.cc` file.
+-```cpp
+- freeflowProblem->updateStaticWallProperties();
+-```
+-
+-In addition, there is also a non-local solution-dependent aspect of the turbulence models which is to be updated at the end of each step.
+-An example of this is a stencil extended velocity gradient, and other flow properties in relation to the wall.
+-These dynamic interactions are to be initialized in the mainfile directly after the `updateStaticWallProperties()` call,
+-as well as in the time loop after `// Update dynamic wall properties`:
+-```cpp
+-freeflowProblem->updateDynamicWallProperties(freeflowSol);
+-```
+-
+-In addition to designating the locations of walls,
+-additional boundary conditions and initial conditions need to be set for the two new primary variables, the turbulent kinetic energy $k$ and the specific tubulent dissipation rate/turbulence frequency $\omega$.
+-
+-In the `boundaryTypes` function, set both variables ($k$ and $\omega$) on all free-flow domain boundaries to be dirichlet, except for the right boundary, which should have outflow conditions. The name of the indices for those two variables that have to be used for setting the boundary conditions types can be found in `dumux/freeflow/rans/twoeq/indices.hh`.
+-
+-Within the dirichlet function for cell faces (`dirichlet(element, scvf)`),
+-we also need to specify that these variables ($k$ and $\omega$) should be fixed to 0 at the wall.
+-In addition, dirichlet cell constraints for the dissipation or $\omega$ variable will be set for all wall adjacent cells.
+-This is done in the `isDirichletCell` function, as well as the `dirichlet` function already, and requires no further changes.
+-
+-For the initial conditions, Reynolds number specific base conditions should be applied everywhere.
+-In the problem constructor, uncomment the code calculating these terms,
+-then apply the `turbulentKineticEnergy_`and `dissipation_` variables to their primary variables in the initial state for all spatial locations.
+-
+-Before starting to run your modified code make sure to include the parameter `Rans.IsFlatWallBounded` in your input file and set this parameter to `true`.
+-
+-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 diffusivity should look like this:
+-
+-
+-
+-### Task B: Use symmetry boundary conditions
+-
+-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 at the top with
+-```c++
+-values.setAllSymmetry();
+-```
+-
+-In addition, you have to remove the condition `onUpperBoundary_(globalPos)` from the `initialAtPos(globalPos)` method and the `dirichlet(element, scvf)` method.
+-
+-Now it is sufficient to only calulate half of the domain height for the free-flow domain. For this, adapt the grid coordinates as well as the number of cells in y-direction in your `params.input` - file accordingly.
+-After recompiling and running your simulation with the new boundary condition, check your results with paraview.
+-Does your velocity profile at the symmetry boundary look reasonable?
+-
+-Before closing paraview, Choose `Surface With Edges` instead of `Surface` in the Paraview toolbar to see the discretization grid. Until now we used a regular structured lattice grid.
+-
+-### Task C: Refine grid towards the interface
+-
+-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
+-```
+-
+-Rerun your simulation with the new grid and check your refined grid with paraview.
+-
+-### Task D: Perform a grid convergence study
+-
+-For the grid convergence study, run various simulations with the following grading parameters:
+-```c++
+-* [Freeflow.Grid] Grading1 = 1.2, [PorousMedium.Grid] Grading1 = -1.2
+-* [Freeflow.Grid] Grading1 = 1.3, [PorousMedium.Grid] Grading1 = -1.3
+-* [Freeflow.Grid] Grading1 = 1.4, [PorousMedium.Grid] Grading1 = -1.4
+-* [Freeflow.Grid] Grading1 = 1.5, [PorousMedium.Grid] Grading1 = -1.5
+-* [Freeflow.Grid] Grading1 = 1.6, [PorousMedium.Grid] Grading1 = -1.6
+-```
+-
+-By changing the parameter `Problem.Name` for each grading factor, you avoid losing the `.vtu` and `.pvd` files of the previous simulation runs.
+-Check the first lines of the output in your terminal to see how the grading factors change the height of your grid cells.
+-
+-Compare the velocity fields with Paraview. In Paraview you can e.g. use the tool `PlotOverLine` to plot the velocity in $x$-direction ($v_x$) over the oulet on the right hand side. After choosing `PlotOverLine` make sure to give the Point 1 as (0.25, 0.25, 0) and the Point 2 as (0.25, 0.5, 0) to define the line a the output. Apply this and choose to only visualize `velocity_gas(m/s)_X`. You can do this for all your refinements and plot them in the same diagram to compare.
diff --git a/exercises/exercise-coupling-ff-pm/README.md b/exercises/exercise-coupling-ff-pm/README.md
index 279dde34..7cc5a61d 100644
--- a/exercises/exercise-coupling-ff-pm/README.md
+++ b/exercises/exercise-coupling-ff-pm/README.md
@@ -4,435 +4,122 @@ The aim of this exercise is to get familiar with setting up coupled free flow/po
## 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.
+> Note: The following problem setup holds for all sub-exercises 1-3.
+>
+The model domain consists of two non-overlapping two-dimensional subdomains.
+Free flow is modeled in the upper subdomain, while the lower subdomain models a flow within a porous medium.
Both single-phase flow and two-phase flow will be considered in the porous domain.
-### 0. Getting familiar with the code
+## 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).
+There are three sub folders: `interface` (contains Exercise 1), `models` (contains Exercise 2) and `turbulence` (contains Exercise 3).
-The task-related files for the simulation exercises are as follows:
-* Three __main files__ for the three sub-tasks :`interface/main.cc`, `models/main.cc`, `turbulence/main.cc`,
-* Three __free flow problem files__: `interface/freeflowsubproblem.hh`, `models/freeflowsubproblem.hh`, `turbulence/freeflowsubproblem.hh`
-* Three __porous medium flow problem files__: `interface/porousmediumsubproblem.hh`, `models/porousmediumsubproblem.hh`, `turbulence/porousmediumsubproblem.hh`
-* Three __properties files__: `interface/properties.hh`, `models/properties.hh`, `turbulence/properties.hh`
-* The __input files__: `interface/params.input`, `models/params.input`, `turbulence/params.input`,
-* The __spatial parameters files__: `1pspatialparams.hh`, `2pspatialparams.hh`
+The folders of the three exercises contain the following files:
+* A __main file__ (`interface/main.cc`, `models/main.cc`, `turbulence/main.cc`),
+* one __problem file for the free-flow domain__ (`interface/freeflowsubproblem.hh`, `models/freeflowsubproblem.hh`, `turbulence/freeflowsubproblem.hh`),
+* one __problem file for the porous medium domain__ (`interface/porousmediumsubproblem.hh`, `models/porousmediumsubproblem.hh`, `turbulence/porousmediumsubproblem.hh`),
+* one __properties file__ (`interface/properties.hh`, `models/properties.hh`, `turbulence/properties.hh`),
+* and one __input files__ (`interface/params.input`, `models/params.input`, `turbulence/params.input`).
+Moreover all the exercises share
+* the __spatial parameters files__ (`1pspatialparams.hh` and `2pspatialparams.hh`)
-In the main file, `TypeTags` for both submodels are defined, `FreeflowTypeTag` and `PorousMediumTypeTag`. These `TypeTags` collect all of the properties associated with each subproblem.
-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`, which help to assemble and solve the full coupled problem.
-
-The problem files very much look like the "regular", uncoupled problem files seen in previous exercises, with the exception that they hold a pointer to the `CouplingManager`.
-This allows them to evaluate the coupling conditions and to exchange information between the coupled models.
-The coupling conditions are realized technically in terms of boundary conditions. For instance,
-in `interface/freeflowsubproblem.hh`, `couplingNeumann` boundary conditions are set, which means that the free flow model evaluates the
-mass and momentum fluxes coming from the porous domain and uses these values as boundary conditions at the interface.
-
-Note that 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 to discretize the free-flow domain and a cell-centered finite volume method for the porous medium domain.
-Keep in mind that the staggered grid implementation distinguishes between face variables (velocity components) and
-cell center variables (all other variables).
-For this reason one distinguishes between `CouplingManager::stokesCellCenterIdx` and `CouplingManager::stokesFaceIdx` indices (see `main.cc`), while for the porous medium all variables can be accessed with `CouplingManager::darcyIdx`.
-
-__Task__:
-Take a closer look at the above listed 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 `interface/freeflowsubproblem.hh` and navigate to the part, 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.
-You can also delete the initial vertical velocity from `initialAtPos()`, to be consistent with the new boundary conditions.
-
-__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.setBeaversJoseph(Indices::momentumXBalanceIdx);
-```
-
-at the position where the coupling boundary conditions are set in `interface/freeflowsubproblem.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 the following line in `main.cc` in the subfolder `interface`.
-```cpp
-freeflowVtkWriter.addField(freeflowProblem->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: Change the shape of interface__
-
-Now we want to include a non-flat interface between the two domains.
-We use [`dune-subgrid`](https://doi.org/10.1007/s00607-009-0067-2) to construct two grids for the two domains from one common host grid. Thus, for the following tasks subgrid needs to be correctly installed.
-Our hostgrid will be a Dune-Yasp grid (`Dune::YaspGrid<2, Dune::TensorProductCoordinates<double, dim> >`)
-and the bounds and resolution of the domain will be set in the `params.input`file under the group `[Grid]`.
-This hostgrid, along with `elementSelector` functions defining some spatial cut of this domain, are passed to the grid manager to create each subdomain grid.
-
-To introduce this, open `main.cc` in the subfolder `interface` again and search for `TODO: dumux-course-task 1.C`.
-Comment out the first code block and uncomment the second.
-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 halves, separated by a sinusoidal interface.
-
-```cpp
-auto elementSelectorFreeflow = [&](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 elementSelectorPorousMedium = [&](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 line for including the grid manager in the properties file, i.e.
-```cpp
-#include <dumux/io/grid/gridmanager_sub.hh>
-```
-
-and do the changes in the respective lines for the `Grid` property.
-
-The problem should now compile. However, a runtime 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());
- ```
-in `freeflowsubproblem.hh` in the subfolder `interface`.
-
-The same is true for the BJS condition, however, here we need to consider the tangential direction:
-```cpp
-values.setBeaversJoseph(1 - scvf.directionIndex());
-```
-
-The final result should look something like this:
-
-
-
-*Extra Points:*
-Rather than enforcing a pressure difference across the domain, an inflow velocity profile could be set.
-What changes to the left boundary conditions in the free-flow domain would you make to introduce this? What conditions can be enforced on the right boundary?
-Hint: A relation between velocity and position is used for the vertical velocity component in the original form of the `initialAtPos` method.
-
-### 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 domains is included and the presence of a second phase in the porous medium domain 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 (`models/properties.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.
-* Since two phases are involved now, we do not need to use the `OnePAdapter` anymore. Change the property of the `FluidSystem` such that `H2OAir` is used directly.
- Afterwards, set the `transportCompIdx` to `Indices::switchIdx` in `porousmediumsubproblem.hh`.
-
-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::PorousMediumOnePNC>
-{ static constexpr auto value = TwoPFormulation::p1s0; };
+In the [Exercise Mainfiles](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux-course/-/blob/feature/updateCouplingFFPMExercise/exercises/exercise-mainfile/README.md) the overall structure of a main-file was already introduced. For the coupled setup we need now to define properties related to one subproblem (free-flow or porous medium flow) for each of the subproblems.
+E.g. in the main file, `TypeTags` for both submodels are defined, `FreeflowTypeTag` and `PorousMediumTypeTag`. These `TypeTags` collect all of the properties associated with each subproblem.
+```c++
+ // Define the sub problem type tags
+ using FreeflowTypeTag = Properties::TTag::FreeflowNC;
+ using PorousMediumTypeTag = Properties::TTag::PorousMediumOnePNC;
```
- in the properties 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 under investigation, 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()`. `MY_PHASE_PRESENCE` should be replaced with the correct value for the case where only a gas-phase (second phase of the fluid system) is present.
- Have a look at the `indices.hh` in the `2pnc` model (as the 2p2c model is a special case of the 2pnc model) in the subfolder `porousmediumflow` in your DuMuX directory. (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. The standard method for providing simulation output for visualization is via a `VtkWriter`. These tools take the grid geometries of each domain, as well as the solutions and write spatial output that one can view in visualization tools such as paraview.
-
-Although this Vtk output is very useful, some output is more suited for other forms of visualization.
-Two examples of data output formats are `.csv` files and `.json` files.
-From these files, data can be visualized using many different tools.
-In this exercise, data exported to a `.csv` file will be plotted using gnuplot, and data exported to a `.json` output will be plotted using the matplotlib python library.
-
-First as example output data, we want to investigate the water mass in the (porous-medium) system.
-There are two ways to evaluate this: (1) the total storage of water mass in the porous medium domain, and (2) the water vapor flux along the interface.
-The total storage can be plotted over time, whereas the water vapor flux will vary spatially along the interface as well as over time.
-
-First, we evaluate the storage term of the water component.
-* Have a look at the function `evaluateWaterMassStorageTerm()` in the `porousmediumsubproblem.hh`.
- Then implement a calculation of the total water mass:
- $`M^\text{w}:=\sum_{\alpha \in \textrm{g,l}} \left( \phi S_\alpha X^\text{w}_\alpha \varrho_\alpha V_\textrm{scv} \right)`$.
-* Afterwards, adapt the method `startStorageEvaluation()` such that the variable `initialWaterContent_` is
- initialized correctly using the `evaluateWaterMassStorageTerm()` method and assign that value also to the variable `lastWaterMass_`.
-
-After these functions are correctly implemented, the evaporation in [mm/d] is calculated in the code as
-$`e = \frac{\textrm{d}\, M^\text{w}}{\textrm{d}\, t} / A_\textrm{interface}`$.
-(all boundaries have Neumann no-flow conditions and no sources are present, meaning the water can only leave the system
-via the porous-medium free-flow interface.)
-
-During the simulation a `.csv` file is created (enable by setting `[Problem.ExportStorage] = true` in the `params.input` file). In addition, a gnuplot script `StorageOverTime.gp` is written and the mass loss and cumulative mass are plotted over time (`StorageOverTime.png`), when plotting is enabled with `[Problem.PlotStorage] = true`.
-
-Second, 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. Multiply this with the relevant face areas, extrusion factor, mass fraction, and seconds per day to get the [mm/d] units as seen in the storage evaluation.
-* When the `[Problem.ExportFluxes] = true` parameter is enabled, simulation data will be exported for this simulation to a `.json` file.
- This file can flexibly handle data of different types, which in this case is helpful as we have both temporal and spatial data.
-* Use the python plotting script `plotFluxes.py` to visualize the flux distribution across the surface at different times.
- This script uses `matplotlib`, a very popular python based visualization library.
-
-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 adapting
- 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.
-
-
-### 3. Use a turbulence model in the free flow domain
-
-Several RANS turbulence models are implemented in DuMu<sup>x</sup>.
-This part of the exercise consists of the following steps:
-* replacing the Navier-Stokes model by the K-Omega SST 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 same applies for types such as `GridManager`,
+```c++
+ // try to create a grid (from the given grid file or the input file)
+ // for both sub-domains
+ using PorousMediumGridManager = Dumux::GridManager<GetPropType<PorousMediumTypeTag, Properties::Grid>>;
+ PorousMediumGridManager porousMediumGridManager;
+ porousMediumGridManager.init("PorousMedium"); // pass parameter group
-The file `freeflowsubproblem.hh` is your free flow problem file and `properties.hh` is the properties file within this exercise.
+ using FreeflowGridManager = Dumux::GridManager<GetPropType<FreeflowTypeTag, Properties::Grid>>;
+ FreeflowGridManager freeflowGridManager;
+ freeflowGridManager.init("Freeflow"); // pass parameter group
-For using the compositional zero equation turbulence model, the following header files need to be included
- in properties file:
-```
-#include <dumux/freeflow/compositional/sstncmodel.hh>
-```
-and in problem file:
+ // we compute on the leaf grid view
+ const auto& porousMediumGridView = porousMediumGridManager.grid().leafGridView();
+ const auto& freeflowGridView = freeflowGridManager.grid().leafGridView();
```
-#include <dumux/freeflow/turbulenceproperties.hh>
-#include <dumux/freeflow/rans/twoeq/sst/problem.hh>
-#include <dumux/freeflow/rans/boundarytypes.hh>
+`FVGridGeometry`,
+```c++
+ // create the finite volume grid geometry
+ using FreeflowFVGridGeometry = GetPropType<FreeflowTypeTag, Properties::GridGeometry>;
+ auto freeflowFvGridGeometry = std::make_shared<FreeflowFVGridGeometry>(freeflowGridView);
+ using PorousMediumFVGridGeometry = GetPropType<PorousMediumTypeTag, Properties::GridGeometry>;
+ auto porousMediumFvGridGeometry = std::make_shared<PorousMediumFVGridGeometry>(porousMediumGridView);
```
-
-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 `FreeflowModel` definition accordingly (multi-component non-isothermal K-Omega SST model, `SSTNCNI`) in the properties file,
-* Adapt the inheritance of the problem class in problem file (Use `RANSProblem<TypeTag>` rather than `NavierStokesStaggeredProblem<TypeTag>`).
-
-Here, the turbulent free flow is wall-bounded, meaning shear stress and turbulence in the flow develop primarily due to the walls.
-The porous medium at the bottom and also the channel wall (upper boundary) act here as such walls.
-Within the problem the location of boundaries representing walls need to be defined.
-To do this, add the following function to the upper and lower boundaries within the `boundaryTypes` function in `freeflowsubproblem.hh`:
-```cpp
- values.setWall();
+`Problem`, etc...
+```c++
+ // the problem (initial and boundary conditions)
+ using FreeflowProblem = GetPropType<FreeflowTypeTag, Properties::Problem>;
+ auto freeflowProblem = std::make_shared<FreeflowProblem>(freeflowFvGridGeometry, couplingManager);
+ using PorousMediumProblem = GetPropType<PorousMediumTypeTag, Properties::Problem>;
+ auto porousMediumProblem = std::make_shared<PorousMediumProblem>(porousMediumFvGridGeometry, couplingManager);
```
-With the locations of the wall designated, and the problem initialized, a series of constant spatial properties,
-in particular the distance to the nearest wall from each cell center, should be initialized.
-To do this, add the following to the `main.cc` file.
-```cpp
- freeflowProblem->updateStaticWallProperties();
-```
+Since we use a monolithic coupling scheme, there is only one `Assembler` and one `NewtonSolver`, which help to assemble and solve the full coupled problem.
-In addition, there is also a non-local solution-dependent aspect of the turbulence models which is to be updated at the end of each step.
-An example of this is a stencil extended velocity gradient, and other flow properties in relation to the wall.
-These dynamic interactions are to be initialized in the mainfile directly after the `updateStaticWallProperties()` call,
-as well as in the time loop after `// Update dynamic wall properties`:
-```cpp
-freeflowProblem->updateDynamicWallProperties(freeflowSol);
+The problem files very much look like the "regular", uncoupled problem files seen in previous exercises, with the exception that they hold a pointer to the `CouplingManager`. E.g.
+```c++
+ FreeFlowSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry,std::shared_ptr<CouplingManager> couplingManager)
+ : ParentType(fvGridGeometry, "Freeflow"),
+ eps_(1e-6),
+ couplingManager_(couplingManager)
+ {
+ //...
+ }
```
-
-In addition to designating the locations of walls,
-additional boundary conditions and initial conditions need to be set for the two new primary variables $k$ and $\omega$.
-In the `boundaryTypes` function, set both variables on all walls to be dirichlet, except for the right boundary, which should have outflow conditions.
-
-For the initial conditions, Reynolds number specific base conditions should be applied everywhere.
-In the problem constructor, uncomment the code calculating these terms,
-then apply the `turbulentKineticEnergy_`and `dissipation_` variables to their primary variables in all locations.
-Within the dirichlet function for cell faces (`dirichlet(element, scvf)`),
-we also need to specify that these variables should be fixed to 0 at the wall.
-
-In addition, dirichlet cell constraints for the dissipation or $\omega$ variable will be set for all wall adjacent cells.
-This is done in the `isDirichletCell` function, as well as the `dirichlet` function already, and requires no further changes.
-
-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 at the top with
+This allows them to evaluate the coupling conditions and to exchange information between the coupled models.
+The coupling conditions are realized technically in terms of boundary conditions. For instance,
+in `interface/freeflowsubproblem.hh`, `couplingNeumann` boundary conditions are set, which means that the free flow model evaluates the
+mass and momentum fluxes coming from the porous domain and uses these values as boundary conditions at the interface.
```c++
-values.setAllSymmetry();
+ BoundaryTypes boundaryTypes(const Element& element,
+ const SubControlVolumeFace& scvf) const
+ {
+ BoundaryTypes values;
+ const auto& globalPos = scvf.center();
+
+ //...
+ // 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
+ }
+
+ return values;
+ }
```
+Note that 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.
-In addition, you have to remove the condition `onUpperBoundary_(globalPos)` from the `initialAtPos(globalPos)` method and the `dirichlet(element, scvf)` 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.
+We will use a staggered grid (also calles Marker-and-Cell method - MAC)to discretize the free-flow domain and a cell-centered finite volume (CCFV) method for the porous medium domain.
+Keep in mind that the staggered grid implementation distinguishes between face variables (velocity components) and cell center variables (all other variables).
+For this reason one distinguishes between `CouplingManager::stokesCellCenterIdx` and `CouplingManager::stokesFaceIdx` indices (see `main.cc`), while for the porous medium all variables can be accessed with `CouplingManager::darcyIdx`.
-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
+ // the indices
+ constexpr auto freeflowCellCenterIdx = CouplingManager::stokesCellCenterIdx;
+ constexpr auto freeflowFaceIdx = CouplingManager::stokesFaceIdx;
+ constexpr auto porousMediumIdx = CouplingManager::darcyIdx;
```
-__Task D__:
+__Task__:
+Take a closer look at the above listed files before moving to the three exercises below.
-For the grid convergence study, run various simulations with the following grading parameters:
-```c++
-* [Freeflow.Grid] Grading1 = 1.2, [PorousMedium.Grid] Grading1 = -1.2
-* [Freeflow.Grid] Grading1 = 1.3, [PorousMedium.Grid] Grading1 = -1.3
-* [Freeflow.Grid] Grading1 = 1.4, [PorousMedium.Grid] Grading1 = -1.4
-* [Freeflow.Grid] Grading1 = 1.5, [PorousMedium.Grid] Grading1 = -1.5
-* [Freeflow.Grid] Grading1 = 1.6, [PorousMedium.Grid] Grading1 = -1.6
-```
+## Sub-Exercises
-By changing the parameter `Problem.Name` for each grading factor, you avoid losing 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.
+* [**Exercise 1:** Changing the interface between the free- and the porous medium domain](./interface/readme.md)
+* [**Exercise 2:** Changing the porous medium model](./models/readme.md)
+* [**Exercise 3:** Introducing a turbulence model in the free flow domain](./turbulence/readme.md)
\ No newline at end of file
diff --git a/exercises/exercise-coupling-ff-pm/interface/freeflowsubproblem.hh b/exercises/exercise-coupling-ff-pm/interface/freeflowsubproblem.hh
index 0f7d590f..fdbd4c18 100644
--- a/exercises/exercise-coupling-ff-pm/interface/freeflowsubproblem.hh
+++ b/exercises/exercise-coupling-ff-pm/interface/freeflowsubproblem.hh
@@ -124,7 +124,9 @@ public:
PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
{
PrimaryVariables values(0.0);
- values = initialAtPos(globalPos);
+ values[Indices::velocityYIdx] = -1e-6 * globalPos[0] * (this->gridGeometry().bBoxMax()[0] - globalPos[0]);
+ // TODO: dumux-course-task 1.A
+ // set fixed pressures on the left and right boundary
return values;
}
@@ -164,20 +166,6 @@ public:
NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
{ return NumEqVector(0.0); }
- /*!
- * \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->gridGeometry().bBoxMax()[0] - globalPos[0]);
- // TODO: dumux-course-task 1.A
- // Set fixed pressures on the left and right boundary
-
- return values;
- }
/*!
* \brief Returns the intrinsic permeability of required as input parameter for the Beavers-Joseph-Saffman boundary condition
diff --git a/exercises/exercise-coupling-ff-pm/interface/main.cc b/exercises/exercise-coupling-ff-pm/interface/main.cc
index b869b4b9..1de0405b 100644
--- a/exercises/exercise-coupling-ff-pm/interface/main.cc
+++ b/exercises/exercise-coupling-ff-pm/interface/main.cc
@@ -163,9 +163,6 @@ int main(int argc, char** argv)
auto freeflowSol = partial(sol, freeflowFaceIdx, freeflowCellCenterIdx);
- freeflowProblem->applyInitialSolution(freeflowSol);
- porousMediumProblem->applyInitialSolution(sol[porousMediumIdx]);
-
couplingManager->init(freeflowProblem, porousMediumProblem, sol);
// the grid variables
@@ -186,7 +183,6 @@ int main(int argc, char** argv)
//TODO: dumux-course-task 1.B
//****** uncomment the add analytical solution of v_x *****//
// freeflowVtkWriter.addField(freeflowProblem->getAnalyticalVelocityX(), "analyticalV_x");
- freeflowVtkWriter.write(0.0);
using PorousMediumSolutionVector = GetPropType<PorousMediumTypeTag, Properties::SolutionVector>;
VtkOutputModule<PorousMediumGridVariables, PorousMediumSolutionVector> porousMediumVtkWriter(*porousMediumGridVariables,
@@ -196,7 +192,6 @@ int main(int argc, char** argv)
using PorousMediumVelocityOutput = GetPropType<PorousMediumTypeTag, Properties::VelocityOutput>;
porousMediumVtkWriter.addVelocityOutput(std::make_shared<PorousMediumVelocityOutput>(*porousMediumGridVariables));
GetPropType<PorousMediumTypeTag, Properties::IOFields>::initOutputModule(porousMediumVtkWriter);
- porousMediumVtkWriter.write(0.0);
// the assembler for a stationary problem
using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
@@ -217,6 +212,8 @@ int main(int argc, char** argv)
using NewtonSolver = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
+ std::cout << "\nSolving the coupled problem..." << std::endl;
+
// solve the non-linear system
nonLinearSolver.solve(sol);
diff --git a/exercises/exercise-coupling-ff-pm/interface/porousmediumsubproblem.hh b/exercises/exercise-coupling-ff-pm/interface/porousmediumsubproblem.hh
index 68967766..fe40fd96 100644
--- a/exercises/exercise-coupling-ff-pm/interface/porousmediumsubproblem.hh
+++ b/exercises/exercise-coupling-ff-pm/interface/porousmediumsubproblem.hh
@@ -103,7 +103,6 @@ public:
{
// set p = 0 at the bottom
PrimaryVariables values(0.0);
- values = initial(element);
return values;
}
@@ -151,16 +150,6 @@ public:
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)
diff --git a/exercises/exercise-coupling-ff-pm/interface/properties.hh b/exercises/exercise-coupling-ff-pm/interface/properties.hh
index 3f5d388d..6c5a19d3 100644
--- a/exercises/exercise-coupling-ff-pm/interface/properties.hh
+++ b/exercises/exercise-coupling-ff-pm/interface/properties.hh
@@ -134,6 +134,9 @@ struct EnableGridFluxVariablesCache<TypeTag, TTag::FreeflowOneP> { static conste
template<class TypeTag>
struct EnableGridVolumeVariablesCache<TypeTag, TTag::FreeflowOneP> { static constexpr bool value = true; };
+template<class TypeTag>
+struct NormalizePressure<TypeTag, TTag::FreeflowOneP> { static constexpr bool value = false; };
+
} //end namespace Dumux::Properties
#endif
diff --git a/exercises/exercise-coupling-ff-pm/interface/readme.md b/exercises/exercise-coupling-ff-pm/interface/readme.md
new file mode 100644
index 00000000..2b1bf9b7
--- /dev/null
+++ b/exercises/exercise-coupling-ff-pm/interface/readme.md
@@ -0,0 +1,169 @@
+## 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:
+
+
+
+Note that we neglect the influence of gravity and only solve for the stationary problem in this exercise.
+
+You can also compile and run the test without changing anything yet and investigate the above mentioned velocity distribution with paraview.
+
+* 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 `interface/freeflowsubproblem.hh` and navigate to the part, 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.
+You can also delete the vertical velocity from `dirichletAtPos`, as Dirichlet-type boundary conditions should be set for the pressure.
+
+Recompile and rerun your code. Then check your results with the applied changes by visualizing them e.g. with paraview.
+
+### 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.setBeaversJoseph(Indices::momentumXBalanceIdx);
+```
+
+at the position where the coupling boundary conditions are set in `interface/freeflowsubproblem.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 the following line in `main.cc` in the subfolder `interface`.
+```cpp
+freeflowVtkWriter.addField(freeflowProblem->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 in the input file to see how that affects the resulting velocity profile.
+
+### Task C: Change the shape of interface
+
+Now we want to include a non-flat interface between the two domains.
+We use [`dune-subgrid`](https://doi.org/10.1007/s00607-009-0067-2) to construct two grids for the two domains from one common host grid. Thus, for the following tasks subgrid needs to be correctly installed.
+Our hostgrid will be a Dune-Yasp grid (`Dune::YaspGrid<2, Dune::TensorProductCoordinates<double, dim> >`)
+and the bounds and resolution of the grid will be set in the `params.input`file under the group `[Grid]`.
+This hostgrid, along with `elementSelector` functions defining some spatial cut of this domain, are passed to the grid manager to create each subdomain grid.
+
+To introduce this, open `main.cc` in the subfolder `interface` again and search for `TODO: dumux-course-task 1.C`.
+Comment out the first code block and uncomment the second.
+This will instantiate a host grid and define two helper lambda functions that are used to choose elements from the host grid for the respective sub grid.
+In the given case, the domain is split in two halves, separated by a sinusoidal interface.
+
+```cpp
+auto elementSelectorFreeflow = [&](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 elementSelectorPorousMedium = [&](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 line for including the grid manager in the properties file, i.e.
+```cpp
+#include <dumux/io/grid/gridmanager_sub.hh>
+```
+
+and do the changes in the respective lines for the `Grid` property.
+
+The problem should now compile. However, a runtime 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());
+ ```
+in `freeflowsubproblem.hh` in the subfolder `interface`.
+
+The same is true for the BJS condition, however, here we need to consider the tangential direction:
+```cpp
+values.setBeaversJoseph(1 - scvf.directionIndex());
+```
+
+Recompile and rerun to check if the final result looks something like this:
+
+
+
+*Extra Points:*
+Rather than enforcing a pressure difference across the domain, an inflow velocity profile could be set.
+What changes to the left boundary conditions in the free-flow domain would you make to introduce this? What conditions can be enforced on the right boundary?
+Hint: A relation between velocity and position is used for the vertical velocity component in the original form of the `dirichletAtPos` method.
\ No newline at end of file
diff --git a/exercises/exercise-coupling-ff-pm/models/CMakeLists.txt b/exercises/exercise-coupling-ff-pm/models/CMakeLists.txt
index b44a8c7a..a11f26f5 100644
--- a/exercises/exercise-coupling-ff-pm/models/CMakeLists.txt
+++ b/exercises/exercise-coupling-ff-pm/models/CMakeLists.txt
@@ -1,4 +1,4 @@
-# executables for ex_interface_coupling_ff-pm
+# executables for exercise_models_coupling_ff-pm
dumux_add_test(NAME exercise_models_coupling_ff-pm
SOURCES main.cc
LABELS ffpm)
diff --git a/exercises/exercise-coupling-ff-pm/models/params.input b/exercises/exercise-coupling-ff-pm/models/params.input
index 1fa6c1c2..d036bf85 100644
--- a/exercises/exercise-coupling-ff-pm/models/params.input
+++ b/exercises/exercise-coupling-ff-pm/models/params.input
@@ -36,15 +36,15 @@ VanGenuchtenN = 8.0
VanGenuchtenAlpha = 6.5e-4
Swr = 0.005
Snr = 0.01
-PcLowSw = Swr * 5.0
-PcHighSw = 1.0 - Snr * 5.0
+VanGenuchtenPcLowSweThreshold = 0.01
+VanGenuchtenPcHighSweThreshold = 0.99
Temperature = 293.15
[Problem]
Name = models_coupling
-ExportStorage = true
-PlotStorage = true
-ExportFluxes = true
+ExportStorage = false
+PlotStorage = false
+ExportFluxes = false
[Newton]
MaxSteps = 12
diff --git a/exercises/exercise-coupling-ff-pm/models/porousmediumsubproblem.hh b/exercises/exercise-coupling-ff-pm/models/porousmediumsubproblem.hh
index 2f06ee94..1250de48 100644
--- a/exercises/exercise-coupling-ff-pm/models/porousmediumsubproblem.hh
+++ b/exercises/exercise-coupling-ff-pm/models/porousmediumsubproblem.hh
@@ -348,6 +348,8 @@ public:
// TODO: dumux-course-task 2.A
// Declare here which phases are present.
+ // TODO: dumux-course-task 2.C
+ // Change initial condition to 2p system with liquid saturation of 0.1
values[Indices::pressureIdx] = freeflowPressure;
values[transportCompIdx] = moleFraction_;
diff --git a/exercises/exercise-coupling-ff-pm/models/readme.md b/exercises/exercise-coupling-ff-pm/models/readme.md
new file mode 100644
index 00000000..6e123170
--- /dev/null
+++ b/exercises/exercise-coupling-ff-pm/models/readme.md
@@ -0,0 +1,118 @@
+## 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 domains is included and the presence of a second phase in the porous medium domain is considered.
+This enables the simulation of the drying of a porous medium (however no energy balance is included yet).
+
+Compared to exercise `1. Changing the interface` account now for the gravity effect as well as solving for the transient problem.
+
+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 used, 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 (`models/properties.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.
+* Since two phases are involved now, we do not need to use the `OnePAdapter` anymore. Change the property of the `FluidSystem` such that `H2OAir` is used directly.
+ Afterwards, set the `transportCompIdx` to `Indices::switchIdx` in `porousmediumsubproblem.hh`.
+
+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::PorousMediumOnePNC>
+{ static constexpr auto value = TwoPFormulation::p1s0; };
+```
+ in the properties file.
+
+In contrast to the formulation, which remains the same during one simulation,
+the primary variables may vary during one simulation.
+In the case under investigation, 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 the 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()`. `MY_PHASE_PRESENCE` should be replaced with the correct value for the case where only a gas-phase (second phase of the fluid system) is present. For finding out if a gas or liquid phase is the first of second phase of a fluidsystem (here `h20air`), take a look at your fluidsystem (e.g. [h2oair](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/blob/master/dumux/material/fluidsystems/h2oair.hh?ref_type=heads#L75) ).
+ Moreover, have a look at the [indices.hh](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/blob/master/dumux/porousmediumflow/2pnc/indices.hh?ref_type=heads#L22) in the `2pnc` model (as the 2p2c model is a special case of the 2pnc model) in the subfolder `porousmediumflow` in your DuMuX directory to see which value to set for the gas-phase (here: second phase of our fluidsystem).
+
+Make sure your code is compiling and running. You can have a look at the results with paraview. How do the mass or mole fractions of water in air change over time?
+
+### Task B: Add output
+
+In the next step, we want to add some output to the simulation. The standard method for providing simulation output for visualization is via a `VtkWriter`. These tools take the grid geometries of each domain, as well as the solutions and write spatial output that one can view in visualization tools such as paraview.
+
+Although this Vtk output is very useful, some output is more suited for other forms of visualization.
+Two examples of data output formats are `.csv` files and `.json` files.
+From these files, data can be visualized using many different tools.
+In this exercise, data exported to a `.csv` file will be plotted using gnuplot, and data exported to a `.json` output will be plotted using the matplotlib python library.
+
+First as example output data, we want to investigate the water mass in the (porous-medium) system.
+There are two ways to evaluate this: (1) the total storage of water mass in the porous medium domain, and (2) the water vapor flux along the interface.
+The total storage can be plotted over time, whereas the water vapor flux will vary spatially along the interface as well as over time.
+
+First, we evaluate the storage term of the water component.
+* Have a look at the function `evaluateWaterMassStorageTerm()` in the `porousmediumsubproblem.hh`.
+ Then implement a calculation of the total water mass:
+ $`M^\text{w}:=\sum_{\alpha \in \textrm{g,l}} \left( \phi S_\alpha X^\text{w}_\alpha \varrho_\alpha V_\textrm{scv} \right)`$.
+ Hint: the volume is a property of a subcontrolvolume (`scv`), the other properties can be obtained through the volume variables (`volVars`). The functions for the `volVars` can be found in the header [dumux/porousmediumflow/2pnc/volumevariables.hh](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/blob/master/dumux/porousmediumflow/2pnc/volumevariables.hh?ref_type=heads)
+* Afterwards, adapt the method `startStorageEvaluation()` such that the variable `initialWaterContent_` is
+ initialized correctly using the `evaluateWaterMassStorageTerm()` method and assign that value also to the variable `lastWaterMass_`.
+
+After these functions are correctly implemented, the evaporation in [mm/d] is calculated in the code as
+$`e = \frac{\textrm{d}\, M^\text{w}}{\textrm{d}\, t} / A_\textrm{interface}`$.
+(all boundaries have Neumann no-flow conditions and no sources are present, meaning the water can only leave the system
+via the porous-medium free-flow interface.)
+
+During the simulation a `.csv` file is created (enable by setting `[Problem.ExportStorage] = true` in the `params.input` file). In addition, a gnuplot script `StorageOverTime.gp` is written and the mass loss and cumulative mass are plotted over time (`StorageOverTime.png`), when plotting is enabled with `[Problem.PlotStorage] = true`.
+
+Second, 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. Multiply this with the relevant face areas, extrusion factor, molar Mass, and seconds per day to get the [mm/d] units as seen in the storage evaluation.
+* When the `[Problem.ExportFluxes] = true` parameter is enabled, simulation data will be exported for this simulation to a `.json` file.
+ This file can flexibly handle data of different types, which in this case is helpful as we have both temporal and spatial data.
+* Use the python plotting script `plotFluxes.py` to visualize the flux distribution across the surface at different times.
+ This script uses `matplotlib`, a very popular python based visualization library and reads-in the data stored in the file `flux_models_coupling.json`.
+ You can execute the script by
+ ```
+ python3 ./plotFluxes.py
+ ```
+
+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). The regularization of the capillary pressure-saturation relationship with using the vanGenuchten model can be found in the file [dumux/material/fluidmatrixinteractions/2p/vangenuchten.hh](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/-/blob/master/dumux/material/fluidmatrixinteractions/2p/vangenuchten.hh?ref_type=heads#L304).
+* See how adapting the parameter for `VanGenuchtenPcLowSweThreshold` and `VanGenuchtenPcHighSweThreshold` 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. Have a look the resulting liquid saturation distribution within the porous medium (with using paraview).
\ No newline at end of file
diff --git a/exercises/exercise-coupling-ff-pm/turbulence/params.input b/exercises/exercise-coupling-ff-pm/turbulence/params.input
index b64b0558..c9d26a57 100644
--- a/exercises/exercise-coupling-ff-pm/turbulence/params.input
+++ b/exercises/exercise-coupling-ff-pm/turbulence/params.input
@@ -5,10 +5,13 @@ TEnd = 864000 # [s] (6 days)
[Freeflow.Grid]
Positions0 = 0.0 0.25
+# TODO: dumux-course-task 3.B - use only half of FF domain height
Positions1 = 0.25 0.5
+# TODO: dumux-course-task 3.C - refine towards interface
Grading0 = 1.0
Grading1 = 1.0
Cells0 = 15
+# TODO: dumux-course-task 3.B - use only half of FF domain height and adapt cell number
Cells1 = 20
Verbosity = true
@@ -17,6 +20,7 @@ Positions0 = 0.0 0.25
Positions1 = 0.0 0.25
Cells0 = 15
Cells1 = 10
+# TODO: dumux-course-task 3.C - refine towards interface
Grading0 = 1.0
Grading1 = 1.0
Verbosity = true
@@ -44,8 +48,6 @@ VanGenuchtenN = 6.9
VanGenuchtenAlpha = 6.371e-4
Swr = 0.005
Snr = 0.01
-PcLowSw = Swr * 5.0
-PcHighSw = 1.0 - Snr * 5.0
Temperature = 298.15 # [K]
[Problem]
diff --git a/exercises/exercise-coupling-ff-pm/turbulence/readme.md b/exercises/exercise-coupling-ff-pm/turbulence/readme.md
new file mode 100644
index 00000000..23b2106d
--- /dev/null
+++ b/exercises/exercise-coupling-ff-pm/turbulence/readme.md
@@ -0,0 +1,126 @@
+## 3. Use a turbulence model in the free flow domain
+As in the previous exercise, we account for the gravity effect as well as solving for the transient problem hereafter.
+
+Several RANS (Reynolds Averaged Navier Stokes) turbulence models are implemented in DuMu<sup>x</sup>.
+This part of the exercise consists of the following steps:
+* replacing the Navier-Stokes model by the [K-Omega SST (shear stress transport) turbulence model](https://www.cfd-online.com/Wiki/SST_k-omega_model) of [Menter (1994)](https://doi.org/10.2514/3.12149), which is a two-equation RANS 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: Use a $k$-$\omega$ turbulence model
+
+For using the compositional $k$-$\omega$ turbulence model, the following header files need to be included in properties file (`properties.hh`):
+```
+#include <dumux/freeflow/compositional/sstncmodel.hh>
+```
+and in free-flow problem file (`freeflowsubproblem.hh`):
+```
+#include <dumux/freeflow/turbulenceproperties.hh>
+#include <dumux/freeflow/rans/twoeq/sst/problem.hh>
+#include <dumux/freeflow/rans/boundarytypes.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 `FreeflowModel` definition accordingly (multi-component non-isothermal K-Omega SST model, `SSTNCNI`) in the properties file,
+* Adapt the inheritance of the problem class in problem file (Use `RANSProblem<TypeTag>` rather than `NavierStokesStaggeredProblem<TypeTag>`).
+
+Additionally make sure that the alias `BoundaryTypes` uses the new boundary types for the RANS model (as included by the respective header above) instead of `NavierStokesBoundaryTypes`.
+
+Here, the turbulent free flow is wall-bounded, meaning shear stress and turbulence in the flow develop primarily due to the walls.
+The porous medium at the bottom and also the channel wall (upper boundary) act here as such walls.
+Within the problem the location of boundaries representing walls need to be defined.
+To do this, add the following function to the upper and lower boundaries within the `boundaryTypes` function in `freeflowsubproblem.hh`:
+```cpp
+ values.setWall();
+```
+
+With the locations of the wall designated, and the problem initialized, a series of constant spatial properties,
+in particular the distance to the nearest wall from each cell center, should be initialized.
+To do this, add the following to the `main.cc` file.
+```cpp
+ freeflowProblem->updateStaticWallProperties();
+```
+
+In addition, there is also a non-local solution-dependent aspect of the turbulence models which is to be updated at the end of each step.
+An example of this is a stencil extended velocity gradient, and other flow properties in relation to the wall.
+These dynamic interactions are to be initialized in the mainfile directly after the `updateStaticWallProperties()` call,
+as well as in the time loop after `// Update dynamic wall properties`:
+```cpp
+freeflowProblem->updateDynamicWallProperties(freeflowSol);
+```
+
+In addition to designating the locations of walls,
+additional boundary conditions and initial conditions need to be set for the two new primary variables, the turbulent kinetic energy $k$ and the specific tubulent dissipation rate/turbulence frequency $\omega$.
+
+In the `boundaryTypes` function, set both variables ($k$ and $\omega$) on all free-flow domain boundaries to be dirichlet, except for the right boundary, which should have outflow conditions. The name of the indices for those two variables that have to be used for setting the boundary conditions types can be found in `dumux/freeflow/rans/twoeq/indices.hh`.
+
+Within the dirichlet function for cell faces (`dirichlet(element, scvf)`),
+we also need to specify that these variables ($k$ and $\omega$) should be fixed to 0 at the wall.
+In addition, dirichlet cell constraints for the dissipation or $\omega$ variable will be set for all wall adjacent cells.
+This is done in the `isDirichletCell` function, as well as the `dirichlet` function already, and requires no further changes.
+
+For the initial conditions, Reynolds number specific base conditions should be applied everywhere.
+In the problem constructor, uncomment the code calculating these terms,
+then apply the `turbulentKineticEnergy_`and `dissipation_` variables to their primary variables in the initial state for all spatial locations.
+
+Before starting to run your modified code make sure to include the parameter `Rans.IsFlatWallBounded` in your input file and set this parameter to `true`.
+
+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 diffusivity should look like this:
+
+
+
+### Task B: Use symmetry boundary conditions
+
+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 at the top with
+```c++
+values.setAllSymmetry();
+```
+
+In addition, you have to remove the condition `onUpperBoundary_(globalPos)` from the `initialAtPos(globalPos)` method and the `dirichlet(element, scvf)` method.
+
+Now it is sufficient to only calulate half of the domain height for the free-flow domain. For this, adapt the grid coordinates as well as the number of cells in y-direction in your `params.input` - file accordingly.
+After recompiling and running your simulation with the new boundary condition, check your results with paraview.
+Does your velocity profile at the symmetry boundary look reasonable?
+
+Before closing paraview, Choose `Surface With Edges` instead of `Surface` in the Paraview toolbar to see the discretization grid. Until now we used a regular structured lattice grid.
+
+### Task C: Refine grid towards the interface
+
+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
+```
+
+Rerun your simulation with the new grid and check your refined grid with paraview.
+
+### Task D: Perform a grid convergence study
+
+For the grid convergence study, run various simulations with the following grading parameters:
+```c++
+* [Freeflow.Grid] Grading1 = 1.2, [PorousMedium.Grid] Grading1 = -1.2
+* [Freeflow.Grid] Grading1 = 1.3, [PorousMedium.Grid] Grading1 = -1.3
+* [Freeflow.Grid] Grading1 = 1.4, [PorousMedium.Grid] Grading1 = -1.4
+* [Freeflow.Grid] Grading1 = 1.5, [PorousMedium.Grid] Grading1 = -1.5
+* [Freeflow.Grid] Grading1 = 1.6, [PorousMedium.Grid] Grading1 = -1.6
+```
+
+By changing the parameter `Problem.Name` for each grading factor, you avoid losing the `.vtu` and `.pvd` files of the previous simulation runs.
+Check the first lines of the output in your terminal to see how the grading factors change the height of your grid cells.
+
+Compare the velocity fields with Paraview. In Paraview you can e.g. use the tool `PlotOverLine` to plot the velocity in $x$-direction ($v_x$) over the oulet on the right hand side. After choosing `PlotOverLine` make sure to give the Point 1 as (0.25, 0.25, 0) and the Point 2 as (0.25, 0.5, 0) to define the line a the output. Apply this and choose to only visualize `velocity_gas(m/s)_X`. You can do this for all your refinements and plot them in the same diagram to compare.
diff --git a/exercises/solution/exercise-coupling-ff-pm/interface/freeflowsubproblem.hh b/exercises/solution/exercise-coupling-ff-pm/interface/freeflowsubproblem.hh
index 25606fdb..8f5f4e07 100644
--- a/exercises/solution/exercise-coupling-ff-pm/interface/freeflowsubproblem.hh
+++ b/exercises/solution/exercise-coupling-ff-pm/interface/freeflowsubproblem.hh
@@ -148,7 +148,18 @@ public:
PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
{
PrimaryVariables values(0.0);
- values = initialAtPos(globalPos);
+#if EXNUMBER == 0
+ values[Indices::velocityYIdx] = -1e-6 * globalPos[0] * (this->gridGeometry().bBoxMax()[0] - globalPos[0]);
+#elif EXNUMBER == 4
+ values[Indices::velocityXIdx] = 1e-6 * (globalPos[1] - this->gridGeometry().bBoxMin()[1])
+ * (this->gridGeometry().bBoxMax()[1] - globalPos[1]);
+#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;
}
@@ -192,28 +203,6 @@ public:
NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
{ return NumEqVector(0.0); }
- /*!
- * \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->gridGeometry().bBoxMax()[0] - globalPos[0]);
-#elif EXNUMBER == 4
- values[Indices::velocityXIdx] = 1e-6 * (globalPos[1] - this->gridGeometry().bBoxMin()[1])
- * (this->gridGeometry().bBoxMax()[1] - globalPos[1]);
-#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
diff --git a/exercises/solution/exercise-coupling-ff-pm/interface/main.cc b/exercises/solution/exercise-coupling-ff-pm/interface/main.cc
index 94536d74..b3ee5f0a 100644
--- a/exercises/solution/exercise-coupling-ff-pm/interface/main.cc
+++ b/exercises/solution/exercise-coupling-ff-pm/interface/main.cc
@@ -155,9 +155,6 @@ int main(int argc, char** argv)
auto freeflowSol = partial(sol, freeflowFaceIdx, freeflowCellCenterIdx);
- freeflowProblem->applyInitialSolution(freeflowSol);
- porousMediumProblem->applyInitialSolution(sol[porousMediumIdx]);
-
couplingManager->init(freeflowProblem, porousMediumProblem, sol);
// the grid variables
@@ -179,8 +176,6 @@ int main(int argc, char** argv)
freeflowVtkWriter.addField(freeflowProblem->getAnalyticalVelocityX(), "analyticalV_x");
#endif
- freeflowVtkWriter.write(0.0);
-
using PorousMediumSolutionVector = GetPropType<PorousMediumTypeTag, Properties::SolutionVector>;
VtkOutputModule<PorousMediumGridVariables, PorousMediumSolutionVector> porousMediumVtkWriter(*porousMediumGridVariables,
sol[porousMediumIdx],
@@ -189,7 +184,6 @@ int main(int argc, char** argv)
using PorousMediumVelocityOutput = GetPropType<PorousMediumTypeTag, Properties::VelocityOutput>;
porousMediumVtkWriter.addVelocityOutput(std::make_shared<PorousMediumVelocityOutput>(*porousMediumGridVariables));
GetPropType<PorousMediumTypeTag, Properties::IOFields>::initOutputModule(porousMediumVtkWriter);
- porousMediumVtkWriter.write(0.0);
// the assembler for a stationary problem
using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
@@ -210,6 +204,8 @@ int main(int argc, char** argv)
using NewtonSolver = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
+ std::cout << "\nSolving the coupled problem..." << std::endl;
+
// solve the non-linear system
nonLinearSolver.solve(sol);
diff --git a/exercises/solution/exercise-coupling-ff-pm/interface/porousmediumsubproblem.hh b/exercises/solution/exercise-coupling-ff-pm/interface/porousmediumsubproblem.hh
index e7aa0ba3..deb4a9e7 100644
--- a/exercises/solution/exercise-coupling-ff-pm/interface/porousmediumsubproblem.hh
+++ b/exercises/solution/exercise-coupling-ff-pm/interface/porousmediumsubproblem.hh
@@ -103,7 +103,6 @@ public:
{
// set p = 0 at the bottom
PrimaryVariables values(0.0);
- values = initial(element);
return values;
}
@@ -151,16 +150,6 @@ public:
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)
diff --git a/exercises/solution/exercise-coupling-ff-pm/interface/properties.hh b/exercises/solution/exercise-coupling-ff-pm/interface/properties.hh
index 36da2f0d..ce959b8d 100644
--- a/exercises/solution/exercise-coupling-ff-pm/interface/properties.hh
+++ b/exercises/solution/exercise-coupling-ff-pm/interface/properties.hh
@@ -132,6 +132,9 @@ struct EnableGridFluxVariablesCache<TypeTag, TTag::FreeflowOneP> { static conste
template<class TypeTag>
struct EnableGridVolumeVariablesCache<TypeTag, TTag::FreeflowOneP> { static constexpr bool value = true; };
+template<class TypeTag>
+struct NormalizePressure<TypeTag, TTag::FreeflowOneP> { static constexpr bool value = false; };
+
} //end namespace Dumux::Properties
#endif
diff --git a/exercises/solution/exercise-coupling-ff-pm/models/CMakeLists.txt b/exercises/solution/exercise-coupling-ff-pm/models/CMakeLists.txt
index a7f4ba40..1a529025 100644
--- a/exercises/solution/exercise-coupling-ff-pm/models/CMakeLists.txt
+++ b/exercises/solution/exercise-coupling-ff-pm/models/CMakeLists.txt
@@ -3,27 +3,27 @@ dumux_add_test(NAME exercise_models_coupling_ff-pm_original
LABELS ffpm
COMPILE_DEFINITIONS EXNUMBER=0
COMMAND ${CMAKE_CURRENT_BINARY_DIR}/exercise_models_coupling_ff-pm_original
- CMD_ARGS -Problem.ExportStorage false -Problem.PlotStorage 0 -Problem.ExportFluxes 0)
+ CMD_ARGS params_orig_a.input)
dumux_add_test(NAME exercise_models_coupling_ff-pm_a_solution
SOURCES main.cc
LABELS ffpm
COMPILE_DEFINITIONS EXNUMBER=1
COMMAND ${CMAKE_CURRENT_BINARY_DIR}/exercise_models_coupling_ff-pm_a_solution
- CMD_ARGS -Problem.ExportStorage false -Problem.PlotStorage 0 -Problem.ExportFluxes 0)
+ CMD_ARGS params_orig_a.input)
dumux_add_test(NAME exercise_models_coupling_ff-pm_b_solution
SOURCES main.cc
LABELS ffpm
COMPILE_DEFINITIONS EXNUMBER=2
COMMAND ${CMAKE_CURRENT_BINARY_DIR}/exercise_models_coupling_ff-pm_b_solution
- CMD_ARGS -Problem.ExportStorage false -Problem.PlotStorage 0 -Problem.ExportFluxes 0)
+ CMD_ARGS params_b_c.input)
dumux_add_test(NAME exercise_models_coupling_ff-pm_c_solution
SOURCES main.cc
LABELS ffpm
COMPILE_DEFINITIONS EXNUMBER=3
COMMAND ${CMAKE_CURRENT_BINARY_DIR}/exercise_models_coupling_ff-pm_c_solution
- CMD_ARGS -Problem.ExportStorage false -Problem.PlotStorage 0 -Problem.ExportFluxes 0)
+ CMD_ARGS params_b_c.input)
-dune_symlink_to_source_files(FILES "params.input" plotFluxes.py)
+dune_symlink_to_source_files(FILES "params_orig_a.input" "params_b_c.input" plotFluxes.py)
diff --git a/exercises/solution/exercise-coupling-ff-pm/models/params.input b/exercises/solution/exercise-coupling-ff-pm/models/params_b_c.input
similarity index 92%
rename from exercises/solution/exercise-coupling-ff-pm/models/params.input
rename to exercises/solution/exercise-coupling-ff-pm/models/params_b_c.input
index 54e20321..1eba77e4 100644
--- a/exercises/solution/exercise-coupling-ff-pm/models/params.input
+++ b/exercises/solution/exercise-coupling-ff-pm/models/params_b_c.input
@@ -23,9 +23,8 @@ Velocity = 1 # m/s
[PorousMedium.Problem]
Name = porousmedium
EnableGravity = true
-Saturation = 0.1 # -
MoleFraction = 0.1 # -
-Pressure = 1e5 # Pa
+Saturation = 0.1 # -
[SpatialParams]
Permeability = 2.65e-10 # m^2
@@ -36,8 +35,8 @@ VanGenuchtenN = 8.0
VanGenuchtenAlpha = 6.5e-4
Swr = 0.005
Snr = 0.01
-PcLowSw = Swr * 5.0
-PcHighSw = 1.0 - Snr * 5.0
+VanGenuchtenPcLowSweThreshold = 0.01
+VanGenuchtenPcHighSweThreshold = 0.99
Temperature = 293.15
[Problem]
diff --git a/exercises/solution/exercise-coupling-ff-pm/models/params_orig_a.input b/exercises/solution/exercise-coupling-ff-pm/models/params_orig_a.input
new file mode 100644
index 00000000..b8984ede
--- /dev/null
+++ b/exercises/solution/exercise-coupling-ff-pm/models/params_orig_a.input
@@ -0,0 +1,56 @@
+[TimeLoop]
+DtInitial = 100 # s
+EpisodeLength = -360 # s # 0.25 days
+TEnd = 256000 # s # 2 days
+
+[Freeflow.Grid]
+LowerLeft = 0 1
+UpperRight = 1 2
+Cells = 16 16
+
+[PorousMedium.Grid]
+UpperRight = 1 1
+Cells = 16 16
+
+[Freeflow.Problem]
+Name = freeflow
+EnableGravity = true
+EnableInertiaTerms = true
+MoleFraction = 0.0 # -
+Pressure = 1e5 # Pa
+Velocity = 1 # m/s
+
+[PorousMedium.Problem]
+Name = porousmedium
+EnableGravity = true
+MoleFraction = 0.1 # -
+
+[SpatialParams]
+Permeability = 2.65e-10 # m^2
+Porosity = 0.4 # -
+AlphaBeaversJoseph = 1.0 # -
+# EXNUMBER >= 1
+VanGenuchtenN = 8.0
+VanGenuchtenAlpha = 6.5e-4
+Swr = 0.005
+Snr = 0.01
+VanGenuchtenPcLowSweThreshold = 0.01
+VanGenuchtenPcHighSweThreshold = 0.99
+Temperature = 293.15
+
+[Problem]
+Name = models_coupling
+ExportStorage = false
+PlotStorage = false
+ExportFluxes = false
+
+[Newton]
+MaxSteps = 12
+MaxRelativeShift = 1e-7
+TargetSteps = 5
+
+[Vtk]
+AddVelocity = 1
+
+[Assembly]
+NumericDifference.BaseEpsilon = 1e-8
diff --git a/exercises/solution/exercise-coupling-ff-pm/turbulence/CMakeLists.txt b/exercises/solution/exercise-coupling-ff-pm/turbulence/CMakeLists.txt
index 2b9a7257..ca29c21e 100644
--- a/exercises/solution/exercise-coupling-ff-pm/turbulence/CMakeLists.txt
+++ b/exercises/solution/exercise-coupling-ff-pm/turbulence/CMakeLists.txt
@@ -1,17 +1,30 @@
dumux_add_test(NAME exercise_turbulence_coupling_ff-pm_original
SOURCES main.cc
LABELS ffpm
- COMPILE_DEFINITIONS EXNUMBER=0)
+ COMPILE_DEFINITIONS EXNUMBER=0
+ COMMAND ./exercise_turbulence_coupling_ff-pm_original
+ CMD_ARGS params_orig_a.input)
dumux_add_test(NAME exercise_turbulence_coupling_ff-pm_a_solution
SOURCES main.cc
LABELS ffpm
- COMPILE_DEFINITIONS EXNUMBER=1)
+ COMPILE_DEFINITIONS EXNUMBER=1
+ COMMAND ./exercise_turbulence_coupling_ff-pm_a_solution
+ CMD_ARGS params_orig_a.input)
dumux_add_test(NAME exercise_turbulence_coupling_ff-pm_b_solution
SOURCES main.cc
LABELS ffpm
- COMPILE_DEFINITIONS EXNUMBER=2)
+ COMPILE_DEFINITIONS EXNUMBER=2
+ COMMAND ./exercise_turbulence_coupling_ff-pm_b_solution
+ CMD_ARGS params_b.input)
+
+dumux_add_test(NAME exercise_turbulence_coupling_ff-pm_c_d_solution
+ SOURCES main.cc
+ LABELS ffpm
+ COMPILE_DEFINITIONS EXNUMBER=3
+ COMMAND ./exercise_turbulence_coupling_ff-pm_c_d_solution
+ CMD_ARGS params_c_d.input)
# add a symlink for each input file
add_input_file_links()
diff --git a/exercises/solution/exercise-coupling-ff-pm/turbulence/params_b.input b/exercises/solution/exercise-coupling-ff-pm/turbulence/params_b.input
new file mode 100644
index 00000000..09f9cf8f
--- /dev/null
+++ b/exercises/solution/exercise-coupling-ff-pm/turbulence/params_b.input
@@ -0,0 +1,72 @@
+[TimeLoop]
+DtInitial = 1e-1 # [s]
+MaxTimeStepSize = 43200 # [s] (12 hours)
+TEnd = 864000 # [s] (6 days)
+
+[Freeflow.Grid]
+Positions0 = 0.0 0.25
+Positions1 = 0.25 0.375
+Grading0 = 1.0
+Grading1 = 1.0
+Cells0 = 15
+Cells1 = 10
+Verbosity = true
+
+[PorousMedium.Grid]
+Positions0 = 0.0 0.25
+Positions1 = 0.0 0.25
+Cells0 = 15
+Cells1 = 10
+Grading0 = 1.0
+Grading1 = 1.0
+Verbosity = true
+
+[Freeflow.Problem]
+Name = freeflow
+RefVelocity = 3.5 # [m/s]
+RefPressure = 1e5 # [Pa]
+refMoleFrac = 0 # [-]
+RefTemperature = 298.15 # [K]
+EnableInertiaTerms = true
+
+[PorousMedium.Problem]
+Name = porousmedium
+Pressure = 1e5
+Saturation = 0.5 # initial Sw
+Temperature = 298.15 # [K]
+InitPhasePresence = 3 # bothPhases
+
+[SpatialParams]
+Porosity = 0.41
+Permeability = 2.65e-10
+AlphaBeaversJoseph = 1.0
+VanGenuchtenN = 6.9
+VanGenuchtenAlpha = 6.371e-4
+Swr = 0.005
+Snr = 0.01
+Temperature = 298.15 # [K]
+
+[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
+
+[RANS]
+IsFlatWallBounded = True
diff --git a/exercises/solution/exercise-coupling-ff-pm/turbulence/params_c_d.input b/exercises/solution/exercise-coupling-ff-pm/turbulence/params_c_d.input
new file mode 100644
index 00000000..eb1f47cf
--- /dev/null
+++ b/exercises/solution/exercise-coupling-ff-pm/turbulence/params_c_d.input
@@ -0,0 +1,72 @@
+[TimeLoop]
+DtInitial = 1e-1 # [s]
+MaxTimeStepSize = 43200 # [s] (12 hours)
+TEnd = 864000 # [s] (6 days)
+
+[Freeflow.Grid]
+Positions0 = 0.0 0.25
+Positions1 = 0.25 0.375
+Grading0 = 1.0
+Grading1 = 1.2
+Cells0 = 15
+Cells1 = 10
+Verbosity = true
+
+[PorousMedium.Grid]
+Positions0 = 0.0 0.25
+Positions1 = 0.0 0.25
+Cells0 = 15
+Cells1 = 10
+Grading0 = 1.0
+Grading1 = -1.2
+Verbosity = true
+
+[Freeflow.Problem]
+Name = freeflow
+RefVelocity = 3.5 # [m/s]
+RefPressure = 1e5 # [Pa]
+refMoleFrac = 0 # [-]
+RefTemperature = 298.15 # [K]
+EnableInertiaTerms = true
+
+[PorousMedium.Problem]
+Name = porousmedium
+Pressure = 1e5
+Saturation = 0.5 # initial Sw
+Temperature = 298.15 # [K]
+InitPhasePresence = 3 # bothPhases
+
+[SpatialParams]
+Porosity = 0.41
+Permeability = 2.65e-10
+AlphaBeaversJoseph = 1.0
+VanGenuchtenN = 6.9
+VanGenuchtenAlpha = 6.371e-4
+Swr = 0.005
+Snr = 0.01
+Temperature = 298.15 # [K]
+
+[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
+
+[RANS]
+IsFlatWallBounded = True
diff --git a/exercises/solution/exercise-coupling-ff-pm/turbulence/params.input b/exercises/solution/exercise-coupling-ff-pm/turbulence/params_orig_a.input
similarity index 96%
rename from exercises/solution/exercise-coupling-ff-pm/turbulence/params.input
rename to exercises/solution/exercise-coupling-ff-pm/turbulence/params_orig_a.input
index a4d7be6e..b7947b88 100644
--- a/exercises/solution/exercise-coupling-ff-pm/turbulence/params.input
+++ b/exercises/solution/exercise-coupling-ff-pm/turbulence/params_orig_a.input
@@ -44,8 +44,6 @@ VanGenuchtenN = 6.9
VanGenuchtenAlpha = 6.371e-4
Swr = 0.005
Snr = 0.01
-PcLowSw = Swr * 5.0
-PcHighSw = 1.0 - Snr * 5.0
Temperature = 298.15 # [K]
[Problem]
--
GitLab