From 7ffe99edd182496ee340306e56124e9de156f010 Mon Sep 17 00:00:00 2001
From: Andreas Lauser <and@poware.org>
Date: Fri, 16 Dec 2011 13:23:23 +0000
Subject: [PATCH] some adaptions to the stuff decided today

- GenericFluidState -> NonEquilibriumFluidState (I don't remember the new name of EquilibriumFluidState)
- removed updateAverageMolarMass from NonEquilibriumFluidState and EquilibriumFluidState
- removed sum{Mole|Mass}Fractions from NonEquilibriumFluidState and EquilibriumFluidState
- all fluid states and systems now use enthalpy instead of internal energy

On a completely different front, this also makes the MpNc model
compile with energy enabled, although the implementation is still
incorrect. (as long as nobody is unhappy about it I won't change
this in this branch!)

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@7059 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
 .../1p2cvs2p/watercontaminantfluidsystem.hh   | 14 +--
 dumux/boxmodels/1p2c/1p2cvolumevariables.hh   |  1 -
 .../boxmodels/2p2cni/2p2cnivolumevariables.hh |  4 +-
 dumux/boxmodels/2pni/2pnivolumevariables.hh   |  4 +-
 dumux/boxmodels/MpNc/MpNcvolumevariables.hh   |  1 -
 .../MpNc/energy/MpNcfluxvariablesenergy.hh    | 74 ++++++++-------
 .../MpNc/energy/MpNclocalresidualenergy.hh    |  2 +-
 .../MpNc/energy/MpNcvolumevariablesenergy.hh  |  6 +-
 .../compositionfromfugacities.hh              |  6 --
 .../computefromreferencephase.hh              | 24 +++--
 .../misciblemultiphasecomposition.hh          |  6 +-
 .../MpNcfluidstates/equilibriumfluidstate.hh  | 93 +++++++------------
 .../MpNcfluidstates/immisciblefluidstate.hh   | 33 ++-----
 ...idstate.hh => nonequilibriumfluidstate.hh} | 91 ++++++++----------
 .../MpNcfluidsystems/1pfluidsystem.hh         | 10 +-
 .../2pimmisciblefluidsystem.hh                | 15 +--
 .../MpNcfluidsystems/h2on2fluidsystem.hh      | 17 ++--
 .../1p2c/interstitialfluidtrailfluidsystem.hh |  4 +-
 .../MpNc/obstaclespatialparameters.hh         | 12 +--
 19 files changed, 177 insertions(+), 240 deletions(-)
 rename dumux/material/MpNcfluidstates/{genericfluidstate.hh => nonequilibriumfluidstate.hh} (83%)

diff --git a/appl/lecture/msm/1p2cvs2p/watercontaminantfluidsystem.hh b/appl/lecture/msm/1p2cvs2p/watercontaminantfluidsystem.hh
index 91169019d6..99f490b5ed 100644
--- a/appl/lecture/msm/1p2cvs2p/watercontaminantfluidsystem.hh
+++ b/appl/lecture/msm/1p2cvs2p/watercontaminantfluidsystem.hh
@@ -282,7 +282,7 @@ public:
 
     /*!
      * \brief Given a phase's composition, temperature, pressure and
-     *        density, calculate its specific internal energy [J/kg].
+     *        density, calculate its specific enthalpy [J/kg].
      *
      *  \todo This fluid system neglects the contribution of
      *        gas-molecules in the liquid phase. This contribution is
@@ -294,12 +294,12 @@ public:
      * \param phaseIdx The index of the fluid phase to consider
      */
     template <class FluidState>
-    static Scalar internalEnergy(const FluidState &fluidState,
-                                 const ParameterCache &paramCache,
-                                 int phaseIdx)
+    static Scalar enthalpy(const FluidState &fluidState,
+                           const ParameterCache &paramCache,
+                           int phaseIdx)
     {
         // TODO!
-        DUNE_THROW(Dune::NotImplemented, "Internal energies!");
+        DUNE_THROW(Dune::NotImplemented, "Enthalpy");
     };
 
     /*!
@@ -319,7 +319,7 @@ public:
                                       int phaseIdx)
     {
         // TODO!
-        DUNE_THROW(Dune::NotImplemented, "Thermal conductivity!");
+        DUNE_THROW(Dune::NotImplemented, "Thermal conductivity");
     }
     
       /*!
@@ -336,7 +336,7 @@ public:
                                int phaseIdx)
     {
         // TODO!
-        DUNE_THROW(Dune::NotImplemented, "Heat capacity!");
+        DUNE_THROW(Dune::NotImplemented, "Heat capacity");
     }
 
 };
diff --git a/dumux/boxmodels/1p2c/1p2cvolumevariables.hh b/dumux/boxmodels/1p2c/1p2cvolumevariables.hh
index 55da13154b..4d41328539 100644
--- a/dumux/boxmodels/1p2c/1p2cvolumevariables.hh
+++ b/dumux/boxmodels/1p2c/1p2cvolumevariables.hh
@@ -120,7 +120,6 @@ public:
         }
         fluidState_.setMoleFraction(/*phaseIdx=*/0, /*compIdx=*/0, 1 - x1);
         fluidState_.setMoleFraction(/*phaseIdx=*/0, /*compIdx=*/1, x1);
-        fluidState_.updateAverageMolarMass(/*phaseIdx=*/0);
 
         typename FluidSystem::ParameterCache paramCache;
         paramCache.updatePhase(fluidState_, /*phaseIdx=*/0);
diff --git a/dumux/boxmodels/2p2cni/2p2cnivolumevariables.hh b/dumux/boxmodels/2p2cni/2p2cnivolumevariables.hh
index f7ab0f41d4..b3d040276c 100644
--- a/dumux/boxmodels/2p2cni/2p2cnivolumevariables.hh
+++ b/dumux/boxmodels/2p2cni/2p2cnivolumevariables.hh
@@ -132,9 +132,9 @@ protected:
     {
         // copmute and set the internal energies of the fluid phases
         for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
-            Scalar u = FluidSystem::internalEnergy(this->fluidState_, paramCache, phaseIdx);
+            Scalar h = FluidSystem::enthalpy(this->fluidState_, paramCache, phaseIdx);
 
-            this->fluidState_.setInternalEnergy(phaseIdx, u);
+            this->fluidState_.setEnthalpy(phaseIdx, h);
         }
 
         // copmute and set the heat capacity of the solid phase
