Commit bfaf4a63 authored by Bernd Flemisch's avatar Bernd Flemisch
Browse files

update properties exercise

parent 7b7340da
# the compositional two-phase simulation program
# the property exercise simulation program
dune_add_test(NAME exercise_properties
SOURCES exercise_properties.cc
COMPILE_DEFINITIONS TYPETAG=Injection2p2pcCCProblem
CMD_ARGS exercise_properties.input)
# add tutorial to the common target
......
......@@ -2,13 +2,13 @@
<br>
## Problem set-up
The problem setup is identical to the previous [_Exercise Basic_](../exercise-basic/README.md).
The problem setup is identical to the two-phase incompressible test from DuMu<sup>x</sup>.
## Preparing the exercise
* Navigate to the directory `exercise-properties`
_Exercise Properties_ deals with a two-phase compositional problem (__2p2c__). Goal is to learn how to use compile and runtime parameters and the _DuMu<sup>x</sup> property system_.
_Exercise Properties_ deals with a two-phase immiscible incompressible problem (__2p__). The goal is to learn how to adapt compile-time parameters by employing the _DuMu<sup>x</sup> property system_.
<br><br>
### Task 1: Getting familiar with the code
......@@ -16,12 +16,11 @@ _Exercise Properties_ deals with a two-phase compositional problem (__2p2c__). G
Locate all the files you will need for this exercise
* The __main file__: `exercise_properties.cc`
* The __problem file__: `injection2p2cproblem.hh`
* The __spatial parameters file__: `injection2p2cspatialparams.hh`
* The __problem file__: `problem.hh`
* The __spatial parameters file__: `spatialparams.hh`
* The __input file__: `exercise_properties.input`
* Two header files containing:
* One header file containing:
* a custom __local residual__ in: `mylocalresidual.hh`
* a custom __material law__ in: `mymateriallaw.hh`
<hr><br><br>
......@@ -46,88 +45,32 @@ make exercise_properties
./exercise_properties
```
Note: Because the input file has the same name as the
executable, DuMuX will find it automatically.
If gnuplot is installed on your system, you should see a plot of the capillary pressure - saturation relationship.
<hr><br><br>
### Task 3: Implement and use a different material law
<hr>
DuMuX uses the term _material law_ to describe the law used to compute
* pc-Sw relations
* kr-Sw relations
* their inverse relations
The file `mymateriallaw.hh` contains a custom implementation of such a material law.
* Implement the method `Scalar pc(const Params &params, Scalar swe)` by implementing your own capillary pressure relationship, e.g. a simple linear relationship $`p_C(S_w) = 1\cdot 10^5 \cdot (1-S_w) + p_e`$.
Note: `MyMaterialLaw` uses the `BrooksCoreyParams` class as parameter input. You can get the entry pressure that is set in the spatial params as follows
```c++
const auto pe = params.pe();
```
The type (i.e. C++ type) of the material law is set in the file `injection2p2cspatialparams.hh` by declaring the following alias in the public section of the spatial parameters class:
```c++
using MaterialLaw = EffToAbsLaw<RegularizedBrooksCorey<Scalar>>;
```
* Make DuMuX use your own material law by including the header `mymateriallaw.hh` and changing the alias `MaterialLaw`. This will make sure that your material law is used everywhere else in the code.
Note: Also use the wrapper class `EffToAbsLaw`. It takes care of converting absolute to effective saturations considering residual saturations. `MyMaterialLaw`
as other material laws (like Brooks-Corey, VanGenuchten, ...) in DuMuX only deals with effective saturations.
* Verify your changes by recompiling and running the program. You should see a plot of your new function.
For the next task, disable the plotting feature by changing the settings in the input file `exercise2.input`
```ini
[Problem]
OnlyPlotMaterialLaws = false
```
executable, DuMu<sup>x</sup> will find it automatically.
<hr><br><br>
### Task 4: Implement a custom local residual
### Task 3: Implement a custom local residual
<hr>
Most types in DuMuX are properties that can be changed just like the material law. In the following task we implement our own 2p2c local residual, i.e. the class that computes the element residual in every Newton iteration. The file `mylocalresidual.hh` contains a copy of the similar to the original local residual class used for all compositional models renamed to `template<class TypeTag> class MyCompositionalLocalResidual`.
Many types in DuMu<sup>x</sup> are properties that can be changed on the problem level by using the property system. In the following task we implement our own 2p local residual, i.e. the class that computes the element residual in every Newton iteration. The file `mylocalresidual.hh` contains a copy of the original local residual class used for all immiscible models renamed to `template<class TypeTag> class MyLocalResidual`.
* Make DuMuX use this new local residual by inluding the header `mylocalresidual.hh` and setting the corresponding property in the `Property` namespace in the file `injection2p2cproblem.hh`
* Make DuMu<sup>x</sup> use this new local residual by including the header `mylocalresidual.hh` and setting the corresponding property in the `Property` namespace in the file `problem.hh`
```c++
// note that every property struct knows about TypeTag
SET_PROP(Injection2p2cTypeTag, LocalResidual)
SET_PROP(TwoPIncompressible, LocalResidual)
{
using type = MyCompositionalLocalResidual<TypeTag>;
using type = MyLocalResidual<TypeTag>;
};
// or using the convenience macro
SET_TYPE_PROP(Injection2p2cTypeTag, LocalResidual,
MyCompositionalLocalResidual<TypeTag>);
SET_TYPE_PROP(TwoPIncompressible, LocalResidual,
MyLocalResidual<TypeTag>);
```
You want to make the new local residual special by adding a switch enabling / disabling diffusion. We will achieve this with a DuMuX parameter which is read from the input file and defaults to a property value if the input file doesn't contain the parameter.
You want to make the simplify the original local residual special by using the assumption that incompressible fluid phases are present. As a consequence, one can balance phase volume instead of phase mass.
* Modify the `computeFlux` method to only call the `diffusiveFlux` method if diffusion is enabled. You can get the new parameter by adding the lines
```c++
// ... in the computeFlux method of MyCompositionalLocalResidual
auto enableDiffusion = getParam<bool>("Problem.EnableDiffusion", true);
```
You can now enable and disable diffusion through the input file
```ini
[Problem]
EnableDiffusion = true / false
```
* Eliminate the density from the `computeStorage` and the `computeFlux` methods.
* Verify the difference in the parameter $`x_w^{N2}`$, i.e. the mole fraction of nitrogen in the
water phase, with and without diffusion.
* Adapt the relevant routine in the problem file such that volume instead of mass is injected. The density of the employed NAPL is 1460 kg/m^3. For a first try, you can use the hardcoded value.
Note that due to diffusion being a slow process you
can only see the difference in later times.
* Generalize your approach by using the component's `liquidDensity(t, p)` function instead of the hardcoded value.
......@@ -19,7 +19,7 @@
/*!
* \file
*
* \brief Test for the two-phase two-component CC model.
* \brief test for the two-phase porousmedium flow model
*/
#include <config.h>
......@@ -30,14 +30,18 @@
#include <dune/common/timer.hh>
#include <dune/grid/io/file/dgfparser/dgfexception.hh>
#include <dune/grid/io/file/vtk.hh>
#include <dune/istl/io.hh>
#include "problem.hh"
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/valgrind.hh>
#include <dumux/common/dumuxmessage.hh>
#include <dumux/common/defaultusagemessage.hh>
#include <dumux/linear/amgbackend.hh>
#include <dumux/nonlinear/privarswitchnewtonsolver.hh>
#include <dumux/nonlinear/newtonsolver.hh>
#include <dumux/assembly/fvassembler.hh>
#include <dumux/assembly/diffmethod.hh>
......@@ -47,14 +51,44 @@
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/io/grid/gridmanager.hh>
#include "injection2p2cproblem.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 arguments for this program is:\n"
"\t-TimeManager.TEnd End of the simulation [s] \n"
"\t-TimeManager.DtInitial Initial timestep size [s] \n"
"\t-Grid.LowerLeft Lower left corner coordinates\n"
"\t-Grid.UpperRight Upper right corner coordinates\n"
"\t-Grid.Cells Number of cells in respective coordinate directions\n"
"\t definition in DGF format\n"
"\t-SpatialParams.LensLowerLeft coordinates of the lower left corner of the lens [m] \n"
"\t-SpatialParams.LensUpperRight coordinates of the upper right corner of the lens [m] \n"
"\t-SpatialParams.Permeability Permeability of the domain [m^2] \n"
"\t-SpatialParams.PermeabilityLens Permeability of the lens [m^2] \n";
std::cout << errorMessageOut
<< "\n";
}
}
int main(int argc, char** argv)try
int main(int argc, char** argv) try
{
using namespace Dumux;
// define the type tag for this problem
using TypeTag = TTAG(Injection2p2pcCCTypeTag);
using TypeTag = TTAG(TwoPIncompressibleTpfa);
// initialize MPI, finalize is done automatically on exit
const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
......@@ -64,7 +98,7 @@ int main(int argc, char** argv)try
DumuxMessage::print(/*firstCall=*/true);
// parse command line arguments and input file
Parameters::init(argc, argv);
Parameters::init(argc, argv, usage);
// try to create a grid (from the given grid file or the input file)
GridManager<typename GET_PROP_TYPE(TypeTag, Grid)> gridManager;
......@@ -103,18 +137,19 @@ int main(int argc, char** argv)try
const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
auto dt = getParam<Scalar>("TimeLoop.DtInitial");
// check if we are about to restart a previously interrupted simulation
Scalar restartTime = 0;
if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart"))
restartTime = getParam<Scalar>("TimeLoop.Restart");
// intialize the vtk output module
using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
VtkOutputModule<TypeTag> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
VtkOutputFields::init(vtkWriter); //! Add model specific output fields
// use non-conforming output for the test with interface solver
const auto ncOutput = getParam<bool>("Problem.UseNonConformingOutput", false);
VtkOutputModule<TypeTag> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name(), "",
(ncOutput ? Dune::VTK::nonconforming : Dune::VTK::conforming));
VtkOutputFields::init(vtkWriter); //!< Add model specific output fields
vtkWriter.write(0.0);
// instantiate time loop
auto timeLoop = std::make_shared<TimeLoop<Scalar>>(restartTime, dt, tEnd);
auto timeLoop = std::make_shared<TimeLoop<Scalar>>(0.0, dt, tEnd);
timeLoop->setMaxTimeStepSize(maxDt);
// the assembler with time loop for instationary problem
......@@ -126,8 +161,7 @@ int main(int argc, char** argv)try
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>;
using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
NewtonSolver nonLinearSolver(assembler, linearSolver);
// time loop
......@@ -136,9 +170,6 @@ int main(int argc, char** argv)try
// 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);
......@@ -149,16 +180,20 @@ int main(int argc, char** argv)try
// advance to the time loop to the next step
timeLoop->advanceTimeStep();
// write vtk output
vtkWriter.write(timeLoop->time());
// report statistics of this time step
timeLoop->reportTimeStep();
// set new dt as suggested by the newton solver
// set new dt as suggested by the Newton solver
timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
vtkWriter.write(timeLoop->time());
} while (!timeLoop->finished());
// output some Newton statistics
nonLinearSolver.report();
timeLoop->finalize(leafGridView.comm());
////////////////////////////////////////////////////////////
......
[TimeLoop]
DtInitial = 3600 # in seconds
TEnd = 3.154e9 # in seconds, i.e ten years
DtInitial = 250 # [s]
TEnd = 30000 # [s]
[Grid]
LowerLeft = 0 0
UpperRight = 60 40
Cells = 24 16
UpperRight = 6 4
Cells = 48 32
[Problem]
Name = infiltration
OnlyPlotMaterialLaws = true
AquiferDepth = 2700.0 # m
TotalAreaSpecificInflow = 1e-4 # kg / (s*m^2)
InjectionDuration = 2.628e6 # in seconds, i.e. one month
[SpatialParams]
LensLowerLeft = 1.0 2.0 # [m] coordinates of the lower left lens corner
LensUpperRight = 4.0 3.0 # [m] coordinates of the upper right lens corner
# TODO: dumux-course-task
# Set Problem.EnableGravity
# Set Problem.EnableDiffusion
[Problem]
Name = exercise_properties
[SpatialParams]
PermeabilityAquitard = 1e-15 # m^2
EntryPressureAquitard = 4.5e4 # Pa
PermeabilityAquifer = 1e-12 # m^2
EntryPressureAquifer = 1e4 # Pa
[Newton]
MaxRelativeShift = 1e-5
......@@ -18,127 +18,113 @@
*****************************************************************************/
/*!
* \file
*
* \brief Element-wise calculation of the local residual for problems
* using compositional fully implicit model.
* \ingroup PorousmediumImmiscible
* \brief Element-wise calculation of the residual for problems
* using the n-phase immiscible fully implicit models.
*/
#ifndef DUMUX_MY_COMPOSITIONAL_LOCAL_RESIDUAL_HH
#define DUMUX_MY_COMPOSITIONAL_LOCAL_RESIDUAL_HH
#ifndef DUMUX_MY_LOCAL_RESIDUAL_HH
#define DUMUX_MY_LOCAL_RESIDUAL_HH
#include <dumux/common/properties.hh>
namespace Dumux {
namespace Dumux
{
/*!
* \ingroup Implicit
* \ingroup ImplicitLocalResidual
* \brief Element-wise calculation of the local residual for problems
* using compositional fully implicit model.
*
* \ingroup PorousmediumImmiscible
* \brief Element-wise calculation of the residual for problems
* using the n-phase immiscible fully implicit models.
*/
template<class TypeTag>
class MyCompositionalLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
class MyLocalResidual : public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
{
using ParentType = typename GET_PROP_TYPE(TypeTag, BaseLocalResidual);
using Implementation = typename GET_PROP_TYPE(TypeTag, LocalResidual);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, GridFluxVariablesCache)::LocalView;
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
using SubControlVolume = typename FVElementGeometry::SubControlVolume;
using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, GridFluxVariablesCache)::LocalView;
using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Element = typename GridView::template Codim<0>::Entity;
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables)::LocalView;
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using EnergyLocalResidual = typename GET_PROP_TYPE(TypeTag, EnergyLocalResidual);
static constexpr int numPhases = GET_PROP_TYPE(TypeTag, ModelTraits)::numPhases();
static constexpr int numComponents = GET_PROP_TYPE(TypeTag, ModelTraits)::numComponents();
enum { conti0EqIdx = Indices::conti0EqIdx };
using ModelTraits = typename GET_PROP_TYPE(TypeTag, ModelTraits);
static constexpr int numPhases = ModelTraits::numPhases();
static constexpr int conti0EqIdx = ModelTraits::Indices::conti0EqIdx; //!< first index for the mass balance
public:
using ParentType::ParentType;
/*!
* \brief Evaluate the amount of all conservation quantities
* (e.g. phase mass) within a sub-control volume.
*
* The result should be averaged over the volume (e.g. phase mass
* inside a sub control volume divided by the volume)
*
* \param storage The mass of the component within the sub-control volume
* \param scvIdx The SCV (sub-control-volume) index
* \param usePrevSol Evaluate function with solution of current or previous time step
* \brief Evaluate the rate of change of all conservation
* quantites (e.g. phase mass) within a sub-control
* volume of a finite volume element for the immiscible models.
* \param problem The problem
* \param scv The sub control volume
* \param volVars The current or previous volVars
* \note This function should not include the source and sink terms.
* \note The volVars can be different to allow computing
* the implicit euler time derivative here
*/
ResidualVector computeStorage(const Problem& problem,
const SubControlVolume& scv,
const VolumeVariables& volVars) const
// TODO: eliminate density from the storage term
NumEqVector computeStorage(const Problem& problem,
const SubControlVolume& scv,
const VolumeVariables& volVars) const
{
ResidualVector storage(0.0);
// compute storage term of all components within all phases
// partial time derivative of the phase mass
NumEqVector storage;
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
{
auto eqIdx = conti0EqIdx + compIdx;
storage[eqIdx] += volVars.porosity()
* volVars.saturation(phaseIdx)
* volVars.molarDensity(phaseIdx)
* volVars.moleFraction(phaseIdx, compIdx);
}
auto eqIdx = conti0EqIdx + phaseIdx;
storage[eqIdx] = volVars.porosity()
* volVars.density(phaseIdx)
* volVars.saturation(phaseIdx);
//! The energy storage in the fluid phase with index phaseIdx
EnergyLocalResidual::fluidPhaseStorage(storage, scv, volVars, phaseIdx);
}
//! The energy storage in the solid matrix
EnergyLocalResidual::solidPhaseStorage(storage, scv, volVars);
return storage;
}
/*!
* \brief Evaluates the total flux of all conservation quantities
* over a face of a sub-control volume.
* \brief Evaluate the mass flux over a face of a sub control volume
*
* \param flux The flux over the SCV (sub-control-volume) face for each component
* \param fIdx The index of the SCV face
* \param onBoundary A boolean variable to specify whether the flux variables
* are calculated for interior SCV faces or boundary faces, default=false
* \param problem The problem
* \param element The element
* \param fvGeometry The finite volume geometry context
* \param elemVolVars The volume variables for all flux stencil elements
* \param scvf The sub control volume face to compute the flux on
* \param elemFluxVarsCache The cache related to flux compuation
*/
ResidualVector computeFlux(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf,
const ElementFluxVariablesCache& elemFluxVarsCache) const
// TODO: eliminate the density from the flux term
NumEqVector computeFlux(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace& scvf,
const ElementFluxVariablesCache& elemFluxVarsCache) const
{
FluxVariables fluxVars;
fluxVars.init(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
// get upwind weights into local scope
ResidualVector flux(0.0);
// TODO: dumux-course-task
// get parameter Problem.EnableDiffusion
// advective fluxes
NumEqVector flux;
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
const auto diffusiveFluxes = fluxVars.molecularDiffusionFlux(phaseIdx);
for (int compIdx = 0; compIdx < numComponents; ++compIdx)
{
// get equation index
const auto eqIdx = conti0EqIdx + compIdx;
// the physical quantities for which we perform upwinding
const auto upwindTerm = [phaseIdx,compIdx] (const auto& volVars)
{ return volVars.molarDensity(phaseIdx)*volVars.moleFraction(phaseIdx, compIdx)*volVars.mobility(phaseIdx); };
flux[eqIdx] += fluxVars.advectiveFlux(phaseIdx, upwindTerm);
// TODO: dumux-course-task
// same here: only add diffusive fluxes if diffusion is enabled
flux[eqIdx] += diffusiveFluxes[compIdx];
}
// the physical quantities for which we perform upwinding
auto upwindTerm = [phaseIdx](const auto& volVars)
{ return volVars.density(phaseIdx)*volVars.mobility(phaseIdx); };
auto eqIdx = conti0EqIdx + phaseIdx;
flux[eqIdx] = fluxVars.advectiveFlux(phaseIdx, upwindTerm);
//! Add advective phase energy fluxes. For isothermal model the contribution is zero.
EnergyLocalResidual::heatConvectionFlux(flux, fluxVars, phaseIdx);
......@@ -149,13 +135,6 @@ public:
return flux;
}
protected:
Implementation *asImp_()
{ return static_cast<Implementation *> (this); }
const Implementation *asImp_() const
{ return static_cast<const Implementation *> (this); }
};
} // end namespace Dumux
......
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
*
* \brief Implementation of the capillary pressure and
* relative permeability <-> saturation relations.
*
*/
#ifndef DUMUX_MY_MATERIAL_LAW_HH
#define DUMUX_MY_MATERIAL_LAW_HH
#include <dumux/material/fluidmatrixinteractions/2p/brookscoreyparams.hh>
#include <cmath>
#include <algorithm>
namespace Dumux
{
/*!
* \ingroup fluidmatrixinteractionslaws
* \note a simple material law using the BrooksCoreyParams
*/
template <class ScalarT, class ParamsT = BrooksCoreyParams<ScalarT> >
class MyMaterialLaw
{
public:
typedef ParamsT Params;
typedef typename Params::Scalar Scalar;
/*!
* \brief The capillary pressure-saturation curve
* \param swe saturation of the wetting phase \f$\mathrm{[\overline{S}_w]}\f$
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen,
and then the params container is constructed accordingly. Afterwards the values are set there, too.
* \return capillary pressure
* TODO: dumux-course-task
* Implement the pc(swe) function
*/
static Scalar pc(const Params &params, Scalar swe)
{
return 0.0;
}