From 9c05b51813d2f33802ab02ee44b8ae2beea6dd09 Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Tue, 31 Jan 2017 20:04:28 +0100
Subject: [PATCH] [richards] Use single phase fluidsystem

---
 .../richards/implicit/indices.hh              |   4 +-
 .../richards/implicit/localresidual.hh        | 133 ---------
 .../richards/implicit/model.hh                |  15 +-
 .../richards/implicit/propertydefaults.hh     |  65 +----
 .../richards/implicit/volumevariables.hh      | 100 +++----
 .../richards/implicit/CMakeLists.txt          |   4 +-
 .../implicit/richardslensspatialparams.hh     |   2 +-
 .../richardslensbox-reference-parallel.vtu    | 170 ++++-------
 test/references/richardslensbox-reference.vtu | 270 +++++------------
 test/references/richardslenscc-reference.vtu  | 274 +++++++-----------
 10 files changed, 305 insertions(+), 732 deletions(-)
 delete mode 100644 dumux/porousmediumflow/richards/implicit/localresidual.hh

diff --git a/dumux/porousmediumflow/richards/implicit/indices.hh b/dumux/porousmediumflow/richards/implicit/indices.hh
index 32be851c89..0216062c0c 100644
--- a/dumux/porousmediumflow/richards/implicit/indices.hh
+++ b/dumux/porousmediumflow/richards/implicit/indices.hh
@@ -55,8 +55,8 @@ struct RichardsIndices
     //////////
     // phase indices
     //////////
-    static const int wPhaseIdx = FluidSystem::wPhaseIdx; //!< Index of the wetting phase;
-    static const int nPhaseIdx = FluidSystem::nPhaseIdx; //!< Index of the non-wetting phase;
+    static const int wPhaseIdx = 0; //!< Index of the wetting phase;
+    static const int nPhaseIdx = 1; //!< Index of the non-wetting phase;
 };
 // \}
 
