diff --git a/exercises/CMakeLists.txt b/exercises/CMakeLists.txt
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..13f2f3a1b001c7891b0928c54ea4bed4d86d08e7 100644
--- a/exercises/CMakeLists.txt
+++ b/exercises/CMakeLists.txt
@@ -0,0 +1,7 @@
+# add a target that builds all exercise solutions
+add_custom_target(test_exercises)
+
+add_subdirectory(exercise-basic)
+add_subdirectory(exercise-properties)
+add_subdirectory(exercise-fluidsystem)
+add_subdirectory(solution)
diff --git a/exercises/exercise-basic/CMakeLists.txt b/exercises/exercise-basic/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..4b6f85690a08832c9b219fce8f35b5bd75296449
--- /dev/null
+++ b/exercises/exercise-basic/CMakeLists.txt
@@ -0,0 +1,17 @@
+# the immiscible two-phase simulation program
+dune_add_test(NAME exercise1_2p
+ SOURCES exercise1_2p.cc
+ CMD_ARGS exercise1.input)
+
+# the compositional two-phase two-component simulation program
+dune_add_test(NAME exercise1_2p2c
+ SOURCES exercise1_2p2c.cc
+ CMD_ARGS exercise1.input)
+
+# here, add the two-phase non-isothermal simulation program
+
+# add tutorial to the common target
+add_dependencies(test_exercises exercise1_2p exercise1_2p2c)
+
+# add a symlink for the input file
+dune_symlink_to_source_files(FILES "exercise1.input")
diff --git a/exercises/exercise-basic/README.md b/exercises/exercise-basic/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..e81ca62821af5bc60947ff259d848187825b22f1
--- /dev/null
+++ b/exercises/exercise-basic/README.md
@@ -0,0 +1,151 @@
+# Exercise #1 (DuMuX course)
+<br>
+## Problem set-up
+
+N$`_2`$ is injected in an aquifer previously saturated with water with an injection rate of 0.0001 kg/(s*m$`^2`$).
+The aquifer is situated 2700 m below sea level and the domain size is 60 m x 40 m. It consists of two layers, a moderately permeable one ($`\Omega_1`$) and a lower permeable one ($`\Omega_2`$).
+
+<img src="https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/raw/master/tutorial/extradoc/exercise1_setup.png" width="1000">
+
+## Preparing the exercise
+
+* Navigate to the directory `dumux/tutorial/ex1`
+
+_Exercise 1_ deals with two problems: a two-phase immiscible problem (__2p__) and a two-phase compositional problem (__2p2c__). They both set up the same scenario with the difference that the 2p2c assumes a miscible fluid state for the two fluids (water and gaseous N$`_2`$) and the 2p model assumes an immiscible fluid state.
+
+<br><br>
+### Task 1: Getting familiar with the code
+<hr>
+
+Locate all the files you will need for this exercise
+* The __main file__ for the __2p__ problem : `exercise1_2p.cc`
+* The __main file__ for the __2p2c__ problem : `exercise1_2p2c.cc`
+* The __problem file__ for the __2p__ problem: `injection2pproblem.hh`
+* The __problem file__ for the __2p2c__ problem: `injection2p2cproblem.hh`
+* The shared __spatial parameters file__: `injection2pspatialparams.hh`
+* The shared __input file__: `exercise1.input`
+
+<hr><br><br>
+### Task 2: Compiling and running an executable
+<hr>
+
+* Change to the build-directory
+
+```bash
+cd ../../build-cmake/tutorial/ex1
+```
+
+* Compile both executables `exercise1_2p` and `exercise1_2p2c`
+
+```bash
+make exercise1_2p exercise1_2p2c
+```
+
+* Execute the two problems and inspect the result
+
+```bash
+./exercise1_2p exercise1.input
+./exercise1_2p2c exercise1.input
+```
+
+* you can look at the results with paraview
+
+```bash
+paraview injection-2p2c.pvd
+```
+
+<hr><br><br>
+### Task 3: Changing input parameters
+<hr>
+
+In the input file `exercise1.input` you can find the following section
+
+```ini
+[SpatialParams]
+EntryPressureAquitard = 4.5e4
+EntryPressureAquifer = 1e4
+```
+
+* Change the values for the aquitard entry pressure in the input file to a lower value and compare the results with the previous solution. You don't need to recompile the executable.
+
+<hr><br><br>
+### Task 4: Runtime parameters
+<hr>
+
+The injection rate is currently hard-coded in `injection2p2cproblem.hh` to $`1e-4 kg/(s m^2)`$.
+
+```c++
+ // set the Neumann values for the Nitrogen component balance
+ // convert from units kg/(s*m^2) to mole/(s*m^2)
+values[Indices::contiNEqIdx] = -1e-4/FluidSystem::molarMass(FluidSystem::nCompIdx);
+values[Indices::contiWEqIdx] = 0.0;
+```
+
+We want to be able to set it at runtime. To this end,
+* use the following DuMuX macro to read a runtime parameter from the input file
+
+```c++
+// read the injection rate from the input file at run time
+totalAreaSpecificInflow_ = getParam<TYPE>("GROUPNAME.PARAMNAME");
+```
+
+* Replace
+`<TYPE>`,`<GROUPNAME>`,`<PARAMNAME>` by what is appropriate for your case:
+ * `<TYPE>` is the type of the parameter to read
+ * `<GROUPNAME>` is the group in the input file
+ * `<PARAMNAME>` is the name of the parameter in the input file
+
+Note that due to the way the macro works, the names are specified as plain text within the "quotation marks".`<GROUPNAME>` and `<PARAMNAME>` need to be separated by a dot (.).
+Follow the instructions given as a
+
+```c++
+// TODO: dumux-course-task
+```
+in the `injection2p2cproblem.hh` file and also remember to also set the parameter totalAreaSpecificInflow in the input file.
+
+* Check the influence of that parameter on the simulation result by rerunning the simulation with different injection rates. Remember to also set the parameter totalAreaSpecificInflow in the input file.
+
+Since you have changed your header file, you have to recompile the program.
+
+<hr><br><br>
+### 5. Setting up a new executable (for a non-isothermal simulation)
+<hr>
+
+* Copy the main file `exercise1_2p.cc` and rename it to `exercise1_2pni.cc`
+* In `exercise1_2pni.cc`, include the header `injection2pniproblem.hh` instead of `injection2pproblem.hh`.
+* In `exercise1_2pni.cc`, change `Injection2pCCTypeTag` to `Injection2pNICCTypeTag` in the line `using TypeTag = TTAG(Injection2pCCTypeTag);`
+* Add a new executable in `CMakeLists.txt` by adding the lines
+
+```cmake
+# the two-phase non-isothermal simulation program
+dune_add_test(NAME exercise1_2pni
+ SOURCES exercise1_2pni.cc
+ CMD_ARGS exercise1.input)
+```
+
+* Test that everything compiles without error
+
+```bash
+make # should rerun cmake
+make exercise1_2pni # builds new executable
+```
+
+<hr><br><br>
+### 6. Setting up a non-isothermal __2pni__ test problem
+<hr>
+
+* Open the file `injection2pniproblem.hh`. It is a copy of the `injection2pproblem.hh` with some useful comments on how to implement a non-isothermal model. Look for comments containing
+
+```c++
+// TODO: dumux-course-task
+```
+
+* The following set-up should be realized:
+
+ __Boundary conditions:__ Dirichlet conditions for the temperature with a temperature gradient of 0.03 K/m and a starting temperature of 283 K.
+
+ __Initial conditions:__ The same temperature gradient as in the boundary conditions with an additional lens (with position: 20 < x < 30, 5 < y < 35), which has an initial temperature of 380 K.
+
+<img src="https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/raw/master/tutorial/extradoc/exercise1_nonisothermal.png" width="800">
+
+The non-isothermal model requires additional parameters like the thermal conductivity of the solid component. They are already implemented and set in `exercise1.input`, you just need to _uncomment_ them.
diff --git a/exercises/exercise-basic/exercise1.cc b/exercises/exercise-basic/exercise1.cc
new file mode 100644
index 0000000000000000000000000000000000000000..3bf3e8d87c4cc32d2f4db34d9ecef6ade0aee790
--- /dev/null
+++ b/exercises/exercise-basic/exercise1.cc
@@ -0,0 +1,197 @@
+// -*- 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 The main file for the two-phase porousmediumflow problem of exercise 1
+ */
+#include <config.h>
+
+#include <ctime>
+#include <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/timer.hh>
+#include <dune/grid/io/file/dgfparser/dgfexception.hh>
+#include <dune/grid/io/file/vtk.hh>
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/dumuxmessage.hh>
+#include <dumux/common/defaultusagemessage.hh>
+
+#include <dumux/linear/amgbackend.hh>
+#include <dumux/nonlinear/newtonsolver.hh>
+
+#include <dumux/assembly/fvassembler.hh>
+#include <dumux/assembly/diffmethod.hh>
+
+#include <dumux/discretization/methods.hh>
+
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/grid/gridmanager.hh>
+
+#include "injection2pproblem.hh"
+
+////////////////////////
+// the main function
+////////////////////////
+int main(int argc, char** argv) try
+{
+ using namespace Dumux;
+
+ // define the type tag for this problem
+ using TypeTag = TTAG(Injection2pCCTypeTag);
+
+ // initialize MPI, finalize is done automatically on exit
+ const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+ // print dumux start message
+ if (mpiHelper.rank() == 0)
+ DumuxMessage::print(/*firstCall=*/true);
+
+ // parse command line arguments and input file
+ Parameters::init(argc, argv);
+
+ // try to create a grid (from the given grid file or the input file)
+ GridManager<typename GET_PROP_TYPE(TypeTag, Grid)> gridManager;
+ gridManager.init();
+
+ ////////////////////////////////////////////////////////////
+ // run instationary non-linear problem on this grid
+ ////////////////////////////////////////////////////////////
+
+ // we compute on the leaf grid view
+ const auto& leafGridView = gridManager.grid().leafGridView();
+
+ // create the finite volume grid geometry
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
+ fvGridGeometry->update();
+
+ // the problem (initial and boundary conditions)
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ auto problem = std::make_shared<Problem>(fvGridGeometry);
+
+ // the solution vector
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ SolutionVector x(fvGridGeometry->numDofs());
+ problem->applyInitialSolution(x);
+ auto xOld = x;
+
+ // the grid variables
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+ auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
+ gridVariables->init(x, xOld);
+
+ // get some time loop parameters
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+ const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+ auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+
+ // 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
+
+ // instantiate time loop
+ auto timeLoop = std::make_shared<TimeLoop<Scalar>>(0.0, dt, tEnd);
+ timeLoop->setMaxTimeStepSize(maxDt);
+
+ // the assembler with time loop for instationary problem
+ using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
+ auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop);
+
+ // the linear solver
+ using LinearSolver = AMGBackend<TypeTag>;
+ auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->dofMapper());
+
+ // the non-linear solver
+ using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
+ NewtonSolver nonLinearSolver(assembler, linearSolver);
+
+ // time loop
+ timeLoop->start();
+ while (!timeLoop->finished())
+ {
+ // 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);
+
+ // make the new solution the old solution
+ xOld = x;
+ gridVariables->advanceTimeStep();
+
+ // advance to the time loop to the next step
+ timeLoop->advanceTimeStep();
+
+ // report statistics of this time step
+ timeLoop->reportTimeStep();
+
+ // set new dt as suggested by the newton solver
+ timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
+
+ // output to vtk
+ vtkWriter.write(timeLoop->time());
+ }
+
+ timeLoop->finalize(leafGridView.comm());
+
+ ////////////////////////////////////////////////////////////
+ // finalize, print dumux message to say goodbye
+ ////////////////////////////////////////////////////////////
+
+ // print dumux end message
+ if (mpiHelper.rank() == 0)
+ {
+ Parameters::print();
+ DumuxMessage::print(/*firstCall=*/false);
+ }
+
+ return 0;
+} // end main
+catch (Dumux::ParameterException &e)
+{
+ std::cerr << std::endl << e << " ---> Abort!" << std::endl;
+ return 1;
+}
+catch (Dune::DGFException & e)
+{
+ std::cerr << "DGF exception thrown (" << e <<
+ "). Most likely, the DGF file name is wrong "
+ "or the DGF file is corrupted, "
+ "e.g. missing hash at end of file or wrong number (dimensions) of entries."
+ << " ---> Abort!" << std::endl;
+ return 2;
+}
+catch (Dune::Exception &e)
+{
+ std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
+ return 3;
+}
+catch (...)
+{
+ std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
+ return 4;
+}
diff --git a/exercises/exercise-basic/exercise1.input b/exercises/exercise-basic/exercise1.input
new file mode 100644
index 0000000000000000000000000000000000000000..792afb3c7b4f9373c4d327572790a4b13534a19a
--- /dev/null
+++ b/exercises/exercise-basic/exercise1.input
@@ -0,0 +1,29 @@
+[TimeLoop]
+DtInitial = 3600 # in seconds
+TEnd = 3.154e9 # in seconds, i.e ten years
+
+[Grid]
+LowerLeft = 0 0
+UpperRight = 60 40
+Cells = 24 16
+
+[Problem]
+Name = injection
+OnlyPlotMaterialLaws = true
+AquiferDepth = 2700.0 # m
+InjectionDuration = 2.628e6 # in seconds, i.e. one month
+
+#TODO: dumux-course-task:
+#set totalAreaSpecificInflow
+
+[SpatialParams]
+PermeabilityAquitard = 1e-15 # m^2
+EntryPressureAquitard = 4.5e4 # Pa
+PermeabilityAquifer = 1e-12 # m^2
+EntryPressureAquifer = 1e4 # Pa
+
+# these parameters are only used in the nonisothermal model. Uncomment them for that
+#[Component]
+#SolidDensity = 2700 # solid density of granite
+#SolidThermalConductivity = 2.8 # solid thermal conducitivity of granite
+#SolidHeatCapacity = 790 # solid heat capacity of granite
diff --git a/exercises/exercise-basic/exercise1_2p.cc b/exercises/exercise-basic/exercise1_2p.cc
new file mode 100644
index 0000000000000000000000000000000000000000..558f2d0f3373ed23df1fefbb0af5c296a698484a
--- /dev/null
+++ b/exercises/exercise-basic/exercise1_2p.cc
@@ -0,0 +1,200 @@
+// -*- 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 The main file for the two-phase porousmediumflow problem of exercise 1
+ */
+#include <config.h>
+
+#include <ctime>
+#include <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/timer.hh>
+#include <dune/grid/io/file/dgfparser/dgfexception.hh>
+#include <dune/grid/io/file/vtk.hh>
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/dumuxmessage.hh>
+#include <dumux/common/defaultusagemessage.hh>
+
+#include <dumux/linear/amgbackend.hh>
+#include <dumux/nonlinear/newtonsolver.hh>
+
+#include <dumux/assembly/fvassembler.hh>
+#include <dumux/assembly/diffmethod.hh>
+
+#include <dumux/discretization/methods.hh>
+
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/grid/gridmanager.hh>
+
+// The problem file, where setup-specific boundary and initial conditions are defined.
+#include "injection2pproblem.hh"
+
+////////////////////////
+// the main function
+////////////////////////
+int main(int argc, char** argv) try
+{
+ using namespace Dumux;
+
+ // define the type tag for this problem
+ using TypeTag = TTAG(Injection2pCCTypeTag);
+
+ // initialize MPI, finalize is done automatically on exit
+ const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+ // print dumux start message
+ if (mpiHelper.rank() == 0)
+ DumuxMessage::print(/*firstCall=*/true);
+
+ // parse command line arguments and input file
+ Parameters::init(argc, argv);
+
+ // try to create a grid (from the given grid file or the input file)
+ GridManager<typename GET_PROP_TYPE(TypeTag, Grid)> gridManager;
+ gridManager.init();
+
+ ////////////////////////////////////////////////////////////
+ // run instationary non-linear problem on this grid
+ ////////////////////////////////////////////////////////////
+
+ // we compute on the leaf grid view
+ const auto& leafGridView = gridManager.grid().leafGridView();
+
+ // create the finite volume grid geometry
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
+ fvGridGeometry->update();
+
+ // the problem (initial and boundary conditions)
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ auto problem = std::make_shared<Problem>(fvGridGeometry);
+
+ // the solution vector
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ SolutionVector x(fvGridGeometry->numDofs());
+ problem->applyInitialSolution(x);
+ auto xOld = x;
+
+ // the grid variables
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+ auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
+ gridVariables->init(x, xOld);
+
+ // get some time loop parameters
+ // getParam<TYPE>("GROUPNAME.PARAMNAME") reads and sets parameter PARAMNAME
+ // of type TYPE given in the group GROUPNAME from the input file
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+ const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+ auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+
+ // 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
+
+ // instantiate time loop
+ auto timeLoop = std::make_shared<TimeLoop<Scalar>>(0.0, dt, tEnd);
+ timeLoop->setMaxTimeStepSize(maxDt);
+
+ // the assembler with time loop for instationary problem
+ using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
+ auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop);
+
+ // the linear solver
+ using LinearSolver = AMGBackend<TypeTag>;
+ auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->dofMapper());
+
+ // the non-linear solver
+ using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
+ NewtonSolver nonLinearSolver(assembler, linearSolver);
+
+ // time loop
+ timeLoop->start();
+ while (!timeLoop->finished())
+ {
+ // 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);
+
+ // make the new solution the old solution
+ xOld = x;
+ gridVariables->advanceTimeStep();
+
+ // advance to the time loop to the next step
+ timeLoop->advanceTimeStep();
+
+ // report statistics of this time step
+ timeLoop->reportTimeStep();
+
+ // set new dt as suggested by the newton solver
+ timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
+
+ // output to vtk
+ vtkWriter.write(timeLoop->time());
+ }
+
+ timeLoop->finalize(leafGridView.comm());
+
+ ////////////////////////////////////////////////////////////
+ // finalize, print dumux message to say goodbye
+ ////////////////////////////////////////////////////////////
+
+ // print dumux end message
+ if (mpiHelper.rank() == 0)
+ {
+ Parameters::print();
+ DumuxMessage::print(/*firstCall=*/false);
+ }
+
+ return 0;
+} // end main
+catch (Dumux::ParameterException &e)
+{
+ std::cerr << std::endl << e << " ---> Abort!" << std::endl;
+ return 1;
+}
+catch (Dune::DGFException & e)
+{
+ std::cerr << "DGF exception thrown (" << e <<
+ "). Most likely, the DGF file name is wrong "
+ "or the DGF file is corrupted, "
+ "e.g. missing hash at end of file or wrong number (dimensions) of entries."
+ << " ---> Abort!" << std::endl;
+ return 2;
+}
+catch (Dune::Exception &e)
+{
+ std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
+ return 3;
+}
+catch (...)
+{
+ std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
+ return 4;
+}
diff --git a/exercises/exercise-basic/exercise1_2p2c.cc b/exercises/exercise-basic/exercise1_2p2c.cc
new file mode 100644
index 0000000000000000000000000000000000000000..c2b46ea15fcab30a078f85acf0f8d07b122b8209
--- /dev/null
+++ b/exercises/exercise-basic/exercise1_2p2c.cc
@@ -0,0 +1,201 @@
+// -*- 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 The main file for the 2p2c porousmediumflow problem in exercise1
+ */
+#include <config.h>
+
+#include <ctime>
+#include <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/timer.hh>
+#include <dune/grid/io/file/dgfparser/dgfexception.hh>
+#include <dune/grid/io/file/vtk.hh>
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/valgrind.hh>
+#include <dumux/common/dumuxmessage.hh>
+
+#include <dumux/linear/amgbackend.hh>
+#include <dumux/nonlinear/privarswitchnewtonsolver.hh>
+
+#include <dumux/assembly/fvassembler.hh>
+#include <dumux/assembly/diffmethod.hh>
+
+#include <dumux/discretization/methods.hh>
+
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/grid/gridmanager.hh>
+
+// The problem file, where setup-specific boundary and initial conditions are defined.
+#include "injection2p2cproblem.hh"
+
+////////////////////////
+// the main function
+////////////////////////
+int main(int argc, char** argv) try
+{
+ using namespace Dumux;
+
+ // define the type tag for this problem
+ using TypeTag = TTAG(Injection2p2cCCTypeTag);
+
+ // initialize MPI, finalize is done automatically on exit
+ const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+ // print dumux start message
+ if (mpiHelper.rank() == 0)
+ DumuxMessage::print(/*firstCall=*/true);
+
+ // parse command line arguments and input file
+ Parameters::init(argc, argv);
+
+ // try to create a grid (from the given grid file or the input file)
+ GridManager<typename GET_PROP_TYPE(TypeTag, Grid)> gridManager;
+ gridManager.init();
+
+ ////////////////////////////////////////////////////////////
+ // run instationary non-linear problem on this grid
+ ////////////////////////////////////////////////////////////
+
+ // we compute on the leaf grid view
+ const auto& leafGridView = gridManager.grid().leafGridView();
+
+ // create the finite volume grid geometry
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
+ fvGridGeometry->update();
+
+ // the problem (initial and boundary conditions)
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ auto problem = std::make_shared<Problem>(fvGridGeometry);
+
+ // the solution vector
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ SolutionVector x(fvGridGeometry->numDofs());
+ problem->applyInitialSolution(x);
+ auto xOld = x;
+
+ // the grid variables
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+ auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
+ gridVariables->init(x, xOld);
+
+ // get some time loop parameters
+ // getParam<TYPE>("GROUPNAME.PARAMNAME") reads and sets parameter PARAMNAME
+ // of type TYPE given in the group GROUPNAME from the input file
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+ const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+ auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+
+ // 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
+
+ // instantiate time loop
+ auto timeLoop = std::make_shared<TimeLoop<Scalar>>(0.0, dt, tEnd);
+ timeLoop->setMaxTimeStepSize(maxDt);
+
+ // the assembler with time loop for instationary problem
+ using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
+ auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop);
+
+ // the linear solver
+ using LinearSolver = AMGBackend<TypeTag>;
+ 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>;
+ NewtonSolver nonLinearSolver(assembler, linearSolver);
+
+ // time loop
+ timeLoop->start();
+ while (!timeLoop->finished())
+ {
+ // 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);
+
+ // make the new solution the old solution
+ xOld = x;
+ gridVariables->advanceTimeStep();
+
+ // advance to the time loop to the next step
+ timeLoop->advanceTimeStep();
+
+ // report statistics of this time step
+ timeLoop->reportTimeStep();
+
+ // set new dt as suggested by the newton solver
+ timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
+
+ // output to vtk
+ vtkWriter.write(timeLoop->time());
+ }
+
+ timeLoop->finalize(leafGridView.comm());
+
+ ////////////////////////////////////////////////////////////
+ // finalize, print dumux message to say goodbye
+ ////////////////////////////////////////////////////////////
+
+ // print dumux end message
+ if (mpiHelper.rank() == 0)
+ {
+ Parameters::print();
+ DumuxMessage::print(/*firstCall=*/false);
+ }
+
+ return 0;
+} // end main
+catch (Dumux::ParameterException &e)
+{
+ std::cerr << std::endl << e << " ---> Abort!" << std::endl;
+ return 1;
+}
+catch (Dune::DGFException & e)
+{
+ std::cerr << "DGF exception thrown (" << e <<
+ "). Most likely, the DGF file name is wrong "
+ "or the DGF file is corrupted, "
+ "e.g. missing hash at end of file or wrong number (dimensions) of entries."
+ << " ---> Abort!" << std::endl;
+ return 2;
+}
+catch (Dune::Exception &e)
+{
+ std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
+ return 3;
+}
+catch (...)
+{
+ std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
+ return 4;
+}
diff --git a/exercises/exercise-basic/injection2p2cproblem.hh b/exercises/exercise-basic/injection2p2cproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..3c30c2a2bf971657cb8da4c06dc360bacabb1acd
--- /dev/null
+++ b/exercises/exercise-basic/injection2p2cproblem.hh
@@ -0,0 +1,275 @@
+// -*- 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 The two-phase porousmediumflow problem for exercise 1
+ */
+#ifndef DUMUX_EX1_INJECTION_2P2C_PROBLEM_HH
+#define DUMUX_EX1_INJECTION_2P2C_PROBLEM_HH
+
+#include <dumux/discretization/cellcentered/tpfa/properties.hh>
+#include <dumux/porousmediumflow/2p2c/model.hh>
+#include <dumux/porousmediumflow/problem.hh>
+#include <dumux/material/fluidsystems/h2on2.hh>
+
+#include "injection2pspatialparams.hh"
+
+namespace Dumux {
+
+// forward declare problem
+template <class TypeTag>
+class Injection2p2cProblem;
+
+namespace Properties {
+NEW_TYPE_TAG(Injection2p2cTypeTag, INHERITS_FROM(TwoPTwoC));
+NEW_TYPE_TAG(Injection2p2cCCTypeTag, 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 fluid configuration
+SET_TYPE_PROP(Injection2p2cTypeTag, FluidSystem, FluidSystems::H2ON2<typename GET_PROP_TYPE(TypeTag, Scalar), /*useComplexRelations=*/ false>);
+
+// Define whether mole(true) or mass (false) fractions are used
+SET_BOOL_PROP(Injection2p2cTypeTag, UseMoles, true);
+} // end namespace Properties
+
+/*!
+ * \ingroup TwoPTwoCModel
+ * \ingroup ImplicitTestProblems
+ * \brief Gas injection problem where a gas (here nitrogen) is injected into a fully
+ * water saturated medium. During buoyancy driven upward migration the gas
+ * passes a high temperature area.
+ *
+ * The domain is sized 60 m times 40 m.
+ *
+ * For the mass conservation equation neumann boundary conditions are used on
+ * the top, on the bottom and on the right of the domain, while dirichlet conditions
+ * apply on the left boundary.
+ *
+ * Gas is injected at the right boundary from 7 m to 15 m at a rate of
+ * 0.001 kg/(s m), the remaining neumann boundaries are no-flow
+ * boundaries.
+ *
+ * At the dirichlet boundaries a hydrostatic pressure and a gas saturation of zero a
+ *
+ * This problem uses the \ref TwoPModel model.
+ */
+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 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;
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+
+ 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
+ // getParam<TYPE>("GROUPNAME.PARAMNAME") reads and sets parameter PARAMNAME
+ // of type TYPE given in the group GROUPNAME from the input file
+ name_ = getParam<std::string>("Problem.Name");
+ // depth of the aquifer, units: m
+ aquiferDepth_ = getParam<Scalar>("Problem.AquiferDepth");
+ // the duration of the injection, units: second
+ injectionDuration_ = getParam<Scalar>("Problem.InjectionDuration");
+
+ // TODO: dumux-course-task
+ // Get the specific inflow of 1e-4 kg/(s m^2) from the input file (totalAreaSpecificInflow_) here as it is done for the injectionDuration_.
+
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief Returns the problem name
+ *
+ * This is used as a prefix for files generated by the simulation.
+ */
+ std::string name() const
+ { return name_+"-2p2c"; }
+
+ /*!
+ * \brief Returns the temperature \f$ K \f$
+ */
+ Scalar temperature() const
+ {
+ return 273.15 + 30; // [K]
+ }
+
+ // \}
+
+ /*!
+ * \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();
+ 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])
+ {
+ // TODO: dumux-course-task
+ //instead of setting -1e-4 here directly use totalAreaSpecificInflow_ in the computation
+
+ // inject nitrogen. negative values mean injection
+ // convert from units kg/(s*m^2) to mole/(s*m^2)
+ values[Indices::conti0EqIdx + FluidSystem::N2Idx] = -1e-4/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::pressureIdx] = pw;
+ values[Indices::switchIdx] = moleFracLiquidN2;
+
+ return values;
+ }
+
+ // \}
+
+ //! set the time for the time dependent boundary conditions (called from main)
+ 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 injectionDuration_; //! Duration of the injection in seconds
+ Scalar time_;
+ //TODO: dumux-course-task
+ //define the Scalar totalAreaSpecificInflow_ here
+
+};
+
+} //end namespace Dumux
+
+#endif
diff --git a/exercises/exercise-basic/injection2pniproblem.hh b/exercises/exercise-basic/injection2pniproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..ad241c88f28575338cd3cdc54579b2a37ad11d0b
--- /dev/null
+++ b/exercises/exercise-basic/injection2pniproblem.hh
@@ -0,0 +1,266 @@
+// -*- 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 The two-phase nonisothermal porousmediumflow problem for exercise 1
+ */
+
+#ifndef DUMUX_EX1_INJECTION_PROBLEM_2PNI_HH
+#define DUMUX_EX1_INJECTION_PROBLEM_2PNI_HH
+
+#include <dumux/discretization/cellcentered/tpfa/properties.hh>
+#include <dumux/porousmediumflow/2p/model.hh>
+#include <dumux/porousmediumflow/problem.hh>
+#include <dumux/material/fluidsystems/h2on2.hh>
+
+#include "injection2pspatialparams.hh"
+
+namespace Dumux {
+
+// forward declare problem
+template <class TypeTag>
+class InjectionProblem2PNI;
+
+namespace Properties
+{
+ /*!
+* TODO:dumux-course-task:
+* inherit from the TwoPNI model instead of TwoP here
+*/
+NEW_TYPE_TAG(Injection2pNITypeTag, INHERITS_FROM(TwoP));
+NEW_TYPE_TAG(Injection2pNICCTypeTag, INHERITS_FROM(CCTpfaModel, Injection2pNITypeTag));
+
+// Set the grid type
+SET_TYPE_PROP(Injection2pNITypeTag, Grid, Dune::YaspGrid<2>);
+
+// Set the problem property
+SET_TYPE_PROP(Injection2pNITypeTag, Problem, InjectionProblem2PNI<TypeTag>);
+
+// Set the spatial parameters
+SET_TYPE_PROP(Injection2pNITypeTag, SpatialParams,
+ InjectionSpatialParams<typename GET_PROP_TYPE(TypeTag, FVGridGeometry),
+ typename GET_PROP_TYPE(TypeTag, Scalar)>);
+
+// Set fluid configuration
+SET_TYPE_PROP(Injection2pNITypeTag, FluidSystem, FluidSystems::H2ON2<typename GET_PROP_TYPE(TypeTag, Scalar), false>);
+} // end namespace Properties
+
+/*!
+ * \ingroup TwoPModel
+ * \ingroup ImplicitTestProblems
+ * \brief Gas injection problem where a gas (here nitrogen) is injected into a fully
+ * water saturated medium. During buoyancy driven upward migration the gas
+ * passes a high temperature area.
+ *
+ * The domain is sized 60 m times 40 m.
+ *
+ * For the mass conservation equation neumann boundary conditions are used on
+ * the top, on the bottom and on the right of the domain, while dirichlet conditions
+ * apply on the left boundary.
+ *
+ * Gas is injected at the right boundary from 7 m to 15 m at a rate of
+ * 0.001 kg/(s m), the remaining neumann boundaries are no-flow
+ * boundaries.
+ *
+ * At the dirichlet boundaries a hydrostatic pressure and a gas saturation of zero a
+ *
+ * This problem uses the \ref TwoPModel model.
+ */
+template<class TypeTag>
+class InjectionProblem2PNI : public PorousMediumFlowProblem<TypeTag>
+{
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ 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;
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+
+ enum { dimWorld = GridView::dimensionworld };
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+public:
+ InjectionProblem2PNI(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
+ // getParam<TYPE>("GROUPNAME.PARAMNAME") reads and sets parameter PARAMNAME
+ // of type TYPE given in the group GROUPNAME from the input file
+ name_ = getParam<std::string>("Problem.Name");
+ // depth of the aquifer, units: m
+ aquiferDepth_ = getParam<Scalar>("Problem.AquiferDepth");
+ // 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.
+ */
+ std::string name() const
+ { return "injection-2pni"; }
+
+ // \}
+
+ /*!
+ * \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();
+ 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])
+ {
+ // inject nitrogen. negative values mean injection
+ // units kg/(s*m^2)
+ values[Indices::conti0EqIdx + FluidSystem::N2Idx] = -1e-4;
+ values[Indices::conti0EqIdx + FluidSystem::H2OIdx] = 0.0;
+
+ /*!
+ * TODO:dumux-course-task:
+ * dumux-course-task:
+ * set Neumann noflow conditions for the energy equation everywhere
+ * hint: use Indices::energyEqIdx) for that
+ */
+ }
+
+ 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);
+
+ // get the water density at atmospheric conditions
+ const Scalar densityW = FluidSystem::H2O::liquidDensity(283.15, 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]);
+
+ values[Indices::pressureIdx] = pw;
+ values[Indices::saturationIdx] = 0.0;
+
+ /*!
+ * TODO:dumux-course-task:
+ * set a temperature gradient of 0.03 K per m beginning at 283 K here.
+ * Hint: you can use aquiferDepth_ and the globalPos similar to the pressure gradient
+ * use globalPos[0] and globalPos[1] to implement the high temperature lens with 380 K
+ * Hint : use Indices::temperatureIdx to address the initial values for temperature
+ */
+ return values;
+ }
+
+ // \}
+
+ //! set the time for the time dependent boundary conditions (called from main)
+ 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 injectionDuration_; //! Duration of the injection in seconds
+ Scalar time_;
+};
+
+} //end namespace Dumux
+
+#endif
diff --git a/exercises/exercise-basic/injection2pproblem.hh b/exercises/exercise-basic/injection2pproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..46b5d38930906614e8d279085bb5fdd346fd744f
--- /dev/null
+++ b/exercises/exercise-basic/injection2pproblem.hh
@@ -0,0 +1,260 @@
+// -*- 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 The two-phase porousmediumflow problem for exercise 1
+ */
+
+#ifndef DUMUX_EX1_INJECTION_PROBLEM_2P_HH
+#define DUMUX_EX1_INJECTION_PROBLEM_2P_HH
+
+#include <dumux/discretization/cellcentered/tpfa/properties.hh>
+#include <dumux/porousmediumflow/2p/model.hh>
+#include <dumux/porousmediumflow/problem.hh>
+#include <dumux/material/fluidsystems/h2on2.hh>
+
+#include "injection2pspatialparams.hh"
+
+namespace Dumux {
+
+// forward declare problem
+template <class TypeTag>
+class InjectionProblem2P;
+
+namespace Properties {
+// define the TypeTag for this problem with a cell-centered two-point flux approximation spatial discretization.
+NEW_TYPE_TAG(Injection2pTypeTag, INHERITS_FROM(TwoP));
+NEW_TYPE_TAG(Injection2pCCTypeTag, INHERITS_FROM(CCTpfaModel, Injection2pTypeTag));
+
+// Set the grid type
+SET_TYPE_PROP(Injection2pTypeTag, Grid, Dune::YaspGrid<2>);
+
+// Set the problem property
+SET_TYPE_PROP(Injection2pTypeTag, Problem, InjectionProblem2P<TypeTag>);
+
+// Set the spatial parameters
+SET_TYPE_PROP(Injection2pTypeTag, SpatialParams,
+ InjectionSpatialParams<typename GET_PROP_TYPE(TypeTag, FVGridGeometry),
+ typename GET_PROP_TYPE(TypeTag, Scalar)>);
+
+// Set fluid configuration
+SET_TYPE_PROP(Injection2pTypeTag, FluidSystem, FluidSystems::H2ON2<typename GET_PROP_TYPE(TypeTag, Scalar), /*useComplexRelations=*/ false>);
+} // end namespace Properties
+
+/*!
+ * \ingroup TwoPModel
+ * \ingroup ImplicitTestProblems
+ * \brief Gas injection problem where a gas (here nitrogen) is injected into a fully
+ * water saturated medium. During buoyancy driven upward migration the gas
+ * passes a high temperature area.
+ *
+ * The domain is sized 60 m times 40 m.
+ *
+ * For the mass conservation equation neumann boundary conditions are used on
+ * the top, on the bottom and on the right of the domain, while dirichlet conditions
+ * apply on the left boundary.
+ *
+ * Gas is injected at the right boundary from 7 m to 15 m at a rate of
+ * 0.001 kg/(s m), the remaining neumann boundaries are no-flow
+ * boundaries.
+ *
+ * At the dirichlet boundaries a hydrostatic pressure and a gas saturation of zero a
+ *
+ * This problem uses the \ref TwoPModel model.
+ */
+template<class TypeTag>
+class InjectionProblem2P : public PorousMediumFlowProblem<TypeTag>
+{
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ 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;
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+
+ enum { dimWorld = GridView::dimensionworld };
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+public:
+ InjectionProblem2P(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
+ // getParam<TYPE>("GROUPNAME.PARAMNAME") reads and sets parameter PARAMNAME
+ // of type TYPE given in the group GROUPNAME from the input file
+ name_ = getParam<std::string>("Problem.Name");
+ // depth of the aquifer, units: m
+ aquiferDepth_ = getParam<Scalar>("Problem.AquiferDepth");
+ // 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.
+ */
+ std::string name() const
+ { return name_+"-2p"; }
+
+ /*!
+ * \brief Returns the temperature \f$ K \f$
+ */
+ Scalar temperature() const
+ {
+ return 273.15 + 30; // [K]
+ }
+
+ // \}
+
+ /*!
+ * \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;
+ // set the left of the domain (with the global position in "0 = x" direction as a Dirichlet boundary
+ if (globalPos[0] < eps_)
+ bcTypes.setAllDirichlet();
+ // set all other as Neumann boundaries
+ 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
+ // using < boundary + eps_ or > boundary - eps_ is safer for floating point comparisons
+ // than using <= or >= as it is robust with regard to imprecision introduced by rounding errors.
+ if (time_ < injectionDuration_
+ && globalPos[1] < 15 + eps_ && globalPos[1] > 7 - eps_ && globalPos[0] > 0.9*this->fvGridGeometry().bBoxMax()[0])
+ {
+ // inject nitrogen. negative values mean injection
+ // units kg/(s*m^2)
+ values[Indices::conti0EqIdx + FluidSystem::N2Idx] = -1e-4;
+ 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);
+
+ // 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]);
+
+ values[Indices::pressureIdx] = pw;
+ values[Indices::saturationIdx] = 0.0;
+
+ return values;
+ }
+
+ // \}
+
+ //! set the time for the time dependent boundary conditions (called from main)
+ 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 injectionDuration_; //! Duration of the injection in seconds
+ Scalar time_;
+};
+
+} //end namespace Dumux
+
+#endif
diff --git a/exercises/exercise-basic/injection2pspatialparams.hh b/exercises/exercise-basic/injection2pspatialparams.hh
new file mode 100644
index 0000000000000000000000000000000000000000..fac5a1da9896f747dc94b9b9787383976c063f86
--- /dev/null
+++ b/exercises/exercise-basic/injection2pspatialparams.hh
@@ -0,0 +1,171 @@
+// -*- 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_EX1_INJECTION_SPATIAL_PARAMS_HH
+#define DUMUX_EX1_INJECTION_SPATIAL_PARAMS_HH
+
+#include <dumux/material/spatialparams/fv.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.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;
+
+ // get the dimensions of the simulation domain from 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;
+
+ using MaterialLaw = EffToAbsLaw<RegularizedBrooksCorey<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);
+ }
+
+ /*!
+ * \brief Define the intrinsic permeability \f$\mathrm{[m^2]}\f$.
+ *
+ * \param globalPos The global position
+ */
+ PermeabilityType permeabilityAtPos(const GlobalPosition& globalPos) const
+ {
+ // here, either aquitard or aquifer permeability are returned, depending on the global position
+ 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
+ {
+ // here, either aquitard or aquifer porosity are returned, depending on the global position
+ 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; }
+
+private:
+
+ static constexpr Scalar eps_ = 1e-6;
+
+ // provides a convenient way distinguishing whether a given location is inside the aquitard
+ bool isInAquitard_(const GlobalPosition &globalPos) const
+ {
+ // globalPos[dimWorld-1] is the y direction for 2D grids or the z direction for 3D grids
+ 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-dunemodule/README.md b/exercises/exercise-dunemodule/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..dff7efc942fec3089b220b3f539223c06836bf49
--- /dev/null
+++ b/exercises/exercise-dunemodule/README.md
@@ -0,0 +1,92 @@
+# Exercise #4 (DuMuX course)
+
+This exercise describes how to create a new DuMuX module
+and how to create a corresponding GitLab project.
+
+This is the suggested
+workflow to develop code on top of DuMuX.
+
+### Task 1: Create new dune module
+<hr>
+
+* Execute the following command (bash environment) in the top-folder, i.e. above the dumux folder
+
+```bash
+./dune-common/bin/duneproject
+```
+
+* Follow the introductions and specify
+ * as name of the new module: `dumux-example`
+ * as module dependencies: `dumux`
+ * a version at your choice
+ * your email address
+
+<hr><br><br>
+### Task 2: Rerun dunecontrol to configure your new project
+<hr>
+
+The following command will configure your new module
+
+```bash
+./dune-common/bin/dunecontrol --opts=<opts file> --only=dumux-example all
+```
+
+<hr><br><br>
+### Task 3: Create a new test case within your new DuMuX module
+<hr>
+
+* Create a new folder (in your module folder), e.g. `appl`
+
+```bash
+mkdir appl
+```
+
+* Copy some test case from the dumux module, e.g. test_box1p
+ * Copy the problem, spatialparams, cc source file, input file
+
+* Adjust the CMakeLists.txt file to include your new subdirectory
+
+* Add a new CMakeLists.txt in the folder `appl` with the content
+
+```cmake
+# add a new box 1p test
+add_dumux_test(test_box1p test_box1p test_box1p.cc test_box1p)
+
+# link the input file to the build folder
+dune_symlink_to_source_files(FILES test_box1p.input)
+```
+
+* Reconfigure your module by running in the topmost directory of your new module
+
+```bash
+cmake build-cmake
+```
+
+* Build and execute the test problem
+
+```bash
+cd build-cmake
+make build_tests
+cd appl
+./test_box1p
+```
+
+<hr><br><br>
+### Task 4: Create a new GitLab project
+<hr>
+
+* Login with your username and password at https://git.iws.uni-stuttgart.de/
+
+Note: If you don't have an account create one. We allow anyone to host repositories
+on our GitLab instance as long as it is DuMuX related.
+
+* Click the **New project** button
+
+* Specify your project name and click the **Create project** button
+
+* Follow the given instructions for an *existing folder*
+
+**Important**: Before executing the `git add .` command, you should add your cmake build folder to `.gitignore`.
+The easiest way to do so is to copy the `.gitignore` file from the dumux module into your module path. If everything
+worked, executing `git status` should not show `build-cmake` anymore. Never put your executables or other build files
+under version control. Only source files (`*.hh`, `*.cc`, `*.input`, `CMakeLists.txt`) should be under version control.
\ No newline at end of file
diff --git a/exercises/exercise-fluidsystem/2p2cproblem.hh b/exercises/exercise-fluidsystem/2p2cproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..0553b406b48ca1344bdba5f9c0de925261bd7da6
--- /dev/null
+++ b/exercises/exercise-fluidsystem/2p2cproblem.hh
@@ -0,0 +1,269 @@
+// -*- 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 Tutorial problem for a fully coupled twophase box model.
+ */
+#ifndef DUMUX_EXERCISE_THREE_B_PROBLEM_HH
+#define DUMUX_EXERCISE_THREE_B_PROBLEM_HH
+
+// The numerical model
+#include <dumux/porousmediumflow/2p2c/model.hh>
+
+//The box discretization
+#include <dumux/discretization/box/properties.hh>
+
+// The base porous media box problem
+#include <dumux/porousmediumflow/problem.hh>
+
+// Spatially dependent parameters
+#include "spatialparams.hh"
+
+// The fluid system that is created in this exercise
+#include "fluidsystems/h2omycompressiblecomponent.hh"
+
+namespace Dumux{
+// Forward declaration of the problem class
+template <class TypeTag> class ExerciseThreeProblemTwoPTwoC;
+
+namespace Properties {
+// Create a new type tag for the problem
+NEW_TYPE_TAG(ExerciseThreeTwoPTwoCTypeTag, INHERITS_FROM(TwoPTwoC));
+NEW_TYPE_TAG(ExerciseThreeBoxTwoPTwoCTypeTag, INHERITS_FROM(BoxModel, ExerciseThreeTwoPTwoCTypeTag));
+
+// Set the "Problem" property
+SET_TYPE_PROP(ExerciseThreeTwoPTwoCTypeTag, Problem, ExerciseThreeProblemTwoPTwoC<TypeTag>);
+
+// Set the spatial parameters
+SET_TYPE_PROP(ExerciseThreeTwoPTwoCTypeTag, SpatialParams,
+ ExerciseThreeSpatialParams<typename GET_PROP_TYPE(TypeTag, FVGridGeometry),
+ typename GET_PROP_TYPE(TypeTag, Scalar)>);
+
+// Set grid and the grid creator to be used
+#if HAVE_DUNE_ALUGRID
+SET_TYPE_PROP(ExerciseThreeTwoPTwoCTypeTag, Grid, Dune::ALUGrid</*dim=*/2, 2, Dune::cube, Dune::nonconforming>);
+#elif HAVE_UG
+SET_TYPE_PROP(ExerciseThreeTwoPTwoCTypeTag, Grid, Dune::UGGrid<2>);
+#else
+SET_TYPE_PROP(ExerciseThreeTwoPTwoCTypeTag, Grid, Dune::YaspGrid<2>);
+#endif // HAVE_DUNE_ALUGRID
+
+ // The fluid system property
+SET_PROP(ExerciseThreeTwoPTwoCTypeTag, FluidSystem)
+{
+private:
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+public:
+ using type = FluidSystems::H2OMyCompressibleComponent<Scalar>;
+};
+
+}
+
+/*!
+ * \ingroup TwoPBoxModel
+ *
+ * \brief Tutorial problem for a fully coupled twophase box model.
+ */
+template <class TypeTag>
+class ExerciseThreeProblemTwoPTwoC : public PorousMediumFlowProblem<TypeTag>
+{
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+
+ // Grid dimension
+ enum { dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld
+ };
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ // Dumux specific types
+ 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;
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+
+public:
+ ExerciseThreeProblemTwoPTwoC(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+ : ParentType(fvGridGeometry)
+ , eps_(3e-6)
+ {
+#if !(HAVE_DUNE_ALUGRID || HAVE_UG)
+ std::cout << "If you want to use simplices instead of cubes, install and use dune-ALUGrid or UGGrid." << std::endl;
+#endif // !(HAVE_DUNE_ALUGRID || HAVE_UG)
+
+ // initialize the fluid system
+ FluidSystem::init();
+
+ // set the depth of the bottom of the reservoir
+ depthBOR_ = this->fvGridGeometry().bBoxMax()[dimWorld-1];
+
+ // name of the problem and output file
+ name_ = getParam<std::string>("Problem.Name");
+ }
+
+ /*!
+ * \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 \f$ K \f$
+ */
+ Scalar temperature() const
+ { return 283.15; }
+
+ /*!
+ * \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_ || globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_)
+ bcTypes.setAllDirichlet();
+ 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
+ {
+ PrimaryVariables priVars;
+ priVars.setState(Indices::firstPhaseOnly);
+ priVars[Indices::pressureIdx] = 200.0e3 + 9.81*1000*(depthBOR_ - globalPos[dimWorld-1]);
+ priVars[Indices::switchIdx] = 0.0; // 0 % oil saturation on left boundary
+ return priVars;
+ }
+
+ /*!
+ * \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);
+ Scalar up = this->fvGridGeometry().bBoxMax()[dimWorld-1];
+ // extraction of oil on the right boundary for approx. 1.e6 seconds
+ if (globalPos[dimWorld-1] > up - eps_ && globalPos[0] > 20 && globalPos[0] < 40) {
+ // oil outflux of 30 g/(m * s) on the right boundary.
+ // we solve for the mole balance, so we have to divide by the molar mass
+ values[Indices::conti0EqIdx + FluidSystem::H2OIdx] = 0;
+ values[Indices::conti0EqIdx + FluidSystem::NAPLIdx] = -3e-2/FluidSystem::MyCompressibleComponent::molarMass();
+ } else {
+ // no-flow on the remaining Neumann-boundaries.
+ values[Indices::conti0EqIdx + FluidSystem::H2OIdx] = 0;
+ values[Indices::conti0EqIdx + FluidSystem::NAPLIdx] = 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);
+
+ values[Indices::pressureIdx] = 200.0e3 + 9.81*1000*(depthBOR_ - globalPos[dimWorld-1]); // 200 kPa = 2 bar
+ values[Indices::switchIdx] = 0.0;
+
+ return values;
+ }
+
+ // \}
+
+ /*!
+ * \brief Returns the source term
+ *
+ * \param values Stores the source values for the conservation equations in
+ * \f$ [ \textnormal{unit of primary variable} / (m^\textrm{dim} \cdot s )] \f$
+ * \param globalPos The global position
+ */
+ PrimaryVariables sourceAtPos(const GlobalPosition &globalPos) const
+ {
+ PrimaryVariables values(0.0);
+
+ return values;
+ }
+
+
+private:
+ // small epsilon value
+ Scalar eps_;
+
+ // depth at the bottom of the reservoir
+ Scalar depthBOR_;
+ std::string name_;
+};
+}
+
+#endif
diff --git a/exercises/exercise-fluidsystem/2pproblem.hh b/exercises/exercise-fluidsystem/2pproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..fd3e3767401b226f44b95529db23603ab58c827f
--- /dev/null
+++ b/exercises/exercise-fluidsystem/2pproblem.hh
@@ -0,0 +1,308 @@
+// -*- 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 Tutorial problem for a fully coupled twophase box model.
+ */
+#ifndef DUMUX_EXERCISE_THREE_A_PROBLEM_HH
+#define DUMUX_EXERCISE_THREE_A_PROBLEM_HH
+
+// The numerical model
+#include <dumux/porousmediumflow/2p/model.hh>
+
+//The box discretization
+#include <dumux/discretization/box/properties.hh>
+
+//The grid managers
+#if HAVE_DUNE_ALUGRID
+#include <dune/alugrid/grid.hh>
+#elif HAVE_UG
+#include <dune/grid/uggrid.hh>
+#else
+#include <dune/grid/yaspgrid.hh>
+#endif
+
+// The base porous media box problem
+#include <dumux/porousmediumflow/problem.hh>
+
+// Spatially dependent parameters
+#include "spatialparams.hh"
+
+// The water component
+#include <dumux/material/components/tabulatedcomponent.hh>
+#include <dumux/material/components/h2o.hh>
+
+// The components that will be created in this exercise
+#include "components/myincompressiblecomponent.hh"
+// #include "components/mycompressiblecomponent.hh"
+
+// We will only have liquid phases here
+#include <dumux/material/fluidsystems/1pliquid.hh>
+
+// The two-phase immiscible fluid system
+#include <dumux/material/fluidsystems/2pimmiscible.hh>
+
+namespace Dumux{
+// Forward declaration of the problem class
+template <class TypeTag> class ExerciseThreeProblemTwoP;
+
+namespace Properties {
+// Create a new type tag for the problem
+NEW_TYPE_TAG(ExerciseThreeTwoPTypeTag, INHERITS_FROM(TwoP));
+NEW_TYPE_TAG(ExerciseThreeBoxTwoPTypeTag, INHERITS_FROM(BoxModel, ExerciseThreeTwoPTypeTag));
+
+// Set the "Problem" property
+SET_TYPE_PROP(ExerciseThreeTwoPTypeTag, Problem, ExerciseThreeProblemTwoP<TypeTag>);
+
+// Set the spatial parameters
+SET_TYPE_PROP(ExerciseThreeTwoPTypeTag, SpatialParams,
+ ExerciseThreeSpatialParams<typename GET_PROP_TYPE(TypeTag, FVGridGeometry),
+ typename GET_PROP_TYPE(TypeTag, Scalar)>);
+
+// Set grid and the grid creator to be used
+#if HAVE_DUNE_ALUGRID
+SET_TYPE_PROP(ExerciseThreeTwoPTypeTag, Grid, Dune::ALUGrid</*dim=*/2, 2, Dune::cube, Dune::nonconforming>);
+#elif HAVE_UG
+SET_TYPE_PROP(ExerciseThreeTwoPTypeTag, Grid, Dune::UGGrid<2>);
+#else
+SET_TYPE_PROP(ExerciseThreeTwoPTypeTag, Grid, Dune::YaspGrid<2>);
+#endif // HAVE_DUNE_ALUGRID
+
+// we use the immiscible fluid system here
+SET_PROP(ExerciseThreeTwoPTypeTag, FluidSystem)
+{
+private:
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using TabulatedH2O = Components::TabulatedComponent<Components::H2O<Scalar>>;
+ using WettingPhase = typename FluidSystems::OnePLiquid<Scalar, TabulatedH2O>;
+ /*!
+ * Uncomment first line and comment second line for using the incompressible component
+ * Uncomment second line and comment first line for using the compressible component
+ */
+ using NonWettingPhase = typename FluidSystems::OnePLiquid<Scalar, MyIncompressibleComponent<Scalar> >;
+ // using NonWettingPhase = typename FluidSystems::OnePLiquid<Scalar, MyCompressibleComponent<Scalar> >;
+
+public:
+ using type = typename FluidSystems::TwoPImmiscible<Scalar, WettingPhase, NonWettingPhase>;
+};
+
+}
+
+/*!
+ * \ingroup TwoPBoxModel
+ *
+ * \brief Tutorial problem for a fully coupled twophase box model.
+ */
+template <class TypeTag>
+class ExerciseThreeProblemTwoP : public PorousMediumFlowProblem<TypeTag>
+{
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+
+ // Grid dimension
+ enum { dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld
+ };
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ // Dumux specific types
+ 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;
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+
+ enum {
+ waterPressureIdx = Indices::pressureIdx,
+ naplSaturationIdx = Indices::saturationIdx,
+ contiWEqIdx = Indices::conti0EqIdx + FluidSystem::comp0Idx, // water transport equation index
+ contiNEqIdx = Indices::conti0EqIdx + FluidSystem::comp1Idx // napl transport equation index
+ };
+
+public:
+ ExerciseThreeProblemTwoP(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+ : ParentType(fvGridGeometry)
+ , eps_(3e-6)
+ {
+#if !(HAVE_DUNE_ALUGRID || HAVE_UG)
+ std::cout << "If you want to use simplices instead of cubes, install and use dune-ALUGrid or UGGrid." << std::endl;
+#endif // !(HAVE_DUNE_ALUGRID || HAVE_UG)
+
+ // initialize the tables for the water properties
+ std::cout << "Initializing the tables for the water properties" << std::endl;
+ Components::TabulatedComponent<Components::H2O<Scalar>>::init(/*tempMin=*/273.15,
+ /*tempMax=*/623.15,
+ /*numTemp=*/100,
+ /*pMin=*/0.0,
+ /*pMax=*/20e6,
+ /*numP=*/200);
+
+ // set the depth of the bottom of the reservoir
+ depthBOR_ = this->fvGridGeometry().bBoxMax()[dimWorld-1];
+
+ // name of the problem and output file
+ name_ = getParam<std::string>("Problem.Name");
+ }
+
+ /*!
+ * \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 \f$ K \f$
+ */
+ Scalar temperature() const
+ { return 283.15; }
+
+ /*!
+ * \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_ || globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_)
+ bcTypes.setAllDirichlet();
+ 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
+ {
+ PrimaryVariables priVars;
+
+ priVars[waterPressureIdx] = 200.0e3 + 9.81*1000*(depthBOR_ - globalPos[dimWorld-1]); // 200 kPa = 2 bar
+ priVars[naplSaturationIdx] = 0.0; // 0 % oil saturation on left boundary
+ return priVars;
+ }
+
+ /*!
+ * \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);
+
+ Scalar up = this->fvGridGeometry().bBoxMax()[dimWorld-1];
+ // extraction of oil on the right boundary for approx. 1.e6 seconds
+ if (globalPos[dimWorld-1] > up - eps_ && globalPos[0] > 20 && globalPos[0] < 40) {
+ // oil outflux of 30 g/(m * s) on the right boundary.
+ values[contiWEqIdx] = 0;
+ values[contiNEqIdx] = -3e-2;
+ } else {
+ // no-flow on the remaining Neumann-boundaries.
+ values[contiWEqIdx] = 0;
+ values[contiNEqIdx] = 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[waterPressureIdx] = 200.0e3 + 9.81*1000*(depthBOR_ - globalPos[dimWorld-1]); // 200 kPa = 2 bar (pw)
+ values[naplSaturationIdx] = 0.0; // (sn)
+
+ return values;
+ }
+ // \}
+
+ /*!
+ * \brief Returns the source term
+ *
+ * \param values Stores the source values for the conservation equations in
+ * \f$ [ \textnormal{unit of primary variable} / (m^\textrm{dim} \cdot s )] \f$
+ * \param globalPos The global position
+ */
+ PrimaryVariables sourceAtPos(const GlobalPosition &globalPos) const
+ {
+ PrimaryVariables values(0.0);
+
+ return values;
+ }
+
+private:
+ // small epsilon value
+ Scalar eps_;
+ std::string name_; //! Problem name
+
+ // depth at the bottom of the reservoir
+ Scalar depthBOR_;
+};
+}
+
+#endif
diff --git a/exercises/exercise-fluidsystem/CMakeLists.txt b/exercises/exercise-fluidsystem/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..01a4fd79d4c9b90fae2192fe33be4401f0c297b0
--- /dev/null
+++ b/exercises/exercise-fluidsystem/CMakeLists.txt
@@ -0,0 +1,20 @@
+# executables for exercise part a & b
+#part a: 2pproblem
+dune_add_test(NAME exercise3_a
+ SOURCES exercise3.cc
+ COMPILE_DEFINITIONS TYPETAG=ExerciseThreeBoxTwoPTypeTag
+ COMPILE_ONLY # for testing purposes, ignore for the exercise
+ CMD_ARGS exercise3_a.input)
+
+#part b: 2p2cproblem
+dune_add_test(NAME exercise3_b
+ SOURCES exercise3.cc
+ COMPILE_DEFINITIONS TYPETAG=ExerciseThreeBoxTwoPTwoCTypeTag
+ COMPILE_ONLY # for testing purposes, ignore for the exercise
+ CMD_ARGS exercise3_b.input)
+
+# add tutorial to the common target
+add_dependencies(test_exercises exercise3_a exercise3_b)
+
+# add symlinks for the input files
+add_input_file_links()
diff --git a/exercises/exercise-fluidsystem/README.md b/exercises/exercise-fluidsystem/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..612e88ec4a3d904c78ba4f493536925d8ebcdd53
--- /dev/null
+++ b/exercises/exercise-fluidsystem/README.md
@@ -0,0 +1,216 @@
+# Exercise #3 (DuMuX course)
+
+The aim of this exercise is to get familiar with the _DuMuX_ way of implementing new components (fluids) and fluid systems (mixtures). In the scope of this exercise, a new fictitious component is implemented (exercise _3a_) as well as its mixture with water (exercise _3b_).
+
+## Problem set-up
+
+The domain has a size of 60 x 60 m and contains two low-permeable lenses. Initially, the domain is fully water saturated and the fictitious component is injected through the middle portion of the upper boundary by means of a Neumann boundary condition. The remaining parts of the upper and the entire lower boundary are Neumann no-flow while on the two lateral sides Dirichlet boundary conditions are applied (hydrostatic conditions for the pressure and zero saturation).
+
+
+
+
+## Preparing the exercise
+
+* Navigate to the directory `dumux/tutorial/ex3`
+
+### 1. Getting familiar with the code
+
+Locate all the files you will need for this exercise
+* The shared __main file__ : `exercise3.cc`
+* The __input file__ for part a: `exercise3_a.input`
+* The __problem file__ for part a: `2pproblem.hh`
+* The __input file__ for part b: `exercise3_b.input`
+* The __problem file__ for part b: `2p2cproblem.hh`
+* The __spatial parameters file__: `spatialparams.hh`
+
+Furthermore you will find the following folders:
+* `binarycoefficients`: Stores headers containing data/methods on binary mixtures
+* `components`: Stores headers containing data/methods on pure components
+* `fluidsystems`: Stores headers containing data/methods on mixtures of pure components. Uses methods from `binarycoefficients`.
+
+To see more components, fluidsystems and binarycoefficients implementations, have a look at the folder `dumux/material`.
+
+### 2. Implement a new component
+
+In the following, the basic steps required to set the desired fluid system are outlined. Here, this is done in the __problem file__, i.e. for this part of the exercise the code shown below is taken from the `2pproblem.hh` file.
+
+In this part of the exercise we will consider a system consisting of two immiscible phases. Therefore, the _TypeTag_ for this problem (`ExerciseThreeBoxTwoPTypeTag`) derives from a base _TypeTag_ (`ExerciseThreeTwoPTypeTag`) that itself derives from the `TwoP` _TypeTag_ (immiscible two-phase model properties).
+
+```c++
+NEW_TYPE_TAG(ExerciseThreeTwoPTypeTag, INHERITS_FROM(TwoP));
+```
+
+In order to be able to derive from this _TypeTag_, the declaration of the `TwoP` _TypeTag_ has to be included. It can be found in the `2p/model.hh` header:
+
+```c++
+// The numerical model
+#include <dumux/porousmediumflow/2p/model.hh>
+```
+
+Additionally, the _TypeTag_ for this problem (`ExerciseThreeBoxTwoPTypeTag`) derives from the `BoxModel` _TypeTag_, to specify properties of the discretization scheme. For a cell-centered scheme, you could derive from `CCTpfaModel` or `CCMpfaModel` instead. Again the corresponding header has to be included
+
+```c++
+// The discretization
+#include <dumux/discretization/box/properties.hh>
+```
+
+As wetting phase we want to use water and we want to precompute tables on which the properties are then interpolated in order to save computational time. Thus, in a first step we have to include the following headers:
+
+```c++
+// The water component
+#include <dumux/material/components/tabulatedcomponent.hh>
+#include <dumux/material/components/h2o.hh>
+```
+The non-wetting phase will be our new component, where we want to implement an incompressible and a compressible variant. The respective headers are prepared, but still incomplete. The compressible variant is still commented so that compilation does not fail when finishing the incompressible variant.
+
+```c++
+// The components that will be created in this exercise
+#include "components/myincompressiblecomponent.hh"
+// #include "components/mycompressiblecomponent.hh"
+```
+As mentioned above, we want to simulate two non-mixing components. The respective fluid system is found in:
+
+```c++
+// The two-phase immiscible fluid system
+#include <dumux/material/fluidsystems/2pimmiscible.hh>
+```
+
+This fluid system expects __phases__ as input and so far we have only included the components, which contain data on the pure component for all physical states. Thus, we need to include
+
+```c++
+// We will only have liquid phases here
+#include <dumux/material/fluidsystems/1pliquid.hh>
+```
+
+which creates a _liquid phase_ from a given component. Finally, using all of the included classes we set the fluid system property by choosing that the non-wetting phase is a one-phase liquid (OnePLiquid) consisting of the incompressible fictitious component and that the wetting-phase consists of tabulated water in the immiscible fluid system:
+
+
+```c++
+// we use the immiscible fluid system here
+SET_PROP(ExerciseThreeTwoPTypeTag, FluidSystem)
+{
+private:
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using TabulatedH2O = Components::TabulatedComponent<Components::H2O<Scalar>>;
+ using WettingPhase = typename FluidSystems::OnePLiquid<Scalar, TabulatedH2O>;
+ /*!
+ * Uncomment first line and comment second line for using the incompressible component
+ * Uncomment second line and comment first line for using the compressible component
+ */
+ using NonWettingPhase = typename FluidSystems::OnePLiquid<Scalar, MyIncompressibleComponent<Scalar> >;
+ // using NonWettingPhase = typename FluidSystems::OnePLiquid<Scalar, MyCompressibleComponent<Scalar> >;
+
+public:
+ using type = typename FluidSystems::TwoPImmiscible<Scalar, WettingPhase, NonWettingPhase>;
+};
+```
+
+### 2.1. Incompressible component
+
+Open the file `myincompressiblecomponent.hh`. You can see in line 42 that a component should always derive from the _Base_ class (see `dumux/material/components/base.hh`), which defines the interface of a _DuMuX_ component with possibly required functions to be overloaded by the actual implementation. Additionally it is required for liquids to derive from the _Liquid_ class (see `dumux/material/components/liquid.hh`), for gases to derive from the _Gas_ class (see `dumux/material/components/gas.hh`) and for solids to derive from the _Solid_ class (see `dumux/material/components/solid.hh`), with functions specific to liquid, gas or solid.
+
+```c++
+/*!
+ * \ingroup Components
+ * \brief A ficitious component to be implemented in exercise 3.
+ *
+ * \tparam Scalar The type used for scalar values
+ */
+template <class Scalar>
+class MyIncompressibleComponent
+: public Components::Base<Scalar, MyIncompressibleComponent<Scalar> >
+, public Components::Liquid<Scalar, MyIncompressibleComponent<Scalar> >
+```
+
+__Task__:
+
+Implement an incompressible component into the file `myincompressiblecomponent.hh`, which has the following specifications:
+
+| Parameter | unit | value |
+| -----| --------| -------- |
+| $`M`$ | $`Kg/mol`$ | $`131.39 \cdot 10^{-3}`$ |
+| $`\rho_{liquid}`$ | $`Kg/m^3`$ | $`1460`$ |
+| $`\mu_{liquid}`$ | $`Pa \cdot s`$ | $`5.7 \cdot 10^{-4}`$ |
+
+In order to do so, have a look at the files `dumux/material/components/base.hh` and `dumux/material/components/liquid.hh` to see how the interfaces are defined and overload them accordingly.
+
+In order to execute the program, change to the build directory and compile and execute the program by typing
+
+```bash
+cd build-cmake/tutorial/ex3
+make exercise3_a
+./exercise3_a exercise3_a.input
+```
+
+The saturation distribution of the nonwetting phase at the final simulation time should look like this:
+
+
+
+### 2.2. Compressible component
+
+We now want to implement a pressure-dependent density for our component. Open the file `mycompressiblecomponent.hh` and copy in the functions you implemented for the incompressible variant. Now substitute the method that returns the density by the following expression:
+
+$`\displaystyle \rho_{MyComp} = \rho_{min} + \frac{ \rho_{max} - \rho_{min} }{ 1 + \rho_{min}*e^{-1.0*k*(\rho_{max} - \rho_{min})*p} } `$
+
+where $`p`$ is the pressure and $`\rho_{min} = 1440 `$, $`\rho_{max} = 1480 `$ and $`k = 5 \cdot 10^{-7} `$. Also, make sure the header is included in the `2pproblem.hh` file by uncommenting line 45. Furthermore, the new component has to be set as the non-wetting phase in the fluid system, i.e. comment line 90 and uncomment line 91. The non-wetting density distribution at the final simulation time should look like this:
+
+
+
+### 3. Implement a new fluid system
+
+The problem file for this part of the exercise is `2p2cproblem.hh`. We now want to implement a new fluid system consisting of two liquid phases, which are water and the previously implemented compressible component. We will consider compositional effects, which is why we now have to derive our _TypeTag_ (`ExerciseThreeBoxTwoPTwoCTypeTag`) from a _TypeTag_ (`ExerciseThreeTwoPTwoCTypeTag`) that derives from the `TwoPTwoC` model _TypeTag_:
+
+```c++
+// The numerical model
+#include <dumux/porousmediumflow/2p2c/model.hh>
+```
+
+```c++
+// Create a new type tag for the problem
+NEW_TYPE_TAG(ExerciseThreeTwoPTwoCTypeTag, INHERITS_FROM(TwoPTwoC));
+NEW_TYPE_TAG(ExerciseThreeBoxTwoPTwoCTypeTag, INHERITS_FROM(BoxModel, ExerciseThreeTwoPTwoCTypeTag));
+```
+
+The new fluid system is to be implemented in the file `fluidsystems/h2omycompressiblecomponent.hh`. This is already included in the problem and the fluid system property is set accordingly.
+
+```c++
+// The fluid system that is created in this exercise
+#include "fluidsystems/h2omycompressiblecomponent.hh"
+```
+
+```c++
+// The fluid system property
+SET_PROP(ExerciseThreeTwoPTwoCTypeTag, FluidSystem)
+{
+private:
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+public:
+ using type = FluidSystems::H2OMyCompressibleComponent<Scalar>;
+};
+```
+
+In the `fluidsystems/h2omycompressiblecomponent.hh` file, your implemented component and the binary coefficient files are already included.
+
+```c++
+// the ficitious component that was created in exercise 3a
+#include <tutorial/ex3/components/mycompressiblecomponent.hh>
+
+// the binary coefficients corresponding to this fluid system
+#include <tutorial/ex3/binarycoefficients/h2omycompressiblecomponent.hh>
+```
+
+__Task__:
+
+Under the assumption that one molecule of `MyCompressibleComponent` displaces exactly one molecule of water, the water phase density can be expressed as follows:
+
+$` \rho_{w} = \frac{ \rho_{w, pure} }{ M_{H_2O} }*(M_{H_2O}*x_{H_2O} + M_{MyComponent}*x_{MyComponent}) `$
+
+Implement this dependency in the `density()` method in the fluid system. In order to compile and execute the program run
+
+```bash
+cd build-cmake/tutorial/ex3
+make exercise3_b
+./exercise3_b exercise3_b.input
+```
+
+You will observe an error message and an abortion of the program. This is due to the fact that in order for the constraint solver and other mechanisms in the two-phase two-component model to work, an additional functionality in the component has to be implemented. The model has to know the vapour pressure. As in the previous exercise, check the `dumux/material/components/base.hh` file for this function and implement it into `mycompressiblecomponent.hh`. For the vapour pressure, use a value of $`3900`$ Pa.
diff --git a/exercises/exercise-fluidsystem/binarycoefficients/h2omycompressiblecomponent.hh b/exercises/exercise-fluidsystem/binarycoefficients/h2omycompressiblecomponent.hh
new file mode 100644
index 0000000000000000000000000000000000000000..843a69a675bc1b8095a850c94f07a8ca65232228
--- /dev/null
+++ b/exercises/exercise-fluidsystem/binarycoefficients/h2omycompressiblecomponent.hh
@@ -0,0 +1,74 @@
+// -*- 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 Binary coefficients for water and a ficticious component implemented in tutorial exercise 3a.
+ */
+#ifndef DUMUX_BINARY_COEFF_H2O_MYCOMPRESSIBLECOMPONENT_HH
+#define DUMUX_BINARY_COEFF_H2O_MYCOMPRESSIBLECOMPONENT_HH
+
+namespace Dumux
+{
+namespace BinaryCoeff
+{
+
+/*!
+ * \brief Binary coefficients for water and a ficticious component implemented in tutorial exercise 3a
+ * The implementation of the missing methods in this file is part of exercise 3b.
+ */
+class H2O_MyCompressibleComponent
+{
+public:
+ /*!
+ * \brief Henry coefficient \f$[N/m^2]\f$ for the fictitous component in liquid water.
+ */
+ template <class Scalar>
+ static Scalar henryMyCompressibleComponentInWater(Scalar temperature)
+ {
+ Scalar dumuxH = 1.5e-1 / 101.325; // unit [(mol/m^3)/Pa]
+ dumuxH *= 18.02e-6; //multiplied by molar volume of reference phase = water
+ return 1.0/dumuxH; // [Pa]
+ }
+
+ /*!
+ * \brief Henry coefficient \f$[N/m^2]\f$ for water in the ficticious component.
+ */
+ template <class Scalar>
+ static Scalar henryWaterInMyCompressibleComponent(Scalar temperature)
+ {
+ // arbitrary
+ return 1.0e8; // [Pa]
+ }
+
+ /*!
+ * \brief Diffusion coefficient [m^2/s] for my ficticious component in liquid water or vice versa.
+ */
+ template <class Scalar>
+ static Scalar liquidDiffCoeff(Scalar temperature, Scalar pressure)
+ {
+ // arbitrary
+ return 1.e-9;
+ }
+};
+
+}
+} // end namespace
+
+#endif
diff --git a/exercises/exercise-fluidsystem/components/mycompressiblecomponent.hh b/exercises/exercise-fluidsystem/components/mycompressiblecomponent.hh
new file mode 100644
index 0000000000000000000000000000000000000000..0d18c90a0c309f5bc7c395fe516f03e853e81e6f
--- /dev/null
+++ b/exercises/exercise-fluidsystem/components/mycompressiblecomponent.hh
@@ -0,0 +1,126 @@
+// -*- 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
+ * \ingroup Components
+ * \brief A fictitious component to be implemented in exercise 3.
+ */
+#ifndef DUMUX_MYCOMPRESSIBLECOMPONENT_HH
+#define DUMUX_MYCOMPRESSIBLECOMPONENT_HH
+
+#include <dumux/material/idealgas.hh>
+
+#include <dumux/material/components/base.hh>
+#include <dumux/material/components/liquid.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup Components
+ * \brief A fictitious component to be implemented in exercise 3.
+ *
+ * \tparam Scalar The type used for scalar values
+ */
+template <class Scalar>
+class MyCompressibleComponent
+: public Components::Base<Scalar, MyCompressibleComponent<Scalar> >
+, public Components::Liquid<Scalar, MyCompressibleComponent<Scalar> >
+{
+
+public:
+ /*!
+ * \brief A human readable name for MyCompressibleComponent.
+ */
+ static std::string name()
+ { return "MyCompressibleComponent"; }
+
+ /*!
+ * TODO: Copy the methods implemented in MyIncompressibleComponent and substitute
+ * the density calculation by the expression given in the exercise description.
+ */
+
+ /*!
+ * \brief Returns true if the liquid phase is assumed to be compressible
+ */
+ static constexpr bool liquidIsCompressible()
+ { return true; }
+
+ /*!
+ * \brief The molar mass in \f$\mathrm{[kg/mol]}\f$ of the component.
+ */
+ static Scalar molarMass()
+ {
+ // TODO: replace the line below by a meaningful return statement
+ DUNE_THROW(Dune::NotImplemented, "Todo: implement molar mass");
+ }
+
+ /*!
+ * \brief The density \f$\mathrm{[kg/m^3]}\f$ of the liquid component at a given pressure in
+ * \f$\mathrm{[Pa]}\f$ and temperature in \f$\mathrm{[K]}\f$.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static Scalar liquidDensity(Scalar temperature, Scalar pressure)
+ {
+ // TODO: replace the line below by a meaningful return statement
+ DUNE_THROW(Dune::NotImplemented, "Todo: implement liquid density");
+ }
+
+ /*!
+ * \brief The molar density of MyCompressibleComponent in \f$\mathrm{[mol/m^3]}\f$ at a given pressure and temperature.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ *
+ */
+ static Scalar liquidMolarDensity(Scalar temperature, Scalar pressure)
+ {
+ // TODO: replace the line below by a meaningful return statement
+ DUNE_THROW(Dune::NotImplemented, "Todo: implement liquid molar density");
+ }
+
+ /*!
+ * \brief The dynamic liquid viscosity \f$\mathrm{[Pa*s]}\f$ of the pure component.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static Scalar liquidViscosity(Scalar temperature, Scalar pressure)
+ {
+ // TODO: replace the line below by a meaningful return statement
+ DUNE_THROW(Dune::NotImplemented, "Todo: implement liquid viscosity");
+ }
+
+ /*!
+ * \brief The vapor pressure in \f$\mathrm{[Pa]}\f$ of the component at a given
+ * temperature in \f$\mathrm{[K]}\f$.
+ *
+ * \param T temperature of the component in \f$\mathrm{[K]}\f$
+ */
+ static Scalar vaporPressure(Scalar t)
+ {
+ // TODO: replace the line below by a meaningful return statement
+ DUNE_THROW(Dune::NotImplemented, "Todo: implement vapour pressure");
+ }
+};
+
+} // end namespace
+
+#endif
diff --git a/exercises/exercise-fluidsystem/components/myincompressiblecomponent.hh b/exercises/exercise-fluidsystem/components/myincompressiblecomponent.hh
new file mode 100644
index 0000000000000000000000000000000000000000..dae8fdde10f9fc72ac6108c2902877f0ccd0bb09
--- /dev/null
+++ b/exercises/exercise-fluidsystem/components/myincompressiblecomponent.hh
@@ -0,0 +1,112 @@
+// -*- 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
+ * \ingroup Components
+ * \brief A fictitious component to be implemented in tutorial exercise 3.
+ */
+#ifndef DUMUX_MYINCOMPRESSIBLECOMPONENT_HH
+#define DUMUX_MYINCOMPRESSIBLECOMPONENT_HH
+
+#include <dumux/material/idealgas.hh>
+
+#include <dumux/material/components/base.hh>
+#include <dumux/material/components/liquid.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup Components
+ * \brief A fictitious component to be implemented in exercise 3.
+ *
+ * \tparam Scalar The type used for scalar values
+ */
+template <class Scalar>
+class MyIncompressibleComponent
+: public Components::Base<Scalar, MyIncompressibleComponent<Scalar> >
+, public Components::Liquid<Scalar, MyIncompressibleComponent<Scalar> >
+{
+public:
+ /*!
+ * \brief A human readable name for MyIncompressibleComponent.
+ */
+ static std::string name()
+ { return "MyIncompressibleComponent"; }
+
+ /*!
+ * TODO: Implement the methods for the component data given in the exercise description.
+ */
+
+ /*!
+ * \brief Returns true if the liquid phase is assumed to be compressible
+ */
+ static constexpr bool liquidIsCompressible()
+ { return false; }
+
+ /*!
+ * \brief The molar mass in \f$\mathrm{[kg/mol]}\f$ of the component.
+ */
+ static Scalar molarMass()
+ {
+ // TODO: replace the line below by a meaningful return statement
+ DUNE_THROW(Dune::NotImplemented, "Todo: implement molarMass()");
+ }
+
+ /*!
+ * \brief The density \f$\mathrm{[kg/m^3]}\f$ of the liquid component at a given pressure in
+ * \f$\mathrm{[Pa]}\f$ and temperature in \f$\mathrm{[K]}\f$.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static Scalar liquidDensity(Scalar temperature, Scalar pressure)
+ {
+ // TODO: replace the line below by a meaningful return statement
+ DUNE_THROW(Dune::NotImplemented, "Todo: implement liquidDensity()");
+ }
+
+ /*!
+ * \brief The molar density of MyIncompressibleComponent in \f$\mathrm{[mol/m^3]}\f$ at a given pressure and temperature.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ *
+ */
+ static Scalar liquidMolarDensity(Scalar temperature, Scalar pressure)
+ {
+ // TODO: replace the line below by a meaningful return statement
+ DUNE_THROW(Dune::NotImplemented, "Todo: implement liquidMolarDensity()");
+ }
+
+ /*!
+ * \brief The dynamic liquid viscosity \f$\mathrm{[Pa*s]}\f$ of the pure component.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static Scalar liquidViscosity(Scalar temperature, Scalar pressure)
+ {
+ // TODO: replace the line below by a meaningful return statement
+ DUNE_THROW(Dune::NotImplemented, "Todo: implement liquidViscosity()");
+ }
+};
+
+} // end namespace
+
+#endif
diff --git a/exercises/exercise-fluidsystem/components/plotdensityfunction.py b/exercises/exercise-fluidsystem/components/plotdensityfunction.py
new file mode 100644
index 0000000000000000000000000000000000000000..2b256689339b8c094d4bd3afb235be2b57fbb214
--- /dev/null
+++ b/exercises/exercise-fluidsystem/components/plotdensityfunction.py
@@ -0,0 +1,19 @@
+#!usr/bin/env python
+import numpy as np
+import matplotlib.pyplot as plt
+
+# function to calculate rho dependent on pressure
+rho_min = 1440;
+rho_max = 1480;
+k = 5e-7;
+
+def rho(p):
+ return rho_min + (rho_max - rho_min)/(1 + rho_min*np.exp(-1.0*k*(rho_max - rho_min)*p));
+
+# sample pressure in range (1e4, 1e7) and compute corresponding densities
+p = np.logspace(4, 7, 100)
+r = rho(p)
+
+# plot density vs. pressure
+plt.semilogx(p, r)
+plt.show()
diff --git a/exercises/exercise-fluidsystem/exercise3.cc b/exercises/exercise-fluidsystem/exercise3.cc
new file mode 100644
index 0000000000000000000000000000000000000000..acfc230b5ff53cb4de42a955a3bdba04799d5a1e
--- /dev/null
+++ b/exercises/exercise-fluidsystem/exercise3.cc
@@ -0,0 +1,236 @@
+// -*- 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 DOC ME!
+ */
+#include <config.h>
+
+#include <ctime>
+#include <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#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 "2pproblem.hh"
+#include "2p2cproblem.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/newtonsolver.hh>
+
+#include <dumux/assembly/fvassembler.hh>
+#include <dumux/assembly/diffmethod.hh>
+
+#include <dumux/discretization/methods.hh>
+
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/grid/gridmanager.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) //TODO check usage
+{
+ if (errorMsg.size() > 0) {
+ std::string errorMessageOut = "\nUsage: ";
+ errorMessageOut += progName;
+ errorMessageOut += " [options]\n";
+ errorMessageOut += errorMsg;
+ errorMessageOut += "\n\nThe list of mandatory options for this program is:\n"
+ "\t-TimeManager.TEnd End of the simulation [s] \n"
+ "\t-TimeManager.DtInitial Initial timestep size [s] \n"
+ "\t-Grid.File Name of the file containing the grid \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-Problem.Name String for naming of the output files \n"
+ "\n";
+ std::cout << errorMessageOut
+ << "\n";
+ }
+}
+
+////////////////////////
+// the main function
+////////////////////////
+int main(int argc, char** argv) try
+{
+ using namespace Dumux;
+
+ // define the type tag for this problem
+ // TYPETAG is set in CMakeLists.txt as compile time definition
+ // alternatively you could write `using TypeTag = TTAG(ExerciseThreeBoxTwoPTypeTag);`
+ // then, for the 2p2c problem you have to change this line to `using TypeTag = TTAG(ExerciseThreeBoxTwoPTypeTag);`
+ // and recompile the executable
+ using TypeTag = TTAG(TYPETAG);
+
+ // initialize MPI, finalize is done automatically on exit
+ const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+ // print dumux start message
+ if (mpiHelper.rank() == 0)
+ DumuxMessage::print(/*firstCall=*/true);
+
+ // parse command line arguments and input file
+ 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;
+ gridManager.init();
+
+ ////////////////////////////////////////////////////////////
+ // run instationary non-linear problem on this grid
+ ////////////////////////////////////////////////////////////
+
+ // we compute on the leaf grid view
+ const auto& leafGridView = gridManager.grid().leafGridView();
+
+ // create the finite volume grid geometry
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
+ fvGridGeometry->update();
+
+ // the problem (initial and boundary conditions)
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ auto problem = std::make_shared<Problem>(fvGridGeometry);
+
+ // the solution vector
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ SolutionVector x(fvGridGeometry->numDofs());
+ problem->applyInitialSolution(x);
+ auto xOld = x;
+
+ // the grid variables
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+ auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
+ gridVariables->init(x, xOld);
+
+ // get some time loop parameters
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+ const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+ auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+
+ // check if we are about to restart a previously interrupted simulation
+ Scalar restartTime = 0;
+ if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart"))
+ restartTime = getParam<Scalar>("TimeLoop.Restart");
+
+ // 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
+
+ // instantiate time loop
+ auto timeLoop = std::make_shared<TimeLoop<Scalar>>(restartTime, dt, tEnd);
+ timeLoop->setMaxTimeStepSize(maxDt);
+
+ // the assembler with time loop for instationary problem
+ using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
+ auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop);
+
+ // the linear solver
+ using LinearSolver = AMGBackend<TypeTag>;
+ auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->dofMapper());
+
+ // the non-linear solver
+ using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
+ NewtonSolver nonLinearSolver(assembler, linearSolver);
+
+ // time loop
+ timeLoop->start(); do
+ {
+ // set previous solution for storage evaluations
+ assembler->setPreviousSolution(xOld);
+
+ // solve the non-linear system with time step control
+ nonLinearSolver.solve(x, *timeLoop);
+
+ // make the new solution the old solution
+ xOld = x;
+ gridVariables->advanceTimeStep();
+
+ // advance to the time loop to the next step
+ timeLoop->advanceTimeStep();
+
+ // report statistics of this time step
+ timeLoop->reportTimeStep();
+
+ // set new dt as suggested by the newton solver
+ timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
+
+ // write solution and grid to vtk
+ vtkWriter.write(timeLoop->time());
+
+ } while (!timeLoop->finished());
+
+ timeLoop->finalize(leafGridView.comm());
+
+ ////////////////////////////////////////////////////////////
+ // finalize, print dumux message to say goodbye
+ ////////////////////////////////////////////////////////////
+
+ // print dumux end message
+ if (mpiHelper.rank() == 0)
+ {
+ Parameters::print();
+ DumuxMessage::print(/*firstCall=*/false);
+ }
+
+ return 0;
+} // end main
+catch (Dumux::ParameterException &e)
+{
+ std::cerr << std::endl << e << " ---> Abort!" << std::endl;
+ return 1;
+}
+catch (Dune::DGFException & e)
+{
+ std::cerr << "DGF exception thrown (" << e <<
+ "). Most likely, the DGF file name is wrong "
+ "or the DGF file is corrupted, "
+ "e.g. missing hash at end of file or wrong number (dimensions) of entries."
+ << " ---> Abort!" << std::endl;
+ return 2;
+}
+catch (Dune::Exception &e)
+{
+ std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
+ return 3;
+}
+catch (...)
+{
+ std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
+ return 4;
+}
diff --git a/exercises/exercise-fluidsystem/exercise3_a.input b/exercises/exercise-fluidsystem/exercise3_a.input
new file mode 100644
index 0000000000000000000000000000000000000000..4aa7757254eed862114f79834ac813398bd72f39
--- /dev/null
+++ b/exercises/exercise-fluidsystem/exercise3_a.input
@@ -0,0 +1,10 @@
+[TimeLoop]
+TEnd = 20000 # duration of the simulation [s]
+DtInitial = 10 # initial time step size [s]
+
+[Problem]
+Name = exercise3_a
+
+[Grid]
+UpperRight = 60 60 # x-/y-coordinates of the upper-right corner of the grid [m]
+Cells = 60 60 # x-/y-resolution of the grid
diff --git a/exercises/exercise-fluidsystem/exercise3_b.input b/exercises/exercise-fluidsystem/exercise3_b.input
new file mode 100644
index 0000000000000000000000000000000000000000..1e8340de2b34b3774ccf8c276587043d5cbe1bd6
--- /dev/null
+++ b/exercises/exercise-fluidsystem/exercise3_b.input
@@ -0,0 +1,10 @@
+[TimeLoop]
+TEnd = 20000 # duration of the simulation [s]
+DtInitial = 10 # initial time step size [s]
+
+[Problem]
+Name = exercise3_b
+
+[Grid]
+UpperRight = 60 60 # x-/y-coordinates of the upper-right corner of the grid [m]
+Cells = 60 60 # x-/y-resolution of the grid
diff --git a/exercises/exercise-fluidsystem/fluidsystems/h2omycompressiblecomponent.hh b/exercises/exercise-fluidsystem/fluidsystems/h2omycompressiblecomponent.hh
new file mode 100644
index 0000000000000000000000000000000000000000..0fd412e1be08535c3329e2c39b1dfa934ae19ff9
--- /dev/null
+++ b/exercises/exercise-fluidsystem/fluidsystems/h2omycompressiblecomponent.hh
@@ -0,0 +1,485 @@
+// -*- 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 A fluid system with water and a fictitious component, which is to be
+ * implemented in tutorial exercise 3a, as phases and components. This
+ * fluid system is to be implemented in exercise 3b.
+ */
+#ifndef DUMUX_H2O_MYCOMPRESSIBLECOMPONENT_FLUID_SYSTEM_HH
+#define DUMUX_H2O_MYCOMPRESSIBLECOMPONENT_FLUID_SYSTEM_HH
+
+#include <dumux/material/components/h2o.hh>
+#include <dumux/material/components/tabulatedcomponent.hh>
+
+#include <dumux/material/fluidsystems/base.hh>
+
+// the fictitious component that was created in exercise 3a
+#include <exercises/exercise-fluidsystem/components/mycompressiblecomponent.hh>
+
+// the binary coefficients corresponding to this fluid system
+#include <exercises/exercise-fluidsystem/binarycoefficients/h2omycompressiblecomponent.hh>
+
+namespace Dumux {
+namespace FluidSystems {
+
+/*!
+ * \brief A compositional fluid consisting of two liquid phases,
+ * which are water and a fictitious component from tutorial exercise 3a.
+ */
+template <class Scalar,
+ class H2OType = Dumux::Components::TabulatedComponent<Dumux::Components::H2O<Scalar> > >
+class H2OMyCompressibleComponent
+ : public BaseFluidSystem< Scalar, H2OMyCompressibleComponent<Scalar, H2OType> >
+{
+ typedef H2OMyCompressibleComponent<Scalar, H2OType> ThisType;
+ typedef BaseFluidSystem<Scalar, ThisType> Base;
+
+public:
+ typedef Dumux::MyCompressibleComponent<Scalar> MyCompressibleComponent;
+ typedef H2OType H2O;
+
+ static constexpr int numPhases = 2;
+ static constexpr int numComponents = 2;
+
+ static constexpr int phase0Idx = 0; // index of the first phase
+ static constexpr int phase1Idx = 1; // index of the second phase
+
+ static constexpr int H2OIdx = 0;
+ static constexpr int NAPLIdx = 1;
+ // export component indices to indicate the main component
+ // of the corresponding phase at atmospheric pressure 1 bar
+ // and room temperature 20°C:
+ static constexpr int comp0Idx = H2OIdx;
+ static constexpr int comp1Idx = NAPLIdx;
+
+ /*!
+ * \brief Initialize the fluid system's static parameters generically
+ *
+ * If a tabulated H2O component is used, we do our best to create
+ * tables that always work.
+ */
+ static void init()
+ {
+ init(/*tempMin=*/273.15,
+ /*tempMax=*/623.15,
+ /*numTemp=*/100,
+ /*pMin=*/0.0,
+ /*pMax=*/20e6,
+ /*numP=*/200);
+ }
+
+ /*!
+ * \brief Initialize the fluid system's static parameters using
+ * problem specific temperature and pressure ranges
+ *
+ * \param tempMin The minimum temperature used for tabulation of water [K]
+ * \param tempMax The maximum temperature used for tabulation of water [K]
+ * \param nTemp The number of ticks on the temperature axis of the table of water
+ * \param pressMin The minimum pressure used for tabulation of water [Pa]
+ * \param pressMax The maximum pressure used for tabulation of water [Pa]
+ * \param nPress The number of ticks on the pressure axis of the table of water
+ */
+ static void init(Scalar tempMin, Scalar tempMax, unsigned nTemp,
+ Scalar pressMin, Scalar pressMax, unsigned nPress)
+ {
+ if (H2O::isTabulated) {
+ std::cout << "Initializing tables for the H2O fluid properties ("
+ << nTemp*nPress
+ << " entries).\n";
+
+ H2O::init(tempMin, tempMax, nTemp,
+ pressMin, pressMax, nPress);
+ }
+ }
+
+
+ /*!
+ * \brief Return whether a phase is liquid
+ *
+ * \param phaseIdx The index of the fluid phase to consider
+ */
+ static constexpr bool isLiquid(int phaseIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+ return true;
+ }
+
+ static constexpr bool isIdealGas(int phaseIdx)
+ { return H2O::gasIsIdeal() && MyCompressibleComponent::gasIsIdeal(); }
+
+ /*!
+ * \brief Returns true if and only if a fluid phase is assumed to
+ * be an ideal mixture.
+ *
+ * We define an ideal mixture as a fluid phase where the fugacity
+ * coefficients of all components times the pressure of the phase
+ * are indepent on the fluid composition. This assumtion is true
+ * if Henry's law and Raoult's law apply. If you are unsure what
+ * this function should return, it is safe to return false. The
+ * only damage done will be (slightly) increased computation times
+ * in some cases.
+ *
+ * \param phaseIdx The index of the fluid phase to consider
+ */
+ static constexpr bool isIdealMixture(int phaseIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+ return true;
+ }
+
+ /*!
+ * \brief Returns true if and only if a fluid phase is assumed to
+ * be compressible.
+ *
+ * Compressible means that the partial derivative of the density
+ * to the fluid pressure is always larger than zero.
+ *
+ * \param phaseIdx The index of the fluid phase to consider
+ */
+ static constexpr bool isCompressible(int phaseIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+ if (phaseIdx == phase0Idx)
+ // the water component decides for the water phase...
+ return H2O::liquidIsCompressible();
+
+ // the NAPL component decides for the napl phase...
+ return MyCompressibleComponent::liquidIsCompressible();
+ }
+
+ /*!
+ * \brief Return the human readable name of a phase (used in indices)
+ */
+ static std::string phaseName(int phaseIdx)
+ {
+ switch (phaseIdx) {
+ case phase0Idx: return "w";
+ case phase1Idx: return "n";
+ };
+ DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
+ }
+
+ /*!
+ * \brief Return the human readable name of a component (used in indices)
+ */
+ static std::string componentName(int compIdx)
+ {
+ switch (compIdx) {
+ case H2OIdx: return H2O::name();
+ case NAPLIdx: return MyCompressibleComponent::name();
+ };
+ DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
+ }
+
+ /*!
+ * \brief Return the molar mass of a component in [kg/mol].
+ */
+ static Scalar molarMass(int compIdx)
+ {
+ switch (compIdx) {
+ case H2OIdx: return H2O::molarMass();
+ case NAPLIdx: return MyCompressibleComponent::molarMass();
+ };
+ DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
+ }
+
+ /*!
+ * \brief Given all mole fractions in a phase, return the phase
+ * density [kg/m^3].
+ */
+ using Base::density;
+ template <class FluidState>
+ static Scalar density(const FluidState &fluidState, int phaseIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+ if (phaseIdx == phase0Idx) {
+ // See: doctoral thesis of Steffen Ochs 2007
+ // Steam injection into saturated porous media : process analysis including experimental and numerical investigations
+ // http://elib.uni-stuttgart.de/bitstream/11682/271/1/Diss_Ochs_OPUS.pdf
+
+ // Scalar rholH2O = H2O::liquidDensity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx));
+ // Scalar clH2O = rholH2O/H2O::molarMass();
+ // Scalar x_H2O = fluidState.moleFraction (phase0Idx, H2OIdx);
+ // Scalar x_myComp = fluidState.moleFraction (phase0Idx, NAPLIdx);
+
+ /*!
+ * TODO: implement the composition-dependent water density from the exercise sheet.
+ */
+ return /*???*/1000.0;
+ }
+ else {
+ // assume the density of the fictious component to be independent of the composition
+ Scalar pressure = MyCompressibleComponent::liquidIsCompressible()?fluidState.pressure(phaseIdx):1e100;
+ return MyCompressibleComponent::liquidDensity(fluidState.temperature(phaseIdx), pressure);
+ }
+ }
+
+ using Base::molarDensity;
+ /*!
+ * \brief The molar density \f$\rho_{mol,\alpha}\f$
+ * of a fluid phase \f$\alpha\f$ in \f$\mathrm{[mol/m^3]}\f$
+ *
+ * The molar density for the simple relation is defined by the
+ * mass density \f$\rho_\alpha\f$ and the molar mass of the main component
+ *
+ * The molar density for the complrex relation is defined by the
+ * mass density \f$\rho_\alpha\f$ and the mean molar mass \f$\overline M_\alpha\f$:
+ *
+ * \f[\rho_{mol,\alpha} = \frac{\rho_\alpha}{\overline M_\alpha} \;.\f]
+ */
+ template <class FluidState>
+ static Scalar molarDensity(const FluidState &fluidState, int phaseIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+
+ Scalar T = fluidState.temperature(phaseIdx);
+ Scalar p = fluidState.pressure(phaseIdx);
+
+ // liquid phase
+ if (phaseIdx == phase0Idx)
+ {
+ // assume pure water or that each gas molecule displaces exactly one
+ // molecule in the liquid.
+ return H2O::liquidMolarDensity(T, p);
+ }
+ else
+ {
+ // assume the molar density of the fictious component to be independent of the composition
+ Scalar pressure = MyCompressibleComponent::liquidIsCompressible()?fluidState.pressure(phaseIdx):1e100;
+ return MyCompressibleComponent::liquidMolarDensity(fluidState.temperature(phaseIdx), pressure);
+ }
+ }
+
+ /*!
+ * \brief Return the viscosity of a phase.
+ */
+ using Base::viscosity;
+ template <class FluidState>
+ static Scalar viscosity(const FluidState &fluidState,
+ int phaseIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+ if (phaseIdx == phase0Idx) {
+ // assume pure water viscosity
+ return H2O::liquidViscosity(fluidState.temperature(phaseIdx),
+ fluidState.pressure(phaseIdx));
+ }
+ else {
+ // assume pure NAPL viscosity
+ return MyCompressibleComponent::liquidViscosity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx));
+ }
+ }
+
+ using Base::diffusionCoefficient;
+ template <class FluidState>
+ static Scalar diffusionCoefficient(const FluidState &fluidState, int phaseIdx, int compIdx)
+ {
+ DUNE_THROW(Dune::NotImplemented, "Diffusion coefficients");
+ }
+
+ /*!
+ * \brief Given a phase's composition, temperature and pressure,
+ * return the binary diffusion coefficient \f$\mathrm{[m^2/s]}\f$ for components
+ * \f$\mathrm{i}\f$ and \f$\mathrm{j}\f$ in this phase.
+ * \param fluidState The fluid state
+ * \param paramCache mutable parameters
+ * \param phaseIdx Index of the fluid phase
+ * \param compIIdx Index of the component i
+ * \param compJIdx Index of the component j
+ */
+ using Base::binaryDiffusionCoefficient;
+ template <class FluidState>
+ static Scalar binaryDiffusionCoefficient(const FluidState &fluidState,
+ int phaseIdx,
+ int compIIdx,
+ int compJIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+ assert(0 <= compIIdx && compIIdx < numComponents);
+ assert(0 <= compJIdx && compJIdx < numComponents);
+
+ const Scalar T = fluidState.temperature(phaseIdx);
+ const Scalar p = fluidState.pressure(phaseIdx);
+
+ // we assume the diffusion coefficient to be the same in both phases
+ return Dumux::BinaryCoeff::H2O_MyCompressibleComponent::liquidDiffCoeff(T, p);
+ }
+
+ /* Henry coefficients
+ */
+ template <class FluidState>
+ static Scalar henryCoefficient(const FluidState &fluidState,
+ int phaseIdx,
+ int compIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+ assert(0 <= compIdx && compIdx < numComponents);
+
+ const Scalar T = fluidState.temperature(phaseIdx);
+ const Scalar p = fluidState.pressure(phaseIdx);
+
+ if (compIdx == NAPLIdx && phaseIdx == phase0Idx)
+ return Dumux::BinaryCoeff::H2O_MyCompressibleComponent::henryMyCompressibleComponentInWater(T)/p;
+
+ else if (phaseIdx == phase1Idx && compIdx == H2OIdx)
+ return Dumux::BinaryCoeff::H2O_MyCompressibleComponent::henryWaterInMyCompressibleComponent(T)/p;
+
+ else
+ DUNE_THROW(Dune::InvalidStateException, "non-existent henry coefficient for phase index " << phaseIdx
+ << " and component index " << compIdx);
+ }
+
+ using Base::fugacityCoefficient;
+ /*!
+ * \brief Returns the fugacity coefficient \f$\mathrm{[-]}\f$ of a component in a
+ * phase.
+ *
+ * In this case, things are actually pretty simple. We have an ideal
+ * solution. Thus, the fugacity coefficient is 1 in the gas phase
+ * (fugacity equals the partial pressure of the component in the gas phase
+ * respectively in the liquid phases it is the inverse of the
+ * Henry coefficients scaled by pressure
+ * \param fluidState The fluid state
+ * \param phaseIdx The index of the phase
+ * \param compIdx The index of the component
+ */
+ template <class FluidState>
+ static Scalar fugacityCoefficient(const FluidState &fluidState,
+ int phaseIdx,
+ int compIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+ assert(0 <= compIdx && compIdx < numComponents);
+
+ Scalar T = fluidState.temperature(phaseIdx);
+ Scalar p = fluidState.pressure(phaseIdx);
+
+ if (phaseIdx == phase0Idx) {
+ if (compIdx == H2OIdx)
+ return H2O::vaporPressure(T)/p;
+ else if (compIdx == NAPLIdx)
+ return Dumux::BinaryCoeff::H2O_MyCompressibleComponent::henryMyCompressibleComponentInWater(T)/p;
+ }
+
+ // for the NAPL phase, we assume currently that nothing is
+ // dissolved. this means that the affinity of the NAPL
+ // component to the NAPL phase is much higher than for the
+ // other components, i.e. the fugacity coefficient is much
+ // smaller.
+ Scalar phiNapl = MyCompressibleComponent::vaporPressure(T)/p;
+ if (compIdx == NAPLIdx)
+ return phiNapl;
+ else
+ return 1e6*phiNapl;
+ }
+
+ template <class FluidState>
+ static Scalar kelvinVaporPressure(const FluidState &fluidState,
+ const int phaseIdx,
+ const int compIdx)
+ {
+ DUNE_THROW(Dune::NotImplemented, "FluidSystems::H2OMyCompressibleComponent::kelvinVaporPressure()");
+ }
+
+ /* partial pressures in the gas phase, taken from saturation vapor pressures
+ */
+ template <class FluidState>
+ static Scalar partialPressureGas(const FluidState &fluidState, int phaseIdx,
+ int compIdx)
+ {
+ assert(0 <= compIdx && compIdx < numComponents);
+
+ const Scalar T = fluidState.temperature(phaseIdx);
+ if (compIdx == NAPLIdx)
+ return MyCompressibleComponent::vaporPressure(T);
+ else if (compIdx == H2OIdx)
+ return H2O::vaporPressure(T);
+ else
+ DUNE_THROW(Dune::InvalidStateException, "non-existent component index " << compIdx);
+ }
+
+ /* inverse vapor pressures, taken from inverse saturation vapor pressures
+ */
+ template <class FluidState>
+ static Scalar inverseVaporPressureCurve(const FluidState &fluidState,
+ int phaseIdx,
+ int compIdx)
+ {
+ assert(0 <= compIdx && compIdx < numComponents);
+
+ const Scalar pressure = fluidState.pressure(phaseIdx);
+ if (compIdx == NAPLIdx)
+ return MyCompressibleComponent::vaporTemperature(pressure);
+ else if (compIdx == H2OIdx)
+ return H2O::vaporTemperature(pressure);
+ else
+ DUNE_THROW(Dune::InvalidStateException, "non-existent component index " << compIdx);
+ }
+
+
+
+ /*!
+ * \brief Given all mole fractions in a phase, return the specific
+ * phase enthalpy [J/kg].
+ */
+ using Base::enthalpy;
+ template <class FluidState>
+ static Scalar enthalpy(const FluidState &fluidState,
+ int phaseIdx)
+ {
+ DUNE_THROW(Dune::NotImplemented, "FluidSystems::H2OMyCompressibleComponent::enthalpy()");
+ }
+
+ using Base::heatCapacity;
+ template <class FluidState>
+ static Scalar heatCapacity(const FluidState &fluidState,
+ int phaseIdx)
+ {
+ DUNE_THROW(Dune::NotImplemented, "FluidSystems::H2OMyCompressibleComponent::heatCapacity()");
+ }
+
+ using Base::thermalConductivity;
+ template <class FluidState>
+ static Scalar thermalConductivity(const FluidState &fluidState,
+ int phaseIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+
+ const Scalar temperature = fluidState.temperature(phaseIdx) ;
+ const Scalar pressure = fluidState.pressure(phaseIdx);
+ if (phaseIdx == phase0Idx)
+ {
+ return H2O::liquidThermalConductivity(temperature, pressure);
+ }
+ else
+ {
+ return MyCompressibleComponent::liquidThermalConductivity(temperature, pressure);
+ }
+ DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
+ }
+
+private:
+
+};
+} // end namespace FluidSystems
+} // end namespace Dumux
+
+#endif
diff --git a/exercises/exercise-fluidsystem/spatialparams.hh b/exercises/exercise-fluidsystem/spatialparams.hh
new file mode 100644
index 0000000000000000000000000000000000000000..36297290e39597bf906f30495a2b6b181a128ae6
--- /dev/null
+++ b/exercises/exercise-fluidsystem/spatialparams.hh
@@ -0,0 +1,163 @@
+// -*- 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 The spatial parameters for the fully coupled tutorial problem
+ * which uses the twophase box model.
+ */
+#ifndef DUMUX_EXERCISE_THREE_SPATIAL_PARAMS_HH
+#define DUMUX_EXERCISE_THREE_SPATIAL_PARAMS_HH
+
+// include parent spatialparameters
+#include <dumux/material/spatialparams/fv.hh>
+
+// include material laws
+#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup TwoPBoxModel
+ *
+ * \brief The spatial parameters for the fully coupled tutorial problem
+ * which uses the twophase box model.
+ */
+template<class FVGridGeometry, class Scalar>
+class ExerciseThreeSpatialParams
+: public FVSpatialParams<FVGridGeometry, Scalar, ExerciseThreeSpatialParams<FVGridGeometry, Scalar>>
+{
+ using ThisType = ExerciseThreeSpatialParams<FVGridGeometry, Scalar>;
+ using ParentType = FVSpatialParams<FVGridGeometry, Scalar, ThisType>;
+ using GridView = typename FVGridGeometry::GridView;
+
+ static constexpr int dim = GridView::dimension;
+ static constexpr int dimWorld = GridView::dimensionworld;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+public:
+ // export permeability type
+ using PermeabilityType = Dune::FieldMatrix<Scalar, dim, dim>;
+
+ using MaterialLaw = EffToAbsLaw<RegularizedBrooksCorey<Scalar>>;
+ using MaterialLawParams = typename MaterialLaw::Params;
+
+ /*!
+ * \brief The constructor
+ *
+ * \param fvGridGeometry The finite volume grid geometry
+ */
+ ExerciseThreeSpatialParams(std::shared_ptr<const FVGridGeometry>& fvGridGeometry)
+ : ParentType(fvGridGeometry)
+ , K_(0)
+ , KLens_(0)
+ {
+ //set main diagonal entries of the permeability tensor to a value
+ //setting to one value means: isotropic, homogeneous
+ for (int i = 0; i < dim; i++)
+ {
+ K_[i][i] = 1e-7;
+ KLens_[i][i] = 1e-10;
+ }
+
+ //set residual saturations
+ materialParams_.setSwr(0.0);
+ materialParamsLens_.setSwr(0.1);
+ materialParams_.setSnr(0.0);
+ materialParamsLens_.setSnr(0.1);
+
+ //parameters of Brooks & Corey Law
+ materialParams_.setPe(500.0);
+ materialParamsLens_.setPe(1000.0);
+ materialParams_.setLambda(2);
+ materialParamsLens_.setLambda(2);
+ }
+
+
+ /*!
+ * \brief Define the intrinsic permeability \f$\mathrm{[m^2]}\f$.
+ *
+ * \param globalPos The global position
+ */
+ PermeabilityType permeabilityAtPos(const GlobalPosition& globalPos) const
+
+ {
+ if (isInLens(globalPos))
+ return KLens_;
+ return K_;
+ }
+
+ /*!
+ * \brief Define the porosity \f$\mathrm{[-]}\f$.
+ *
+ * \param globalPos The global position
+ */
+ Scalar porosityAtPos(const GlobalPosition& globalPos) const
+ {
+ if (isInLens(globalPos))
+ return 0.1;
+ return 0.2;
+ }
+
+ /*!
+ * \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 (isInLens(globalPos))
+ return materialParamsLens_;
+ return materialParams_;
+ }
+
+ /*!
+ * \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::phase0Idx; }
+
+ //! if we are in the lens
+ bool isInLens(const GlobalPosition& globalPos) const
+ {
+ const auto x = globalPos[0];
+ const auto y = globalPos[1];
+ return (x < 40 && x > 20 && y > 35 && y < 45) ||
+ (x < 50 && x > 30 && y < 30 && y > 15);
+ }
+
+private:
+
+ Dune::FieldMatrix<Scalar, dim, dim> K_;
+ Dune::FieldMatrix<Scalar, dim, dim> KLens_;
+ // Object that holds the values/parameters of the selected material law.
+ MaterialLawParams materialParams_;
+ MaterialLawParams materialParamsLens_;
+};
+} // end namespace Dumux
+#endif
diff --git a/exercises/exercise-properties/CMakeLists.txt b/exercises/exercise-properties/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..d700163f5401ed4ba4adf7d9d434e7a969a50b5f
--- /dev/null
+++ b/exercises/exercise-properties/CMakeLists.txt
@@ -0,0 +1,11 @@
+# the compositional two-phase simulation program
+dune_add_test(NAME exercise2
+ SOURCES exercise2.cc
+ COMPILE_DEFINITIONS TYPETAG=Injection2p2pcCCProblem
+ CMD_ARGS exercise2.input)
+
+# add tutorial to the common target
+add_dependencies(test_exercises exercise2)
+
+# add a symlink for the input file
+dune_symlink_to_source_files(FILES "exercise2.input")
diff --git a/exercises/exercise-properties/README.md b/exercises/exercise-properties/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..aff0fbefb85e0f014122b3899e4b0ae82feb7569
--- /dev/null
+++ b/exercises/exercise-properties/README.md
@@ -0,0 +1,151 @@
+# Exercise #2 (DuMuX course)
+<br>
+## Problem set-up
+
+The problem setup is identical to the previous [_exercise 1_](../ex1/README.md).
+
+## Preparing the exercise
+
+* Navigate to the directory `dumux/tutorial/ex2`
+
+_Exercise 2_ deals with a two-phase compositional problem (__2p2c__). Goal is to learn how to use compile and runtime parameters and the _DuMuX property system_.
+
+<br><br>
+### Task 1: Getting familiar with the code
+<hr>
+
+Locate all the files you will need for this exercise
+* The __main file__: `exercise2.cc`
+* The __problem file__: `injection2p2cproblem.hh`
+* The __spatial parameters file__: `injection2p2cspatialparams.hh`
+* The __input file__: `exercise2.input`
+* Two header files containing:
+ * a custom __local residual__ in: `mylocalresidual.hh`
+ * a custom __material law__ in: `mymateriallaw.hh`
+
+
+<hr><br><br>
+### Task 2: Compiling and running the program
+<hr>
+
+* Change to the build-directory
+
+```bash
+cd ../../build-cmake/tutorial/ex2
+```
+
+* Compile the executable `exercise2`
+
+```bash
+make exercise2
+```
+
+* Run the problem and inspect the result
+
+```bash
+./exercise2
+```
+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
+```
+
+<hr><br><br>
+### Task 4: Enable/Disable Gravity -> DuMuX parameters
+<hr>
+
+DuMuX has many parameters that have default values. For example, all simulations consider gravity effects by default.
+You can disable gravity for a study, simply by setting the parameter in the input file
+
+```ini
+[Problem]
+EnableGravity = false
+```
+
+Run the simulation with and without gravity. Change the `Problem.Name` parameter to create output files with different
+names. Compare the results using paraview. You should immediately see the difference.
+
+A list of parameters that can be set through the input file is given [here](http://www.dumux.org/doxygen-stable/html-2.11/a06387.php).
+
+
+<hr><br><br>
+### Task 5: Implement your own 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`.
+
+* 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`
+
+```c++
+// note that every property struct knows about TypeTag
+SET_PROP(Injection2p2cTypeTag, LocalResidual)
+{
+ using type = MyCompositionalLocalResidual<TypeTag>;
+};
+
+// or using the convenience macro
+SET_TYPE_PROP(Injection2p2cTypeTag, LocalResidual,
+ MyCompositionalLocalResidual<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.
+
+* 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
+```
+
+* Verify the difference in the parameter $`x_w^{N2}`$, i.e. the mole fraction of nitrogen in the
+water phase, with and without diffusion.
+
+Note that due to diffusion being a slow process you
+can only see the difference in later times.
diff --git a/exercises/exercise-properties/exercise2.cc b/exercises/exercise-properties/exercise2.cc
new file mode 100644
index 0000000000000000000000000000000000000000..2a7493e6771bb5620efc20f21eb95683485ab42e
--- /dev/null
+++ b/exercises/exercise-properties/exercise2.cc
@@ -0,0 +1,200 @@
+// -*- 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 Test for the two-phase two-component CC model.
+ */
+#include <config.h>
+
+#include <ctime>
+#include <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/timer.hh>
+#include <dune/grid/io/file/dgfparser/dgfexception.hh>
+#include <dune/grid/io/file/vtk.hh>
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/dumuxmessage.hh>
+#include <dumux/common/defaultusagemessage.hh>
+
+#include <dumux/linear/amgbackend.hh>
+#include <dumux/nonlinear/privarswitchnewtonsolver.hh>
+
+#include <dumux/assembly/fvassembler.hh>
+#include <dumux/assembly/diffmethod.hh>
+
+#include <dumux/discretization/methods.hh>
+
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/grid/gridmanager.hh>
+
+#include "injection2p2cproblem.hh"
+
+int main(int argc, char** argv)try
+{
+ using namespace Dumux;
+
+ // define the type tag for this problem
+ using TypeTag = TTAG(Injection2p2pcCCTypeTag);
+
+ // initialize MPI, finalize is done automatically on exit
+ const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+ // print dumux start message
+ if (mpiHelper.rank() == 0)
+ DumuxMessage::print(/*firstCall=*/true);
+
+ // parse command line arguments and input file
+ Parameters::init(argc, argv);
+
+ // try to create a grid (from the given grid file or the input file)
+ GridManager<typename GET_PROP_TYPE(TypeTag, Grid)> gridManager;
+ gridManager.init();
+
+ ////////////////////////////////////////////////////////////
+ // run instationary non-linear problem on this grid
+ ////////////////////////////////////////////////////////////
+
+ // we compute on the leaf grid view
+ const auto& leafGridView = gridManager.grid().leafGridView();
+
+ // create the finite volume grid geometry
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
+ fvGridGeometry->update();
+
+ // the problem (initial and boundary conditions)
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ auto problem = std::make_shared<Problem>(fvGridGeometry);
+
+ // the solution vector
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ SolutionVector x(fvGridGeometry->numDofs());
+ problem->applyInitialSolution(x);
+ auto xOld = x;
+
+ // the grid variables
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+ auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
+ gridVariables->init(x, xOld);
+
+ // get some time loop parameters
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+ const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+ auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+
+ // check if we are about to restart a previously interrupted simulation
+ Scalar restartTime = 0;
+ if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart"))
+ restartTime = getParam<Scalar>("TimeLoop.Restart");
+
+ // 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
+
+ // instantiate time loop
+ auto timeLoop = std::make_shared<TimeLoop<Scalar>>(restartTime, dt, tEnd);
+ timeLoop->setMaxTimeStepSize(maxDt);
+
+ // the assembler with time loop for instationary problem
+ using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
+ auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop);
+
+ // the linear solver
+ using LinearSolver = AMGBackend<TypeTag>;
+ 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>;
+ NewtonSolver nonLinearSolver(assembler, linearSolver);
+
+ // time loop
+ timeLoop->start(); do
+ {
+ // 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);
+
+ // make the new solution the old solution
+ xOld = x;
+ gridVariables->advanceTimeStep();
+
+ // advance to the time loop to the next step
+ timeLoop->advanceTimeStep();
+
+ // report statistics of this time step
+ timeLoop->reportTimeStep();
+
+ // set new dt as suggested by the newton solver
+ timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
+
+ vtkWriter.write(timeLoop->time());
+
+ } while (!timeLoop->finished());
+
+ timeLoop->finalize(leafGridView.comm());
+
+ ////////////////////////////////////////////////////////////
+ // finalize, print dumux message to say goodbye
+ ////////////////////////////////////////////////////////////
+
+ // print dumux end message
+ if (mpiHelper.rank() == 0)
+ {
+ Parameters::print();
+ DumuxMessage::print(/*firstCall=*/false);
+ }
+
+ return 0;
+} // end main
+catch (Dumux::ParameterException &e)
+{
+ std::cerr << std::endl << e << " ---> Abort!" << std::endl;
+ return 1;
+}
+catch (Dune::DGFException & e)
+{
+ std::cerr << "DGF exception thrown (" << e <<
+ "). Most likely, the DGF file name is wrong "
+ "or the DGF file is corrupted, "
+ "e.g. missing hash at end of file or wrong number (dimensions) of entries."
+ << " ---> Abort!" << std::endl;
+ return 2;
+}
+catch (Dune::Exception &e)
+{
+ std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
+ return 3;
+}
+catch (...)
+{
+ std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
+ return 4;
+}
diff --git a/exercises/exercise-properties/exercise2.input b/exercises/exercise-properties/exercise2.input
new file mode 100644
index 0000000000000000000000000000000000000000..85ea966262be41ab57d58bc3e3376d1eb5334426
--- /dev/null
+++ b/exercises/exercise-properties/exercise2.input
@@ -0,0 +1,25 @@
+[TimeLoop]
+DtInitial = 3600 # in seconds
+TEnd = 3.154e9 # in seconds, i.e ten years
+
+[Grid]
+LowerLeft = 0 0
+UpperRight = 60 40
+Cells = 24 16
+
+[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
+
+# TODO: dumux-course-task
+# Set Problem.EnableGravity
+# Set Problem.EnableDiffusion
+
+[SpatialParams]
+PermeabilityAquitard = 1e-15 # m^2
+EntryPressureAquitard = 4.5e4 # Pa
+PermeabilityAquifer = 1e-12 # m^2
+EntryPressureAquifer = 1e4 # Pa
diff --git a/exercises/exercise-properties/injection2p2cproblem.hh b/exercises/exercise-properties/injection2p2cproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..a765f8bb4676f84a949bbeddbb1253994d02cbea
--- /dev/null
+++ b/exercises/exercise-properties/injection2p2cproblem.hh
@@ -0,0 +1,279 @@
+// -*- 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), false /*useComplexRelations*/>);
+
+// 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
new file mode 100644
index 0000000000000000000000000000000000000000..ab9951e28c9b296f67aa95ed0fd8843f35ab4f28
--- /dev/null
+++ b/exercises/exercise-properties/injection2p2cspatialparams.hh
@@ -0,0 +1,193 @@
+// -*- 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
new file mode 100644
index 0000000000000000000000000000000000000000..5c79e24b00230a1fb8beee824436fce0e72c39ab
--- /dev/null
+++ b/exercises/exercise-properties/mylocalresidual.hh
@@ -0,0 +1,163 @@
+// -*- 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 Element-wise calculation of the local residual for problems
+ * using compositional fully implicit model.
+ */
+#ifndef DUMUX_MY_COMPOSITIONAL_LOCAL_RESIDUAL_HH
+#define DUMUX_MY_COMPOSITIONAL_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.
+ *
+ */
+template<class TypeTag>
+class MyCompositionalLocalResidual: 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 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 };
+
+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
+ */
+ ResidualVector computeStorage(const Problem& problem,
+ const SubControlVolume& scv,
+ const VolumeVariables& volVars) const
+ {
+ ResidualVector storage(0.0);
+
+ // compute storage term of all components within all phases
+ 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);
+ }
+ }
+
+ return storage;
+ }
+
+ /*!
+ * \brief Evaluates the total flux of all conservation quantities
+ * 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
+ */
+ ResidualVector 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
+ 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];
+ }
+
+ //! Add advective phase energy fluxes. For isothermal model the contribution is zero.
+ EnergyLocalResidual::heatConvectionFlux(flux, fluxVars, phaseIdx);
+ }
+
+ //! Add diffusive energy fluxes. For isothermal model the contribution is zero.
+ EnergyLocalResidual::heatConductionFlux(flux, fluxVars);
+
+ return flux;
+ }
+
+protected:
+ Implementation *asImp_()
+ { return static_cast<Implementation *> (this); }
+
+ const Implementation *asImp_() const
+ { return static_cast<const Implementation *> (this); }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/exercises/exercise-properties/mymateriallaw.hh b/exercises/exercise-properties/mymateriallaw.hh
new file mode 100644
index 0000000000000000000000000000000000000000..e49c39a83c0afb776453782ac299e9cd1d3882c3
--- /dev/null
+++ b/exercises/exercise-properties/mymateriallaw.hh
@@ -0,0 +1,116 @@
+// -*- 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/extradoc/exercise1_nonisothermal.png b/exercises/extradoc/exercise1_nonisothermal.png
new file mode 100644
index 0000000000000000000000000000000000000000..938ba158eb97e83d1ae80145161b94d2ab6137cf
Binary files /dev/null and b/exercises/extradoc/exercise1_nonisothermal.png differ
diff --git a/exercises/extradoc/exercise1_setup.png b/exercises/extradoc/exercise1_setup.png
new file mode 100644
index 0000000000000000000000000000000000000000..db9c1b4acdfe9368781cd7593b0def3963146844
Binary files /dev/null and b/exercises/extradoc/exercise1_setup.png differ
diff --git a/exercises/extradoc/exercise2_properties.png b/exercises/extradoc/exercise2_properties.png
new file mode 100644
index 0000000000000000000000000000000000000000..28c700ed18443d78995d890cb7fb085a84146ff7
Binary files /dev/null and b/exercises/extradoc/exercise2_properties.png differ
diff --git a/exercises/extradoc/exercise3_a_solution.png b/exercises/extradoc/exercise3_a_solution.png
new file mode 100644
index 0000000000000000000000000000000000000000..fdf3ef23c81a0bc414554a9bc2e4f1e784e134c4
Binary files /dev/null and b/exercises/extradoc/exercise3_a_solution.png differ
diff --git a/exercises/extradoc/exercise3_a_solution2.png b/exercises/extradoc/exercise3_a_solution2.png
new file mode 100644
index 0000000000000000000000000000000000000000..ba90e6fcfa783fcc5a4eaf1c5fbe252dc46ea8d9
Binary files /dev/null and b/exercises/extradoc/exercise3_a_solution2.png differ
diff --git a/exercises/extradoc/exercise3_setup.png b/exercises/extradoc/exercise3_setup.png
new file mode 100644
index 0000000000000000000000000000000000000000..eed8371e6a96f9b4d4eee22849f4b3bbbcd83fa1
Binary files /dev/null and b/exercises/extradoc/exercise3_setup.png differ
diff --git a/exercises/solution/CMakeLists.txt b/exercises/solution/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..ac5b07f6334f00562b2464e9e07226335d1cdf86
--- /dev/null
+++ b/exercises/solution/CMakeLists.txt
@@ -0,0 +1,2 @@
+add_subdirectory(ex1)
+add_subdirectory(ex2)
diff --git a/exercises/solution/ex1/CMakeLists.txt b/exercises/solution/ex1/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..bf559fc2ddcf5f555b2972adb0e884584920ecff
--- /dev/null
+++ b/exercises/solution/ex1/CMakeLists.txt
@@ -0,0 +1,21 @@
+# the immiscible two-phase simulation program
+#dune_add_test(NAME exercise1_2p_solution
+# SOURCES exercise1_2p_solution.cc
+# CMD_ARGS exercise1.input)
+
+# the compositional two-phase two-component simulation program
+dune_add_test(NAME exercise1_2p2c_solution
+ SOURCES exercise1_2p2c.cc
+ CMD_ARGS exercise1.input)
+
+# the two-phase non-isothermal simulation program
+dune_add_test(NAME exercise1_2pni_solution
+ SOURCES exercise1_2pni_solution.cc
+ CMD_ARGS exercise1.input)
+
+
+# add tutorial to the common target
+add_dependencies(test_exercises exercise1_2p2c_solution exercise1_2pni_solution)
+
+# add a symlink for the input file
+dune_symlink_to_source_files(FILES "exercise1.input")
diff --git a/exercises/solution/ex1/exercise1.input b/exercises/solution/ex1/exercise1.input
new file mode 100644
index 0000000000000000000000000000000000000000..545a170fb594b33f1a0600be18d50495c3953ca0
--- /dev/null
+++ b/exercises/solution/ex1/exercise1.input
@@ -0,0 +1,27 @@
+[TimeLoop]
+DtInitial = 3600 # in seconds
+TEnd = 3.154e9 # in seconds, i.e ten years
+
+[Grid]
+LowerLeft = 0 0
+UpperRight = 60 40
+Cells = 24 16
+
+[Problem]
+Name = injection
+OnlyPlotMaterialLaws = true
+AquiferDepth = 2700.0 # m
+InjectionDuration = 2.628e6 # in seconds, i.e. one month
+TotalAreaSpecificInflow = -1e-4 # kg/(s*m^2)
+
+[SpatialParams]
+PermeabilityAquitard = 1e-15 # m^2
+EntryPressureAquitard = 4.5e4 # Pa
+PermeabilityAquifer = 1e-12 # m^2
+EntryPressureAquifer = 1e4 # Pa
+
+# these parameters are only used in the nonisothermal model. Uncomment them for that
+[Component]
+SolidDensity = 2700 # solid density of granite
+SolidThermalConductivity = 2.8 # solid thermal conducitivity of granite
+SolidHeatCapacity = 790 # solid heat capacity of granite
diff --git a/exercises/solution/ex1/exercise1_2p2c.cc b/exercises/solution/ex1/exercise1_2p2c.cc
new file mode 100644
index 0000000000000000000000000000000000000000..55496aad06bae940e473fa8ad013a27d23be5444
--- /dev/null
+++ b/exercises/solution/ex1/exercise1_2p2c.cc
@@ -0,0 +1,198 @@
+// -*- 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 The main file for the 2p2c porousmediumflow problem in exercise1
+ */
+#include <config.h>
+
+#include <ctime>
+#include <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/timer.hh>
+#include <dune/grid/io/file/dgfparser/dgfexception.hh>
+#include <dune/grid/io/file/vtk.hh>
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/valgrind.hh>
+#include <dumux/common/dumuxmessage.hh>
+
+#include <dumux/linear/amgbackend.hh>
+#include <dumux/nonlinear/privarswitchnewtonsolver.hh>
+
+#include <dumux/assembly/fvassembler.hh>
+#include <dumux/assembly/diffmethod.hh>
+
+#include <dumux/discretization/methods.hh>
+
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/grid/gridmanager.hh>
+
+#include "injection2p2cproblem.hh"
+
+////////////////////////
+// the main function
+////////////////////////
+int main(int argc, char** argv) try
+{
+ using namespace Dumux;
+
+ // define the type tag for this problem
+ using TypeTag = TTAG(Injection2p2cCCTypeTag);
+
+ // initialize MPI, finalize is done automatically on exit
+ const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+ // print dumux start message
+ if (mpiHelper.rank() == 0)
+ DumuxMessage::print(/*firstCall=*/true);
+
+ // parse command line arguments and input file
+ Parameters::init(argc, argv);
+
+ // try to create a grid (from the given grid file or the input file)
+ GridManager<typename GET_PROP_TYPE(TypeTag, Grid)> gridManager;
+ gridManager.init();
+
+ ////////////////////////////////////////////////////////////
+ // run instationary non-linear problem on this grid
+ ////////////////////////////////////////////////////////////
+
+ // we compute on the leaf grid view
+ const auto& leafGridView = gridManager.grid().leafGridView();
+
+ // create the finite volume grid geometry
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
+ fvGridGeometry->update();
+
+ // the problem (initial and boundary conditions)
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ auto problem = std::make_shared<Problem>(fvGridGeometry);
+
+ // the solution vector
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ SolutionVector x(fvGridGeometry->numDofs());
+ problem->applyInitialSolution(x);
+ auto xOld = x;
+
+ // the grid variables
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+ auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
+ gridVariables->init(x, xOld);
+
+ // get some time loop parameters
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+ const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+ auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+
+ // 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
+
+ // instantiate time loop
+ auto timeLoop = std::make_shared<TimeLoop<Scalar>>(0.0, dt, tEnd);
+ timeLoop->setMaxTimeStepSize(maxDt);
+
+ // the assembler with time loop for instationary problem
+ using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
+ auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop);
+
+ // the linear solver
+ using LinearSolver = AMGBackend<TypeTag>;
+ 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>;
+ NewtonSolver nonLinearSolver(assembler, linearSolver);
+
+ // time loop
+ timeLoop->start();
+ while (!timeLoop->finished())
+ {
+ // 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);
+
+ // make the new solution the old solution
+ xOld = x;
+ gridVariables->advanceTimeStep();
+
+ // advance to the time loop to the next step
+ timeLoop->advanceTimeStep();
+
+ // report statistics of this time step
+ timeLoop->reportTimeStep();
+
+ // set new dt as suggested by the newton solver
+ timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
+
+ // output to vtk
+ vtkWriter.write(timeLoop->time());
+ }
+
+ timeLoop->finalize(leafGridView.comm());
+
+ ////////////////////////////////////////////////////////////
+ // finalize, print dumux message to say goodbye
+ ////////////////////////////////////////////////////////////
+
+ // print dumux end message
+ if (mpiHelper.rank() == 0)
+ {
+ Parameters::print();
+ DumuxMessage::print(/*firstCall=*/false);
+ }
+
+ return 0;
+} // end main
+catch (Dumux::ParameterException &e)
+{
+ std::cerr << std::endl << e << " ---> Abort!" << std::endl;
+ return 1;
+}
+catch (Dune::DGFException & e)
+{
+ std::cerr << "DGF exception thrown (" << e <<
+ "). Most likely, the DGF file name is wrong "
+ "or the DGF file is corrupted, "
+ "e.g. missing hash at end of file or wrong number (dimensions) of entries."
+ << " ---> Abort!" << std::endl;
+ return 2;
+}
+catch (Dune::Exception &e)
+{
+ std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
+ return 3;
+}
+catch (...)
+{
+ std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
+ return 4;
+}
diff --git a/exercises/solution/ex1/exercise1_2pni_solution.cc b/exercises/solution/ex1/exercise1_2pni_solution.cc
new file mode 100644
index 0000000000000000000000000000000000000000..68f9d3b50ffa07c4aedf2e9843174da557e1544d
--- /dev/null
+++ b/exercises/solution/ex1/exercise1_2pni_solution.cc
@@ -0,0 +1,198 @@
+// -*- 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 The solution main file for the two-phase porousmediumflow problem of exercise 1
+ */
+#include <config.h>
+
+#include <ctime>
+#include <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/timer.hh>
+#include <dune/grid/io/file/dgfparser/dgfexception.hh>
+#include <dune/grid/io/file/vtk.hh>
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/dumuxmessage.hh>
+#include <dumux/common/defaultusagemessage.hh>
+
+#include <dumux/linear/amgbackend.hh>
+#include <dumux/nonlinear/newtonsolver.hh>
+
+#include <dumux/assembly/fvassembler.hh>
+#include <dumux/assembly/diffmethod.hh>
+
+#include <dumux/discretization/methods.hh>
+
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/grid/gridmanager.hh>
+
+// The problem file, where setup-specific boundary and initial conditions are defined.
+#include "injection2pniproblem.hh"
+
+////////////////////////
+// the main function
+////////////////////////
+int main(int argc, char** argv) try
+{
+ using namespace Dumux;
+
+ // define the type tag for this problem
+ using TypeTag = TTAG(Injection2pNICCTypeTag);
+
+ // initialize MPI, finalize is done automatically on exit
+ const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+ // print dumux start message
+ if (mpiHelper.rank() == 0)
+ DumuxMessage::print(/*firstCall=*/true);
+
+ // parse command line arguments and input file
+ Parameters::init(argc, argv);
+
+ // try to create a grid (from the given grid file or the input file)
+ GridManager<typename GET_PROP_TYPE(TypeTag, Grid)> gridManager;
+ gridManager.init();
+
+ ////////////////////////////////////////////////////////////
+ // run instationary non-linear problem on this grid
+ ////////////////////////////////////////////////////////////
+
+ // we compute on the leaf grid view
+ const auto& leafGridView = gridManager.grid().leafGridView();
+
+ // create the finite volume grid geometry
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
+ fvGridGeometry->update();
+
+ // the problem (initial and boundary conditions)
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ auto problem = std::make_shared<Problem>(fvGridGeometry);
+
+ // the solution vector
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ SolutionVector x(fvGridGeometry->numDofs());
+ problem->applyInitialSolution(x);
+ auto xOld = x;
+
+ // the grid variables
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+ auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
+ gridVariables->init(x, xOld);
+
+ // get some time loop parameters
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+ const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+ auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+
+ // 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
+
+ // instantiate time loop
+ auto timeLoop = std::make_shared<TimeLoop<Scalar>>(0.0, dt, tEnd);
+ timeLoop->setMaxTimeStepSize(maxDt);
+
+ // the assembler with time loop for instationary problem
+ using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
+ auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop);
+
+ // the linear solver
+ using LinearSolver = AMGBackend<TypeTag>;
+ auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->dofMapper());
+
+ // the non-linear solver
+ using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
+ NewtonSolver nonLinearSolver(assembler, linearSolver);
+
+ // time loop
+ timeLoop->start();
+ while (!timeLoop->finished())
+ {
+ // 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);
+
+ // make the new solution the old solution
+ xOld = x;
+ gridVariables->advanceTimeStep();
+
+ // advance to the time loop to the next step
+ timeLoop->advanceTimeStep();
+
+ // report statistics of this time step
+ timeLoop->reportTimeStep();
+
+ // set new dt as suggested by the newton solver
+ timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
+
+ // output to vtk
+ vtkWriter.write(timeLoop->time());
+ }
+
+ timeLoop->finalize(leafGridView.comm());
+
+ ////////////////////////////////////////////////////////////
+ // finalize, print dumux message to say goodbye
+ ////////////////////////////////////////////////////////////
+
+ // print dumux end message
+ if (mpiHelper.rank() == 0)
+ {
+ Parameters::print();
+ DumuxMessage::print(/*firstCall=*/false);
+ }
+
+ return 0;
+} // end main
+catch (Dumux::ParameterException &e)
+{
+ std::cerr << std::endl << e << " ---> Abort!" << std::endl;
+ return 1;
+}
+catch (Dune::DGFException & e)
+{
+ std::cerr << "DGF exception thrown (" << e <<
+ "). Most likely, the DGF file name is wrong "
+ "or the DGF file is corrupted, "
+ "e.g. missing hash at end of file or wrong number (dimensions) of entries."
+ << " ---> Abort!" << std::endl;
+ return 2;
+}
+catch (Dune::Exception &e)
+{
+ std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
+ return 3;
+}
+catch (...)
+{
+ std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
+ return 4;
+}
diff --git a/exercises/solution/ex1/injection2p2cproblem.hh b/exercises/solution/ex1/injection2p2cproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..d33e5a46bb25ea3bffec3f54a804af1947d5f7d9
--- /dev/null
+++ b/exercises/solution/ex1/injection2p2cproblem.hh
@@ -0,0 +1,271 @@
+// -*- 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 The two-phase porousmediumflow problem for exercise 1
+ */
+#ifndef DUMUX_EX1_INJECTION_2P2C_PROBLEM_HH
+#define DUMUX_EX1_INJECTION_2P2C_PROBLEM_HH
+
+#include <dumux/discretization/cellcentered/tpfa/properties.hh>
+#include <dumux/porousmediumflow/2p2c/model.hh>
+#include <dumux/porousmediumflow/problem.hh>
+#include <dumux/material/fluidsystems/h2on2.hh>
+
+#include "injection2pspatialparams.hh"
+
+namespace Dumux {
+
+// forward declare problem
+template <class TypeTag>
+class Injection2p2cProblem;
+
+namespace Properties {
+NEW_TYPE_TAG(Injection2p2cTypeTag, INHERITS_FROM(TwoPTwoC));
+NEW_TYPE_TAG(Injection2p2cCCTypeTag, 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 fluid configuration
+SET_TYPE_PROP(Injection2p2cTypeTag, FluidSystem, FluidSystems::H2ON2<typename GET_PROP_TYPE(TypeTag, Scalar), /*useComplexRelations=*/ false>);
+
+// Define whether mole(true) or mass (false) fractions are used
+SET_BOOL_PROP(Injection2p2cTypeTag, UseMoles, true);
+} // end namespace Properties
+
+/*!
+ * \ingroup TwoPTwoCModel
+ * \ingroup ImplicitTestProblems
+ * \brief Gas injection problem where a gas (here nitrogen) is injected into a fully
+ * water saturated medium. During buoyancy driven upward migration the gas
+ * passes a high temperature area.
+ *
+ * The domain is sized 60 m times 40 m.
+ *
+ * For the mass conservation equation neumann boundary conditions are used on
+ * the top, on the bottom and on the right of the domain, while dirichlet conditions
+ * apply on the left boundary.
+ *
+ * Gas is injected at the right boundary from 7 m to 15 m at a rate of
+ * 0.001 kg/(s m), the remaining neumann boundaries are no-flow
+ * boundaries.
+ *
+ * At the dirichlet boundaries a hydrostatic pressure and a gas saturation of zero a
+ *
+ * This problem uses the \ref TwoPModel model.
+ */
+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 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;
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+
+ 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
+ // getParam<TYPE>("GROUPNAME.PARAMNAME") reads and sets parameter PARAMNAME
+ // of type TYPE given in the group GROUPNAME from the input file
+ name_ = getParam<std::string>("Problem.Name");
+ // depth of the aquifer, units: m
+ aquiferDepth_ = getParam<Scalar>("Problem.AquiferDepth");
+ // the duration of the injection, units: second
+ injectionDuration_ = getParam<Scalar>("Problem.InjectionDuration");
+ totalAreaSpecificInflow_ = getParam<Scalar>("Problem.TotalAreaSpecificInflow");
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief Returns the problem name
+ *
+ * This is used as a prefix for files generated by the simulation.
+ */
+ std::string name() const
+ { return name_+"-2p2c"; }
+
+ /*!
+ * \brief Returns the temperature \f$ K \f$
+ */
+ Scalar temperature() const
+ {
+ return 273.15 + 30; // [K]
+ }
+
+ // \}
+
+ /*!
+ * \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();
+ 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])
+ {
+ // TODO: dumux-course-task
+ //instead of setting -1e-4 here directly use totalAreaSpecificInflow_ in the computation
+
+ // inject nitrogen. negative values mean injection
+ // 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::pressureIdx] = pw;
+ values[Indices::switchIdx] = moleFracLiquidN2;
+
+ return values;
+ }
+
+ // \}
+
+ //! set the time for the time dependent boundary conditions (called from main)
+ 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 injectionDuration_; //! Duration of the injection in seconds
+ Scalar time_;
+ Scalar totalAreaSpecificInflow_;
+
+};
+
+} //end namespace Dumux
+
+#endif
diff --git a/exercises/solution/ex1/injection2pniproblem.hh b/exercises/solution/ex1/injection2pniproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..70135414a985fcfa0bbf4efdd75ccf0e9306763e
--- /dev/null
+++ b/exercises/solution/ex1/injection2pniproblem.hh
@@ -0,0 +1,252 @@
+// -*- 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 The two-phase nonisothermal porousmediumflow problem for exercise 1
+ */
+
+#ifndef DUMUX_EX1_INJECTION_PROBLEM_2PNI_HH
+#define DUMUX_EX1_INJECTION_PROBLEM_2PNI_HH
+
+#include <dumux/discretization/cellcentered/tpfa/properties.hh>
+#include <dumux/porousmediumflow/2p/model.hh>
+#include <dumux/porousmediumflow/problem.hh>
+#include <dumux/material/fluidsystems/h2on2.hh>
+
+#include "injection2pspatialparams.hh"
+
+namespace Dumux {
+
+// forward declare problem
+template <class TypeTag>
+class InjectionProblem2PNI;
+
+namespace Properties {
+NEW_TYPE_TAG(Injection2pNITypeTag, INHERITS_FROM(TwoPNI));
+NEW_TYPE_TAG(Injection2pNICCTypeTag, INHERITS_FROM(CCTpfaModel, Injection2pNITypeTag));
+
+// Set the grid type
+SET_TYPE_PROP(Injection2pNITypeTag, Grid, Dune::YaspGrid<2>);
+
+// Set the problem property
+SET_TYPE_PROP(Injection2pNITypeTag, Problem, InjectionProblem2PNI<TypeTag>);
+
+// Set the spatial parameters
+SET_TYPE_PROP(Injection2pNITypeTag, SpatialParams,
+ InjectionSpatialParams<typename GET_PROP_TYPE(TypeTag, FVGridGeometry),
+ typename GET_PROP_TYPE(TypeTag, Scalar)>);
+
+// Set fluid configuration
+SET_TYPE_PROP(Injection2pNITypeTag, FluidSystem, FluidSystems::H2ON2<typename GET_PROP_TYPE(TypeTag, Scalar), false>);
+} // end namespace Properties
+
+/*!
+ * \ingroup TwoPModel
+ * \ingroup ImplicitTestProblems
+ * \brief Gas injection problem where a gas (here nitrogen) is injected into a fully
+ * water saturated medium. During buoyancy driven upward migration the gas
+ * passes a high temperature area.
+ *
+ * The domain is sized 60 m times 40 m.
+ *
+ * For the mass conservation equation neumann boundary conditions are used on
+ * the top, on the bottom and on the right of the domain, while dirichlet conditions
+ * apply on the left boundary.
+ *
+ * Gas is injected at the right boundary from 7 m to 15 m at a rate of
+ * 0.001 kg/(s m), the remaining neumann boundaries are no-flow
+ * boundaries.
+ *
+ * At the dirichlet boundaries a hydrostatic pressure and a gas saturation of zero a
+ *
+ * This problem uses the \ref TwoPModel model.
+ */
+template<class TypeTag>
+class InjectionProblem2PNI : public PorousMediumFlowProblem<TypeTag>
+{
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ 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;
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+
+ enum { dimWorld = GridView::dimensionworld };
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+public:
+ InjectionProblem2PNI(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
+ // getParam<TYPE>("GROUPNAME.PARAMNAME") reads and sets parameter PARAMNAME
+ // of type TYPE given in the group GROUPNAME from the input file
+ name_ = getParam<std::string>("Problem.Name");
+ // depth of the aquifer, units: m
+ aquiferDepth_ = getParam<Scalar>("Problem.AquiferDepth");
+ // 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.
+ */
+ std::string name() const
+ { return "injection-2pni"; }
+
+ // \}
+
+ /*!
+ * \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();
+ 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])
+ {
+ // inject nitrogen. negative values mean injection
+ // units kg/(s*m^2)
+ values[Indices::conti0EqIdx + FluidSystem::N2Idx] = -1e-4;
+ values[Indices::conti0EqIdx + FluidSystem::H2OIdx] = 0.0;
+ values[Indices::energyEqIdx] = 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);
+
+ // get the water density at atmospheric conditions
+ const Scalar densityW = FluidSystem::H2O::liquidDensity(283.15, 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]);
+
+ values[Indices::pressureIdx] = pw;
+ values[Indices::saturationIdx] = 0.0;
+
+ values[Indices::temperatureIdx] = 283.0 + (aquiferDepth_ - globalPos[1])*0.03;
+ if (globalPos[0] > 20 - eps_ && globalPos[0] < 30 + eps_ && globalPos[1] > 5 - eps_ && globalPos[1] < 35 + eps_)
+ values[Indices::temperatureIdx] = 380;
+
+ return values;
+ }
+
+ // \}
+
+ //! set the time for the time dependent boundary conditions (called from main)
+ 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 injectionDuration_; //! Duration of the injection in seconds
+ Scalar time_;
+};
+
+} //end namespace Dumux
+
+#endif
diff --git a/exercises/solution/ex1/injection2pspatialparams.hh b/exercises/solution/ex1/injection2pspatialparams.hh
new file mode 100644
index 0000000000000000000000000000000000000000..979f14ce9782f08b1a08010b0f8053506265ebf3
--- /dev/null
+++ b/exercises/solution/ex1/injection2pspatialparams.hh
@@ -0,0 +1,164 @@
+// -*- 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_EX1_INJECTION_SPATIAL_PARAMS_HH
+#define DUMUX_EX1_INJECTION_SPATIAL_PARAMS_HH
+
+#include <dumux/material/spatialparams/fv.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.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;
+
+ using MaterialLaw = EffToAbsLaw<RegularizedBrooksCorey<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);
+ }
+
+ /*!
+ * \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; }
+
+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/ex2/CMakeLists.txt b/exercises/solution/ex2/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..fe1e87464e3d126605f090e6b72ad21ae636f920
--- /dev/null
+++ b/exercises/solution/ex2/CMakeLists.txt
@@ -0,0 +1,10 @@
+# the compositional two-phase simulation program
+dune_add_test(NAME exercise2_solution
+ SOURCES exercise2_solution.cc
+ CMD_ARGS exercise2.input -Problem.OnlyPlotMaterialLaws false)
+
+# add tutorial to the common target
+add_dependencies(test_exercises exercise2_solution)
+
+# add a symlink for the input file
+dune_symlink_to_source_files(FILES "exercise2.input")
diff --git a/exercises/solution/ex2/exercise2.input b/exercises/solution/ex2/exercise2.input
new file mode 100644
index 0000000000000000000000000000000000000000..df7737b869ce9359d9f4cbfdfcdf5edbb8c4e5a3
--- /dev/null
+++ b/exercises/solution/ex2/exercise2.input
@@ -0,0 +1,24 @@
+[TimeLoop]
+DtInitial = 3600 # in seconds
+TEnd = 3.154e9 # in seconds, i.e ten years
+
+[Grid]
+LowerLeft = 0 0
+UpperRight = 60 40
+Cells = 24 16
+
+[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
+
+EnableGravity = true
+EnableDiffusion = true
+
+[SpatialParams]
+PermeabilityAquitard = 1e-15 # m^2
+EntryPressureAquitard = 4.5e4 # Pa
+PermeabilityAquifer = 1e-12 # m^2
+EntryPressureAquifer = 1e4 # Pa
diff --git a/exercises/solution/ex2/exercise2_solution.cc b/exercises/solution/ex2/exercise2_solution.cc
new file mode 100644
index 0000000000000000000000000000000000000000..2a7493e6771bb5620efc20f21eb95683485ab42e
--- /dev/null
+++ b/exercises/solution/ex2/exercise2_solution.cc
@@ -0,0 +1,200 @@
+// -*- 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 Test for the two-phase two-component CC model.
+ */
+#include <config.h>
+
+#include <ctime>
+#include <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/timer.hh>
+#include <dune/grid/io/file/dgfparser/dgfexception.hh>
+#include <dune/grid/io/file/vtk.hh>
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/dumuxmessage.hh>
+#include <dumux/common/defaultusagemessage.hh>
+
+#include <dumux/linear/amgbackend.hh>
+#include <dumux/nonlinear/privarswitchnewtonsolver.hh>
+
+#include <dumux/assembly/fvassembler.hh>
+#include <dumux/assembly/diffmethod.hh>
+
+#include <dumux/discretization/methods.hh>
+
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/grid/gridmanager.hh>
+
+#include "injection2p2cproblem.hh"
+
+int main(int argc, char** argv)try
+{
+ using namespace Dumux;
+
+ // define the type tag for this problem
+ using TypeTag = TTAG(Injection2p2pcCCTypeTag);
+
+ // initialize MPI, finalize is done automatically on exit
+ const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+ // print dumux start message
+ if (mpiHelper.rank() == 0)
+ DumuxMessage::print(/*firstCall=*/true);
+
+ // parse command line arguments and input file
+ Parameters::init(argc, argv);
+
+ // try to create a grid (from the given grid file or the input file)
+ GridManager<typename GET_PROP_TYPE(TypeTag, Grid)> gridManager;
+ gridManager.init();
+
+ ////////////////////////////////////////////////////////////
+ // run instationary non-linear problem on this grid
+ ////////////////////////////////////////////////////////////
+
+ // we compute on the leaf grid view
+ const auto& leafGridView = gridManager.grid().leafGridView();
+
+ // create the finite volume grid geometry
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
+ fvGridGeometry->update();
+
+ // the problem (initial and boundary conditions)
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ auto problem = std::make_shared<Problem>(fvGridGeometry);
+
+ // the solution vector
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ SolutionVector x(fvGridGeometry->numDofs());
+ problem->applyInitialSolution(x);
+ auto xOld = x;
+
+ // the grid variables
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+ auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
+ gridVariables->init(x, xOld);
+
+ // get some time loop parameters
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+ const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+ auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+
+ // check if we are about to restart a previously interrupted simulation
+ Scalar restartTime = 0;
+ if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart"))
+ restartTime = getParam<Scalar>("TimeLoop.Restart");
+
+ // 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
+
+ // instantiate time loop
+ auto timeLoop = std::make_shared<TimeLoop<Scalar>>(restartTime, dt, tEnd);
+ timeLoop->setMaxTimeStepSize(maxDt);
+
+ // the assembler with time loop for instationary problem
+ using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
+ auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop);
+
+ // the linear solver
+ using LinearSolver = AMGBackend<TypeTag>;
+ 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>;
+ NewtonSolver nonLinearSolver(assembler, linearSolver);
+
+ // time loop
+ timeLoop->start(); do
+ {
+ // 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);
+
+ // make the new solution the old solution
+ xOld = x;
+ gridVariables->advanceTimeStep();
+
+ // advance to the time loop to the next step
+ timeLoop->advanceTimeStep();
+
+ // report statistics of this time step
+ timeLoop->reportTimeStep();
+
+ // set new dt as suggested by the newton solver
+ timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
+
+ vtkWriter.write(timeLoop->time());
+
+ } while (!timeLoop->finished());
+
+ timeLoop->finalize(leafGridView.comm());
+
+ ////////////////////////////////////////////////////////////
+ // finalize, print dumux message to say goodbye
+ ////////////////////////////////////////////////////////////
+
+ // print dumux end message
+ if (mpiHelper.rank() == 0)
+ {
+ Parameters::print();
+ DumuxMessage::print(/*firstCall=*/false);
+ }
+
+ return 0;
+} // end main
+catch (Dumux::ParameterException &e)
+{
+ std::cerr << std::endl << e << " ---> Abort!" << std::endl;
+ return 1;
+}
+catch (Dune::DGFException & e)
+{
+ std::cerr << "DGF exception thrown (" << e <<
+ "). Most likely, the DGF file name is wrong "
+ "or the DGF file is corrupted, "
+ "e.g. missing hash at end of file or wrong number (dimensions) of entries."
+ << " ---> Abort!" << std::endl;
+ return 2;
+}
+catch (Dune::Exception &e)
+{
+ std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
+ return 3;
+}
+catch (...)
+{
+ std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
+ return 4;
+}
diff --git a/exercises/solution/ex2/injection2p2cproblem.hh b/exercises/solution/ex2/injection2p2cproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..8c62fdf59d81d85b551f7a3e1d38846657a70106
--- /dev/null
+++ b/exercises/solution/ex2/injection2p2cproblem.hh
@@ -0,0 +1,278 @@
+// -*- 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), false /*useComplexRelations*/>);
+
+// 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/ex2/injection2p2cspatialparams.hh b/exercises/solution/ex2/injection2p2cspatialparams.hh
new file mode 100644
index 0000000000000000000000000000000000000000..4c89071a25ab8dd1964a44b4a1093dbd86e91c97
--- /dev/null
+++ b/exercises/solution/ex2/injection2p2cspatialparams.hh
@@ -0,0 +1,190 @@
+// -*- 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/ex2/mylocalresidual.hh b/exercises/solution/ex2/mylocalresidual.hh
new file mode 100644
index 0000000000000000000000000000000000000000..ac6685555de9f509779b0e972700d5c802f82409
--- /dev/null
+++ b/exercises/solution/ex2/mylocalresidual.hh
@@ -0,0 +1,163 @@
+// -*- 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 Element-wise calculation of the local residual for problems
+ * using compositional fully implicit model.
+ */
+#ifndef DUMUX_MY_COMPOSITIONAL_LOCAL_RESIDUAL_HH
+#define DUMUX_MY_COMPOSITIONAL_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.
+ *
+ */
+template<class TypeTag>
+class MyCompositionalLocalResidual: 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 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 };
+
+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
+ */
+ ResidualVector computeStorage(const Problem& problem,
+ const SubControlVolume& scv,
+ const VolumeVariables& volVars) const
+ {
+ ResidualVector storage(0.0);
+
+ // compute storage term of all components within all phases
+ 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);
+ }
+ }
+
+ return storage;
+ }
+
+ /*!
+ * \brief Evaluates the total flux of all conservation quantities
+ * 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
+ */
+ ResidualVector 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
+ 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];
+ }
+
+ //! Add advective phase energy fluxes. For isothermal model the contribution is zero.
+ EnergyLocalResidual::heatConvectionFlux(flux, fluxVars, phaseIdx);
+ }
+
+ //! Add diffusive energy fluxes. For isothermal model the contribution is zero.
+ EnergyLocalResidual::heatConductionFlux(flux, fluxVars);
+
+ return flux;
+ }
+
+protected:
+ Implementation *asImp_()
+ { return static_cast<Implementation *> (this); }
+
+ const Implementation *asImp_() const
+ { return static_cast<const Implementation *> (this); }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/exercises/solution/ex2/mymateriallaw.hh b/exercises/solution/ex2/mymateriallaw.hh
new file mode 100644
index 0000000000000000000000000000000000000000..fe36bec694fe2d89363ffa9693274ac2efe4b5a7
--- /dev/null
+++ b/exercises/solution/ex2/mymateriallaw.hh
@@ -0,0 +1,114 @@
+// -*- 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/ex3/h2omycompressiblecomponent.hh b/exercises/solution/ex3/h2omycompressiblecomponent.hh
new file mode 100644
index 0000000000000000000000000000000000000000..7ffb47b3aa7b7c0b63333a3007b35b165adfa09e
--- /dev/null
+++ b/exercises/solution/ex3/h2omycompressiblecomponent.hh
@@ -0,0 +1,489 @@
+// -*- 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 A fluid system with water and a fictitious component, which is to be
+ * implemented in tutorial exercise 3a, as phases and components. This
+ * fluid system is to be implemented in exercise 3b.
+ */
+#ifndef DUMUX_H2O_MYCOMPRESSIBLECOMPONENT_FLUID_SYSTEM_HH
+#define DUMUX_H2O_MYCOMPRESSIBLECOMPONENT_FLUID_SYSTEM_HH
+
+#include <dumux/material/components/h2o.hh>
+#include <dumux/material/components/tabulatedcomponent.hh>
+
+#include <dumux/material/fluidsystems/base.hh>
+
+// the fictitious component that was created in exercise 3a
+#include <tutorial/ex3/components/mycompressiblecomponent.hh>
+
+// the binary coefficients corresponding to this fluid system
+#include <tutorial/ex3/binarycoefficients/h2omycompressiblecomponent.hh>
+
+namespace Dumux
+{
+namespace FluidSystems
+{
+
+/*!
+ * \brief A compositional fluid consisting of two liquid phases,
+ * which are water and a fictitious component from tutorial exercise 3a.
+ */
+template <class TypeTag, class Scalar,
+ class H2OType = Dumux::Components::TabulatedComponent<Dumux::Components::H2O<Scalar> > >
+class H2OMyCompressibleComponent
+ : public BaseFluidSystem< Scalar, H2OMyCompressibleComponent<TypeTag, Scalar, H2OType> >
+{
+ typedef H2OMyCompressibleComponent<TypeTag, Scalar, H2OType> ThisType;
+ typedef BaseFluidSystem<Scalar, ThisType> Base;
+
+public:
+ typedef Dumux::MyCompressibleComponent<Scalar> MyCompressibleComponent;
+ typedef H2OType H2O;
+
+ static constexpr int numPhases = 2;
+ static constexpr int numComponents = 2;
+
+ static constexpr int phase0Idx = 0; // index of the first phase
+ static constexpr int phase1Idx = 1; // index of the second phase
+
+ static constexpr int H2OIdx = 0;
+ static constexpr int NAPLIdx = 1;
+ // export component indices to indicate the main component
+ // of the corresponding phase at atmospheric pressure 1 bar
+ // and room temperature 20°C:
+ static constexpr int comp0Idx = H2OIdx;
+ static constexpr int comp1Idx = NAPLIdx;
+
+ /*!
+ * \brief Initialize the fluid system's static parameters generically
+ *
+ * If a tabulated H2O component is used, we do our best to create
+ * tables that always work.
+ */
+ static void init()
+ {
+ init(/*tempMin=*/273.15,
+ /*tempMax=*/623.15,
+ /*numTemp=*/100,
+ /*pMin=*/0.0,
+ /*pMax=*/20e6,
+ /*numP=*/200);
+ }
+
+ /*!
+ * \brief Initialize the fluid system's static parameters using
+ * problem specific temperature and pressure ranges
+ *
+ * \param tempMin The minimum temperature used for tabulation of water [K]
+ * \param tempMax The maximum temperature used for tabulation of water [K]
+ * \param nTemp The number of ticks on the temperature axis of the table of water
+ * \param pressMin The minimum pressure used for tabulation of water [Pa]
+ * \param pressMax The maximum pressure used for tabulation of water [Pa]
+ * \param nPress The number of ticks on the pressure axis of the table of water
+ */
+ static void init(Scalar tempMin, Scalar tempMax, unsigned nTemp,
+ Scalar pressMin, Scalar pressMax, unsigned nPress)
+ {
+ if (H2O::isTabulated) {
+ std::cout << "Initializing tables for the H2O fluid properties ("
+ << nTemp*nPress
+ << " entries).\n";
+
+ H2O::init(tempMin, tempMax, nTemp,
+ pressMin, pressMax, nPress);
+ }
+ }
+
+
+ /*!
+ * \brief Return whether a phase is liquid
+ *
+ * \param phaseIdx The index of the fluid phase to consider
+ */
+ static constexpr bool isLiquid(int phaseIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+ return true;
+ }
+
+ static constexpr bool isIdealGas(int phaseIdx)
+ { return H2O::gasIsIdeal() && MyCompressibleComponent::gasIsIdeal(); }
+
+ /*!
+ * \brief Returns true if and only if a fluid phase is assumed to
+ * be an ideal mixture.
+ *
+ * We define an ideal mixture as a fluid phase where the fugacity
+ * coefficients of all components times the pressure of the phase
+ * are indepent on the fluid composition. This assumtion is true
+ * if Henry's law and Raoult's law apply. If you are unsure what
+ * this function should return, it is safe to return false. The
+ * only damage done will be (slightly) increased computation times
+ * in some cases.
+ *
+ * \param phaseIdx The index of the fluid phase to consider
+ */
+ static constexpr bool isIdealMixture(int phaseIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+ return true;
+ }
+
+ /*!
+ * \brief Returns true if and only if a fluid phase is assumed to
+ * be compressible.
+ *
+ * Compressible means that the partial derivative of the density
+ * to the fluid pressure is always larger than zero.
+ *
+ * \param phaseIdx The index of the fluid phase to consider
+ */
+ static constexpr bool isCompressible(int phaseIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+ if (phaseIdx == phase0Idx)
+ // the water component decides for the water phase...
+ return H2O::liquidIsCompressible();
+
+ // the NAPL component decides for the napl phase...
+ return MyCompressibleComponent::liquidIsCompressible();
+ }
+
+ /*!
+ * \brief Return the human readable name of a phase (used in indices)
+ */
+ static std::string phaseName(int phaseIdx)
+ {
+ switch (phaseIdx) {
+ case phase0Idx: return "w";
+ case phase1Idx: return "n";
+ };
+ DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
+ }
+
+ /*!
+ * \brief Return the human readable name of a component (used in indices)
+ */
+ static std::string componentName(int compIdx)
+ {
+ switch (compIdx) {
+ case H2OIdx: return H2O::name();
+ case NAPLIdx: return MyCompressibleComponent::name();
+ };
+ DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
+ }
+
+ /*!
+ * \brief Return the molar mass of a component in [kg/mol].
+ */
+ static Scalar molarMass(int compIdx)
+ {
+ switch (compIdx) {
+ case H2OIdx: return H2O::molarMass();
+ case NAPLIdx: return MyCompressibleComponent::molarMass();
+ };
+ DUNE_THROW(Dune::InvalidStateException, "Invalid component index " << compIdx);
+ }
+
+ /*!
+ * \brief Given all mole fractions in a phase, return the phase
+ * density [kg/m^3].
+ */
+ using Base::density;
+ template <class FluidState>
+ static Scalar density(const FluidState &fluidState, int phaseIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+ if (phaseIdx == phase0Idx) {
+ // See: doctoral thesis of Steffen Ochs 2007
+ // Steam injection into saturated porous media : process analysis including experimental and numerical investigations
+ // http://elib.uni-stuttgart.de/bitstream/11682/271/1/Diss_Ochs_OPUS.pdf
+ Scalar rholH2O = H2O::liquidDensity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx));
+ Scalar clH2O = rholH2O/H2O::molarMass();
+ Scalar x_H2O = fluidState.moleFraction(phase0Idx, H2OIdx);
+ Scalar x_myComp = fluidState.moleFraction(phase0Idx, NAPLIdx);
+
+ // return composition-dependent water phase density
+ return clH2O*(H2O::molarMass()*x_H2O + MyCompressibleComponent::molarMass()*x_myComp);
+ }
+ else {
+ // assume the density of the fictious component to be independent of the composition
+ Scalar pressure = MyCompressibleComponent::liquidIsCompressible()?fluidState.pressure(phaseIdx):1e100;
+ return MyCompressibleComponent::liquidDensity(fluidState.temperature(phaseIdx), pressure);
+ }
+ }
+
+ using Base::molarDensity;
+ /*!
+ * \brief The molar density \f$\rho_{mol,\alpha}\f$
+ * of a fluid phase \f$\alpha\f$ in \f$\mathrm{[mol/m^3]}\f$
+ *
+ * The molar density for the simple relation is defined by the
+ * mass density \f$\rho_\alpha\f$ and the molar mass of the main component
+ *
+ * The molar density for the complrex relation is defined by the
+ * mass density \f$\rho_\alpha\f$ and the mean molar mass \f$\overline M_\alpha\f$:
+ *
+ * \f[\rho_{mol,\alpha} = \frac{\rho_\alpha}{\overline M_\alpha} \;.\f]
+ */
+ template <class FluidState>
+ static Scalar molarDensity(const FluidState &fluidState, int phaseIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+
+ Scalar T = fluidState.temperature(phaseIdx);
+ Scalar p = fluidState.pressure(phaseIdx);
+
+ // liquid phase
+ if (phaseIdx == phase0Idx)
+ {
+ // assume pure water or that each gas molecule displaces exactly one
+ // molecule in the liquid.
+ return H2O::liquidMolarDensity(T, p);
+ }
+ else
+ {
+ // assume the molar density of the fictious component to be independent of the composition
+ Scalar pressure = MyCompressibleComponent::liquidIsCompressible()?fluidState.pressure(phaseIdx):1e100;
+ return MyCompressibleComponent::liquidMolarDensity(fluidState.temperature(phaseIdx), pressure);
+ }
+ }
+
+ /*!
+ * \brief Return the viscosity of a phase.
+ */
+ using Base::viscosity;
+ template <class FluidState>
+ static Scalar viscosity(const FluidState &fluidState,
+ int phaseIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+ if (phaseIdx == phase0Idx) {
+ // assume pure water viscosity
+ return H2O::liquidViscosity(fluidState.temperature(phaseIdx),
+ fluidState.pressure(phaseIdx));
+ }
+ else {
+ // assume pure NAPL viscosity
+ return MyCompressibleComponent::liquidViscosity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx));
+ }
+ }
+
+ using Base::diffusionCoefficient;
+ template <class FluidState>
+ static Scalar diffusionCoefficient(const FluidState &fluidState, int phaseIdx, int compIdx)
+ {
+ DUNE_THROW(Dune::NotImplemented, "Diffusion coefficients");
+ }
+
+ /*!
+ * \brief Given a phase's composition, temperature and pressure,
+ * return the binary diffusion coefficient \f$\mathrm{[m^2/s]}\f$ for components
+ * \f$\mathrm{i}\f$ and \f$\mathrm{j}\f$ in this phase.
+ * \param fluidState The fluid state
+ * \param paramCache mutable parameters
+ * \param phaseIdx Index of the fluid phase
+ * \param compIIdx Index of the component i
+ * \param compJIdx Index of the component j
+ */
+ using Base::binaryDiffusionCoefficient;
+ template <class FluidState>
+ static Scalar binaryDiffusionCoefficient(const FluidState &fluidState,
+ int phaseIdx,
+ int compIIdx,
+ int compJIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+ assert(0 <= compIIdx && compIIdx < numComponents);
+ assert(0 <= compJIdx && compJIdx < numComponents);
+
+ const Scalar T = fluidState.temperature(phaseIdx);
+ const Scalar p = fluidState.pressure(phaseIdx);
+
+ // we assume the diffusion coefficient to be the same in both phases
+ return Dumux::BinaryCoeff::H2O_MyCompressibleComponent::liquidDiffCoeff(T, p);
+ }
+
+ /* Henry coefficients
+ */
+ template <class FluidState>
+ static Scalar henryCoefficient(const FluidState &fluidState,
+ int phaseIdx,
+ int compIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+ assert(0 <= compIdx && compIdx < numComponents);
+
+ const Scalar T = fluidState.temperature(phaseIdx);
+ const Scalar p = fluidState.pressure(phaseIdx);
+
+ if (compIdx == NAPLIdx && phaseIdx == phase0Idx)
+ return Dumux::BinaryCoeff::H2O_MyCompressibleComponent::henryMyCompressibleComponentInWater(T)/p;
+
+ else if (phaseIdx == phase1Idx && compIdx == H2OIdx)
+ return Dumux::BinaryCoeff::H2O_MyCompressibleComponent::henryWaterInMyCompressibleComponent(T)/p;
+
+ else
+ DUNE_THROW(Dune::InvalidStateException, "non-existent henry coefficient for phase index " << phaseIdx
+ << " and component index " << compIdx);
+ }
+
+ using Base::fugacityCoefficient;
+ /*!
+ * \brief Returns the fugacity coefficient \f$\mathrm{[-]}\f$ of a component in a
+ * phase.
+ *
+ * In this case, things are actually pretty simple. We have an ideal
+ * solution. Thus, the fugacity coefficient is 1 in the gas phase
+ * (fugacity equals the partial pressure of the component in the gas phase
+ * respectively in the liquid phases it is the inverse of the
+ * Henry coefficients scaled by pressure
+ * \param fluidState The fluid state
+ * \param phaseIdx The index of the phase
+ * \param compIdx The index of the component
+ */
+ template <class FluidState>
+ static Scalar fugacityCoefficient(const FluidState &fluidState,
+ int phaseIdx,
+ int compIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+ assert(0 <= compIdx && compIdx < numComponents);
+
+ Scalar T = fluidState.temperature(phaseIdx);
+ Scalar p = fluidState.pressure(phaseIdx);
+
+ if (phaseIdx == phase0Idx) {
+ if (compIdx == H2OIdx)
+ return H2O::vaporPressure(T)/p;
+ else if (compIdx == NAPLIdx)
+ return Dumux::BinaryCoeff::H2O_MyCompressibleComponent::henryMyCompressibleComponentInWater(T)/p;
+ }
+
+ // for the NAPL phase, we assume currently that nothing is
+ // dissolved. this means that the affinity of the NAPL
+ // component to the NAPL phase is much higher than for the
+ // other components, i.e. the fugacity coefficient is much
+ // smaller.
+ Scalar phiNapl = MyCompressibleComponent::vaporPressure(T)/p;
+ if (compIdx == NAPLIdx)
+ return phiNapl;
+ else
+ return 1e6*phiNapl;
+ }
+
+ template <class FluidState>
+ static Scalar kelvinVaporPressure(const FluidState &fluidState,
+ const int phaseIdx,
+ const int compIdx)
+ {
+ DUNE_THROW(Dune::NotImplemented, "FluidSystems::H2OMyCompressibleComponent::kelvinVaporPressure()");
+ }
+
+ /* partial pressures in the gas phase, taken from saturation vapor pressures
+ */
+ template <class FluidState>
+ static Scalar partialPressureGas(const FluidState &fluidState, int phaseIdx,
+ int compIdx)
+ {
+ assert(0 <= compIdx && compIdx < numComponents);
+
+ const Scalar T = fluidState.temperature(phaseIdx);
+ if (compIdx == NAPLIdx)
+ return MyCompressibleComponent::vaporPressure(T);
+ else if (compIdx == H2OIdx)
+ return H2O::vaporPressure(T);
+ else
+ DUNE_THROW(Dune::InvalidStateException, "non-existent component index " << compIdx);
+ }
+
+ /* inverse vapor pressures, taken from inverse saturation vapor pressures
+ */
+ template <class FluidState>
+ static Scalar inverseVaporPressureCurve(const FluidState &fluidState,
+ int phaseIdx,
+ int compIdx)
+ {
+ assert(0 <= compIdx && compIdx < numComponents);
+
+ const Scalar pressure = fluidState.pressure(phaseIdx);
+ if (compIdx == NAPLIdx)
+ return MyCompressibleComponent::vaporTemperature(pressure);
+ else if (compIdx == H2OIdx)
+ return H2O::vaporTemperature(pressure);
+ else
+ DUNE_THROW(Dune::InvalidStateException, "non-existent component index " << compIdx);
+ }
+
+ /*!
+ * \brief Given all mole fractions in a phase, return the specific
+ * phase enthalpy [J/kg].
+ */
+ using Base::enthalpy;
+ template <class FluidState>
+ static Scalar enthalpy(const FluidState &fluidState,
+ int phaseIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+ if (phaseIdx == phase0Idx) {
+ return H2O::liquidEnthalpy(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx));
+ }
+ else {
+ return MyCompressibleComponent::liquidEnthalpy(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx));
+ }
+ DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
+ }
+
+ using Base::heatCapacity;
+ template <class FluidState>
+ static Scalar heatCapacity(const FluidState &fluidState,
+ int phaseIdx)
+ {
+ DUNE_THROW(Dune::NotImplemented, "FluidSystems::H2ONAPL::heatCapacity()");
+ }
+
+ using Base::thermalConductivity;
+ template <class FluidState>
+ static Scalar thermalConductivity(const FluidState &fluidState,
+ int phaseIdx)
+ {
+ assert(0 <= phaseIdx && phaseIdx < numPhases);
+
+ const Scalar temperature = fluidState.temperature(phaseIdx) ;
+ const Scalar pressure = fluidState.pressure(phaseIdx);
+ if (phaseIdx == phase0Idx)
+ {
+ return H2O::liquidThermalConductivity(temperature, pressure);
+ }
+ else
+ {
+ return MyCompressibleComponent::liquidThermalConductivity(temperature, pressure);
+ }
+ DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
+ }
+
+private:
+
+};
+} // end namespace FluidSystems
+} // end namespace Dumux
+
+#endif
diff --git a/exercises/solution/ex3/mycompressiblecomponent.hh b/exercises/solution/ex3/mycompressiblecomponent.hh
new file mode 100644
index 0000000000000000000000000000000000000000..c19fae6d638510b14bcfb08cd7f924aa6dc3030c
--- /dev/null
+++ b/exercises/solution/ex3/mycompressiblecomponent.hh
@@ -0,0 +1,124 @@
+// -*- 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
+ * \ingroup Components
+ * \brief A fictitious component to be implemented in exercise 3.
+ */
+#ifndef DUMUX_MYCOMPRESSIBLECOMPONENT_HH
+#define DUMUX_MYCOMPRESSIBLECOMPONENT_HH
+
+#include <dumux/material/idealgas.hh>
+
+#include <dumux/material/components/base.hh>
+#include <dumux/material/components/liquid.hh>
+
+#include "myincompressiblecomponent.hh"
+
+namespace Dumux
+{
+/*!
+ * \ingroup Components
+ * \brief A fictitious component to be implemented in exercise 3.
+ *
+ * \tparam Scalar The type used for scalar values
+ */
+template <class Scalar>
+class MyCompressibleComponent
+: public Components::Base<Scalar, MyIncompressibleComponent<Scalar> >
+, public Components::Liquid<Scalar, MyIncompressibleComponent<Scalar> >
+{
+public:
+ /*!
+ * \brief A human readable name for MyCompressibleComponent.
+ */
+ static std::string name()
+ { return "MyCompressibleComponent"; }
+
+ /*!
+ * \brief The molar mass in \f$\mathrm{[kg/mol]}\f$ of MyCompressibleComponent.
+ */
+ static Scalar molarMass()
+ {
+ return 131.39e-3; // [kg/mol]
+ }
+
+ /*!
+ * \brief The density of MyCompressibleComponent at a given pressure and temperature \f$\mathrm{[kg/m^3]}\f$.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static Scalar liquidDensity(Scalar temperature, Scalar pressure)
+ {
+ static const Scalar rho_min = 1440;
+ static const Scalar rho_max = 1480;
+ static const Scalar k = 5e-7;
+
+ using std::exp;
+ return rho_min + (rho_max - rho_min)/(1 + rho_min*exp(-1.0*k*(rho_max - rho_min)*pressure)); // [kg/m^3]
+ }
+
+ /*!
+ * \brief The molar density of MyCompressibleComponent in \f$\mathrm{[mol/m^3]}\f$ at a given pressure and temperature.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ *
+ */
+ static Scalar liquidMolarDensity(Scalar temperature, Scalar pressure)
+ { return liquidDensity(temperature, pressure)/molarMass(); }
+
+
+ /*!
+ * \brief The dynamic viscosity \f$\mathrm{[Pa*s]}\f$ of MyCompressibleComponent.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static Scalar liquidViscosity(Scalar temperature, Scalar pressure)
+ {
+ return 5.7e-4;// [Pa*s]
+ }
+
+ /*****************************************************************
+ * The function below is implemented in the scope of exercise 3b *
+ *****************************************************************/
+
+ /*!
+ * \brief The vapor pressure in \f$\mathrm{[Pa]}\f$ of MyCompressibleComponent
+ * at a given temperature.
+ *
+ * \param T temperature of component in \f$\mathrm{[K]}\f$
+ */
+ static Scalar vaporPressure(Scalar T)
+ {
+ return 3900; // [Pa] (at 20C)
+ }
+
+ /*!
+ * \brief Returns true if the liquid phase is assumed to be compressible
+ */
+ static constexpr bool liquidIsCompressible()
+ { return true; }
+};
+
+} // end namespace
+
+#endif
diff --git a/exercises/solution/ex3/myincompressiblecomponent.hh b/exercises/solution/ex3/myincompressiblecomponent.hh
new file mode 100644
index 0000000000000000000000000000000000000000..24fe65567fb3b4d66acd2bb8bbe9d458fce9e892
--- /dev/null
+++ b/exercises/solution/ex3/myincompressiblecomponent.hh
@@ -0,0 +1,101 @@
+// -*- 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
+ * \ingroup Components
+ * \brief A fictitious component to be implemented in tutorial exercise 3.
+ */
+#ifndef DUMUX_MYINCOMPRESSIBLECOMPONENT_HH
+#define DUMUX_MYINCOMPRESSIBLECOMPONENT_HH
+
+#include <dumux/material/idealgas.hh>
+
+#include <dumux/material/components/base.hh>
+#include <dumux/material/components/liquid.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup Components
+ * \brief A fictitious component to be implemented in exercise 3.
+ *
+ * \tparam Scalar The type used for scalar values
+ */
+template <class Scalar>
+class MyIncompressibleComponent
+: public Components::Base<Scalar, MyIncompressibleComponent<Scalar> >
+, public Components::Liquid<Scalar, MyIncompressibleComponent<Scalar> >
+{
+public:
+ /*!
+ * \brief A human readable name for MyIncompressibleComponent.
+ */
+ static std::string name()
+ { return "MyIncompressibleComponent"; }
+
+ /*!
+ * \brief The molar mass in \f$\mathrm{[kg/mol]}\f$ of MyIncompressibleComponent.
+ */
+ static Scalar molarMass()
+ {
+ return 131.39e-3; // [kg/mol]
+ }
+
+ /*!
+ * \brief The density of MyIncompressibleComponent at a given pressure and temperature \f$\mathrm{[kg/m^3]}\f$.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static Scalar liquidDensity(Scalar temperature, Scalar pressure)
+ {
+ return 1460.0; // [kg/m^3]
+ }
+
+ /*!
+ * \brief The molar density of MyIncompressibleComponent in \f$\mathrm{[mol/m^3]}\f$ at a given pressure and temperature.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ *
+ */
+ static Scalar liquidMolarDensity(Scalar temperature, Scalar pressure)
+ { return liquidDensity(temperature, pressure)/molarMass(); }
+
+ /*!
+ * \brief The dynamic viscosity \f$\mathrm{[Pa*s]}\f$ of MyIncompressibleComponent.
+ *
+ * \param temperature temperature of component in \f$\mathrm{[K]}\f$
+ * \param pressure pressure of component in \f$\mathrm{[Pa]}\f$
+ */
+ static Scalar liquidViscosity(Scalar temperature, Scalar pressure)
+ {
+ return 5.7e-4;// [Pa*s]
+ }
+
+ /*!
+ * \brief Returns true if the liquid phase is assumed to be compressible
+ */
+ static constexpr bool liquidIsCompressible()
+ { return false; }
+};
+
+} // end namespace
+
+#endif
diff --git a/scripts/README.md b/scripts/README.md
index 8bf9b20837e8ad92a988272f15f59776eb762177..6f87492760a9f7936137cc9bf9ce99f23a4d130e 100644
--- a/scripts/README.md
+++ b/scripts/README.md
@@ -13,21 +13,19 @@ You need to have the following REQUIREMENTS installed:
* git
* pkg-config
* paraview (to visualize the results)
+ * gnuplot (to plot some curves)
+ * wget (to download some config files during the installation)
On debian-based system you can use this:
- apt-get install build-essential pkg-config cmake git paraview
+ apt-get install build-essential pkg-config cmake git paraview wget gnuplot
Then, you can the execute the script and it will download the dune repositories and dumux
and configure all modules with CMake
./install.sh
-This will clone the necessary repositories (in a subfolder 'dune'),
+This will clone the necessary repositories (in a subfolder 'dumux'),
build all libaries.
-The dune core modules are also available as debian packages
-(see http://www.dune-project.org/download.html#binary) where
-the version 2.4 is required.
-
-Run the script test_dumux.sh in the same folder you ran this script
+Run the script test_dumux.sh in the newly created dumux folder
to test your installation of dumux.
./test_dumux.sh
It will compile and run a simple one-phase ground water flow example
diff --git a/scripts/install.sh b/scripts/install.sh
index bdad6681f74d0a9586e33bbb62aed88720bfb6c3..c1c7b23392ab4aedb300a6017dc149664eeabc0c 100755
--- a/scripts/install.sh
+++ b/scripts/install.sh
@@ -2,7 +2,7 @@
DUNE_VERSION=2.6
# check some prerequistes
-for PRGRM in git cmake gcc g++ wget paraview pkg-config; do
+for PRGRM in git cmake gcc g++ wget paraview pkg-config gnuplot; do
if ! [ -x "$(command -v $PRGRM)" ]; then
echo "Error: $PRGRM is not installed." >&2
exit 1
@@ -29,9 +29,14 @@ echo "(0/2) Downloading supplementary files. Make sure to be connected to the in
echo "*********************************************************************************************"
# download the install.opts and the test script
-wget https://git.iws.uni-stuttgart.de/dumux-repositories/dumux-course/raw/master/scripts/test_dumux.sh
-chmod +x test_dumux.sh
-wget https://git.iws.uni-stuttgart.de/dumux-repositories/dumux-course/raw/master/cmake.opts
+if [ ! -f "test_dumux.sh" ]; then
+ wget https://git.iws.uni-stuttgart.de/dumux-repositories/dumux-course/raw/master/scripts/test_dumux.sh
+ chmod +x test_dumux.sh
+fi
+if [ ! -f "cmake.opts" ]; then
+ wget https://git.iws.uni-stuttgart.de/dumux-repositories/dumux-course/raw/master/cmake.opts
+ chmod +x test_dumux.sh
+fi
# get the testing script
@@ -64,11 +69,11 @@ done
# dune-subgrid
if [ ! -d "dune-subgrid" ]; then
- git clone -b releases/$DUNE_VERSION https://git.imp.fu-berlin.de/agnumpde/dune-subgrid.git
+ git clone -b releases/$DUNE_VERSION-1 https://git.imp.fu-berlin.de/agnumpde/dune-subgrid.git
else
echo "Skip cloning dune-subgrid because the folder already exists."
cd dune-subgrid
- git checkout releases/$DUNE_VERSION
+ git checkout releases/$DUNE_VERSION-1
cd ..
fi