diff --git a/dumux/common/deprecated.hh b/dumux/common/deprecated.hh
index c8e75e6c2ccbd721468cabed15603d5c3ca1eb92..7e3c82c40f3e1acfed9cd7b63a8264e75e9a255c 100644
--- a/dumux/common/deprecated.hh
+++ b/dumux/common/deprecated.hh
@@ -27,7 +27,10 @@
#include <utility>
+#include <dune/common/ftraits.hh>
+#include <dune/common/exceptions.hh>
#include <dune/common/std/type_traits.hh>
+#include <dumux/common/numeqvector.hh>
namespace Dumux {
@@ -63,8 +66,92 @@ constexpr inline bool hasSCVFGeometryInterface()
#pragma clang diagnostic pop
#endif // __clang__
+template<bool enableWaterDiffusionInAir>
+struct ExtendedRichardsHelper
+{
+ template<bool b>
+ [[deprecated("Enabling the extended Richards model through a template parameter/properties is deprecated and will be removed after release (3.6). Use the new model ExtendedRichards instead.")]]
+ static constexpr void extendedRichardsWithTemplateParameter() {}
+
+ static constexpr bool isExtendedRichards()
+ {
+ if constexpr(enableWaterDiffusionInAir)
+ extendedRichardsWithTemplateParameter<enableWaterDiffusionInAir>();
+ return enableWaterDiffusionInAir;
+ }
+};
+
+/*!
+ * \ingroup RichardsModel
+ * \brief A primary variable vector with a state to allow variable switches.
+ */
+template<class PrimaryVariables, class StateType>
+class RichardsSwitchablePrimaryVariables : public PrimaryVariables
+{
+ using ParentType = PrimaryVariables;
+public:
+ //! Inherit the constructors
+ using ParentType::ParentType;
+ //! Use the assignment operators from the field vector
+ using ParentType::operator=;
+
+ [[deprecated("Will be removed after release (3.6). Normal Richards does not require setting a state. Use ExtendedRichards for extended model. "
+ "After removing all .setState()/.state() calls in the user files the remaining deprecation warnings can be removed by changing the PrimaryVariables property to Dune::FieldVector. "
+ "template<class TypeTag> "
+ "struct PrimaryVariables<TypeTag, TTag::Richards> "
+ "{ using type = Dune::FieldVector<GetPropType<TypeTag, Properties::Scalar>, GetPropType<TypeTag, Properties::ModelTraits>::numEq()>; }")]]
+ StateType state() const
+ {
+ return state_;
+ }
+
+ //! Set the state of this primary variable object, e.g. the phase presence.
+ [[deprecated("Will be removed after release (3.6). Normal Richards does not require setting a state. Use ExtendedRichards for extended model."
+ "After removing all .setState()/.state() calls in the user files the remaining deprecation warnings can be removed by changing the PrimaryVariables property to Dune::FieldVector. "
+ "template<class TypeTag> "
+ "struct PrimaryVariables<TypeTag, TTag::Richards> "
+ "{ using type = Dune::FieldVector<GetPropType<TypeTag, Properties::Scalar>, GetPropType<TypeTag, Properties::ModelTraits>::numEq()>; }")]]
+ void setState(StateType state)
+ {
+ state_ = std::move(state);
+ stateIsSet_ = true;
+ }
+
+ //! Invalidate the state
+ void invalidateState()
+ {
+ stateIsSet_ = false;
+ }
+
+private:
+ StateType state_;
+ bool stateIsSet_{false};
+};
+
+
} // end namespace Deprecated
#endif
+/*!
+ * \ingroup PorousmediumflowModels
+ * \brief The corresponding NumEqVectorTraits for the primary variables with switchable state
+ */
+template<class PrimaryVariables, class StateType>
+struct NumEqVectorTraits<Deprecated::RichardsSwitchablePrimaryVariables<PrimaryVariables, StateType>>
+{
+ static constexpr std::size_t numEq = PrimaryVariables::size();
+ using type = PrimaryVariables;
+};
+
} // end namespace Dumux
+
+// specialize field traits for this type
+namespace Dune {
+
+template <class PrimaryVariables, class StateType>
+struct FieldTraits<Dumux::Deprecated::RichardsSwitchablePrimaryVariables<PrimaryVariables, StateType>>
+: public FieldTraits<PrimaryVariables>
+{};
+
+} // end namespace Dune
#endif
diff --git a/dumux/porousmediumflow/CMakeLists.txt b/dumux/porousmediumflow/CMakeLists.txt
index 4b82b2cdb7bae10abdb64499e4fba12c662f6f50..43fb49ffcf6e124d5c674e5f3c0643c8e1b68205 100644
--- a/dumux/porousmediumflow/CMakeLists.txt
+++ b/dumux/porousmediumflow/CMakeLists.txt
@@ -18,6 +18,7 @@ add_subdirectory(mpnc)
add_subdirectory(nonequilibrium)
add_subdirectory(nonisothermal)
add_subdirectory(richards)
+add_subdirectory(richardsextended)
add_subdirectory(richardsnc)
add_subdirectory(solidenergy)
add_subdirectory(tracer)
diff --git a/dumux/porousmediumflow/richards/model.hh b/dumux/porousmediumflow/richards/model.hh
index 6da72cc3838a3f376e5e942097308d047cf114f8..93b53a661ebb7d8b9d1610984c1116c41a0f8b79 100644
--- a/dumux/porousmediumflow/richards/model.hh
+++ b/dumux/porousmediumflow/richards/model.hh
@@ -90,8 +90,11 @@
#ifndef DUMUX_RICHARDS_MODEL_HH
#define DUMUX_RICHARDS_MODEL_HH
+#include <type_traits>
+
#include <dune/common/fvector.hh>
+#include <dumux/common/deprecated.hh>
#include <dumux/common/properties.hh>
#include <dumux/porousmediumflow/immiscible/localresidual.hh>
@@ -132,7 +135,7 @@ struct RichardsModelTraits
static constexpr int numFluidComponents() { return 1; }
static constexpr bool enableAdvection() { return true; }
- static constexpr bool enableMolecularDiffusion() { return enableDiff; }
+ static constexpr bool enableMolecularDiffusion() { return Dumux::Deprecated::ExtendedRichardsHelper<enableDiff>::isExtendedRichards(); }
static constexpr bool enableEnergyBalance() { return false; }
//! The Richards model has some assumptions on the fluid systems
@@ -268,7 +271,10 @@ private:
using PrimaryVariablesVector = Dune::FieldVector<GetPropType<TypeTag, Properties::Scalar>,
GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
public:
- using type = SwitchablePrimaryVariables<PrimaryVariablesVector, int>;
+ // TODO: After release 3.6, use the default (FieldVector) for Richards and remove the whole property specialization
+ using type = std::conditional_t<getPropValue<TypeTag, Properties::EnableWaterDiffusionInAir>(),
+ SwitchablePrimaryVariables<PrimaryVariablesVector, int>,
+ Deprecated::RichardsSwitchablePrimaryVariables<PrimaryVariablesVector, int>>;
};
/*!
diff --git a/dumux/porousmediumflow/richards/volumevariables.hh b/dumux/porousmediumflow/richards/volumevariables.hh
index 076b5cba53cc46f42bdda3113087be680e0bea3e..23a329164c7f81354ecc9c7578a05c352d3f492b 100644
--- a/dumux/porousmediumflow/richards/volumevariables.hh
+++ b/dumux/porousmediumflow/richards/volumevariables.hh
@@ -30,14 +30,14 @@
#include <dune/common/exceptions.hh>
#include <dumux/common/typetraits/state.hh>
+#include <dumux/common/deprecated.hh>
#include <dumux/porousmediumflow/volumevariables.hh>
#include <dumux/porousmediumflow/nonisothermal/volumevariables.hh>
#include <dumux/material/idealgas.hh>
-#include <dumux/material/constants.hh>
#include <dumux/material/solidstates/updatesolidvolumefractions.hh>
-#include "primaryvariableswitch.hh"
+#include <dumux/porousmediumflow/richardsextended/primaryvariableswitch.hh>
namespace Dumux {
@@ -92,7 +92,7 @@ public:
using SolidSystem = typename Traits::SolidSystem;
using Indices = typename Traits::ModelTraits::Indices;
//! If water diffusion in air is enabled
- static constexpr bool enableWaterDiffusionInAir() { return ModelTraits::enableMolecularDiffusion(); };
+ static constexpr bool enableWaterDiffusionInAir() { return Dumux::Deprecated::ExtendedRichardsHelper<ModelTraits::enableMolecularDiffusion()>::isExtendedRichards(); }
//! Export phase indices
static constexpr auto liquidPhaseIdx = Traits::FluidSystem::phase0Idx;
@@ -117,9 +117,6 @@ public:
const auto fluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element, scv, elemSol);
- const auto& priVars = elemSol[scv.localDofIndex()];
- const auto phasePresence = getState_(priVars);
-
// precompute the minimum capillary pressure (entry pressure)
// needed to make sure we don't compute unphysical capillary pressures and thus saturations
minPc_ = fluidMatrixInteraction.endPointPc();
@@ -127,102 +124,112 @@ public:
//update porosity before calculating the effective properties depending on it
updateSolidVolumeFractions(elemSol, problem, element, scv, solidState_, numFluidComps);
- typename FluidSystem::ParameterCache paramCache;
- auto getEffectiveDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx)
+ if constexpr (enableWaterDiffusionInAir())
{
- return EffDiffModel::effectiveDiffusionCoefficient(*this, phaseIdx, compIIdx, compJIdx);
- };
+ const auto& priVars = elemSol[scv.localDofIndex()];
+ const auto phasePresence = priVars.state();
- if (phasePresence == Indices::gasPhaseOnly)
- {
- moleFraction_[liquidPhaseIdx] = 1.0;
- massFraction_[liquidPhaseIdx] = 1.0;
- moleFraction_[gasPhaseIdx] = priVars[Indices::switchIdx];
+ typename FluidSystem::ParameterCache paramCache;
+ auto getEffectiveDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx)
+ {
+ return EffDiffModel::effectiveDiffusionCoefficient(*this, phaseIdx, compIIdx, compJIdx);
+ };
- const auto averageMolarMassGasPhase = (moleFraction_[gasPhaseIdx]*FluidSystem::molarMass(liquidPhaseIdx)) +
- ((1-moleFraction_[gasPhaseIdx])*FluidSystem::molarMass(gasPhaseIdx));
+ if (phasePresence == Indices::gasPhaseOnly)
+ {
+ moleFraction_[liquidPhaseIdx] = 1.0;
+ massFraction_[liquidPhaseIdx] = 1.0;
+ moleFraction_[gasPhaseIdx] = priVars[Indices::switchIdx];
- //X_w = x_w* M_w/ M_avg
- massFraction_[gasPhaseIdx] = priVars[Indices::switchIdx]*FluidSystem::molarMass(liquidPhaseIdx)/averageMolarMassGasPhase;
+ const auto averageMolarMassGasPhase = (moleFraction_[gasPhaseIdx]*FluidSystem::molarMass(liquidPhaseIdx)) +
+ ((1-moleFraction_[gasPhaseIdx])*FluidSystem::molarMass(gasPhaseIdx));
- fluidState_.setSaturation(liquidPhaseIdx, 0.0);
- fluidState_.setSaturation(gasPhaseIdx, 1.0);
+ //X_w = x_w* M_w/ M_avg
+ massFraction_[gasPhaseIdx] = priVars[Indices::switchIdx]*FluidSystem::molarMass(liquidPhaseIdx)/averageMolarMassGasPhase;
- EnergyVolVars::updateTemperature(elemSol, problem, element, scv, fluidState_, solidState_);
+ fluidState_.setSaturation(liquidPhaseIdx, 0.0);
+ fluidState_.setSaturation(gasPhaseIdx, 1.0);
- // get pc for sw = 0.0
- const Scalar pc = fluidMatrixInteraction.pc(0.0);
+ EnergyVolVars::updateTemperature(elemSol, problem, element, scv, fluidState_, solidState_);
- // set the wetting pressure
- fluidState_.setPressure(liquidPhaseIdx, problem.nonwettingReferencePressure() - pc);
- fluidState_.setPressure(gasPhaseIdx, problem.nonwettingReferencePressure());
+ // get pc for sw = 0.0
+ const Scalar pc = fluidMatrixInteraction.pc(0.0);
- // set molar densities
- if constexpr (enableWaterDiffusionInAir())
- {
- molarDensity_[liquidPhaseIdx] = FluidSystem::H2O::liquidDensity(temperature(), pressure(liquidPhaseIdx))/FluidSystem::H2O::molarMass();
- molarDensity_[gasPhaseIdx] = IdealGas<Scalar>::molarDensity(temperature(), problem.nonwettingReferencePressure());
- }
+ // set the wetting pressure
+ fluidState_.setPressure(liquidPhaseIdx, problem.nonwettingReferencePressure() - pc);
+ fluidState_.setPressure(gasPhaseIdx, problem.nonwettingReferencePressure());
- // density and viscosity
- paramCache.updateAll(fluidState_);
- fluidState_.setDensity(liquidPhaseIdx, FluidSystem::density(fluidState_, paramCache, liquidPhaseIdx));
- fluidState_.setDensity(gasPhaseIdx, FluidSystem::density(fluidState_, paramCache, gasPhaseIdx));
- fluidState_.setViscosity(liquidPhaseIdx, FluidSystem::viscosity(fluidState_, paramCache, liquidPhaseIdx));
-
- // compute and set the enthalpy
- fluidState_.setEnthalpy(liquidPhaseIdx, EnergyVolVars::enthalpy(fluidState_, paramCache, liquidPhaseIdx));
- fluidState_.setEnthalpy(gasPhaseIdx, EnergyVolVars::enthalpy(fluidState_, paramCache, gasPhaseIdx));
-
- //binary diffusion coefficients
- paramCache.updateAll(fluidState_);
- effectiveDiffCoeff_ = getEffectiveDiffusionCoefficient(gasPhaseIdx,
- FluidSystem::comp1Idx,
- FluidSystem::comp0Idx);
- }
- else if (phasePresence == Indices::bothPhases)
- {
- completeFluidState(elemSol, problem, element, scv, fluidState_, solidState_);
+ // set molar densities
+ if constexpr (enableWaterDiffusionInAir())
+ {
+ molarDensity_[liquidPhaseIdx] = FluidSystem::H2O::liquidDensity(temperature(), pressure(liquidPhaseIdx))/FluidSystem::H2O::molarMass();
+ molarDensity_[gasPhaseIdx] = IdealGas<Scalar>::molarDensity(temperature(), problem.nonwettingReferencePressure());
+ }
- // if we want to account for diffusion in the air phase
- // use Raoult to compute the water mole fraction in air
- if constexpr (enableWaterDiffusionInAir())
- {
- molarDensity_[liquidPhaseIdx] = FluidSystem::H2O::liquidDensity(temperature(), pressure(liquidPhaseIdx))/FluidSystem::H2O::molarMass();
- molarDensity_[gasPhaseIdx] = IdealGas<Scalar>::molarDensity(temperature(), problem.nonwettingReferencePressure());
- moleFraction_[liquidPhaseIdx] = 1.0;
+ // density and viscosity
+ paramCache.updateAll(fluidState_);
+ fluidState_.setDensity(liquidPhaseIdx, FluidSystem::density(fluidState_, paramCache, liquidPhaseIdx));
+ fluidState_.setDensity(gasPhaseIdx, FluidSystem::density(fluidState_, paramCache, gasPhaseIdx));
+ fluidState_.setViscosity(liquidPhaseIdx, FluidSystem::viscosity(fluidState_, paramCache, liquidPhaseIdx));
- moleFraction_[gasPhaseIdx] = FluidSystem::H2O::vaporPressure(temperature()) / problem.nonwettingReferencePressure();
+ // compute and set the enthalpy
+ fluidState_.setEnthalpy(liquidPhaseIdx, EnergyVolVars::enthalpy(fluidState_, paramCache, liquidPhaseIdx));
+ fluidState_.setEnthalpy(gasPhaseIdx, EnergyVolVars::enthalpy(fluidState_, paramCache, gasPhaseIdx));
- const auto averageMolarMassGasPhase = (moleFraction_[gasPhaseIdx]*FluidSystem::molarMass(liquidPhaseIdx)) +
- ((1-moleFraction_[gasPhaseIdx])*FluidSystem::molarMass(gasPhaseIdx));
-
- //X_w = x_w* M_w/ M_avg
- massFraction_[gasPhaseIdx] = moleFraction_[gasPhaseIdx]*FluidSystem::molarMass(liquidPhaseIdx)/averageMolarMassGasPhase;
-
- // binary diffusion coefficients
+ //binary diffusion coefficients
paramCache.updateAll(fluidState_);
effectiveDiffCoeff_ = getEffectiveDiffusionCoefficient(gasPhaseIdx,
FluidSystem::comp1Idx,
FluidSystem::comp0Idx);
}
+ else if (phasePresence == Indices::bothPhases)
+ {
+ completeFluidState(elemSol, problem, element, scv, fluidState_, solidState_);
+
+ // if we want to account for diffusion in the air phase
+ // use Raoult to compute the water mole fraction in air
+ if constexpr (enableWaterDiffusionInAir())
+ {
+ molarDensity_[liquidPhaseIdx] = FluidSystem::H2O::liquidDensity(temperature(), pressure(liquidPhaseIdx))/FluidSystem::H2O::molarMass();
+ molarDensity_[gasPhaseIdx] = IdealGas<Scalar>::molarDensity(temperature(), problem.nonwettingReferencePressure());
+ moleFraction_[liquidPhaseIdx] = 1.0;
+
+ moleFraction_[gasPhaseIdx] = FluidSystem::H2O::vaporPressure(temperature()) / problem.nonwettingReferencePressure();
+
+ const auto averageMolarMassGasPhase = (moleFraction_[gasPhaseIdx]*FluidSystem::molarMass(liquidPhaseIdx)) +
+ ((1-moleFraction_[gasPhaseIdx])*FluidSystem::molarMass(gasPhaseIdx));
+
+ //X_w = x_w* M_w/ M_avg
+ massFraction_[gasPhaseIdx] = moleFraction_[gasPhaseIdx]*FluidSystem::molarMass(liquidPhaseIdx)/averageMolarMassGasPhase;
+
+ // binary diffusion coefficients
+ paramCache.updateAll(fluidState_);
+ effectiveDiffCoeff_ = getEffectiveDiffusionCoefficient(gasPhaseIdx,
+ FluidSystem::comp1Idx,
+ FluidSystem::comp0Idx);
+ }
+ }
+ else if (phasePresence == Indices::liquidPhaseOnly)
+ {
+ completeFluidState(elemSol, problem, element, scv, fluidState_, solidState_);
+
+ if constexpr (enableWaterDiffusionInAir())
+ {
+ molarDensity_[liquidPhaseIdx] = FluidSystem::H2O::liquidDensity(temperature(), pressure(liquidPhaseIdx))/FluidSystem::H2O::molarMass();
+ molarDensity_[gasPhaseIdx] = IdealGas<Scalar>::molarDensity(temperature(), problem.nonwettingReferencePressure());
+ moleFraction_[liquidPhaseIdx] = 1.0;
+ moleFraction_[gasPhaseIdx] = 0.0;
+ massFraction_[liquidPhaseIdx] = 1.0;
+ massFraction_[gasPhaseIdx] = 0.0;
+
+ // binary diffusion coefficients (none required for liquid phase only)
+ effectiveDiffCoeff_ = 0.0;
+ }
+ }
}
- else if (phasePresence == Indices::liquidPhaseOnly)
+ else
{
completeFluidState(elemSol, problem, element, scv, fluidState_, solidState_);
-
- if constexpr (enableWaterDiffusionInAir())
- {
- molarDensity_[liquidPhaseIdx] = FluidSystem::H2O::liquidDensity(temperature(), pressure(liquidPhaseIdx))/FluidSystem::H2O::molarMass();
- molarDensity_[gasPhaseIdx] = IdealGas<Scalar>::molarDensity(temperature(), problem.nonwettingReferencePressure());
- moleFraction_[liquidPhaseIdx] = 1.0;
- moleFraction_[gasPhaseIdx] = 0.0;
- massFraction_[liquidPhaseIdx] = 1.0;
- massFraction_[gasPhaseIdx] = 0.0;
-
- // binary diffusion coefficients (none required for liquid phase only)
- effectiveDiffCoeff_ = 0.0;
- }
}
//////////
diff --git a/dumux/porousmediumflow/richardsextended/CMakeLists.txt b/dumux/porousmediumflow/richardsextended/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..8f8e866fc5acc387f42abae50e178746e52ef18f
--- /dev/null
+++ b/dumux/porousmediumflow/richardsextended/CMakeLists.txt
@@ -0,0 +1,3 @@
+file(GLOB DUMUX_POROUSMEDIUMFLOW_RICHARDSEXTENDED_HEADERS *.hh *.inc)
+install(FILES ${DUMUX_POROUSMEDIUMFLOW_RICHARDSEXTENDED_HEADERS}
+ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/porousmediumflow/richardsextended)
diff --git a/dumux/porousmediumflow/richardsextended/indices.hh b/dumux/porousmediumflow/richardsextended/indices.hh
new file mode 100644
index 0000000000000000000000000000000000000000..8a0fc76475d62077526d6b8c58bd20422b03dfa6
--- /dev/null
+++ b/dumux/porousmediumflow/richardsextended/indices.hh
@@ -0,0 +1,53 @@
+// -*- 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 3 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
+ * \ingroup ExtendedRichardsModel
+ * \brief Index names for the extended Richards model.
+ */
+
+#ifndef DUMUX_RICHARDSEXTENDED_INDICES_HH
+#define DUMUX_RICHARDSEXTENDED_INDICES_HH
+
+namespace Dumux {
+
+/*!
+ * \ingroup ExtendedRichardsModel
+ * \brief Index names for the extended Richards model.
+ */
+
+struct ExtendedRichardsIndices
+{
+ //! Primary variable index for the wetting phase pressure
+ static constexpr int pressureIdx = 0;
+ static constexpr int switchIdx = 0;
+
+ //! Equation index for the mass conservation of the wetting phase
+ static constexpr int conti0EqIdx = 0;
+
+ //TODO: After release (3.6) replace with inheritance from base Richards model plus phase states
+ // present phases (-> 'pseudo' primary variable)
+ static constexpr int liquidPhaseOnly = 1; //!< Only the liquid phase is present
+ static constexpr int gasPhaseOnly = 2; //!< Only the gas phase is present
+ static constexpr int bothPhases = 3; //!< Both phases are present
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/porousmediumflow/richardsextended/iofields.hh b/dumux/porousmediumflow/richardsextended/iofields.hh
new file mode 100644
index 0000000000000000000000000000000000000000..eb1e6db6be9729ef691a8848e007373ba65765f6
--- /dev/null
+++ b/dumux/porousmediumflow/richardsextended/iofields.hh
@@ -0,0 +1,83 @@
+// -*- 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 3 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
+ * \ingroup ExtendedRichardsModel
+ * \brief Adds I/O fields specific to the extended Richards model.
+ */
+
+#ifndef DUMUX_RICHARDSEXTENDED_IO_FIELDS_HH
+#define DUMUX_RICHARDSEXTENDED_IO_FIELDS_HH
+
+#include <dumux/io/name.hh>
+#include <dumux/porousmediumflow/richards/iofields.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup ExtendedRichardsModel
+ * \brief Adds I/O fields specific to the extended Richards model.
+ */
+class ExtendedRichardsIOFields : public RichardsIOFields<false>
+{
+ using ParentType = RichardsIOFields<false>;
+public:
+ template <class OutputModule>
+ static void initOutputModule(OutputModule& out)
+ {
+ using VV = typename OutputModule::VolumeVariables;
+ using FS = typename VV::FluidSystem;
+
+ // TODO: replace by a call to the base class plus code specific to extended model after release (3.6)
+ out.addVolumeVariable([](const auto& v){ return v.saturation(FS::phase0Idx); },
+ IOName::saturation<FS>(FS::phase0Idx));
+ out.addVolumeVariable([](const auto& v){ return v.saturation(FS::phase1Idx); },
+ IOName::saturation<FS>(FS::phase1Idx));
+ out.addVolumeVariable([](const auto& v){ return v.pressure(FS::phase0Idx); },
+ IOName::pressure<FS>(FS::phase0Idx));
+ out.addVolumeVariable([](const auto& v){ return v.pressure(FS::phase1Idx); },
+ IOName::pressure<FS>(FS::phase1Idx));
+ out.addVolumeVariable([](const auto& v){ return v.capillaryPressure(); },
+ IOName::capillaryPressure());
+ out.addVolumeVariable([](const auto& v){ return v.density(FS::phase0Idx); },
+ IOName::density<FS>(FS::phase0Idx));
+ out.addVolumeVariable([](const auto& v){ return v.mobility(FS::phase0Idx); },
+ IOName::mobility<FS>(FS::phase0Idx));
+ out.addVolumeVariable([](const auto& v){ return v.relativePermeability(FS::phase0Idx); },
+ IOName::relativePermeability<FS>(FS::phase0Idx));
+ out.addVolumeVariable([](const auto& v){ return v.porosity(); },
+ IOName::porosity());
+
+ static const bool gravity = getParamFromGroup<bool>(out.paramGroup(), "Problem.EnableGravity");
+
+ if(gravity)
+ out.addVolumeVariable([](const auto& v){ return v.pressureHead(FS::phase0Idx); },
+ IOName::pressureHead());
+ out.addVolumeVariable([](const auto& v){ return v.moleFraction(FS::phase1Idx, FS::comp0Idx); },
+ IOName::moleFraction<FS>(FS::phase1Idx, FS::comp0Idx));
+ out.addVolumeVariable([](const auto& v){ return v.waterContent(FS::phase0Idx); },
+ IOName::waterContent());
+ out.addVolumeVariable([](const auto& v){ return v.priVars().state(); },
+ IOName::phasePresence());
+ }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/porousmediumflow/richardsextended/localresidual.hh b/dumux/porousmediumflow/richardsextended/localresidual.hh
new file mode 100644
index 0000000000000000000000000000000000000000..8003ce641ebcc2585c00fd4396a99522c9d3fed0
--- /dev/null
+++ b/dumux/porousmediumflow/richardsextended/localresidual.hh
@@ -0,0 +1,160 @@
+// -*- 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 3 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
+ * \ingroup ExtendedRichardsModel
+ * \brief Element-wise calculation of the Jacobian matrix for problems
+ * using the extended Richards fully implicit models.
+ */
+
+#ifndef DUMUX_RICHARDSEXTENDED_LOCAL_RESIDUAL_HH
+#define DUMUX_RICHARDSEXTENDED_LOCAL_RESIDUAL_HH
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/typetraits/typetraits.hh>
+#include <dumux/common/numeqvector.hh>
+#include <dumux/discretization/method.hh>
+#include <dumux/discretization/extrusion.hh>
+#include <dumux/flux/referencesystemformulation.hh>
+#include <dumux/porousmediumflow/richards/localresidual.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup ExtendedRichardsModel
+ * \brief Element-wise calculation of the Jacobian matrix for problems
+ * using the extended Richards fully implicit models.
+ */
+template<class TypeTag>
+class ExtendedRichardsLocalResidual : public RichardsLocalResidual<TypeTag>
+{
+ using Implementation = GetPropType<TypeTag, Properties::LocalResidual>;
+
+ using ParentType = RichardsLocalResidual<TypeTag>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using Problem = GetPropType<TypeTag, Properties::Problem>;
+ using NumEqVector = Dumux::NumEqVector<GetPropType<TypeTag, Properties::PrimaryVariables>>;
+ using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
+ using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
+ using FluxVariables = GetPropType<TypeTag, Properties::FluxVariables>;
+ using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView;
+ using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+ using FVElementGeometry = typename GridGeometry::LocalView;
+ using SubControlVolume = typename GridGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
+ using GridView = typename GridGeometry::GridView;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using EnergyLocalResidual = GetPropType<TypeTag, Properties::EnergyLocalResidual>;
+ using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+ using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+ // first index for the mass balance
+ enum { conti0EqIdx = Indices::conti0EqIdx };
+
+ // phase & component indices
+ static constexpr auto liquidPhaseIdx = FluidSystem::phase0Idx;
+ static constexpr auto gasPhaseIdx = FluidSystem::phase1Idx;
+ static constexpr auto liquidCompIdx = FluidSystem::comp0Idx;
+
+public:
+ using ParentType::ParentType;
+
+ /*!
+ * \brief Evaluates 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
+ */
+ NumEqVector computeStorage(const Problem& problem,
+ const SubControlVolume& scv,
+ const VolumeVariables& volVars) const
+ {
+ // TODO: replace by a call to the bass class plus code specific to extended model after release (3.6)
+ // partial time derivative of the phase mass
+ NumEqVector storage(0.0);
+ storage[conti0EqIdx] = volVars.porosity()
+ * volVars.density(liquidPhaseIdx)
+ * volVars.saturation(liquidPhaseIdx);
+
+ // for extended Richards we consider water in air
+ storage[conti0EqIdx] += volVars.porosity()
+ * volVars.molarDensity(gasPhaseIdx)
+ * volVars.moleFraction(gasPhaseIdx, liquidCompIdx)
+ * FluidSystem::molarMass(liquidCompIdx)
+ * volVars.saturation(gasPhaseIdx);
+
+ //! The energy storage in the water, air and solid phase
+ EnergyLocalResidual::fluidPhaseStorage(storage, scv, volVars, liquidPhaseIdx);
+ EnergyLocalResidual::fluidPhaseStorage(storage, scv, volVars, gasPhaseIdx);
+ EnergyLocalResidual::solidPhaseStorage(storage, scv, volVars);
+
+ return storage;
+ }
+
+
+ /*!
+ * \brief Evaluates the mass flux over a face of a sub control volume.
+ *
+ * \param problem The problem
+ * \param element The current element.
+ * \param fvGeometry The finite-volume geometry
+ * \param elemVolVars The volume variables of the current element
+ * \param scvf The sub control volume face to compute the flux on
+ * \param elemFluxVarsCache The cache related to flux computation
+ */
+ NumEqVector computeFlux(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf,
+ const ElementFluxVariablesCache& elemFluxVarsCache) const
+ {
+ // TODO: replace by a call to the bass class plus code specific to extended model after release (3.6)
+ FluxVariables fluxVars;
+ fluxVars.init(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
+
+ NumEqVector flux(0.0);
+ // the physical quantities for which we perform upwinding
+ auto upwindTerm = [](const auto& volVars)
+ { return volVars.density(liquidPhaseIdx)*volVars.mobility(liquidPhaseIdx); };
+
+ flux[conti0EqIdx] = fluxVars.advectiveFlux(liquidPhaseIdx, upwindTerm);
+
+ // for extended Richards we consider water vapor diffusion in air
+ //check for the reference system and adapt units of the diffusive flux accordingly.
+ if (FluxVariables::MolecularDiffusionType::referenceSystemFormulation() == ReferenceSystemFormulation::massAveraged)
+ flux[conti0EqIdx] += fluxVars.molecularDiffusionFlux(gasPhaseIdx)[liquidCompIdx];
+ else
+ flux[conti0EqIdx] += fluxVars.molecularDiffusionFlux(gasPhaseIdx)[liquidCompIdx]*FluidSystem::molarMass(liquidCompIdx);
+
+ //! Add advective phase energy fluxes for the water phase only. For isothermal model the contribution is zero.
+ EnergyLocalResidual::heatConvectionFlux(flux, fluxVars, liquidPhaseIdx);
+
+ //! Add diffusive energy fluxes. For isothermal model the contribution is zero.
+ //! The effective lambda is averaged over both fluid phases and the solid phase
+ EnergyLocalResidual::heatConductionFlux(flux, fluxVars);
+
+ return flux;
+ }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/porousmediumflow/richardsextended/model.hh b/dumux/porousmediumflow/richardsextended/model.hh
new file mode 100644
index 0000000000000000000000000000000000000000..58c4b161b79c29c7039394a26f2c575faf984b1e
--- /dev/null
+++ b/dumux/porousmediumflow/richardsextended/model.hh
@@ -0,0 +1,316 @@
+// -*- 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 3 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
+ * \ingroup ExtendedRichardsModel
+ * \brief This model implements a variant of the extended Richards'
+ * equation for quasi-twophase flow (see e.g. Vanderborght et al. 2017).
+ *
+ * In the unsaturated zone, Richards' equation
+ \f[
+ \frac{\partial\;\phi S_w \varrho_w}{\partial t}
+ +
+ \frac{\partial\;\phi (1-S_w)\varrho_n X_n^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)
+ +
+ {\bf D_{n, pm}^w} \varrho_n \text{grad}\, X^w_n
+ \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 nonwetting 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 infinitely 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 chooses 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!
+ */
+
+#ifndef DUMUX_RICHARDSEXTENDED_MODEL_HH
+#define DUMUX_RICHARDSEXTENDED_MODEL_HH
+
+#include <dune/common/fvector.hh>
+
+#include <dumux/common/properties.hh>
+
+#include <dumux/porousmediumflow/immiscible/localresidual.hh>
+#include <dumux/porousmediumflow/compositional/switchableprimaryvariables.hh>
+#include <dumux/material/fluidmatrixinteractions/diffusivitymillingtonquirk.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/thermalconductivity/somerton.hh>
+#include <dumux/material/components/simpleh2o.hh>
+#include <dumux/material/fluidsystems/h2oair.hh>
+#include <dumux/material/fluidstates/immiscible.hh>
+
+#include <dumux/porousmediumflow/properties.hh>
+#include <dumux/porousmediumflow/nonisothermal/model.hh>
+#include <dumux/porousmediumflow/nonisothermal/indices.hh>
+#include <dumux/porousmediumflow/nonisothermal/iofields.hh>
+#include <dumux/porousmediumflow/richards/balanceequationopts.hh>
+#include <dumux/porousmediumflow/richards/model.hh>
+#include <dumux/porousmediumflow/richards/velocityoutput.hh>
+
+#include "indices.hh"
+#include "volumevariables.hh"
+#include "iofields.hh"
+#include "localresidual.hh"
+
+namespace Dumux {
+
+/*!
+ * \ingroup ExtendedRichardsModel
+ * \brief Specifies a number properties of the extended Richards model.
+ */
+struct ExtendedRichardsModelTraits : public RichardsModelTraits<false>
+{
+ using Indices = ExtendedRichardsIndices;
+
+ static constexpr int numEq() { return 1; }
+ static constexpr int numFluidPhases() { return 2; }
+ static constexpr int numFluidComponents() { return 1; }
+
+ static constexpr bool enableAdvection() { return true; }
+ static constexpr bool enableMolecularDiffusion() { return true; }
+ static constexpr bool enableEnergyBalance() { return false; }
+
+ //! The Richards model has some assumptions on the fluid systems
+ //! that can be verified with this trait
+ template<class FluidSystem>
+ static constexpr bool fluidSystemIsCompatible()
+ {
+ return !FluidSystem::isGas(FluidSystem::phase0Idx)
+ && FluidSystem::isGas(FluidSystem::phase1Idx);
+ }
+
+ //! The Richards model has some assumptions on the fluid systems
+ //! that can be verified with this trait
+ template<class FluidSystem>
+ static constexpr auto checkFluidSystem(const FluidSystem& fs)
+ {
+ struct FluidSystemCheck {
+ static_assert(fluidSystemIsCompatible<FluidSystem>(),
+ "Richards model currently assumes the first phase to be liquid and the second phase to be gaseous.");
+ };
+ return FluidSystemCheck{};
+ }
+};
+
+/*!
+ * \ingroup RichardsModel
+ * \brief Traits class for the Richards model.
+ *
+ * \tparam PV The type used for primary variables
+ * \tparam FSY The fluid system type
+ * \tparam FST The fluid state type
+ * \tparam PT The type used for permeabilities
+ * \tparam MT The model traits
+ */
+template<class PV, class FSY, class FST, class SSY, class SST, class PT, class MT, class DT, class EDM>
+struct ExtendedRichardsVolumeVariablesTraits
+{
+ using PrimaryVariables = PV;
+ using FluidSystem = FSY;
+ using FluidState = FST;
+ using SolidSystem = SSY;
+ using SolidState = SST;
+ using PermeabilityType = PT;
+ using ModelTraits = MT;
+ using DiffusionType = DT;
+ using EffectiveDiffusivityModel = EDM;
+};
+
+// \{
+///////////////////////////////////////////////////////////////////////////
+// properties for the isothermal Richards model.
+///////////////////////////////////////////////////////////////////////////
+namespace Properties {
+
+//////////////////////////////////////////////////////////////////
+// Type tags
+//////////////////////////////////////////////////////////////////
+
+//! The type tags for the implicit isothermal one-phase two-component problems
+// Create new type tags
+namespace TTag {
+struct ExtendedRichards { using InheritsFrom = std::tuple<Richards>; };
+struct ExtendedRichardsNI { using InheritsFrom = std::tuple<ExtendedRichards>; };
+} // end namespace TTag
+
+//////////////////////////////////////////////////////////////////
+// Properties values
+//////////////////////////////////////////////////////////////////
+
+//! The local residual operator
+template<class TypeTag>
+struct LocalResidual<TypeTag, TTag::ExtendedRichards> { using type = ExtendedRichardsLocalResidual<TypeTag>; };
+
+//! Set the vtk output fields specific to this model
+template<class TypeTag>
+struct IOFields<TypeTag, TTag::ExtendedRichards>
+{
+ using type = ExtendedRichardsIOFields;
+};
+
+//! The model traits
+template<class TypeTag>
+struct ModelTraits<TypeTag, TTag::ExtendedRichards> { using type = ExtendedRichardsModelTraits; };
+
+//! Set the volume variables property
+template<class TypeTag>
+struct VolumeVariables<TypeTag, TTag::ExtendedRichards>
+{
+private:
+ using PV = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ using FSY = GetPropType<TypeTag, Properties::FluidSystem>;
+ using FST = GetPropType<TypeTag, Properties::FluidState>;
+ using SSY = GetPropType<TypeTag, Properties::SolidSystem>;
+ using SST = GetPropType<TypeTag, Properties::SolidState>;
+ using PT = typename GetPropType<TypeTag, Properties::SpatialParams>::PermeabilityType;
+ using MT = GetPropType<TypeTag, Properties::ModelTraits>;
+ using DT = GetPropType<TypeTag, Properties::MolecularDiffusionType>;
+ using EDM = GetPropType<TypeTag, Properties::EffectiveDiffusivityModel>;
+ using Traits = ExtendedRichardsVolumeVariablesTraits<PV, FSY, FST, SSY, SST, PT, MT, DT, EDM>;
+public:
+ using type = ExtendedRichardsVolumeVariables<Traits>;
+};
+
+//! Use the model after Millington (1961) for the effective diffusivity
+template<class TypeTag>
+struct EffectiveDiffusivityModel<TypeTag, TTag::ExtendedRichards>
+{ using type = DiffusivityMillingtonQuirk<GetPropType<TypeTag, Properties::Scalar>>; };
+
+//! The primary variables vector for the richards model
+template<class TypeTag>
+struct PrimaryVariables<TypeTag, TTag::ExtendedRichards>
+{
+private:
+ using PrimaryVariablesVector = Dune::FieldVector<GetPropType<TypeTag, Properties::Scalar>,
+ GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
+public:
+ using type = SwitchablePrimaryVariables<PrimaryVariablesVector, int>;
+};
+
+/////////////////////////////////////////////////////
+// Property values for non-isothermal Richars model
+/////////////////////////////////////////////////////
+
+//! set non-isothermal model traits
+template<class TypeTag>
+struct ModelTraits<TypeTag, TTag::ExtendedRichardsNI>
+{
+private:
+ using IsothermalTraits = ExtendedRichardsModelTraits;
+public:
+ using type = PorousMediumFlowNIModelTraits<IsothermalTraits>;
+};
+
+//! Set the vtk output fields specific to th non-isothermal model
+template<class TypeTag>
+struct IOFields<TypeTag, TTag::ExtendedRichardsNI>
+{
+ using type = EnergyIOFields<ExtendedRichardsIOFields>;
+};
+
+//! Set the volume variables property
+template<class TypeTag>
+struct VolumeVariables<TypeTag, TTag::ExtendedRichardsNI>
+{
+private:
+ using PV = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ using FSY = GetPropType<TypeTag, Properties::FluidSystem>;
+ using FST = GetPropType<TypeTag, Properties::FluidState>;
+ using SSY = GetPropType<TypeTag, Properties::SolidSystem>;
+ using SST = GetPropType<TypeTag, Properties::SolidState>;
+ using PT = typename GetPropType<TypeTag, Properties::SpatialParams>::PermeabilityType;
+ using MT = GetPropType<TypeTag, Properties::ModelTraits>;
+ using DT = GetPropType<TypeTag, Properties::MolecularDiffusionType>;
+ using EDM = GetPropType<TypeTag, Properties::EffectiveDiffusivityModel>;
+ using BaseTraits = ExtendedRichardsVolumeVariablesTraits<PV, FSY, FST, SSY, SST, PT, MT, DT, EDM>;
+
+ using ETCM = GetPropType< TypeTag, Properties::ThermalConductivityModel>;
+ template<class BaseTraits, class ETCM>
+ struct NITraits : public BaseTraits { using EffectiveThermalConductivityModel = ETCM; };
+
+public:
+ using type = ExtendedRichardsVolumeVariables<NITraits<BaseTraits, ETCM>>;
+};
+
+//! Somerton is used as default model to compute the effective thermal heat conductivity
+template<class TypeTag>
+struct ThermalConductivityModel<TypeTag, TTag::ExtendedRichardsNI>
+{
+private:
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+public:
+ using type = ThermalConductivitySomerton<Scalar>;
+};
+
+// \}
+} // end namespace Properties
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/porousmediumflow/richards/primaryvariableswitch.hh b/dumux/porousmediumflow/richardsextended/primaryvariableswitch.hh
similarity index 94%
rename from dumux/porousmediumflow/richards/primaryvariableswitch.hh
rename to dumux/porousmediumflow/richardsextended/primaryvariableswitch.hh
index 0c8ea2307689a315337691e0a123fe46ed26f42f..dfbc638c1956f45ade1c283c68504c7ea3efd897 100644
--- a/dumux/porousmediumflow/richards/primaryvariableswitch.hh
+++ b/dumux/porousmediumflow/richardsextended/primaryvariableswitch.hh
@@ -18,12 +18,12 @@
*****************************************************************************/
/*!
* \file
- * \ingroup RichardsModel
+ * \ingroup ExtenndedRichardsModel
* \brief The primary variable switch for the extended Richards model.
*/
-#ifndef DUMUX_RICHARDS_PRIMARY_VARIABLE_SWITCH_HH
-#define DUMUX_RICHARDS_PRIMARY_VARIABLE_SWITCH_HH
+#ifndef DUMUX_RICHARDSEXTENDED_PRIMARY_VARIABLE_SWITCH_HH
+#define DUMUX_RICHARDSEXTENDED_PRIMARY_VARIABLE_SWITCH_HH
#include <dumux/common/exceptions.hh>
#include <dumux/common/parameters.hh>
@@ -62,9 +62,6 @@ protected:
if (!usePriVarSwitch)
return false;
- if (!VolumeVariables::enableWaterDiffusionInAir())
- DUNE_THROW(Dune::InvalidStateException, "The Richards primary variable switch only works with water diffusion in air enabled!");
-
static constexpr int liquidCompIdx = FluidSystem::liquidPhaseIdx;
// evaluate primary variable switch
diff --git a/dumux/porousmediumflow/richardsextended/volumevariables.hh b/dumux/porousmediumflow/richardsextended/volumevariables.hh
new file mode 100644
index 0000000000000000000000000000000000000000..d69ad5a256f6c3718c511d9a591f0a0ba961f22b
--- /dev/null
+++ b/dumux/porousmediumflow/richardsextended/volumevariables.hh
@@ -0,0 +1,498 @@
+// -*- 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 3 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
+ * \ingroup ExtendedRichardsModel
+ * \brief Volume averaged quantities required by the extended Richards model.
+ */
+
+#ifndef DUMUX_RICHARDSEXTENDED_VOLUME_VARIABLES_HH
+#define DUMUX_RICHARDSEXTENDED_VOLUME_VARIABLES_HH
+
+#include <cassert>
+
+#include <dune/common/exceptions.hh>
+
+#include <dumux/porousmediumflow/volumevariables.hh>
+#include <dumux/porousmediumflow/nonisothermal/volumevariables.hh>
+#include <dumux/material/idealgas.hh>
+#include <dumux/material/constants.hh>
+#include <dumux/material/solidstates/updatesolidvolumefractions.hh>
+
+#include "primaryvariableswitch.hh"
+
+namespace Dumux {
+
+/*!
+ * \ingroup ExtendedRichardsModel
+ * \brief Volume averaged quantities required by the extended Richards model.
+ *
+ * This contains the quantities which are are constant within a finite
+ * volume in the Richards model.
+ */
+template <class Traits>
+class ExtendedRichardsVolumeVariables
+: public PorousMediumFlowVolumeVariables<Traits>
+, public EnergyVolumeVariables<Traits, ExtendedRichardsVolumeVariables<Traits> >
+{
+ using ParentType = PorousMediumFlowVolumeVariables<Traits>;
+ using EnergyVolVars = EnergyVolumeVariables<Traits, ExtendedRichardsVolumeVariables<Traits> >;
+ using Scalar = typename Traits::PrimaryVariables::value_type;
+ using PermeabilityType = typename Traits::PermeabilityType;
+ using ModelTraits = typename Traits::ModelTraits;
+
+ static constexpr int numFluidComps = ParentType::numFluidComponents();
+ static constexpr int numPhases = ParentType::numFluidPhases();
+
+ using EffDiffModel = typename Traits::EffectiveDiffusivityModel;
+
+ // checks if the fluid system uses the Richards model index convention
+ static constexpr auto fsCheck = ModelTraits::checkFluidSystem(typename Traits::FluidSystem{});
+
+public:
+ //! Export type of the fluid system
+ using FluidSystem = typename Traits::FluidSystem;
+ //! Export type of the fluid state
+ using FluidState = typename Traits::FluidState;
+ //! Export type of the fluid state
+ //! Export type of solid state
+ using SolidState = typename Traits::SolidState;
+ //! Export type of solid system
+ using SolidSystem = typename Traits::SolidSystem;
+ using Indices = typename Traits::ModelTraits::Indices;
+ using PrimaryVariableSwitch = ExtendedRichardsPrimaryVariableSwitch;
+
+ //! Export phase indices
+ static constexpr auto liquidPhaseIdx = Traits::FluidSystem::phase0Idx;
+ static constexpr auto gasPhaseIdx = Traits::FluidSystem::phase1Idx;
+
+ /*!
+ * \brief Updates all quantities for a given control volume.
+ *
+ * \param elemSol A vector containing all primary variables connected to the element
+ * \param problem The object specifying the problem which ought to
+ * be simulated
+ * \param element An element which contains part of the control volume
+ * \param scv The sub-control volume
+ */
+ template<class ElemSol, class Problem, class Element, class Scv>
+ void update(const ElemSol &elemSol,
+ const Problem &problem,
+ const Element &element,
+ const Scv& scv)
+ {
+ ParentType::update(elemSol, problem, element, scv);
+
+ const auto fluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element, scv, elemSol);
+
+ const auto& priVars = elemSol[scv.localDofIndex()];
+ const auto phasePresence = priVars.state();
+
+ // precompute the minimum capillary pressure (entry pressure)
+ // needed to make sure we don't compute unphysical capillary pressures and thus saturations
+ minPc_ = fluidMatrixInteraction.endPointPc();
+
+ //update porosity before calculating the effective properties depending on it
+ updateSolidVolumeFractions(elemSol, problem, element, scv, solidState_, numFluidComps);
+
+ typename FluidSystem::ParameterCache paramCache;
+ auto getEffectiveDiffusionCoefficient = [&](int phaseIdx, int compIIdx, int compJIdx)
+ {
+ return EffDiffModel::effectiveDiffusionCoefficient(*this, phaseIdx, compIIdx, compJIdx);
+ };
+
+ if (phasePresence == Indices::gasPhaseOnly)
+ {
+ moleFraction_[liquidPhaseIdx] = 1.0;
+ massFraction_[liquidPhaseIdx] = 1.0;
+ moleFraction_[gasPhaseIdx] = priVars[Indices::switchIdx];
+
+ const auto averageMolarMassGasPhase = (moleFraction_[gasPhaseIdx]*FluidSystem::molarMass(liquidPhaseIdx)) +
+ ((1-moleFraction_[gasPhaseIdx])*FluidSystem::molarMass(gasPhaseIdx));
+
+ //X_w = x_w* M_w/ M_avg
+ massFraction_[gasPhaseIdx] = priVars[Indices::switchIdx]*FluidSystem::molarMass(liquidPhaseIdx)/averageMolarMassGasPhase;
+
+ fluidState_.setSaturation(liquidPhaseIdx, 0.0);
+ fluidState_.setSaturation(gasPhaseIdx, 1.0);
+
+ EnergyVolVars::updateTemperature(elemSol, problem, element, scv, fluidState_, solidState_);
+
+ // get pc for sw = 0.0
+ const Scalar pc = fluidMatrixInteraction.pc(0.0);
+
+ // set the wetting pressure
+ fluidState_.setPressure(liquidPhaseIdx, problem.nonwettingReferencePressure() - pc);
+ fluidState_.setPressure(gasPhaseIdx, problem.nonwettingReferencePressure());
+
+ molarDensity_[liquidPhaseIdx] = FluidSystem::H2O::liquidDensity(temperature(), pressure(liquidPhaseIdx))/FluidSystem::H2O::molarMass();
+ molarDensity_[gasPhaseIdx] = IdealGas<Scalar>::molarDensity(temperature(), problem.nonwettingReferencePressure());
+
+ // density and viscosity
+ paramCache.updateAll(fluidState_);
+ fluidState_.setDensity(liquidPhaseIdx, FluidSystem::density(fluidState_, paramCache, liquidPhaseIdx));
+ fluidState_.setDensity(gasPhaseIdx, FluidSystem::density(fluidState_, paramCache, gasPhaseIdx));
+ fluidState_.setViscosity(liquidPhaseIdx, FluidSystem::viscosity(fluidState_, paramCache, liquidPhaseIdx));
+
+ // compute and set the enthalpy
+ fluidState_.setEnthalpy(liquidPhaseIdx, EnergyVolVars::enthalpy(fluidState_, paramCache, liquidPhaseIdx));
+ fluidState_.setEnthalpy(gasPhaseIdx, EnergyVolVars::enthalpy(fluidState_, paramCache, gasPhaseIdx));
+
+ //binary diffusion coefficients
+ paramCache.updateAll(fluidState_);
+ effectiveDiffCoeff_ = getEffectiveDiffusionCoefficient(gasPhaseIdx,
+ FluidSystem::comp1Idx,
+ FluidSystem::comp0Idx);
+ }
+ else if (phasePresence == Indices::bothPhases)
+ {
+ completeFluidState_(elemSol, problem, element, scv, fluidState_, solidState_);
+
+ // we want to account for diffusion in the air phase
+ // use Raoult to compute the water mole fraction in air
+ molarDensity_[liquidPhaseIdx] = FluidSystem::H2O::liquidDensity(temperature(), pressure(liquidPhaseIdx))/FluidSystem::H2O::molarMass();
+ molarDensity_[gasPhaseIdx] = IdealGas<Scalar>::molarDensity(temperature(), problem.nonwettingReferencePressure());
+ moleFraction_[liquidPhaseIdx] = 1.0;
+
+ moleFraction_[gasPhaseIdx] = FluidSystem::H2O::vaporPressure(temperature()) / problem.nonwettingReferencePressure();
+
+ const auto averageMolarMassGasPhase = (moleFraction_[gasPhaseIdx]*FluidSystem::molarMass(liquidPhaseIdx)) +
+ ((1-moleFraction_[gasPhaseIdx])*FluidSystem::molarMass(gasPhaseIdx));
+
+ //X_w = x_w* M_w/ M_avg
+ massFraction_[gasPhaseIdx] = moleFraction_[gasPhaseIdx]*FluidSystem::molarMass(liquidPhaseIdx)/averageMolarMassGasPhase;
+
+ // binary diffusion coefficients
+ paramCache.updateAll(fluidState_);
+ effectiveDiffCoeff_ = getEffectiveDiffusionCoefficient(gasPhaseIdx,
+ FluidSystem::comp1Idx,
+ FluidSystem::comp0Idx);
+ }
+ else if (phasePresence == Indices::liquidPhaseOnly)
+ {
+ completeFluidState_(elemSol, problem, element, scv, fluidState_, solidState_);
+
+ molarDensity_[liquidPhaseIdx] = FluidSystem::H2O::liquidDensity(temperature(), pressure(liquidPhaseIdx))/FluidSystem::H2O::molarMass();
+ molarDensity_[gasPhaseIdx] = IdealGas<Scalar>::molarDensity(temperature(), problem.nonwettingReferencePressure());
+ moleFraction_[liquidPhaseIdx] = 1.0;
+ moleFraction_[gasPhaseIdx] = 0.0;
+ massFraction_[liquidPhaseIdx] = 1.0;
+ massFraction_[gasPhaseIdx] = 0.0;
+
+ // binary diffusion coefficients (none required for liquid phase only)
+ effectiveDiffCoeff_ = 0.0;
+ }
+
+ //////////
+ // specify the other parameters
+ //////////
+ relativePermeabilityWetting_ = fluidMatrixInteraction.krw(fluidState_.saturation(liquidPhaseIdx));
+ EnergyVolVars::updateSolidEnergyParams(elemSol, problem, element, scv, solidState_);
+ permeability_ = problem.spatialParams().permeability(element, scv, elemSol);
+ EnergyVolVars::updateEffectiveThermalConductivity();
+ }
+
+ /*!
+ * \brief Returns the fluid configuration at the given primary
+ * variables.
+ */
+ const FluidState &fluidState() const
+ { return fluidState_; }
+
+ /*!
+ * \brief Returns the phase state for the control volume.
+ */
+ const SolidState &solidState() const
+ { return solidState_; }
+
+ /*!
+ * \brief Returns 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 solidState_.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 = liquidPhaseIdx) const
+ { return fluidState_.saturation(phaseIdx); }
+
+ /*!
+ * \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 = liquidPhaseIdx) const
+ { return fluidState_.density(phaseIdx); }
+
+ /*!
+ * \brief Returns the effective pressure \f$\mathrm{[Pa]}\f$ of a given phase within
+ * the control volume.
+ *
+ * For the nonwetting phase (i.e. the gas phase), we assume
+ * infinite mobility, which implies that the nonwetting 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 = liquidPhaseIdx) const
+ { return fluidState_.pressure(phaseIdx); }
+
+ /*!
+ * \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 = liquidPhaseIdx) 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 nonwetting phase is infinitely mobile
+ */
+ Scalar viscosity(const int phaseIdx = liquidPhaseIdx) const
+ { return phaseIdx == liquidPhaseIdx ? fluidState_.viscosity(liquidPhaseIdx) : 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 = liquidPhaseIdx) const
+ { return phaseIdx == liquidPhaseIdx ? 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 nonwetting 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_, pressure(gasPhaseIdx) - pressure(liquidPhaseIdx));
+ }
+
+ /*!
+ * \brief Returns the pressureHead \f$\mathrm{[cm]}\f$ of a given phase within
+ * the control volume.
+ *
+ * For the nonwetting phase (i.e. the gas phase), we assume
+ * infinite mobility, which implies that the nonwetting 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 = liquidPhaseIdx) const
+ { return 100.0 *(pressure(phaseIdx) - pressure(gasPhaseIdx))/density(phaseIdx)/9.81; }
+
+ /*!
+ * \brief Returns the water content of a fluid phase within the finite volume.
+ *
+ * The water content is defined as the fraction of
+ * the saturation divided 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 = liquidPhaseIdx) const
+ { return saturation(phaseIdx) * solidState_.porosity(); }
+
+ /*!
+ * \brief Returns the mole fraction of a given component in a
+ * given phase within the control volume in \f$[-]\f$.
+ *
+ * \param phaseIdx The phase index
+ * \param compIdx The component index
+ */
+ Scalar moleFraction(const int phaseIdx, const int compIdx) const
+ {
+ if (compIdx != FluidSystem::comp0Idx)
+ DUNE_THROW(Dune::InvalidStateException, "There is only one component for Richards!");
+ return moleFraction_[phaseIdx];
+ }
+
+ /*!
+ * \brief Returns the mole fraction of a given component in a
+ * given phase within the control volume in \f$[-]\f$.
+ *
+ * \param phaseIdx The phase index
+ * \param compIdx The component index
+ */
+ Scalar massFraction(const int phaseIdx, const int compIdx) const
+ {
+ if (compIdx != FluidSystem::comp0Idx)
+ DUNE_THROW(Dune::InvalidStateException, "There is only one component for Richards!");
+ return massFraction_[phaseIdx];
+ }
+
+ /*!
+ * \brief Returns the mass density of a given phase within the
+ * control volume in \f$[mol/m^3]\f$.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar molarDensity(const int phaseIdx) const
+ {
+ return molarDensity_[phaseIdx];
+ }
+
+ /*!
+ * \brief Returns the binary diffusion coefficients for a phase in \f$[m^2/s]\f$.
+ */
+ Scalar diffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
+ {
+ assert(phaseIdx == gasPhaseIdx);
+ assert(compIIdx != compJIdx);
+ typename FluidSystem::ParameterCache paramCache;
+ paramCache.updatePhase(fluidState_, phaseIdx);
+ return FluidSystem::binaryDiffusionCoefficient(fluidState_, paramCache, phaseIdx, compIIdx, compJIdx);
+ }
+
+ /*!
+ * \brief Returns the effective diffusion coefficients for a phase in \f$[m^2/s]\f$.
+ */
+ Scalar effectiveDiffusionCoefficient(int phaseIdx, int compIIdx, int compJIdx) const
+ {
+ assert(phaseIdx == gasPhaseIdx);
+ assert(compIIdx != compJIdx);
+ return effectiveDiffCoeff_;
+ }
+private:
+ /*!
+ * \brief Fills the fluid state according to the primary variables.
+ *
+ * Taking the information from the primary variables,
+ * the fluid state is filled with every information that is
+ * necessary to evaluate the model's local residual.
+ *
+ * \param elemSol A vector containing all primary variables connected to the element
+ * \param problem The problem at hand.
+ * \param element The current element.
+ * \param scv The subcontrol volume.
+ * \param fluidState The fluid state to fill.
+ * \param solidState The solid state to fill.
+ */
+ template<class ElemSol, class Problem, class Element, class Scv>
+ void completeFluidState_(const ElemSol& elemSol,
+ const Problem& problem,
+ const Element& element,
+ const Scv& scv,
+ FluidState& fluidState,
+ SolidState& solidState)
+ {
+ EnergyVolVars::updateTemperature(elemSol, problem, element, scv, fluidState, solidState);
+
+ const auto fluidMatrixInteraction = problem.spatialParams().fluidMatrixInteraction(element, scv, elemSol);
+
+ const auto& priVars = elemSol[scv.localDofIndex()];
+
+ // set the wetting pressure
+ using std::max;
+ Scalar minPc = fluidMatrixInteraction.pc(1.0);
+ fluidState.setPressure(liquidPhaseIdx, priVars[Indices::pressureIdx]);
+ fluidState.setPressure(gasPhaseIdx, max(problem.nonwettingReferencePressure(), fluidState.pressure(liquidPhaseIdx) + minPc));
+
+ // compute the capillary pressure to compute the saturation
+ // make sure that we the capillary pressure is not smaller than the minimum pc
+ // this would possibly return unphysical values from regularized material laws
+ using std::max;
+ const Scalar pc = max(fluidMatrixInteraction.endPointPc(),
+ problem.nonwettingReferencePressure() - fluidState.pressure(liquidPhaseIdx));
+ const Scalar sw = fluidMatrixInteraction.sw(pc);
+ fluidState.setSaturation(liquidPhaseIdx, sw);
+ fluidState.setSaturation(gasPhaseIdx, 1.0-sw);
+
+ // density and viscosity
+ typename FluidSystem::ParameterCache paramCache;
+ paramCache.updateAll(fluidState);
+ fluidState.setDensity(liquidPhaseIdx,
+ FluidSystem::density(fluidState, paramCache, liquidPhaseIdx));
+ fluidState.setDensity(gasPhaseIdx,
+ FluidSystem::density(fluidState, paramCache, gasPhaseIdx));
+
+ fluidState.setViscosity(liquidPhaseIdx,
+ FluidSystem::viscosity(fluidState, paramCache, liquidPhaseIdx));
+
+ // compute and set the enthalpy
+ fluidState.setEnthalpy(liquidPhaseIdx, EnergyVolVars::enthalpy(fluidState, paramCache, liquidPhaseIdx));
+ fluidState.setEnthalpy(gasPhaseIdx, EnergyVolVars::enthalpy(fluidState, paramCache, gasPhaseIdx));
+ }
+
+protected:
+ FluidState fluidState_; //!< the fluid state
+ SolidState solidState_;
+ Scalar relativePermeabilityWetting_; //!< the relative permeability of the wetting phase
+ PermeabilityType permeability_; //!< the intrinsic permeability
+ Scalar minPc_; //!< the minimum capillary pressure (entry pressure)
+ Scalar moleFraction_[numPhases]; //!< The water mole fractions in water and air
+ Scalar massFraction_[numPhases]; //!< The water mass fractions in water and air
+ Scalar molarDensity_[numPhases]; //!< The molar density of water and air
+
+ // Effective diffusion coefficients for the phases
+ Scalar effectiveDiffCoeff_;
+
+};
+} // end namespace Dumux
+
+#endif
diff --git a/test/multidomain/embedded/1d3d/1p_richards/problem_soil.hh b/test/multidomain/embedded/1d3d/1p_richards/problem_soil.hh
index 25f00f7c81a19bc2e09b6ff6ee1d1ee6db3d5ae8..8f12a73a38856cce99ec7d1d07cd3794df4ec7a3 100644
--- a/test/multidomain/embedded/1d3d/1p_richards/problem_soil.hh
+++ b/test/multidomain/embedded/1d3d/1p_richards/problem_soil.hh
@@ -198,7 +198,6 @@ public:
PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{
PrimaryVariables priVars({initPressure_});
- priVars.setState(Indices::bothPhases);
return priVars;
}
diff --git a/test/multidomain/embedded/1d3d/1p_richards/properties.hh b/test/multidomain/embedded/1d3d/1p_richards/properties.hh
index 71f7e0c2622d9bdfcb70574a40f435e5fa98c0bf..928c06070b9ee4abb9af4588f03b55626e05e9e1 100644
--- a/test/multidomain/embedded/1d3d/1p_richards/properties.hh
+++ b/test/multidomain/embedded/1d3d/1p_richards/properties.hh
@@ -103,6 +103,12 @@ struct SpatialParams<TypeTag, TTag::Soil>
GetPropType<TypeTag, Properties::Scalar>>;
};
+// TODO: remove after release (3.6)
+// Set the primary variables type
+template<class TypeTag>
+struct PrimaryVariables<TypeTag, TTag::Soil>
+{ using type = Dune::FieldVector<GetPropType<TypeTag, Properties::Scalar>, GetPropType<TypeTag, Properties::ModelTraits>::numEq()>; };
+
// Create new type tags
namespace TTag {
struct Root { using InheritsFrom = std::tuple<OneP, CCTpfaModel>; };
@@ -149,6 +155,12 @@ struct SpatialParams<TypeTag, TTag::Root>
GetPropType<TypeTag, Properties::Scalar>>;
};
+// TODO: remove after release (3.6)
+// Set the primary variables type
+template<class TypeTag>
+struct PrimaryVariables<TypeTag, TTag::Root>
+{ using type = Dune::FieldVector<GetPropType<TypeTag, Properties::Scalar>, GetPropType<TypeTag, Properties::ModelTraits>::numEq()>; };
+
template<class TypeTag>
struct CouplingManager<TypeTag, TTag::SOILTYPETAG>
{
diff --git a/test/porousmediumflow/richards/analytical/problem.hh b/test/porousmediumflow/richards/analytical/problem.hh
index eb28d0b2d881c646a008a4fc05b9fecc4bdbced8..b486dd7d84b39d5e6a656041a7b00704fea7418a 100644
--- a/test/porousmediumflow/richards/analytical/problem.hh
+++ b/test/porousmediumflow/richards/analytical/problem.hh
@@ -184,7 +184,6 @@ public:
PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
{
PrimaryVariables values(0.0);
- values.setState(bothPhases);
const Scalar time = time_;
const Scalar pwTop = 98942.8;
const Scalar pwBottom = 95641.1;
@@ -226,7 +225,6 @@ public:
PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{
PrimaryVariables values(0.0);
- values.setState(bothPhases);
analyticalSolution(values, time_, globalPos);
return values;
}
diff --git a/test/porousmediumflow/richards/analytical/properties.hh b/test/porousmediumflow/richards/analytical/properties.hh
index 8f165a11551871a339504aefe3a689afc7936559..a244235daf1dfb5ac6cda47e3a40f8735ccd2c7f 100644
--- a/test/porousmediumflow/richards/analytical/properties.hh
+++ b/test/porousmediumflow/richards/analytical/properties.hh
@@ -80,6 +80,12 @@ struct FluidSystem<TypeTag, TTag::RichardsAnalytical>
using type = FluidSystems::TwoPImmiscible<Scalar, L, G>;
};
+// TODO: remove after release (3.6)
+// Set the primary variables type
+template<class TypeTag>
+struct PrimaryVariables<TypeTag, TTag::RichardsAnalytical>
+{ using type = Dune::FieldVector<GetPropType<TypeTag, Properties::Scalar>, GetPropType<TypeTag, Properties::ModelTraits>::numEq()>; };
+
} // end namespace Dumux::Properties
#endif
diff --git a/test/porousmediumflow/richards/annulus/problem.hh b/test/porousmediumflow/richards/annulus/problem.hh
index 2104e254109255dfbe58ed786731ae47c252cca5..c6492dbea960e5b781aa7b04021797db9e06b369 100644
--- a/test/porousmediumflow/richards/annulus/problem.hh
+++ b/test/porousmediumflow/richards/annulus/problem.hh
@@ -133,9 +133,7 @@ public:
return exactSolution(globalPos);
else
{
- PrimaryVariables values(initialPressure_);
- values.setState(Indices::bothPhases);
- return values;
+ return { initialPressure_ };
}
}
@@ -167,9 +165,7 @@ public:
- A*mu_/(2*M_PI*k_)*std::log(r/innerRadius_)
+ source_*0.25*(r*r - ri2)*mu_/k_;
- PrimaryVariables priVars(applyInverseKirchhoffTransformation_(psi));
- priVars.setState(Indices::bothPhases);
- return priVars;
+ return { applyInverseKirchhoffTransformation_(psi) };
}
PrimaryVariables exactSolution(const GlobalPosition& globalPos) const
diff --git a/test/porousmediumflow/richards/annulus/properties.hh b/test/porousmediumflow/richards/annulus/properties.hh
index 3d1fc01ac1281083aaddaf0453b5d40d1f28c08c..9b82614753d094467f793a8e02053490b1b64ca3 100644
--- a/test/porousmediumflow/richards/annulus/properties.hh
+++ b/test/porousmediumflow/richards/annulus/properties.hh
@@ -74,6 +74,12 @@ public:
using type = BoxFVGridGeometry<Scalar, GridView, enableCache, GGTraits>;
};
+// TODO: remove after release (3.6)
+// Set the primary variables type
+template<class TypeTag>
+struct PrimaryVariables<TypeTag, TTag::RichardsAnnulus>
+{ using type = Dune::FieldVector<GetPropType<TypeTag, Properties::Scalar>, GetPropType<TypeTag, Properties::ModelTraits>::numEq()>; };
+
} // end namespace Dumux::Properties
#endif
diff --git a/test/porousmediumflow/richards/benchmarks/problem.hh b/test/porousmediumflow/richards/benchmarks/problem.hh
index 8abc56cfb54aa71a0d6c225ff403db56d89020cf..05abd250f2b2659731a38433ca8e9bd3bcc60b32 100644
--- a/test/porousmediumflow/richards/benchmarks/problem.hh
+++ b/test/porousmediumflow/richards/benchmarks/problem.hh
@@ -135,7 +135,6 @@ public:
{
PrimaryVariables values(0.0);
values[Indices::pressureIdx] = initialPressure_;
- values.setState(Indices::bothPhases);
return values;
}
diff --git a/test/porousmediumflow/richards/lens/problem.hh b/test/porousmediumflow/richards/lens/problem.hh
index 7d56055fa14b6da26e650d1170fddc4f20231cd3..6837075a12d0926944484e63394d347b5710a083 100644
--- a/test/porousmediumflow/richards/lens/problem.hh
+++ b/test/porousmediumflow/richards/lens/problem.hh
@@ -184,7 +184,6 @@ private:
const Scalar sw = 0.0;
const Scalar pc = this->spatialParams().fluidMatrixInteractionAtPos(globalPos).pc(sw);
values[pressureIdx] = nonwettingReferencePressure() - pc;
- values.setState(bothPhases);
return values;
}
diff --git a/test/porousmediumflow/richards/lens/properties.hh b/test/porousmediumflow/richards/lens/properties.hh
index 8966d7e86399345e7d626384d6ea9b703e1a761f..4052eaaa97322d092c46f2259fe342960ed30f7b 100644
--- a/test/porousmediumflow/richards/lens/properties.hh
+++ b/test/porousmediumflow/richards/lens/properties.hh
@@ -78,6 +78,12 @@ struct SpatialParams<TypeTag, TTag::RichardsLens>
GetPropType<TypeTag, Properties::Scalar>>;
};
+// TODO: remove after release (3.6)
+// Set the primary variables type
+template<class TypeTag>
+struct PrimaryVariables<TypeTag, TTag::RichardsLens>
+{ using type = Dune::FieldVector<GetPropType<TypeTag, Properties::Scalar>, GetPropType<TypeTag, Properties::ModelTraits>::numEq()>; };
+
} // end namespace Dumux::Properties
#endif
diff --git a/test/porousmediumflow/richards/nonisothermal/conduction/problem.hh b/test/porousmediumflow/richards/nonisothermal/conduction/problem.hh
index 5ee04f58fad2dd87bdd1d3642bc202e97bab8603..6c850cf12c91e4ea81431f1731f1beb3b04d16df 100644
--- a/test/porousmediumflow/richards/nonisothermal/conduction/problem.hh
+++ b/test/porousmediumflow/richards/nonisothermal/conduction/problem.hh
@@ -255,7 +255,6 @@ private:
PrimaryVariables initial_(const GlobalPosition &globalPos) const
{
PrimaryVariables priVars(0.0);
- priVars.setState(liquidPhaseOnly);
priVars[pressureIdx] = 1e5; // initial condition for the pressure
priVars[temperatureIdx] = 290.;
diff --git a/test/porousmediumflow/richards/nonisothermal/conduction/properties.hh b/test/porousmediumflow/richards/nonisothermal/conduction/properties.hh
index 85f49de8f527b82c3e79161430a5833d9b86a05d..8ea69e39ed3eff902dec83bbc2b4ac104d257229 100644
--- a/test/porousmediumflow/richards/nonisothermal/conduction/properties.hh
+++ b/test/porousmediumflow/richards/nonisothermal/conduction/properties.hh
@@ -67,6 +67,12 @@ struct SpatialParams<TypeTag, TTag::RichardsNIConduction>
using type = RichardsNISpatialParams<GridGeometry, Scalar>;
};
+// TODO: remove after release (3.6)
+// Set the primary variables type
+template<class TypeTag>
+struct PrimaryVariables<TypeTag, TTag::RichardsNIConduction>
+{ using type = Dune::FieldVector<GetPropType<TypeTag, Properties::Scalar>, GetPropType<TypeTag, Properties::ModelTraits>::numEq()>; };
+
} // end namespace Dumux::Properties
#endif
diff --git a/test/porousmediumflow/richards/nonisothermal/convection/problem.hh b/test/porousmediumflow/richards/nonisothermal/convection/problem.hh
index fa0d78f7c38ad24f86827f71df0c14edb45ee338..1c78cac578427396e4e74768d0a8787a65da6185 100644
--- a/test/porousmediumflow/richards/nonisothermal/convection/problem.hh
+++ b/test/porousmediumflow/richards/nonisothermal/convection/problem.hh
@@ -279,7 +279,6 @@ private:
PrimaryVariables initial_(const GlobalPosition &globalPos) const
{
PrimaryVariables priVars(0.0);
- priVars.setState(liquidPhaseOnly);
priVars[pressureIdx] = pressureLow_; // initial condition for the pressure
priVars[temperatureIdx] = temperatureLow_;
return priVars;
diff --git a/test/porousmediumflow/richards/nonisothermal/convection/properties.hh b/test/porousmediumflow/richards/nonisothermal/convection/properties.hh
index 28512a4733751b634b957fed552ff9708969e7b2..5992161831acae176ba540e969aa30e9a9be95c3 100644
--- a/test/porousmediumflow/richards/nonisothermal/convection/properties.hh
+++ b/test/porousmediumflow/richards/nonisothermal/convection/properties.hh
@@ -68,6 +68,12 @@ struct SpatialParams<TypeTag, TTag::RichardsNIConvection>
using type = RichardsNISpatialParams<GridGeometry, Scalar>;
};
+// TODO: remove after release (3.6)
+// Set the primary variables type
+template<class TypeTag>
+struct PrimaryVariables<TypeTag, TTag::RichardsNIConvection>
+{ using type = Dune::FieldVector<GetPropType<TypeTag, Properties::Scalar>, GetPropType<TypeTag, Properties::ModelTraits>::numEq()>; };
+
} // end namespace Dumux::Properties
#endif
diff --git a/test/porousmediumflow/richards/nonisothermal/evaporation/properties.hh b/test/porousmediumflow/richards/nonisothermal/evaporation/properties.hh
index 1f1b341e169d2309ae74cde77e3f4abe002e5ea3..f7f34c5d1ee1f5986c111904fb07896f7507118a 100644
--- a/test/porousmediumflow/richards/nonisothermal/evaporation/properties.hh
+++ b/test/porousmediumflow/richards/nonisothermal/evaporation/properties.hh
@@ -32,7 +32,7 @@
#include <dumux/discretization/cctpfa.hh>
#include <dumux/discretization/box.hh>
-#include <dumux/porousmediumflow/richards/model.hh>
+#include <dumux/porousmediumflow/richardsextended/model.hh>
#include <dumux/material/fluidmatrixinteractions/2p/thermalconductivity/somerton.hh>
#include <dumux/material/fluidsystems/h2on2.hh>
@@ -43,7 +43,7 @@ namespace Dumux::Properties {
// Create new type tags
namespace TTag {
-struct RichardsNIEvaporation { using InheritsFrom = std::tuple<RichardsNI>; };
+struct RichardsNIEvaporation { using InheritsFrom = std::tuple<ExtendedRichardsNI>; };
struct RichardsNIEvaporationBox { using InheritsFrom = std::tuple<RichardsNIEvaporation, BoxModel>; };
struct RichardsNIEvaporationCC { using InheritsFrom = std::tuple<RichardsNIEvaporation, CCTpfaModel>; };
} // end namespace TTag
@@ -69,9 +69,6 @@ struct SpatialParams<TypeTag, TTag::RichardsNIEvaporation>
using type = RichardsNISpatialParams<GridGeometry, Scalar>;
};
-template<class TypeTag>
-struct EnableWaterDiffusionInAir<TypeTag, TTag::RichardsNIEvaporation> { static constexpr bool value = true; };
-
} // end namespace Dumux::Properties
#endif
diff --git a/test/references/test_md_embedded_1d3d_1p_richards_proj_3d-reference.vtu b/test/references/test_md_embedded_1d3d_1p_richards_proj_3d-reference.vtu
index b5ac1fc5e979d9c6eddaf3d3549a3fb0d22a6f16..aef45e62b62c9309a9b6dd3562dcee39d2e9229d 100644
--- a/test/references/test_md_embedded_1d3d_1p_richards_proj_3d-reference.vtu
+++ b/test/references/test_md_embedded_1d3d_1p_richards_proj_3d-reference.vtu
@@ -1613,167 +1613,6 @@
0.0497021 0.0522536 0.0518026 0.0532135 0.0529031 0.0529637 0.0528625 0.0529667 0.0533109 0.0532539 0.0529564 0.0528012
0.0530829 0.0537059 0.0538073 0.0536238 0.05325 0.0534762
</DataArray>
- <DataArray type="Float32" Name="phase presence" NumberOfComponents="1" format="ascii">
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3
- </DataArray>
</PointData>
<CellData Scalars="process rank">
<DataArray type="Float32" Name="process rank" NumberOfComponents="1" format="ascii">
diff --git a/test/references/test_md_embedded_1d3d_1p_richards_tpfabox_3d-reference.vtu b/test/references/test_md_embedded_1d3d_1p_richards_tpfabox_3d-reference.vtu
index 08e2d5f633b127867a7b1c1d5ca356e8e016c35a..5a5cc88de5e2ca9b71fba28a178e5fa22fca1458 100644
--- a/test/references/test_md_embedded_1d3d_1p_richards_tpfabox_3d-reference.vtu
+++ b/test/references/test_md_embedded_1d3d_1p_richards_tpfabox_3d-reference.vtu
@@ -1863,192 +1863,6 @@
0.061555 0.0615276 0.061451 0.0613407 0.0612206 0.0611194 0.061064 0.061071 0.0611381 0.0612426 0.0613508 0.0614303
0.0614594
</DataArray>
- <DataArray type="Float32" Name="phase presence" NumberOfComponents="1" format="ascii">
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3
- </DataArray>
</PointData>
<CellData Scalars="process rank">
<DataArray type="Float32" Name="process rank" NumberOfComponents="1" format="ascii">
diff --git a/test/references/test_md_embedded_1d3d_1p_richards_tpfatpfa_3d-reference.vtu b/test/references/test_md_embedded_1d3d_1p_richards_tpfatpfa_3d-reference.vtu
index 5972a86c7a8f8f3faae0bb9ba11ce15fc86d9047..2dfbdb227d94973d1de03f9e4ff6a3563de33789 100644
--- a/test/references/test_md_embedded_1d3d_1p_richards_tpfatpfa_3d-reference.vtu
+++ b/test/references/test_md_embedded_1d3d_1p_richards_tpfatpfa_3d-reference.vtu
@@ -1463,152 +1463,6 @@
0.0614909 0.0614254 0.0613047 0.0611489 0.0609875 0.0608593 0.0608076 0.0608684 0.0610141 0.0611769 0.0613093 0.0613825
0.0615633 0.0615065 0.0614028 0.0612712 0.0611389 0.0610388 0.0610033 0.0610504 0.0611616 0.0612921 0.0614016 0.0614632
</DataArray>
- <DataArray type="Float32" Name="phase presence" NumberOfComponents="1" format="ascii">
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- </DataArray>
<DataArray type="Float32" Name="process rank" NumberOfComponents="1" format="ascii">
0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0
diff --git a/test/references/test_richards_analytical_tpfa-reference.vtu b/test/references/test_richards_analytical_tpfa-reference.vtu
index 9bcbeef80d7814f386a1d1c929706bae80bf8a3b..e4e57bb0135914cdf94562eb58ac76de828ce291 100644
--- a/test/references/test_richards_analytical_tpfa-reference.vtu
+++ b/test/references/test_richards_analytical_tpfa-reference.vtu
@@ -453,51 +453,6 @@
0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4
0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4
</DataArray>
- <DataArray type="Float32" Name="phase presence" NumberOfComponents="1" format="ascii">
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3
- </DataArray>
<DataArray type="Float32" Name="process rank" NumberOfComponents="1" format="ascii">
0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0
diff --git a/test/references/test_richards_lens_box-reference.vtu b/test/references/test_richards_lens_box-reference.vtu
index 32ce87a8dc41fa210ba3a9e504da695ee4999199..9a399f43f95d8cb4a7ab3f30e1859a3b2805c522 100644
--- a/test/references/test_richards_lens_box-reference.vtu
+++ b/test/references/test_richards_lens_box-reference.vtu
@@ -421,44 +421,6 @@
0.285214 0.373147 0.374603 0.373147 0.285214 0.0325891 1.57507e-06 3.60822e-17 3.60822e-17 3.60822e-17 3.60822e-17 3.60822e-17
3.60822e-17 3.60822e-17 3.60822e-17 3.60822e-17 3.60822e-17
</DataArray>
- <DataArray type="Float32" Name="phase presence" NumberOfComponents="1" format="ascii">
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3
- </DataArray>
</PointData>
<CellData Scalars="process rank">
<DataArray type="Float32" Name="process rank" NumberOfComponents="1" format="ascii">
diff --git a/test/references/test_richards_lens_box_parallel-reference.vtu b/test/references/test_richards_lens_box_parallel-reference.vtu
index 901ca580caaa6296f63d6de186495e7f0dbf74d6..2eab6b679eb4215f810aae90af9f3544e20d3395 100644
--- a/test/references/test_richards_lens_box_parallel-reference.vtu
+++ b/test/references/test_richards_lens_box_parallel-reference.vtu
@@ -234,27 +234,6 @@
0.24609 0.254316 0.24609 0.0912612 3.60822e-17 3.60822e-17 3.60822e-17 3.60822e-17 3.60822e-17 3.60822e-17 1.57507e-06 0.0325891
0.285214 0.373147 0.374603 0.373147 0.285214
</DataArray>
- <DataArray type="Float32" Name="phase presence" NumberOfComponents="1" format="ascii">
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3
- </DataArray>
</PointData>
<CellData Scalars="process rank">
<DataArray type="Float32" Name="process rank" NumberOfComponents="1" format="ascii">
diff --git a/test/references/test_richards_lens_tpfa-reference.vtu b/test/references/test_richards_lens_tpfa-reference.vtu
index fa6c1e1b7ff7fb38be07c087ff9d65cc4fa2e1ba..5e90ef05ada4d552f7edb18ba3be40f7a082353d 100644
--- a/test/references/test_richards_lens_tpfa-reference.vtu
+++ b/test/references/test_richards_lens_tpfa-reference.vtu
@@ -412,40 +412,6 @@
0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0
</DataArray>
- <DataArray type="Float32" Name="phase presence" NumberOfComponents="1" format="ascii">
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- </DataArray>
</CellData>
<Points>
<DataArray type="Float32" Name="Coordinates" NumberOfComponents="3" format="ascii">
diff --git a/test/references/test_richards_lens_tpfa_parallel-reference.vtu b/test/references/test_richards_lens_tpfa_parallel-reference.vtu
index 7a8a5e3c2a909beb910f9e35789c96e1760b676c..130943b47fa1af5803dbc1ca38f59dbf3533bc7e 100644
--- a/test/references/test_richards_lens_tpfa_parallel-reference.vtu
+++ b/test/references/test_richards_lens_tpfa_parallel-reference.vtu
@@ -201,24 +201,6 @@
3.74335e-17 3.74335e-17 3.74335e-17 3.74335e-17 3.74335e-17 3.74335e-17 3.74335e-17 0.00133304 0.118767 0.135959 0.135959 0.118767
3.74335e-17 3.74335e-17 3.74335e-17 3.74335e-17 3.74335e-17 3.74335e-17 1.14143e-06 0.0316293 0.325895 0.338906 0.338906 0.325895
</DataArray>
- <DataArray type="Float32" Name="phase presence" NumberOfComponents="1" format="ascii">
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- 3 3 3 3 3 3 3 3 3 3 3 3
- </DataArray>
<DataArray type="Float32" Name="process rank" NumberOfComponents="1" format="ascii">
0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0
diff --git a/test/references/test_richardsni_conduction_box-reference.vtu b/test/references/test_richardsni_conduction_box-reference.vtu
index a822b334c890d16e693f22fb7270a20c62879430..ce677a27483b4cd4cb5742f16e696f9c9c450a77 100644
--- a/test/references/test_richardsni_conduction_box-reference.vtu
+++ b/test/references/test_richardsni_conduction_box-reference.vtu
@@ -327,42 +327,6 @@
0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4
0.4 0.4 0.4 0.4 0.4 0.4
</DataArray>
- <DataArray type="Float32" Name="phase presence" NumberOfComponents="1" format="ascii">
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1
- </DataArray>
<DataArray type="Float32" Name="T" NumberOfComponents="1" format="ascii">
300 299.327 300 299.327 298.659 298.659 298.003 298.003 297.362 297.362 296.742 296.742
296.147 296.147 295.58 295.58 295.043 295.043 294.539 294.539 294.068 294.068 293.632 293.632
diff --git a/test/references/test_richardsni_conduction_tpfa-reference.vtu b/test/references/test_richardsni_conduction_tpfa-reference.vtu
index 2b6f63a245a0b0f50a3eb2a4db6b9c85e954b47a..50566938138ba2d5fd654d97185031b5938f6ac4 100644
--- a/test/references/test_richardsni_conduction_tpfa-reference.vtu
+++ b/test/references/test_richardsni_conduction_tpfa-reference.vtu
@@ -212,25 +212,6 @@
0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4
0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4
</DataArray>
- <DataArray type="Float32" Name="phase presence" NumberOfComponents="1" format="ascii">
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1
- </DataArray>
<DataArray type="Float32" Name="velocity_liq (m/s)" NumberOfComponents="3" format="ascii">
-1.0497e-09 0 0 -1.0375e-09 0 0 -1.01365e-09 0 0 -9.81329e-10 0 0
-9.42115e-10 0 0 -8.95233e-10 0 0 -8.41905e-10 0 0 -7.83441e-10 0 0
diff --git a/test/references/test_richardsni_convection_box-reference.vtu b/test/references/test_richardsni_convection_box-reference.vtu
index b751ca7007ff5e256f56a9812f1ce83927096996..244d42b596c1f9f9f888874cc6eccc6a4657694c 100644
--- a/test/references/test_richardsni_convection_box-reference.vtu
+++ b/test/references/test_richardsni_convection_box-reference.vtu
@@ -179,22 +179,6 @@
0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4
0.4 0.4 0.4 0.4 0.4 0.4
</DataArray>
- <DataArray type="Float32" Name="phase presence" NumberOfComponents="1" format="ascii">
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1
- </DataArray>
<DataArray type="Float32" Name="temperatureExact" NumberOfComponents="1" format="ascii">
291 291 291 291 291 291 291 291 291 291 291 291
291 291 291 291 291 291 291 291 291 291 291 291
diff --git a/test/references/test_richardsni_convection_tpfa-reference.vtu b/test/references/test_richardsni_convection_tpfa-reference.vtu
index fdf46b11a0a8a1a6c54620b436fdf18b8d34a980..7494831bc616d0c856dd7764a7aa0b7cc7324910 100644
--- a/test/references/test_richardsni_convection_tpfa-reference.vtu
+++ b/test/references/test_richardsni_convection_tpfa-reference.vtu
@@ -102,15 +102,6 @@
0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4
0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4
</DataArray>
- <DataArray type="Float32" Name="phase presence" NumberOfComponents="1" format="ascii">
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1 1 1 1 1
- 1 1 1 1 1 1 1 1
- </DataArray>
<DataArray type="Float32" Name="process rank" NumberOfComponents="1" format="ascii">
0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0