diff --git a/dumux/porousmediumflow/richards/implicit/localresidual.hh b/dumux/porousmediumflow/richards/implicit/localresidual.hh
deleted file mode 100644
index 5659a08f07..0000000000
--- a/dumux/porousmediumflow/richards/implicit/localresidual.hh
+++ /dev/null
@@ -1,133 +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 2 of the License, or       *
- *   (at your option) any later version.                                     *
- *                                                                           *
- *   This program is distributed in the hope that it will be useful,         *
- *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
- *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
- *   GNU General Public License for more details.                            *
- *                                                                           *
- *   You should have received a copy of the GNU General Public License       *
- *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
- *****************************************************************************/
-/*!
- * \file
- *
- * \brief Element-wise calculation of the Jacobian matrix for problems
- *        using the n-phase immiscible fully implicit models.
- */
-#ifndef DUMUX_RICHARDS_LOCAL_RESIDUAL_HH
-#define DUMUX_RICHARDS_LOCAL_RESIDUAL_HH
-
-namespace Dumux
-{
-/*!
- * \ingroup OnePModel
- * \ingroup RichardsLocalResidual
- * \brief Element-wise calculation of the Jacobian matrix for problems
- *        using the n-phase immiscible fully implicit models.
- */
-template<class TypeTag>
-class RichardsLocalResidual : public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
-{
-    typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation;
-
-    using ParentType = typename GET_PROP_TYPE(TypeTag, BaseLocalResidual);
-    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
-    using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
-    using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
-    using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
-    using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
-    using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
-    using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
-    using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
-    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
-    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
-    using Element = typename GridView::template Codim<0>::Entity;
-    using EnergyLocalResidual = typename GET_PROP_TYPE(TypeTag, EnergyLocalResidual);
-    using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
-    // first index for the mass balance
-    enum { conti0EqIdx = Indices::conti0EqIdx };
-
-    enum {
-        nPhaseIdx = Indices::nPhaseIdx,
-        wPhaseIdx = Indices::wPhaseIdx
-    };
-
-public:
-
-    /*!
-     * \brief Evaluate the rate of change of all conservation
-     *        quantites (e.g. phase mass) within a sub-control
-     *        volume of a finite volume element for the immiscible models.
-     * \param scv The sub control volume
-     * \param volVars The current or previous volVars
-     * \note This function should not include the source and sink terms.
-     * \note The volVars can be different to allow computing
-     *       the implicit euler time derivative here
-     */
-    PrimaryVariables computeStorage(const SubControlVolume& scv,
-                                    const VolumeVariables& volVars) const
-    {
-        // partial time derivative of the phase mass
-        PrimaryVariables storage;
-        storage[conti0EqIdx] = volVars.porosity()
-                               * volVars.density(wPhaseIdx)
-                               * volVars.saturation(wPhaseIdx);
-
-        //! The energy storage in the fluid phase with index phaseIdx
-        EnergyLocalResidual::fluidPhaseStorage(storage, scv, volVars, wPhaseIdx);
-
-        //! The energy storage in the solid matrix
-        EnergyLocalResidual::solidPhaseStorage(storage, scv, volVars);
-
-        return storage;
-    }
-
-
-    /*!
-     * \brief Evaluate the mass flux over a face of a sub control volume
-     * \param scvf The sub control volume face to compute the flux on
-     */
-    PrimaryVariables computeFlux(const Element& element,
-                                 const FVElementGeometry& fvGeometry,
-                                 const ElementVolumeVariables& elemVolVars,
-                                 const SubControlVolumeFace& scvf,
-                                 const ElementFluxVariablesCache& elemFluxVarsCache) const
-    {
-        FluxVariables fluxVars;
-        fluxVars.init(this->problem(), element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
-
-        PrimaryVariables flux;
-        // the physical quantities for which we perform upwinding
-        auto upwindTerm = [](const VolumeVariables& volVars)
-                          { return volVars.density(wPhaseIdx)*volVars.mobility(wPhaseIdx); };
-
-        flux[conti0EqIdx] = fluxVars.advectiveFlux(wPhaseIdx, upwindTerm);
-
-        //! Add advective phase energy fluxes. For isothermal model the contribution is zero.
-        EnergyLocalResidual::heatConvectionFlux(flux, fluxVars, wPhaseIdx);
-
-        //! Add diffusive energy fluxes. For isothermal model the contribution is zero.
-        EnergyLocalResidual::heatConductionFlux(flux, fluxVars);
-
-        return flux;
-    }
-
-private:
-    Implementation *asImp_()
-    { return static_cast<Implementation *> (this); }
-
-    const Implementation *asImp_() const
-    { return static_cast<const Implementation *> (this); }
-};
-
-} // end namespace Dumux
-
-#endif
diff --git a/dumux/porousmediumflow/richards/implicit/model.hh b/dumux/porousmediumflow/richards/implicit/model.hh
index b27bc52cfc..d71307fc9b 100644
--- a/dumux/porousmediumflow/richards/implicit/model.hh
+++ b/dumux/porousmediumflow/richards/implicit/model.hh
@@ -26,7 +26,7 @@
 #define DUMUX_RICHARDS_MODEL_HH
 
 #include <dumux/implicit/model.hh>
-#include <dumux/porousmediumflow/implicit/velocityoutput.hh>
+#include <dumux/porousmediumflow/nonisothermal/implicit/model.hh>
 #include "properties.hh"
 
 namespace Dumux
@@ -108,6 +108,9 @@ class RichardsModel : public GET_PROP_TYPE(TypeTag, BaseModel)
     using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
     using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
     using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+
+    using NonIsothermalModel = Dumux::NonIsothermalModel<TypeTag>;
+
     enum {
         nPhaseIdx = Indices::nPhaseIdx,
         wPhaseIdx = Indices::wPhaseIdx
@@ -132,15 +135,15 @@ public:
         vtkOutputModule.addSecondaryVariable("pw", [](const VolumeVariables& v){ return v.pressure(wPhaseIdx); });
         vtkOutputModule.addSecondaryVariable("pn", [](const VolumeVariables& v){ return v.pressure(nPhaseIdx); });
         vtkOutputModule.addSecondaryVariable("pc", [](const VolumeVariables& v){ return v.capillaryPressure(); });
-        vtkOutputModule.addSecondaryVariable("rhoW", [](const VolumeVariables& v){ return v.density(wPhaseIdx); });
-        vtkOutputModule.addSecondaryVariable("rhoN", [](const VolumeVariables& v){ return v.density(nPhaseIdx); });
-        vtkOutputModule.addSecondaryVariable("mobW", [](const VolumeVariables& v){ return v.mobility(wPhaseIdx); });
-        vtkOutputModule.addSecondaryVariable("mobN", [](const VolumeVariables& v){ return v.mobility(nPhaseIdx); });
-        vtkOutputModule.addSecondaryVariable("temperature", [](const VolumeVariables& v){ return v.temperature(); });
+        vtkOutputModule.addSecondaryVariable("density", [](const VolumeVariables& v){ return v.density(wPhaseIdx); });
+        vtkOutputModule.addSecondaryVariable("mobility", [](const VolumeVariables& v){ return v.mobility(wPhaseIdx); });
+        vtkOutputModule.addSecondaryVariable("kr", [](const VolumeVariables& v){ return v.relativePermeability(wPhaseIdx); });
         vtkOutputModule.addSecondaryVariable("porosity", [](const VolumeVariables& v){ return v.porosity(); });
         if(GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity));
             vtkOutputModule.addSecondaryVariable("pressure head", [](const VolumeVariables& v){ return v.pressureHead(wPhaseIdx); });
         vtkOutputModule.addSecondaryVariable("water content", [](const VolumeVariables& v){ return v.waterContent(wPhaseIdx); });
+
+        NonIsothermalModel::maybeAddTemperature(vtkOutputModule);
     }
 };
 