diff --git a/dumux/boxmodels/2pni/2pnivolumevariables.hh b/dumux/boxmodels/2pni/2pnivolumevariables.hh
index 297c6d53b8..5399ac08ac 100644
--- a/dumux/boxmodels/2pni/2pnivolumevariables.hh
+++ b/dumux/boxmodels/2pni/2pnivolumevariables.hh
@@ -129,9 +129,9 @@ protected:
     {
         // copmute and set the internal energies of the fluid phases
         for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
-            Scalar u = FluidSystem::internalEnergy(this->fluidState_, paramCache, phaseIdx);
+            Scalar h = FluidSystem::enthalpy(this->fluidState_, paramCache, phaseIdx);
 
-            this->fluidState_.setInternalEnergy(phaseIdx, u);
+            this->fluidState_.setEnthalpy(phaseIdx, h);
         }
 
         // copmute and set the heat capacity of the solid phase
diff --git a/dumux/boxmodels/MpNc/MpNcvolumevariables.hh b/dumux/boxmodels/MpNc/MpNcvolumevariables.hh
index 37038c8c4b..163225029b 100644
--- a/dumux/boxmodels/MpNc/MpNcvolumevariables.hh
+++ b/dumux/boxmodels/MpNc/MpNcvolumevariables.hh
@@ -219,7 +219,6 @@ public:
         // update the remaining parts of the energy module
         EnergyVolumeVariables::update(fluidState_,
                                       paramCache,
-                                      priVars,
                                       element,
                                       elemGeom,
                                       scvIdx,
diff --git a/dumux/boxmodels/MpNc/energy/MpNcfluxvariablesenergy.hh b/dumux/boxmodels/MpNc/energy/MpNcfluxvariablesenergy.hh
index d4589e3b4c..8233476427 100644
--- a/dumux/boxmodels/MpNc/energy/MpNcfluxvariablesenergy.hh
+++ b/dumux/boxmodels/MpNc/energy/MpNcfluxvariablesenergy.hh
@@ -112,7 +112,7 @@ public:
         for (int scvIdx = 0; scvIdx < fvElemGeom.numVertices; scvIdx++)
         {
             tmp = fvElemGeom.subContVolFace[scvfIdx].grad[scvIdx];
-            tmp *= elemVolVars[scvIdx].temperature();
+            tmp *= elemVolVars[scvIdx].fluidState().temperature(/*phaseIdx=*/0);
             temperatureGradient += tmp;
         }
 
@@ -121,14 +121,16 @@ public:
 
 
         lambdaPm_ = lumpedLambdaPm(problem,
+                                   element,
                                     fvElemGeom,
                                     scvfIdx,
                                     elemVolVars) ;
 
     }
 
-     Scalar lumpedLambdaPm(const Problem & problem,
-                              const FVElementGeometry & fvElemGeom,
+    Scalar lumpedLambdaPm(const Problem &problem,
+                          const Element &element,
+                           const FVElementGeometry & fvElemGeom,
                            int faceIdx,
                            const ElementVolumeVariables & elemVolVars)
     {
@@ -136,38 +138,46 @@ public:
          const int i = fvElemGeom.subContVolFace[faceIdx].i;
          const int j = fvElemGeom.subContVolFace[faceIdx].j;
 
-         const Scalar Sli = elemVolVars[i].fluidState().saturation(lPhaseIdx);
-         const Scalar Slj = elemVolVars[j].fluidState().saturation(lPhaseIdx);
+         typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables))::FluidState FluidState;
+         const FluidState &fsI = elemVolVars[i].fluidState();
+         const FluidState &fsJ = elemVolVars[j].fluidState();
+         const Scalar Sli = fsI.saturation(lPhaseIdx);
+         const Scalar Slj = fsJ.saturation(lPhaseIdx);
 
-         const Scalar Sl = std::max<Scalar>(0.0, 0.5*(Sli + Slj));
+         typename FluidSystem::ParameterCache paramCacheI, paramCacheJ;
+         paramCacheI.updateAll(fsI);
+         paramCacheJ.updateAll(fsJ);
 
-        //        const Scalar lambdaDry = 0.583; // W / (K m) // works, orig
-        //        const Scalar lambdaWet = 1.13; // W / (K m) // works, orig
-
-        const typename FluidSystem::MutableParameters mutParams; //dummy
-        const Scalar lambdaDry = 0.5 * (problem.spatialParameters().soilThermalConductivity() + FluidSystem::thermalConductivity(mutParams, gPhaseIdx) ); // W / (K m)
-        const Scalar lambdaWet = 0.5 * (problem.spatialParameters().soilThermalConductivity() + FluidSystem::thermalConductivity(mutParams, lPhaseIdx)) ; // W / (K m)
-
-        // the heat conductivity of the matrix. in general this is a
-        // tensorial value, but we assume isotropic heat conductivity.
-        // This is the Sommerton approach with lambdaDry =
-        // lambdaSn100%.  Taken from: H. Class: "Theorie und
-        // numerische Modellierung nichtisothermer Mehrphasenprozesse
-        // in NAPL-kontaminierten poroesen Medien", PhD Thesis, University of
-        // Stuttgart, Institute of Hydraulic Engineering, p. 57
-
-        Scalar result;
-        if (Sl < 0.1) {
-            // regularization
-            Dumux::Spline<Scalar> sp(0, 0.1, // x1, x2
-                                    0, sqrt(0.1), // y1, y2
-                                    5*0.5/sqrt(0.1), 0.5/sqrt(0.1)); // m1, m2
-            result = lambdaDry + sp.eval(Sl)*(lambdaWet - lambdaDry);
-        }
-        else
-            result = lambdaDry + std::sqrt(Sl)*(lambdaWet - lambdaDry);
+         const Scalar Sl = std::max<Scalar>(0.0, 0.5*(Sli + Slj));
 
