diff --git a/exercises/exercise-properties/CMakeLists.txt b/exercises/exercise-properties/CMakeLists.txt
index 0cc953b240e1125797b72a761dfda59febb70588..0db052ca95a2eaff0ae3b69a58b1fbd582720ae4 100644
--- a/exercises/exercise-properties/CMakeLists.txt
+++ b/exercises/exercise-properties/CMakeLists.txt
@@ -1,7 +1,6 @@
-# the compositional two-phase simulation program
+# the property exercise simulation program
dune_add_test(NAME exercise_properties
SOURCES exercise_properties.cc
- COMPILE_DEFINITIONS TYPETAG=Injection2p2pcCCProblem
CMD_ARGS exercise_properties.input)
# add tutorial to the common target
diff --git a/exercises/exercise-properties/README.md b/exercises/exercise-properties/README.md
index c6e521bdfed0fcb3ed54d115c231e653e1ee6772..9c263bf572d35657d7c5c0bac24f41acbc5ac839 100644
--- a/exercises/exercise-properties/README.md
+++ b/exercises/exercise-properties/README.md
@@ -2,13 +2,13 @@
<br>
## Problem set-up
-The problem setup is identical to the previous [_Exercise Basic_](../exercise-basic/README.md).
+The problem setup is identical to the two-phase incompressible test from DuMu<sup>x</sup>.
## Preparing the exercise
* Navigate to the directory `exercise-properties`
-_Exercise Properties_ deals with a two-phase compositional problem (__2p2c__). Goal is to learn how to use compile and runtime parameters and the _DuMu<sup>x</sup> property system_.
+_Exercise Properties_ deals with a two-phase immiscible incompressible problem (__2p__). The goal is to learn how to adapt compile-time parameters by employing the _DuMu<sup>x</sup> property system_.
<br><br>
### Task 1: Getting familiar with the code
@@ -16,12 +16,11 @@ _Exercise Properties_ deals with a two-phase compositional problem (__2p2c__). G
Locate all the files you will need for this exercise
* The __main file__: `exercise_properties.cc`
-* The __problem file__: `injection2p2cproblem.hh`
-* The __spatial parameters file__: `injection2p2cspatialparams.hh`
+* The __problem file__: `problem.hh`
+* The __spatial parameters file__: `spatialparams.hh`
* The __input file__: `exercise_properties.input`
-* Two header files containing:
+* One header file containing:
* a custom __local residual__ in: `mylocalresidual.hh`
- * a custom __material law__ in: `mymateriallaw.hh`
<hr><br><br>
@@ -46,88 +45,32 @@ make exercise_properties
./exercise_properties
```
Note: Because the input file has the same name as the
-executable, DuMuX will find it automatically.
-
-If gnuplot is installed on your system, you should see a plot of the capillary pressure - saturation relationship.
-
-<hr><br><br>
-### Task 3: Implement and use a different material law
-<hr>
-
-DuMuX uses the term _material law_ to describe the law used to compute
-* pc-Sw relations
-* kr-Sw relations
-* their inverse relations
-
-The file `mymateriallaw.hh` contains a custom implementation of such a material law.
-
-* Implement the method `Scalar pc(const Params ¶ms, Scalar swe)` by implementing your own capillary pressure relationship, e.g. a simple linear relationship $`p_C(S_w) = 1\cdot 10^5 \cdot (1-S_w) + p_e`$.
-
-Note: `MyMaterialLaw` uses the `BrooksCoreyParams` class as parameter input. You can get the entry pressure that is set in the spatial params as follows
-
-```c++
-const auto pe = params.pe();
-```
-
-The type (i.e. C++ type) of the material law is set in the file `injection2p2cspatialparams.hh` by declaring the following alias in the public section of the spatial parameters class:
-
-```c++
-using MaterialLaw = EffToAbsLaw<RegularizedBrooksCorey<Scalar>>;
-```
-
-* Make DuMuX use your own material law by including the header `mymateriallaw.hh` and changing the alias `MaterialLaw`. This will make sure that your material law is used everywhere else in the code.
-
-Note: Also use the wrapper class `EffToAbsLaw`. It takes care of converting absolute to effective saturations considering residual saturations. `MyMaterialLaw`
-as other material laws (like Brooks-Corey, VanGenuchten, ...) in DuMuX only deals with effective saturations.
-
-* Verify your changes by recompiling and running the program. You should see a plot of your new function.
-
-For the next task, disable the plotting feature by changing the settings in the input file `exercise2.input`
-
-```ini
-[Problem]
-OnlyPlotMaterialLaws = false
-```
+executable, DuMu<sup>x</sup> will find it automatically.
<hr><br><br>
-### Task 4: Implement a custom local residual
+### Task 3: Implement a custom local residual
<hr>
-Most types in DuMuX are properties that can be changed just like the material law. In the following task we implement our own 2p2c local residual, i.e. the class that computes the element residual in every Newton iteration. The file `mylocalresidual.hh` contains a copy of the similar to the original local residual class used for all compositional models renamed to `template<class TypeTag> class MyCompositionalLocalResidual`.
+Many types in DuMu<sup>x</sup> are properties that can be changed on the problem level by using the property system. In the following task we implement our own 2p local residual, i.e. the class that computes the element residual in every Newton iteration. The file `mylocalresidual.hh` contains a copy of the original local residual class used for all immiscible models renamed to `template<class TypeTag> class MyLocalResidual`.
-* Make DuMuX use this new local residual by inluding the header `mylocalresidual.hh` and setting the corresponding property in the `Property` namespace in the file `injection2p2cproblem.hh`
+* Make DuMu<sup>x</sup> use this new local residual by including the header `mylocalresidual.hh` and setting the corresponding property in the `Property` namespace in the file `problem.hh`
```c++
// note that every property struct knows about TypeTag
-SET_PROP(Injection2p2cTypeTag, LocalResidual)
+SET_PROP(TwoPIncompressible, LocalResidual)
{
- using type = MyCompositionalLocalResidual<TypeTag>;
+ using type = MyLocalResidual<TypeTag>;
};
// or using the convenience macro
-SET_TYPE_PROP(Injection2p2cTypeTag, LocalResidual,
- MyCompositionalLocalResidual<TypeTag>);
+SET_TYPE_PROP(TwoPIncompressible, LocalResidual,
+ MyLocalResidual<TypeTag>);
```
-You want to make the new local residual special by adding a switch enabling / disabling diffusion. We will achieve this with a DuMuX parameter which is read from the input file and defaults to a property value if the input file doesn't contain the parameter.
+You want to make the simplify the original local residual special by using the assumption that incompressible fluid phases are present. As a consequence, one can balance phase volume instead of phase mass.
-* Modify the `computeFlux` method to only call the `diffusiveFlux` method if diffusion is enabled. You can get the new parameter by adding the lines
-
-```c++
-// ... in the computeFlux method of MyCompositionalLocalResidual
-auto enableDiffusion = getParam<bool>("Problem.EnableDiffusion", true);
-
-```
-
-You can now enable and disable diffusion through the input file
-
-```ini
-[Problem]
-EnableDiffusion = true / false
-```
+* Eliminate the density from the `computeStorage` and the `computeFlux` methods.
-* Verify the difference in the parameter $`x_w^{N2}`$, i.e. the mole fraction of nitrogen in the
-water phase, with and without diffusion.
+* Adapt the relevant routine in the problem file such that volume instead of mass is injected. The density of the employed NAPL is 1460 kg/m^3. For a first try, you can use the hardcoded value.
-Note that due to diffusion being a slow process you
-can only see the difference in later times.
+* Generalize your approach by using the component's `liquidDensity(t, p)` function instead of the hardcoded value.
diff --git a/exercises/exercise-properties/exercise_properties.cc b/exercises/exercise-properties/exercise_properties.cc
index 2a7493e6771bb5620efc20f21eb95683485ab42e..3583260ed9bfd2083e546a8771dca79488c4cae5 100644
--- a/exercises/exercise-properties/exercise_properties.cc
+++ b/exercises/exercise-properties/exercise_properties.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/exercise-properties/exercise_properties.input b/exercises/exercise-properties/exercise_properties.input
index 85ea966262be41ab57d58bc3e3376d1eb5334426..7ae1761e98b6048ded46586398d5d406f2c8fe80 100644
--- a/exercises/exercise-properties/exercise_properties.input
+++ b/exercises/exercise-properties/exercise_properties.input
@@ -1,25 +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
-# TODO: dumux-course-task
-# Set Problem.EnableGravity
-# Set Problem.EnableDiffusion
+[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/exercise-properties/injection2p2cproblem.hh b/exercises/exercise-properties/injection2p2cproblem.hh
deleted file mode 100644
index 3a9023f7a530ab9b9c6da6dacbd12c8e860e640c..0000000000000000000000000000000000000000
--- a/exercises/exercise-properties/injection2p2cproblem.hh
+++ /dev/null
@@ -1,279 +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"
-
-// TODO: dumux-course-task
-// Include the local residual header
-
-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)>);
-
-// TODO: dumux-course-task
-// change the local residual type to MyTwoPTwoCLocalResidual<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 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.
- */
- PrimaryVariables neumannAtPos(const GlobalPosition &globalPos) const
- {
- // initialize values to zero, i.e. no-flow Neumann boundary conditions
- PrimaryVariables 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/exercise-properties/injection2p2cspatialparams.hh b/exercises/exercise-properties/injection2p2cspatialparams.hh
deleted file mode 100644
index ab9951e28c9b296f67aa95ed0fd8843f35ab4f28..0000000000000000000000000000000000000000
--- a/exercises/exercise-properties/injection2p2cspatialparams.hh
+++ /dev/null
@@ -1,193 +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 <dumux/io/gnuplotinterface.hh>
-#include <dumux/io/plotmateriallaw.hh>
-
-namespace Dumux {
-
-/*!
- * \ingroup TwoPTwoCModel
- * \ingroup ImplicitTestProblems
- * \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<RegularizedBrooksCorey<Scalar>>;
- using MaterialLawParams = typename MaterialLaw::Params;
-
-public:
-
- /*!
- * \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/exercise-properties/mylocalresidual.hh b/exercises/exercise-properties/mylocalresidual.hh
index 5c79e24b00230a1fb8beee824436fce0e72c39ab..1a35cf0faa26b8f179d24719f8c5f6c155d04d8d 100644
--- a/exercises/exercise-properties/mylocalresidual.hh
+++ b/exercises/exercise-properties/mylocalresidual.hh
@@ -18,127 +18,113 @@
*****************************************************************************/
/*!
* \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 {
-
+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
+ // TODO: eliminate density from the storage term
+ 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.density(phaseIdx)
+ * 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
+ // TODO: eliminate the density from the flux term
+ 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);
-
- // TODO: dumux-course-task
- // get parameter Problem.EnableDiffusion
- // 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
- flux[eqIdx] += diffusiveFluxes[compIdx];
- }
+ // the physical quantities for which we perform upwinding
+ auto upwindTerm = [phaseIdx](const auto& volVars)
+ { return volVars.density(phaseIdx)*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 +135,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/exercise-properties/mymateriallaw.hh b/exercises/exercise-properties/mymateriallaw.hh
deleted file mode 100644
index e49c39a83c0afb776453782ac299e9cd1d3882c3..0000000000000000000000000000000000000000
--- a/exercises/exercise-properties/mymateriallaw.hh
+++ /dev/null
@@ -1,116 +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
- * TODO: dumux-course-task
- * Implement the pc(swe) function
- */
- static Scalar pc(const Params ¶ms, Scalar swe)
- {
- return 0.0;
- }
-
- /*!
- * \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/exercise-properties/problem.hh b/exercises/exercise-properties/problem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..020272c32ad3ee0f94e11a7db32b462a4c3dfdd0
--- /dev/null
+++ b/exercises/exercise-properties/problem.hh
@@ -0,0 +1,249 @@
+// -*- 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"
+
+// TODO: include the local residual header
+
+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>);
+
+// TODO: use MyLocalResidual as LocalResidual
+
+// 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))
+ values[contiDNAPLEqIdx] = -0.04; // kg / (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/exercise-properties/spatialparams.hh b/exercises/exercise-properties/spatialparams.hh
new file mode 100644
index 0000000000000000000000000000000000000000..275984d814916b4c744251f4c104ed5eb375da1d
--- /dev/null
+++ b/exercises/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