diff --git a/dumux/porousmediumflow/richards/implicit/propertydefaults.hh b/dumux/porousmediumflow/richards/implicit/propertydefaults.hh
index c1948c41c2..859ca52789 100644
--- a/dumux/porousmediumflow/richards/implicit/propertydefaults.hh
+++ b/dumux/porousmediumflow/richards/implicit/propertydefaults.hh
@@ -32,14 +32,15 @@
 #include "indices.hh"
 #include "volumevariables.hh"
 #include "properties.hh"
-#include "localresidual.hh"
 #include "newtoncontroller.hh"
 
+#include <dumux/porousmediumflow/immiscible/localresidual.hh>
 #include <dumux/porousmediumflow/nonisothermal/implicit/propertydefaults.hh>
-#include <dumux/material/components/nullcomponent.hh>
 #include <dumux/material/fluidmatrixinteractions/2p/thermalconductivitysomerton.hh>
-#include <dumux/material/fluidsystems/2pimmiscible.hh>
 #include <dumux/material/spatialparams/implicit.hh>
+#include <dumux/material/components/simpleh2o.hh>
+#include <dumux/material/fluidsystems/liquidphase.hh>
+#include <dumux/material/fluidstates/immiscible.hh>
 
 namespace Dumux
 {
@@ -51,11 +52,15 @@ namespace Properties {
 //////////////////////////////////////////////////////////////////
 //! Number of equations required by the model
 SET_INT_PROP(Richards, NumEq, 1);
+
 //! Number of fluid phases considered
-SET_INT_PROP(Richards, NumPhases, 2);
+//! Although the number of phases is two for the Richards model the
+//! non-wetting phase is not balanced and thus we only have one
+//! phase with a balance equation
+SET_INT_PROP(Richards, NumPhases, 1);
 
 //! The local residual operator
-SET_TYPE_PROP(Richards, LocalResidual, RichardsLocalResidual<TypeTag>);
+SET_TYPE_PROP(Richards, LocalResidual, ImmiscibleLocalResidual<TypeTag>);
 
 //! The global model used
 SET_TYPE_PROP(Richards, Model, RichardsModel<TypeTag>);
@@ -89,59 +94,15 @@ SET_TYPE_PROP(Richards, SpatialParams, ImplicitSpatialParams<TypeTag>);
  */
 SET_TYPE_PROP(Richards, MaterialLawParams, typename GET_PROP_TYPE(TypeTag, MaterialLaw)::Params);
 
-/*!
- * \brief The wetting phase used.
- *
- * By default we use the null-phase, i.e. this has to be defined by
- * the problem for the program to work. Please be aware that you
- * should be careful to use the Richards model in conjunction with
- * liquid non-wetting phases. This is only meaningful if the viscosity
- * of the liquid phase is _much_ lower than the viscosity of the
- * wetting phase.
- */
-
-SET_PROP(Richards, WettingPhase)
-{
-private:
-    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
-public:
-    using type = FluidSystems::LiquidPhase<Scalar, NullComponent<Scalar>>;
-};
-
-/*!
- * \brief The wetting phase used.
- *
- * By default we use the null-phase, i.e. this has to be defined by
- * the problem for the program to work. This doed not need to be
- * specified by the problem for the Richards model to work because the
- * Richards model does not conserve the non-wetting phase.
- */
-SET_PROP(Richards, NonwettingPhase)
-{
-private:
-    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
-public:
-    using type = FluidSystems::GasPhase<Scalar, NullComponent<Scalar>>;
-};
-
 /*!
  *\brief The fluid system used by the model.
  *
- * By default this uses the immiscible twophase fluid system. The
- * actual fluids used are specified using in the problem definition by
- * the WettingPhase and NonwettingPhase properties. Be aware that
- * using different fluid systems in conjunction with the Richards
- * model only makes very limited sense.
+ * By default this uses the liquid phase fluid system with simple H2O.
  */
 SET_PROP(Richards, FluidSystem)
 {
-private:
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
-    using WettingPhase = typename GET_PROP_TYPE(TypeTag, WettingPhase);
-    using NonwettingPhase = typename GET_PROP_TYPE(TypeTag, NonwettingPhase);
-
-public:
-    using type = FluidSystems::TwoPImmiscible<Scalar, WettingPhase, NonwettingPhase>;
+    using type = FluidSystems::LiquidPhase<Scalar, SimpleH2O<Scalar>>;
 };
 
 /*!
@@ -190,7 +151,7 @@ SET_TYPE_PROP(RichardsNI, IsothermalModel, RichardsModel<TypeTag>);
 SET_TYPE_PROP(RichardsNI, IsothermalVolumeVariables, RichardsVolumeVariables<TypeTag>);
 
 //set isothermal LocalResidual
-SET_TYPE_PROP(RichardsNI, IsothermalLocalResidual, RichardsLocalResidual<TypeTag>);
+SET_TYPE_PROP(RichardsNI, IsothermalLocalResidual, ImmiscibleLocalResidual<TypeTag>);
 
 //set isothermal Indices
 SET_TYPE_PROP(RichardsNI, IsothermalIndices, RichardsIndices<TypeTag>);
diff --git a/dumux/porousmediumflow/richards/implicit/volumevariables.hh b/dumux/porousmediumflow/richards/implicit/volumevariables.hh
index 2c5eca5fea..ddf3d0372e 100644
--- a/dumux/porousmediumflow/richards/implicit/volumevariables.hh
+++ b/dumux/porousmediumflow/richards/implicit/volumevariables.hh
@@ -26,7 +26,6 @@
 
 #include "properties.hh"
 
-#include <dumux/material/fluidstates/immiscible.hh>
 #include <dumux/discretization/volumevariables.hh>
 
 namespace Dumux
@@ -79,8 +78,6 @@ public:
                 const Element &element,
                 const SubControlVolume& scv)
     {
-        assert(!FluidSystem::isLiquid(nPhaseIdx));
-
         ParentType::update(elemSol, problem, element, scv);
 
         completeFluidState(elemSol, problem, element, scv, fluidState_);
@@ -93,9 +90,6 @@ public:
         pnRef_ = problem.nonWettingReferencePressure();
         porosity_ = problem.spatialParams().porosity(element, scv, elemSol);
         permeability_ = problem.spatialParams().permeability(element, scv, elemSol);
-
-        // energy related quantities not belonging to the fluid state
-        asImp_().updateEnergy_(elemSol, problem, element, scv);
     }
 
     /*!
@@ -107,7 +101,7 @@ public:
                                    const SubControlVolume& scv,
                                    FluidState& fluidState)
     {
-        Scalar t = Implementation::temperature_(elemSol, problem, element, scv);
+        Scalar t = ParentType::temperature(elemSol, problem, element, scv);
         fluidState.setTemperature(t);
 
         const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
@@ -115,26 +109,23 @@ public:
 
         const Scalar minPc = MaterialLaw::pc(materialParams, 1.0);
         fluidState.setPressure(wPhaseIdx, priVars[pressureIdx]);
-        fluidState.setPressure(nPhaseIdx, std::max(problem.nonWettingReferencePressure(), priVars[pressureIdx] + minPc));
+
+        using std::max;
+        pn_ = max(problem.nonWettingReferencePressure(), priVars[pressureIdx] + minPc);
 
         // saturations
-        const Scalar sw = MaterialLaw::sw(materialParams, fluidState.pressure(nPhaseIdx) - fluidState.pressure(wPhaseIdx));
+        const Scalar sw = MaterialLaw::sw(materialParams, pn_ - fluidState.pressure(wPhaseIdx));
         fluidState.setSaturation(wPhaseIdx, sw);
-        fluidState.setSaturation(nPhaseIdx, 1 - sw);
 
         // density and viscosity
         typename FluidSystem::ParameterCache paramCache;
         paramCache.updateAll(fluidState);
         fluidState.setDensity(wPhaseIdx, FluidSystem::density(fluidState, paramCache, wPhaseIdx));
-        fluidState.setDensity(nPhaseIdx, 1e-10);
 
         fluidState.setViscosity(wPhaseIdx, FluidSystem::viscosity(fluidState, paramCache, wPhaseIdx));
-        fluidState.setViscosity(nPhaseIdx, 1e-10);
 
         // compute and set the enthalpy
-        fluidState.setEnthalpy(wPhaseIdx, Implementation::enthalpy_(fluidState, paramCache, wPhaseIdx));
-        fluidState.setEnthalpy(nPhaseIdx, Implementation::enthalpy_(fluidState, paramCache, nPhaseIdx));
-
+        fluidState.setEnthalpy(wPhaseIdx, Implementation::enthalpy(fluidState, paramCache, wPhaseIdx));
     }
 
     /*!
@@ -144,6 +135,12 @@ public:
     const FluidState &fluidState() const
     { return fluidState_; }
 
+    /*!
+     * \brief Return the temperature
+     */
+    Scalar temperature() const
+    { return fluidState_.temperature(); }
+
     /*!
      * \brief Returns the average porosity [] within the control volume.
      *
@@ -170,7 +167,7 @@ public:
      * \param phaseIdx The index of the fluid phase
      */
     Scalar saturation(const int phaseIdx) const
-    { return fluidState_.saturation(phaseIdx); }
+    { return phaseIdx == wPhaseIdx ? fluidState_.saturation(wPhaseIdx) : 1.0-fluidState_.saturation(wPhaseIdx); }
 
     /*!
      * \brief Returns the average mass density \f$\mathrm{[kg/m^3]}\f$ of a given
@@ -179,7 +176,7 @@ public:
      * \param phaseIdx The index of the fluid phase
      */
     Scalar density(const int phaseIdx) const
-    { return fluidState_.density(phaseIdx); }
+    { return phaseIdx == wPhaseIdx ? fluidState_.density(phaseIdx) : 0.0; }
 
     /*!
      * \brief Returns the effective pressure \f$\mathrm{[Pa]}\f$ of a given phase within
@@ -193,17 +190,7 @@ public:
      * \param phaseIdx The index of the fluid phase
      */
     Scalar pressure(const int phaseIdx) const
-    { return fluidState_.pressure(phaseIdx); }
-
-    /*!
-     * \brief Returns average temperature \f$\mathrm{[K]}\f$ inside the control volume.
-     *
-     * Note that we assume thermodynamic equilibrium, i.e. the
-     * temperature of the rock matrix and of all fluid phases are
-     * identical.
-     */
-    Scalar temperature() const
-    { return fluidState_.temperature(); }
+    { return phaseIdx == wPhaseIdx ? fluidState_.pressure(phaseIdx) : pn_; }
 
     /*!
      * \brief Returns the effective mobility \f$\mathrm{[1/(Pa*s)]}\f$ of a given phase within
@@ -219,6 +206,16 @@ public:
     Scalar mobility(const int phaseIdx) const
     { return relativePermeability(phaseIdx)/fluidState_.viscosity(phaseIdx); }
 
+    /*!
+     * \brief Returns the dynamic viscosity \f$\mathrm{[Pa*s]}\f$ of a given phase within
+     *        the control volume.
+     *
+     * \param phaseIdx The index of the fluid phase
+     * \note The non-wetting phase is infinitely mobile
+     */
+    Scalar viscosity(const int phaseIdx) const
+    { return phaseIdx == wPhaseIdx ? fluidState_.viscosity(wPhaseIdx) : 0.0; }
+
     /*!
      * \brief Returns relative permeability [-] of a given phase within
      *        the control volume.
@@ -226,11 +223,7 @@ public:
      * \param phaseIdx The index of the fluid phase
      */
     Scalar relativePermeability(const int phaseIdx) const
-    {
-        if (phaseIdx == wPhaseIdx)
-            return relativePermeabilityWetting_;
-        return 1;
-    }
+    { return phaseIdx == wPhaseIdx ? relativePermeabilityWetting_ : 1.0; }
 
     /*!
      * \brief Returns the effective capillary pressure \f$\mathrm{[Pa]}\f$ within the
@@ -241,9 +234,7 @@ public:
      * \f[ p_c = p_n - p_w \f]
      */
     Scalar capillaryPressure() const
-    {
-        return fluidState_.pressure(nPhaseIdx) - fluidState_.pressure(wPhaseIdx);
-    }
+    { return pn_ - fluidState_.pressure(wPhaseIdx); }
 
     /*!
      * \brief Returns the pressureHead \f$\mathrm{[cm]}\f$ of a given phase within
@@ -260,9 +251,7 @@ public:
      *       or the gravity different
      */
     Scalar pressureHead(const int phaseIdx) const
-    {
-        return 100.0 *(fluidState_.pressure(phaseIdx) - pnRef_)/fluidState_.density(phaseIdx)/9.81;
-    }
+    { return 100.0 *(pressure(phaseIdx) - pnRef_)/density(phaseIdx)/9.81; }
 
     /*!
      * \brief Returns the water content
@@ -276,41 +265,15 @@ public:
      *       manually do a conversion.
      */
     Scalar waterContent(const int phaseIdx) const
-    {
-        return fluidState_.saturation(phaseIdx)* porosity_;
-    }
+    { return saturation(phaseIdx) * porosity_; }
 
 protected:
-    static Scalar temperature_(const ElementSolutionVector &elemSol,
-                               const Problem& problem,
-                               const Element &element,
-                               const SubControlVolume &scv)
-    {
-        return problem.temperatureAtPos(scv.dofPosition());
-    }
-
-    template<class ParameterCache>
-    static Scalar enthalpy_(const FluidState& fluidState,
-                            const ParameterCache& paramCache,
-                            const int phaseIdx)
-    {
-        return 0;
-    }
-
-    /*!
-     * \brief Called by update() to compute the energy related quantities.
-     */
-    void updateEnergy_(const ElementSolutionVector &elemSol,
-                       const Problem &problem,
-                       const Element &element,
-                       const SubControlVolume& scv)
-    {}
-
     FluidState fluidState_;
     Scalar relativePermeabilityWetting_;
     Scalar porosity_;
     PermeabilityType permeability_;
     Scalar pnRef_;
+    static Scalar pn_;
 
 private:
     Implementation &asImp_()
@@ -320,6 +283,9 @@ private:
     { return *static_cast<const Implementation*>(this); }
 };
 
+template <class TypeTag>
+typename GET_PROP_TYPE(TypeTag, Scalar) RichardsVolumeVariables<TypeTag>::pn_;
+
 } // end namespace Dumux
 
 #endif