-        return result;
+         //        const Scalar lambdaDry = 0.583; // W / (K m) // works, orig
+         //        const Scalar lambdaWet = 1.13; // W / (K m) // works, orig
+         
+         Scalar lambdaSoilI = problem.spatialParameters().soilThermalConductivity(element, fvElemGeom, i);
+         Scalar lambdaSoilJ = problem.spatialParameters().soilThermalConductivity(element, fvElemGeom, i);
+         const Scalar lambdaDry = 0.5 * (lambdaSoilI + FluidSystem::thermalConductivity(fsI, paramCacheI, gPhaseIdx)); // W / (K m)
+         const Scalar lambdaWet = 0.5 * (lambdaSoilJ + FluidSystem::thermalConductivity(fsJ, paramCacheJ, lPhaseIdx)) ; // W / (K m)
+         
+         // the heat conductivity of the matrix. in general this is a
+         // tensorial value, but we assume isotropic heat conductivity.
+         // This is the Sommerton approach with lambdaDry =
+         // lambdaSn100%.  Taken from: H. Class: "Theorie und
+         // numerische Modellierung nichtisothermer Mehrphasenprozesse
+         // in NAPL-kontaminierten poroesen Medien", PhD Thesis, University of
+         // Stuttgart, Institute of Hydraulic Engineering, p. 57
+         
+         Scalar result;
+         if (Sl < 0.1) {
+             // regularization
+             Dumux::Spline<Scalar> sp(0, 0.1, // x1, x2
+                                      0, sqrt(0.1), // y1, y2
+                                      5*0.5/sqrt(0.1), 0.5/sqrt(0.1)); // m1, m2
+             result = lambdaDry + sp.eval(Sl)*(lambdaWet - lambdaDry);
+         }
+         else
+             result = lambdaDry + std::sqrt(Sl)*(lambdaWet - lambdaDry);
+         
+         return result;
     }
 
     /*!
diff --git a/dumux/boxmodels/MpNc/energy/MpNclocalresidualenergy.hh b/dumux/boxmodels/MpNc/energy/MpNclocalresidualenergy.hh
index cae84ef117..45e668dd49 100644
--- a/dumux/boxmodels/MpNc/energy/MpNclocalresidualenergy.hh
+++ b/dumux/boxmodels/MpNc/energy/MpNclocalresidualenergy.hh
@@ -142,7 +142,7 @@ public:
 
         // heat stored in the rock matrix
         result[energyEqIdx] +=
-            volVars.temperature()
+            volVars.fluidState().temperature(/*phaseIdx=*/0)
             * volVars.soilDensity()
             * (1.0 - volVars.porosity())
             * volVars.heatCapacity();
diff --git a/dumux/boxmodels/MpNc/energy/MpNcvolumevariablesenergy.hh b/dumux/boxmodels/MpNc/energy/MpNcvolumevariablesenergy.hh
index 7cd6b7681d..e9bc865e62 100644
--- a/dumux/boxmodels/MpNc/energy/MpNcvolumevariablesenergy.hh
+++ b/dumux/boxmodels/MpNc/energy/MpNcvolumevariablesenergy.hh
@@ -88,7 +88,6 @@ public:
      */
     void update(FluidState &fs,
                 ParameterCache &paramCache,
-                const PrimaryVariables &sol,
                 const Element &element,
                 const FVElementGeometry &elemGeom,
                 int scvIdx,
@@ -174,9 +173,8 @@ public:
 
         // set the enthalpies
         for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
-            Scalar u =
-                FluidSystem::internalEnergy(fs, paramCache, phaseIdx);
-            fs.setInternalEnergy(phaseIdx, u);
+            Scalar h = FluidSystem::enthalpy(fs, paramCache, phaseIdx);
+            fs.setEnthalpy(phaseIdx, h);
         }
     }
 
diff --git a/dumux/material/MpNcconstraintsolvers/compositionfromfugacities.hh b/dumux/material/MpNcconstraintsolvers/compositionfromfugacities.hh
index 594a966679..79cf8f60ed 100644
--- a/dumux/material/MpNcconstraintsolvers/compositionfromfugacities.hh
+++ b/dumux/material/MpNcconstraintsolvers/compositionfromfugacities.hh
@@ -26,7 +26,6 @@
 #ifndef DUMUX_COMPOSITION_FROM_FUGACITIES_HH
 #define DUMUX_COMPOSITION_FROM_FUGACITIES_HH
 
-#include "../MpNcfluidstates/nonequilibriumfluidstate.hh"
 #include <dumux/common/exceptions.hh>
 
 namespace Dumux {
@@ -105,7 +104,6 @@ public:
         // right hand side
         Dune::FieldVector<Scalar, numComponents> b;
 
-        fluidState.updateAverageMolarMass(phaseIdx);
         paramCache.updatePhase(fluidState, phaseIdx);
 
         // maximum number of iterations
@@ -195,7 +193,6 @@ protected:
             fluidState.setMoleFraction(phaseIdx, i, fugacities[i]/gamma);
         };
 
-        fluidState.updateAverageMolarMass(phaseIdx);
         paramCache.updatePhase(fluidState, phaseIdx);
 
         Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
@@ -240,7 +237,6 @@ protected:
             // deviate the mole fraction of the i-th component
             Scalar x_i = fluidState.moleFraction(phaseIdx, i);
             fluidState.setMoleFraction(phaseIdx, i, x_i + eps);
-            fluidState.updateAverageMolarMass(phaseIdx);
             paramCache.updatePhaseSingleMoleFraction(fluidState, phaseIdx, i);
 
             // compute new defect and derivative for all component
@@ -266,7 +262,6 @@ protected:
 
             // reset composition to original value
             fluidState.setMoleFraction(phaseIdx, i, x_i);
-            fluidState.updateAverageMolarMass(phaseIdx);
             paramCache.updatePhaseSingleMoleFraction(fluidState, phaseIdx, i);
 
             // end forward differences
@@ -325,7 +320,6 @@ protected:
                 //sumMoleFrac += std::abs(newx);
             }
 
