diff --git a/exercises/exercise-basic/exercise1_2p.cc b/exercises/exercise-basic/exercise1_2p.cc
index 558f2d0f3373ed23df1fefbb0af5c296a698484a..fbefdd1eae0db5ea4225dd120be224f7cf86f615 100644
--- a/exercises/exercise-basic/exercise1_2p.cc
+++ b/exercises/exercise-basic/exercise1_2p.cc
@@ -101,8 +101,6 @@ int main(int argc, char** argv) try
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");
diff --git a/exercises/exercise-basic/exercise1_2p2c.cc b/exercises/exercise-basic/exercise1_2p2c.cc
index c2b46ea15fcab30a078f85acf0f8d07b122b8209..f24009436f37aa1199b4cddacefa4189bdbb12b8 100644
--- a/exercises/exercise-basic/exercise1_2p2c.cc
+++ b/exercises/exercise-basic/exercise1_2p2c.cc
@@ -101,8 +101,6 @@ int main(int argc, char** argv) try
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");
diff --git a/exercises/exercise-basic/injection2p2cproblem.hh b/exercises/exercise-basic/injection2p2cproblem.hh
index 26e111092c6ff2d19ed587c3c2609a6a40389eb8..12535cbe8f06d2247e79993414171414740f7151 100644
--- a/exercises/exercise-basic/injection2p2cproblem.hh
+++ b/exercises/exercise-basic/injection2p2cproblem.hh
@@ -118,10 +118,6 @@ public:
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_.
-
}
/*!
diff --git a/exercises/exercise-basic/injection2pspatialparams.hh b/exercises/exercise-basic/injection2pspatialparams.hh
index fac5a1da9896f747dc94b9b9787383976c063f86..fa328bc1a82ecded90d95e7d67419c924f8a71ae 100644
--- a/exercises/exercise-basic/injection2pspatialparams.hh
+++ b/exercises/exercise-basic/injection2pspatialparams.hh
@@ -126,7 +126,7 @@ public:
*
* \return the material parameters object
*/
- const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const
+ const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const
{
if (isInAquitard_(globalPos))
return aquitardMaterialParams_;
diff --git a/exercises/solution/ex1/exercise1.input b/exercises/solution/ex1/exercise1.input
index 545a170fb594b33f1a0600be18d50495c3953ca0..ada32b39acafe33819dfff1518deb878b2398790 100644
--- a/exercises/solution/ex1/exercise1.input
+++ b/exercises/solution/ex1/exercise1.input
@@ -12,7 +12,6 @@ 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
diff --git a/exercises/solution/ex1/exercise1_2p2c.cc b/exercises/solution/ex1/exercise1_2p2c.cc
deleted file mode 100644
index 55496aad06bae940e473fa8ad013a27d23be5444..0000000000000000000000000000000000000000
--- a/exercises/solution/ex1/exercise1_2p2c.cc
+++ /dev/null
@@ -1,198 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * See the file COPYING for full copying permissions. *
- * *
- * This program is free software: you can redistribute it and/or modify *
- * it under the terms of the GNU General Public License as published by *
- * the Free Software Foundation, either version 2 of the License, or *
- * (at your option) any later version. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * *
- * You should have received a copy of the GNU General Public License *
- * along with this program. If not, see <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- * \brief 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/injection2p2cproblem.hh b/exercises/solution/ex1/injection2p2cproblem.hh
deleted file mode 100644
index d33e5a46bb25ea3bffec3f54a804af1947d5f7d9..0000000000000000000000000000000000000000
--- a/exercises/solution/ex1/injection2p2cproblem.hh
+++ /dev/null
@@ -1,271 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * See the file COPYING for full copying permissions. *
- * *
- * This program is free software: you can redistribute it and/or modify *
- * it under the terms of the GNU General Public License as published by *
- * the Free Software Foundation, either version 2 of the License, or *
- * (at your option) any later version. *
- * *
- * This program is distributed in the hope that it will be useful, *
- * but WITHOUT ANY WARRANTY; without even the implied warranty of *
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
- * GNU General Public License for more details. *
- * *
- * You should have received a copy of the GNU General Public License *
- * along with this program. If not, see <http://www.gnu.org/licenses/>. *
- *****************************************************************************/
-/*!
- * \file
- *
- * \brief 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
index 70135414a985fcfa0bbf4efdd75ccf0e9306763e..6cc7132b280e6408639f2eeffd114e39a2bd624d 100644
--- a/exercises/solution/ex1/injection2pniproblem.hh
+++ b/exercises/solution/ex1/injection2pniproblem.hh
@@ -90,6 +90,7 @@ class InjectionProblem2PNI : public PorousMediumFlowProblem<TypeTag>
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);
+ using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
enum { dimWorld = GridView::dimensionworld };
using Element = typename GridView::template Codim<0>::Entity;
@@ -204,6 +205,22 @@ public:
*/
// \{
+ /*!
+ * \brief Evaluate the source term for all phases within a given
+ * sub-control-volume.
+ *
+ * For this method, the \a priVars parameter stores the rate mass
+ * of a component is generated or annihilate per volume
+ * unit. Positive values mean that mass is created, negative ones
+ * mean that it vanishes.
+ *
+ * The units must be according to either using mole or mass fractions. (mole/(m^3*s) or kg/(m^3*s))
+ */
+ NumEqVector sourceAtPos(const GlobalPosition &globalPos) const
+ {
+ return NumEqVector(0.0);
+ }
+
/*!
* \brief Evaluate the initial value for a control volume.
*
diff --git a/exercises/solution/ex1/injection2pspatialparams.hh b/exercises/solution/ex1/injection2pspatialparams.hh
index 979f14ce9782f08b1a08010b0f8053506265ebf3..fa328bc1a82ecded90d95e7d67419c924f8a71ae 100644
--- a/exercises/solution/ex1/injection2pspatialparams.hh
+++ b/exercises/solution/ex1/injection2pspatialparams.hh
@@ -49,6 +49,8 @@ class InjectionSpatialParams
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;
@@ -65,7 +67,7 @@ public:
*
* \param fvGridGeometry The finite volume grid geometry
*/
- InjectionSpatialParams(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+ InjectionSpatialParams(std::shared_ptr<const FVGridGeometry>& fvGridGeometry)
: ParentType(fvGridGeometry)
{
aquiferHeightFromBottom_ = 30.0;
@@ -97,8 +99,8 @@ public:
* \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_;
@@ -111,6 +113,7 @@ public:
*/
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_;
@@ -144,8 +147,12 @@ 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
- { return globalPos[dimWorld-1] > aquiferHeightFromBottom_ + eps_; }
+ {
+ // 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_;