From e398bd2f51a1196b126926959b333700ca771576 Mon Sep 17 00:00:00 2001
From: Mathis Kelm <mathis.kelm@iws.uni-stuttgart.de>
Date: Fri, 2 Sep 2022 14:15:12 +0200
Subject: [PATCH] [richards] Deprecate old extended Richards model

---
 dumux/common/deprecated.hh                    |  87 ++++++++
 dumux/porousmediumflow/richards/model.hh      |  10 +-
 .../richards/primaryvariableswitch.hh         | 137 -------------
 .../richards/volumevariables.hh               | 169 ++++++++--------
 .../embedded/1d3d/1p_richards/problem_soil.hh |   1 -
 .../embedded/1d3d/1p_richards/properties.hh   |  12 ++
 .../richards/analytical/problem.hh            |   2 -
 .../richards/analytical/properties.hh         |   6 +
 .../richards/annulus/problem.hh               |   8 +-
 .../richards/annulus/properties.hh            |   6 +
 .../richards/benchmarks/problem.hh            |   1 -
 .../porousmediumflow/richards/lens/problem.hh |   1 -
 .../richards/lens/properties.hh               |   6 +
 .../nonisothermal/conduction/problem.hh       |   1 -
 .../nonisothermal/conduction/properties.hh    |   6 +
 .../nonisothermal/convection/problem.hh       |   1 -
 .../nonisothermal/convection/properties.hh    |   6 +
 ...ded_1d3d_1p_richards_proj_3d-reference.vtu | 161 ---------------
 ..._1d3d_1p_richards_tpfabox_3d-reference.vtu | 186 ------------------
 ...1d3d_1p_richards_tpfatpfa_3d-reference.vtu | 146 --------------
 ...est_richards_analytical_tpfa-reference.vtu |  45 -----
 .../test_richards_lens_box-reference.vtu      |  38 ----
 ...t_richards_lens_box_parallel-reference.vtu |  21 --
 .../test_richards_lens_tpfa-reference.vtu     |  34 ----
 ..._richards_lens_tpfa_parallel-reference.vtu |  18 --
 ...st_richardsni_conduction_box-reference.vtu |  36 ----
 ...t_richardsni_conduction_tpfa-reference.vtu |  19 --
 ...st_richardsni_convection_box-reference.vtu |  16 --
 ...t_richardsni_convection_tpfa-reference.vtu |   9 -
 29 files changed, 227 insertions(+), 962 deletions(-)
 delete mode 100644 dumux/porousmediumflow/richards/primaryvariableswitch.hh