-            fluidState.updateAverageMolarMass(phaseIdx);
             paramCache.updatePhaseComposition(fluidState, phaseIdx);
 
             /*
diff --git a/dumux/material/MpNcconstraintsolvers/computefromreferencephase.hh b/dumux/material/MpNcconstraintsolvers/computefromreferencephase.hh
index 866b9d7c25..f251dc06c0 100644
--- a/dumux/material/MpNcconstraintsolvers/computefromreferencephase.hh
+++ b/dumux/material/MpNcconstraintsolvers/computefromreferencephase.hh
@@ -30,7 +30,6 @@
 #ifndef DUMUX_COMPUTE_FROM_REFERENCE_PHASE_HH
 #define DUMUX_COMPUTE_FROM_REFERENCE_PHASE_HH
 
-#include "../MpNcfluidstates/genericfluidstate.hh"
 #include <dumux/material/MpNcconstraintsolvers/compositionfromfugacities.hh>
 
 namespace Dumux {
@@ -107,24 +106,23 @@ public:
                       ParameterCache &paramCache,
                       int refPhaseIdx,
                       bool setViscosity,
-                      bool setInternalEnergy)
+                      bool setEnthalpy)
     {
         ComponentVector fugVec;
 
         // compute the density and enthalpy of the
         // reference phase
-        fluidState.updateAverageMolarMass(refPhaseIdx);
         paramCache.updatePhase(fluidState, refPhaseIdx);
         fluidState.setDensity(refPhaseIdx,
                               FluidSystem::density(fluidState, 
                                                    paramCache,
                                                    refPhaseIdx));
 
-        if (setInternalEnergy)
-            fluidState.setInternalEnergy(refPhaseIdx,
-                                         FluidSystem::internalEnergy(fluidState, 
-                                                                     paramCache,
-                                                                     refPhaseIdx));
+        if (setEnthalpy)
+            fluidState.setEnthalpy(refPhaseIdx,
+                                   FluidSystem::enthalpy(fluidState, 
+                                                         paramCache,
+                                                         refPhaseIdx));
 
         if (setViscosity)
             fluidState.setViscosity(refPhaseIdx,
@@ -157,11 +155,11 @@ public:
                                                                paramCache, 
                                                                phaseIdx));
             
-            if (setInternalEnergy)
-                fluidState.setInternalEnergy(phaseIdx,
-                                             FluidSystem::internalEnergy(fluidState,
-                                                                         paramCache,
-                                                                         phaseIdx));
+            if (setEnthalpy)
+                fluidState.setEnthalpy(phaseIdx,
+                                       FluidSystem::enthalpy(fluidState,
+                                                             paramCache,
+                                                             phaseIdx));
         }
     };
 };
diff --git a/dumux/material/MpNcconstraintsolvers/misciblemultiphasecomposition.hh b/dumux/material/MpNcconstraintsolvers/misciblemultiphasecomposition.hh
index 68b97fdf45..d7db5a5512 100644
--- a/dumux/material/MpNcconstraintsolvers/misciblemultiphasecomposition.hh
+++ b/dumux/material/MpNcconstraintsolvers/misciblemultiphasecomposition.hh
@@ -167,8 +167,6 @@ public:
                 int rowIdx = phaseIdx*numComponents + compIdx;
                 fluidState.setMoleFraction(phaseIdx, compIdx, x[rowIdx]);
             }
-            
-            fluidState.updateAverageMolarMass(phaseIdx);
             paramCache.updatePhaseComposition(fluidState, phaseIdx);
         
             Scalar value = FluidSystem::density(fluidState, paramCache, phaseIdx);
@@ -180,8 +178,8 @@ public:
             }
             
             if (setInternalEnergy) {
-                value = FluidSystem::internalEnergy(fluidState, paramCache, phaseIdx);
-                fluidState.setInternalEnergy(phaseIdx, value);
+                value = FluidSystem::enthalpy(fluidState, paramCache, phaseIdx);
+                fluidState.setEnthalpy(phaseIdx, value);
             }
         }
     };
diff --git a/dumux/material/MpNcfluidstates/equilibriumfluidstate.hh b/dumux/material/MpNcfluidstates/equilibriumfluidstate.hh
index c44b79d1be..6568ec3317 100644
--- a/dumux/material/MpNcfluidstates/equilibriumfluidstate.hh
+++ b/dumux/material/MpNcfluidstates/equilibriumfluidstate.hh
@@ -44,7 +44,19 @@ public:
     enum { numComponents = FluidSystem::numComponents };
 
     EquilibriumFluidState()
-    { Valgrind::SetUndefined(*this); }
+    {
+        // set the composition to 0
+        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)  {
+            for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+                moleFraction_[phaseIdx][compIdx] = 0;
+            
+            averageMolarMass_[phaseIdx] = 0;
+            sumMoleFractions_[phaseIdx] = 0;
+        }
+        
+        // make everything undefined so that valgrind will complain
+        Valgrind::SetUndefined(*this);
+    }
 
     template <class FluidState>
     EquilibriumFluidState(FluidState &fs)
@@ -72,29 +84,12 @@ public:
     Scalar massFraction(int phaseIdx, int compIdx) const
     {
         return
-            sumMassFractions(phaseIdx)*
-            molarity(phaseIdx, compIdx)*
-            FluidSystem::molarMass(compIdx)
-            /
-            density(phaseIdx);
+            sumMoleFractions_[phaseIdx]
+            * moleFraction_[phaseIdx][compIdx]
+            * FluidSystem::molarMass(compIdx)
+            / averageMolarMass_[phaseIdx];
     }
 
-    /*!
-     * \brief The sum of all component mole fractions in a phase []
-     *
-     * We define this to be the same as the sum of all mass fractions.
-     */
-    Scalar sumMoleFractions(int phaseIdx) const
-    { return sumMoleFractions_[phaseIdx]; }
-
-    /*!
-     * \brief The sum of all component mass fractions in a phase []
-     *
-     * We define this to be the same as the sum of all mole fractions.
-     */
-    Scalar sumMassFractions(int phaseIdx) const
-    { return sumMoleFractions_[phaseIdx]; }
-
     /*!
      * \brief The average molar mass of a fluid phase [kg/mol]
      */
