From b69ade30fe984957606d8ca6544d35fd750f1421 Mon Sep 17 00:00:00 2001
From: Bernd Flemisch <bernd@iws.uni-stuttgart.de>
Date: Thu, 19 Jul 2018 03:16:24 +0200
Subject: [PATCH] update solution to properties exercise
---
.../exercise-properties/CMakeLists.txt | 4 +-
.../exercise_properties_solution.cc | 79 +++--
.../exercise_properties_solution.input | 28 +-
.../injection2p2cproblem.hh | 278 ------------------
.../injection2p2cspatialparams.hh | 190 ------------
.../exercise-properties/mylocalresidual.hh | 146 ++++-----
.../exercise-properties/mymateriallaw.hh | 114 -------
.../solution/exercise-properties/problem.hh | 251 ++++++++++++++++
.../exercise-properties/spatialparams.hh | 161 ++++++++++
9 files changed, 543 insertions(+), 708 deletions(-)
delete mode 100644 exercises/solution/exercise-properties/injection2p2cproblem.hh
delete mode 100644 exercises/solution/exercise-properties/injection2p2cspatialparams.hh
delete mode 100644 exercises/solution/exercise-properties/mymateriallaw.hh
create mode 100644 exercises/solution/exercise-properties/problem.hh
create mode 100644 exercises/solution/exercise-properties/spatialparams.hh
diff --git a/exercises/solution/exercise-properties/CMakeLists.txt b/exercises/solution/exercise-properties/CMakeLists.txt
index fec78d7b..0d8fe365 100644
--- a/exercises/solution/exercise-properties/CMakeLists.txt
+++ b/exercises/solution/exercise-properties/CMakeLists.txt
@@ -1,7 +1,7 @@
-# the compositional two-phase simulation program
+# the solution to the properties exercise
dune_add_test(NAME exercise_properties_solution
SOURCES exercise_properties_solution.cc
- CMD_ARGS exercise_properties.input -Problem.OnlyPlotMaterialLaws false)
+ CMD_ARGS exercise_properties_solution.input)
# add tutorial to the common target
add_dependencies(test_exercises exercise_properties_solution)
diff --git a/exercises/solution/exercise-properties/exercise_properties_solution.cc b/exercises/solution/exercise-properties/exercise_properties_solution.cc
index 2a7493e6..3583260e 100644
--- a/exercises/solution/exercise-properties/exercise_properties_solution.cc
+++ b/exercises/solution/exercise-properties/exercise_properties_solution.cc
@@ -19,7 +19,7 @@
/*!
* \file
*
- * \brief Test for the two-phase two-component CC model.
+ * \brief test for the two-phase porousmedium flow model
*/
#include <config.h>
@@ -30,14 +30,18 @@
#include <dune/common/timer.hh>
#include <dune/grid/io/file/dgfparser/dgfexception.hh>
#include <dune/grid/io/file/vtk.hh>
+#include <dune/istl/io.hh>
+
+#include "problem.hh"
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
+#include <dumux/common/valgrind.hh>
#include <dumux/common/dumuxmessage.hh>
#include <dumux/common/defaultusagemessage.hh>
#include <dumux/linear/amgbackend.hh>
-#include <dumux/nonlinear/privarswitchnewtonsolver.hh>
+#include <dumux/nonlinear/newtonsolver.hh>
#include <dumux/assembly/fvassembler.hh>
#include <dumux/assembly/diffmethod.hh>
@@ -47,14 +51,44 @@
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/io/grid/gridmanager.hh>
-#include "injection2p2cproblem.hh"
+/*!
+ * \brief Provides an interface for customizing error messages associated with
+ * reading in parameters.
+ *
+ * \param progName The name of the program, that was tried to be started.
+ * \param errorMsg The error message that was issued by the start function.
+ * Comprises the thing that went wrong and a general help message.
+ */
+void usage(const char *progName, const std::string &errorMsg)
+{
+ if (errorMsg.size() > 0) {
+ std::string errorMessageOut = "\nUsage: ";
+ errorMessageOut += progName;
+ errorMessageOut += " [options]\n";
+ errorMessageOut += errorMsg;
+ errorMessageOut += "\n\nThe list of mandatory arguments for this program is:\n"
+ "\t-TimeManager.TEnd End of the simulation [s] \n"
+ "\t-TimeManager.DtInitial Initial timestep size [s] \n"
+ "\t-Grid.LowerLeft Lower left corner coordinates\n"
+ "\t-Grid.UpperRight Upper right corner coordinates\n"
+ "\t-Grid.Cells Number of cells in respective coordinate directions\n"
+ "\t definition in DGF format\n"
+ "\t-SpatialParams.LensLowerLeft coordinates of the lower left corner of the lens [m] \n"
+ "\t-SpatialParams.LensUpperRight coordinates of the upper right corner of the lens [m] \n"
+ "\t-SpatialParams.Permeability Permeability of the domain [m^2] \n"
+ "\t-SpatialParams.PermeabilityLens Permeability of the lens [m^2] \n";
+
+ std::cout << errorMessageOut
+ << "\n";
+ }
+}
-int main(int argc, char** argv)try
+int main(int argc, char** argv) try
{
using namespace Dumux;
// define the type tag for this problem
- using TypeTag = TTAG(Injection2p2pcCCTypeTag);
+ using TypeTag = TTAG(TwoPIncompressibleTpfa);
// initialize MPI, finalize is done automatically on exit
const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
@@ -64,7 +98,7 @@ int main(int argc, char** argv)try
DumuxMessage::print(/*firstCall=*/true);
// parse command line arguments and input file
- Parameters::init(argc, argv);
+ Parameters::init(argc, argv, usage);
// try to create a grid (from the given grid file or the input file)
GridManager<typename GET_PROP_TYPE(TypeTag, Grid)> gridManager;
@@ -103,18 +137,19 @@ int main(int argc, char** argv)try
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");
-
// intialize the vtk output module
using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
- VtkOutputModule<TypeTag> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
- VtkOutputFields::init(vtkWriter); //! Add model specific output fields
+
+ // use non-conforming output for the test with interface solver
+ const auto ncOutput = getParam<bool>("Problem.UseNonConformingOutput", false);
+ VtkOutputModule<TypeTag> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name(), "",
+ (ncOutput ? Dune::VTK::nonconforming : Dune::VTK::conforming));
+
+ VtkOutputFields::init(vtkWriter); //!< Add model specific output fields
+ vtkWriter.write(0.0);
// instantiate time loop
- auto timeLoop = std::make_shared<TimeLoop<Scalar>>(restartTime, dt, tEnd);
+ auto timeLoop = std::make_shared<TimeLoop<Scalar>>(0.0, dt, tEnd);
timeLoop->setMaxTimeStepSize(maxDt);
// the assembler with time loop for instationary problem
@@ -126,8 +161,7 @@ int main(int argc, char** argv)try
auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->dofMapper());
// the non-linear solver
- using PrimaryVariableSwitch = typename GET_PROP_TYPE(TypeTag, PrimaryVariableSwitch);
- using NewtonSolver = Dumux::PriVarSwitchNewtonSolver<Assembler, LinearSolver, PrimaryVariableSwitch>;
+ using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
NewtonSolver nonLinearSolver(assembler, linearSolver);
// time loop
@@ -136,9 +170,6 @@ int main(int argc, char** argv)try
// set previous solution for storage evaluations
assembler->setPreviousSolution(xOld);
- //set time in problem (is used in time-dependent Neumann boundary condition)
- problem->setTime(timeLoop->time()+timeLoop->timeStepSize());
-
// solve the non-linear system with time step control
nonLinearSolver.solve(x, *timeLoop);
@@ -149,16 +180,20 @@ int main(int argc, char** argv)try
// advance to the time loop to the next step
timeLoop->advanceTimeStep();
+ // write vtk output
+ vtkWriter.write(timeLoop->time());
+
// report statistics of this time step
timeLoop->reportTimeStep();
- // set new dt as suggested by the newton solver
+ // set new dt as suggested by the Newton solver
timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
- vtkWriter.write(timeLoop->time());
-
} while (!timeLoop->finished());
+ // output some Newton statistics
+ nonLinearSolver.report();
+
timeLoop->finalize(leafGridView.comm());
////////////////////////////////////////////////////////////
diff --git a/exercises/solution/exercise-properties/exercise_properties_solution.input b/exercises/solution/exercise-properties/exercise_properties_solution.input
index df7737b8..7ae1761e 100644
--- a/exercises/solution/exercise-properties/exercise_properties_solution.input
+++ b/exercises/solution/exercise-properties/exercise_properties_solution.input
@@ -1,24 +1,18 @@
[TimeLoop]
-DtInitial = 3600 # in seconds
-TEnd = 3.154e9 # in seconds, i.e ten years
+DtInitial = 250 # [s]
+TEnd = 30000 # [s]
[Grid]
LowerLeft = 0 0
-UpperRight = 60 40
-Cells = 24 16
+UpperRight = 6 4
+Cells = 48 32
-[Problem]
-Name = infiltration
-OnlyPlotMaterialLaws = true
-AquiferDepth = 2700.0 # m
-TotalAreaSpecificInflow = 1e-4 # kg / (s*m^2)
-InjectionDuration = 2.628e6 # in seconds, i.e. one month
+[SpatialParams]
+LensLowerLeft = 1.0 2.0 # [m] coordinates of the lower left lens corner
+LensUpperRight = 4.0 3.0 # [m] coordinates of the upper right lens corner
-EnableGravity = true
-EnableDiffusion = true
+[Problem]
+Name = exercise_properties
-[SpatialParams]
-PermeabilityAquitard = 1e-15 # m^2
-EntryPressureAquitard = 4.5e4 # Pa
-PermeabilityAquifer = 1e-12 # m^2
-EntryPressureAquifer = 1e4 # Pa
+[Newton]
+MaxRelativeShift = 1e-5
diff --git a/exercises/solution/exercise-properties/injection2p2cproblem.hh b/exercises/solution/exercise-properties/injection2p2cproblem.hh
deleted file mode 100644
index e6512829..00000000
--- a/exercises/solution/exercise-properties/injection2p2cproblem.hh
+++ /dev/null
@@ -1,278 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * See the file COPYING for full copying permissions. *
- * *
- * This program is free software: you can redistribute it and/or modify *
- * it under the terms of the GNU General Public License as published by *
- * the Free Software Foundation, either version 2 of the License, or *
- * (at your option) any later version. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * *
- * You should have received a copy of the GNU General Public License *
- * along with this program. If not, see <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- *
- * \brief Problem where air is injected under a low permeable layer in a depth of 2700m.
- */
-#ifndef DUMUX_INJECTION_2P2C_PROBLEM_HH
-#define DUMUX_INJECTION_2P2C_PROBLEM_HH
-
-#include <dumux/porousmediumflow/2p2c/model.hh>
-#include <dumux/porousmediumflow/problem.hh>
-#include <dumux/material/fluidsystems/h2on2.hh>
-#include <dumux/discretization/cellcentered/tpfa/properties.hh>
-
-#include "injection2p2cspatialparams.hh"
-
-#include "mylocalresidual.hh"
-
-namespace Dumux {
-
-// foward declaration
-template <class TypeTag>
-class Injection2p2cProblem;
-
-// setup property TypeTag
-namespace Properties {
-NEW_TYPE_TAG(Injection2p2cTypeTag, INHERITS_FROM(TwoPTwoC));
-NEW_TYPE_TAG(Injection2p2pcCCTypeTag, INHERITS_FROM(CCTpfaModel, Injection2p2cTypeTag));
-
-// Set the grid type
-SET_TYPE_PROP(Injection2p2cTypeTag, Grid, Dune::YaspGrid<2>);
-
-// Set the problem property
-SET_TYPE_PROP(Injection2p2cTypeTag, Problem, Injection2p2cProblem<TypeTag>);
-
-// Set the spatial parameters
-SET_TYPE_PROP(Injection2p2cTypeTag, SpatialParams,
- InjectionSpatialParams<typename GET_PROP_TYPE(TypeTag, FVGridGeometry),
- typename GET_PROP_TYPE(TypeTag, Scalar)>);
-
-SET_TYPE_PROP(Injection2p2cTypeTag, LocalResidual, MyCompositionalLocalResidual<TypeTag>);
-
-// Set fluid configuration
-SET_TYPE_PROP(Injection2p2cTypeTag,
- FluidSystem,
- FluidSystems::H2ON2<typename GET_PROP_TYPE(TypeTag, Scalar), FluidSystems::H2ON2DefaultPolicy</*fastButSimplifiedRelations=*/ true>>);
-
-// Define whether mole(true) or mass (false) fractions are used
-SET_BOOL_PROP(Injection2p2cTypeTag, UseMoles, true);
-} // end namespace Properties
-
-/*!
- * \ingroup TwoPTwoCModel
- * \ingroup ImplicitTestProblems
- * \brief Problem where air is injected under a low permeable layer in a depth of 2700m.
- *
- * The domain is sized 60m times 40m and consists of two layers, a moderately
- * permeable one for \f$ y<22m\f$ and one with a lower permeablility
- * in the rest of the domain.
- *
- * Nitrogen is injected into a water-filled aquifer through a well. First, we inject for one month.
- * Then, we continue simulating the development of the nitrogen plume for 10 years.
- * This is realized with a Neumann boundary condition at the right boundary
- * (\f$ 7m<y<15m\f$). The aquifer is situated 2700m below sea level (the depth can be changed through the input file).
- * The injected fluid phase migrates upwards due to buoyancy.
- * It accumulates and partially enters the top layer lower permeable aquitard.
- *
- * The default setting for useMoles is true, i.e. each component is balaced in units of mole.
- * The property useMoles can be set to either true or false in the
- * problem file. If you change this, make sure that the according units are used in the problem setup.
- *
- * This problem uses the \ref TwoPTwoCModel.
- */
-template <class TypeTag>
-class Injection2p2cProblem : public PorousMediumFlowProblem<TypeTag>
-{
- using ParentType = PorousMediumFlowProblem<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 NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
- using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
- using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
- using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
-
- enum { dimWorld = GridView::dimensionworld };
- using Element = typename GridView::template Codim<0>::Entity;
- using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
-
-public:
- Injection2p2cProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
- : ParentType(fvGridGeometry)
- {
- // initialize the tables of the fluid system
- FluidSystem::init(/*tempMin=*/273.15,
- /*tempMax=*/423.15,
- /*numTemp=*/50,
- /*pMin=*/0.0,
- /*pMax=*/30e6,
- /*numP=*/300);
-
- // name of the problem and output file
- name_ = getParam<std::string>("Problem.Name");
- // depth of the aquifer, units: m
- aquiferDepth_ = getParam<Scalar>("Problem.AquiferDepth");
- // inflow rate of nitrogen water vapor mixture, units: kg/(s m^2)
- totalAreaSpecificInflow_ = getParam<Scalar>("Problem.TotalAreaSpecificInflow");
- // the duration of the injection, units: second
- injectionDuration_ = getParam<Scalar>("Problem.InjectionDuration");
- }
-
- /*!
- * \name Problem parameters
- */
- // \{
-
- /*!
- * \brief Returns the problem name
- *
- * This is used as a prefix for files generated by the simulation.
- */
- const std::string& name() const
- { return name_; }
-
- /*!
- * \brief Returns the temperature in \f$ K \f$
- */
- Scalar temperature() const
- { return 273.15 + 30; }
-
- // \}
-
- /*!
- * \name Boundary conditions
- */
- // \{
-
- /*!
- * \brief Specifies which kind of boundary condition should be
- * used for which equation on a given boundary segment.
- *
- * \param bcTypes The boundary types for the conservation equations
- * \param globalPos The position for which the bc type should be evaluated
- */
- BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
- {
- BoundaryTypes bcTypes;
- if (globalPos[0] < eps_)
- bcTypes.setAllDirichlet();
-
- // and Neuman boundary conditions everywhere else
- // note that we don't differentiate between Neumann and Robin boundary types
- else
- bcTypes.setAllNeumann();
-
- return bcTypes;
- }
-
- /*!
- * \brief Evaluates the boundary conditions for a Dirichlet
- * boundary segment
- *
- * \param globalPos The global position
- */
- PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
- {
- return initialAtPos(globalPos);
- }
-
- /*!
- * \brief Evaluate the boundary conditions for a neumann
- * boundary segment.
- *
- * \param values Stores the Neumann values for the conservation equations in
- * \f$ [ \textnormal{unit of conserved quantity} / (m^(dim-1) \cdot s )] \f$
- * \param globalPos The position of the integration point of the boundary segment.
- *
- * For this method, the \a values parameter stores the mass flux
- * in normal direction of each phase. Negative values mean influx.
- */
- NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
- {
- // initialize values to zero, i.e. no-flow Neumann boundary conditions
- NumEqVector values(0.0);
-
- //if we are inside the injection zone set inflow Neumann boundary conditions
- if (time_ < injectionDuration_
- && globalPos[1] < 15 + eps_ && globalPos[1] > 7 - eps_ && globalPos[0] > 0.9*this->fvGridGeometry().bBoxMax()[0])
- {
- // set the Neumann values for the Nitrogen component balance
- // convert from units kg/(s*m^2) to mole/(s*m^2)
- values[Indices::conti0EqIdx + FluidSystem::N2Idx] = -totalAreaSpecificInflow_/FluidSystem::molarMass(FluidSystem::N2Idx);
- values[Indices::conti0EqIdx + FluidSystem::H2OIdx] = 0.0;
- }
-
- return values;
- }
-
- // \}
-
- /*!
- * \name Volume terms
- */
- // \{
-
- /*!
- * \brief Evaluate the initial value for a control volume.
- *
- * \param globalPos The position for which the initial condition should be evaluated
- *
- * For this method, the \a values parameter stores primary
- * variables.
- */
- PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
- {
- PrimaryVariables values(0.0);
- values.setState(Indices::firstPhaseOnly);
-
- // get the water density at atmospheric conditions
- const Scalar densityW = FluidSystem::H2O::liquidDensity(temperature(), 1.0e5);
-
- // assume an intially hydrostatic liquid pressure profile
- // note: we subtract rho_w*g*h because g is defined negative
- const Scalar pw = 1.0e5 - densityW*this->gravity()[dimWorld-1]*(aquiferDepth_ - globalPos[dimWorld-1]);
-
- // initially we have some nitrogen dissolved
- // saturation mole fraction would be
- // moleFracLiquidN2 = (pw + pc + p_vap^sat)/henry;
- const Scalar moleFracLiquidN2 = pw*0.95/BinaryCoeff::H2O_N2::henry(temperature());
-
- // note that because we start with a single phase system the primary variables
- // are pl and x^w_N2. This will switch as soon after we start injecting to a two
- // phase system so the primary variables will be pl and Sn (non-wetting saturation).
- values[Indices::switchIdx] = moleFracLiquidN2;
- values[Indices::pressureIdx] = pw;
-
- return values;
- }
-
- // \}
-
- void setTime(Scalar time)
- {
- time_ = time;
- }
-
-private:
- static constexpr Scalar eps_ = 1e-6;
- std::string name_; //! Problem name
- Scalar aquiferDepth_; //! Depth of the aquifer in m
- Scalar totalAreaSpecificInflow_; //! Area specific inflow rate in mole/(s*m^2)
- Scalar time_;
- Scalar injectionDuration_; //! Duration of the injection in seconds
-};
-
-} // end namespace Dumux
-
-#endif
diff --git a/exercises/solution/exercise-properties/injection2p2cspatialparams.hh b/exercises/solution/exercise-properties/injection2p2cspatialparams.hh
deleted file mode 100644
index 4c89071a..00000000
--- a/exercises/solution/exercise-properties/injection2p2cspatialparams.hh
+++ /dev/null
@@ -1,190 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * See the file COPYING for full copying permissions. *
- * *
- * This program is free software: you can redistribute it and/or modify *
- * it under the terms of the GNU General Public License as published by *
- * the Free Software Foundation, either version 2 of the License, or *
- * (at your option) any later version. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * *
- * You should have received a copy of the GNU General Public License *
- * along with this program. If not, see <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- *
- * \brief Definition of the spatial parameters for the injection problem
- * which uses the isothermal two-phase two-component
- * fully implicit model.
- */
-
-#ifndef DUMUX_INJECTION_SPATIAL_PARAMS_HH
-#define DUMUX_INJECTION_SPATIAL_PARAMS_HH
-
-#include <dumux/material/spatialparams/fv.hh>
-#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
-#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
-// TODO: dumux-course-task
-// Inlcude your own material law
-#include "mymateriallaw.hh"
-
-#include <dumux/io/gnuplotinterface.hh>
-#include <dumux/io/plotmateriallaw.hh>
-
-namespace Dumux {
-
-/*!
- * \ingroup TwoPTwoCModel
- * \brief Definition of the spatial parameters for the injection problem
- * which uses the isothermal two-phase two-component
- * fully implicit model.
- */
-template<class FVGridGeometry, class Scalar>
-class InjectionSpatialParams
-: public FVSpatialParams<FVGridGeometry, Scalar, InjectionSpatialParams<FVGridGeometry, Scalar>>
-{
- using ThisType = InjectionSpatialParams<FVGridGeometry, Scalar>;
- using ParentType = FVSpatialParams<FVGridGeometry, Scalar, ThisType>;
- using GridView = typename FVGridGeometry::GridView;
- static const int dimWorld = GridView::dimensionworld;
- using Element = typename GridView::template Codim<0>::Entity;
- using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
-
-public:
- // export permeability type
- using PermeabilityType = Scalar;
-
- // TODO: dumux-course-task
- // Use your own material law instead
- // Set the material law parameterized by absolute saturations
- using MaterialLaw = EffToAbsLaw<MyMaterialLaw<Scalar>>;
- using MaterialLawParams = typename MaterialLaw::Params;
-
- /*!
- * \brief The constructor
- *
- * \param fvGridGeometry The finite volume grid geometry
- */
- InjectionSpatialParams(std::shared_ptr<const FVGridGeometry>& fvGridGeometry)
- : ParentType(fvGridGeometry)
- {
- aquiferHeightFromBottom_ = 30.0;
-
- // intrinsic permeabilities
- aquitardK_ = getParam<Scalar>("SpatialParams.PermeabilityAquitard");
- aquiferK_ = getParam<Scalar>("SpatialParams.PermeabilityAquifer");
-
- // porosities
- aquitardPorosity_ = 0.2;
- aquiferPorosity_ = 0.4;
-
- // residual saturations
- aquitardMaterialParams_.setSwr(0.2);
- aquitardMaterialParams_.setSnr(0.0);
- aquiferMaterialParams_.setSwr(0.2);
- aquiferMaterialParams_.setSnr(0.0);
-
- // parameters for the Brooks-Corey law
- aquitardMaterialParams_.setPe(getParam<Scalar>("SpatialParams.EntryPressureAquitard"));
- aquiferMaterialParams_.setPe(getParam<Scalar>("SpatialParams.EntryPressureAquifer"));
- aquitardMaterialParams_.setLambda(2.0);
- aquiferMaterialParams_.setLambda(2.0);
-
- // plot the material laws using gnuplot and exit
- if (getParam<bool>("Problem.OnlyPlotMaterialLaws"))
- {
- plotMaterialLaws();
- exit(0);
- }
- }
-
- /*!
- * \brief Define the intrinsic permeability \f$\mathrm{[m^2]}\f$.
- *
- * \param globalPos The global position
- */
- PermeabilityType permeabilityAtPos(const GlobalPosition& globalPos) const
-
- {
- if (isInAquitard_(globalPos))
- return aquitardK_;
- return aquiferK_;
- }
-
- /*!
- * \brief Define the porosity \f$\mathrm{[-]}\f$.
- *
- * \param globalPos The global position
- */
- Scalar porosityAtPos(const GlobalPosition& globalPos) const
- {
- if (isInAquitard_(globalPos))
- return aquitardPorosity_;
- return aquiferPorosity_;
- }
-
- /*!
- * \brief Function for defining the parameters needed by constitutive relationships (kr-sw, pc-sw, etc.).
- *
- * \param globalPos The global position
- *
- * \return the material parameters object
- */
- const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const
- {
- if (isInAquitard_(globalPos))
- return aquitardMaterialParams_;
- return aquiferMaterialParams_;
- }
-
- /*!
- * \brief Function for defining which phase is to be considered as the wetting phase.
- *
- * \return the wetting phase index
- * \param globalPos The position of the center of the element
- */
- template<class FluidSystem>
- int wettingPhaseAtPos(const GlobalPosition& globalPos) const
- { return FluidSystem::H2OIdx; }
-
- /*!
- * \brief Creates a gnuplot output of the pc-Sw curve
- */
- void plotMaterialLaws()
- {
- PlotMaterialLaw<Scalar, MaterialLaw> plotMaterialLaw;
- GnuplotInterface<Scalar> gnuplot;
- plotMaterialLaw.addpcswcurve(gnuplot, aquitardMaterialParams_, 0.2, 1.0, "upper layer (fine, aquitard)", "w lp");
- plotMaterialLaw.addpcswcurve(gnuplot, aquiferMaterialParams_, 0.2, 1.0, "lower layer (coarse, aquifer)", "w l");
- gnuplot.setOption("set xrange [0:1]");
- gnuplot.setOption("set label \"residual\\nsaturation\" at 0.1,100000 center");
- gnuplot.plot("pc-Sw");
- }
-
-private:
-
- static constexpr Scalar eps_ = 1e-6;
-
- bool isInAquitard_(const GlobalPosition &globalPos) const
- { return globalPos[dimWorld-1] > aquiferHeightFromBottom_ + eps_; }
-
- Scalar aquitardK_;
- Scalar aquiferK_;
- Scalar aquiferHeightFromBottom_;
-
- Scalar aquitardPorosity_;
- Scalar aquiferPorosity_;
-
- MaterialLawParams aquitardMaterialParams_;
- MaterialLawParams aquiferMaterialParams_;
-};
-
-} // end namespace Dumux
-
-#endif
diff --git a/exercises/solution/exercise-properties/mylocalresidual.hh b/exercises/solution/exercise-properties/mylocalresidual.hh
index ac668555..485f33a8 100644
--- a/exercises/solution/exercise-properties/mylocalresidual.hh
+++ b/exercises/solution/exercise-properties/mylocalresidual.hh
@@ -18,127 +18,110 @@
*****************************************************************************/
/*!
* \file
- *
- * \brief Element-wise calculation of the local residual for problems
- * using compositional fully implicit model.
+ * \ingroup PorousmediumImmiscible
+ * \brief Element-wise calculation of the residual for problems
+ * using the n-phase immiscible fully implicit models.
*/
-#ifndef DUMUX_MY_COMPOSITIONAL_LOCAL_RESIDUAL_HH
-#define DUMUX_MY_COMPOSITIONAL_LOCAL_RESIDUAL_HH
+#ifndef DUMUX_MY_LOCAL_RESIDUAL_HH
+#define DUMUX_MY_LOCAL_RESIDUAL_HH
#include <dumux/common/properties.hh>
namespace Dumux
{
/*!
- * \ingroup Implicit
- * \ingroup ImplicitLocalResidual
- * \brief Element-wise calculation of the local residual for problems
- * using compositional fully implicit model.
- *
+ * \ingroup PorousmediumImmiscible
+ * \brief Element-wise calculation of the residual for problems
+ * using the n-phase immiscible fully implicit models.
*/
template<class TypeTag>
-class MyCompositionalLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
+class MyLocalResidual : public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
{
using ParentType = typename GET_PROP_TYPE(TypeTag, BaseLocalResidual);
- using Implementation = typename GET_PROP_TYPE(TypeTag, LocalResidual);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
+ using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
+ using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, GridFluxVariablesCache)::LocalView;
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
- using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
- using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
- using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, GridFluxVariablesCache)::LocalView;
- using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Element = typename GridView::template Codim<0>::Entity;
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
- using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using EnergyLocalResidual = typename GET_PROP_TYPE(TypeTag, EnergyLocalResidual);
- static constexpr int numPhases = GET_PROP_TYPE(TypeTag, ModelTraits)::numPhases();
- static constexpr int numComponents = GET_PROP_TYPE(TypeTag, ModelTraits)::numComponents();
-
- enum { conti0EqIdx = Indices::conti0EqIdx };
+ using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
+ static constexpr int numPhases = ModelTraits::numPhases();
+ static constexpr int conti0EqIdx = ModelTraits::Indices::conti0EqIdx; //!< first index for the mass balance
public:
using ParentType::ParentType;
/*!
- * \brief Evaluate the amount of all conservation quantities
- * (e.g. phase mass) within a sub-control volume.
- *
- * The result should be averaged over the volume (e.g. phase mass
- * inside a sub control volume divided by the volume)
- *
- * \param storage The mass of the component within the sub-control volume
- * \param scvIdx The SCV (sub-control-volume) index
- * \param usePrevSol Evaluate function with solution of current or previous time step
+ * \brief Evaluate the rate of change of all conservation
+ * quantites (e.g. phase mass) within a sub-control
+ * volume of a finite volume element for the immiscible models.
+ * \param problem The problem
+ * \param scv The sub control volume
+ * \param volVars The current or previous volVars
+ * \note This function should not include the source and sink terms.
+ * \note The volVars can be different to allow computing
+ * the implicit euler time derivative here
*/
- ResidualVector computeStorage(const Problem& problem,
- const SubControlVolume& scv,
- const VolumeVariables& volVars) const
+ NumEqVector computeStorage(const Problem& problem,
+ const SubControlVolume& scv,
+ const VolumeVariables& volVars) const
{
- ResidualVector storage(0.0);
-
- // compute storage term of all components within all phases
+ // partial time derivative of the phase mass
+ NumEqVector storage;
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
- for (int compIdx = 0; compIdx < numComponents; ++compIdx)
- {
- auto eqIdx = conti0EqIdx + compIdx;
- storage[eqIdx] += volVars.porosity()
- * volVars.saturation(phaseIdx)
- * volVars.molarDensity(phaseIdx)
- * volVars.moleFraction(phaseIdx, compIdx);
- }
+ auto eqIdx = conti0EqIdx + phaseIdx;
+ storage[eqIdx] = volVars.porosity()
+ * volVars.saturation(phaseIdx);
+
+ //! The energy storage in the fluid phase with index phaseIdx
+ EnergyLocalResidual::fluidPhaseStorage(storage, scv, volVars, phaseIdx);
}
+ //! The energy storage in the solid matrix
+ EnergyLocalResidual::solidPhaseStorage(storage, scv, volVars);
+
return storage;
}
+
/*!
- * \brief Evaluates the total flux of all conservation quantities
- * over a face of a sub-control volume.
+ * \brief Evaluate the mass flux over a face of a sub control volume
*
- * \param flux The flux over the SCV (sub-control-volume) face for each component
- * \param fIdx The index of the SCV face
- * \param onBoundary A boolean variable to specify whether the flux variables
- * are calculated for interior SCV faces or boundary faces, default=false
+ * \param problem The problem
+ * \param element The element
+ * \param fvGeometry The finite volume geometry context
+ * \param elemVolVars The volume variables for all flux stencil elements
+ * \param scvf The sub control volume face to compute the flux on
+ * \param elemFluxVarsCache The cache related to flux compuation
*/
- ResidualVector computeFlux(const Problem& problem,
- const Element& element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& elemVolVars,
- const SubControlVolumeFace& scvf,
- const ElementFluxVariablesCache& elemFluxVarsCache) const
+ NumEqVector computeFlux(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf,
+ const ElementFluxVariablesCache& elemFluxVarsCache) const
{
FluxVariables fluxVars;
fluxVars.init(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
- // get upwind weights into local scope
- ResidualVector flux(0.0);
- auto enableDiffusion = getParam<bool>("Problem.EnableDiffusion", true);
-
- // advective fluxes
+ NumEqVector flux;
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
- const auto diffusiveFluxes = fluxVars.molecularDiffusionFlux(phaseIdx);
- for (int compIdx = 0; compIdx < numComponents; ++compIdx)
- {
- // get equation index
- const auto eqIdx = conti0EqIdx + compIdx;
-
- // the physical quantities for which we perform upwinding
- const auto upwindTerm = [phaseIdx,compIdx] (const auto& volVars)
- { return volVars.molarDensity(phaseIdx)*volVars.moleFraction(phaseIdx, compIdx)*volVars.mobility(phaseIdx); };
- flux[eqIdx] += fluxVars.advectiveFlux(phaseIdx, upwindTerm);
-
- // TODO: dumux-course-task
- // same here: only add diffusive fluxes if diffusion is enabled
- if (enableDiffusion)
- flux[eqIdx] += diffusiveFluxes[compIdx];
- }
+ // the physical quantities for which we perform upwinding
+ auto upwindTerm = [phaseIdx](const auto& volVars)
+ { return volVars.mobility(phaseIdx); };
+
+ auto eqIdx = conti0EqIdx + phaseIdx;
+ flux[eqIdx] = fluxVars.advectiveFlux(phaseIdx, upwindTerm);
//! Add advective phase energy fluxes. For isothermal model the contribution is zero.
EnergyLocalResidual::heatConvectionFlux(flux, fluxVars, phaseIdx);
@@ -149,13 +132,6 @@ public:
return flux;
}
-
-protected:
- Implementation *asImp_()
- { return static_cast<Implementation *> (this); }
-
- const Implementation *asImp_() const
- { return static_cast<const Implementation *> (this); }
};
} // end namespace Dumux
diff --git a/exercises/solution/exercise-properties/mymateriallaw.hh b/exercises/solution/exercise-properties/mymateriallaw.hh
deleted file mode 100644
index fe36bec6..00000000
--- a/exercises/solution/exercise-properties/mymateriallaw.hh
+++ /dev/null
@@ -1,114 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * See the file COPYING for full copying permissions. *
- * *
- * This program is free software: you can redistribute it and/or modify *
- * it under the terms of the GNU General Public License as published by *
- * the Free Software Foundation, either version 2 of the License, or *
- * (at your option) any later version. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * *
- * You should have received a copy of the GNU General Public License *
- * along with this program. If not, see <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- *
- * \brief Implementation of the capillary pressure and
- * relative permeability <-> saturation relations.
- *
- */
-#ifndef DUMUX_MY_MATERIAL_LAW_HH
-#define DUMUX_MY_MATERIAL_LAW_HH
-
-#include <dumux/material/fluidmatrixinteractions/2p/brookscoreyparams.hh>
-#include <cmath>
-#include <algorithm>
-
-namespace Dumux
-{
-/*!
- * \ingroup fluidmatrixinteractionslaws
- * \note a simple material law using the BrooksCoreyParams
- */
-template <class ScalarT, class ParamsT = BrooksCoreyParams<ScalarT> >
-class MyMaterialLaw
-{
-public:
- typedef ParamsT Params;
- typedef typename Params::Scalar Scalar;
-
- /*!
- * \brief The capillary pressure-saturation curve
- * \param swe saturation of the wetting phase \f$\mathrm{[\overline{S}_w]}\f$
- * \param params A container object that is populated with the appropriate coefficients for the respective law.
- * Therefore, in the (problem specific) spatialParameters first, the material law is chosen,
- and then the params container is constructed accordingly. Afterwards the values are set there, too.
- * \return capillary pressure
- */
- static Scalar pc(const Params ¶ms, Scalar swe)
- {
- return 1.0e5*(1.0-swe) + params.pe();
- }
-
- /*!
- * \brief The relative permeability for the wetting phase of
- * the medium implied by the Brooks-Corey
- * parameterization.
- *
- * \param swe The mobile saturation of the wetting phase.
- * \param params A container object that is populated with the appropriate coefficients for the respective law.
- * Therefore, in the (problem specific) spatialParameters first, the material law is chosen,
- * and then the params container is constructed accordingly. Afterwards the values are set there, too.
- * \return Relative permeability of the wetting phase calculated as implied by Brooks & Corey.
- *
- * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number,
- * by clamping the input.
- */
- static Scalar krw(const Params ¶ms, Scalar swe)
- {
- using std::pow;
- using std::min;
- using std::max;
-
- swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= swe <= 1.0
-
- return pow(swe, 2.0/params.lambda() + 3);
- }
-
- /*!
- * \brief The relative permeability for the non-wetting phase of
- * the medium as implied by the Brooks-Corey
- * parameterization.
- *
- * \param swe The mobile saturation of the wetting phase.
- * \param params A container object that is populated with the appropriate coefficients for the respective law.
- * Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
- * is constructed accordingly. Afterwards the values are set there, too.
- * \return Relative permeability of the non-wetting phase calculated as implied by Brooks & Corey.
- *
- * \note Instead of undefined behaviour if pc is not in the valid range, we return a valid number,
- * by clamping the input.
- */
- static Scalar krn(const Params ¶ms, Scalar swe)
- {
- using std::pow;
- using std::min;
- using std::max;
-
- swe = min(max(swe, 0.0), 1.0); // the equation below is only defined for 0.0 <= swe <= 1.0
-
- const Scalar exponent = 2.0/params.lambda() + 1;
- const Scalar tmp = 1.0 - swe;
- return tmp*tmp*(1.0 - pow(swe, exponent));
- }
-};
-
-} // end namespace Dumux
-
-#endif
diff --git a/exercises/solution/exercise-properties/problem.hh b/exercises/solution/exercise-properties/problem.hh
new file mode 100644
index 00000000..4e1a3b40
--- /dev/null
+++ b/exercises/solution/exercise-properties/problem.hh
@@ -0,0 +1,251 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \ingroup TwoPTests
+ * \brief The properties for the incompressible 2p test
+ */
+#ifndef DUMUX_INCOMPRESSIBLE_TWOP_TEST_PROBLEM_HH
+#define DUMUX_INCOMPRESSIBLE_TWOP_TEST_PROBLEM_HH
+
+#include <dune/grid/yaspgrid.hh>
+
+#include <dumux/discretization/box/properties.hh>
+#include <dumux/discretization/cellcentered/tpfa/properties.hh>
+#include <dumux/discretization/cellcentered/mpfa/properties.hh>
+
+#include <dumux/material/components/trichloroethene.hh>
+#include <dumux/material/components/simpleh2o.hh>
+#include <dumux/material/fluidsystems/1pliquid.hh>
+#include <dumux/material/fluidsystems/2pimmiscible.hh>
+
+#include <dumux/porousmediumflow/problem.hh>
+#include <dumux/porousmediumflow/2p/model.hh>
+#include <dumux/porousmediumflow/2p/incompressiblelocalresidual.hh>
+
+#include "spatialparams.hh"
+#include "mylocalresidual.hh"
+
+namespace Dumux {
+// forward declarations
+template<class TypeTag> class TwoPTestProblem;
+
+namespace Properties {
+NEW_TYPE_TAG(TwoPIncompressible, INHERITS_FROM(TwoP));
+NEW_TYPE_TAG(TwoPIncompressibleTpfa, INHERITS_FROM(CCTpfaModel, TwoPIncompressible));
+
+// Set the grid type
+SET_TYPE_PROP(TwoPIncompressible, Grid, Dune::YaspGrid<2>);
+
+// Set the problem type
+SET_TYPE_PROP(TwoPIncompressible, Problem, TwoPTestProblem<TypeTag>);
+
+SET_TYPE_PROP(TwoPIncompressible, LocalResidual, MyLocalResidual<TypeTag>);
+
+// Set the fluid system
+SET_PROP(TwoPIncompressible, FluidSystem)
+{
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using WettingPhase = FluidSystems::OnePLiquid<Scalar, Components::SimpleH2O<Scalar> >;
+ using NonwettingPhase = FluidSystems::OnePLiquid<Scalar, Components::Trichloroethene<Scalar> >;
+ using type = FluidSystems::TwoPImmiscible<Scalar, WettingPhase, NonwettingPhase>;
+};
+
+// Set the spatial parameters
+SET_PROP(TwoPIncompressible, SpatialParams)
+{
+private:
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+public:
+ using type = TwoPTestSpatialParams<FVGridGeometry, Scalar>;
+};
+} // end namespace Properties
+
+/*!
+ * \ingroup TwoPTests
+ * \brief The incompressible 2p test problem.
+ */
+template<class TypeTag>
+class TwoPTestProblem : public PorousMediumFlowProblem<TypeTag>
+{
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+ using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+ using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
+ enum {
+ pressureH2OIdx = Indices::pressureIdx,
+ saturationDNAPLIdx = Indices::saturationIdx,
+ contiDNAPLEqIdx = Indices::conti0EqIdx + FluidSystem::comp1Idx,
+ waterPhaseIdx = FluidSystem::phase0Idx,
+ dnaplPhaseIdx = FluidSystem::phase1Idx
+ };
+
+public:
+ TwoPTestProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+ : ParentType(fvGridGeometry) {}
+
+ /*!
+ * \brief Specifies which kind of boundary condition should be
+ * used for which equation on a given boundary segment
+ *
+ * \param values Stores the value of the boundary type
+ * \param globalPos The global position
+ */
+ BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
+ {
+ BoundaryTypes values;
+ if (onLeftBoundary_(globalPos) || onRightBoundary_(globalPos))
+ values.setAllDirichlet();
+ else
+ values.setAllNeumann();
+ return values;
+ }
+
+ /*!
+ * \brief Evaluates the boundary conditions for a Dirichlet
+ * boundary segment
+ *
+ * \param values Stores the Dirichlet values for the conservation equations in
+ * \f$ [ \textnormal{unit of primary variable} ] \f$
+ * \param globalPos The global position
+ */
+ PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
+ {
+ PrimaryVariables values;
+ typename GET_PROP_TYPE(TypeTag, FluidState) fluidState;
+ fluidState.setTemperature(temperature());
+ fluidState.setPressure(waterPhaseIdx, /*pressure=*/1e5);
+ fluidState.setPressure(dnaplPhaseIdx, /*pressure=*/1e5);
+
+ Scalar densityW = FluidSystem::density(fluidState, waterPhaseIdx);
+
+ Scalar height = this->fvGridGeometry().bBoxMax()[1] - this->fvGridGeometry().bBoxMin()[1];
+ Scalar depth = this->fvGridGeometry().bBoxMax()[1] - globalPos[1];
+ Scalar alpha = 1 + 1.5/height;
+ Scalar width = this->fvGridGeometry().bBoxMax()[0] - this->fvGridGeometry().bBoxMin()[0];
+ Scalar factor = (width*alpha + (1.0 - alpha)*globalPos[0])/width;
+
+ // hydrostatic pressure scaled by alpha
+ values[pressureH2OIdx] = 1e5 - factor*densityW*this->gravity()[1]*depth;
+ values[saturationDNAPLIdx] = 0.0;
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a neumann
+ * boundary segment.
+ *
+ * \param values Stores the Neumann values for the conservation equations in
+ * \f$ [ \textnormal{unit of conserved quantity} / (m^(dim-1) \cdot s )] \f$
+ * \param globalPos The position of the integration point of the boundary segment.
+ *
+ * For this method, the \a values parameter stores the mass flux
+ * in normal direction of each phase. Negative values mean influx.
+ */
+ NumEqVector neumannAtPos(const GlobalPosition &globalPos) const
+ {
+ NumEqVector values(0.0);
+ if (onInlet_(globalPos))
+ {
+ using TCE = Components::Trichloroethene<Scalar>;
+ values[contiDNAPLEqIdx] = -0.04/TCE::liquidDensity(temperature(), /*pressure=*/1e5);// m/s
+ }
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluates the initial values for a control volume
+ *
+ * \param values Stores the initial values for the conservation equations in
+ * \f$ [ \textnormal{unit of primary variables} ] \f$
+ * \param globalPos The global position
+ */
+ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+ {
+ PrimaryVariables values;
+ typename GET_PROP_TYPE(TypeTag, FluidState) fluidState;
+ fluidState.setTemperature(temperature());
+ fluidState.setPressure(waterPhaseIdx, /*pressure=*/1e5);
+ fluidState.setPressure(dnaplPhaseIdx, /*pressure=*/1e5);
+
+ Scalar densityW = FluidSystem::density(fluidState, waterPhaseIdx);
+
+ Scalar depth = this->fvGridGeometry().bBoxMax()[1] - globalPos[1];
+
+ // hydrostatic pressure
+ values[pressureH2OIdx] = 1e5 - densityW*this->gravity()[1]*depth;
+ values[saturationDNAPLIdx] = 0;
+ return values;
+ }
+
+ /*!
+ * \brief Returns the temperature \f$\mathrm{[K]}\f$ for an isothermal problem.
+ *
+ * This is not specific to the discretization. By default it just
+ * throws an exception so it must be overloaded by the problem if
+ * no energy equation is used.
+ */
+ Scalar temperature() const
+ {
+ return 293.15; // 10°C
+ }
+
+private:
+ bool onLeftBoundary_(const GlobalPosition &globalPos) const
+ {
+ return globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_;
+ }
+
+ bool onRightBoundary_(const GlobalPosition &globalPos) const
+ {
+ return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_;
+ }
+
+ bool onLowerBoundary_(const GlobalPosition &globalPos) const
+ {
+ return globalPos[1] < this->fvGridGeometry().bBoxMin()[1] + eps_;
+ }
+
+ bool onUpperBoundary_(const GlobalPosition &globalPos) const
+ {
+ return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_;
+ }
+
+ bool onInlet_(const GlobalPosition &globalPos) const
+ {
+ Scalar width = this->fvGridGeometry().bBoxMax()[0] - this->fvGridGeometry().bBoxMin()[0];
+ Scalar lambda = (this->fvGridGeometry().bBoxMax()[0] - globalPos[0])/width;
+ return onUpperBoundary_(globalPos) && 0.5 < lambda && lambda < 2.0/3.0;
+ }
+
+ static constexpr Scalar eps_ = 1e-6;
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/exercises/solution/exercise-properties/spatialparams.hh b/exercises/solution/exercise-properties/spatialparams.hh
new file mode 100644
index 00000000..275984d8
--- /dev/null
+++ b/exercises/solution/exercise-properties/spatialparams.hh
@@ -0,0 +1,161 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \ingroup TwoPTests
+ * \brief The spatial params for the incompressible 2p test
+ */
+#ifndef DUMUX_INCOMPRESSIBLE_TWOP_TEST_SPATIAL_PARAMS_HH
+#define DUMUX_INCOMPRESSIBLE_TWOP_TEST_SPATIAL_PARAMS_HH
+
+#include <dumux/material/spatialparams/fv.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup TwoPTests
+ * \brief The spatial params for the incompressible 2p test
+ */
+template<class FVGridGeometry, class Scalar>
+class TwoPTestSpatialParams
+: public FVSpatialParams<FVGridGeometry, Scalar, TwoPTestSpatialParams<FVGridGeometry, Scalar>>
+{
+ using GridView = typename FVGridGeometry::GridView;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using ThisType = TwoPTestSpatialParams<FVGridGeometry, Scalar>;
+ using ParentType = FVSpatialParams<FVGridGeometry, Scalar, ThisType>;
+
+ static constexpr int dimWorld = GridView::dimensionworld;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ using EffectiveLaw = RegularizedVanGenuchten<Scalar>;
+
+public:
+ using MaterialLaw = EffToAbsLaw<EffectiveLaw>;
+ using MaterialLawParams = typename MaterialLaw::Params;
+ using PermeabilityType = Scalar;
+
+ TwoPTestSpatialParams(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+ : ParentType(fvGridGeometry)
+ {
+ lensLowerLeft_ = getParam<GlobalPosition>("SpatialParams.LensLowerLeft");
+ lensUpperRight_ = getParam<GlobalPosition>("SpatialParams.LensUpperRight");
+
+ // residual saturations
+ lensMaterialParams_.setSwr(0.18);
+ lensMaterialParams_.setSnr(0.0);
+ outerMaterialParams_.setSwr(0.05);
+ outerMaterialParams_.setSnr(0.0);
+
+ // parameters for the Van Genuchten law
+ // alpha and n
+ lensMaterialParams_.setVgAlpha(0.00045);
+ lensMaterialParams_.setVgn(7.3);
+ outerMaterialParams_.setVgAlpha(0.0037);
+ outerMaterialParams_.setVgn(4.7);
+
+ lensK_ = 9.05e-12;
+ outerK_ = 4.6e-10;
+ }
+
+ /*!
+ * \brief Function for defining the (intrinsic) permeability \f$[m^2]\f$.
+ * In this test, we use element-wise distributed permeabilities.
+ *
+ * \param element The current element
+ * \param scv The sub-control volume inside the element.
+ * \param elemSol The solution at the dofs connected to the element.
+ * \return permeability
+ */
+ template<class ElementSolution>
+ PermeabilityType permeability(const Element& element,
+ const SubControlVolume& scv,
+ const ElementSolution& elemSol) const
+ {
+ if (isInLens_(element.geometry().center()))
+ return lensK_;
+ return outerK_;
+ }
+
+ /*!
+ * \brief Returns the porosity \f$[-]\f$
+ *
+ * \param globalPos The global position
+ */
+ Scalar porosityAtPos(const GlobalPosition& globalPos) const
+ { return 0.4; }
+
+ /*!
+ * \brief Returns the parameter object for the Brooks-Corey material law.
+ * In this test, we use element-wise distributed material parameters.
+ *
+ * \param element The current element
+ * \param scv The sub-control volume inside the element.
+ * \param elemSol The solution at the dofs connected to the element.
+ * \return the material parameters object
+ */
+ template<class ElementSolution>
+ const MaterialLawParams& materialLawParams(const Element& element,
+ const SubControlVolume& scv,
+ const ElementSolution& elemSol) const
+ {
+ if (isInLens_(element.geometry().center()))
+ return lensMaterialParams_;
+ return outerMaterialParams_;
+ }
+
+ /*!
+ * \brief Function for defining which phase is to be considered as the wetting phase.
+ *
+ * \return the wetting phase index
+ * \param globalPos The global position
+ */
+ template<class FluidSystem>
+ int wettingPhaseAtPos(const GlobalPosition& globalPos) const
+ {
+ return FluidSystem::phase0Idx;
+ }
+
+private:
+ bool isInLens_(const GlobalPosition &globalPos) const
+ {
+ for (int i = 0; i < dimWorld; ++i) {
+ if (globalPos[i] < lensLowerLeft_[i] + eps_ || globalPos[i] > lensUpperRight_[i] - eps_)
+ return false;
+ }
+ return true;
+ }
+
+ GlobalPosition lensLowerLeft_;
+ GlobalPosition lensUpperRight_;
+
+ Scalar lensK_;
+ Scalar outerK_;
+ MaterialLawParams lensMaterialParams_;
+ MaterialLawParams outerMaterialParams_;
+
+ static constexpr Scalar eps_ = 1.5e-7;
+};
+
+} // end namespace Dumux
+
+#endif
--
GitLab