From 5ec5aaa0780c25bfb1a71d376172e20255b4ccca Mon Sep 17 00:00:00 2001
From: Andreas Lauser <and@poware.org>
Date: Mon, 31 Oct 2011 14:59:52 +0000
Subject: [PATCH] new fluid systems: implement the changes discussed last
 friday

If you have any major objections, speak at meeting on friday or
forever hold your peace.  (Where "forever" means "until the next major
fluid framework refactoring". [Which won't be initiated by me.])

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@6807 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
 dumux/boxmodels/MpNc/MpNcproperties.hh        |   2 +-
 dumux/boxmodels/MpNc/MpNcpropertydefaults.hh  |   2 +-
 dumux/boxmodels/MpNc/MpNcvolumevariables.hh   |  61 ++-
 dumux/boxmodels/MpNc/MpNcvolumevariablesia.hh |   5 +-
 dumux/boxmodels/MpNc/MpNcvtkwritercommon.hh   |  18 +-
 .../MpNc/diffusion/volumevariables.hh         |  11 +-
 .../MpNc/energy/MpNcvolumevariablesenergy.hh  |  92 ++--
 .../MpNc/mass/MpNcvolumevariablesmass.hh      |  30 +-
 .../compositionfromfugacities.hh              | 138 +++---
 .../MpNcfluidstates/equilibriumfluidstate.hh  | 124 +++---
 .../MpNcfluidstates/genericfluidstate.hh      | 132 +++---
 .../MpNcfluidstates/immisciblefluidstate.hh   |  70 ++-
 .../trackinggenericfluidstate.hh              | 398 ------------------
 .../MpNcfluidsystems/h2on2fluidsystem.hh      | 329 +++++++--------
 .../MpNcfluidsystems/nullparametercache.hh    |  56 +++
 .../MpNcfluidsystems/parametercachebase.hh    | 120 ++++++
 dumux/material/idealgas.hh                    |  17 +-
 17 files changed, 673 insertions(+), 932 deletions(-)
 delete mode 100644 dumux/material/MpNcfluidstates/trackinggenericfluidstate.hh
 create mode 100644 dumux/material/MpNcfluidsystems/nullparametercache.hh
 create mode 100644 dumux/material/MpNcfluidsystems/parametercachebase.hh

diff --git a/dumux/boxmodels/MpNc/MpNcproperties.hh b/dumux/boxmodels/MpNc/MpNcproperties.hh
index 7cc4bf98cb..7050f9c65a 100644
--- a/dumux/boxmodels/MpNc/MpNcproperties.hh
+++ b/dumux/boxmodels/MpNc/MpNcproperties.hh
@@ -67,7 +67,7 @@ NEW_PROP_TAG(MPNCVtkAddVarPressures);
 NEW_PROP_TAG(MPNCVtkAddVelocities);
 NEW_PROP_TAG(MPNCVtkAddDensities);
 NEW_PROP_TAG(MPNCVtkAddMobilities);
-NEW_PROP_TAG(MPNCVtkAddMeanMolarMass);
+NEW_PROP_TAG(MPNCVtkAddAverageMolarMass);
 NEW_PROP_TAG(MPNCVtkAddMassFractions);
 NEW_PROP_TAG(MPNCVtkAddMoleFractions);
 NEW_PROP_TAG(MPNCVtkAddMolarities);
diff --git a/dumux/boxmodels/MpNc/MpNcpropertydefaults.hh b/dumux/boxmodels/MpNc/MpNcpropertydefaults.hh
index b4006423ec..e67ab9defb 100644
--- a/dumux/boxmodels/MpNc/MpNcpropertydefaults.hh
+++ b/dumux/boxmodels/MpNc/MpNcpropertydefaults.hh
@@ -214,7 +214,7 @@ SET_BOOL_PROP(BoxMPNC, MPNCVtkAddVarPressures, false);
 SET_BOOL_PROP(BoxMPNC, MPNCVtkAddVelocities, false);
 SET_BOOL_PROP(BoxMPNC, MPNCVtkAddDensities, true);
 SET_BOOL_PROP(BoxMPNC, MPNCVtkAddMobilities, true);
-SET_BOOL_PROP(BoxMPNC, MPNCVtkAddMeanMolarMass, false);
+SET_BOOL_PROP(BoxMPNC, MPNCVtkAddAverageMolarMass, false);
 SET_BOOL_PROP(BoxMPNC, MPNCVtkAddMassFractions, false);
 SET_BOOL_PROP(BoxMPNC, MPNCVtkAddMoleFractions, true);
 SET_BOOL_PROP(BoxMPNC, MPNCVtkAddMolarities, false);
diff --git a/dumux/boxmodels/MpNc/MpNcvolumevariables.hh b/dumux/boxmodels/MpNc/MpNcvolumevariables.hh
index 38253ce229..55e846d944 100644
--- a/dumux/boxmodels/MpNc/MpNcvolumevariables.hh
+++ b/dumux/boxmodels/MpNc/MpNcvolumevariables.hh
@@ -90,7 +90,6 @@ class MPNCVolumeVariables
 
     typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
 
-
 public:
     //! The return type of the fluidState() method
     typedef typename MassVolumeVariables::FluidState FluidState;
@@ -124,8 +123,9 @@ public:
                            elemGeom,
                            scvIdx,
                            isOldSol);
+        ParentType::checkDefined();
 
-        typename FluidSystem::MutableParameters mutParams;
+        typename FluidSystem::ParameterCache paramCache;
 
         /////////////
         // set the phase saturations
@@ -133,24 +133,21 @@ public:
         Scalar sumSat = 0;
         for (int i = 0; i < numPhases - 1; ++i) {
             sumSat += priVars[S0Idx + i];
-            mutParams.setSaturation(i, priVars[S0Idx + i]);
+            fluidState_.setSaturation(i, priVars[S0Idx + i]);
         }
         Valgrind::CheckDefined(sumSat);
-        mutParams.setSaturation(numPhases - 1, 1.0 - sumSat);
+        fluidState_.setSaturation(numPhases - 1, 1.0 - sumSat);
 
         /////////////
         // set the fluid phase temperatures
         /////////////
-        // update the temperature part of the energy module
-        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx){
-            Scalar T = EnergyVolumeVariables::getTemperature(priVars,
-                                                             element,
-                                                             elemGeom,
-                                                             scvIdx,
-                                                             problem,
-                                                             phaseIdx);
-            mutParams.setTemperature(phaseIdx, T);
-        }
+        EnergyVolumeVariables::updateTemperatures(fluidState_,
+                                                  paramCache,
+                                                  priVars,
+                                                  element,
+                                                  elemGeom,
+                                                  scvIdx,
+                                                  problem);
 
         /////////////
         // set the phase pressures
@@ -161,16 +158,17 @@ public:
             problem.spatialParameters().materialLawParams(element, elemGeom, scvIdx);
         // capillary pressures
         Scalar capPress[numPhases];
-        MaterialLaw::capillaryPressures(capPress, materialParams, mutParams);
+        MaterialLaw::capillaryPressures(capPress, materialParams, fluidState_);
         // add to the pressure of the first fluid phase
         Scalar p0 = priVars[p0Idx];
         for (int i = 0; i < numPhases; ++ i)
-            mutParams.setPressure(i, p0 - capPress[0] + capPress[i]);
+            fluidState_.setPressure(i, p0 - capPress[0] + capPress[i]);
 
         /////////////
         // set the fluid compositions
         /////////////
-        MassVolumeVariables::update(mutParams,
+        MassVolumeVariables::update(fluidState_,
+                                    paramCache,
                                     priVars,
                                     hint_,
                                     problem,
@@ -178,7 +176,6 @@ public:
                                     elemGeom,
                                     scvIdx);
         MassVolumeVariables::checkDefined();
-        ParentType::checkDefined();
 
         /////////////
         // Porosity
@@ -189,7 +186,6 @@ public:
                                                          elemGeom,
                                                          scvIdx);
         Valgrind::CheckDefined(porosity_);
-        ParentType::checkDefined();
 
         /////////////
         // Phase mobilities
@@ -198,31 +194,30 @@ public:
         // relative permeabilities
         MaterialLaw::relativePermeabilities(relativePermeability_,
                                             materialParams, 
-                                            mutParams);
+                                            fluidState_);
 
         // dynamic viscosities
         for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
             // viscosities
-            Scalar mu = FluidSystem::computeViscosity(mutParams, phaseIdx);
-            mutParams.setViscosity(phaseIdx, mu);
+            Scalar mu = FluidSystem::viscosity(fluidState_, paramCache, phaseIdx);
+            fluidState_.setViscosity(phaseIdx, mu);
         }
-        ParentType::checkDefined();
 
         /////////////
         // diffusion
         /////////////
 
         // update the diffusion part of the volume data
-        DiffusionVolumeVariables::update(mutParams, *this, problem);
+        DiffusionVolumeVariables::update(fluidState_, paramCache, *this, problem);
         DiffusionVolumeVariables::checkDefined();
-        ParentType::checkDefined();
 
         /////////////
         // energy
         /////////////
 
         // update the remaining parts of the energy module
-        EnergyVolumeVariables::update(mutParams,
+        EnergyVolumeVariables::update(fluidState_,
+                                      paramCache,
                                       priVars,
                                       element,
                                       elemGeom,
@@ -230,18 +225,16 @@ public:
                                       problem);
 
         EnergyVolumeVariables::checkDefined();
-        ParentType::checkDefined();
-
 
-        // assign the fluid state to be stored (needed in the next update)
-        fluidState_.assign(mutParams);
+        // make sure the quantities in the fluid state are well-defined
         fluidState_.checkDefined();
 
         // specific interfacial area,
         // well also all the dimensionless numbers :-)
         // well, also the mass transfer rate
         IAVolumeVariables::update(*this,
-                                  mutParams,
+                                  fluidState_,
+                                  paramCache,
                                   priVars,
                                   problem,
                                   element,
@@ -249,11 +242,7 @@ public:
                                   scvIdx);
         IAVolumeVariables::checkDefined();
 
-        ParentType::checkDefined();
-
         checkDefined();
-
-        ParentType::checkDefined();
     }
 
     /*!
@@ -330,7 +319,7 @@ public:
         // difference of sum of mole fractions in the phase from 100%
         Scalar a = 1;
         for (int i = 0; i < numComponents; ++i)
-            a -= fluidState.moleFrac(phaseIdx, i);
+            a -= fluidState.moleFraction(phaseIdx, i);
         return a;
     }
 
diff --git a/dumux/boxmodels/MpNc/MpNcvolumevariablesia.hh b/dumux/boxmodels/MpNc/MpNcvolumevariablesia.hh
index 23764d1a10..ffc2b3e210 100644
--- a/dumux/boxmodels/MpNc/MpNcvolumevariablesia.hh
+++ b/dumux/boxmodels/MpNc/MpNcvolumevariablesia.hh
@@ -61,9 +61,10 @@ public:
     /*!
      * \brief Updates the volume specific interfacial area [m^2 / m^3] between the phases.
      */
-    template <class MutableParams>
+    template <class FluidState, class ParameterCache>
     void update(const VolumeVariables & volVars,
-                const MutableParams & mutParams,
+                const FluidState &fluidState,
+                const ParameterCache &paramCache,
                 const PrimaryVariables &priVars,
                 const Problem &problem,
                 const Element & element,
diff --git a/dumux/boxmodels/MpNc/MpNcvtkwritercommon.hh b/dumux/boxmodels/MpNc/MpNcvtkwritercommon.hh
index a6b53eedae..9d751c5a67 100644
--- a/dumux/boxmodels/MpNc/MpNcvtkwritercommon.hh
+++ b/dumux/boxmodels/MpNc/MpNcvtkwritercommon.hh
@@ -85,7 +85,7 @@ public:
         velocityOutput_ = GET_PARAM(TypeTag, bool, MPNC, VtkAddVelocities);
         densityOutput_ = GET_PARAM(TypeTag, bool, MPNC, VtkAddDensities);
         mobilityOutput_ = GET_PARAM(TypeTag, bool, MPNC, VtkAddMobilities);
-        meanMolarMassOutput_ = GET_PARAM(TypeTag, bool, MPNC, VtkAddMeanMolarMass);
+        averageMolarMassOutput_ = GET_PARAM(TypeTag, bool, MPNC, VtkAddAverageMolarMass);
         massFracOutput_ = GET_PARAM(TypeTag, bool, MPNC, VtkAddMassFractions);
         moleFracOutput_ = GET_PARAM(TypeTag, bool, MPNC, VtkAddMoleFractions);
         molarityOutput_ = GET_PARAM(TypeTag, bool, MPNC, VtkAddMolarities);
@@ -116,7 +116,7 @@ public:
         if (pressureOutput_) this->resizePhaseBuffer_(pressure_);
         if (densityOutput_) this->resizePhaseBuffer_(density_);
         if (mobilityOutput_) this->resizePhaseBuffer_(mobility_);
-        if (meanMolarMassOutput_) this->resizePhaseBuffer_(meanMolarMass_);
+        if (averageMolarMassOutput_) this->resizePhaseBuffer_(averageMolarMass_);
 
         if (moleFracOutput_) this->resizePhaseComponentBuffer_(moleFrac_);
         if (massFracOutput_) this->resizePhaseComponentBuffer_(massFrac_);
@@ -154,10 +154,10 @@ public:
                 if (pressureOutput_) pressure_[phaseIdx][I] = volVars.fluidState().pressure(phaseIdx);
                 if (densityOutput_) density_[phaseIdx][I] = volVars.fluidState().density(phaseIdx);
                 if (mobilityOutput_) mobility_[phaseIdx][I] = volVars.mobility(phaseIdx);
-                if (meanMolarMassOutput_) meanMolarMass_[phaseIdx][I] = volVars.fluidState().meanMolarMass(phaseIdx);
+                if (averageMolarMassOutput_) averageMolarMass_[phaseIdx][I] = volVars.fluidState().averageMolarMass(phaseIdx);
                 for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
-                    if (moleFracOutput_) moleFrac_[phaseIdx][compIdx][I] = volVars.fluidState().moleFrac(phaseIdx, compIdx);
-                    if (massFracOutput_) massFrac_[phaseIdx][compIdx][I] = volVars.fluidState().massFrac(phaseIdx, compIdx);
+                    if (moleFracOutput_) moleFrac_[phaseIdx][compIdx][I] = volVars.fluidState().moleFraction(phaseIdx, compIdx);
+                    if (massFracOutput_) massFrac_[phaseIdx][compIdx][I] = volVars.fluidState().massFraction(phaseIdx, compIdx);
                     if (molarityOutput_) molarity_[phaseIdx][compIdx][I] = volVars.fluidState().molarity(phaseIdx, compIdx);
                 }
             }