@@ -159,13 +154,13 @@ public:
      * \brief The specific enthalpy of a fluid phase [J/kg]
      */
     Scalar enthalpy(int phaseIdx) const
-    { return internalEnergy_[phaseIdx] + pressure(phaseIdx)/(density(phaseIdx)); }
+    { return enthalpy_[phaseIdx]; }
 
     /*!
      * \brief The specific internal energy of a fluid phase [J/kg]
      */
     Scalar internalEnergy(int phaseIdx) const
-    { return internalEnergy_[phaseIdx]; }
+    { return enthalpy_[phaseIdx] - pressure(phaseIdx)/density(phaseIdx); }
 
     /*!
      * \brief The dynamic viscosity of a fluid phase [Pa s]
@@ -185,16 +180,6 @@ public:
     Scalar temperature() const
     { return temperature_; }
 
-    /*!
-     * \brief The capillary pressure [Pa] between a phase and a
-     *        reference phase.
-     *
-     * In order to make the term "capillary pressure" meaningful in a
-     * physical sense, mechanic equilibrium needs to be assumed.
-     */
-    Scalar capillaryPressure(int refPhaseIdx, int phaseIdx) const
-    { return pressure_[phaseIdx] - pressure_[refPhaseIdx]; }
-
     /*!
      * \brief The fugacity of a component
      *
@@ -230,7 +215,7 @@ public:
             pressure_[phaseIdx] = fs.pressure(phaseIdx);
             saturation_[phaseIdx] = fs.saturation(phaseIdx);
             density_[phaseIdx] = fs.density(phaseIdx);
-            internalEnergy_[phaseIdx] = fs.internalEnergy(phaseIdx);
+            enthalpy_[phaseIdx] = fs.enthalpy(phaseIdx);
             viscosity_[phaseIdx] = fs.viscosity(phaseIdx);
         }
         temperature_ = fs.temperature(0);
@@ -258,7 +243,19 @@ public:
      * \brief Set the mole fraction of a component in a phase []
      */
     void setMoleFraction(int phaseIdx, int compIdx, Scalar value)
-    { moleFraction_[phaseIdx][compIdx] = value; }
+    { 
+        Valgrind::CheckDefined(value);
+
+        Scalar delta = value - moleFraction_[phaseIdx][compIdx];
+
+        moleFraction_[phaseIdx][compIdx] = value;
+
+        sumMoleFractions_[phaseIdx] += delta;
+        averageMolarMass_[phaseIdx] += delta*FluidSystem::molarMass(compIdx);
+        
+        Valgrind::SetDefined(sumMoleFractions_[phaseIdx]);
+        Valgrind::SetDefined(averageMolarMass_[phaseIdx]);
+    };
 
     /*!
      * \brief Set the fugacity of a component in a phase []
@@ -275,8 +272,8 @@ public:
     /*!
      * \brief Set the specific enthalpy of a phase [J/m^3]
      */
-    void setInternalEnergy(int phaseIdx, Scalar value)
-    { internalEnergy_[phaseIdx] = value; }
+    void setEnthalpy(int phaseIdx, Scalar value)
+    { enthalpy_[phaseIdx] = value; }
 
     /*!
      * \brief Set the dynamic viscosity of a phase [Pa s]
@@ -284,24 +281,6 @@ public:
     void setViscosity(int phaseIdx, Scalar value)
     { viscosity_[phaseIdx] = value; }
 
-    /*!
-     * \brief Calculatate the mean molar mass of a phase given that
-     *        all mole fractions have been set
-     */
-    void updateAverageMolarMass(int phaseIdx)
-    {
-        averageMolarMass_[phaseIdx] = 0;
-        sumMoleFractions_[phaseIdx] = 0;
-
-        for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
-            sumMoleFractions_[phaseIdx] += moleFraction_[phaseIdx][compIdx];
-            averageMolarMass_[phaseIdx] += moleFraction_[phaseIdx][compIdx]*FluidSystem::molarMass(compIdx);
-        }
-        sumMoleFractions_[phaseIdx] = std::max(1e-15, sumMoleFractions_[phaseIdx]);
-        Valgrind::CheckDefined(averageMolarMass_[phaseIdx]);
-        Valgrind::CheckDefined(sumMoleFractions_[phaseIdx]);
-    }
-
     /*!
      * \brief Make sure that all attributes are defined.
      *
@@ -339,7 +318,7 @@ protected:
     Scalar pressure_[numPhases];
     Scalar saturation_[numPhases];
     Scalar density_[numPhases];
-    Scalar internalEnergy_[numPhases];
+    Scalar enthalpy_[numPhases];
     Scalar viscosity_[numPhases];
     Scalar temperature_;
 };
diff --git a/dumux/material/MpNcfluidstates/immisciblefluidstate.hh b/dumux/material/MpNcfluidstates/immisciblefluidstate.hh
index 69836d04a5..eb087e6e70 100644
--- a/dumux/material/MpNcfluidstates/immisciblefluidstate.hh
+++ b/dumux/material/MpNcfluidstates/immisciblefluidstate.hh
@@ -74,22 +74,6 @@ public:
     Scalar massFraction(int phaseIdx, int compIdx) const
     { return (phaseIdx == compIdx)?1.0:0.0; }
 
-    /*!
-     * \brief The sum of all component mole fractions in a phase []
-     *
-     * We define this to be the same as the sum of all mass fractions.
-     */
-    Scalar sumMoleFractions(int phaseIdx) const
-    { return 1.0; }
-
-    /*!
-     * \brief The sum of all component mass fractions in a phase []
-     *
-     * We define this to be the same as the sum of all mole fractions.
-     */
-    Scalar sumMassFractions(int phaseIdx) const
-    { return 1.0; }
-
     /*!
      * \brief The average molar mass of a fluid phase [kg/mol]
      */
@@ -175,13 +159,13 @@ public:
      * \brief The specific enthalpy of a fluid phase [J/kg]
      */
     Scalar enthalpy(int phaseIdx) const