diff --git a/test/porousmediumflow/richards/implicit/CMakeLists.txt b/test/porousmediumflow/richards/implicit/CMakeLists.txt
index 96adf2a404..0726f6c2a8 100644
--- a/test/porousmediumflow/richards/implicit/CMakeLists.txt
+++ b/test/porousmediumflow/richards/implicit/CMakeLists.txt
@@ -5,7 +5,7 @@ if(MPI_FOUND)
                  python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
                    --script fuzzy
                    --files ${CMAKE_SOURCE_DIR}/test/references/richardslensbox-reference-parallel.vtu
-                           ${CMAKE_CURRENT_BINARY_DIR}/s0002-p0000-richardslensbox-00007.vtu
+                           ${CMAKE_CURRENT_BINARY_DIR}/s0002-p0000-richardslensbox-00008.vtu
                    --command "${MPIEXEC} -np 2 ${CMAKE_CURRENT_BINARY_DIR}/test_boxrichards")
 else()
 # isothermal tests
@@ -13,7 +13,7 @@ else()
                  python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
                    --script fuzzy
                    --files ${CMAKE_SOURCE_DIR}/test/references/richardslensbox-reference.vtu
-                           ${CMAKE_CURRENT_BINARY_DIR}/richardslensbox-00007.vtu
+                           ${CMAKE_CURRENT_BINARY_DIR}/richardslensbox-00008.vtu
                    --command "${CMAKE_CURRENT_BINARY_DIR}/test_boxrichards")
 endif()
 