@@ -213,8 +213,8 @@ public:
         if (densityOutput_)
             this->commitPhaseBuffer_(writer, "rho_%s", density_);
 
-        if (meanMolarMassOutput_)
-            this->commitPhaseBuffer_(writer, "M_%s", meanMolarMass_);
+        if (averageMolarMassOutput_)
+            this->commitPhaseBuffer_(writer, "M_%s", averageMolarMass_);
 
         if (mobilityOutput_)
             this->commitPhaseBuffer_(writer, "lambda_%s", mobility_);
@@ -263,13 +263,13 @@ private:
     bool massFracOutput_;
     bool moleFracOutput_;
     bool molarityOutput_;
-    bool meanMolarMassOutput_;
+    bool averageMolarMassOutput_;
 
     PhaseBuffer saturation_;
     PhaseBuffer pressure_;
     PhaseBuffer density_;
     PhaseBuffer mobility_;
-    PhaseBuffer meanMolarMass_;
+    PhaseBuffer averageMolarMass_;
 
     PhaseVelocityField velocity_;
     ScalarBuffer boxSurface_;
diff --git a/dumux/boxmodels/MpNc/diffusion/volumevariables.hh b/dumux/boxmodels/MpNc/diffusion/volumevariables.hh
index 096724bc9b..b8be021382 100644
--- a/dumux/boxmodels/MpNc/diffusion/volumevariables.hh
+++ b/dumux/boxmodels/MpNc/diffusion/volumevariables.hh
@@ -129,18 +129,17 @@ class MPNCVolumeVariablesDiffusion<TypeTag, false>
     typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
     typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
     typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
-    typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
-    typedef typename FluidSystem::MutableParameters MutableParameters;
 
 public:
     MPNCVolumeVariablesDiffusion()
     {}
 