-    { return internalEnergy_[phaseIdx] + pressure(phaseIdx)/(density(phaseIdx)); }
+    { return enthalpy_[phaseIdx]; }
 
     /*!
      * \brief The specific internal energy of a fluid phase [J/kg]
      */
     Scalar internalEnergy(int phaseIdx) const
-    { return internalEnergy_[phaseIdx]; }
+    { return enthalpy_[phaseIdx] - pressure(phaseIdx)/density(phaseIdx); }
 
     /*!
      * \brief The dynamic viscosity of a fluid phase [Pa s]
@@ -189,7 +173,6 @@ public:
     Scalar viscosity(int phaseIdx) const
     { return viscosity_[phaseIdx]; }
 
-
     /*****************************************************
      * Access to fluid properties which only make sense
      * if assuming thermodynamic equilibrium
@@ -231,7 +214,7 @@ public:
             pressure_[phaseIdx] = fs.pressure(phaseIdx);
             saturation_[phaseIdx] = fs.saturation(phaseIdx);
             density_[phaseIdx] = fs.density(phaseIdx);
-            internalEnergy_[phaseIdx] = fs.internalEnergy(phaseIdx);
+            enthalpy_[phaseIdx] = fs.enthalpy(phaseIdx);
             viscosity_[phaseIdx] = fs.viscosity(phaseIdx);
         }
         temperature_ = fs.temperature(0);
@@ -262,10 +245,10 @@ public:
     { density_[phaseIdx] = value; }
 
     /*!
-     * \brief Set the specific internal energy of a phase [J/m^3]
+     * \brief Set the specific enthalpy of a phase [J/m^3]
      */
-    void setInternalEnergy(int phaseIdx, Scalar value)
-    { internalEnergy_[phaseIdx] = value; }
+    void setEnthalpy(int phaseIdx, Scalar value)
+    { enthalpy_[phaseIdx] = value; }
 
     /*!
      * \brief Set the dynamic viscosity of a phase [Pa s]
@@ -288,7 +271,7 @@ public:
             //for (int j = 0; j < numComponents; ++j) {
             //    Valgrind::CheckDefined(fugacityCoefficient_[i][j]);
             //}
-            Valgrind::CheckDefined(pressure_[i]);
+            Valgrind:CheckDefined(pressure_[i]);
             Valgrind::CheckDefined(saturation_[i]);
             Valgrind::CheckDefined(density_[i]);
             //Valgrind::CheckDefined(internalEnergy_[i]);
@@ -303,7 +286,7 @@ protected:
     Scalar pressure_[numPhases];
     Scalar saturation_[numPhases];
     Scalar density_[numPhases];
-    Scalar internalEnergy_[numPhases];
+    Scalar enthalpy_[numPhases];
     Scalar viscosity_[numPhases];
     Scalar temperature_;
 };
diff --git a/dumux/material/MpNcfluidstates/genericfluidstate.hh b/dumux/material/MpNcfluidstates/nonequilibriumfluidstate.hh
similarity index 83%
rename from dumux/material/MpNcfluidstates/genericfluidstate.hh
rename to dumux/material/MpNcfluidstates/nonequilibriumfluidstate.hh
index 9b4778c7c3..c6e64b4bdf 100644
--- a/dumux/material/MpNcfluidstates/genericfluidstate.hh
+++ b/dumux/material/MpNcfluidstates/nonequilibriumfluidstate.hh
@@ -40,14 +40,26 @@ namespace Dumux
  *        any assumptions.
  */
 template <class Scalar, class FluidSystem>
-class GenericFluidState
+class NonEquilibriumFluidState
 {
 public:
     enum { numPhases = FluidSystem::numPhases };
     enum { numComponents = FluidSystem::numComponents };
 
-    GenericFluidState()
-    { Valgrind::SetUndefined(*this); }
+    NonEquilibriumFluidState()
+    { 
+        // set the composition to 0
+        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)  {
+            for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+                moleFractions_[phaseIdx][compIdx] = 0;
+            
+            averageMolarMass_[phaseIdx] = 0;
+            sumMoleFractions_[phaseIdx] = 0;
+        }
+        
+        // make everything undefined so that valgrind will complain
+        Valgrind::SetUndefined(*this);
+    }
 
     /*****************************************************
      * Generic access to fluid properties (No assumptions
@@ -72,28 +84,12 @@ public:
     Scalar massFraction(int phaseIdx, int compIdx) const
     {
         return
-            sumMassFraction(phaseIdx)
+            sumMoleFractions_[phaseIdx]
             * moleFraction_[phaseIdx][compIdx]
             * FluidSystem::molarMass(compIdx)
             / averageMolarMass_[phaseIdx];
     }
 
-    /*!
-     * \brief The sum of all component mole fractions in a phase []
-     *
-     * We define this to be the same as the sum of all mass fractions.
-     */
-    Scalar sumMoleFractions(int phaseIdx) const
-    { return sumMoleFractions_[phaseIdx]; }
-
-    /*!
-     * \brief The sum of all component mass fractions in a phase []
-     *
-     * We define this to be the same as the sum of all mole fractions.
-     */
-    Scalar sumMassFraction(int phaseIdx) const
-    { return sumMoleFractions_[phaseIdx]; }
-
     /*!
      * \brief The mean molar mass of a fluid phase [kg/mol]
      */
@@ -153,18 +149,17 @@ public:
      */
     Scalar pressure(int phaseIdx) const
     { return pressure_[phaseIdx]; };
-
     /*!
-     * \brief The specific internal energy of a fluid phase [J/kg]
+     * \brief The specific enthalpy of a fluid phase [J/kg]
      */
-    Scalar internalEnergy(int phaseIdx) const
-    { return internalEnergy_[phaseIdx]; };
+    Scalar enthalpy(int phaseIdx) const
+    { return enthalpy_[phaseIdx]; }
 
     /*!
-     * \brief The specific enthalpy of a fluid phase [J/kg]
+     * \brief The specific internal energy of a fluid phase [J/kg]
      */