diff --git a/test/porousmediumflow/richards/implicit/richardslensspatialparams.hh b/test/porousmediumflow/richards/implicit/richardslensspatialparams.hh
index 6c589ceac4..19a6914b5a 100644
--- a/test/porousmediumflow/richards/implicit/richardslensspatialparams.hh
+++ b/test/porousmediumflow/richards/implicit/richardslensspatialparams.hh
@@ -95,7 +95,7 @@ public:
     {
 
         lensLowerLeft_ = {1.0, 2.0};
-        lensUpperRight_ = {4.0, 3.0};
+        lensUpperRight_ = {4.1, 3.1};
 
         // residual saturations
         lensMaterialParams_.setSwr(0.18);
diff --git a/test/references/richardslensbox-reference-parallel.vtu b/test/references/richardslensbox-reference-parallel.vtu
index acca99ab72..d64d207eca 100644
--- a/test/references/richardslensbox-reference-parallel.vtu
+++ b/test/references/richardslensbox-reference-parallel.vtu
@@ -2,28 +2,7 @@
 <VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
   <UnstructuredGrid>
     <Piece NumberOfCells="192" NumberOfPoints="221">
-      <PointData Scalars="Sn">
-        <DataArray type="Float32" Name="Sn" 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 0.999972 0.999959
-          0.999972 1 1 1 1 1 1 1 1 1 0.997392 0.909005
-          0.904424 0.909005 0.997392 1 1 1 1 1 1 1 0.998941 0.771659
-          0.385726 0.365066 0.385726 0.771659 1 1 1 1 1 1 0.999996 0.918448
-          0.287481 0.0673969 0.0636792 0.0673969 0.287481
-        </DataArray>
+      <PointData Scalars="Sw">
         <DataArray type="Float32" Name="Sw" NumberOfComponents="1" format="ascii">
           9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17
           9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17
@@ -45,26 +24,26 @@
           0.614274 0.634934 0.614274 0.228341 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 4.31542e-06 0.0815522
           0.712519 0.932603 0.936321 0.932603 0.712519
         </DataArray>
-        <DataArray type="Float32" Name="pn" NumberOfComponents="1" format="ascii">
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000
+        <DataArray type="Float32" Name="Sn" 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 0.999972 0.999959
+          0.999972 1 1 1 1 1 1 1 1 1 0.997392 0.909005
+          0.904424 0.909005 0.997392 1 1 1 1 1 1 1 0.998941 0.771659
+          0.385726 0.365066 0.385726 0.771659 1 1 1 1 1 1 0.999996 0.918448
+          0.287481 0.0673969 0.0636792 0.0673969 0.287481
         </DataArray>
         <DataArray type="Float32" Name="pw" NumberOfComponents="1" format="ascii">
           97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4
@@ -87,6 +66,27 @@
           99733.4 99738.8 99733.4 99586.6 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.5 99323.6
           99759.2 99835.1 99837.3 99835.1 99759.2
         </DataArray>
+        <DataArray type="Float32" Name="pn" NumberOfComponents="1" format="ascii">
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000
+        </DataArray>
         <DataArray type="Float32" Name="pc" NumberOfComponents="1" format="ascii">
           2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62
           2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62
@@ -108,7 +108,7 @@
           266.621 261.187 266.621 413.416 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.51 676.417
           240.748 164.873 162.729 164.873 240.748
         </DataArray>
-        <DataArray type="Float32" Name="rhoW" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="density" NumberOfComponents="1" format="ascii">
           1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
           1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
           1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
@@ -129,28 +129,7 @@
           1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
           1000 1000 1000 1000 1000
         </DataArray>
-        <DataArray type="Float32" Name="rhoN" NumberOfComponents="1" format="ascii">
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10
-        </DataArray>
-        <DataArray type="Float32" Name="mobW" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="mobility" 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
           0 0 0 0 0 0 0 0 0 0 0 0
@@ -171,26 +150,26 @@
           145.944 164.193 145.944 3.93436 0 0 0 0 0 0 0 0.0198338
           248.432 732.351 752.186 732.351 248.432
         </DataArray>
-        <DataArray type="Float32" Name="mobN" NumberOfComponents="1" format="ascii">
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10
+        <DataArray type="Float32" Name="kr" 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
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 1.57172e-05
+          2.3053e-05 1.57172e-05 0 0 0 0 0 0 0 0 0 0.00287218
+          0.120485 0.135536 0.120485 0.00287218 0 0 0 0 0 0 0 1.05648e-05
+          0.237214 0.71697 0.737982 0.71697 0.237214
         </DataArray>
         <DataArray type="Float32" Name="porosity" NumberOfComponents="1" format="ascii">
           0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4
@@ -213,27 +192,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
         </DataArray>
-        <DataArray type="Float32" Name="temperature" NumberOfComponents="1" format="ascii">
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15
-        </DataArray>
         <DataArray type="Float32" Name="pressure head" NumberOfComponents="1" format="ascii">
           -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861
           -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861
diff --git a/test/references/richardslensbox-reference.vtu b/test/references/richardslensbox-reference.vtu
index 9dff729c81..2d03ab8ffe 100644
--- a/test/references/richardslensbox-reference.vtu
+++ b/test/references/richardslensbox-reference.vtu
@@ -2,45 +2,7 @@
 <VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
   <UnstructuredGrid>
     <Piece NumberOfCells="384" NumberOfPoints="425">
-      <PointData Scalars="Sn">
-        <DataArray type="Float32" Name="Sn" 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 0.999976 0.999964
-          0.999976 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 0.997534 0.909917
-          0.905361 0.909917 0.997534 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 0.999 0.771847
-          0.384776 0.36421 0.384776 0.771847 0.999 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 0.999996 0.918527
-          0.286966 0.0671314 0.0634913 0.0671314 0.286966 0.918527 0.999996 1 1 1 1 1
-          1 1 1 1 1
-        </DataArray>
+      <PointData Scalars="Sw">
         <DataArray type="Float32" Name="Sw" NumberOfComponents="1" format="ascii">
           9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17
           9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17
@@ -79,43 +41,43 @@
           0.713034 0.932869 0.936509 0.932869 0.713034 0.0814728 3.93768e-06 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17
           9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17
         </DataArray>
-        <DataArray type="Float32" Name="pn" NumberOfComponents="1" format="ascii">
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000
+        <DataArray type="Float32" Name="Sn" 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 0.999976 0.999964
+          0.999976 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 0.997534 0.909917
+          0.905361 0.909917 0.997534 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 0.999 0.771847
+          0.384776 0.36421 0.384776 0.771847 0.999 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 0.999996 0.918527
+          0.286966 0.0671314 0.0634913 0.0671314 0.286966 0.918527 0.999996 1 1 1 1 1
+          1 1 1 1 1
         </DataArray>
         <DataArray type="Float32" Name="pw" NumberOfComponents="1" format="ascii">
           97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4
@@ -155,6 +117,44 @@
           99759.4 99835.3 99837.4 99835.3 99759.4 99323.1 97470.5 97470.4 97470.4 97470.4 97470.4 97470.4
           97470.4 97470.4 97470.4 97470.4 97470.4
         </DataArray>
+        <DataArray type="Float32" Name="pn" NumberOfComponents="1" format="ascii">
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000
+        </DataArray>
         <DataArray type="Float32" Name="pc" NumberOfComponents="1" format="ascii">
           2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62
           2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62
@@ -193,7 +193,7 @@
           240.611 164.722 162.618 164.722 240.611 676.884 2529.52 2529.62 2529.62 2529.62 2529.62 2529.62
           2529.62 2529.62 2529.62 2529.62 2529.62
         </DataArray>
-        <DataArray type="Float32" Name="rhoW" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="density" NumberOfComponents="1" format="ascii">
           1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
           1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
           1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
@@ -231,45 +231,7 @@
           1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
           1000 1000 1000 1000 1000
         </DataArray>
-        <DataArray type="Float32" Name="rhoN" NumberOfComponents="1" format="ascii">
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10
-        </DataArray>
-        <DataArray type="Float32" Name="mobW" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="mobility" 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
           0 0 0 0 0 0 0 0 0 0 0 0
@@ -307,44 +269,6 @@
           249.082 733.751 753.199 733.751 249.082 0.0196823 0 0 0 0 0 0
           0 0 0 0 0
         </DataArray>
-        <DataArray type="Float32" Name="mobN" NumberOfComponents="1" format="ascii">
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10
-        </DataArray>
         <DataArray type="Float32" Name="porosity" NumberOfComponents="1" format="ascii">
           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 0.4 0.4 0.4 0.4
@@ -383,44 +307,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
         </DataArray>
-        <DataArray type="Float32" Name="temperature" NumberOfComponents="1" format="ascii">
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15
-        </DataArray>
         <DataArray type="Float32" Name="pressure head" NumberOfComponents="1" format="ascii">
           -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861
           -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861
diff --git a/test/references/richardslenscc-reference.vtu b/test/references/richardslenscc-reference.vtu
index cca30ff245..fde4bda5f6 100644
--- a/test/references/richardslenscc-reference.vtu
+++ b/test/references/richardslenscc-reference.vtu
@@ -2,41 +2,7 @@
 <VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
   <UnstructuredGrid>
     <Piece NumberOfCells="384" NumberOfPoints="425">
-      <CellData Scalars="Sn">
-        <DataArray type="Float32" Name="Sn" 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 0.992854 0.988373 0.988373 0.992854
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 0.996667 0.703082 0.660102 0.660102 0.703082
-          0.996667 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 0.999997 0.920927 0.185263 0.152734 0.152734 0.185263
-          0.920927 0.999997 1 1 1 1 1 1 1 1 1 1
-        </DataArray>
+      <CellData Scalars="Sw">
         <DataArray type="Float32" Name="Sw" NumberOfComponents="1" format="ascii">
           9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17
           9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17
@@ -71,39 +37,39 @@
           9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 2.85357e-06 0.0790732 0.814737 0.847266 0.847266 0.814737
           0.0790732 2.85357e-06 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17
         </DataArray>
-        <DataArray type="Float32" Name="pn" NumberOfComponents="1" format="ascii">
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+        <DataArray type="Float32" Name="Sn" 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 0.992854 0.988373 0.988373 0.992854
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 0.996667 0.703082 0.660102 0.660102 0.703082
+          0.996667 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 0.999997 0.920927 0.185263 0.152734 0.152734 0.185263
+          0.920927 0.999997 1 1 1 1 1 1 1 1 1 1
         </DataArray>
         <DataArray type="Float32" Name="pw" NumberOfComponents="1" format="ascii">
           97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4
@@ -139,6 +105,40 @@
           97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.5 99308.3 99788.3 99798.9 99798.9 99788.3
           99308.3 97470.5 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4
         </DataArray>
+        <DataArray type="Float32" Name="pn" NumberOfComponents="1" format="ascii">
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+        </DataArray>
         <DataArray type="Float32" Name="pc" NumberOfComponents="1" format="ascii">
           2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62
           2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62
@@ -173,7 +173,7 @@
           2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.54 691.737 211.699 201.127 201.127 211.699
           691.737 2529.54 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62
         </DataArray>
-        <DataArray type="Float32" Name="rhoW" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="density" NumberOfComponents="1" format="ascii">
           1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
           1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
           1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
@@ -207,41 +207,7 @@
           1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
           1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
         </DataArray>
-        <DataArray type="Float32" Name="rhoN" NumberOfComponents="1" format="ascii">
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-        </DataArray>
-        <DataArray type="Float32" Name="mobW" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="mobility" 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
           0 0 0 0 0 0 0 0 0 0 0 0
@@ -275,39 +241,39 @@
           0 0 0 0 0 0 0 0.0154609 407.511 473.279 473.279 407.511
           0.0154609 0 0 0 0 0 0 0 0 0 0 0
         </DataArray>
-        <DataArray type="Float32" Name="mobN" NumberOfComponents="1" format="ascii">
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
+        <DataArray type="Float32" Name="kr" 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
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0.0107357 0.0176661 0.0176661 0.0107357
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 1.54609e-05 0.407511 0.473279 0.473279 0.407511
+          1.54609e-05 0 0 0 0 0 0 0 0 0 0 0
         </DataArray>
         <DataArray type="Float32" Name="porosity" NumberOfComponents="1" format="ascii">
           0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4
@@ -343,40 +309,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 0.4 0.4 0.4 0.4
         </DataArray>
-        <DataArray type="Float32" Name="temperature" NumberOfComponents="1" format="ascii">
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-        </DataArray>
         <DataArray type="Float32" Name="pressure head" NumberOfComponents="1" format="ascii">
           -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861
           -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861
-- 
GitLab