Commit 6695f293 authored by Sina Ackermann's avatar Sina Ackermann Committed by Kilian Weishaupt
Browse files

[coupling-ff-pm] Add turbulunce exercise

parent 5df03184
......@@ -206,7 +206,7 @@ In the first task, the porous-medium model will be changed from a 1p2c system to
Although a 2p2c system is plugged in, we still want to simulate the same situation as before.
The following changes have to be made in the porous-medium model (`ex_models_pmproblem.hh`):
* Include the 2pnc model: include the respective headers and inherit from the new model `TwoPNC`
* Exchange the spatial parameters for the 1-phase system by those for a 2-phase system (hint: two occurrences).
* Exchange the spatial parameters for the 1-phase system by those for a 2-phase system (hint: two occurrences).
One big difference between the 1p and 2p model is, that the primary variables for which
the problem is solved, are not fixed.
......@@ -285,7 +285,7 @@ The regularization has a strong influence on how long liquid water is present at
The porous-medium model can now reach a liquid saturation of zero. As already explained above,
for this case the liquid saturation cannot serve as primary variable anymore. However,
manually adapting the primary variable states and values is inconvenient.
[Class et al. (2007)](
[Class et al. (2002)](
describe an algorithm to switch the primary variables, if phases should appear or disappear during a simulation.
- Replace the current implementation of the Newton solver with the version which can handle
primary variable switches: `dumux/multidomain/privarswitchnewtonsolver.hh`.
......@@ -295,4 +295,95 @@ describe an algorithm to switch the primary variables, if phases should appear o
Now you are able to simulate a complete drying of the porous medium.
### 4. ...
### 4. Use a turbulence model in the free flow domain
Several RANS turbulence models are implemented in DuMuX.
This part of the exercise consists of the following steps:
* replacing the Navier-Stokes model by the zero equation turbulence model,
* switching to a symmetry boundary condition,
* applying a grid refinement towards the interface,
* subsequently refining the grid (convergence study).
We will work with a `1p2cni/2p2cni` coupled problem, where `ni` stands for non-isothermal, and take the inertial forces into account.
All the prepared files can be found in the subfolder `exercise-coupling-ff-pm/turbulence`.
__Task A__:
The file `ex_turbulence_ffproblem.hh` is your free flow problem file within this exercise.
For using the compositional zero equation turbulence model, the following header files need to be included:
#include <dumux/freeflow/compositional/zeroeqncmodel.hh>
#include <dumux/freeflow/rans/zeroeq/problem.hh>
The includes for the NavierStokesNC model and the NavierStokesProblem are no longer needed and can be removed.
Make sure your free flow problem inherits from the correct parent type:
* Change the last entry in the `NEW_TYPE_TAG` definition accordingly (non-isothermal zero equation model)
* Adapt the inheritance of the problem class (hint: two occurrences)
Take a look into the two headers included above to see how the correct TypeTag and the Problem class the inherit from are called.
Additionally, you have to update the static and dynamic wall properties with
after applying the initial solution to the Stokes problem in the main file.
In order to update these properties, the Stokes problem file needs to be extended to provide an isOnWall() method:
bool isOnWall(const GlobalPosition& globalPos) const
return (onLowerBoundary_(globalPos) || onUpperBoundary_(globalPos));
Since the dynamic wall properties can change during the simulation and depend on the solution, they have to be updated in each time step.
Include the following lines into your main file (after `// update dynamic wall properties`):
Compile and run your new coupled problem and take a look at the results in Paraview.
In addition to the standard variables and parameters, you can now choose turbulence model specific quantities (e.g. the turbulent viscosity `nu_t`) for the free flow domain.
The result for the turbulent viscosity should look like this:
__Task B__:
Instead of computing the whole cross-section of a channel, you can use symmetric boundary conditions at the top boundary of your free flow domain by replacing all previous boundary conditions with
In addition, you have to remove the condition `onUpperBoundary_(globalPos)` from the `isOnWall(globalPos)` method.
__Task C__:
Choose `Surface With Edges` instead of `Surface` in the Paraview toolbar to see the discretization grid.
We will refine the grid now in order to better resolve the processes at the coupling interface.
Since not much is happening at the upper and lower boundaries of the whole domain, we want to keep the resolution low in these areas to save some computation time.
A grid refinement towards the interface is called _grading_.
Try different gradings by changing the values in the respective sections in the input file:
Grading0 = 1.0
Grading1 = 1.0
__Task D__:
For the grid convergence study, run various simulations with the following grading parameters:
* [Stokes.Grid] Grading1 = 1.2, [Darcy.Grid] Grading1 = -1.2
* [Stokes.Grid] Grading1 = 1.3, [Darcy.Grid] Grading1 = -1.3
* [Stokes.Grid] Grading1 = 1.4, [Darcy.Grid] Grading1 = -1.4
* [Stokes.Grid] Grading1 = 1.5, [Darcy.Grid] Grading1 = -1.5
* [Stokes.Grid] Grading1 = 1.6, [Darcy.Grid] Grading1 = -1.6
By changing the parameter `Problem.Name` for each grading factor, you avoid loosing the `.vtu` and `.pvd` files of the previous simulation runs.
Check the first lines of the output to see how the grading factors change the height of your grid cells.
Compare the velocity fields with Paraview.
# executables for ex_interface_coupling_ff-pm
dune_add_test(NAME ex_turbulence_coupling_ff-pm
CMD_ARGS ex_turbulence_coupling_ff-pm.input)
# add tutorial to the common target
add_dependencies(test_exercises ex_turbulence_coupling_ff-pm)
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <>. *
* \file
* \brief A test problem for the isothermal coupled Stokes/Darcy problem (1p2c/2p2c)
#include <config.h>
#include <ctime>
#include <iostream>
#include <fstream>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/timer.hh>
#include <dune/istl/io.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/dumuxmessage.hh>
#include <dumux/common/geometry/diameter.hh>
#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/assembly/fvassembler.hh>
#include <dumux/assembly/diffmethod.hh>
#include <dumux/discretization/methods.hh>
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/io/staggeredvtkoutputmodule.hh>
#include <dumux/io/grid/gridmanager.hh>
#include <dumux/multidomain/staggeredtraits.hh>
#include <dumux/multidomain/fvassembler.hh>
#include <dumux/multidomain/privarswitchnewtonsolver.hh>
#include <dumux/multidomain/boundary/stokesdarcy/couplingmanager.hh>
#include "ex_turbulence_pmproblem.hh"
#include "ex_turbulence_ffproblem.hh"
namespace Dumux {
namespace Properties {
SET_PROP(ZeroEqTypeTag, CouplingManager)
using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, TTAG(DarcyTwoPTwoCTypeTag)>;
using type = Dumux::StokesDarcyCouplingManager<Traits>;
SET_PROP(DarcyTwoPTwoCTypeTag, CouplingManager)
using Traits = StaggeredMultiDomainTraits<TTAG(ZeroEqTypeTag), TTAG(ZeroEqTypeTag), TypeTag>;
using type = Dumux::StokesDarcyCouplingManager<Traits>;
} // end namespace Properties
} // end namespace Dumux
int main(int argc, char** argv) try
using namespace Dumux;
// initialize MPI, finalize is done automatically on exit
const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
// print dumux start message
if (mpiHelper.rank() == 0)
// parse command line arguments and input file
Parameters::init(argc, argv);
// Define the sub problem type tags
using StokesTypeTag = TTAG(ZeroEqTypeTag);
using DarcyTypeTag = TTAG(DarcyTwoPTwoCTypeTag);
// try to create a grid (from the given grid file or the input file)
// for both sub-domains
using DarcyGridManager = Dumux::GridManager<typename GET_PROP_TYPE(DarcyTypeTag, Grid)>;
DarcyGridManager darcyGridManager;
darcyGridManager.init("Darcy"); // pass parameter group
using StokesGridManager = Dumux::GridManager<typename GET_PROP_TYPE(StokesTypeTag, Grid)>;
StokesGridManager stokesGridManager;
stokesGridManager.init("Stokes"); // pass parameter group
// we compute on the leaf grid view
const auto& darcyGridView = darcyGridManager.grid().leafGridView();
const auto& stokesGridView = stokesGridManager.grid().leafGridView();
// create the finite volume grid geometry
using StokesFVGridGeometry = typename GET_PROP_TYPE(StokesTypeTag, FVGridGeometry);
auto stokesFvGridGeometry = std::make_shared<StokesFVGridGeometry>(stokesGridView);
using DarcyFVGridGeometry = typename GET_PROP_TYPE(DarcyTypeTag, FVGridGeometry);
auto darcyFvGridGeometry = std::make_shared<DarcyFVGridGeometry>(darcyGridView);
using Traits = StaggeredMultiDomainTraits<StokesTypeTag, StokesTypeTag, DarcyTypeTag>;
// the coupling manager
using CouplingManager = StokesDarcyCouplingManager<Traits>;
auto couplingManager = std::make_shared<CouplingManager>(stokesFvGridGeometry, darcyFvGridGeometry);
// the indices
constexpr auto stokesCellCenterIdx = CouplingManager::stokesCellCenterIdx;
constexpr auto stokesFaceIdx = CouplingManager::stokesFaceIdx;
constexpr auto darcyIdx = CouplingManager::darcyIdx;
// the problem (initial and boundary conditions)
using StokesProblem = typename GET_PROP_TYPE(StokesTypeTag, Problem);
auto stokesProblem = std::make_shared<StokesProblem>(stokesFvGridGeometry, couplingManager);
using DarcyProblem = typename GET_PROP_TYPE(DarcyTypeTag, Problem);
auto darcyProblem = std::make_shared<DarcyProblem>(darcyFvGridGeometry, couplingManager);
// initialize the fluidsystem (tabulation)
GET_PROP_TYPE(StokesTypeTag, FluidSystem)::init();
// get some time loop parameters
using Scalar = typename GET_PROP_TYPE(StokesTypeTag, Scalar);
const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
auto dt = getParam<Scalar>("TimeLoop.DtInitial");
// check if we are about to restart a previously interrupted simulation
Scalar restartTime = 0;
if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart"))
restartTime = getParam<Scalar>("TimeLoop.Restart");
// instantiate time loop
auto timeLoop = std::make_shared<TimeLoop<Scalar>>(restartTime, dt, tEnd);
// set timeloop for the subproblems, needed for boundary value variations
// the solution vector
Traits::SolutionVector sol;
// apply initial solution for instationary problems
// auxiliary free flow solution vector
typename GET_PROP_TYPE(StokesTypeTag, SolutionVector) stokesSol;
auto solStokesOld = stokesSol;
sol[stokesCellCenterIdx] = stokesSol[stokesCellCenterIdx];
sol[stokesFaceIdx] = stokesSol[stokesFaceIdx];
// TODO: update static wall properties
// TODO: update dynamic wall properties
auto solDarcyOld = sol[darcyIdx];
auto solOld = sol;
couplingManager->init(stokesProblem, darcyProblem, sol);
// the grid variables
using StokesGridVariables = typename GET_PROP_TYPE(StokesTypeTag, GridVariables);
auto stokesGridVariables = std::make_shared<StokesGridVariables>(stokesProblem, stokesFvGridGeometry);
stokesGridVariables->init(stokesSol, solStokesOld);
using DarcyGridVariables = typename GET_PROP_TYPE(DarcyTypeTag, GridVariables);
auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyFvGridGeometry);
darcyGridVariables->init(sol[darcyIdx], solDarcyOld);
// intialize the vtk output module
const auto stokesName = getParam<std::string>("Problem.Name") + "_" + stokesProblem->name();
const auto darcyName = getParam<std::string>("Problem.Name") + "_" + darcyProblem->name();
StaggeredVtkOutputModule<StokesTypeTag, GET_PROP_VALUE(StokesTypeTag, PhaseIdx)> stokesVtkWriter(*stokesProblem, *stokesFvGridGeometry, *stokesGridVariables, stokesSol, stokesName);
GET_PROP_TYPE(StokesTypeTag, VtkOutputFields)::init(stokesVtkWriter);
VtkOutputModule<DarcyTypeTag> darcyVtkWriter(*darcyProblem, *darcyFvGridGeometry, *darcyGridVariables, sol[darcyIdx], darcyName);
GET_PROP_TYPE(DarcyTypeTag, VtkOutputFields)::init(darcyVtkWriter);
// the assembler with time loop for instationary problem
using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
auto assembler = std::make_shared<Assembler>(std::make_tuple(stokesProblem, stokesProblem, darcyProblem),
// the linear solver
using LinearSolver = UMFPackBackend;
auto linearSolver = std::make_shared<LinearSolver>();
// the primary variable switches used by the sub models
using PriVarSwitchTuple = std::tuple<NoPrimaryVariableSwitch, NoPrimaryVariableSwitch, typename GET_PROP_TYPE(DarcyTypeTag, PrimaryVariableSwitch)>;
// the non-linear solver
using NewtonSolver = MultiDomainPriVarSwitchNewtonSolver<Assembler, LinearSolver, CouplingManager, PriVarSwitchTuple>;
NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
// time loop
timeLoop->start(); do
// set previous solution for storage evaluations
// solve the non-linear system with time step control
nonLinearSolver.solve(sol, *timeLoop);
// make the new solution the old solution
solOld = sol;
// update the auxiliary free flow solution vector
stokesSol[stokesCellCenterIdx] = sol[stokesCellCenterIdx];
stokesSol[stokesFaceIdx] = sol[stokesFaceIdx];
// TODO: update dynamic wall properties
// post time step treatment of Darcy problem
darcyProblem->postTimeStep(sol[darcyIdx], *darcyGridVariables, timeLoop->timeStepSize());
// advance grid variables to the next time step
// advance to the time loop to the next step
// write vtk output
// report statistics of this time step
// set new dt as suggested by newton solver
} while (!timeLoop->finished());
// finalize, print dumux message to say goodbye
// print dumux end message
if (mpiHelper.rank() == 0)
return 0;
} // end main
catch (Dumux::ParameterException &e)
std::cerr << std::endl << e << " ---> Abort!" << std::endl;
return 1;
catch (Dune::DGFException & e)
std::cerr << "DGF exception thrown (" << e <<
"). Most likely, the DGF file name is wrong "
"or the DGF file is corrupted, "
"e.g. missing hash at end of file or wrong number (dimensions) of entries."
<< " ---> Abort!" << std::endl;
return 2;
catch (Dune::Exception &e)
std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
return 3;
catch (...)
std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
return 4;
DtInitial = 1e-1 # [s]
MaxTimeStepSize = 43200 # [s] (12 hours)
TEnd = 864000 # [s] (6 days)
Positions0 = 0.0 0.25
Positions1 = 0.25 0.5
Grading0 = 1.0
Grading1 = 1.0
Cells0 = 15
Cells1 = 20
Verbosity = true
Positions0 = 0.0 0.25
Positions1 = 0.0 0.25
Cells0 = 15
Cells1 = 10
Grading0 = 1.0
Grading1 = 1.0
Verbosity = true
Name = stokes
RefVelocity = 3.5 # [m/s]
RefPressure = 1e5 # [Pa]
refMoleFrac = 0 # [-]
RefTemperature = 298.15 # [K]
Name = darcy
Pressure = 1e5
Saturation = 0.5 # initial Sw
Temperature = 298.15 # [K]
InitPhasePresence = 3 # bothPhases
Porosity = 0.41
Permeability = 2.65e-10
AlphaBeaversJoseph = 1.0
Swr = 0.005
Snr = 0.01
VgAlpha = 6.371e-4
VgN = 6.9
Name = ex_coupling_turbulence_ff-pm
EnableGravity = true
InterfaceDiffusionCoefficientAvg = Harmonic
AddVelocity = true
WriteFaceData = false
MaxSteps = 12
MaxRelativeShift = 1e-5
NumericDifferenceMethod = 0
NumericDifference.BaseEpsilon = 1e-8
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <>. *
* \file
* \brief The free-flow sub problem
#include <dune/grid/yaspgrid.hh>
#include <dumux/material/fluidsystems/h2oair.hh>
#include <dumux/discretization/staggered/freeflow/properties.hh>
#include <dumux/freeflow/compositional/navierstokesncmodel.hh>
#include <dumux/freeflow/navierstokes/problem.hh>
namespace Dumux
template <class TypeTag>
class FreeFlowSubProblem;
namespace Properties
NEW_TYPE_TAG(ZeroEqTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, NavierStokesNCNI));
// Set the grid type
SET_TYPE_PROP(ZeroEqTypeTag, Grid, Dune::YaspGrid<2, Dune::TensorProductCoordinates<typename GET_PROP_TYPE(TypeTag, Scalar), 2> >);
// set the fluid system
SET_TYPE_PROP(ZeroEqTypeTag, FluidSystem, FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>);
// set phase index (air)
SET_INT_PROP(ZeroEqTypeTag, PhaseIdx, GET_PROP_TYPE(TypeTag, FluidSystem)::gasPhaseIdx);
SET_INT_PROP(ZeroEqTypeTag, ReplaceCompEqIdx, 3);
// Use formulation based on mass fractions
SET_BOOL_PROP(ZeroEqTypeTag, UseMoles, true);
// Set the problem property
SET_TYPE_PROP(ZeroEqTypeTag, Problem, Dumux::FreeFlowSubProblem<TypeTag> );
SET_BOOL_PROP(ZeroEqTypeTag, EnableFVGridGeometryCache, true);
SET_BOOL_PROP(ZeroEqTypeTag, EnableGridFluxVariablesCache, true);
SET_BOOL_PROP(ZeroEqTypeTag, EnableGridVolumeVariablesCache, true);
SET_BOOL_PROP(ZeroEqTypeTag, EnableInertiaTerms, true);
* \brief The free-flow sub problem
template <class TypeTag>
class FreeFlowSubProblem : public NavierStokesProblem<TypeTag>
using ParentType = NavierStokesProblem<TypeTag>;
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using FVElementGeometry = typename FVGridGeometry::LocalView;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using Element = typename GridView::template Codim<0>::Entity;
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
using ElementFaceVariables = typename GET_PROP_TYPE(TypeTag, GridFaceVariables)::LocalView;
using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);