From dfac295498722a5c8b6379c69b7034ef76fec215 Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Wed, 20 Sep 2017 19:52:18 +0200
Subject: [PATCH] Add some stuff for exercise 2
---
tutorial/ex2/CMakeLists.txt | 6 +
tutorial/ex2/README.md | 128 ++++++++
tutorial/ex2/exercise2.cc | 67 +++++
tutorial/ex2/exercise2.input | 18 ++
tutorial/ex2/injection2p2cproblem.hh | 367 +++++++++++++++++++++++
tutorial/ex2/injection2pspatialparams.hh | 253 ++++++++++++++++
6 files changed, 839 insertions(+)
create mode 100644 tutorial/ex2/CMakeLists.txt
create mode 100644 tutorial/ex2/README.md
create mode 100644 tutorial/ex2/exercise2.cc
create mode 100755 tutorial/ex2/exercise2.input
create mode 100644 tutorial/ex2/injection2p2cproblem.hh
create mode 100644 tutorial/ex2/injection2pspatialparams.hh
diff --git a/tutorial/ex2/CMakeLists.txt b/tutorial/ex2/CMakeLists.txt
new file mode 100644
index 0000000000..d53f775c2a
--- /dev/null
+++ b/tutorial/ex2/CMakeLists.txt
@@ -0,0 +1,6 @@
+# a compositional two-phase simulation program
+dune_add_test(NAME exercise2
+ SOURCES exercise2.cc)
+
+# add a symlink for the input file
+dune_symlink_to_source_files(FILES "exercise2.input")
diff --git a/tutorial/ex2/README.md b/tutorial/ex2/README.md
new file mode 100644
index 0000000000..49de2542c7
--- /dev/null
+++ b/tutorial/ex2/README.md
@@ -0,0 +1,128 @@
+# Exercise #2 (DuMuX course)
+
+## Problem set-up
+
+The problem setup is identical to the previous [_exercise 1_](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux/blob/6c2ca5dd000599f793f89fe5e3fc69ef9c9d8b73/tutorial/ex1/README.md) with a lower injection rate of 1e-6 kg/(m*s) so that diffusion plays a more dominant role in the transport process.
+
+## Preparing the exercise
+
+* Navigate to the directory `dumux/tutorial/dumux-course`
+
+_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_.
+
+### 1. Getting familiar with the code
+
+Locate all the files you will need for this exercise
+* The __main file__: `exercise2.cc`
+* The __problem file__: `injection2p2cproblem.hh`
+* The __spatial parameters file__: `injection2pspatialparams.hh`
+* The __input file__: `exercise2.input`
+* Two header files containing:
+ * a custom __local residual__ in: `mylocalresidual.hh`
+ * a custom __material law__ in:: `mymateriallaw.hh`
+
+### 2. Compiling and running the program
+
+* Change to the build-directory
+
+```bash
+cd build-cmake/tutorial/exercise2
+```
+
+* Compile the executable `exercise2`
+
+```bash
+make exercise2
+```
+
+* Execute the two problems 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.
+
+### 3. Implement and use a different material law
+
+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(Scalar sw)` by implementing your own capillary pressure relationship, e.g. pc(Sw) = 1e5*(1-Sw).
+
+The type (i.e. C++ type) of the material law is set in the file `injection2pspatialparams.hh` by using the DuMuX property system
+
+```c++
+SET_PROP(InjectionSpatialParams, MaterialLaw)
+{
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using type = EffToAbsLaw<RegularizedBrooksCorey<Scalar>>;
+};
+```
+
+* Make DuMuX use your own material law by including the header `mymateriallaw.hh` and changing the alias `type`. This will make sure that your material law is used everywhere else in the code.
+
+Verify your changes by recompiling and running the program. You should see a plot of your new function.
+
+### 4. Implement your own local residual
+
+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 step. The file `mylocalresidual.hh` contains a copy of the original local residual class used for the 2p2c model renamed to `template<class TypeTag> class MyTwoPTwoCLocalResidual`.
+
+* Make DuMuX use this new local residual by setting the corresponding property in the `Property` namespace in the file `injection2p2cproblem.hh`
+
+```c++
+// note that every property struct knows about TypeTag
+SET_PROP(Injection2p2cProblem, LocalResidual)
+{
+ using type = MyTwoPTwoCLocalResidualocal<TypeTag>;
+};
+
+// or using the convenience macro
+SET_TYPE_PROP(Injection2p2cProblem, LocalResidual,
+ MyTwoPTwoCLocalResidualocal<TypeTag>);
+```
+
+* Implement an output to the terminal in the constructor of `MyTwoPTwoCLocalResidual` e.g.
+
+```c++
+MyTwoPTwoCLocalResidual()
+{
+ std::cout << "Using MyTwoPTwoCLocalResidual." << std::endl;
+}
+```
+
+* Verify you are using the new class by compiling and running the new program and inspecting the terminal output.
+
+You want to make the new local residual special by adding a switch enabling / disabling diffusion. We will achieve this with a DuMuX parameter, a parameter read from the input file that defaults to a property value if the input file doesn't contain the parameter.
+
+* Create a new `TypeTag` node, a new `PropertyTag`, and set a default in the `mylocalresidual.hh` file by adding
+
+```c++
+namespace Dumux {
+
+namespace Properties
+{
+ NEW_TYPE_TAG(MyLocalResidualParams); // creates a new TypeTag node
+ NEW_PROP_TAG(ProblemEnableDiffusion); // creates a new property
+ SET_BOOL_PROP(MyLocalResidualParams,
+ ProblemEnableDiffusion, true); // set a default value
+}
+...
+```
+
+* 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 constructor of MyTwoPTwoCLocalResidual
+ enableDiffusion_ = GET_PARAM_FROM_GROUP(TypeTag, bool,
+ Problem, EnableDiffusion);
+
+// ... in the private member section of MyTwoPTwoCLocalResidual
+private:
+ bool enableDiffusion_;
+```
diff --git a/tutorial/ex2/exercise2.cc b/tutorial/ex2/exercise2.cc
new file mode 100644
index 0000000000..aba0b0e8c7
--- /dev/null
+++ b/tutorial/ex2/exercise2.cc
@@ -0,0 +1,67 @@
+// -*- 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 "injection2p2cproblem.hh"
+#include <dumux/common/start.hh>
+
+/*!
+ * \brief Provides an interface for customizing error messages associated with
+ * reading in parameters.
+ *
+ * \param progName The name of the program, that was tried to be started.
+ * \param errorMsg The error message that was issued by the start function.
+ * Comprises the thing that went wrong and a general help message.
+ */
+void usage(const char *progName, const std::string &errorMsg)
+{
+ if (errorMsg.size() > 0) {
+ std::string errorMessageOut = "\nUsage: ";
+ errorMessageOut += progName;
+ errorMessageOut += " [options]\n";
+ errorMessageOut += errorMsg;
+ errorMessageOut += "\n\nThe list of mandatory 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-FluidSystem.NTemperature Number of tabularization entries [-] \n"
+ "\t-FluidSystem.NPressure Number of tabularization entries [-] \n"
+ "\t-FluidSystem.PressureLow Low end for tabularization of fluid properties [Pa] \n"
+ "\t-FluidSystem.PressureHigh High end for tabularization of fluid properties [Pa] \n"
+ "\t-FluidSystem.TemperatureLow Low end for tabularization of fluid properties [Pa] \n"
+ "\t-FluidSystem.TemperatureHigh High end for tabularization of fluid properties [Pa] \n"
+ "\t-SimulationControl.Name The name of the output files [-] \n"
+ "\t-InitialConditions.Temperature Initial temperature in the reservoir [K] \n"
+ "\t-InitialConditions.DepthBOR Depth below ground surface [m] \n";
+
+ std::cout << errorMessageOut
+ << "\n";
+ }
+}
+
+int main(int argc, char** argv)
+{
+ typedef TTAG(Injection2p2pcCCProblem) ProblemTypeTag;
+ return Dumux::start<ProblemTypeTag>(argc, argv, usage);
+}
diff --git a/tutorial/ex2/exercise2.input b/tutorial/ex2/exercise2.input
new file mode 100755
index 0000000000..a4e055be8c
--- /dev/null
+++ b/tutorial/ex2/exercise2.input
@@ -0,0 +1,18 @@
+[TimeManager]
+DtInitial = 250 # [s]
+TEnd = 1e7 # [s]
+
+[Grid]
+LowerLeft = 0 0
+UpperRight = 60 40
+Cells = 24 16
+
+[Problem]
+MaxDepth = 2700.0
+
+[Newton]
+WriteConvergence = 1 # write convergence behaviour to disk?
+
+[SpatialParams]
+EntryPressureFine = 4.5e4
+EntryPressureCoarse = 1e4
diff --git a/tutorial/ex2/injection2p2cproblem.hh b/tutorial/ex2/injection2p2cproblem.hh
new file mode 100644
index 0000000000..a63bc73d1c
--- /dev/null
+++ b/tutorial/ex2/injection2p2cproblem.hh
@@ -0,0 +1,367 @@
+// -*- 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/implicit/model.hh>
+#include <dumux/porousmediumflow/implicit/problem.hh>
+#include <dumux/material/fluidsystems/h2on2.hh>
+
+#include "injection2pspatialparams.hh"
+
+namespace Dumux
+{
+
+template <class TypeTag>
+class Injection2p2cProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(Injection2p2cProblem, INHERITS_FROM(TwoPTwoC, InjectionSpatialParams));
+NEW_TYPE_TAG(Injection2p2cBoxProblem, INHERITS_FROM(BoxModel, Injection2p2cProblem));
+NEW_TYPE_TAG(Injection2p2pcCCProblem, INHERITS_FROM(CCModel, Injection2p2cProblem));
+
+// Set the grid type
+SET_TYPE_PROP(Injection2p2cProblem, Grid, Dune::YaspGrid<2>);
+
+// Set the problem property
+SET_TYPE_PROP(Injection2p2cProblem, Problem, Injection2p2cProblem<TypeTag>);
+
+// Set fluid configuration
+SET_TYPE_PROP(Injection2p2cProblem,
+ FluidSystem,
+ FluidSystems::H2ON2<typename GET_PROP_TYPE(TypeTag, Scalar), false /*useComplexRelations*/>);
+
+// Define whether mole(true) or mass (false) fractions are used
+SET_BOOL_PROP(Injection2p2cProblem, UseMoles, true);
+}
+
+
+/*!
+ * \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 (\f$ K=10e-12\f$) for \f$ y<22m\f$ and one with a lower permeablility (\f$ K=10e-13\f$)
+ * in the rest of the domain.
+ *
+ * A mixture of Nitrogen and Water vapor, which is composed according to the prevailing conditions (temperature, pressure)
+ * enters a water-filled aquifer. This is realized with a solution-dependent Neumann boundary condition at the right boundary
+ * (\f$ 7m<y<15m\f$). The aquifer is situated 2700m below sea level. The injected fluid phase migrates upwards due to buoyancy.
+ * It accumulates and partially enters the lower permeable aquitard.
+ *
+ * The model is able to use either mole or mass fractions. The property useMoles can be set to either true or false in the
+ * problem file. Make sure that the according units are used in the problem setup. The default setting for useMoles is true.
+ *
+ * This problem uses the \ref TwoPTwoCModel.
+ *
+ * To run the simulation execute the following line in shell:
+ * <tt>./exercise1_2p2c -parameterFile exercise1.input>
+ */
+template <class TypeTag>
+class Injection2p2cProblem : public ImplicitPorousMediaProblem<TypeTag>
+{
+ typedef ImplicitPorousMediaProblem<TypeTag> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+
+ enum {
+ // Grid and world dimension
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld
+ };
+
+ // copy some indices for convenience
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+ enum {
+ wPhaseIdx = Indices::wPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx,
+
+
+ wCompIdx = FluidSystem::wCompIdx,
+ nCompIdx = FluidSystem::nCompIdx,
+
+ contiH2OEqIdx = Indices::contiWEqIdx,
+ contiN2EqIdx = Indices::contiNEqIdx
+ };
+
+
+ typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
+ typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::template Codim<dim>::Entity Vertex;
+ typedef typename GridView::Intersection Intersection;
+
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+ enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+
+ //! property that defines whether mole or mass fractions are used
+ static const bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
+
+public:
+ /*!
+ * \brief The constructor
+ *
+ * \param timeManager The time manager
+ * \param gridView The grid view
+ */
+ Injection2p2cProblem(TimeManager &timeManager,
+ const GridView &gridView)
+ : ParentType(timeManager, gridView)
+ {
+ maxDepth_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.MaxDepth);
+
+ // initialize the tables of the fluid system
+ FluidSystem::init(/*tempMin=*/273.15,
+ /*tempMax=*/423.15,
+ /*numTemp=*/50,
+ /*pMin=*/0.0,
+ /*pMax=*/30e6,
+ /*numP=*/300);
+
+ //stateing in the console whether mole or mass fractions are used
+ if(useMoles)
+ {
+ std::cout<<"problem uses mole-fractions"<<std::endl;
+ }
+ else
+ {
+ std::cout<<"problem uses mass-fractions"<<std::endl;
+ }
+ }
+
+ /*!
+ * \brief User defined output after the time integration
+ *
+ * Will be called diretly after the time integration.
+ */
+ void postTimeStep()
+ {
+ // Calculate storage terms
+ PrimaryVariables storageW, storageN;
+ this->model().globalPhaseStorage(storageW, wPhaseIdx);
+ this->model().globalPhaseStorage(storageN, nPhaseIdx);
+
+ // Write mass balance information for rank 0
+ if (this->gridView().comm().rank() == 0) {
+ std::cout<<"Storage: wetting=[" << storageW << "]"
+ << " nonwetting=[" << storageN << "]\n";
+ }
+ }
+
+
+ /*!
+ * \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 "injection-2p2c"; }
+
+ /*!
+ * \brief Returns the temperature \f$ K \f$
+ */
+ Scalar temperature() const
+ {
+ return 273.15 + 30; // [K]
+ }
+
+ /*!
+ * \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
+ */
+ void sourceAtPos(PrimaryVariables &values,
+ const GlobalPosition &globalPos) const
+ {
+ values = 0;
+ }
+
+ // \}
+
+ /*!
+ * \name Boundary conditions
+ */
+ // \{
+
+ /*!
+ * \brief Specifies which kind of boundary condition should be
+ * used for which equation on a given boundary segment
+ *
+ * \param values Stores the value of the boundary type
+ * \param globalPos The global position
+ */
+ void boundaryTypesAtPos(BoundaryTypes &values,
+ const GlobalPosition &globalPos) const
+ {
+ if (globalPos[0] < eps_)
+ values.setAllDirichlet();
+ else
+ values.setAllNeumann();
+ }
+
+ /*!
+ * \brief Evaluates the boundary conditions for a Dirichlet
+ * boundary segment
+ *
+ * \param values Stores the Dirichlet values for the conservation equations in
+ * \f$ [ \textnormal{unit of primary variable} ] \f$
+ * \param globalPos The global position
+ */
+ void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
+ {
+ initial_(values, globalPos);
+ }
+
+ /*!
+ * \brief Evaluates the boundary conditions for a Neumann
+ * boundary segment in dependency on the current solution.
+ *
+ * \param values Stores the Neumann values for the conservation equations in
+ * \f$ [ \textnormal{unit of conserved quantity} / (m^(dim-1) \cdot s )] \f$
+ * \param element The finite element
+ * \param fvGeometry The finite volume geometry of the element
+ * \param intersection The intersection between element and boundary
+ * \param scvIdx The local index of the sub-control volume
+ * \param boundaryFaceIdx The index of the boundary face
+ * \param elemVolVars All volume variables for the element
+ *
+ * This method is used for cases, when the Neumann condition depends on the
+ * solution and requires some quantities that are specific to the fully-implicit method.
+ * The \a values store the mass flux of each phase normal to the boundary.
+ * Negative values indicate an inflow.
+ */
+ void solDependentNeumann(PrimaryVariables &values,
+ const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const Intersection &intersection,
+ const int scvIdx,
+ const int boundaryFaceIdx,
+ const ElementVolumeVariables &elemVolVars) const
+ {
+ values = 0;
+
+ GlobalPosition globalPos;
+ if (isBox)
+ globalPos = element.geometry().corner(scvIdx);
+ else
+ globalPos = intersection.geometry().center();
+
+ Scalar injectedPhaseMass = 1e-3;
+ Scalar moleFracW = elemVolVars[scvIdx].moleFraction(nPhaseIdx, wCompIdx);
+ if (globalPos[1] < 15 + eps_ && globalPos[1] > 7 -eps_) {
+ values[contiN2EqIdx] = -(1-moleFracW)*injectedPhaseMass/FluidSystem::molarMass(nCompIdx); //mole/(m^2*s) -> kg/(s*m^2)
+ values[contiH2OEqIdx] = -moleFracW*injectedPhaseMass/FluidSystem::molarMass(wCompIdx); //mole/(m^2*s) -> kg/(s*m^2)
+ }
+ }
+
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ /*!
+ * \brief Evaluates the initial values for a control volume
+ *
+ * \param values Stores the initial values for the conservation equations in
+ * \f$ [ \textnormal{unit of primary variables} ] \f$
+ * \param globalPos The global position
+ */
+ void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
+ {
+ initial_(values, globalPos);
+ }
+
+ /*!
+ * \brief Return the initial phase state inside a control volume.
+ *
+ * \param globalPos The global position
+ */
+ int initialPhasePresenceAtPos(const GlobalPosition &globalPos) const
+ { return Indices::wPhaseOnly; }
+
+ // \}
+
+private:
+ /*!
+ * \brief Evaluates the initial values for a control volume
+ *
+ * The internal method for the initial condition
+ *
+ * \param values Stores the initial values for the conservation equations in
+ * \f$ [ \textnormal{unit of primary variables} ] \f$
+ * \param globalPos The global position
+ */
+ void initial_(PrimaryVariables &values,
+ const GlobalPosition &globalPos) const
+ {
+ Scalar densityW = FluidSystem::H2O::liquidDensity(temperature(), 1e5);
+
+ Scalar pl = 1e5 - densityW*this->gravity()[1]*(maxDepth_ - globalPos[1]);
+ Scalar moleFracLiquidN2 = pl*0.95/BinaryCoeff::H2O_N2::henry(temperature());
+ Scalar moleFracLiquidH2O = 1.0 - moleFracLiquidN2;
+
+ Scalar meanM =
+ FluidSystem::molarMass(wCompIdx)*moleFracLiquidH2O +
+ FluidSystem::molarMass(nCompIdx)*moleFracLiquidN2;
+ if(useMoles)
+ {
+ //mole-fraction formulation
+ values[Indices::switchIdx] = moleFracLiquidN2;
+ }
+ else
+ {
+ //mass fraction formulation
+ Scalar massFracLiquidN2 = moleFracLiquidN2*FluidSystem::molarMass(nCompIdx)/meanM;
+ values[Indices::switchIdx] = massFracLiquidN2;
+ }
+ values[Indices::pressureIdx] = pl;
+ }
+
+
+ static constexpr Scalar eps_ = 1e-6;
+ Scalar maxDepth_;
+
+};
+} //end namespace
+
+#endif
diff --git a/tutorial/ex2/injection2pspatialparams.hh b/tutorial/ex2/injection2pspatialparams.hh
new file mode 100644
index 0000000000..fc85aa4c21
--- /dev/null
+++ b/tutorial/ex2/injection2pspatialparams.hh
@@ -0,0 +1,253 @@
+// -*- 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/implicit.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
+
+#include <dumux/porousmediumflow/2p2c/implicit/properties.hh>
+
+namespace Dumux
+{
+
+//forward declaration
+template<class TypeTag>
+class InjectionSpatialParams;
+
+namespace Properties
+{
+// The spatial parameters TypeTag
+NEW_TYPE_TAG(InjectionSpatialParams);
+
+// Set the spatial parameters
+SET_TYPE_PROP(InjectionSpatialParams, SpatialParams, InjectionSpatialParams<TypeTag>);
+
+// Set the material law parameterized by absolute saturations
+SET_TYPE_PROP(InjectionSpatialParams,
+ MaterialLaw,
+ EffToAbsLaw<RegularizedBrooksCorey<typename GET_PROP_TYPE(TypeTag, Scalar)> >);
+}
+
+/*!
+ * \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 TypeTag>
+class InjectionSpatialParams : public ImplicitSpatialParams<TypeTag>
+{
+ typedef ImplicitSpatialParams<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename Grid::ctype CoordScalar;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+
+ enum {
+ dim=GridView::dimension,
+ dimWorld=GridView::dimensionworld
+ };
+
+ typedef Dune::FieldVector<CoordScalar,dimWorld> GlobalPosition;
+
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename GridView::template Codim<0>::Entity Element;
+
+public:
+ typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
+ typedef typename MaterialLaw::Params MaterialLawParams;
+
+ /*!
+ * \brief The constructor
+ *
+ * \param gridView The grid view
+ */
+ InjectionSpatialParams(const GridView &gridView)
+ : ParentType(gridView)
+ {
+ layerBottom_ = 25.0;
+
+ // intrinsic permeabilities
+ fineK_ = 1e-13;
+ coarseK_ = 1e-12;
+
+ // porosities
+ finePorosity_ = 0.2;
+ coarsePorosity_ = 0.4;
+
+ //materialLawParams
+ fineEntryPressure_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.EntryPressureFine);
+ coarseEntryPressure_ = GET_RUNTIME_PARAM(TypeTag, Scalar, SpatialParams.EntryPressureCoarse);
+
+ // heat conductivity of granite
+ lambdaSolid_ = 2.8;
+
+ // residual saturations
+ fineMaterialParams_.setSwr(0.2);
+ fineMaterialParams_.setSnr(0.0);
+ coarseMaterialParams_.setSwr(0.2);
+ coarseMaterialParams_.setSnr(0.0);
+
+ // parameters for the Brooks-Corey law
+ fineMaterialParams_.setPe(fineEntryPressure_);
+ coarseMaterialParams_.setPe(coarseEntryPressure_);
+ fineMaterialParams_.setLambda(2.0);
+ coarseMaterialParams_.setLambda(2.0);
+ }
+
+ /*!
+ * \brief Returns the intrinsic permeability tensor \f$[m^2]\f$
+ *
+ * \param element The finite element
+ * \param fvGeometry The finite volume geometry of the element
+ * \param scvIdx The local index of the sub-control volume
+ */
+ const Scalar intrinsicPermeability(const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const int scvIdx) const
+ {
+ const GlobalPosition &globalPos = fvGeometry.subContVol[scvIdx].global;
+ if (isFineMaterial_(globalPos))
+ return fineK_;
+ return coarseK_;
+ }
+
+ /*!
+ * \brief Returns the porosity \f$[-]\f$
+ *
+ * \param element The finite element
+ * \param fvGeometry The finite volume geometry of the element
+ * \param scvIdx The local index of the sub-control volume
+ */
+ Scalar porosity(const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const int scvIdx) const
+ {
+ const GlobalPosition &globalPos = fvGeometry.subContVol[scvIdx].global;
+ if (isFineMaterial_(globalPos))
+ return finePorosity_;
+ return coarsePorosity_;
+ }
+
+
+ /*!
+ * \brief Returns the parameter object for the capillary-pressure/
+ * saturation material law
+ *
+ * \param element The finite element
+ * \param fvGeometry The finite volume geometry of the element
+ * \param scvIdx The local index of the sub-control volume
+ */
+ const MaterialLawParams& materialLawParams(const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const int scvIdx) const
+ {
+ const GlobalPosition &globalPos = fvGeometry.subContVol[scvIdx].global;
+ if (isFineMaterial_(globalPos))
+ return fineMaterialParams_;
+ return coarseMaterialParams_;
+ }
+
+ /*!
+ * These parameters are only needed for nonisothermal models. Comment them in if you want to implement the 2pni model.
+ */
+
+ /*!
+ * \brief Returns the heat capacity \f$[J / (kg K)]\f$ of the rock matrix.
+ *
+ * This is only required for non-isothermal models.
+ *
+ * \param element The finite element
+ * \param fvGeometry The finite volume geometry
+ * \param scvIdx The local index of the sub-control volume
+ */
+// Scalar solidHeatCapacity(const Element &element,
+// const FVElementGeometry &fvGeometry,
+// const int scvIdx) const
+// {
+// return 790; // specific heat capacity of granite [J / (kg K)]
+// }
+
+ /*!
+ * \brief Returns the mass density \f$[kg / m^3]\f$ of the rock matrix.
+ *
+ * This is only required for non-isothermal models.
+ *
+ * \param element The finite element
+ * \param fvGeometry The finite volume geometry
+ * \param scvIdx The local index of the sub-control volume
+ */
+// Scalar solidDensity(const Element &element,
+// const FVElementGeometry &fvGeometry,
+// const int scvIdx) const
+// {
+// return 2700; // density of granite [kg/m^3]
+// }
+
+ /*!
+ * \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the solid
+ *
+ * This is only required for non-isothermal models.
+ *
+ * \param element The finite element
+ * \param fvGeometry The finite volume geometry of the element
+ * \param scvIdx The local index of the sub-control volume
+ */
+// Scalar solidThermalConductivity(const Element &element,
+// const FVElementGeometry &fvGeometry,
+// const int scvIdx) const
+// {
+// return lambdaSolid_;
+// }
+
+private:
+ bool isFineMaterial_(const GlobalPosition &globalPos) const
+ { return globalPos[dimWorld-1] > layerBottom_; }
+
+ Scalar fineK_;
+ Scalar coarseK_;
+ Scalar layerBottom_;
+
+ Scalar finePorosity_;
+ Scalar coarsePorosity_;
+
+ Scalar lambdaSolid_;
+
+ Scalar fineEntryPressure_;
+ Scalar coarseEntryPressure_;
+
+ MaterialLawParams fineMaterialParams_;
+ MaterialLawParams coarseMaterialParams_;
+};
+
+}
+
+#endif
--
GitLab