From 65c6463920ae5aa50fb462e38bc47eb4ebd1a34e Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Fri, 27 Jan 2017 18:28:34 +0100
Subject: [PATCH] [richardsnc] Implement a Richards n-component model
The components are only in the water phase. The number of components
is statically extracted from the fluidsystem. Simply changing the
fluid system and boundary and initial conditions is sufficient to setup
a new model with n tracers.
---
dumux/porousmediumflow/CMakeLists.txt | 1 +
.../richardsnc/CMakeLists.txt | 1 +
.../richardsnc/implicit/CMakeLists.txt | 10 +
.../richardsnc/implicit/indices.hh | 63 +++
.../richardsnc/implicit/model.hh | 156 +++++++
.../richardsnc/implicit/properties.hh | 79 ++++
.../richardsnc/implicit/propertydefaults.hh | 194 +++++++++
.../richardsnc/implicit/volumevariables.hh | 392 ++++++++++++++++++
test/porousmediumflow/CMakeLists.txt | 1 +
.../richardsnc/CMakeLists.txt | 1 +
.../richardsnc/implicit/CMakeLists.txt | 24 ++
11 files changed, 922 insertions(+)
create mode 100644 dumux/porousmediumflow/richardsnc/CMakeLists.txt
create mode 100644 dumux/porousmediumflow/richardsnc/implicit/CMakeLists.txt
create mode 100644 dumux/porousmediumflow/richardsnc/implicit/indices.hh
create mode 100644 dumux/porousmediumflow/richardsnc/implicit/model.hh
create mode 100644 dumux/porousmediumflow/richardsnc/implicit/properties.hh
create mode 100644 dumux/porousmediumflow/richardsnc/implicit/propertydefaults.hh
create mode 100644 dumux/porousmediumflow/richardsnc/implicit/volumevariables.hh
create mode 100644 test/porousmediumflow/richardsnc/CMakeLists.txt
create mode 100644 test/porousmediumflow/richardsnc/implicit/CMakeLists.txt
diff --git a/dumux/porousmediumflow/CMakeLists.txt b/dumux/porousmediumflow/CMakeLists.txt
index 76eeae0efb..ff813a915a 100644
--- a/dumux/porousmediumflow/CMakeLists.txt
+++ b/dumux/porousmediumflow/CMakeLists.txt
@@ -16,5 +16,6 @@ add_subdirectory("implicit")
add_subdirectory("mpnc")
add_subdirectory("nonisothermal")
add_subdirectory("richards")
+add_subdirectory("richardsnc")
add_subdirectory("sequential")
add_subdirectory("tracer")
diff --git a/dumux/porousmediumflow/richardsnc/CMakeLists.txt b/dumux/porousmediumflow/richardsnc/CMakeLists.txt
new file mode 100644
index 0000000000..ba8341c614
--- /dev/null
+++ b/dumux/porousmediumflow/richardsnc/CMakeLists.txt
@@ -0,0 +1 @@
+add_subdirectory("implicit")
diff --git a/dumux/porousmediumflow/richardsnc/implicit/CMakeLists.txt b/dumux/porousmediumflow/richardsnc/implicit/CMakeLists.txt
new file mode 100644
index 0000000000..c19a6e9a86
--- /dev/null
+++ b/dumux/porousmediumflow/richardsnc/implicit/CMakeLists.txt
@@ -0,0 +1,10 @@
+
+#install headers
+install(FILES
+indices.hh
+model.hh
+localresidual.hh
+properties.hh
+propertydefaults.hh
+volumevariables.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/porousmediumflow/richardsnc/implicit)
diff --git a/dumux/porousmediumflow/richardsnc/implicit/indices.hh b/dumux/porousmediumflow/richardsnc/implicit/indices.hh
new file mode 100644
index 0000000000..f3912564cc
--- /dev/null
+++ b/dumux/porousmediumflow/richardsnc/implicit/indices.hh
@@ -0,0 +1,63 @@
+// -*- 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 Defines the primary variable and equation indices used by
+ * the richardsnc model
+ */
+
+#ifndef DUMUX_RICHARDSNC_INDICES_HH
+#define DUMUX_RICHARDSNC_INDICES_HH
+
+#include "properties.hh"
+
+namespace Dumux
+{
+// \{
+
+/*!
+ * \ingroup RichardsNCModel
+ * \ingroup ImplicitIndices
+ * \brief The indices for the isothermal Richards, n-component model.
+ */
+template <class TypeTag, int PVOffset = 0>
+struct RichardsNCIndices
+{
+
+ //! Set the index of the phases for accessing the volvars
+ static const int wPhaseIdx = 0;
+ static const int nPhaseIdx = 1;
+
+ //! Component indices
+ static const int compMainIdx = PVOffset + 0; //!< main component index
+
+ //! primary variable indices
+ static const int pressureIdx = PVOffset + 0; //!< pressure
+
+ //! \note These indices make sense if the first balance is replaced by the
+ //! total mass balance.
+
+ //! Equation indices
+ static const int conti0EqIdx = PVOffset + 0; //!< continuity equation index
+};
+
+// \}
+}
+
+#endif
diff --git a/dumux/porousmediumflow/richardsnc/implicit/model.hh b/dumux/porousmediumflow/richardsnc/implicit/model.hh
new file mode 100644
index 0000000000..51e576823c
--- /dev/null
+++ b/dumux/porousmediumflow/richardsnc/implicit/model.hh
@@ -0,0 +1,156 @@
+// -*- 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 Base class for all models which use the Richards,
+ * n-component fully implicit model.
+ */
+
+#ifndef DUMUX_RICHARDSNC_MODEL_HH
+#define DUMUX_RICHARDSNC_MODEL_HH
+
+#include "properties.hh"
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup RichardsModel
+ *
+ * \brief This model which implements a variant of the Richards'
+ * equation for quasi-twophase flow.
+ *
+ * In the unsaturated zone, Richards' equation
+ \f[
+ \frac{\partial\;\phi S_w \varrho_w}{\partial t}
+ -
+ \text{div} \left\lbrace
+ \varrho_w \frac{k_{rw}}{\mu_w} \; \mathbf{K} \;
+ \left( \text{\textbf{grad}}
+ p_w - \varrho_w \textbf{g}
+ \right)
+ \right\rbrace
+ =
+ q_w,
+ \f]
+ * is frequently used to
+ * approximate the water distribution above the groundwater level.
+ *
+ * It can be derived from the two-phase equations, i.e.
+ \f[
+ \phi\frac{\partial S_\alpha \varrho_\alpha}{\partial t}
+ -
+ \text{div} \left\lbrace
+ \varrho_\alpha \frac{k_{r\alpha}}{\mu_\alpha}\; \mathbf{K} \;
+ \left( \text{\textbf{grad}}
+ p_\alpha - \varrho_\alpha \textbf{g}
+ \right)
+ \right\rbrace
+ =
+ q_\alpha,
+ \f]
+ * where \f$\alpha \in \{w, n\}\f$ is the fluid phase,
+ * \f$\kappa \in \{ w, a \}\f$ are the components,
+ * \f$\rho_\alpha\f$ is the fluid density, \f$S_\alpha\f$ is the fluid
+ * saturation, \f$\phi\f$ is the porosity of the soil,
+ * \f$k_{r\alpha}\f$ is the relative permeability for the fluid,
+ * \f$\mu_\alpha\f$ is the fluid's dynamic viscosity, \f$\mathbf{K}\f$ is the
+ * intrinsic permeability, \f$p_\alpha\f$ is the fluid pressure and
+ * \f$g\f$ is the potential of the gravity field.
+ *
+ * In contrast to the full two-phase model, the Richards model assumes
+ * gas as the non-wetting fluid and that it exhibits a much lower
+ * viscosity than the (liquid) wetting phase. (For example at
+ * atmospheric pressure and at room temperature, the viscosity of air
+ * is only about \f$1\%\f$ of the viscosity of liquid water.) As a
+ * consequence, the \f$\frac{k_{r\alpha}}{\mu_\alpha}\f$ term
+ * typically is much larger for the gas phase than for the wetting
+ * phase. For this reason, the Richards model assumes that
+ * \f$\frac{k_{rn}}{\mu_n}\f$ is infinitly large. This implies that
+ * the pressure of the gas phase is equivalent to the static pressure
+ * distribution and that therefore, mass conservation only needs to be
+ * considered for the wetting phase.
+ *
+ * The model thus choses the absolute pressure of the wetting phase
+ * \f$p_w\f$ as its only primary variable. The wetting phase
+ * saturation is calculated using the inverse of the capillary
+ * pressure, i.e.
+ \f[
+ S_w = p_c^{-1}(p_n - p_w)
+ \f]
+ * holds, where \f$p_n\f$ is a given reference pressure. Nota bene,
+ * that the last step is assumes that the capillary
+ * pressure-saturation curve can be uniquely inverted, so it is not
+ * possible to set the capillary pressure to zero when using the
+ * Richards model!
+ */
+
+template<class TypeTag >
+class RichardsNCModel : public GET_PROP_TYPE(TypeTag, BaseModel)
+{
+ using ParentType = typename GET_PROP_TYPE(TypeTag, BaseModel);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+
+ static const int wPhaseIdx = Indices::wPhaseIdx;
+ static const int nPhaseIdx = Indices::nPhaseIdx;
+ static const int numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
+
+public:
+
+ /*!
+ * \brief Apply the initial conditions to the model.
+ *
+ * \param problem The object representing the problem which needs to
+ * be simulated.
+ */
+ void init(Problem& problem)
+ {
+ ParentType::init(problem);
+
+ // register standardized vtk output fields
+ auto& vtkOutputModule = problem.vtkOutputModule();
+ vtkOutputModule.addSecondaryVariable("Sw", [](const VolumeVariables& v){ return v.saturation(wPhaseIdx); });
+ vtkOutputModule.addSecondaryVariable("Sn", [](const VolumeVariables& v){ return v.saturation(nPhaseIdx); });
+ vtkOutputModule.addSecondaryVariable("pw", [](const VolumeVariables& v){ return v.pressure(wPhaseIdx); });
+ vtkOutputModule.addSecondaryVariable("pn", [](const VolumeVariables& v){ return v.pressure(nPhaseIdx); });
+ vtkOutputModule.addSecondaryVariable("pc", [](const VolumeVariables& v){ return v.capillaryPressure(); });
+ vtkOutputModule.addSecondaryVariable("density", [](const VolumeVariables& v){ return v.density(wPhaseIdx); });
+ vtkOutputModule.addSecondaryVariable("mobility", [](const VolumeVariables& v){ return v.mobility(wPhaseIdx); });
+ vtkOutputModule.addSecondaryVariable("kr", [](const VolumeVariables& v){ return v.relativePermeability(wPhaseIdx); });
+ vtkOutputModule.addSecondaryVariable("temperature", [](const VolumeVariables& v){ return v.temperature(); });
+ vtkOutputModule.addSecondaryVariable("porosity", [](const VolumeVariables& v){ return v.porosity(); });
+ if(GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity));
+ vtkOutputModule.addSecondaryVariable("pressure head", [](const VolumeVariables& v){ return v.pressureHead(wPhaseIdx); });
+ vtkOutputModule.addSecondaryVariable("water content", [](const VolumeVariables& v){ return v.waterContent(wPhaseIdx); });
+
+ for (int k = 0; k < numComponents; ++k)
+ vtkOutputModule.addSecondaryVariable("x^" + FluidSystem::phaseName(wPhaseIdx) + "_" + FluidSystem::componentName(k),
+ [k](const VolumeVariables& v){ return v.moleFraction(wPhaseIdx, k); });
+ }
+};
+
+} // end namespace Dumux
+
+#include "propertydefaults.hh"
+
+#endif
diff --git a/dumux/porousmediumflow/richardsnc/implicit/properties.hh b/dumux/porousmediumflow/richardsnc/implicit/properties.hh
new file mode 100644
index 0000000000..386bb184a6
--- /dev/null
+++ b/dumux/porousmediumflow/richardsnc/implicit/properties.hh
@@ -0,0 +1,79 @@
+// -*- 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/>. *
+ *****************************************************************************/
+/*!
+ * \ingroup Properties
+ * \ingroup ImplicitProperties
+ * \ingroup RichardsNCModel
+ * \file
+ *
+ * \brief Defines the properties required for the Richards,
+ * n-component fully implicit model.
+ */
+
+#ifndef DUMUX_RICHARDSNC_PROPERTIES_HH
+#define DUMUX_RICHARDSNC_PROPERTIES_HH
+
+
+#include <dumux/implicit/box/properties.hh>
+#include <dumux/implicit/cellcentered/tpfa/properties.hh>
+#include <dumux/porousmediumflow/nonisothermal/implicit/properties.hh>
+
+namespace Dumux
+{
+// \{
+namespace Properties
+{
+
+//////////////////////////////////////////////////////////////////
+// Type tags
+//////////////////////////////////////////////////////////////////
+
+//! The type tags for the implicit isothermal one-phase two-component problems
+NEW_TYPE_TAG(RichardsNC);
+NEW_TYPE_TAG(BoxRichardsNC, INHERITS_FROM(BoxModel, RichardsNC));
+NEW_TYPE_TAG(CCRichardsNC, INHERITS_FROM(CCTpfaModel, RichardsNC));
+
+//! The type tags for the corresponding non-isothermal problems
+NEW_TYPE_TAG(RichardsNCNI, INHERITS_FROM(RichardsNC, NonIsothermal));
+NEW_TYPE_TAG(BoxRichardsNCNI, INHERITS_FROM(BoxModel, RichardsNCNI));
+NEW_TYPE_TAG(CCRichardsNCNI, INHERITS_FROM(CCTpfaModel, RichardsNCNI));
+
+//////////////////////////////////////////////////////////////////
+// Property tags
+//////////////////////////////////////////////////////////////////
+
+NEW_PROP_TAG(NumPhases); //!< Number of fluid phases in the system
+NEW_PROP_TAG(NumComponents); //!< Number of fluid components in the system
+NEW_PROP_TAG(MaterialLaw); //!< The material law which ought to be used (by default extracted from the spatial parameters)
+NEW_PROP_TAG(MaterialLawParams); //!< The type of the parameter object for the material law (by default extracted from the spatial parameters)
+NEW_PROP_TAG(Indices); //!< Enumerations for the model
+NEW_PROP_TAG(SpatialParams); //!< The type of the spatial parameters
+NEW_PROP_TAG(EffectiveDiffusivityModel); //!< The employed model for the computation of the effective diffusivity
+NEW_PROP_TAG(FluidSystem); //!< Type of the multi-component relations
+NEW_PROP_TAG(FluidState); //!< Type of the fluid state to be used
+NEW_PROP_TAG(ProblemEnableGravity); //!< Returns whether gravity is considered in the problem
+NEW_PROP_TAG(UseMoles); //!< Defines whether mole (true) or mass (false) fractions are used
+NEW_PROP_TAG(SpatialParamsForchCoeff); //!< Property for the forchheimer coefficient
+NEW_PROP_TAG(TauTortuosity); //!< Tortuosity value (tau) used in macroscopic diffusion
+
+}
+// \}
+}
+
+#endif
diff --git a/dumux/porousmediumflow/richardsnc/implicit/propertydefaults.hh b/dumux/porousmediumflow/richardsnc/implicit/propertydefaults.hh
new file mode 100644
index 0000000000..e5f51c1361
--- /dev/null
+++ b/dumux/porousmediumflow/richardsnc/implicit/propertydefaults.hh
@@ -0,0 +1,194 @@
+// -*- 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/>. *
+ *****************************************************************************/
+/*!
+ * \ingroup Properties
+ * \ingroup ImplicitProperties
+ * \ingroup RichardsNCModel
+ * \file
+ *
+ * \brief Defines some default values for the properties of the
+ * Richards, n-component fully implicit model.
+ */
+
+#ifndef DUMUX_RICHARDSNC_PROPERTY_DEFAULTS_HH
+#define DUMUX_RICHARDSNC_PROPERTY_DEFAULTS_HH
+
+#include "properties.hh"
+#include "model.hh"
+#include "volumevariables.hh"
+#include "indices.hh"
+
+#include <dumux/porousmediumflow/compositional/localresidual.hh>
+#include <dumux/porousmediumflow/richards/implicit/newtoncontroller.hh>
+#include <dumux/porousmediumflow/nonisothermal/implicit/propertydefaults.hh>
+#include <dumux/material/spatialparams/implicit1p.hh>
+#include <dumux/material/fluidmatrixinteractions/diffusivitymillingtonquirk.hh>
+#include <dumux/material/fluidmatrixinteractions/1p/thermalconductivityaverage.hh>
+#include <dumux/material/components/simpleh2o.hh>
+#include <dumux/material/components/constant.hh>
+#include <dumux/material/fluidsystems/liquidphase2c.hh>
+#include <dumux/material/fluidstates/compositional.hh>
+
+namespace Dumux
+{
+// \{
+namespace Properties
+{
+//////////////////////////////////////////////////////////////////
+// Property values
+//////////////////////////////////////////////////////////////////
+SET_PROP(RichardsNC, NumEq)
+{
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ static const int value = FluidSystem::numComponents;
+};
+
+SET_PROP(RichardsNC, NumPhases)
+{
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ static const int value = FluidSystem::numPhases;
+ static_assert(value == 1, "You can only use one-phasic fluid systems with the Richards n-component model!");
+};
+
+SET_PROP(RichardsNC, NumComponents)
+{
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ static const int value = FluidSystem::numComponents;
+};
+
+SET_BOOL_PROP(RichardsNC, UseMoles, true); //!< Define that per default mole fractions are used in the balance equations
+
+//! Use the dedicated local residual
+SET_TYPE_PROP(RichardsNC, LocalResidual, CompositionalLocalResidual<TypeTag>);
+
+//! We set the replaceCompIdx to 0, i.e. the first equation is substituted with
+//! the total mass balance, i.e. the phase balance
+SET_INT_PROP(RichardsNC, ReplaceCompEqIdx, 0);
+
+//! define the model
+SET_TYPE_PROP(RichardsNC, Model, RichardsNCModel<TypeTag>);
+
+//! define the VolumeVariables
+SET_TYPE_PROP(RichardsNC, VolumeVariables, RichardsNCVolumeVariables<TypeTag>);
+
+//! Smarter newton controller
+SET_TYPE_PROP(RichardsNC, NewtonController, RichardsNewtonController<TypeTag>);
+
+/*!
+ *\brief The fluid system used by the model.
+ *
+ * By default this uses the liquid phase fluid system with simple H2O.
+ */
+SET_PROP(RichardsNC, FluidSystem)
+{
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using type = FluidSystems::LiquidPhaseTwoC<TypeTag, SimpleH2O<Scalar>, Constant<TypeTag, Scalar>>;
+};
+
+/*!
+ * \brief The fluid state which is used by the volume variables to
+ * store the thermodynamic state. This should be chosen
+ * appropriately for the model ((non-)isothermal, equilibrium, ...).
+ * This can be done in the problem.
+ */
+SET_PROP(RichardsNC, FluidState)
+{
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using type = CompositionalFluidState<Scalar, FluidSystem>;
+};
+
+/*!
+ * \brief Set type of the parameter objects for the material law
+ *
+ * By default this is just retrieved from the material law.
+ */
+SET_TYPE_PROP(RichardsNC, MaterialLawParams, typename GET_PROP_TYPE(TypeTag, MaterialLaw)::Params);
+
+//! Set the indices used
+SET_TYPE_PROP(RichardsNC, Indices, RichardsNCIndices<TypeTag>);
+//! The spatial parameters to be employed.
+//! Use ImplicitSpatialParamsOneP by default.
+SET_TYPE_PROP(RichardsNC, SpatialParams, ImplicitSpatialParamsOneP<TypeTag>);
+
+//! The model after Millington (1961) is used for the effective diffusivity
+SET_PROP(RichardsNC, EffectiveDiffusivityModel)
+{
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using type = DiffusivityMillingtonQuirk<Scalar>;
+};
+
+// enable gravity by default
+SET_BOOL_PROP(RichardsNC, ProblemEnableGravity, true);
+
+//physical processes to be considered by the isothermal model
+SET_BOOL_PROP(RichardsNC, EnableAdvection, true);
+SET_BOOL_PROP(RichardsNC, EnableMolecularDiffusion, true);
+SET_BOOL_PROP(RichardsNC, EnableEnergyBalance, false);
+
+//! default value for the forchheimer coefficient
+// Source: Ward, J.C. 1964 Turbulent flow in porous media. ASCE J. Hydraul. Div 90.
+// Actually the Forchheimer coefficient is also a function of the dimensions of the
+// porous medium. Taking it as a constant is only a first approximation
+// (Nield, Bejan, Convection in porous media, 2006, p. 10)
+SET_SCALAR_PROP(RichardsNC, SpatialParamsForchCoeff, 0.55);
+
+/*!
+ * \brief default value for tortuosity value (tau) used in macroscopic diffusion
+ *
+ * Value is 0.5 according to Carman 1937: <i>Fluid flow through granular beds</i>
+ * \cite carman1937
+ */
+SET_SCALAR_PROP(RichardsNC, TauTortuosity, 0.5);
+
+//! average is used as default model to compute the effective thermal heat conductivity
+SET_PROP(RichardsNCNI, ThermalConductivityModel)
+{
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using type = ThermalConductivityAverage<Scalar>;
+};
+
+//////////////////////////////////////////////////////////////////
+// Property values for isothermal model required for the general non-isothermal model
+//////////////////////////////////////////////////////////////////
+
+// set isothermal Model
+SET_TYPE_PROP(RichardsNCNI, IsothermalModel, RichardsNCModel<TypeTag>);
+
+//set isothermal VolumeVariables
+SET_TYPE_PROP(RichardsNCNI, IsothermalVolumeVariables, RichardsNCVolumeVariables<TypeTag>);
+
+//set isothermal LocalResidual
+SET_TYPE_PROP(RichardsNCNI, IsothermalLocalResidual, CompositionalLocalResidual<TypeTag>);
+
+//set isothermal Indices
+SET_TYPE_PROP(RichardsNCNI, IsothermalIndices, RichardsNCIndices<TypeTag>);
+
+//set isothermal NumEq
+SET_PROP(RichardsNCNI, IsothermalNumEq)
+{
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ static const int value = FluidSystem::numComponents;
+};
+
+} // end namespace Properties
+// \}
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/porousmediumflow/richardsnc/implicit/volumevariables.hh b/dumux/porousmediumflow/richardsnc/implicit/volumevariables.hh
new file mode 100644
index 0000000000..e1f325215b
--- /dev/null
+++ b/dumux/porousmediumflow/richardsnc/implicit/volumevariables.hh
@@ -0,0 +1,392 @@
+// -*- 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 @copybrief Dumux::RichardsNCVolumeVariables
+ */
+#ifndef DUMUX_RICHARDSNC_VOLUME_VARIABLES_HH
+#define DUMUX_RICHARDSNC_VOLUME_VARIABLES_HH
+
+#include <dumux/discretization/volumevariables.hh>
+
+#include "properties.hh"
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup RichardsNCModel
+ * \ingroup ImplicitVolumeVariables
+ * \brief Contains the quantities which are constant within a
+ * finite volume in the Richards, n-component model.
+ */
+template <class TypeTag>
+class RichardsNCVolumeVariables : public ImplicitVolumeVariables<TypeTag>
+{
+ using ParentType = ImplicitVolumeVariables<TypeTag>;
+
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Implementation = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+ using SpatialParams = typename GET_PROP_TYPE(TypeTag, SpatialParams);
+ using PermeabilityType = typename SpatialParams::PermeabilityType;
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+
+ enum
+ {
+ wPhaseIdx = Indices::wPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx,
+ pressureIdx = Indices::pressureIdx
+ };
+
+ static const bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
+ static const int dim = GridView::dimension;
+ static const int dimWorld = GridView::dimensionworld;
+ static const int numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
+
+ using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
+ using GlobalPosition = Dune::FieldVector<Scalar,dimWorld>;
+ using Element = typename GridView::template Codim<0>::Entity;
+
+public:
+
+ using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
+
+ /*!
+ * \copydoc ImplicitVolumeVariables::update
+ */
+ void update(const ElementSolutionVector &elemSol,
+ const Problem &problem,
+ const Element &element,
+ const SubControlVolume &scv)
+ {
+ ParentType::update(elemSol, problem, element, scv);
+
+ //calculate all secondary variables from the primary variables and store results in fluidstate
+ completeFluidState(elemSol, problem, element, scv, fluidState_);
+
+ const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
+ relativePermeabilityWetting_ = MaterialLaw::krw(materialParams, fluidState_.saturation(wPhaseIdx));
+
+ //! precompute the minimum capillary pressure (entry pressure)
+ minPc_ = MaterialLaw::pc(materialParams, /*Sw=*/ 1.0);
+ pn_ = problem.nonWettingReferencePressure();
+ porosity_ = problem.spatialParams().porosity(element, scv, elemSol);
+ permeability_ = problem.spatialParams().permeability(element, scv, elemSol);
+
+ // Second instance of a parameter cache.
+ // Could be avoided if diffusion coefficients also
+ // became part of the fluid state.
+ typename FluidSystem::ParameterCache paramCache;
+ paramCache.updatePhase(fluidState_, wPhaseIdx);
+
+ const int compIIdx = wPhaseIdx;
+ for (unsigned int compJIdx = 0; compJIdx < numComponents; ++compJIdx)
+ if(compIIdx != compJIdx)
+ setDiffusionCoefficient_(compJIdx,
+ FluidSystem::binaryDiffusionCoefficient(fluidState_,
+ paramCache,
+ wPhaseIdx,
+ compIIdx,
+ compJIdx));
+ }
+
+ /*!
+ * \copydoc ImplicitModel::completeFluidState
+ */
+ static void completeFluidState(const ElementSolutionVector &elemSol,
+ const Problem& problem,
+ const Element& element,
+ const SubControlVolume &scv,
+ FluidState& fluidState)
+ {
+ Scalar t = ParentType::temperature(elemSol, problem, element, scv);
+ fluidState.setTemperature(t);
+
+ const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
+ const auto& priVars = ParentType::extractDofPriVars(elemSol, scv);
+
+ // set the wetting pressure
+ fluidState.setPressure(wPhaseIdx, priVars[pressureIdx]);
+
+ // compute the capillary pressure to compute the saturation
+ // the material law takes care of cases where pc gets negative
+ const Scalar pc = problem.nonWettingReferencePressure() - fluidState.pressure(wPhaseIdx);
+ const Scalar sw = MaterialLaw::sw(materialParams, pc);
+ fluidState.setSaturation(wPhaseIdx, sw);
+
+ // set the mole/mass fractions
+ if(useMoles)
+ {
+ Scalar sumSecondaryFractions = 0.0;
+ for (int compIdx = 1; compIdx < numComponents; ++compIdx)
+ {
+ fluidState.setMoleFraction(wPhaseIdx, compIdx, priVars[compIdx]);
+ sumSecondaryFractions += priVars[compIdx];
+ }
+ fluidState.setMoleFraction(wPhaseIdx, 0, 1.0 - sumSecondaryFractions);
+ }
+ else
+ {
+ for (int compIdx = 1; compIdx < numComponents; ++compIdx)
+ fluidState.setMassFraction(wPhaseIdx, compIdx, priVars[compIdx]);
+ }
+
+ // density and viscosity
+ typename FluidSystem::ParameterCache paramCache;
+ paramCache.updateAll(fluidState);
+ fluidState.setDensity(wPhaseIdx, FluidSystem::density(fluidState, paramCache, wPhaseIdx));
+ fluidState.setViscosity(wPhaseIdx, FluidSystem::viscosity(fluidState, paramCache, wPhaseIdx));
+
+ // compute and set the enthalpy
+ fluidState.setEnthalpy(wPhaseIdx, Implementation::enthalpy(fluidState, paramCache, wPhaseIdx));
+ }
+
+ /*!
+ * \brief Return the fluid configuration at the given primary
+ * variables
+ */
+ const FluidState &fluidState() const
+ { return fluidState_; }
+
+ /*!
+ * \brief Return the temperature
+ */
+ Scalar temperature() const
+ { return fluidState_.temperature(); }
+
+ /*!
+ * \brief Returns the average porosity [] within the control volume.
+ *
+ * The porosity is defined as the ratio of the pore space to the
+ * total volume, i.e. \f[ \Phi := \frac{V_{pore}}{V_{pore} + V_{rock}} \f]
+ */
+ Scalar porosity() const
+ { return porosity_; }
+
+ /*!
+ * \brief Returns the permeability within the control volume in \f$[m^2]\f$.
+ */
+ const PermeabilityType& permeability() const
+ { return permeability_; }
+
+ /*!
+ * \brief Returns the average absolute saturation [] of a given
+ * fluid phase within the finite volume.
+ *
+ * The saturation of a fluid phase is defined as the fraction of
+ * the pore volume filled by it, i.e.
+ * \f[ S_\alpha := \frac{V_\alpha}{V_{pore}} = \phi \frac{V_\alpha}{V} \f]
+ *
+ * \param phaseIdx The index of the fluid phase
+ */
+ Scalar saturation(const int phaseIdx) const
+ { return phaseIdx == wPhaseIdx ? fluidState_.saturation(wPhaseIdx) : 1.0-fluidState_.saturation(wPhaseIdx); }
+
+ /*!
+ * \brief Returns the average mass density \f$\mathrm{[kg/m^3]}\f$ of a given
+ * fluid phase within the control volume.
+ *
+ * \param phaseIdx The index of the fluid phase
+ */
+ Scalar density(const int phaseIdx) const
+ { return phaseIdx == wPhaseIdx ? fluidState_.density(phaseIdx) : 0.0; }
+
+ /*!
+ * \brief Returns the effective pressure \f$\mathrm{[Pa]}\f$ of a given phase within
+ * the control volume.
+ *
+ * For the non-wetting phase (i.e. the gas phase), we assume
+ * infinite mobility, which implies that the non-wetting phase
+ * pressure is equal to the finite volume's reference pressure
+ * defined by the problem.
+ *
+ * \param phaseIdx The index of the fluid phase
+ */
+ Scalar pressure(const int phaseIdx) const
+ { return phaseIdx == wPhaseIdx ? fluidState_.pressure(phaseIdx) : pn_; }
+
+ /*!
+ * \brief Returns the effective mobility \f$\mathrm{[1/(Pa*s)]}\f$ of a given phase within
+ * the control volume.
+ *
+ * The mobility of a fluid phase is defined as the relative
+ * permeability of the phase (given by the chosen material law)
+ * divided by the dynamic viscosity of the fluid, i.e.
+ * \f[ \lambda_\alpha := \frac{k_{r\alpha}}{\mu_\alpha} \f]
+ *
+ * \param phaseIdx The index of the fluid phase
+ */
+ Scalar mobility(const int phaseIdx) const
+ { return relativePermeability(phaseIdx)/fluidState_.viscosity(phaseIdx); }
+
+ /*!
+ * \brief Returns the dynamic viscosity \f$\mathrm{[Pa*s]}\f$ of a given phase within
+ * the control volume.
+ *
+ * \param phaseIdx The index of the fluid phase
+ * \note The non-wetting phase is infinitely mobile
+ */
+ Scalar viscosity(const int phaseIdx) const
+ { return phaseIdx == wPhaseIdx ? fluidState_.viscosity(wPhaseIdx) : 0.0; }
+
+ /*!
+ * \brief Returns relative permeability [-] of a given phase within
+ * the control volume.
+ *
+ * \param phaseIdx The index of the fluid phase
+ */
+ Scalar relativePermeability(const int phaseIdx) const
+ { return phaseIdx == wPhaseIdx ? relativePermeabilityWetting_ : 1.0; }
+
+ /*!
+ * \brief Returns the effective capillary pressure \f$\mathrm{[Pa]}\f$ within the
+ * control volume.
+ *
+ * The capillary pressure is defined as the difference in
+ * pressures of the non-wetting and the wetting phase, i.e.
+ * \f[ p_c = p_n - p_w \f]
+ *
+ * \note Capillary pressures are always larger than the entry pressure
+ * This regularization doesn't affect the residual in which pc is not needed.
+ */
+ Scalar capillaryPressure() const
+ {
+ using std::max;
+ return max(minPc_, pn_ - fluidState_.pressure(wPhaseIdx));
+ }
+
+ /*!
+ * \brief Returns the pressureHead \f$\mathrm{[cm]}\f$ of a given phase within
+ * the control volume.
+ *
+ * For the non-wetting phase (i.e. the gas phase), we assume
+ * infinite mobility, which implies that the non-wetting phase
+ * pressure is equal to the finite volume's reference pressure
+ * defined by the problem.
+ *
+ * \param phaseIdx The index of the fluid phase
+ * \note this function is here as a convenience to the user to not have to
+ * manually do a conversion. It is not correct if the density is not constant
+ * or the gravity different
+ */
+ Scalar pressureHead(const int phaseIdx) const
+ { return 100.0 *(pressure(phaseIdx) - pn_)/density(phaseIdx)/9.81; }
+
+ /*!
+ * \brief Returns the water content
+ * fluid phase within the finite volume.
+ *
+ * The water content is defined as the fraction of
+ * the saturation devided by the porosity
+
+ * \param phaseIdx The index of the fluid phase
+ * \note this function is here as a convenience to the user to not have to
+ * manually do a conversion.
+ */
+ Scalar waterContent(const int phaseIdx) const
+ { return saturation(phaseIdx) * porosity_; }
+
+ /*!
+ * \brief Return molar density \f$\mathrm{[mol/m^3]}\f$ the of the fluid phase.
+ *
+ * We always forward to the fluid state with the phaseIdx property (see class description).
+ */
+ Scalar molarDensity(const int phaseIdx) const
+ { return phaseIdx == wPhaseIdx ? fluidState_.molarDensity(phaseIdx) : 0.0; }
+
+ /*!
+ * \brief Return mole fraction \f$\mathrm{[mol/mol]}\f$ of a component in the phase.
+ *
+ * \param compIdx The index of the component
+ *
+ * We always forward to the fluid state with the phaseIdx property (see class description).
+ */
+ Scalar moleFraction(const int phaseIdx, const int compIdx) const
+ { return phaseIdx == wPhaseIdx ? fluidState_.moleFraction(phaseIdx, compIdx) : 0.0; }
+
+ /*!
+ * \brief Return mass fraction \f$\mathrm{[kg/kg]}\f$ of a component in the phase.
+ *
+ * \param compIdx The index of the component
+ *
+ * We always forward to the fluid state with the phaseIdx property (see class description).
+ */
+ Scalar massFraction(const int phaseIdx, const int compIdx) const
+ { return phaseIdx == wPhaseIdx ? fluidState_.massFraction(phaseIdx, compIdx) : 0.0; }
+
+ /*!
+ * \brief Return concentration \f$\mathrm{[mol/m^3]}\f$ of a component in the phase.
+ *
+ * \param compIdx The index of the component
+ *
+ * We always forward to the fluid state with the phaseIdx property (see class description).
+ */
+ Scalar molarity(const int phaseIdx, const int compIdx) const
+ { return phaseIdx == wPhaseIdx ? fluidState_.molarity(phaseIdx, compIdx) : 0.0; }
+
+ /*!
+ * \brief Return the binary diffusion coefficient \f$\mathrm{[m^2/s]}\f$ in the fluid.
+ */
+ Scalar diffusionCoefficient(const int phaseIdx, const int compIdx) const
+ {
+ assert(phaseIdx == wPhaseIdx);
+ assert(compIdx > wPhaseIdx);
+ return diffCoefficient_[compIdx-1];
+ }
+
+ /*!
+ * \brief Returns the dispersivity of the fluid's streamlines.
+ */
+ const GlobalPosition &dispersivity() const
+ { return dispersivity_; }
+
+protected:
+ FluidState fluidState_; //! the fluid state
+ Scalar relativePermeabilityWetting_; //! the relative permeability of the wetting phase
+ Scalar porosity_; //! the porosity
+ PermeabilityType permeability_; //! the instrinsic permeability
+ Scalar pn_; //! the reference non-wetting pressure
+ Scalar minPc_; //! the minimum capillary pressure (entry pressure)
+
+private:
+ void setDiffusionCoefficient_(int compIdx, Scalar d)
+ {
+ assert(compIdx > wPhaseIdx);
+ diffCoefficient_[compIdx-1] = d;
+ }
+
+ std::array<Scalar, numComponents-1> diffCoefficient_;
+ GlobalPosition dispersivity_;
+
+ Implementation &asImp_()
+ { return *static_cast<Implementation*>(this); }
+
+ const Implementation &asImp_() const
+ { return *static_cast<const Implementation*>(this); }
+};
+
+}// end namespace Dumux
+
+#endif
diff --git a/test/porousmediumflow/CMakeLists.txt b/test/porousmediumflow/CMakeLists.txt
index 8b6532dc0c..83c6b1991c 100644
--- a/test/porousmediumflow/CMakeLists.txt
+++ b/test/porousmediumflow/CMakeLists.txt
@@ -13,4 +13,5 @@ add_subdirectory("3pwateroil")
add_subdirectory("co2")
add_subdirectory("mpnc")
add_subdirectory("richards")
+add_subdirectory("richardsnc")
add_subdirectory("tracer")
diff --git a/test/porousmediumflow/richardsnc/CMakeLists.txt b/test/porousmediumflow/richardsnc/CMakeLists.txt
new file mode 100644
index 0000000000..ba8341c614
--- /dev/null
+++ b/test/porousmediumflow/richardsnc/CMakeLists.txt
@@ -0,0 +1 @@
+add_subdirectory("implicit")
diff --git a/test/porousmediumflow/richardsnc/implicit/CMakeLists.txt b/test/porousmediumflow/richardsnc/implicit/CMakeLists.txt
new file mode 100644
index 0000000000..8352fc02e6
--- /dev/null
+++ b/test/porousmediumflow/richardsnc/implicit/CMakeLists.txt
@@ -0,0 +1,24 @@
+add_input_file_links()
+
+
+add_dumux_test(test_boxrichardsnc test_boxrichardsnc test_boxrichardsnc.cc
+ python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/richardslensbox-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/richardswelltracerbox-00008.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_boxrichardsnc")
+
+add_dumux_test(test_ccrichardsnc test_ccrichardsnc test_ccrichardsnc.cc
+ python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/richardslenscc-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/richardswelltracercc-00007.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_ccrichardsnc")
+
+#install sources
+install(FILES
+richardswelltracerproblem.hh
+richardswelltracerspatialparams.hh
+test_boxrichardsnc.cc
+test_ccrichardsnc.cc
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/test/porousmediumflow/richardsnc/implicit)
--
GitLab