-    void update(MutableParameters &mutParams,
-        const VolumeVariables &volVars,
+    template <class FluidState, class ParameterCache>
+    void update(FluidState &fluidState,
+                ParameterCache &paramCache,
+                const VolumeVariables &volVars,
                 const Problem &problem)
-    {
-    };
+    { };
 
     Scalar diffCoeffL(int compIdx) const
     { return 0; };
diff --git a/dumux/boxmodels/MpNc/energy/MpNcvolumevariablesenergy.hh b/dumux/boxmodels/MpNc/energy/MpNcvolumevariablesenergy.hh
index 8d0bba1ecf..5da3e6ab9d 100644
--- a/dumux/boxmodels/MpNc/energy/MpNcvolumevariablesenergy.hh
+++ b/dumux/boxmodels/MpNc/energy/MpNcvolumevariablesenergy.hh
@@ -27,6 +27,8 @@
 #ifndef DUMUX_MPNC_ENERGY_VOLUME_VARIABLES_HH
 #define DUMUX_MPNC_ENERGY_VOLUME_VARIABLES_HH
 
+#include <dumux/material/MpNcfluidstates/equilibriumfluidstate.hh>
+
 namespace Dumux
 {
 /*!
@@ -57,19 +59,24 @@ class MPNCVolumeVariablesEnergy
     enum { numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)) };
     //typedef typename GET_PROP_TYPE(TypeTag, PTAG(MPNCEnergyIndices)) EnergyIndices;
 
-
+    typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+    typedef typename FluidSystem::ParameterCache ParameterCache; 
+    typedef Dumux::EquilibriumFluidState<Scalar, FluidSystem> FluidState;
+    
 public:
     /*!
      * \brief Update the temperature of the sub-control volume.
      */
-    Scalar getTemperature(const PrimaryVariables &sol,
-                          const Element &element,
-                          const FVElementGeometry &elemGeom,
-                          int scvIdx,
-                          const Problem &problem,
-                          const int temperatureIdx) const
+    void updateTemperatures(FluidState &fs,
+                            ParameterCache &paramCache,
+                            const PrimaryVariables &sol,
+                            const Element &element,
+                            const FVElementGeometry &elemGeom,
+                            int scvIdx,
+                            const Problem &problem) const
     {
-        return problem.boxTemperature(element, elemGeom, scvIdx);
+        Scalar T = problem.boxTemperature(element, elemGeom, scvIdx);
+        fs.setTemperature(T);
     }
 
 
@@ -79,39 +86,23 @@ public:
      *
      * Since we are isothermal, we don't need to do anything!
      */
-    template <class MutableParams>
-    void update(MutableParams &mutParams,
+    void update(FluidState &fs,
+                ParameterCache &paramCache,
                 const PrimaryVariables &sol,
                 const Element &element,
                 const FVElementGeometry &elemGeom,
                 int scvIdx,
                 const Problem &problem)
     {
-        temperature_ = problem.boxTemperature(element, elemGeom, scvIdx);
-
-        // set the fluid temperatures
-        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
-            mutParams.setTemperature(phaseIdx, temperature_);
     }
 
-    /*!
-     * \brief Return the temperature of any given phase [K]
-     */
-    Scalar temperature(int phaseIdx = 0) const
-    { return temperature_; }
-
     /*!
      * \brief If running under valgrind this produces an error message
      *        if some of the object's attributes is undefined.
      */
     void checkDefined() const
     {
-        Valgrind::CheckDefined(temperature_);
     }
-
-protected:
-    Scalar temperature_;
-
 };
 
 /*!
@@ -130,7 +121,6 @@ class MPNCVolumeVariablesEnergy<TypeTag, /*enableEnergy=*/true, /*kineticEnergyT
     typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
     typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
 
-    typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
     typedef typename GET_PROP_TYPE(TypeTag, PTAG(MPNCIndices)) Indices;
 
     enum { numPhases        = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)) };
@@ -139,28 +129,33 @@ class MPNCVolumeVariablesEnergy<TypeTag, /*enableEnergy=*/true, /*kineticEnergyT
     enum { numEnergyEqs     = Indices::NumPrimaryEnergyVars};
     enum { temperature0Idx = Indices::temperatureIdx };
 
+    typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+    typedef typename FluidSystem::ParameterCache ParameterCache; 
+    typedef Dumux::EquilibriumFluidState<Scalar, FluidSystem> FluidState;
+    
 public:
     /*!
      * \brief Update the temperature of the sub-control volume.
      */
-    Scalar getTemperature(const PrimaryVariables &sol,
-                          const Element &element,
-                          const FVElementGeometry &elemGeom,
-                          int scvIdx,
-                          const Problem &problem,
-                          const int dummy) const
+    void updateTemperatures(FluidState &fs,
+                            ParameterCache &paramCache,
+                            const PrimaryVariables &sol,
+                            const Element &element,
+                            const FVElementGeometry &elemGeom,
+                            int scvIdx,
+                            const Problem &problem) const
     {
         // retrieve temperature from solution vector
-        return sol[temperatureIdx];
+        Scalar T = sol[temperatureIdx];
+        fs.setTemperature(T);
     }
 
     /*!
      * \brief Update the enthalpy and the internal energy for a given
      *        control volume.
      */
-    template <class MutableParams>
-    void update(MutableParams &mutParams,
-                const PrimaryVariables &sol,
+    void update(FluidState &fs,
+                ParameterCache &paramCache,
                 const Element &element,
                 const FVElementGeometry &elemGeom,
                 int scvIdx,
@@ -168,24 +163,20 @@ public:
     {
         Valgrind::SetUndefined(*this);
 
-        // set the fluid temperatures
-        temperature_ = sol[temperature0Idx];
-        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
-            mutParams.setTemperature(phaseIdx, temperature_);
-
+        // heat capacities of the fluids plus the porous medium
         heatCapacity_ =
             problem.spatialParameters().heatCapacity(element, elemGeom, scvIdx);
         Valgrind::CheckDefined(heatCapacity_);
 
         soilDensity_ =
-                problem.spatialParameters().soilDensity(element, elemGeom, scvIdx);
+            problem.spatialParameters().soilDensity(element, elemGeom, scvIdx);
         Valgrind::CheckDefined(soilDensity_);
 
         // set the enthalpies
         for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
-            Scalar h =
-                FluidSystem::computeEnthalpy(mutParams, phaseIdx);
-            mutParams.setEnthalpy(phaseIdx, h);
+            Scalar u =
+                FluidSystem::internalEnergy(fs, paramCache, phaseIdx);
+            fs.setInternalEnergy(phaseIdx, u);
         }
     }
 
@@ -196,13 +187,6 @@ public:
     Scalar heatCapacity() const
     { return heatCapacity_; };
 
-    /*!
-     * \brief Returns the temperature in fluid / solid phase(s)
-     *        the sub-control volume.
-     */
-    Scalar temperature(int phaseIdx = 0) const
-    { return temperature_; }
-
     /*!
      * \brief Returns the total density of the given soil [kg / m^3] in
      *        the sub-control volume.
@@ -218,13 +202,11 @@ public:
     {
         Valgrind::CheckDefined(heatCapacity_);
         Valgrind::CheckDefined(soilDensity_);
-        Valgrind::CheckDefined(temperature_);
     };
 
 protected:
     Scalar heatCapacity_;
     Scalar soilDensity_;
-    Scalar temperature_;
 };
 
 } // end namepace
diff --git a/dumux/boxmodels/MpNc/mass/MpNcvolumevariablesmass.hh b/dumux/boxmodels/MpNc/mass/MpNcvolumevariablesmass.hh
index 3887239dbf..d67f6574e1 100644
--- a/dumux/boxmodels/MpNc/mass/MpNcvolumevariablesmass.hh
+++ b/dumux/boxmodels/MpNc/mass/MpNcvolumevariablesmass.hh
@@ -60,7 +60,8 @@ class MPNCVolumeVariablesMass
     enum { fug0Idx = Indices::fug0Idx };
 
     typedef Dune::FieldVector<Scalar, numComponents> ComponentVector;
-
+    
+    typedef typename FluidSystem::ParameterCache ParameterCache;
 public:
     /*!
      * \brief The fluid state which is used by the volume variables to
@@ -75,11 +76,10 @@ public:
      * \brief Update composition of all phases in the mutable
      *        parameters from the primary variables.
      */
-    template <class MutableParams>
-    void update(MutableParams &mutParams,
+    void update(FluidState &fs,
+                ParameterCache &paramCache,
                 const PrimaryVariables &priVars,
                 const VolumeVariables *hint,
-
                 const Problem &problem,
                 const Element &element,
                 const FVElementGeometry &elemGeom,
@@ -97,10 +97,10 @@ public:
                 Scalar x_ij = 1.0/numComponents;
                 if (hint)
                     // use the hint for the initial mole fraction!
-                    x_ij = hint->fluidState().moleFrac(phaseIdx, compIdx);
+                    x_ij = hint->fluidState().moleFraction(phaseIdx, compIdx);
 
                 // set initial guess of the component's mole fraction
-                mutParams.setMoleFrac(phaseIdx,
+                fs.setMoleFraction(phaseIdx,
                                       compIdx,
                                       x_ij);
 
@@ -111,18 +111,18 @@ public:
             if (!hint)
                 // if no hint was given, we ask the constraint solver
                 // to make an initial guess
-                CompositionFromFugacitiesSolver::guessInitial(mutParams, phaseIdx, fug);
-            CompositionFromFugacitiesSolver::solve(mutParams, phaseIdx, fug);
+                CompositionFromFugacitiesSolver::guessInitial(fs, paramCache, phaseIdx, fug);
+            CompositionFromFugacitiesSolver::solve(fs, paramCache, phaseIdx, fug);
 
             /*
               std::cout << "final composition: " << FluidSystem::phaseName(phaseIdx) << "="
-              << mutParams.moleFrac(phaseIdx, 0) << " "
-              << mutParams.moleFrac(phaseIdx, 1) << " "
-              << mutParams.moleFrac(phaseIdx, 2) << " "
-              << mutParams.moleFrac(phaseIdx, 3) << " "
-              << mutParams.moleFrac(phaseIdx, 4) << " "
-              << mutParams.moleFrac(phaseIdx, 5) << " "
-              << mutParams.moleFrac(phaseIdx, 6) << "\n";
+              << fs.moleFrac(phaseIdx, 0) << " "
+              << fs.moleFrac(phaseIdx, 1) << " "
+              << fs.moleFrac(phaseIdx, 2) << " "
+              << fs.moleFrac(phaseIdx, 3) << " "
+              << fs.moleFrac(phaseIdx, 4) << " "
+              << fs.moleFrac(phaseIdx, 5) << " "
+              << fs.moleFrac(phaseIdx, 6) << "\n";
             */
 
         }
diff --git a/dumux/material/MpNcconstraintsolvers/compositionfromfugacities.hh b/dumux/material/MpNcconstraintsolvers/compositionfromfugacities.hh
index a3fde6d5c8..1b40b35484 100644
--- a/dumux/material/MpNcconstraintsolvers/compositionfromfugacities.hh
+++ b/dumux/material/MpNcconstraintsolvers/compositionfromfugacities.hh
@@ -38,18 +38,22 @@ namespace Dumux {
 template <class Scalar, class FluidSystem>
 class CompositionFromFugacities
 {
-    typedef typename FluidSystem::MutableParameters MutableParameters;
-
     enum { numPhases = FluidSystem::numPhases };
     enum { numComponents = FluidSystem::numComponents };
-    
+
+    typedef typename FluidSystem::ParameterCache ParameterCache;
+
 public:
     typedef Dune::FieldVector<Scalar, numComponents> ComponentVector;
 
     /*!
      * \brief Guess an initial value for the composition of the phase.
      */
-    static void guessInitial(MutableParameters &mutParams, int phaseIdx, const ComponentVector &fugVec)
+    template <class FluidState>
+    static void guessInitial(FluidState &fluidState,
+                             ParameterCache &paramCache,
+                             int phaseIdx,
+                             const ComponentVector &fugVec)
     {
         if (FluidSystem::isIdealMixture(phaseIdx)) 
             return;
@@ -57,9 +61,9 @@ public:
         // Pure component fugacities
         for (int i = 0; i < numComponents; ++ i) {
             //std::cout << f << " -> " << mutParams.fugacity(phaseIdx, i)/f << "\n";
-            mutParams.setMoleFrac(phaseIdx,
-                                  i, 
-                                  1.0/numComponents);
+            fluidState.setMoleFraction(phaseIdx,
+                                   i, 
+                                   1.0/numComponents);
         }
     };
 
@@ -69,14 +73,16 @@ public:
      *
      * The phase's fugacities must already be set.
      */
-    static void solve(MutableParameters &mutParams, 
+    template <class FluidState>
+    static void solve(FluidState &fluidState,
+                      ParameterCache &paramCache,
                       int phaseIdx, 
                       const ComponentVector &targetFug)
     {
         // use a much more efficient method in case the phase is an
         // ideal mixture
         if (FluidSystem::isIdealMixture(phaseIdx)) {
-            solveIdealMix_(mutParams, phaseIdx, targetFug);
+            solveIdealMix_(fluidState, paramCache, phaseIdx, targetFug);
             return;
         }
 
@@ -85,7 +91,7 @@ public:
         // save initial composition in case something goes wrong
         Dune::FieldVector<Scalar, numComponents> xInit;
         for (int i = 0; i < numComponents; ++i) {
-            xInit[i] = mutParams.moleFrac(phaseIdx, i);
+            xInit[i] = fluidState.moleFraction(phaseIdx, i);
         };
         
         /////////////////////////
@@ -99,22 +105,25 @@ public:
         // right hand side
         Dune::FieldVector<Scalar, numComponents> b;
         
+        fluidState.updateAverageMolarMass(phaseIdx);
+        paramCache.updatePhase(fluidState, phaseIdx);
+        
         // maximum number of iterations
         const int nMax = 25;
         for (int nIdx = 0; nIdx < nMax; ++nIdx) {
             // calculate Jacobian matrix and right hand side
-            linearize_(J, b, mutParams, phaseIdx, targetFug); 
+            linearize_(J, b, fluidState, paramCache, phaseIdx, targetFug); 
             Valgrind::CheckDefined(J);
             Valgrind::CheckDefined(b);
 
             /*
             std::cout << FluidSystem::phaseName(phaseIdx) << "Phase composition: ";
             for (int i = 0; i < FluidSystem::numComponents; ++i)
-                std::cout << mutParams.moleFrac(phaseIdx, i) << " ";
+                std::cout << fluidState.moleFraction(phaseIdx, i) << " ";
             std::cout << "\n";
             std::cout << FluidSystem::phaseName(phaseIdx) << "Phase phi: ";
             for (int i = 0; i < FluidSystem::numComponents; ++i)
-                std::cout << mutParams.fugacityCoeff(phaseIdx, i) << " ";
+                std::cout << fluidState.fugacityCoeff(phaseIdx, i) << " ";
             std::cout << "\n";
             */
 
@@ -131,10 +140,10 @@ public:
             /*
             std::cout << FluidSystem::phaseName(phaseIdx) << "Phase composition: ";
             for (int i = 0; i < FluidSystem::numComponents; ++i)
-                std::cout << mutParams.moleFrac(phaseIdx, i) << " ";
+                std::cout << fluidState.moleFraction(phaseIdx, i) << " ";
             std::cout << "\n";
             std::cout << "J: " << J << "\n";
-            std::cout << "rho: " << mutParams.density(phaseIdx) << "\n";
+            std::cout << "rho: " << fluidState.density(phaseIdx) << "\n";
             std::cout << "delta: " << x << "\n";
             std::cout << "defect: " << b << "\n";
 
@@ -145,14 +154,12 @@ public:
 
             // update the fluid composition. b is also used to store
             // the defect for the next iteration.
-            Scalar relError = update_(mutParams, x, b, phaseIdx, targetFug);
+            Scalar relError = update_(fluidState, paramCache, x, b, phaseIdx, targetFug);
             //std::cout << "relError: " << relError << "\n";
             
             if (relError < 1e-9) {
-                mutParams.updateMeanMolarMass(phaseIdx);
-                
-                Scalar Vm = FluidSystem::computeMolarVolume(mutParams, phaseIdx);
-                mutParams.setMolarVolume(phaseIdx, Vm);
+                Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
+                fluidState.setDensity(phaseIdx, rho);
 
                 //std::cout << "num iterations: " << nIdx << "\n";
                 return;
@@ -163,8 +170,8 @@ public:
                    "Calculating the " << FluidSystem::phaseName(phaseIdx)
                    << "Phase composition failed. Initial {x} = {"
                    << xInit 
-                   << "}, {fug_t} = {" << targetFug << "}, p = " << mutParams.pressure(phaseIdx)
-                   << ", T = " << mutParams.temperature(phaseIdx));
+                   << "}, {fug_t} = {" << targetFug << "}, p = " << fluidState.pressure(phaseIdx)
+                   << ", T = " << fluidState.temperature(phaseIdx));
     };
 
 
@@ -172,26 +179,35 @@ protected:
     // update the phase composition in case the phase is an ideal
     // mixture, i.e. the component's fugacity coefficients are
     // independent of the phase's composition.
-    static void solveIdealMix_(MutableParameters &mutParams, 
+    template <class FluidState>
+    static void solveIdealMix_(FluidState &fluidState, 
+                               ParameterCache &paramCache,
                                int phaseIdx,
                                const ComponentVector &fugacities)
     {        
         for (int i = 0; i < numComponents; ++ i) {
-            Scalar phi = FluidSystem::computeFugacityCoeff(mutParams, phaseIdx, i);
-            Scalar gamma = phi * mutParams.pressure(phaseIdx);
-            mutParams.setFugacityCoeff(phaseIdx, i, phi);
-            mutParams.setMoleFrac(phaseIdx, i, fugacities[i]/gamma);
+            Scalar phi = FluidSystem::fugacityCoefficient(fluidState,
+                                                          paramCache,
+                                                          phaseIdx,
+                                                          i);
+            Scalar gamma = phi * fluidState.pressure(phaseIdx);
+            fluidState.setFugacityCoefficient(phaseIdx, i, phi);
+            fluidState.setMoleFraction(phaseIdx, i, fugacities[i]/gamma);
         };
-        mutParams.updateMeanMolarMass(phaseIdx);
         
-        Scalar Vm = FluidSystem::computeMolarVolume(mutParams, phaseIdx);
-        mutParams.setMolarVolume(phaseIdx, Vm);
+        fluidState.updateAverageMolarMass(phaseIdx);
+        paramCache.updatePhase(fluidState, phaseIdx);
+        
+        Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
+        fluidState.setDensity(phaseIdx, rho);
         return;
     };
 
+    template <class FluidState>
     static void linearize_(Dune::FieldMatrix<Scalar, numComponents, numComponents> &J,
                            Dune::FieldVector<Scalar, numComponents> &defect,
-                           MutableParameters &mutParams, 
+                           FluidState &fluidState, 
+                           ParameterCache &paramCache,
                            int phaseIdx,
                            const ComponentVector &targetFug)
     {       
@@ -201,11 +217,12 @@ protected:
         // calculate the defect (deviation of the current fugacities
         // from the target fugacities)
         for (int i = 0; i < numComponents; ++ i) {
-            Scalar phi = FluidSystem::computeFugacityCoeff(mutParams,
-                                                           phaseIdx,
-                                                           i);
-            Scalar f = phi*mutParams.pressure(phaseIdx)*mutParams.moleFrac(phaseIdx, i);
-            mutParams.setFugacityCoeff(phaseIdx, i, phi);
+            Scalar phi = FluidSystem::fugacityCoefficient(fluidState,
+                                                          paramCache,
+                                                          phaseIdx,
+                                                          i);
+            Scalar f = phi*fluidState.pressure(phaseIdx)*fluidState.moleFraction(phaseIdx, i);
+            fluidState.setFugacityCoefficient(phaseIdx, i, phi);
                 
             defect[i] = targetFug[i] - f;
         }
@@ -221,21 +238,24 @@ protected:
             // forward differences
 
             // deviate the mole fraction of the i-th component
-            Scalar x_i = mutParams.moleFrac(phaseIdx, i);
-            mutParams.setMoleFrac(phaseIdx, i, x_i + eps);
+            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
             // fugacities
             for (int j = 0; j < numComponents; ++j) {
                 // compute the j-th component's fugacity coefficient ...
-                Scalar phi = FluidSystem::computeFugacityCoeff(mutParams,
-                                                               phaseIdx,
-                                                               j);
+                Scalar phi = FluidSystem::fugacityCoefficient(fluidState,
+                                                              paramCache,
+                                                              phaseIdx,
+                                                              j);
                 // ... and its fugacity ...
                 Scalar f = 
                     phi *
-                    mutParams.pressure(phaseIdx) *
-                    mutParams.moleFrac(phaseIdx, j);
+                    fluidState.pressure(phaseIdx) *
+                    fluidState.moleFraction(phaseIdx, j);
                 // as well as the defect for this fugacity
                 Scalar defJPlusEps = targetFug[j] - f;
 
@@ -245,14 +265,18 @@ protected:
             }
             
             // reset composition to original value
-            mutParams.setMoleFrac(phaseIdx, i, x_i);
-            
+            fluidState.setMoleFraction(phaseIdx, i, x_i);
+            fluidState.updateAverageMolarMass(phaseIdx);
+            paramCache.updatePhaseSingleMoleFraction(fluidState, phaseIdx, i);
+
             // end forward differences
             ////////
         }
     }
 
-    static Scalar update_(MutableParameters &mutParams,
+    template <class FluidState>
+    static Scalar update_(FluidState &fluidState,
+                          ParameterCache &paramCache,
                           Dune::FieldVector<Scalar, numComponents> &x,
                           Dune::FieldVector<Scalar, numComponents> &b,
                           int phaseIdx,
@@ -264,10 +288,10 @@ protected:
         Scalar sumDelta = 0;
         Scalar sumx = 0;
         for (int i = 0; i < numComponents; ++i) {
-            origComp[i] = mutParams.moleFrac(phaseIdx, i);
+            origComp[i] = fluidState.moleFraction(phaseIdx, i);
             relError = std::max(relError, std::abs(x[i]));
             
-            sumx += std::abs(mutParams.moleFrac(phaseIdx, i));
+            sumx += std::abs(fluidState.moleFraction(phaseIdx, i));
             sumDelta += std::abs(x[i]);
         };
 
@@ -278,7 +302,7 @@ protected:
             x /= (sumDelta/maxDelta);
 #endif
 
-        //Scalar curDefect = calculateDefect_(mutParams, phaseIdx, targetFug);
+        //Scalar curDefect = calculateDefect_(fluidState, phaseIdx, targetFug);
         //Scalar nextDefect;
         //Scalar sumMoleFrac = 0.0;
         //ComponentVector newB(1e100);
@@ -297,16 +321,19 @@ protected:
                 else
                     newx = 0;
 #endif
-                mutParams.setMoleFrac(phaseIdx, i, newx);
+                fluidState.setMoleFraction(phaseIdx, i, newx);
                 //sumMoleFrac += std::abs(newx);
             }
 
+            fluidState.updateAverageMolarMass(phaseIdx);
+            paramCache.updatePhaseComposition(fluidState, phaseIdx);
+
             /*
             // if the sum of the mole fractions gets 0, we take the
             // original composition divided by 100
             if (sumMoleFrac < 1e-10) {
                 for (int i = 0; i < numComponents; ++i) {
-                    mutParams.setMoleFrac(phaseIdx, i, origComp[i]/100);
+                    fluidState.setMoleFraction(phaseIdx, i, origComp[i]/100);
                 }
                 return relError;
             }
@@ -315,13 +342,13 @@ protected:
             /*
             // calculate new residual
             for (int i = 0; i < numComponents; ++i) {
-                Scalar phi = FluidSystem::computeFugacityCoeff(mutParams,
+                Scalar phi = FluidSystem::computeFugacityCoeff(fluidState,
                                                                phaseIdx,
                                                                i);
-                mutParams.setFugacityCoeff(phaseIdx, i, phi);
+                fluidState.setFugacityCoeff(phaseIdx, i, phi);
             }
 
-            nextDefect = calculateDefect_(mutParams, phaseIdx, targetFug);
+            nextDefect = calculateDefect_(fluidState, phaseIdx, targetFug);
             //std::cout << "try delta: " << x << "\n";
             //std::cout << "defect: old=" << curDefect << " new=" << nextDefect << "\n"; 
             if (nextDefect <= curDefect)
@@ -335,7 +362,8 @@ protected:
         return relError;
     }
 
-    static Scalar calculateDefect_(const MutableParameters &params, 
+    template <class FluidState>
+    static Scalar calculateDefect_(const FluidState &params, 
                                    int phaseIdx, 
                                    const ComponentVector &targetFug)
     {
diff --git a/dumux/material/MpNcfluidstates/equilibriumfluidstate.hh b/dumux/material/MpNcfluidstates/equilibriumfluidstate.hh
index 5f2e2d611c..d5f5eff514 100644
--- a/dumux/material/MpNcfluidstates/equilibriumfluidstate.hh
+++ b/dumux/material/MpNcfluidstates/equilibriumfluidstate.hh
@@ -36,12 +36,12 @@ namespace Dumux
  *        multi-phase, multi-component fluid system assuming
  *        thermodynamic equilibrium.
  */
-template <class Scalar, class StaticParameters>
+template <class Scalar, class FluidSystem>
 class EquilibriumFluidState
 {
 public:
-    enum { numComponents = StaticParameters::numComponents };
-    enum { numPhases = StaticParameters::numPhases };
+    enum { numPhases = FluidSystem::numPhases };
+    enum { numComponents = FluidSystem::numComponents };
 
     EquilibriumFluidState()
     { Valgrind::SetUndefined(*this); }
@@ -63,18 +63,18 @@ public:
     /*!
      * \brief The mole fraction of a component in a phase []
      */
-    Scalar moleFrac(int phaseIdx, int compIdx) const
-    { return moleFrac_[phaseIdx][compIdx]; }
+    Scalar moleFraction(int phaseIdx, int compIdx) const
+    { return moleFraction_[phaseIdx][compIdx]; }
 
     /*!
      * \brief The mass fraction of a component in a phase []
      */
-    Scalar massFrac(int phaseIdx, int compIdx) const
+    Scalar massFraction(int phaseIdx, int compIdx) const
     { 
         return
-            sumMassFrac(phaseIdx)*
+            sumMassFractions(phaseIdx)*
             molarity(phaseIdx, compIdx)*
-            StaticParameters::molarMass(compIdx)
+            FluidSystem::molarMass(compIdx)
             /
             density(phaseIdx);
     }
@@ -84,22 +84,22 @@ public:
      *
      * We define this to be the same as the sum of all mass fractions.
      */
-    Scalar sumMoleFrac(int phaseIdx) const
-    { return sumMoleFrac_[phaseIdx]; }
+    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 sumMassFrac(int phaseIdx) const
-    { return sumMoleFrac_[phaseIdx]; }
+    Scalar sumMassFractions(int phaseIdx) const
+    { return sumMoleFractions_[phaseIdx]; }
 
     /*!
-     * \brief The mean molar mass of a fluid phase [kg/mol]
+     * \brief The average molar mass of a fluid phase [kg/mol]
      */
-    Scalar meanMolarMass(int phaseIdx) const
-    { return meanMolarMass_[phaseIdx]; }
+    Scalar averageMolarMass(int phaseIdx) const
+    { return averageMolarMass_[phaseIdx]; }
 
     /*!
      * \brief The concentration of a component in a phase [mol/m^3]
@@ -111,37 +111,37 @@ public:
      * http://en.wikipedia.org/wiki/Concentration
      */
     Scalar molarity(int phaseIdx, int compIdx) const
-    { return molarDensity(phaseIdx)*moleFrac(phaseIdx, compIdx); }
+    { return molarDensity(phaseIdx)*moleFraction(phaseIdx, compIdx); }
 
     /*!
      * \brief The fugacity of a component in a phase [Pa]
      */
     Scalar fugacity(int phaseIdx, int compIdx) const
-    { return fugacityCoeff(phaseIdx, compIdx)*moleFrac(phaseIdx, compIdx)*pressure(phaseIdx); }
+    { return fugacityCoefficient(phaseIdx, compIdx)*moleFraction(phaseIdx, compIdx)*pressure(phaseIdx); }
 
     /*!
      * \brief The fugacity coefficient of a component in a phase [Pa]
      */
-    Scalar fugacityCoeff(int phaseIdx, int compIdx) const
-    { return fugacityCoeff_[phaseIdx][compIdx]; }
+    Scalar fugacityCoefficient(int phaseIdx, int compIdx) const
+    { return fugacityCoefficient_[phaseIdx][compIdx]; }
 
     /*!
      * \brief The molar volume of a fluid phase [m^3/mol]
      */
     Scalar molarVolume(int phaseIdx) const
-    { return molarVolume_[phaseIdx]; }
+    { return 1/molarDensity(phaseIdx); }
 
     /*!
      * \brief The mass density of a fluid phase [kg/m^3]
      */
     Scalar density(int phaseIdx) const
-    { return molarDensity(phaseIdx)*meanMolarMass(phaseIdx); }
+    { return density_[phaseIdx]; }
 
     /*!
      * \brief The molar density of a fluid phase [mol/m^3]
      */
     Scalar molarDensity(int phaseIdx) const
-    { return 1/molarVolume(phaseIdx); }
+    { return density_[phaseIdx]/averageMolarMass(phaseIdx); }
 
     /*!
      * \brief The temperature of a fluid phase [K]
@@ -159,13 +159,13 @@ public:
      * \brief The specific enthalpy of a fluid phase [J/kg]
      */
     Scalar enthalpy(int phaseIdx) const
-    { return enthalpy_[phaseIdx]; }
+    { return internalEnergy_[phaseIdx] + pressure(phaseIdx)/(density(phaseIdx)); }
 
     /*!
      * \brief The specific internal energy of a fluid phase [J/kg]
      */
     Scalar internalEnergy(int phaseIdx) const
-    { return enthalpy(phaseIdx) - pressure(phaseIdx)/(density(phaseIdx)); }
+    { return internalEnergy_[phaseIdx]; }
 
     /*!
      * \brief The dynamic viscosity of a fluid phase [Pa s]
@@ -223,14 +223,14 @@ public:
     {
         for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
             for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
-                moleFrac_[phaseIdx][compIdx] = fs.moleFrac(phaseIdx, compIdx);
-                fugacityCoeff_[phaseIdx][compIdx] = fs.fugacityCoeff(phaseIdx, compIdx);
+                moleFraction_[phaseIdx][compIdx] = fs.moleFraction(phaseIdx, compIdx);
+                fugacityCoefficient_[phaseIdx][compIdx] = fs.fugacityCoefficient(phaseIdx, compIdx);
             }
-            meanMolarMass_[phaseIdx] = fs.meanMolarMass(phaseIdx);
+            averageMolarMass_[phaseIdx] = fs.averageMolarMass(phaseIdx);
             pressure_[phaseIdx] = fs.pressure(phaseIdx);
             saturation_[phaseIdx] = fs.saturation(phaseIdx);
-            molarVolume_[phaseIdx] = fs.molarVolume(phaseIdx);
-            enthalpy_[phaseIdx] = fs.enthalpy(phaseIdx);
+            density_[phaseIdx] = fs.density(phaseIdx);
+            internalEnergy_[phaseIdx] = fs.internalEnergy(phaseIdx);
             viscosity_[phaseIdx] = fs.viscosity(phaseIdx);
         }
         temperature_ = fs.temperature(0);
@@ -257,26 +257,26 @@ public:
     /*!
      * \brief Set the mole fraction of a component in a phase []
      */
-    void setMoleFrac(int phaseIdx, int compIdx, Scalar value)
-    { moleFrac_[phaseIdx][compIdx] = value; }   
+    void setMoleFraction(int phaseIdx, int compIdx, Scalar value)
+    { moleFraction_[phaseIdx][compIdx] = value; }   
 
     /*!
      * \brief Set the fugacity of a component in a phase []
      */
-    void setFugacityCoeff(int phaseIdx, int compIdx, Scalar value)
-    { fugacityCoeff_[phaseIdx][compIdx] = value; }   
+    void setFugacityCoefficient(int phaseIdx, int compIdx, Scalar value)
+    { fugacityCoefficient_[phaseIdx][compIdx] = value; }   
 
     /*!
-     * \brief Set the molar volume of a phase [m^3/mol]
+     * \brief Set the density of a phase [kg / m^3]
      */
-    void setMolarVolume(int phaseIdx, Scalar value)
-    { molarVolume_[phaseIdx] = value; }   
+    void setDensity(int phaseIdx, Scalar value)
+    { density_[phaseIdx] = value; }   
 
     /*!
      * \brief Set the specific enthalpy of a phase [J/m^3]
      */
-    void setEnthalpy(int phaseIdx, Scalar value)
-    { enthalpy_[phaseIdx] = value; }   
+    void setInternalEnergy(int phaseIdx, Scalar value)
+    { internalEnergy_[phaseIdx] = value; }   
 
     /*!
      * \brief Set the dynamic viscosity of a phase [Pa s]
@@ -288,18 +288,18 @@ public:
      * \brief Calculatate the mean molar mass of a phase given that
      *        all mole fractions have been set
      */
-    void updateMeanMolarMass(int phaseIdx)
+    void updateAverageMolarMass(int phaseIdx)
     {
-        meanMolarMass_[phaseIdx] = 0;
-        sumMoleFrac_[phaseIdx] = 0;
+        averageMolarMass_[phaseIdx] = 0;
+        sumMoleFractions_[phaseIdx] = 0;
      
         for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
-            sumMoleFrac_[phaseIdx] += moleFrac_[phaseIdx][compIdx];
-            meanMolarMass_[phaseIdx] += moleFrac_[phaseIdx][compIdx]*StaticParameters::molarMass(compIdx);
+            sumMoleFractions_[phaseIdx] += moleFraction_[phaseIdx][compIdx];
+            averageMolarMass_[phaseIdx] += moleFraction_[phaseIdx][compIdx]*FluidSystem::molarMass(compIdx);
         }
-        sumMoleFrac_[phaseIdx] = std::max(1e-15, sumMoleFrac_[phaseIdx]);
-        Valgrind::CheckDefined(meanMolarMass_[phaseIdx]);
-        Valgrind::CheckDefined(sumMoleFrac_[phaseIdx]);
+        sumMoleFractions_[phaseIdx] = std::max(1e-15, sumMoleFractions_[phaseIdx]);
+        Valgrind::CheckDefined(averageMolarMass_[phaseIdx]);
+        Valgrind::CheckDefined(sumMoleFractions_[phaseIdx]);
     }
     
     /*!
@@ -313,17 +313,17 @@ public:
     void checkDefined() const
     {
 #if HAVE_VALGRIND && ! defined NDEBUG
-        for (int i = 0; i < numPhases; ++i) {
-            for (int j = 0; j < numComponents; ++j) {
-                Valgrind::CheckDefined(moleFrac_[i][j]);
-                Valgrind::CheckDefined(fugacityCoeff_[i][j]);
+        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+            for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
+                Valgrind::CheckDefined(moleFraction_[phaseIdx][compIdx]);
+                Valgrind::CheckDefined(fugacityCoefficient_[phaseIdx][compIdx]);
             }
-            Valgrind::CheckDefined(meanMolarMass_[i]);
-            Valgrind::CheckDefined(pressure_[i]);
-            Valgrind::CheckDefined(saturation_[i]);
-            Valgrind::CheckDefined(molarVolume_[i]);
-            //Valgrind::CheckDefined(enthalpy_[i]);
-            Valgrind::CheckDefined(viscosity_[i]);
+            Valgrind::CheckDefined(averageMolarMass_[phaseIdx]);
+            Valgrind::CheckDefined(pressure_[phaseIdx]);
+            Valgrind::CheckDefined(saturation_[phaseIdx]);
+            Valgrind::CheckDefined(density_[phaseIdx]);
+            //Valgrind::CheckDefined(enthalpy_[phaseIdx]);
+            Valgrind::CheckDefined(viscosity_[phaseIdx]);
         }
 
         Valgrind::CheckDefined(temperature_);
@@ -331,15 +331,15 @@ public:
     }
 
 protected:
-    Scalar moleFrac_[numPhases][numComponents];
-    Scalar fugacityCoeff_[numPhases][numComponents];
+    Scalar moleFraction_[numPhases][numComponents];
+    Scalar fugacityCoefficient_[numPhases][numComponents];
 
-    Scalar meanMolarMass_[numPhases];
-    Scalar sumMoleFrac_[numPhases];
+    Scalar averageMolarMass_[numPhases];
+    Scalar sumMoleFractions_[numPhases];
     Scalar pressure_[numPhases];
     Scalar saturation_[numPhases];
-    Scalar molarVolume_[numPhases];
-    Scalar enthalpy_[numPhases];
+    Scalar density_[numPhases];
+    Scalar internalEnergy_[numPhases];
     Scalar viscosity_[numPhases];
     Scalar temperature_;
 };
diff --git a/dumux/material/MpNcfluidstates/genericfluidstate.hh b/dumux/material/MpNcfluidstates/genericfluidstate.hh
index 16a00e2038..d19bfce8df 100644
--- a/dumux/material/MpNcfluidstates/genericfluidstate.hh
+++ b/dumux/material/MpNcfluidstates/genericfluidstate.hh
@@ -39,12 +39,12 @@ namespace Dumux
  *        multi-phase, multi-component fluid system without using
  *        any assumptions.
  */
-template <class Scalar, class StaticParameters>
+template <class Scalar, class FluidSystem>
 class GenericFluidState
 {
 public:
-    enum { numComponents    = StaticParameters::numComponents };
-    enum { numPhases        = StaticParameters::numPhases };
+    enum { numPhases = FluidSystem::numPhases };
+    enum { numComponents = FluidSystem::numComponents };
 
     GenericFluidState()
     { Valgrind::SetUndefined(*this); }
@@ -62,20 +62,20 @@ public:
     /*!
      * \brief The mole fraction of a component in a phase []
      */
-    Scalar moleFrac(int phaseIdx, int compIdx) const
-    { return moleFrac_[phaseIdx][compIdx]; }
+    Scalar moleFraction(int phaseIdx, int compIdx) const
+    { return moleFraction_[phaseIdx][compIdx]; }
 
 
     /*!
      * \brief The mass fraction of a component in a phase []
      */
-    Scalar massFrac(int phaseIdx, int compIdx) const
+    Scalar massFraction(int phaseIdx, int compIdx) const
     { 
         return
-            sumMassFrac(phaseIdx) *
-            moleFrac_[phaseIdx][compIdx] *
-            StaticParameters::molarMass(compIdx) /
-            meanMolarMass_[phaseIdx];
+            sumMassFraction(phaseIdx)
+            * moleFraction_[phaseIdx][compIdx]
+            * FluidSystem::molarMass(compIdx)
+            / averageMolarMass_[phaseIdx];
     }
 
     /*!
@@ -83,22 +83,22 @@ public:
      *
      * We define this to be the same as the sum of all mass fractions.
      */
-    Scalar sumMoleFrac(int phaseIdx) const
-    { return sumMoleFrac_[phaseIdx]; }
+    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 sumMassFrac(int phaseIdx) const
-    { return sumMoleFrac_[phaseIdx]; }
+    Scalar sumMassFraction(int phaseIdx) const
+    { return sumMoleFractions_[phaseIdx]; }
 
     /*!
      * \brief The mean molar mass of a fluid phase [kg/mol]
      */
-    Scalar meanMolarMass(int phaseIdx) const
-    { return meanMolarMass_[phaseIdx]; }
+    Scalar averageMolarMass(int phaseIdx) const
+    { return averageMolarMass_[phaseIdx]; }
 
     /*!
      * \brief The molar concentration of a component in a phase [mol/m^3]
@@ -109,40 +109,38 @@ public:
      * http://en.wikipedia.org/wiki/Concentration
      */
     Scalar molarity(int phaseIdx, int compIdx) const
-    { return molarDensity(phaseIdx)*moleFrac(phaseIdx, compIdx); }
+    { return molarDensity(phaseIdx)*moleFraction(phaseIdx, compIdx); }
 
     /*!
      * \brief The fugacity coefficient of a component in a phase []
      */
-    Scalar fugacityCoeff(int phaseIdx, int compIdx) const
-    { return fugacityCoeff_[phaseIdx][compIdx]; }
+    Scalar fugacityCoefficient(int phaseIdx, int compIdx) const
+    { return fugacityCoefficient_[phaseIdx][compIdx]; }
 
     /*!
      * \brief The fugacity of a component in a phase [Pa]
      */
     Scalar fugacity(int phaseIdx, int compIdx) const
-    { return pressure_[phaseIdx]*fugacityCoeff_[phaseIdx][compIdx]*moleFrac_[phaseIdx][compIdx]; }
+    { return pressure_[phaseIdx]*fugacityCoefficient_[phaseIdx][compIdx]*moleFraction_[phaseIdx][compIdx]; }
 
     /*!
      * \brief The molar volume of a fluid phase [m^3/mol]
      */
     Scalar molarVolume(int phaseIdx) const
-    { return molarVolume_[phaseIdx]; }
+    { return 1/molarDensity(phaseIdx); }
 
     /*!
      * \brief The mass density of a fluid phase [kg/m^3]
      */
     Scalar density(int phaseIdx) const
-    {
-        return 
-            meanMolarMass_[phaseIdx]/molarVolume_[phaseIdx];
-    }
+    { return density_[phaseIdx]; }
 
     /*!
      * \brief The molar density of a fluid phase [mol/m^3]
      */
     Scalar molarDensity(int phaseIdx) const
-    { return 1/molarVolume(phaseIdx); }
+    { return density_[phaseIdx]/averageMolarMass(phaseIdx); }
+
 
     /*!
      * \brief The temperature of a fluid phase [K]
@@ -157,16 +155,16 @@ public:
     { return pressure_[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 enthalpy_[phaseIdx]; };
+    Scalar internalEnergy(int phaseIdx) const
+    { return internalEnergy_[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 enthalpy(phaseIdx) - pressure(phaseIdx)/(density(phaseIdx)); };
+    Scalar enthalpy(int phaseIdx) const
+    { return internalEnergy(phaseIdx) + pressure(phaseIdx)/density(phaseIdx); };
 
     /*!
      * \brief The dynamic viscosity of a fluid phase [Pa s]
@@ -200,28 +198,26 @@ public:
     /*!
      * \brief Set the mole fraction of a component in a phase []
      */
-    void setMoleFrac(int phaseIdx, int compIdx, Scalar value)
-    { moleFrac_[phaseIdx][compIdx] = value; }   
+    void setMoleFraction(int phaseIdx, int compIdx, Scalar value)
+    { moleFraction_[phaseIdx][compIdx] = value; }   
 
     /*!
      * \brief Set the fugacity of a component in a phase []
      */
-    void setFugacityCoeff(int phaseIdx, int compIdx, Scalar value)
-    { fugacityCoeff_[phaseIdx][compIdx] = value; }   
+    void setFugacityCoefficient(int phaseIdx, int compIdx, Scalar value)
+    { fugacityCoefficient_[phaseIdx][compIdx] = value; }   
 
     /*!
-     * \brief Set the molar volume of a phase [m^3/mol]
+     * \brief Set the density of a phase [kg / m^3]
      */
-    void setMolarVolume(int phaseIdx, Scalar value)
-    {
-        molarVolume_[phaseIdx] = value;
-    }
+    void setDensity(int phaseIdx, Scalar value)
+    { density_[phaseIdx] = value; }
 
     /*!
-     * \brief Set the specific enthalpy of a phase [J/m^3]
+     * \brief Set the specific internal energy of a phase [J/m^3]
      */
-    void setEnthalpy(int phaseIdx, Scalar value)
-    { enthalpy_[phaseIdx] = value; }
+    void setInternalEnergy(int phaseIdx, Scalar value)
+    { internalEnergy_[phaseIdx] = value; }
 
     /*!
      * \brief Set the dynamic viscosity of a phase [Pa s]
@@ -233,18 +229,18 @@ public:
      * \brief Calculatate the mean molar mass of a phase given that
      *        all mole fractions have been set
      */
-    void updateMeanMolarMass(int phaseIdx)
+    void updateAverageMolarMass(int phaseIdx)
     {
-        meanMolarMass_[phaseIdx] = 0;
-        sumMoleFrac_[phaseIdx] = 0;
+        averageMolarMass_[phaseIdx] = 0;
+        sumMoleFractions_[phaseIdx] = 0;
      
         for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
-            sumMoleFrac_[phaseIdx] += moleFrac_[phaseIdx][compIdx];
-            meanMolarMass_[phaseIdx] += moleFrac_[phaseIdx][compIdx]*StaticParameters::molarMass(compIdx);
+            sumMoleFractions_[phaseIdx] += moleFraction_[phaseIdx][compIdx];
+            averageMolarMass_[phaseIdx] += moleFraction_[phaseIdx][compIdx]*FluidSystem::molarMass(compIdx);
         }
-        sumMoleFrac_[phaseIdx] = std::max(1e-15, sumMoleFrac_[phaseIdx]);
-        Valgrind::CheckDefined(meanMolarMass_[phaseIdx]);
-        Valgrind::CheckDefined(sumMoleFrac_[phaseIdx]);
+        sumMoleFractions_[phaseIdx] = std::max(1e-15, sumMoleFractions_[phaseIdx]);
+        Valgrind::CheckDefined(averageMolarMass_[phaseIdx]);
+        Valgrind::CheckDefined(sumMoleFractions_[phaseIdx]);
     }
     
     /*!
@@ -256,15 +252,15 @@ public:
     {
         for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
             for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
-                moleFrac_[phaseIdx][compIdx] = fs.moleFrac(phaseIdx, compIdx);
-                fugacityCoeff_[phaseIdx][compIdx] = fs.fugacityCoeff(phaseIdx, compIdx);
+                moleFraction_[phaseIdx][compIdx] = fs.moleFraction(phaseIdx, compIdx);
+                fugacityCoefficient_[phaseIdx][compIdx] = fs.fugacityCoefficient(phaseIdx, compIdx);
             }
-            updateMeanMolarMass(phaseIdx);
+            updateAverageMolarMass(phaseIdx);
 
             pressure_[phaseIdx] = fs.pressure(phaseIdx);
             saturation_[phaseIdx] = fs.saturation(phaseIdx);
-            molarVolume_[phaseIdx] = fs.molarVolume(phaseIdx);
-            enthalpy_[phaseIdx] = fs.enthalpy(phaseIdx);
+            density_[phaseIdx] = fs.density(phaseIdx);
+            internalEnergy_[phaseIdx] = fs.internalEnergy(phaseIdx);
             viscosity_[phaseIdx] = fs.viscosity(phaseIdx);
             temperature_[phaseIdx] = fs.temperature(phaseIdx);
         }
@@ -283,30 +279,30 @@ public:
 #if HAVE_VALGRIND && ! defined NDEBUG
         for (int i = 0; i < numPhases; ++i) {
             for (int j = 0; j < numComponents; ++j) {
-                Valgrind::CheckDefined(fugacityCoeff_[i][j]);
-                Valgrind::CheckDefined(moleFrac_[i][j]);
+                Valgrind::CheckDefined(fugacityCoefficient_[i][j]);
+                Valgrind::CheckDefined(moleFraction_[i][j]);
             }
-            Valgrind::CheckDefined(meanMolarMass_[i]);
+            Valgrind::CheckDefined(averageMolarMass_[i]);
             Valgrind::CheckDefined(pressure_[i]);
             Valgrind::CheckDefined(saturation_[i]);
-            Valgrind::CheckDefined(molarVolume_[i]);
+            Valgrind::CheckDefined(density_[i]);
             Valgrind::CheckDefined(temperature_[i]);
-            //Valgrind::CheckDefined(enthalpy_[i]);
+            //Valgrind::CheckDefined(internalEnergy_[i]);
             Valgrind::CheckDefined(viscosity_[i]);
         }
 #endif // HAVE_VALGRIND
     }
 
 protected:
-    Scalar moleFrac_[numPhases][numComponents];
-    Scalar fugacityCoeff_[numPhases][numComponents];
+    Scalar moleFraction_[numPhases][numComponents];
+    Scalar fugacityCoefficient_[numPhases][numComponents];
 
-    Scalar meanMolarMass_[numPhases];
-    Scalar sumMoleFrac_[numPhases];
+    Scalar averageMolarMass_[numPhases];
+    Scalar sumMoleFractions_[numPhases];
     Scalar pressure_[numPhases];
     Scalar saturation_[numPhases];
-    Scalar molarVolume_[numPhases];
-    Scalar enthalpy_[numPhases];
+    Scalar density_[numPhases];
+    Scalar internalEnergy_[numPhases];
     Scalar viscosity_[numPhases];
     Scalar temperature_[numPhases];
 };
diff --git a/dumux/material/MpNcfluidstates/immisciblefluidstate.hh b/dumux/material/MpNcfluidstates/immisciblefluidstate.hh
index 586128df24..b3b548fcd3 100644
--- a/dumux/material/MpNcfluidstates/immisciblefluidstate.hh
+++ b/dumux/material/MpNcfluidstates/immisciblefluidstate.hh
@@ -65,13 +65,13 @@ public:
     /*!
      * \brief The mole fraction of a component in a phase []
      */
-    Scalar moleFrac(int phaseIdx, int compIdx) const
+    Scalar moleFraction(int phaseIdx, int compIdx) const
     { return (phaseIdx == compIdx)?1.0:0.0; }
 
     /*!
      * \brief The mass fraction of a component in a phase []
      */
-    Scalar massFrac(int phaseIdx, int compIdx) const
+    Scalar massFraction(int phaseIdx, int compIdx) const
     { return (phaseIdx == compIdx)?1.0:0.0; }
 
     /*!
@@ -79,7 +79,7 @@ public:
      *
      * We define this to be the same as the sum of all mass fractions.
      */
-    Scalar sumMoleFrac(int phaseIdx) const
+    Scalar sumMoleFractions(int phaseIdx) const
     { return 1.0; }
 
     /*!
@@ -87,13 +87,13 @@ public:
      *
      * We define this to be the same as the sum of all mole fractions.
      */
-    Scalar sumMassFrac(int phaseIdx) const
+    Scalar sumMassFractions(int phaseIdx) const
     { return 1.0; }
 
     /*!
-     * \brief The mean molar mass of a fluid phase [kg/mol]
+     * \brief The average molar mass of a fluid phase [kg/mol]
      */
-    Scalar meanMolarMass(int phaseIdx) const
+    Scalar averageMolarMass(int phaseIdx) const
     { return StaticParameters::molarMass(/*compIdx=*/phaseIdx); }
 
     /*!
@@ -106,7 +106,7 @@ public:
      * http://en.wikipedia.org/wiki/Concentration
      */
     Scalar molarity(int phaseIdx, int compIdx) const
-    { return molarDensity(phaseIdx)*moleFrac(phaseIdx, compIdx); }
+    { return molarDensity(phaseIdx)*moleFraction(phaseIdx, compIdx); }
 
     /*!
      * \brief The fugacity of a component in a phase [Pa]
@@ -120,9 +120,9 @@ public:
     Scalar fugacity(int phaseIdx, int compIdx) const
     { 
         if (phaseIdx == compIdx)
-            return fugacityCoeff(phaseIdx, compIdx)*moleFrac(phaseIdx, compIdx)*pressure(phaseIdx);
+            return pressure(phaseIdx);
         else
-            return 0.0;
+            return 0;
     };
     
     /*!
@@ -133,10 +133,10 @@ public:
      * infinite. Beware that this will very likely break your code if
      * you don't keep that in mind.
      */
-    Scalar fugacityCoeff(int phaseIdx, int compIdx) const
+    Scalar fugacityCoefficient(int phaseIdx, int compIdx) const
     {
         if (phaseIdx == compIdx)
-            return fugacityCoeff_[phaseIdx];
+            return 1.0;
         else
             return std::numeric_limits<Scalar>::infinity();
     }
@@ -145,19 +145,19 @@ public:
      * \brief The molar volume of a fluid phase [m^3/mol]
      */
     Scalar molarVolume(int phaseIdx) const
-    { return molarVolume_[phaseIdx]; }
+    { return 1/molarDensity(phaseIdx); }
 
     /*!
      * \brief The mass density of a fluid phase [kg/m^3]
      */
     Scalar density(int phaseIdx) const
-    { return molarDensity(phaseIdx)*meanMolarMass(phaseIdx); }
+    { return density_[phaseIdx]; }
 
     /*!
      * \brief The molar density of a fluid phase [mol/m^3]
      */
     Scalar molarDensity(int phaseIdx) const
-    { return 1/molarVolume(phaseIdx); }
+    { return density_[phaseIdx]/averageMolarMass(phaseIdx); }
 
     /*!
      * \brief The temperature of a fluid phase [K]
@@ -175,13 +175,13 @@ public:
      * \brief The specific enthalpy of a fluid phase [J/kg]
      */
     Scalar enthalpy(int phaseIdx) const
-    { return enthalpy_[phaseIdx]; }
+    { return internalEnergy_[phaseIdx] + pressure(phaseIdx)/(density(phaseIdx)); }
 
     /*!
      * \brief The specific internal energy of a fluid phase [J/kg]
      */
     Scalar internalEnergy(int phaseIdx) const
-    { return enthalpy(phaseIdx) - pressure(phaseIdx)/density(phaseIdx); }
+    { return internalEnergy_[phaseIdx]; }
 
     /*!
      * \brief The dynamic viscosity of a fluid phase [Pa s]
@@ -228,11 +228,10 @@ public:
     void assign(const FluidState &fs)
     {
         for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
-            fugacityCoeff_[phaseIdx] = fs.fugacityCoeff(phaseIdx, phaseIdx);
             pressure_[phaseIdx] = fs.pressure(phaseIdx);
             saturation_[phaseIdx] = fs.saturation(phaseIdx);
-            molarVolume_[phaseIdx] = fs.molarVolume(phaseIdx);
-            enthalpy_[phaseIdx] = fs.enthalpy(phaseIdx);
+            density_[phaseIdx] = fs.density(phaseIdx);
+            internalEnergy_[phaseIdx] = fs.internalEnergy(phaseIdx);
             viscosity_[phaseIdx] = fs.viscosity(phaseIdx);
         }
         temperature_ = fs.temperature(0);
@@ -257,25 +256,16 @@ public:
     { saturation_[phaseIdx] = value; }   
 
     /*!
-     * \brief Set the fugacity coefficient of a component in a phase []
-     *
-     * The phase is always index the same as the component, since we
-     * assume immiscibility.
+     * \brief Set the density of a phase [kg / m^3]
      */
-    void setFugacityCoeff(int compIdx, Scalar value)
-    { fugacityCoeff_[compIdx] = value; }   
+    void setDensity(int phaseIdx, Scalar value)
+    { density_[phaseIdx] = value; }   
 
     /*!
-     * \brief Set the molar volume of a phase [m^3/mol]
+     * \brief Set the specific internal energy of a phase [J/m^3]
      */
-    void setMolarVolume(int phaseIdx, Scalar value)
-    { molarVolume_[phaseIdx] = value; }   
-
-    /*!
-     * \brief Set the specific enthalpy of a phase [J/m^3]
-     */
-    void setEnthalpy(int phaseIdx, Scalar value)
-    { enthalpy_[phaseIdx] = value; }   
+    void setInternalEnergy(int phaseIdx, Scalar value)
+    { internalEnergy_[phaseIdx] = value; }   
 
     /*!
      * \brief Set the dynamic viscosity of a phase [Pa s]
@@ -296,12 +286,12 @@ public:
 #if HAVE_VALGRIND && ! defined NDEBUG
         for (int i = 0; i < numPhases; ++i) {
             //for (int j = 0; j < numComponents; ++j) {
-            //    Valgrind::CheckDefined(fugacityCoeff_[i][j]);
+            //    Valgrind::CheckDefined(fugacityCoefficient_[i][j]);
             //}
             Valgrind::CheckDefined(pressure_[i]);
             Valgrind::CheckDefined(saturation_[i]);
-            Valgrind::CheckDefined(molarVolume_[i]);
-            //Valgrind::CheckDefined(enthalpy_[i]);
+            Valgrind::CheckDefined(density_[i]);
+            //Valgrind::CheckDefined(internalEnergy_[i]);
             Valgrind::CheckDefined(viscosity_[i]);
         }
 
@@ -310,12 +300,10 @@ public:
     }
 
 protected:
-    Scalar fugacityCoeff_[numComponents];
-
     Scalar pressure_[numPhases];
     Scalar saturation_[numPhases];
-    Scalar molarVolume_[numPhases];
-    Scalar enthalpy_[numPhases];
+    Scalar density_[numPhases];
+    Scalar internalEnergy_[numPhases];
     Scalar viscosity_[numPhases];
     Scalar temperature_;
 };
diff --git a/dumux/material/MpNcfluidstates/trackinggenericfluidstate.hh b/dumux/material/MpNcfluidstates/trackinggenericfluidstate.hh
deleted file mode 100644
index 8960a49dc4..0000000000
--- a/dumux/material/MpNcfluidstates/trackinggenericfluidstate.hh
+++ /dev/null
@@ -1,398 +0,0 @@
-/*****************************************************************************
- *   Copyright (C) 2011 by Andreas Lauser                                    *
- *   Institute of Hydraulic Engineering                                      *
- *   University of Stuttgart, Germany                                        *
- *   email: <givenname>.<name>@iws.uni-stuttgart.de                          *
- *                                                                           *
- *   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 A GenericFluidState which keeps track on which variables changed.
- */
-#ifndef DUMUX_TRACKING_GENERIC_FLUID_STATE_HH
-#define DUMUX_TRACKING_GENERIC_FLUID_STATE_HH
-
-#include "genericfluidstate.hh"
-
-#include <dumux/common/valgrind.hh>
-
-#include <cmath>
-#include <algorithm>
-
-namespace Dumux
-{
-/*!
- * \brief A GenericFluidState which keeps track on which variables changed.
- */
-template <class Scalar, class StaticParameters>
-class TrackingGenericFluidState
-{
-    typedef Dumux::GenericFluidState<Scalar, StaticParameters> GenericFluidState;
-public:
-    enum { numComponents = StaticParameters::numComponents };
-    enum { numPhases = StaticParameters::numPhases };
-
-    TrackingGenericFluidState()
-    {
-        for (int i = 0; i < numPhases; ++i)
-            changed_[i].raw = ~((unsigned) 0x00);
-    }
-
-    /*****************************************************
-     * Generic access to fluid properties (No assumptions 
-     * on thermodynamic equilibrium required)
-     *****************************************************/
-    /*!
-     * \brief Returns the saturation of a phase []
-     */
-    Scalar saturation(int phaseIdx) const
-    { return fs_.saturation(phaseIdx); }
-
-    /*!
-     * \brief The mole fraction of a component in a phase []
-     */
-    Scalar moleFrac(int phaseIdx, int compIdx) const
-    { return fs_.moleFrac(phaseIdx, compIdx); }
-
-    /*!
-     * \brief The mass fraction of a component in a phase []
-     */
-    Scalar massFrac(int phaseIdx, int compIdx) const
-    { return fs_.massFrac(phaseIdx, compIdx); }
-
-    /*!
-     * \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 sumMoleFrac(int phaseIdx) const
-    { return fs_.sumMoleFrac(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 sumMassFrac(int phaseIdx) const
-    { return fs_.sumMassFrac(phaseIdx); }
-
-    /*!
-     * \brief The mean molar mass of a fluid phase [kg/mol]
-     */
-    Scalar meanMolarMass(int phaseIdx) const
-    { return fs_.meanMolarMass(phaseIdx); }
-
-    /*!
-     * \brief The molar concentration of a component in a phase [mol/m^3]
-     *
-     * This is often just called "concentration", but there are many
-     * other (though less common) measures for concentration.
-     *
-     * http://en.wikipedia.org/wiki/Concentration
-     */
-    Scalar molarity(int phaseIdx, int compIdx) const
-    { return fs_.molarity(phaseIdx, compIdx); }
-
-    /*!
-     * \brief The fugacity of a component in a phase [Pa]
-     */
-    Scalar fugacity(int phaseIdx, int compIdx) const
-    { return fs_.fugacity(phaseIdx, compIdx); }
-
-    /*!
-     * \brief The fugacity coefficient of a component in a phase []
-     */
-    Scalar fugacityCoeff(int phaseIdx, int compIdx) const
-    { return fs_.fugacityCoeff(phaseIdx, compIdx); }
-
-    /*!
-     * \brief The molar volume of a fluid phase [m^3/mol]
-     */
-    Scalar molarVolume(int phaseIdx) const
-    { return fs_.molarVolume(phaseIdx); }
-
-    /*!
-     * \brief The mass density of a fluid phase [kg/m^3]
-     */
-    Scalar density(int phaseIdx) const
-    { return fs_.density(phaseIdx); }
-
-    /*!
-     * \brief The molar density of a fluid phase [mol/m^3]
-     */
-    Scalar molarDensity(int phaseIdx) const
-    { return fs_.molarDensity(phaseIdx); }
-
-    /*!
-     * \brief The temperature of a fluid phase [K]
-     */
-    Scalar temperature(int phaseIdx) const
-    { return fs_.temperature(phaseIdx); }
-    
-    /*!
-     * \brief The pressure of a fluid phase [Pa]
-     */
-    Scalar pressure(int phaseIdx) const
-    { return fs_.pressure(phaseIdx); }
-
-    /*!
-     * \brief The specific enthalpy of a fluid phase [J/kg]
-     */
-    Scalar enthalpy(int phaseIdx) const
-    { return fs_.enthalpy(phaseIdx); }
-
-    /*!
-     * \brief The specific internal energy of a fluid phase [J/kg]
-     */
-    Scalar internalEnergy(int phaseIdx) const
-    { return fs_.internalEnergy(phaseIdx); }
-
-    /*!
-     * \brief The dynamic viscosity of a fluid phase [Pa s]
-     */
-    Scalar viscosity(int phaseIdx) const
-    { return fs_.viscosity(phaseIdx); };
-
-    /*****************************************************
-     * Setter methods. Note that these are not part of the 
-     * generic FluidState interface but specific for each
-     * implementation...
-     *****************************************************/
-    /*!
-     * \brief Set the temperature [K] of a fluid phase
-     */
-    void setTemperature(int phaseIdx, Scalar value)
-    { 
-        changed_[phaseIdx].fields.temperature = 1;
-        fs_.setTemperature(phaseIdx, value);
-    };
-
-    /*!
-     * \brief Set the fluid pressure of a phase [Pa]
-     */
-    void setPressure(int phaseIdx, Scalar value)
-    { 
-        changed_[phaseIdx].fields.pressure = 1;
-        fs_.setPressure(phaseIdx, value);
-    };
-
-    /*!
-     * \brief Set the saturation of a phase []
-     */
-    void setSaturation(int phaseIdx, Scalar value)
-    { 
-        changed_[phaseIdx].fields.saturation = 1;
-        fs_.setSaturation(phaseIdx, value);
-    };
-
-    /*!
-     * \brief Set the mole fraction of a component in a phase []
-     */
-    void setMoleFrac(int phaseIdx, int compIdx, Scalar value)
-    { 
-        changed_[phaseIdx].fields.moleFrac |= 1 << compIdx;
-        fs_.setMoleFrac(phaseIdx, compIdx, value);
-    };
-
-
-    /*!
-     * \brief Set the fugacity of a component in a phase []
-     */
-    void setFugacityCoeff(int phaseIdx, int compIdx, Scalar value)
-    { 
-        changed_[phaseIdx].fields.fugacityCoeff |= 1 << compIdx;
-        fs_.setFugacityCoeff(phaseIdx, compIdx, value);
-    }
-
-    /*!
-     * \brief Set the molar volume of a phase [m^3/mol]
-     */
-    void setMolarVolume(int phaseIdx, Scalar value)
-    { 
-        changed_[phaseIdx].fields.molarVolume = 1;
-        fs_.setMolarVolume(phaseIdx, value);
-    }
-
-    /*!
-     * \brief Set the specific enthalpy of a phase [J/m^3]
-     */
-    void setEnthalpy(int phaseIdx, Scalar value)
-    { 
-        changed_[phaseIdx].fields.enthalpy = 1;
-        fs_.setEnthalpy(phaseIdx, value);
-    }   
-
-    /*!
-     * \brief Set the dynamic viscosity of a phase [Pa s]
-     */
-    void setViscosity(int phaseIdx, Scalar value)
-    { 
-        changed_[phaseIdx].fields.viscosity = 1;
-        fs_.setViscosity(phaseIdx, value);
-    }   
-
-    /*!
-     * \brief Calculatate the mean molar mass of a phase given that
-     *        all mole fractions have been set
-     */
-    void updateMeanMolarMass(int phaseIdx)
-    { fs_.updateMeanMolarMass(phaseIdx); }
-    
-
-    /*!
-     * \brief Reset the dirty tracking for a phase.
-     */
-    void setPhaseClean(int phaseIdx)
-    {
-        // set all tracking bits for the phase to zero
-        changed_[phaseIdx].raw = 0;
-    }
-
-    /*!
-     * \brief Reset the dirty tracking of all phases.
-     */
-    void setClean()
-    {
-        for (int i = 0; i < numPhases; ++i)
-            changed_[i].raw = 0;
-    }
-
-    /*!
-     * \brief Return if all phases are clean.
-     */
-    bool isClean() const
-    { 
-        for (int i = 0; i < numPhases; ++i)
-            if (changed_[i].raw != 0)
-                return false;
-        return true;
-    };
-
-    /*!
-     * \brief Find out whether the values for a phase been since the
-     *        last call to cleanPhase().
-     */
-    bool phaseClean(int phaseIdx) const
-    { return changed_[phaseIdx].raw == 0; }
-
-
-    /*!
-     * \brief Find out whether a specific phase's temperature has been
-     *        changed.
-     */
-    bool temperatureClean(int phaseIdx) const
-    { return changed_[phaseIdx].fields.temperature == 0; }
-
-    /*!
-     * \brief Find out whether a specific phase's saturation has been
-     *        changed.
-     */
-    bool saturationClean(int phaseIdx) const
-    { return changed_[phaseIdx].fields.saturation == 0; }
-
-    /*!
-     * \brief Find out whether a specific phase's pressure has been
-     *        changed.
-     */
-    bool pressureClean(int phaseIdx) const
-    { return changed_[phaseIdx].fields.pressure == 0; }
-
-    /*!
-     * \brief Find out whether a specific phase's enthalpy has been
-     *        changed.
-     */
-    bool enthalpyClean(int phaseIdx) const
-    { return changed_[phaseIdx].enthalpy == 0; }
-
-    /*!
-     * \brief Find out whether a specific phase's temperature has been
-     *        changed.
-     */
-    bool molarVolumeClean(int phaseIdx) const
-    { return changed_[phaseIdx].fields.molarVolume == 0; }
-
-    /*!
-     * \brief Find out whether some of the mole fractions of a phase
-     *        have been changed since the last call to
-     *        setCleanPhase().
-     */
-    bool moleFracsClean(int phaseIdx) const
-    { return changed_[phaseIdx].fields.moleFrac == 0; }
-
-    /*!
-     * \brief Find out whether a specific mole fraction of a component
-     *        in a phase has been changed since the last call to
-     *        setCleanPhase().
-     */
-    bool moleFracsClean(int phaseIdx, int compIdx) const
-    { return (changed_[phaseIdx].fields.moleFrac  & (1 << compIdx)) == 0; }
-
-    /*!
-     * \brief Find out whether at least one fugacity of a component in
-     *        a phase has been changed since the last call to
-     *        setCleanPhase().
-     */
-    bool fugacityCoeffsClean(int phaseIdx) const
-    { return changed_[phaseIdx].fields.fugacitCoeffs == 0; }
-
-    /*!
-     * \brief Find out whether a specific fugacity of a component in a
-     *        phase has been changed since the last call to
-     *        setCleanPhase().
-     */
-    bool fugacityCoeffClean(int phaseIdx, int compIdx) const
-    { return (changed_[phaseIdx].fields.fugacityCoeffs  & (1 << compIdx)) == 0; }
-    
-
-    /*!
-     * \brief Make sure that all attributes are defined.
-     *
-     * This method does not do anything if the program is not run
-     * under valgrind. If it is, then valgrind will print an error
-     * message if some attributes of the object have not been properly
-     * defined.
-     */
-    void checkDefined() const
-    {
-#if HAVE_VALGRIND && ! defined NDEBUG
-        Valgrind::CheckDefined(changed_);
-        fs_.checkDefined();
-#endif // HAVE_VALGRIND
-    }
-
-protected:
-    static_assert(sizeof(unsigned)*8 >= 5 + 2*numComponents,
-                  "Too many components to use the bitfield trick");
-    union {
-        unsigned raw;
-        struct {
-            unsigned char moleFrac : numComponents;
-            unsigned char fugacityCoeff : numComponents;
-
-            unsigned char molarVolume : 1;
-            unsigned char temperature : 1;
-            unsigned char saturation : 1;
-            unsigned char pressure : 1;
-            unsigned char enthalpy : 1;
-            unsigned char viscosity : 1;
-        } __attribute__((__packed__)) fields;
-    } changed_[numPhases];
-    GenericFluidState fs_;
-};
-
-} // end namepace Dumux
-
-#endif
diff --git a/dumux/material/MpNcfluidsystems/h2on2fluidsystem.hh b/dumux/material/MpNcfluidsystems/h2on2fluidsystem.hh
index c03aca5fe4..283e9385a3 100644
--- a/dumux/material/MpNcfluidsystems/h2on2fluidsystem.hh
+++ b/dumux/material/MpNcfluidsystems/h2on2fluidsystem.hh
@@ -1,5 +1,5 @@
 /*****************************************************************************
- *   Copyright (C) 2009-2010 by Andreas Lauser                               *
+ *   Copyright (C) 2009-2011 by Andreas Lauser                               *
  *   Institute of Hydraulic Engineering                                      *
  *   University of Stuttgart, Germany                                        *
  *   email: <givenname>.<name>@iws.uni-stuttgart.de                          *
@@ -38,14 +38,32 @@
 
 #include <dumux/common/exceptions.hh>
 
+#include "nullparametercache.hh"
+
 namespace Dumux
 {
+/*!
+ * \brief A twophase fluid system with water and nitrogen as components.
+ */
+template <class Scalar>
+class H2ON2FluidSystem
+{
+    // convenience typedefs
+    typedef Dumux::IdealGas<Scalar> IdealGas;
 
+    typedef Dumux::H2O<Scalar> IapwsH2O;
+    typedef Dumux::SimpleH2O<Scalar> SimpleH2O;
+
+    typedef Dumux::TabulatedComponent<Scalar, IapwsH2O > TabulatedH2O;
+
+    typedef Dumux::N2<Scalar> SimpleN2;
+
+public:
+    //! The type of parameter cache objects
+    typedef Dumux::NullParameterCache ParameterCache;
 
-template <class Scalar>
-struct H2ON2StaticParameters {
     /****************************************
-     * Fluid phase parameters
+     * Fluid phase related static parameters
      ****************************************/
 
     //! Number of phases in the fluid system
@@ -55,38 +73,15 @@ struct H2ON2StaticParameters {
     static constexpr int lPhaseIdx = 0;
     //! Index of the gas phase
     static constexpr int gPhaseIdx = 1;
-    //! Index of the solid phase
-    static constexpr int sPhaseIdx = 2;
     
-    //! The component for pure water 
-    typedef Dumux::SimpleH2O<Scalar> SimpleH2O;
-    typedef Dumux::H2O<Scalar> IapwsH2O;
-    typedef Dumux::TabulatedComponent<Scalar, IapwsH2O> TabulatedH2O;
-
-//    typedef IapwsH2O H2O;
+    //! The components for pure water
     typedef TabulatedH2O H2O;
-//    typedef SimpleH2O H2O;
-    
-    //! The component for pure nitrogen
-    typedef Dumux::N2<Scalar> N2;
-
-    /*!
-     * \brief Initialize the static parameters.
-     */
-    static void init(Scalar tempMin, Scalar tempMax, unsigned nTemp,
-                     Scalar pressMin, Scalar pressMax, unsigned nPress)
-    {
-        if (H2O::isTabulated) {
-            std::cout << "Initializing tables for the H2O fluid properties ("
-                      << nTemp*nPress
-                      << " entries).\n";
-            
-            TabulatedH2O::init(tempMin, tempMax, nTemp,
-                               pressMin, pressMax, nPress);
-        }
+    //typedef SimpleH2O H2O;
+    //typedef IapwsH2O H2O;
 
-    }
-   
+    //! The components for pure nitrogen
+    typedef SimpleN2 N2;
+ 
     /*!
      * \brief Return the human readable name of a fluid phase
      */
@@ -94,11 +89,10 @@ struct H2ON2StaticParameters {
     {
         static const char *name[] = {
             "l",
-            "g",
-            "s"
+            "g"          
         };
 
-        assert(0 <= phaseIdx && phaseIdx < numPhases + 1);
+        assert(0 <= phaseIdx && phaseIdx < numPhases);
         return name[phaseIdx];
     }
 
@@ -133,7 +127,7 @@ struct H2ON2StaticParameters {
     }
 
     /****************************************
-     * Component related parameters
+     * Component related static parameters
      ****************************************/
 
     //! Number of components in the fluid system
@@ -220,81 +214,68 @@ struct H2ON2StaticParameters {
         assert(0 <= compIdx && compIdx < numComponents);
         return accFac[compIdx];
     };
-};
-
-/*!
- * \brief A twophase fluid system with water and nitrogen as components.
- */
-template <class Scalar>
-class H2ON2FluidSystem : public H2ON2StaticParameters<Scalar>
-{
-public:
-    typedef Dumux::H2ON2StaticParameters<Scalar> StaticParameters;
-    typedef Dumux::GenericFluidState<Scalar, StaticParameters> MutableParameters;
-
-private:
-    // convenience typedefs
-    typedef StaticParameters SP;
-    typedef Dumux::IdealGas<Scalar> IdealGas;
 
-public:
+    /*!
+     * \brief Initialize the fluid system's static parameters generically
+     *
+     * If a tabulated H2O component is used, we do our best to create
+     * tables that always work.
+     */
+    static void init()
+    {
+        init(/*tempMin=*/273.15, 
+             /*tempMax=*/623.15, 
+             /*numTemp=*/100,
+             /*pMin=*/-10,
+             /*pMax=*/20e6,
+             /*numP=*/200);
+    }
 
     /*!
-     * \brief Initialize the fluid system's static parameters
+     * \brief Initialize the fluid system's static parameters using 
+     *        problem specific temperature and pressure ranges
      */
     static void init(Scalar tempMin, Scalar tempMax, unsigned nTemp,
-    				 Scalar pressMin, Scalar pressMax, unsigned nPress)
+                     Scalar pressMin, Scalar pressMax, unsigned nPress)
     {
-        SP::init(tempMin, tempMax, nTemp,
-                 pressMin, pressMax, nPress);
+        if (H2O::isTabulated) {
+            std::cout << "Initializing tables for the H2O fluid properties ("
+                      << nTemp*nPress
+                      << " entries).\n";
+            
+            TabulatedH2O::init(tempMin, tempMax, nTemp,
+                               pressMin, pressMax, nPress);
+        }
     }
 
     /*!
      * \brief Calculate the molar volume [m^3/mol] of a fluid phase
      */
-    static Scalar computeMolarVolume(MutableParameters &params,
-                                     int phaseIdx)
+    template <class FluidState>
+    static Scalar density(const FluidState &fluidState,
+                          const ParameterCache &paramCache,
+                          int phaseIdx)
     {
-        assert(0 <= phaseIdx  && phaseIdx <= SP::numPhases);
+        assert(0 <= phaseIdx  && phaseIdx <= numPhases);
 
-        params.updateMeanMolarMass(phaseIdx);
-        Scalar T = params.temperature(phaseIdx);
-        Scalar p = params.pressure(phaseIdx);
+        Scalar T = fluidState.temperature(phaseIdx);
+        Scalar p = fluidState.pressure(phaseIdx);
 
         switch (phaseIdx) {
-        case SP::lPhaseIdx:
+        case lPhaseIdx:
             // assume pure water where one water molecule gets
             // replaced by one nitrogen molecule
-            return SP::H2O::molarMass()/SP::H2O::liquidDensity(T, p);
-        case SP::gPhaseIdx:
+            return H2O::liquidDensity(T, p);
+        case gPhaseIdx:
             // assume ideal gas
-            return 1.0 / IdealGas::concentration(T, p);
+            return
+                IdealGas::molarDensity(T, p) 
+                * fluidState.averageMolarMass(gPhaseIdx);
         }
         
         DUNE_THROW(Dune::InvalidStateException, "Unhandled phase index " << phaseIdx);
     };
 
-    /*!
-     * \brief Calculate the fugacity of a component in a fluid phase
-     *        [Pa]
-     *
-     * The components chemical \f$mu_\kappa\f$ potential is connected
-     * to the component's fugacity \f$f_\kappa\f$ by the relation
-     *
-     * \f[ \mu_\kappa = R T_\alpha \mathrm{ln} \frac{f_\kappa}{p_\alpha} \f]
-     *
-     * where \f$p_\alpha\f$ and \f$T_\alpha\f$ are the fluid phase'
-     * pressure and temperature.
-     */
-    static Scalar computeFugacity(const MutableParameters &params,
-                                  int phaseIdx,
-                                  int compIdx)
-    {
-        Scalar x = params.moleFrac(phaseIdx,compIdx);
-        Scalar p = params.pressure(phaseIdx);
-        return x*p*computeFugacityCoeff(params, phaseIdx, compIdx);
-    };
-
     /*!
      * \brief Calculate the fugacity coefficient [Pa] of an individual
      *        component in a fluid phase
@@ -305,23 +286,24 @@ public:
      *
      * \f[ f_\kappa = \phi_\kappa * x_{\kappa} \f]
      */
-    static Scalar computeFugacityCoeff(MutableParameters &params,
-                                       int phaseIdx,
-                                       int compIdx)
+    template <class FluidState>
+    static Scalar fugacityCoefficient(const FluidState &fluidState,
+                                      const ParameterCache &paramCache,
+                                      int phaseIdx,
+                                      int compIdx)
     {
-        assert(0 <= phaseIdx  && phaseIdx <= SP::numPhases);
-        assert(0 <= compIdx  && compIdx <= SP::numComponents);
+        assert(0 <= phaseIdx  && phaseIdx <= numPhases);
+        assert(0 <= compIdx  && compIdx <= numComponents);
 
-        params.updateMeanMolarMass(phaseIdx);
-        Scalar T = params.temperature(phaseIdx);
-        Scalar p = params.pressure(phaseIdx);
+        Scalar T = fluidState.temperature(phaseIdx);
+        Scalar p = fluidState.pressure(phaseIdx);
         switch (phaseIdx) {
-        case SP::lPhaseIdx: 
+        case lPhaseIdx: 
             switch (compIdx) {
-            case SP::H2OIdx: return SP::H2O::vaporPressure(T)/p;
-            case SP::N2Idx: return BinaryCoeff::H2O_N2::henry(T)/p;
+            case H2OIdx: return H2O::vaporPressure(T)/p;
+            case N2Idx: return BinaryCoeff::H2O_N2::henry(T)/p;
             };
-        case SP::gPhaseIdx:
+        case gPhaseIdx:
             return 1.0; // ideal gas
         };
 
@@ -331,21 +313,22 @@ public:
     /*!
      * \brief Calculate the dynamic viscosity of a fluid phase [Pa*s]
      */
-    static Scalar computeViscosity(MutableParameters &params,
-                                   int phaseIdx)
+    template <class FluidState>
+    static Scalar viscosity(const FluidState &fluidState,
+                            const ParameterCache &paramCache,
+                            int phaseIdx)
     {
-        assert(0 <= phaseIdx  && phaseIdx <= SP::numPhases);
+        assert(0 <= phaseIdx  && phaseIdx <= numPhases);
 
-        params.updateMeanMolarMass(phaseIdx);
-        Scalar T = params.temperature(phaseIdx);
-        Scalar p = params.pressure(phaseIdx);
+        Scalar T = fluidState.temperature(phaseIdx);
+        Scalar p = fluidState.pressure(phaseIdx);
         switch (phaseIdx) {
-        case SP::lPhaseIdx:
+        case lPhaseIdx:
             // assume pure water for the liquid phase
-            return SP::H2O::liquidViscosity(T, p);
-        case SP::gPhaseIdx:
+            return H2O::liquidViscosity(T, p);
+        case gPhaseIdx:
             // assume pure water for the gas phase
-            return SP::N2::gasViscosity(T, p);
+            return N2::gasViscosity(T, p);
         }
         
         DUNE_THROW(Dune::InvalidStateException, "Unhandled phase index " << phaseIdx);
@@ -370,9 +353,11 @@ public:
      * where \f$p_\alpha\f$ and \f$T_\alpha\f$ are the fluid phase'
      * pressure and temperature.
      */
-    static Scalar computeDiffusionCoeff(const MutableParameters &params,
-                                        int phaseIdx,
-                                        int compIdx)
+    template <class FluidState>
+    static Scalar diffusionCoefficient(const FluidState &fluidState,
+                                       const ParameterCache &paramCache,
+                                       int phaseIdx,
+                                       int compIdx)
     {
         // TODO!
         DUNE_THROW(Dune::NotImplemented, "Diffusion coefficients");
@@ -383,20 +368,21 @@ public:
      *        return the binary diffusion coefficient for components
      *        \f$i\f$ and \f$j\f$ in this phase.
      */
-    static Scalar computeBinaryDiffCoeff(MutableParameters &params, 
-                                         int phaseIdx,
-                                         int compIIdx,
-                                         int compJIdx)
+    template <class FluidState>
+    static Scalar binaryDiffusionCoefficient(const FluidState &fluidState, 
+                                             const ParameterCache &paramCache,
+                                             int phaseIdx,
+                                             int compIIdx,
+                                             int compJIdx)
                                   
     {
-        params.updateMeanMolarMass(phaseIdx);
         if (compIIdx > compJIdx)
             std::swap(compIIdx, compJIdx);
 
 #ifndef NDEBUG
         if (compIIdx == compJIdx ||
-            phaseIdx > SP::numPhases - 1 ||
-            compJIdx > SP::numComponents - 1)
+            phaseIdx > numPhases - 1 ||
+            compJIdx > numComponents - 1)
         {
             DUNE_THROW(Dune::InvalidStateException,
                        "Binary diffusion coefficient of components "
@@ -405,26 +391,26 @@ public:
         }
 #endif
 
-        Scalar T = params.temperature(phaseIdx);
-        Scalar p = params.pressure(phaseIdx);       
+        Scalar T = fluidState.temperature(phaseIdx);
+        Scalar p = fluidState.pressure(phaseIdx);       
 
         switch (phaseIdx) {
-        case SP::lPhaseIdx:
+        case lPhaseIdx:
             switch (compIIdx) {
-            case SP::H2OIdx:
+            case H2OIdx:
                 switch (compJIdx) {
-                case SP::N2Idx: return BinaryCoeff::H2O_N2::liquidDiffCoeff(T, p);
+                case N2Idx: return BinaryCoeff::H2O_N2::liquidDiffCoeff(T, p);
                 }
             default:
                 DUNE_THROW(Dune::InvalidStateException,
                            "Binary diffusion coefficients of trace "
                            "substances in liquid phase is undefined!\n");
             }
-        case SP::gPhaseIdx:
+        case gPhaseIdx:
             switch (compIIdx) {
-            case SP::H2OIdx:
+            case H2OIdx:
                 switch (compJIdx) {
-                case SP::N2Idx: return BinaryCoeff::H2O_N2::gasDiffCoeff(T, p);
+                case N2Idx: return BinaryCoeff::H2O_N2::gasDiffCoeff(T, p);
                 }
             }
         }
@@ -438,98 +424,87 @@ public:
     /*!
      * \brief Given a phase's composition, temperature, pressure and
      *        density, calculate its specific enthalpy [J/kg].
+     *
+     *  \todo This fluid system neglects the contribution of
+     *        gas-molecules in the liquid phase. This contribution is
+     *        probably not big. Somebody would have to find out the
+     *        enthalpy of solution for this system. ...
      */
-
-    /*!
-     *  \todo This system neglects the contribution of gas-molecules in the liquid phase.
-     *        This contribution is probably not big. Somebody would have to find out the enthalpy of solution for this system. ...
-     */
-    static Scalar computeEnthalpy(MutableParameters &params, 
-                                  int phaseIdx)
+    template <class FluidState>
+    static Scalar internalEnergy(const FluidState &fluidState, 
+                                 const ParameterCache &paramCache,
+                                 int phaseIdx)
     {
-        params.updateMeanMolarMass(phaseIdx);
-        Scalar T = params.temperature(phaseIdx);
-        Scalar p = params.pressure(phaseIdx);
+        Scalar T = fluidState.temperature(phaseIdx);
+        Scalar p = fluidState.pressure(phaseIdx);
         Valgrind::CheckDefined(T);
         Valgrind::CheckDefined(p);
-        if (phaseIdx == SP::lPhaseIdx) {
+        if (phaseIdx == lPhaseIdx) {
             // TODO: correct way to deal with the solutes???
             return 
-                SP::H2O::liquidEnthalpy(T, p) ;
+                H2O::liquidInternalEnergy(T, p) ;
         }
         else {
             // assume ideal gas
-            Scalar XH2O = params.massFrac(SP::gPhaseIdx, SP::H2OIdx);
-            Scalar XN2 = params.massFrac(SP::gPhaseIdx, SP::N2Idx);          
+            Scalar XH2O = fluidState.massFrac(gPhaseIdx, H2OIdx);
+            Scalar XN2 = fluidState.massFrac(gPhaseIdx, N2Idx);          
             Scalar result = 0;
-            result += XH2O*SP::H2O::gasEnthalpy(T, p);
-            result += XN2*SP::N2::gasEnthalpy(T, p);
+            result += XH2O*H2O::gasInternalEnergy(T, p);
+            result += XN2*N2::gasInternalEnergy(T, p);
             return result;
         }
     }
 
     /*!
-     * \brief Given a phase's composition, temperature, pressure and
-     *        density, calculate its specific internal energy [J/kg].
-     */
-    static Scalar computeInternalEnergy(MutableParameters &params, 
-                                        int phaseIdx)
-    {
-        params.updateMeanMolarMass(phaseIdx);
-        Scalar T = params.temperature(phaseIdx);
-        Scalar p = params.pressure(phaseIdx);
-        Scalar rho = params.density(phaseIdx);
-        return
-            computeEnthalpy(params, phaseIdx) -
-            p/rho;
-    }
-
-    /*!
-     * \brief Thermal conductivity of phases.
+     * \brief Thermal conductivity of a fluid phase.
      *
      * Use the conductivity of air and water as a first approximation.
      * Source:
      * http://en.wikipedia.org/wiki/List_of_thermal_conductivities
      */
-    static Scalar computeThermalConductivity(const MutableParameters &params,
-                                             int phaseIdx)
+    template <class FluidState>
+    static Scalar thermalConductivity(const FluidState &fluidState,
+                                      const ParameterCache &paramCache,
+                                      int phaseIdx)
     {
 //    	TODO thermal conductivity is a function of:
-//        Scalar p = params.pressure(phaseIdx);
-//        Scalar T = params.temperature(phaseIdx);
-//        Scalar x = params.moleFrac(phaseIdx,compIdx);
+//        Scalar p = fluidState.pressure(phaseIdx);
+//        Scalar T = fluidState.temperature(phaseIdx);
+//        Scalar x = fluidState.moleFrac(phaseIdx,compIdx);
 #warning: so far rough estimates from wikipedia
         switch (phaseIdx) {
-        case SP::lPhaseIdx: // use conductivity of pure water
+        case lPhaseIdx: // use conductivity of pure water
             return  0.6;   // conductivity of water[W / (m K ) ]
-        case SP::gPhaseIdx:// use conductivity of pure air
+        case gPhaseIdx:// use conductivity of pure air
             return 0.025; // conductivity of air [W / (m K ) ]
         }
         DUNE_THROW(Dune::InvalidStateException, "Unhandled phase index " << phaseIdx);
     }
     
     /*!
-     * \brief Specific isobaric heat capacity of liquid water / air
+     * \brief Specific isobaric heat capacity of a fluid phase.
      *        \f$\mathrm{[J/kg]}\f$.
      *
      * \param params    mutable parameters
      * \param phaseIdx  for which phase to give back the heat capacity
      */
-    static Scalar computeHeatCapacity(const MutableParameters &params,
-                                      int phaseIdx)
+    template <class FluidState>
+    static Scalar heatCapacity(const FluidState &fluidState,
+                               const ParameterCache &paramCache,
+                               int phaseIdx)
     {
 //        http://en.wikipedia.org/wiki/Heat_capacity
 #warning: so far rough estimates from wikipedia
 //      TODO heatCapacity is a function of composition.
-//        Scalar p = params.pressure(phaseIdx);
-//        Scalar T = params.temperature(phaseIdx);
-//        Scalar x = params.moleFrac(phaseIdx,compIdx);
+//        Scalar p = fluidState.pressure(phaseIdx);
+//        Scalar T = fluidState.temperature(phaseIdx);
+//        Scalar x = fluidState.moleFrac(phaseIdx,compIdx);
         switch (phaseIdx) {
-        case SP::lPhaseIdx: // use heat capacity of pure liquid water
+        case lPhaseIdx: // use heat capacity of pure liquid water
             return  4181.3;  // @(25°C) !!!
             /* [J/(kg K)]*/ /* not working because ddgamma_ddtau is not defined*/ /* Dumux::H2O<Scalar>::liquidHeatCap_p(T,
                                              p); */
-        case SP::gPhaseIdx:
+        case gPhaseIdx:
             return  1003.5 ; // @ (0°C) !!!
             /* [J/(kg K)]*/ /* not working because ddgamma_ddtau is not defined*/ /*Dumux::H2O<Scalar>::gasHeatCap_p(T,
                                           p) ;*/
diff --git a/dumux/material/MpNcfluidsystems/nullparametercache.hh b/dumux/material/MpNcfluidsystems/nullparametercache.hh
new file mode 100644
index 0000000000..2539519780
--- /dev/null
+++ b/dumux/material/MpNcfluidsystems/nullparametercache.hh
@@ -0,0 +1,56 @@
+/*****************************************************************************
+ *   Copyright (C) 2011 by Andreas Lauser                                    *
+ *   Institute of Hydraulic Engineering                                      *
+ *   University of Stuttgart, Germany                                        *
+ *   email: <givenname>.<name>@iws.uni-stuttgart.de                          *
+ *                                                                           *
+ *   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 The a parameter cache which does nothing
+ */
+#ifndef DUMUX_NULL_PARAMETER_CACHE_HH
+#define DUMUX_NULL_PARAMETER_CACHE_HH
+
+#include "parametercachebase.hh"
+
+namespace Dumux
+{
+/*!
+ * \brief The a parameter cache which does nothing
+ */
+class NullParameterCache : public ParameterCacheBase<NullParameterCache>
+{
+public:
+    NullParameterCache() 
+    {};
+    
+    template <class FluidState>
+    void updateAll(const FluidState &fs)
+    {
+    };
+    
+    /*!
+     * \brief Update all cached parameters of a specific fluid phase
+     */
+    template <class FluidState>
+    void updatePhase(const FluidState &fs, int phaseIdx)
+    {};
+};
+
+} // end namepace
+
+#endif
diff --git a/dumux/material/MpNcfluidsystems/parametercachebase.hh b/dumux/material/MpNcfluidsystems/parametercachebase.hh
new file mode 100644
index 0000000000..71b2984f60
--- /dev/null
+++ b/dumux/material/MpNcfluidsystems/parametercachebase.hh
@@ -0,0 +1,120 @@
+/*****************************************************************************
+ *   Copyright (C) 2011 by Andreas Lauser                                    *
+ *   Institute of Hydraulic Engineering                                      *
+ *   University of Stuttgart, Germany                                        *
+ *   email: <givenname>.<name>@iws.uni-stuttgart.de                          *
+ *                                                                           *
+ *   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 The base class of the parameter cache classes for fluid systems
+ */
+#ifndef DUMUX_PARAMETER_CACHE_BASE_HH
+#define DUMUX_PARAMETER_CACHE_BASE_HH
+
+namespace Dumux
+{
+/*!
+ * \brief The base class of the parameter cache classes for fluid systems
+ */
+template <class Implementation>
+class ParameterCacheBase
+{
+public:
+    ParameterCacheBase() 
+    {};
+    
+    template <class FluidState>
+    void updateAll(const FluidState &fs)
+    {
+        for (int phaseIdx = 0; phaseIdx < FluidState::numPhases; ++phaseIdx)
+            updatePhase(fs, phaseIdx);
+    };
+    
+    /*!
+     * \brief Update all cached parameters of a specific fluid phase
+     */
+    template <class FluidState>
+    void updatePhase(const FluidState &fs, int phaseIdx)
+    {};
+
+    /*!
+     * \brief Update all cached parameters of a specific fluid phase
+     *        which depend on temperature
+     *
+     * *Only* use this method if only the temperature of a phase
+     * changed between two update*() calls. If more changed, call
+     * updatePhase()!
+     */
+    template <class FluidState> void
+    updatePhaseTemperature(const FluidState &fs, int phaseIdx)
+    {
+        asImp_().updatePhase(fs, phaseIdx);
+    };
+
+    /*!
+     * \brief Update all cached parameters of a specific fluid phase
+     *        which depend on pressure
+     *
+     * *Only* use this method if only the pressure of a phase changed
+     * between two update*() calls. If more changed, call
+     * updatePhase()!
+     */
+    template <class FluidState>
+    void updatePhasePressure(const FluidState &fs, int phaseIdx)
+    {
+        asImp_().updatePhase(fs, phaseIdx);
+    };
+
+    /*!
+     * \brief Update all cached parameters of a specific fluid phase
+     *        which depend on composition
+     *
+     * *Only* use this method if neither the pressure nor the
+     * temperature of the phase changed between two update*()
+     * calls. If more changed, call updatePhase()!
+     */
+    template <class FluidState>
+    void updatePhaseComposition(const FluidState &fs, int phaseIdx)
+    {
+        asImp_().updatePhase(fs, phaseIdx);
+    };
+
+    /*!
+     * \brief Update all cached parameters of a specific fluid phase
+     *        which depend on the mole fraction of a single component
+     *
+     * *Only* use this method if just a single component's
+     * concentration changed between two update*() calls. If more than
+     * one concentration changed, call updatePhaseComposition() of
+     * updatePhase()!
+     */
+    template <class FluidState>
+    void updatePhaseSingleMoleFraction(const FluidState &fs, 
+                                       int phaseIdx,
+                                       int compIdx)
+    {
+        asImp_().updatePhaseComposition(fs, phaseIdx);
+    };
+
+private:
+    Implementation &asImp_()
+    { return *static_cast<Implementation*>(this); } 
+};
+
+} // end namepace
+
+#endif
diff --git a/dumux/material/idealgas.hh b/dumux/material/idealgas.hh
index 8aaba763a0..66a60313c5 100644
--- a/dumux/material/idealgas.hh
+++ b/dumux/material/idealgas.hh
@@ -51,19 +51,24 @@ public:
 
     /*!
      * \brief The pressure of the gas in \f$\mathrm{[N/m^2]}\f$, depending on
-     *        concentration and temperature.
+     *        the molar density and temperature.
      */
     static Scalar pressure(Scalar temperature,
-                           Scalar concentration)
-    { return R*temperature*concentration; }
+                           Scalar rhoMolar)
+    { return R*temperature*rhoMolar; }
 
     /*!
-     * \brief The molar concentration of the gas in \f$\mathrm{[mol/m^3]}\f$, depending on
-     *        pressure and temperature.
+     * \brief The molar density of the gas \f$\mathrm{[mol/m^3]}\f$,
+     *        depending on pressure and temperature.
      */
-    static Scalar concentration(Scalar temperature,
+    static Scalar molarDensity(Scalar temperature,
                                 Scalar pressure)
     { return pressure/(R*temperature); }
+
+    DUMUX_DEPRECATED_MSG("use molarDensity()")
+    static Scalar concentration(Scalar temperature,
+                                Scalar pressure)
+    { return molarDensity(temperature, pressure); }
 };
 
 } // end namepace
-- 
GitLab