From e53d6fab3235c3b63b6ef8805e203c3e43c0de1e Mon Sep 17 00:00:00 2001
From: Kilian Weishaupt <kilian.weishaupt@iws.uni-stuttgart.de>
Date: Tue, 11 Dec 2018 18:46:33 +0100
Subject: [PATCH] [ff-pnm] Update second exercise
---
exercises/exercise-coupling-ff-pm/README.md | 8 +-
.../models/ex_models_coupling_ff-pm.cc | 73 +++++------
.../models/ex_models_coupling_ff-pm.input | 13 +-
.../models/ex_models_ffproblem.hh | 87 +++++++------
.../models/ex_models_pmproblem.hh | 73 +++++------
.../models/ex_models_coupling_ff-pm.cc | 84 +++++-------
.../models/ex_models_coupling_ff-pm.input | 5 +-
.../models/ex_models_ffproblem.hh | 123 ++++++++++--------
.../models/ex_models_pmproblem.hh | 82 ++++++------
9 files changed, 266 insertions(+), 282 deletions(-)
diff --git a/exercises/exercise-coupling-ff-pm/README.md b/exercises/exercise-coupling-ff-pm/README.md
index ec4e22aa..f3d56c7e 100644
--- a/exercises/exercise-coupling-ff-pm/README.md
+++ b/exercises/exercise-coupling-ff-pm/README.md
@@ -269,7 +269,7 @@ Finally we want to know the distribution of the water mass fluxes across the int
Use the facilities therein to return the values of ...massCouplingCondition... from the `couplingManager`
for each coupling scvf. Then the fluxes are visualized with gnuplot, when setting `Problem.PlotFluxes = true`.
If the simulation is too fast, you can have a look at the flux*.png files after the simulation.
-* You can use the property `Problem.PlotStorage = true` to see the temporal evolution of the evaporation rate
+* You can use the parameter `Problem.PlotStorage = true` to see the temporal evolution of the evaporation rate
and the cumulative water mass loss.
Compile and run the simulation and take a look at the results.
@@ -295,10 +295,6 @@ for this case the liquid saturation cannot serve as primary variable anymore. Ho
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.
-- Replace the current implementation of the Newton solver with the version which can handle
- primary variable switches: `dumux/multidomain/privarswitchnewtonsolver.hh`.
- You also have to uncomment the line containing the `PriVarSwitchTuple` and to overwrite
- the last argument with the `PrimaryVariableSwitch` property from the Darcy model.
Now you are able to simulate a complete drying of the porous medium.
@@ -355,7 +351,7 @@ However, there is also a solution-dependent component of these interactions, e.g
damping of the eddy viscosity toward the wall, the velocity gradient at the wall and inside the
cells is needed.
These dynamic interactions are to be updated by calling
-```cpp
+```cpp
stokesProblem->updateDynamicWallProperties(stokesSol)
```
in the time loop (after `// update dynamic wall properties`).
diff --git a/exercises/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.cc b/exercises/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.cc
index b54ae19e..8dc210c5 100644
--- a/exercises/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.cc
+++ b/exercises/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.cc
@@ -34,11 +34,11 @@
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/dumuxmessage.hh>
-#include <dumux/common/geometry/diameter.hh>
+#include <dumux/common/partial.hh>
#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/assembly/fvassembler.hh>
#include <dumux/assembly/diffmethod.hh>
-#include <dumux/discretization/methods.hh>
+#include <dumux/discretization/method.hh>
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/io/staggeredvtkoutputmodule.hh>
#include <dumux/io/grid/gridmanager.hh>
@@ -55,15 +55,17 @@
namespace Dumux {
namespace Properties {
-SET_PROP(StokesTypeTag, CouplingManager)
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::StokesTypeTag>
{
- using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, TTAG(DarcyTypeTag)>;
+ using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, Properties::TTag::DarcyTypeTag>;
using type = Dumux::StokesDarcyCouplingManager<Traits>;
};
-SET_PROP(DarcyTypeTag, CouplingManager)
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::DarcyTypeTag>
{
- using Traits = StaggeredMultiDomainTraits<TTAG(StokesTypeTag), TTAG(StokesTypeTag), TypeTag>;
+ using Traits = StaggeredMultiDomainTraits<Properties::TTag::StokesTypeTag, Properties::TTag::StokesTypeTag, TypeTag>;
using type = Dumux::StokesDarcyCouplingManager<Traits>;
};
@@ -85,16 +87,16 @@ int main(int argc, char** argv) try
Parameters::init(argc, argv);
// Define the sub problem type tags
- using StokesTypeTag = TTAG(StokesTypeTag);
- using DarcyTypeTag = TTAG(DarcyTypeTag);
+ using StokesTypeTag = Properties::TTag::StokesTypeTag;
+ using DarcyTypeTag = Properties::TTag::DarcyTypeTag;
// try to create a grid (from the given grid file or the input file)
// for both sub-domains
- using DarcyGridManager = Dumux::GridManager<typename GET_PROP_TYPE(DarcyTypeTag, Grid)>;
+ using DarcyGridManager = Dumux::GridManager<GetPropType<DarcyTypeTag, Properties::Grid>>;
DarcyGridManager darcyGridManager;
darcyGridManager.init("Darcy"); // pass parameter group
- using StokesGridManager = Dumux::GridManager<typename GET_PROP_TYPE(StokesTypeTag, Grid)>;
+ using StokesGridManager = Dumux::GridManager<GetPropType<StokesTypeTag, Properties::Grid>>;
StokesGridManager stokesGridManager;
stokesGridManager.init("Stokes"); // pass parameter group
@@ -103,10 +105,10 @@ int main(int argc, char** argv) try
const auto& stokesGridView = stokesGridManager.grid().leafGridView();
// create the finite volume grid geometry
- using StokesFVGridGeometry = typename GET_PROP_TYPE(StokesTypeTag, FVGridGeometry);
+ using StokesFVGridGeometry = GetPropType<StokesTypeTag, Properties::FVGridGeometry>;
auto stokesFvGridGeometry = std::make_shared<StokesFVGridGeometry>(stokesGridView);
stokesFvGridGeometry->update();
- using DarcyFVGridGeometry = typename GET_PROP_TYPE(DarcyTypeTag, FVGridGeometry);
+ using DarcyFVGridGeometry = GetPropType<DarcyTypeTag, Properties::FVGridGeometry>;
auto darcyFvGridGeometry = std::make_shared<DarcyFVGridGeometry>(darcyGridView);
darcyFvGridGeometry->update();
@@ -122,27 +124,22 @@ int main(int argc, char** argv) try
constexpr auto darcyIdx = CouplingManager::darcyIdx;
// the problem (initial and boundary conditions)
- using StokesProblem = typename GET_PROP_TYPE(StokesTypeTag, Problem);
+ using StokesProblem = GetPropType<StokesTypeTag, Properties::Problem>;
auto stokesProblem = std::make_shared<StokesProblem>(stokesFvGridGeometry, couplingManager);
- using DarcyProblem = typename GET_PROP_TYPE(DarcyTypeTag, Problem);
+ using DarcyProblem = GetPropType<DarcyTypeTag, Properties::Problem>;
auto darcyProblem = std::make_shared<DarcyProblem>(darcyFvGridGeometry, couplingManager);
// initialize the fluidsystem (tabulation)
- GET_PROP_TYPE(StokesTypeTag, FluidSystem)::init();
+ GetPropType<StokesTypeTag, Properties::FluidSystem>::init();
// get some time loop parameters
- using Scalar = typename GET_PROP_TYPE(StokesTypeTag, Scalar);
+ using Scalar = GetPropType<StokesTypeTag, Properties::Scalar>;
const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
auto dt = getParam<Scalar>("TimeLoop.DtInitial");
- // check if we are about to restart a previously interrupted simulation
- Scalar restartTime = 0;
- if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart"))
- restartTime = getParam<Scalar>("TimeLoop.Restart");
-
// instantiate time loop
- auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(restartTime, dt, tEnd);
+ auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(0, dt, tEnd);
timeLoop->setMaxTimeStepSize(maxDt);
stokesProblem->setTimeLoop(timeLoop);
darcyProblem->setTimeLoop(timeLoop);
@@ -153,43 +150,35 @@ int main(int argc, char** argv) try
sol[stokesFaceIdx].resize(stokesFvGridGeometry->numFaceDofs());
sol[darcyIdx].resize(darcyFvGridGeometry->numDofs());
- const auto& cellCenterSol = sol[stokesCellCenterIdx];
- const auto& faceSol = sol[stokesFaceIdx];
+ auto stokesSol = partial(sol, stokesCellCenterIdx, stokesFaceIdx);
- // apply initial solution for instationary problems
- typename GET_PROP_TYPE(StokesTypeTag, SolutionVector) stokesSol;
- std::get<0>(stokesSol) = cellCenterSol;
- std::get<1>(stokesSol) = faceSol;
stokesProblem->applyInitialSolution(stokesSol);
- auto solStokesOld = stokesSol;
- sol[stokesCellCenterIdx] = stokesSol[stokesCellCenterIdx];
- sol[stokesFaceIdx] = stokesSol[stokesFaceIdx];
-
darcyProblem->applyInitialSolution(sol[darcyIdx]);
- auto solDarcyOld = sol[darcyIdx];
auto solOld = sol;
couplingManager->init(stokesProblem, darcyProblem, sol);
// the grid variables
- using StokesGridVariables = typename GET_PROP_TYPE(StokesTypeTag, GridVariables);
+ using StokesGridVariables = GetPropType<StokesTypeTag, Properties::GridVariables>;
auto stokesGridVariables = std::make_shared<StokesGridVariables>(stokesProblem, stokesFvGridGeometry);
- stokesGridVariables->init(stokesSol, solStokesOld);
- using DarcyGridVariables = typename GET_PROP_TYPE(DarcyTypeTag, GridVariables);
+ stokesGridVariables->init(stokesSol);
+ using DarcyGridVariables = GetPropType<DarcyTypeTag, Properties::GridVariables>;
auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyFvGridGeometry);
- darcyGridVariables->init(sol[darcyIdx], solDarcyOld);
+ darcyGridVariables->init(sol[darcyIdx]);
// intialize the vtk output module
const auto stokesName = getParam<std::string>("Problem.Name") + "_" + stokesProblem->name();
const auto darcyName = getParam<std::string>("Problem.Name") + "_" + darcyProblem->name();
- StaggeredVtkOutputModule<StokesTypeTag> stokesVtkWriter(*stokesProblem, *stokesFvGridGeometry, *stokesGridVariables, stokesSol, stokesName);
- GET_PROP_TYPE(StokesTypeTag, VtkOutputFields)::init(stokesVtkWriter);
+ StaggeredVtkOutputModule<StokesGridVariables, decltype(stokesSol)> stokesVtkWriter(*stokesGridVariables, stokesSol, stokesName);
+ GetPropType<StokesTypeTag, Properties::VtkOutputFields>::initOutputModule(stokesVtkWriter);
stokesVtkWriter.write(0.0);
- VtkOutputModule<DarcyTypeTag> darcyVtkWriter(*darcyProblem, *darcyFvGridGeometry, *darcyGridVariables, sol[darcyIdx], darcyName);
- GET_PROP_TYPE(DarcyTypeTag, VtkOutputFields)::init(darcyVtkWriter);
+ VtkOutputModule<DarcyGridVariables, GetPropType<DarcyTypeTag, Properties::SolutionVector>> darcyVtkWriter(*darcyGridVariables, sol[darcyIdx], darcyName);
+ using DarcyVelocityOutput = GetPropType<DarcyTypeTag, Properties::VelocityOutput>;
+ darcyVtkWriter.addVelocityOutput(std::make_shared<DarcyVelocityOutput>(*darcyGridVariables));
+ GetPropType<DarcyTypeTag, Properties::IOFields>::initOutputModule(darcyVtkWriter);
darcyVtkWriter.write(0.0);
// intialize the subproblems
@@ -211,8 +200,6 @@ int main(int argc, char** argv) try
using LinearSolver = UMFPackBackend;
auto linearSolver = std::make_shared<LinearSolver>();
- // the primary variable switches used by the sub models and the non-linear solver
-// using PriVarSwitchTuple = std::tuple<NoPrimaryVariableSwitch, NoPrimaryVariableSwitch, NoPrimaryVariableSwitch>;
using NewtonSolver = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
diff --git a/exercises/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.input b/exercises/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.input
index d3f68518..07bb4dd0 100644
--- a/exercises/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.input
+++ b/exercises/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.input
@@ -1,5 +1,5 @@
[TimeLoop]
-DtInitial = 100 # s
+DtInitial = 1 # s
EpisodeLength = -360 # s # 0.25 days
TEnd = 256000 # s # 2 days
@@ -14,10 +14,11 @@ Cells = 16 16
[Stokes.Problem]
Name = stokes
-EnableGravity = false
+EnableGravity = true
+EnableInertiaTerms = true
MoleFraction = 0.0 # -
Pressure = 1e5 # Pa
-Velocity = 1e-3 # m/s
+Velocity = 1 # m/s
[Darcy.Problem]
Name = darcy
@@ -43,7 +44,11 @@ PlotStorage = false
[Newton]
MaxSteps = 12
-MaxRelativeShift = 1e-5
+MaxRelativeShift = 1e-7
+TargetSteps = 5
[Vtk]
AddVelocity = 1
+
+[Assembly]
+NumericDifference.BaseEpsilon = 1e-8
diff --git a/exercises/exercise-coupling-ff-pm/models/ex_models_ffproblem.hh b/exercises/exercise-coupling-ff-pm/models/ex_models_ffproblem.hh
index 1fe98c08..6f00aa10 100644
--- a/exercises/exercise-coupling-ff-pm/models/ex_models_ffproblem.hh
+++ b/exercises/exercise-coupling-ff-pm/models/ex_models_ffproblem.hh
@@ -40,32 +40,41 @@ class StokesSubProblem;
namespace Properties
{
-NEW_TYPE_TAG(StokesTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, NavierStokesNC));
+// Create new type tags
+namespace TTag {
+struct StokesTypeTag { using InheritsFrom = std::tuple<NavierStokesNC, StaggeredFreeFlowModel>; };
+} // end namespace TTag
// Set the grid type
-SET_TYPE_PROP(StokesTypeTag, Grid, Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<typename GET_PROP_TYPE(TypeTag, Scalar), 2> >);
+template<class TypeTag>
+struct Grid<TypeTag, TTag::StokesTypeTag> { using type = Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
// The fluid system
-SET_PROP(StokesTypeTag, FluidSystem)
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::StokesTypeTag>
{
- using H2OAir = FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>;
+ using H2OAir = FluidSystems::H2OAir<GetPropType<TypeTag, Properties::Scalar>>;
using type = FluidSystems::OnePAdapter<H2OAir, H2OAir::gasPhaseIdx>;
};
// Do not replace one equation with a total mass balance
-SET_INT_PROP(StokesTypeTag, ReplaceCompEqIdx, 3);
+template<class TypeTag>
+struct ReplaceCompEqIdx<TypeTag, TTag::StokesTypeTag> { static constexpr int value = 3; };
// Use formulation based on mass fractions
-SET_BOOL_PROP(StokesTypeTag, UseMoles, true);
+template<class TypeTag>
+struct UseMoles<TypeTag, TTag::StokesTypeTag> { static constexpr bool value = true; };
// Set the problem property
-SET_TYPE_PROP(StokesTypeTag, Problem, Dumux::StokesSubProblem<TypeTag> );
-
-SET_BOOL_PROP(StokesTypeTag, EnableFVGridGeometryCache, true);
-SET_BOOL_PROP(StokesTypeTag, EnableGridFluxVariablesCache, true);
-SET_BOOL_PROP(StokesTypeTag, EnableGridVolumeVariablesCache, true);
-
-SET_BOOL_PROP(StokesTypeTag, EnableInertiaTerms, false);
+template<class TypeTag>
+struct Problem<TypeTag, TTag::StokesTypeTag> { using type = Dumux::StokesSubProblem<TypeTag> ; };
+
+template<class TypeTag>
+struct EnableFVGridGeometryCache<TypeTag, TTag::StokesTypeTag> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridFluxVariablesCache<TypeTag, TTag::StokesTypeTag> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridVolumeVariablesCache<TypeTag, TTag::StokesTypeTag> { static constexpr bool value = true; };
}
/*!
@@ -79,29 +88,29 @@ class StokesSubProblem : public NavierStokesProblem<TypeTag>
{
using ParentType = NavierStokesProblem<TypeTag>;
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
- using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
- using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+ using GridView = GetPropType<TypeTag, Properties::GridView>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+ using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+ using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
- using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
using FVElementGeometry = typename FVGridGeometry::LocalView;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using Element = typename GridView::template Codim<0>::Entity;
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
- using ElementFaceVariables = typename GET_PROP_TYPE(TypeTag, GridFaceVariables)::LocalView;
- using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
+ using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
+ using ElementFaceVariables = typename GetPropType<TypeTag, Properties::GridFaceVariables>::LocalView;
+ using FluidState = GetPropType<TypeTag, Properties::FluidState>;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
- using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
- using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+ using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
- using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager);
+ using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
using TimeLoopPtr = std::shared_ptr<TimeLoop<Scalar>>;
- static constexpr bool useMoles = GET_PROP_TYPE(TypeTag, ModelTraits)::useMoles();
+ static constexpr bool useMoles = GetPropType<TypeTag, Properties::ModelTraits>::useMoles();
public:
StokesSubProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry, std::shared_ptr<CouplingManager> couplingManager)
@@ -117,10 +126,6 @@ public:
*/
// \{
-
- bool shouldWriteRestartFile() const
- { return false; }
-
/*!
* \brief Return the temperature within the domain in [K].
*/
@@ -218,9 +223,9 @@ public:
if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
{
- values[Indices::momentumYBalanceIdx] = couplingManager().couplingData().momentumCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf);
+ values[Indices::momentumYBalanceIdx] = couplingManager().couplingData().momentumCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf);
- const auto massFlux = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf);
+ const auto massFlux = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf);
values[Indices::conti0EqIdx] = massFlux[0];
values[Indices::conti0EqIdx + 1] = massFlux[1];
}
@@ -251,14 +256,22 @@ public:
*
* \param globalPos The global position
*/
- PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+ PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
{
+ PrimaryVariables values(0.0);
+
+ // This is only an approximation of the actual hydorostatic pressure gradient.
+ // Air is compressible (the density depends on pressure), thus the actual
+ // vertical pressure profile is non-linear.
+ // This discrepancy can lead to spurious flows at the outlet if the inlet
+ // velocity is chosen too small.
FluidState fluidState;
updateFluidStateForBC_(fluidState, pressure_);
const Scalar density = FluidSystem::density(fluidState, 0);
+ values[Indices::pressureIdx] = pressure_ - density*this->gravity()[1]*(this->fvGridGeometry().bBoxMax()[1] - globalPos[1]);
- PrimaryVariables values(0.0);
- values[Indices::pressureIdx] = pressure_ + density*this->gravity()[1]*(globalPos[1] - this->fvGridGeometry().bBoxMin()[1]);
+
+ // gravity has negative sign
values[Indices::conti0EqIdx + 1] = moleFraction_;
values[Indices::velocityXIdx] = 4.0 * velocity_ * (globalPos[1] - this->fvGridGeometry().bBoxMin()[1])
* (this->fvGridGeometry().bBoxMax()[1] - globalPos[1])
@@ -273,9 +286,9 @@ public:
/*!
* \brief Returns the intrinsic permeability of required as input parameter for the Beavers-Joseph-Saffman boundary condition
*/
- Scalar permeability(const SubControlVolumeFace& scvf) const
+ Scalar permeability(const Element& element, const SubControlVolumeFace& scvf) const
{
- return couplingManager().couplingData().darcyPermeability(scvf);
+ return couplingManager().couplingData().darcyPermeability(element, scvf);
}
/*!
diff --git a/exercises/exercise-coupling-ff-pm/models/ex_models_pmproblem.hh b/exercises/exercise-coupling-ff-pm/models/ex_models_pmproblem.hh
index c20e9e42..fb4f258c 100644
--- a/exercises/exercise-coupling-ff-pm/models/ex_models_pmproblem.hh
+++ b/exercises/exercise-coupling-ff-pm/models/ex_models_pmproblem.hh
@@ -44,52 +44,61 @@ class DarcySubProblem;
namespace Properties
{
-NEW_TYPE_TAG(DarcyTypeTag, INHERITS_FROM(CCTpfaModel, OnePNC));
+// Create new type tags
+namespace TTag {
+struct DarcyTypeTag { using InheritsFrom = std::tuple<OnePNC, CCTpfaModel>; };
+} // end namespace TTag
// Set the problem property
-SET_TYPE_PROP(DarcyTypeTag, Problem, Dumux::DarcySubProblem<TypeTag>);
+template<class TypeTag>
+struct Problem<TypeTag, TTag::DarcyTypeTag> { using type = Dumux::DarcySubProblem<TypeTag>; };
// The fluid system
-SET_PROP(DarcyTypeTag, FluidSystem)
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::DarcyTypeTag>
{
- using H2OAir = FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>;
+ using H2OAir = FluidSystems::H2OAir<GetPropType<TypeTag, Properties::Scalar>>;
using type = FluidSystems::OnePAdapter<H2OAir, H2OAir::gasPhaseIdx>;
};
// Use moles
-SET_BOOL_PROP(DarcyTypeTag, UseMoles, true);
+template<class TypeTag>
+struct UseMoles<TypeTag, TTag::DarcyTypeTag> { static constexpr bool value = true; };
// Do not replace one equation with a total mass balance
-SET_INT_PROP(DarcyTypeTag, ReplaceCompEqIdx, 3);
+template<class TypeTag>
+struct ReplaceCompEqIdx<TypeTag, TTag::DarcyTypeTag> { static constexpr int value = 3; };
//! Use a model with constant tortuosity for the effective diffusivity
SET_TYPE_PROP(DarcyTypeTag, EffectiveDiffusivityModel,
- DiffusivityConstantTortuosity<typename GET_PROP_TYPE(TypeTag, Scalar)>);
+ DiffusivityConstantTortuosity<GetPropType<TypeTag, Properties::Scalar>>);
// Set the grid type
-SET_TYPE_PROP(DarcyTypeTag, Grid, Dune::YaspGrid<2>);
+template<class TypeTag>
+struct Grid<TypeTag, TTag::DarcyTypeTag> { using type = Dune::YaspGrid<2>; };
// Set the spatial paramaters type
-SET_TYPE_PROP(DarcyTypeTag, SpatialParams, OnePSpatialParams<TypeTag>);
+template<class TypeTag>
+struct SpatialParams<TypeTag, TTag::DarcyTypeTag> { using type = OnePSpatialParams<TypeTag>; };
}
template <class TypeTag>
class DarcySubProblem : public PorousMediumFlowProblem<TypeTag>
{
using ParentType = PorousMediumFlowProblem<TypeTag>;
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
- using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
- using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
- using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
+ using GridView = GetPropType<TypeTag, Properties::GridView>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+ using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
+ using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
+ using FVElementGeometry = typename GetPropType<TypeTag, Properties::FVGridGeometry>::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
- using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
// copy some indices for convenience
- using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+ using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
enum {
// grid and world dimension
dim = GridView::dimension,
@@ -105,7 +114,7 @@ class DarcySubProblem : public PorousMediumFlowProblem<TypeTag>
using Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
- using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager);
+ using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
using TimeLoopPtr = std::shared_ptr<TimeLoop<Scalar>>;
public:
@@ -114,7 +123,6 @@ public:
: ParentType(fvGridGeometry, "Darcy"), eps_(1e-7), couplingManager_(couplingManager)
{
moleFraction_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.MoleFraction");
- pressure_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Pressure");
// initialize output file
plotFluxes_ = getParamFromGroup<bool>(this->paramGroup(), "Problem.PlotFluxes", false);
@@ -224,8 +232,6 @@ public:
if (!couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
continue;
- // NOTE: binding the coupling context is necessary
- couplingManager_->bindCouplingContext(CouplingManager::darcyIdx, element);
NumEqVector flux(0.0); // use "massCouplingCondition" from the couplingManager here
x.push_back(scvf.center()[0]);
@@ -246,21 +252,11 @@ public:
gnuplotInterfaceFluxes_.plot("flux_" + std::to_string(timeLoop_->timeStepIndex()));
}
- /*!
- * \brief Returns true if a restart file should be written to
- * disk.
- */
- bool shouldWriteRestartFile() const
- { return false; }
-
/*!
* \name Problem parameters
*/
// \{
- bool shouldWriteOutput() const // define output
- { return true; }
-
/*!
* \brief Return the temperature within the domain in [K].
*
@@ -311,7 +307,7 @@ public:
NumEqVector values(0.0);
if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
- values = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, scvf);
+ values = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, scvf);
return values;
}
@@ -332,10 +328,10 @@ public:
* \param scv The subcontrolvolume
*/
template<class ElementVolumeVariables>
- NumEqVector source(const Element &element,
+ NumEqVector source(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
- const SubControlVolume &scv) const
+ const SubControlVolume& scv) const
{ return NumEqVector(0.0); }
// \}
@@ -348,11 +344,13 @@ public:
* For this method, the \a priVars parameter stores primary
* variables.
*/
- PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+ PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
{
+ static const Scalar stokesPressure = getParamFromGroup<Scalar>("Stokes", "Problem.Pressure");
+
PrimaryVariables values(0.0);
+ values[Indices::pressureIdx] = stokesPressure;
values[transportCompIdx] = moleFraction_;
- values[pressureIdx] = pressure_;
return values;
}
@@ -384,7 +382,6 @@ private:
Scalar eps_;
Scalar moleFraction_;
- Scalar pressure_;
Scalar initialWaterContent_ = 0.0;
Scalar lastWaterMass_ = 0.0;
diff --git a/exercises/solution/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.cc b/exercises/solution/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.cc
index cdb7abdd..8528aeba 100644
--- a/exercises/solution/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.cc
+++ b/exercises/solution/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.cc
@@ -34,22 +34,18 @@
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/dumuxmessage.hh>
-#include <dumux/common/geometry/diameter.hh>
+#include <dumux/common/partial.hh>
#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/assembly/fvassembler.hh>
#include <dumux/assembly/diffmethod.hh>
-#include <dumux/discretization/methods.hh>
+#include <dumux/discretization/method.hh>
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/io/staggeredvtkoutputmodule.hh>
#include <dumux/io/grid/gridmanager.hh>
#include <dumux/multidomain/staggeredtraits.hh>
#include <dumux/multidomain/fvassembler.hh>
-#if EXNUMBER >= 3
-#include <dumux/multidomain/privarswitchnewtonsolver.hh>
-#else
#include <dumux/multidomain/newtonsolver.hh>
-#endif
#include <dumux/multidomain/boundary/stokesdarcy/couplingmanager.hh>
@@ -59,15 +55,17 @@
namespace Dumux {
namespace Properties {
-SET_PROP(StokesTypeTag, CouplingManager)
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::StokesTypeTag>
{
- using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, TTAG(DarcyTypeTag)>;
+ using Traits = StaggeredMultiDomainTraits<TypeTag, TypeTag, Properties::TTag::DarcyTypeTag>;
using type = Dumux::StokesDarcyCouplingManager<Traits>;
};
-SET_PROP(DarcyTypeTag, CouplingManager)
+template<class TypeTag>
+struct CouplingManager<TypeTag, TTag::DarcyTypeTag>
{
- using Traits = StaggeredMultiDomainTraits<TTAG(StokesTypeTag), TTAG(StokesTypeTag), TypeTag>;
+ using Traits = StaggeredMultiDomainTraits<Properties::TTag::StokesTypeTag, Properties::TTag::StokesTypeTag, TypeTag>;
using type = Dumux::StokesDarcyCouplingManager<Traits>;
};
@@ -89,16 +87,16 @@ int main(int argc, char** argv) try
Parameters::init(argc, argv);
// Define the sub problem type tags
- using StokesTypeTag = TTAG(StokesTypeTag);
- using DarcyTypeTag = TTAG(DarcyTypeTag);
+ using StokesTypeTag = Properties::TTag::StokesTypeTag;
+ using DarcyTypeTag = Properties::TTag::DarcyTypeTag;
// try to create a grid (from the given grid file or the input file)
// for both sub-domains
- using DarcyGridManager = Dumux::GridManager<typename GET_PROP_TYPE(DarcyTypeTag, Grid)>;
+ using DarcyGridManager = Dumux::GridManager<GetPropType<DarcyTypeTag, Properties::Grid>>;
DarcyGridManager darcyGridManager;
darcyGridManager.init("Darcy"); // pass parameter group
- using StokesGridManager = Dumux::GridManager<typename GET_PROP_TYPE(StokesTypeTag, Grid)>;
+ using StokesGridManager = Dumux::GridManager<GetPropType<StokesTypeTag, Properties::Grid>>;
StokesGridManager stokesGridManager;
stokesGridManager.init("Stokes"); // pass parameter group
@@ -107,10 +105,10 @@ int main(int argc, char** argv) try
const auto& stokesGridView = stokesGridManager.grid().leafGridView();
// create the finite volume grid geometry
- using StokesFVGridGeometry = typename GET_PROP_TYPE(StokesTypeTag, FVGridGeometry);
+ using StokesFVGridGeometry = GetPropType<StokesTypeTag, Properties::FVGridGeometry>;
auto stokesFvGridGeometry = std::make_shared<StokesFVGridGeometry>(stokesGridView);
stokesFvGridGeometry->update();
- using DarcyFVGridGeometry = typename GET_PROP_TYPE(DarcyTypeTag, FVGridGeometry);
+ using DarcyFVGridGeometry = GetPropType<DarcyTypeTag, Properties::FVGridGeometry>;
auto darcyFvGridGeometry = std::make_shared<DarcyFVGridGeometry>(darcyGridView);
darcyFvGridGeometry->update();
@@ -126,27 +124,22 @@ int main(int argc, char** argv) try
constexpr auto darcyIdx = CouplingManager::darcyIdx;
// the problem (initial and boundary conditions)
- using StokesProblem = typename GET_PROP_TYPE(StokesTypeTag, Problem);
+ using StokesProblem = GetPropType<StokesTypeTag, Properties::Problem>;
auto stokesProblem = std::make_shared<StokesProblem>(stokesFvGridGeometry, couplingManager);
- using DarcyProblem = typename GET_PROP_TYPE(DarcyTypeTag, Problem);
+ using DarcyProblem = GetPropType<DarcyTypeTag, Properties::Problem>;
auto darcyProblem = std::make_shared<DarcyProblem>(darcyFvGridGeometry, couplingManager);
// initialize the fluidsystem (tabulation)
- GET_PROP_TYPE(StokesTypeTag, FluidSystem)::init();
+ GetPropType<StokesTypeTag, Properties::FluidSystem>::init();
// get some time loop parameters
- using Scalar = typename GET_PROP_TYPE(StokesTypeTag, Scalar);
+ using Scalar = GetPropType<StokesTypeTag, Properties::Scalar>;
const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
auto dt = getParam<Scalar>("TimeLoop.DtInitial");
- // check if we are about to restart a previously interrupted simulation
- Scalar restartTime = 0;
- if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart"))
- restartTime = getParam<Scalar>("TimeLoop.Restart");
-
// instantiate time loop
- auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(restartTime, dt, tEnd);
+ auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(0, dt, tEnd);
timeLoop->setMaxTimeStepSize(maxDt);
stokesProblem->setTimeLoop(timeLoop);
darcyProblem->setTimeLoop(timeLoop);
@@ -157,32 +150,22 @@ int main(int argc, char** argv) try
sol[stokesFaceIdx].resize(stokesFvGridGeometry->numFaceDofs());
sol[darcyIdx].resize(darcyFvGridGeometry->numDofs());
- const auto& cellCenterSol = sol[stokesCellCenterIdx];
- const auto& faceSol = sol[stokesFaceIdx];
+ auto stokesSol = partial(sol, stokesCellCenterIdx, stokesFaceIdx);
- // apply initial solution for instationary problems
- typename GET_PROP_TYPE(StokesTypeTag, SolutionVector) stokesSol;
- std::get<0>(stokesSol) = cellCenterSol;
- std::get<1>(stokesSol) = faceSol;
stokesProblem->applyInitialSolution(stokesSol);
- auto solStokesOld = stokesSol;
- sol[stokesCellCenterIdx] = stokesSol[stokesCellCenterIdx];
- sol[stokesFaceIdx] = stokesSol[stokesFaceIdx];
-
darcyProblem->applyInitialSolution(sol[darcyIdx]);
- auto solDarcyOld = sol[darcyIdx];
auto solOld = sol;
couplingManager->init(stokesProblem, darcyProblem, sol);
// the grid variables
- using StokesGridVariables = typename GET_PROP_TYPE(StokesTypeTag, GridVariables);
+ using StokesGridVariables = GetPropType<StokesTypeTag, Properties::GridVariables>;
auto stokesGridVariables = std::make_shared<StokesGridVariables>(stokesProblem, stokesFvGridGeometry);
- stokesGridVariables->init(stokesSol, solStokesOld);
- using DarcyGridVariables = typename GET_PROP_TYPE(DarcyTypeTag, GridVariables);
+ stokesGridVariables->init(stokesSol);
+ using DarcyGridVariables = GetPropType<DarcyTypeTag, Properties::GridVariables>;
auto darcyGridVariables = std::make_shared<DarcyGridVariables>(darcyProblem, darcyFvGridGeometry);
- darcyGridVariables->init(sol[darcyIdx], solDarcyOld);
+ darcyGridVariables->init(sol[darcyIdx]);
// intialize the vtk output module
#if EXNUMBER >= 1
@@ -194,12 +177,14 @@ int main(int argc, char** argv) try
const auto darcyName = "orig_" + getParam<std::string>("Problem.Name") + "_" + darcyProblem->name();
#endif
- StaggeredVtkOutputModule<StokesTypeTag> stokesVtkWriter(*stokesProblem, *stokesFvGridGeometry, *stokesGridVariables, stokesSol, stokesName);
- GET_PROP_TYPE(StokesTypeTag, VtkOutputFields)::init(stokesVtkWriter);
- stokesVtkWriter.write(0.0);
+ StaggeredVtkOutputModule<StokesGridVariables, decltype(stokesSol)> stokesVtkWriter(*stokesGridVariables, stokesSol, stokesName);
+ GetPropType<StokesTypeTag, Properties::VtkOutputFields>::initOutputModule(stokesVtkWriter);
+ stokesVtkWriter.write(0);
- VtkOutputModule<DarcyTypeTag> darcyVtkWriter(*darcyProblem, *darcyFvGridGeometry, *darcyGridVariables, sol[darcyIdx], darcyName);
- GET_PROP_TYPE(DarcyTypeTag, VtkOutputFields)::init(darcyVtkWriter);
+ VtkOutputModule<DarcyGridVariables, GetPropType<DarcyTypeTag, Properties::SolutionVector>> darcyVtkWriter(*darcyGridVariables, sol[darcyIdx], darcyName);
+ using DarcyVelocityOutput = GetPropType<DarcyTypeTag, Properties::VelocityOutput>;
+ darcyVtkWriter.addVelocityOutput(std::make_shared<DarcyVelocityOutput>(*darcyGridVariables));
+ GetPropType<DarcyTypeTag, Properties::IOFields>::initOutputModule(darcyVtkWriter);
darcyVtkWriter.write(0.0);
// intialize the subproblems
@@ -221,14 +206,7 @@ int main(int argc, char** argv) try
using LinearSolver = UMFPackBackend;
auto linearSolver = std::make_shared<LinearSolver>();
- // the primary variable switches used by the sub models and the non-linear solver
-#if EXNUMBER >= 3
- using PriVarSwitchTuple = std::tuple<NoPrimaryVariableSwitch, NoPrimaryVariableSwitch, typename GET_PROP_TYPE(DarcyTypeTag, PrimaryVariableSwitch)>;
- using NewtonSolver = MultiDomainPriVarSwitchNewtonSolver<Assembler, LinearSolver, CouplingManager, PriVarSwitchTuple>;
-#else
-// using PriVarSwitchTuple = std::tuple<NoPrimaryVariableSwitch, NoPrimaryVariableSwitch, NoPrimaryVariableSwitch>;
using NewtonSolver = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
-#endif
NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
// time loop
diff --git a/exercises/solution/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.input b/exercises/solution/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.input
index 8d08fe83..1b0a62ad 100644
--- a/exercises/solution/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.input
+++ b/exercises/solution/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.input
@@ -14,10 +14,11 @@ Cells = 16 16
[Stokes.Problem]
Name = stokes
-EnableGravity = false
+EnableGravity = true
+EnableInertiaTerms = true
MoleFraction = 0.0 # -
Pressure = 1e5 # Pa
-Velocity = 1e-3 # m/s
+Velocity = 1 # m/s
[Darcy.Problem]
Name = darcy
diff --git a/exercises/solution/exercise-coupling-ff-pm/models/ex_models_ffproblem.hh b/exercises/solution/exercise-coupling-ff-pm/models/ex_models_ffproblem.hh
index 3d59b25a..ae63df0d 100644
--- a/exercises/solution/exercise-coupling-ff-pm/models/ex_models_ffproblem.hh
+++ b/exercises/solution/exercise-coupling-ff-pm/models/ex_models_ffproblem.hh
@@ -40,32 +40,41 @@ class StokesSubProblem;
namespace Properties
{
-NEW_TYPE_TAG(StokesTypeTag, INHERITS_FROM(StaggeredFreeFlowModel, NavierStokesNC));
+// Create new type tags
+namespace TTag {
+struct StokesTypeTag { using InheritsFrom = std::tuple<NavierStokesNC, StaggeredFreeFlowModel>; };
+} // end namespace TTag
// Set the grid type
-SET_TYPE_PROP(StokesTypeTag, Grid, Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<typename GET_PROP_TYPE(TypeTag, Scalar), 2> >);
+template<class TypeTag>
+struct Grid<TypeTag, TTag::StokesTypeTag> { using type = Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
// The fluid system
-SET_PROP(StokesTypeTag, FluidSystem)
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::StokesTypeTag>
{
- using H2OAir = FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>;
+ using H2OAir = FluidSystems::H2OAir<GetPropType<TypeTag, Properties::Scalar>>;
using type = FluidSystems::OnePAdapter<H2OAir, H2OAir::gasPhaseIdx>;
};
// Do not replace one equation with a total mass balance
-SET_INT_PROP(StokesTypeTag, ReplaceCompEqIdx, 3);
+template<class TypeTag>
+struct ReplaceCompEqIdx<TypeTag, TTag::StokesTypeTag> { static constexpr int value = 3; };
// Use formulation based on mass fractions
-SET_BOOL_PROP(StokesTypeTag, UseMoles, true);
+template<class TypeTag>
+struct UseMoles<TypeTag, TTag::StokesTypeTag> { static constexpr bool value = true; };
// Set the problem property
-SET_TYPE_PROP(StokesTypeTag, Problem, Dumux::StokesSubProblem<TypeTag> );
-
-SET_BOOL_PROP(StokesTypeTag, EnableFVGridGeometryCache, true);
-SET_BOOL_PROP(StokesTypeTag, EnableGridFluxVariablesCache, true);
-SET_BOOL_PROP(StokesTypeTag, EnableGridVolumeVariablesCache, true);
-
-SET_BOOL_PROP(StokesTypeTag, EnableInertiaTerms, false);
+template<class TypeTag>
+struct Problem<TypeTag, TTag::StokesTypeTag> { using type = Dumux::StokesSubProblem<TypeTag> ; };
+
+template<class TypeTag>
+struct EnableFVGridGeometryCache<TypeTag, TTag::StokesTypeTag> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridFluxVariablesCache<TypeTag, TTag::StokesTypeTag> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridVolumeVariablesCache<TypeTag, TTag::StokesTypeTag> { static constexpr bool value = true; };
}
/*!
@@ -79,31 +88,31 @@ class StokesSubProblem : public NavierStokesProblem<TypeTag>
{
using ParentType = NavierStokesProblem<TypeTag>;
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
- using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
- using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+ using GridView = GetPropType<TypeTag, Properties::GridView>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+ using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+ using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
- using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
using FVElementGeometry = typename FVGridGeometry::LocalView;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using Element = typename GridView::template Codim<0>::Entity;
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
- using ElementFaceVariables = typename GET_PROP_TYPE(TypeTag, GridFaceVariables)::LocalView;
- using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
+ using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
+ using ElementFaceVariables = typename GetPropType<TypeTag, Properties::GridFaceVariables>::LocalView;
+ using FluidState = GetPropType<TypeTag, Properties::FluidState>;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
- using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
- using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+ using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
- using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager);
+ using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
using TimeLoopPtr = std::shared_ptr<TimeLoop<Scalar>>;
- static constexpr bool useMoles = GET_PROP_TYPE(TypeTag, ModelTraits)::useMoles();
+ static constexpr bool useMoles = GetPropType<TypeTag, Properties::ModelTraits>::useMoles();
- static constexpr auto dim = GET_PROP_TYPE(TypeTag, ModelTraits)::dim();
+ static constexpr auto dim = GetPropType<TypeTag, Properties::ModelTraits>::dim();
static constexpr auto transportCompIdx = 1;
public:
@@ -120,10 +129,6 @@ public:
*/
// \{
-
- bool shouldWriteRestartFile() const
- { return false; }
-
/*!
* \brief Return the temperature within the domain in [K].
*/
@@ -221,9 +226,9 @@ public:
if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
{
- values[Indices::momentumYBalanceIdx] = couplingManager().couplingData().momentumCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf);
+ values[Indices::momentumYBalanceIdx] = couplingManager().couplingData().momentumCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf);
- const auto massFlux = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, elemFaceVars, scvf);
+ const auto massFlux = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, elemFaceVars, scvf);
values[Indices::conti0EqIdx] = massFlux[0];
values[Indices::conti0EqIdx + 1] = massFlux[1];
}
@@ -249,26 +254,34 @@ public:
*/
// \{
- /*!
- * \brief Evaluate the initial value for a control volume.
- *
- * \param globalPos The global position
- */
- PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
- {
- FluidState fluidState;
- updateFluidStateForBC_(fluidState, pressure_);
- const Scalar density = FluidSystem::density(fluidState, 0);
-
- PrimaryVariables values(0.0);
- values[Indices::pressureIdx] = pressure_ + density*this->gravity()[1]*(globalPos[1] - this->fvGridGeometry().bBoxMin()[1]);
- values[Indices::conti0EqIdx + 1] = moleFraction_;
- values[Indices::velocityXIdx] = 4.0 * velocity_ * (globalPos[1] - this->fvGridGeometry().bBoxMin()[1])
- * (this->fvGridGeometry().bBoxMax()[1] - globalPos[1])
- / (height_() * height_());
-
- return values;
- }
+ /*!
+ * \brief Evaluate the initial value for a control volume.
+ *
+ * \param globalPos The global position
+ */
+ PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
+ {
+ PrimaryVariables values(0.0);
+
+ // This is only an approximation of the actual hydorostatic pressure gradient.
+ // Air is compressible (the density depends on pressure), thus the actual
+ // vertical pressure profile is non-linear.
+ // This discrepancy can lead to spurious flows at the outlet if the inlet
+ // velocity is chosen too small.
+ FluidState fluidState;
+ updateFluidStateForBC_(fluidState, pressure_);
+ const Scalar density = FluidSystem::density(fluidState, 0);
+ values[Indices::pressureIdx] = pressure_ - density*this->gravity()[1]*(this->fvGridGeometry().bBoxMax()[1] - globalPos[1]);
+
+
+ // gravity has negative sign
+ values[Indices::conti0EqIdx + 1] = moleFraction_;
+ values[Indices::velocityXIdx] = 4.0 * velocity_ * (globalPos[1] - this->fvGridGeometry().bBoxMin()[1])
+ * (this->fvGridGeometry().bBoxMax()[1] - globalPos[1])
+ / (height_() * height_());
+
+ return values;
+ }
void setTimeLoop(TimeLoopPtr timeLoop)
{ timeLoop_ = timeLoop; }
@@ -276,9 +289,9 @@ public:
/*!
* \brief Returns the intrinsic permeability of required as input parameter for the Beavers-Joseph-Saffman boundary condition
*/
- Scalar permeability(const SubControlVolumeFace& scvf) const
+ Scalar permeability(const Element& element, const SubControlVolumeFace& scvf) const
{
- return couplingManager().couplingData().darcyPermeability(scvf);
+ return couplingManager().couplingData().darcyPermeability(element, scvf);
}
/*!
diff --git a/exercises/solution/exercise-coupling-ff-pm/models/ex_models_pmproblem.hh b/exercises/solution/exercise-coupling-ff-pm/models/ex_models_pmproblem.hh
index b2d951bc..cb79f6a3 100644
--- a/exercises/solution/exercise-coupling-ff-pm/models/ex_models_pmproblem.hh
+++ b/exercises/solution/exercise-coupling-ff-pm/models/ex_models_pmproblem.hh
@@ -49,19 +49,24 @@ class DarcySubProblem;
namespace Properties
{
+// Create new type tags
+namespace TTag {
#if EXNUMBER >= 1
-NEW_TYPE_TAG(DarcyTypeTag, INHERITS_FROM(CCTpfaModel, TwoPNC));
+struct DarcyTypeTag { using InheritsFrom = std::tuple<TwoPNC, CCTpfaModel>; };
#else
-NEW_TYPE_TAG(DarcyTypeTag, INHERITS_FROM(CCTpfaModel, OnePNC));
+struct DarcyTypeTag { using InheritsFrom = std::tuple<OnePNC, CCTpfaModel>; };
#endif
+} // end namespace TTag
// Set the problem property
-SET_TYPE_PROP(DarcyTypeTag, Problem, Dumux::DarcySubProblem<TypeTag>);
+template<class TypeTag>
+struct Problem<TypeTag, TTag::DarcyTypeTag> { using type = Dumux::DarcySubProblem<TypeTag>; };
// The fluid system
-SET_PROP(DarcyTypeTag, FluidSystem)
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::DarcyTypeTag>
{
- using H2OAir = FluidSystems::H2OAir<typename GET_PROP_TYPE(TypeTag, Scalar)>;
+ using H2OAir = FluidSystems::H2OAir<GetPropType<TypeTag, Properties::Scalar>>;
#if EXNUMBER == 0
using type = FluidSystems::OnePAdapter<H2OAir, H2OAir::gasPhaseIdx>;
#else
@@ -70,29 +75,35 @@ SET_PROP(DarcyTypeTag, FluidSystem)
};
// Use moles
-SET_BOOL_PROP(DarcyTypeTag, UseMoles, true);
+template<class TypeTag>
+struct UseMoles<TypeTag, TTag::DarcyTypeTag> { static constexpr bool value = true; };
// Do not replace one equation with a total mass balance
-SET_INT_PROP(DarcyTypeTag, ReplaceCompEqIdx, 3);
+template<class TypeTag>
+struct ReplaceCompEqIdx<TypeTag, TTag::DarcyTypeTag> { static constexpr int value = 3; };
//! Use a model with constant tortuosity for the effective diffusivity
SET_TYPE_PROP(DarcyTypeTag, EffectiveDiffusivityModel,
- DiffusivityConstantTortuosity<typename GET_PROP_TYPE(TypeTag, Scalar)>);
+ DiffusivityConstantTortuosity<GetPropType<TypeTag, Properties::Scalar>>);
// Set the grid type
-SET_TYPE_PROP(DarcyTypeTag, Grid, Dune::YaspGrid<2>);
+template<class TypeTag>
+struct Grid<TypeTag, TTag::DarcyTypeTag> { using type = Dune::YaspGrid<2>; };
#if EXNUMBER >= 1
//! Set the default formulation to pw-Sn: This can be over written in the problem.
-SET_PROP(DarcyTypeTag, Formulation)
+template<class TypeTag>
+struct Formulation<TypeTag, TTag::DarcyTypeTag>
{ static constexpr auto value = TwoPFormulation::p1s0; };
#endif
// Set the spatial paramaters type
#if EXNUMBER >= 1
-SET_TYPE_PROP(DarcyTypeTag, SpatialParams, TwoPSpatialParams<TypeTag>);
+template<class TypeTag>
+struct SpatialParams<TypeTag, TTag::DarcyTypeTag> { using type = TwoPSpatialParams<TypeTag>; };
#else
-SET_TYPE_PROP(DarcyTypeTag, SpatialParams, OnePSpatialParams<TypeTag>);
+template<class TypeTag>
+struct SpatialParams<TypeTag, TTag::DarcyTypeTag> { using type = OnePSpatialParams<TypeTag>; };
#endif
}
@@ -100,19 +111,19 @@ template <class TypeTag>
class DarcySubProblem : public PorousMediumFlowProblem<TypeTag>
{
using ParentType = PorousMediumFlowProblem<TypeTag>;
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
- using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
- using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
- using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
+ using GridView = GetPropType<TypeTag, Properties::GridView>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+ using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
+ using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
+ using FVElementGeometry = typename GetPropType<TypeTag, Properties::FVGridGeometry>::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
- using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
// copy some indices for convenience
- using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+ using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
enum {
// grid and world dimension
dim = GridView::dimension,
@@ -135,7 +146,7 @@ class DarcySubProblem : public PorousMediumFlowProblem<TypeTag>
using Element = typename GridView::template Codim<0>::Entity;
using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
- using CouplingManager = typename GET_PROP_TYPE(TypeTag, CouplingManager);
+ using CouplingManager = GetPropType<TypeTag, Properties::CouplingManager>;
using TimeLoopPtr = std::shared_ptr<TimeLoop<Scalar>>;
public:
@@ -148,7 +159,6 @@ public:
#else
moleFraction_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.MoleFraction");
#endif
- pressure_ = getParamFromGroup<Scalar>(this->paramGroup(), "Problem.Pressure");
// initialize output file
plotFluxes_ = getParamFromGroup<bool>(this->paramGroup(), "Problem.PlotFluxes", false);
@@ -275,10 +285,8 @@ public:
if (!couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
continue;
- // NOTE: binding the coupling context is necessary
- couplingManager_->bindCouplingContext(CouplingManager::darcyIdx, element);
#if EXNUMBER >= 2
- NumEqVector flux = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, scvf);
+ NumEqVector flux = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, scvf);
#else
NumEqVector flux(0.0); // add "massCouplingCondition" from the couplingManager here
#endif
@@ -301,21 +309,6 @@ public:
gnuplotInterfaceFluxes_.plot("flux_" + std::to_string(timeLoop_->timeStepIndex()));
}
- /*!
- * \brief Returns true if a restart file should be written to
- * disk.
- */
- bool shouldWriteRestartFile() const
- { return false; }
-
- /*!
- * \name Problem parameters
- */
- // \{
-
- bool shouldWriteOutput() const // define output
- { return true; }
-
/*!
* \brief Return the temperature within the domain in [K].
*
@@ -366,7 +359,7 @@ public:
NumEqVector values(0.0);
if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
- values = couplingManager().couplingData().massCouplingCondition(fvGeometry, elemVolVars, scvf);
+ values = couplingManager().couplingData().massCouplingCondition(element, fvGeometry, elemVolVars, scvf);
return values;
}
@@ -405,6 +398,8 @@ public:
*/
PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{
+ static const Scalar stokesPressure = getParamFromGroup<Scalar>("Stokes", "Problem.Pressure");
+
PrimaryVariables values(0.0);
#if EXNUMBER >= 3
values.setState(3/*bothPhases*/);
@@ -415,7 +410,7 @@ public:
#else
values[transportCompIdx] = moleFraction_;
#endif
- values[pressureIdx] = pressure_;
+ values[pressureIdx] = stokesPressure;
return values;
}
@@ -451,7 +446,6 @@ private:
#else
Scalar moleFraction_;
#endif
- Scalar pressure_;
Scalar initialWaterContent_ = 0.0;
Scalar lastWaterMass_ = 0.0;
--
GitLab