-    Scalar enthalpy(int phaseIdx) const
-    { return internalEnergy(phaseIdx) + pressure(phaseIdx)/density(phaseIdx); };
+    Scalar internalEnergy(int phaseIdx) const
+    { return enthalpy_[phaseIdx] - pressure(phaseIdx)/density(phaseIdx); }
 
     /*!
      * \brief The dynamic viscosity of a fluid phase [Pa s]
@@ -199,7 +194,19 @@ public:
      * \brief Set the mole fraction of a component in a phase []
      */
     void setMoleFraction(int phaseIdx, int compIdx, Scalar value)
-    { moleFraction_[phaseIdx][compIdx] = value; }
+    { 
+        Valgrind::CheckDefined(value);
+
+        Scalar delta = value - moleFraction_[phaseIdx][compIdx];
+
+        moleFraction_[phaseIdx][compIdx] = value;
+
+        sumMoleFractions_[phaseIdx] += delta;
+        averageMolarMass_[phaseIdx] += delta*FluidSystem::molarMass(compIdx);
+        
+        Valgrind::SetDefined(sumMolarMass_[phaseIdx]);
+        Valgrind::SetDefined(averageMolarMass_[phaseIdx]);
+    }
 
     /*!
      * \brief Set the fugacity of a component in a phase []
@@ -214,10 +221,10 @@ public:
     { density_[phaseIdx] = value; }
 
     /*!
-     * \brief Set the specific internal energy of a phase [J/m^3]
+     * \brief Set the specific enthalpy of a phase [J/m^3]
      */
-    void setInternalEnergy(int phaseIdx, Scalar value)
-    { internalEnergy_[phaseIdx] = value; }
+    void setEnthalpy(int phaseIdx, Scalar value)
+    { enthalpy_[phaseIdx] = value; }
 
     /*!
      * \brief Set the dynamic viscosity of a phase [Pa s]
@@ -225,24 +232,6 @@ public:
     void setViscosity(int phaseIdx, Scalar value)
     { viscosity_[phaseIdx] = value; }
 
-    /*!
-     * \brief Calculatate the mean molar mass of a phase given that
-     *        all mole fractions have been set
-     */
-    void updateAverageMolarMass(int phaseIdx)
-    {
-        averageMolarMass_[phaseIdx] = 0;
-        sumMoleFractions_[phaseIdx] = 0;
-
-        for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
-            sumMoleFractions_[phaseIdx] += moleFraction_[phaseIdx][compIdx];
-            averageMolarMass_[phaseIdx] += moleFraction_[phaseIdx][compIdx]*FluidSystem::molarMass(compIdx);
-        }
-        sumMoleFractions_[phaseIdx] = std::max(1e-15, sumMoleFractions_[phaseIdx]);
-        Valgrind::CheckDefined(averageMolarMass_[phaseIdx]);
-        Valgrind::CheckDefined(sumMoleFractions_[phaseIdx]);
-    }
-
     /*!
      * \brief Retrieve all parameters from an arbitrary fluid
      *        state.
@@ -260,7 +249,7 @@ public:
             pressure_[phaseIdx] = fs.pressure(phaseIdx);
             saturation_[phaseIdx] = fs.saturation(phaseIdx);
             density_[phaseIdx] = fs.density(phaseIdx);
-            internalEnergy_[phaseIdx] = fs.internalEnergy(phaseIdx);
+            enthalpy_[phaseIdx] = fs.enthalpy(phaseIdx);
             viscosity_[phaseIdx] = fs.viscosity(phaseIdx);
             temperature_[phaseIdx] = fs.temperature(phaseIdx);
         }
@@ -302,7 +291,7 @@ protected:
     Scalar pressure_[numPhases];
     Scalar saturation_[numPhases];
     Scalar density_[numPhases];
-    Scalar internalEnergy_[numPhases];
+    Scalar enthalpy_[numPhases];
     Scalar viscosity_[numPhases];
     Scalar temperature_[numPhases];
 };
diff --git a/dumux/material/MpNcfluidsystems/1pfluidsystem.hh b/dumux/material/MpNcfluidsystems/1pfluidsystem.hh
index bac653903f..9fcd1a85ec 100644
--- a/dumux/material/MpNcfluidsystems/1pfluidsystem.hh
+++ b/dumux/material/MpNcfluidsystems/1pfluidsystem.hh
@@ -303,22 +303,22 @@ public:
  
     /*!
      * \brief Given a phase's composition, temperature, pressure and
-     *        density, calculate its specific internal energy [J/kg].
+     *        density, calculate its specific enthalpy [J/kg].
      *
      * \param fluidState An abitrary fluid state
      * \param paramCache The fluid system's parameter cache
      * \param phaseIdx The index of the fluid phase to consider
      */
     template <class FluidState>
-    static Scalar internalEnergy(const FluidState &fluidState,
-                                 const ParameterCache &paramCache,
-                                 int phaseIdx)
+    static Scalar enthalpy(const FluidState &fluidState,
+                           const ParameterCache &paramCache,
+                           int phaseIdx)
     {
         assert(0 <= phaseIdx && phaseIdx < numPhases);
 
         Scalar temperature = fluidState.temperature(phaseIdx);
         Scalar pressure = fluidState.pressure(phaseIdx);
-        return Fluid::internalEnergy(temperature, pressure);
+        return Fluid::enthalpy(temperature, pressure);
     }
 
     /*!
diff --git a/dumux/material/MpNcfluidsystems/2pimmisciblefluidsystem.hh b/dumux/material/MpNcfluidsystems/2pimmisciblefluidsystem.hh
index dc36667a27..70196ee0e1 100644
--- a/dumux/material/MpNcfluidsystems/2pimmisciblefluidsystem.hh
+++ b/dumux/material/MpNcfluidsystems/2pimmisciblefluidsystem.hh
@@ -329,18 +329,11 @@ public:
     };
 
     /*!
-     * \brief Return the specific internal of a fluid phase [J/kg].
-     *
-     * \param phaseIdx index of the phase
-     * \param temperature phase temperature in [K]
-     * \param pressure phase pressure in [Pa]
-     * \param fluidState The fluid state of the two-phase model
-     * \tparam FluidState the fluid state class of the two-phase model
-     * \return returns the specific enthalpy of the phase [J/kg]
+     * \brief Return the specific enthalpy of a fluid phase [J/kg].
      */
     using Base::internalEnergy;
     template <class FluidState>
