Skip to content
GitLab
Commit e53d6fab authored by Kilian Weishaupt's avatar Kilian Weishaupt
Browse files

[ff-pnm] Update second exercise

parent 440692db
Hide whitespace changes
Inline Side-by-side
exercises/exercise-coupling-ff-pm/README.md
View file @ e53d6fab
......@@ -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`).
......
exercises/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.cc
View file @ e53d6fab
......@@ -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);
......
exercises/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.input
View file @ e53d6fab
[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
exercises/exercise-coupling-ff-pm/models/ex_models_ffproblem.hh
View file @ e53d6fab
......@@ -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);
}
/*!
......
exercises/exercise-coupling-ff-pm/models/ex_models_pmproblem.hh
View file @ e53d6fab
......@@ -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;
......
exercises/solution/exercise-coupling-ff-pm/models/ex_models_coupling_ff-pm.cc
View file @ e53d6fab
......@@ -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>