exercise-runtimeparams solution files

exercises/solution/exercise-runtimeparams/exercise_runtimeparams_solution.cc

+// -*- 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 runtime parameters
+ */
+#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;
+}

exercises/solution/exercise-runtimeparams/exercise_runtimeparams_solution.input

+[TimeLoop]
+DtInitial = 3600 # in seconds
+TEnd = 3.154e9 # in seconds, i.e ten years
+
+[Grid]
+UpperRight = 60 40
+Cells = 24 16
+
+[Problem]
+Name = runtimeparams_exercise
+OnlyPlotMaterialLaws = true
+AquiferDepth = 2700.0 # m
+InjectionDuration = 2.628e6 # in seconds, i.e. one month
+# TODO: Task 2: Create a parameter called "TotalAreaSpecificInflow"
+TotalAreaSpecificInflow = -1e-4 # kg/(s m^2)
+
+[SpatialParams]
+PermeabilityAquitard = 1e-15 # m^2
+# TODO: Task 1: Change the Aquitard's Entry Pressure
+EntryPressureAquitard = 4.5e3 # Pa
+PermeabilityAquifer = 1e-12 # m^2
+EntryPressureAquifer = 1e4 # Pa

exercises/solution/exercise-runtimeparams/injection2pproblem.hh

+// -*- 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 runtime parameters
+ */
+
+#ifndef DUMUX_EXRUNTIMEPARAMS_INJECTION_PROBLEM_2P_HH
+#define DUMUX_EXRUNTIMEPARAMS_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_ = getParamFromGroup<Scalar>("Problem","InjectionDuration");
+        //TODO: Task 2: Set a variable "TotalAreaSpecificInflow" to read in a value from the parameter tree via the input file
+        totalAreaSpecificInflow_ = getParam<Scalar>("Problem.TotalAreaSpecificInflow");
+        //TODO: Task 3: Set a default value for the above parameter.
+//        totalAreaSpecificInflow_ = getParam<Scalar>("Problem.TotalAreaSpecificInflow", -1e-4);
+        //TODO: Task 4: Provide output describing where the parameter value comes from using parameter bool functions.
+//         if (hasParamInGroup("Problem","TotalAreaSpecificInflow"))
+//             std::cout << "Parameter value is read from the input file." << std::endl;
+//         else
+//             std::cout << "Using the default parameter value." << std::endl;
+    }
+
+    /*!
+     * \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)
+            //TODO: Task 2: incorporate "totalAreaSpecificInflow_"  into the injection boundary condition
+            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);
+
+        // 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
+    //TODO: Task 2: Set a variable "totalAreaSpecificInflow_" to read in a value from the parameter tree via the input file
+    Scalar totalAreaSpecificInflow_; //! Rate of the Injection in kg/(s m^2)
+    Scalar time_;
+};
+
+} //end namespace Dumux
+
+#endif

exercises/solution/exercise-runtimeparams/injection2pspatialparams.hh

+// -*- 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_RUNTIMEPARAMS_INJECTION_SPATIAL_PARAMS_HH
+#define DUMUX_RUNTIMEPARAMS_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)
+    {
+        // Aquifer Height, measured from the bottom
+        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