-    static Scalar internalEnergy(const FluidState &fluidState,
+    static Scalar enthalpy(const FluidState &fluidState,
                                  int phaseIdx)
     {
         assert(0 <= phaseIdx  && phaseIdx < numPhases);
@@ -348,8 +341,8 @@ public:
         Scalar temperature = fluidState.temperature(phaseIdx);
         Scalar pressure = fluidState.pressure(phaseIdx);
         if (phaseIdx == wPhaseIdx)
-            return WettingPhase::internalEnergy(temperature, pressure);
-        return NonWettingPhase::internalEnergy(temperature, pressure);
+            return WettingPhase::enthalpy(temperature, pressure);
+        return NonWettingPhase::enthalpy(temperature, pressure);
     }
 
     /*!
diff --git a/dumux/material/MpNcfluidsystems/h2on2fluidsystem.hh b/dumux/material/MpNcfluidsystems/h2on2fluidsystem.hh
index f643546ed4..e5b4f0ed38 100644
--- a/dumux/material/MpNcfluidsystems/h2on2fluidsystem.hh
+++ b/dumux/material/MpNcfluidsystems/h2on2fluidsystem.hh
@@ -25,8 +25,6 @@
 #ifndef DUMUX_H2O_N2_FLUID_SYSTEM_HH
 #define DUMUX_H2O_N2_FLUID_SYSTEM_HH
 
-#include <dumux/material/MpNcfluidstates/nonequilibriumfluidstate.hh>
-
 #include <dumux/material/components/simpleh2o.hh>
 #include <dumux/material/components/h2o.hh>
 #include <dumux/material/components/n2.hh>
@@ -34,7 +32,6 @@
 #include <dumux/material/idealgas.hh>
 
 #include <dumux/material/binarycoefficients/h2o_n2.hh>
-#include <dumux/material/MpNcfluidstates/nonequilibriumfluidstate.hh>
 
 #include <dumux/common/valgrind.hh>
 #include <dumux/common/exceptions.hh>
@@ -456,7 +453,7 @@ public:
 
     /*!
      * \brief Given a phase's composition, temperature, pressure and
-     *        density, calculate its specific internal energy [J/kg].
+     *        density, calculate its specific enthalpy [J/kg].
      *
      *  \todo This fluid system neglects the contribution of
      *        gas-molecules in the liquid phase. This contribution is
@@ -468,9 +465,9 @@ public:
      * \param phaseIdx The index of the fluid phase to consider
      */
     template <class FluidState>
-    static Scalar internalEnergy(const FluidState &fluidState,
-                                 const ParameterCache &paramCache,
-                                 int phaseIdx)
+    static Scalar enthalpy(const FluidState &fluidState,
+                           const ParameterCache &paramCache,
+                           int phaseIdx)
     {
         Scalar T = fluidState.temperature(phaseIdx);
         Scalar p = fluidState.pressure(phaseIdx);
@@ -479,15 +476,15 @@ public:
         if (phaseIdx == lPhaseIdx) {
             // TODO: correct way to deal with the solutes???
             return
-                H2O::liquidInternalEnergy(T, p) ;
+                H2O::liquidEnthalpy(T, p) ;
         }
         else {
             // assume ideal gas
             Scalar XH2O = fluidState.massFraction(gPhaseIdx, H2OIdx);
             Scalar XN2 = fluidState.massFraction(gPhaseIdx, N2Idx);
             Scalar result = 0;
-            result += XH2O*H2O::gasInternalEnergy(T, p);
-            result += XN2*N2::gasInternalEnergy(T, p);
+            result += XH2O*H2O::gasEnthalpy(T, p);
+            result += XN2*N2::gasEnthalpy(T, p);
             return result;
         }
     }
diff --git a/test/boxmodels/1p2c/interstitialfluidtrailfluidsystem.hh b/test/boxmodels/1p2c/interstitialfluidtrailfluidsystem.hh
index a7aa775159..ab2b6bd91e 100644
--- a/test/boxmodels/1p2c/interstitialfluidtrailfluidsystem.hh
+++ b/test/boxmodels/1p2c/interstitialfluidtrailfluidsystem.hh
@@ -301,13 +301,13 @@ public:
      * \param phaseIdx  for which phase to give back the heat capacity
      */
     template <class FluidState>
-    static Scalar internalEnergy(const FluidState &fluidState,
+    static Scalar enthalpy(const FluidState &fluidState,
                                  const ParameterCache &paramCache,
                                  int phaseIdx)
     {
         assert(0 <= phaseIdx && phaseIdx < numPhases);
 
-        DUNE_THROW(Dune::NotImplemented, "Internal energies");
+        DUNE_THROW(Dune::NotImplemented, "Enthalpies");
     };
 
     /*!
diff --git a/test/boxmodels/MpNc/obstaclespatialparameters.hh b/test/boxmodels/MpNc/obstaclespatialparameters.hh
index 60a2d11703..ddbac95ebd 100644
--- a/test/boxmodels/MpNc/obstaclespatialparameters.hh
+++ b/test/boxmodels/MpNc/obstaclespatialparameters.hh
@@ -272,16 +272,16 @@ public:
             return coarseMaterialParams_;
     }
 
-    const Scalar & soilDensity(const Element &element,
-                                const FVElementGeometry &fvElemGeom,
-                                int scvIdx) const
+    Scalar soilDensity(const Element &element,
+                       const FVElementGeometry &fvElemGeom,
+                       int scvIdx) const
     {
         return 2700. ; // density of granite [kg/m^3]
     }
 
-    const Scalar & soilThermalConductivity(const Element &element,
-                                    const FVElementGeometry &fvElemGeom,
-                                    int scvIdx) const
+    Scalar soilThermalConductivity(const Element &element,
+                                   const FVElementGeometry &fvElemGeom,
+                                   int scvIdx) const
     {
         return 2.8; // conductivity of granite [W / (m K ) ]
     }
-- 
GitLab