diff --git a/exercises/solution/exercise-runtimeparams/exercise_runtimeparams_solution.cc b/exercises/solution/exercise-runtimeparams/exercise_runtimeparams_solution.cc new file mode 100644 index 0000000000000000000000000000000000000000..5356a47778b40b14dedb9e6102b7c6d8b20e3870 --- /dev/null +++ b/exercises/solution/exercise-runtimeparams/exercise_runtimeparams_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 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; +} diff --git a/exercises/solution/exercise-runtimeparams/exercise_runtimeparams_solution.input b/exercises/solution/exercise-runtimeparams/exercise_runtimeparams_solution.input new file mode 100644 index 0000000000000000000000000000000000000000..62cbcc86ab3b49b6471fff0470fa3dbe3e4f04e3 --- /dev/null +++ b/exercises/solution/exercise-runtimeparams/exercise_runtimeparams_solution.input @@ -0,0 +1,22 @@ +[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 diff --git a/exercises/solution/exercise-runtimeparams/injection2pproblem.hh b/exercises/solution/exercise-runtimeparams/injection2pproblem.hh new file mode 100644 index 0000000000000000000000000000000000000000..74c92c3e154dffea18c4ec00e39d8790853ca090 --- /dev/null +++ b/exercises/solution/exercise-runtimeparams/injection2pproblem.hh @@ -0,0 +1,272 @@ +// -*- 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 diff --git a/exercises/solution/exercise-runtimeparams/injection2pspatialparams.hh b/exercises/solution/exercise-runtimeparams/injection2pspatialparams.hh new file mode 100644 index 0000000000000000000000000000000000000000..b42f8114cf4ac0c6ddd07cbb083d868f6c1c0976 --- /dev/null +++ b/exercises/solution/exercise-runtimeparams/injection2pspatialparams.hh @@ -0,0 +1,172 @@ +// -*- 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