diff --git a/dumux/common/deprecated.hh b/dumux/common/deprecated.hh
index c8e75e6c2c..7e3c82c40f 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/richards/model.hh b/dumux/porousmediumflow/richards/model.hh
index 6da72cc383..93b53a661e 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/primaryvariableswitch.hh b/dumux/porousmediumflow/richards/primaryvariableswitch.hh
deleted file mode 100644
index 0c8ea23076..0000000000
--- a/dumux/porousmediumflow/richards/primaryvariableswitch.hh
+++ /dev/null
@@ -1,137 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- *   See the file COPYING for full copying permissions.                      *
- *                                                                           *
- *   This program is free software: you can redistribute it and/or modify    *
- *   it under the terms of the GNU General Public License as published by    *
- *   the Free Software Foundation, either version 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 RichardsModel
- * \brief The primary variable switch for the extended Richards model.
- */
-
-#ifndef DUMUX_RICHARDS_PRIMARY_VARIABLE_SWITCH_HH
-#define DUMUX_RICHARDS_PRIMARY_VARIABLE_SWITCH_HH
-
-#include <dumux/common/exceptions.hh>
-#include <dumux/common/parameters.hh>
-#include <dumux/material/constants.hh>
-#include <dumux/porousmediumflow/compositional/primaryvariableswitch.hh>
-
-namespace Dumux {
-
-/*!
- * \ingroup RichardsModel
- * \brief The primary variable switch controlling the phase presence state variable.
- */
-class ExtendedRichardsPrimaryVariableSwitch
-: public PrimaryVariableSwitch<ExtendedRichardsPrimaryVariableSwitch>
-{
-    using ParentType = PrimaryVariableSwitch<ExtendedRichardsPrimaryVariableSwitch>;
-    friend ParentType;
-
-public:
-    using ParentType::ParentType;
-
-protected:
-
-    // perform variable switch at a degree of freedom location
-    template<class VolumeVariables, class GlobalPosition>
-    bool update_(typename VolumeVariables::PrimaryVariables& priVars,
-                 const VolumeVariables& volVars,
-                 std::size_t dofIdxGlobal,
-                 const GlobalPosition& globalPos)
-    {
-        using Scalar = typename VolumeVariables::PrimaryVariables::value_type;
-        using Indices = typename VolumeVariables::Indices;
-        using FluidSystem = typename VolumeVariables::FluidSystem;
-
-        static const bool usePriVarSwitch = getParam<bool>("Problem.UsePrimaryVariableSwitch");
-        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
-        bool wouldSwitch = false;
-        int phasePresence = priVars.state();
-        int newPhasePresence = phasePresence;
-
-        // check if a primary var switch is necessary
-        if (phasePresence == Indices::gasPhaseOnly)
-        {
-            // if the mole fraction of water is larger than the one
-            // predicted by a liquid-vapor equilibrium
-            Scalar xnw = volVars.moleFraction(FluidSystem::gasPhaseIdx, liquidCompIdx);
-            Scalar xnwPredicted = FluidSystem::H2O::vaporPressure(volVars.temperature())
-                                  / volVars.pressure(FluidSystem::gasPhaseIdx);
-
-            Scalar xwMax = 1.0;
-            if (xnw / xnwPredicted > xwMax)
-                wouldSwitch = true;
-            if (this->wasSwitched_[dofIdxGlobal])
-                xwMax *= 1.01;
-
-            // if the ratio of predicted mole fraction to current mole fraction is larger than
-            // 100%, wetting phase appears
-            if (xnw / xnwPredicted > xwMax)
-            {
-                // wetting phase appears
-                if (this->verbosity() > 1)
-                    std::cout << "Liquid phase appears at dof " << dofIdxGlobal
-                              << ", coordinates: " << globalPos << ", xnw / xnwPredicted * 100: "
-                              << xnw / xnwPredicted * 100 << "%"
-                              << ", at x_n^w: " << priVars[Indices::switchIdx] << std::endl;
-                newPhasePresence = Indices::bothPhases;
-                priVars[Indices::switchIdx] = 0.0;
-            }
-        }
-        else if (phasePresence == Indices::bothPhases)
-        {
-            Scalar Smin = 0.0;
-            if (this->wasSwitched_[dofIdxGlobal])
-                Smin = -0.01;
-
-            if (volVars.saturation(FluidSystem::liquidPhaseIdx) <= Smin)
-            {
-                wouldSwitch = true;
-                // wetting phase disappears
-                newPhasePresence = Indices::gasPhaseOnly;
-                priVars[Indices::switchIdx] = volVars.moleFraction(FluidSystem::gasPhaseIdx, liquidCompIdx);
-
-                if (this->verbosity() > 1)
-                    std::cout << "Liquid phase disappears at dof " << dofIdxGlobal
-                              << ", coordinates: " << globalPos << ", sw: "
-                              << volVars.saturation(FluidSystem::liquidPhaseIdx)
-                              << ", x_n^w: " << priVars[Indices::switchIdx] << std::endl;
-            }
-        }
-        else if (phasePresence == Indices::liquidPhaseOnly)
-        {
-            DUNE_THROW(Dune::NotImplemented, "Water phase only phase presence!");
-        }
-
-        priVars.setState(newPhasePresence);
-        this->wasSwitched_[dofIdxGlobal] = wouldSwitch;
-        return phasePresence != newPhasePresence;
-    }
-};
-
-} // end namespace Dumux
-
-#endif
diff --git a/dumux/porousmediumflow/richards/volumevariables.hh b/dumux/porousmediumflow/richards/volumevariables.hh
index 076b5cba53..23a329164c 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/test/multidomain/embedded/1d3d/1p_richards/problem_soil.hh b/test/multidomain/embedded/1d3d/1p_richards/problem_soil.hh
index 25f00f7c81..8f12a73a38 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 71f7e0c262..928c06070b 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 eb28d0b2d8..b486dd7d84 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 8f165a1155..a244235daf 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 2104e25410..c6492dbea9 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 3d1fc01ac1..9b82614753 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 8abc56cfb5..05abd250f2 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 7d56055fa1..6837075a12 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 8966d7e863..4052eaaa97 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 5ee04f58fa..6c850cf12c 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 85f49de8f5..8ea69e39ed 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 fa0d78f7c3..1c78cac578 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 28512a4733..5992161831 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/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 b5ac1fc5e9..aef45e62b6 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 08e2d5f633..5a5cc88de5 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 5972a86c7a..2dfbdb227d 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 9bcbeef80d..e4e57bb013 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 32ce87a8dc..9a399f43f9 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 901ca580ca..2eab6b679e 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 fa6c1e1b7f..5e90ef05ad 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 7a8a5e3c2a..130943b47f 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 a822b334c8..ce677a2748 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 2b6f63a245..5056693813 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 b751ca7007..244d42b596 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 fdf46b11a0..7494831bc6 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
-- 
GitLab