From a53dd95e3d78ea321c65526b2f1ad35fa0a196f7 Mon Sep 17 00:00:00 2001
From: Katharina Heck <katharina.heck@iws.uni-stuttgart.de>
Date: Wed, 14 Mar 2018 16:37:19 +0100
Subject: [PATCH] [material][fluidstates] add molar density in fluidstate 
 Co-authored-by: Simon Scholz simon.scholz@iws.uni-stuttgart.de 
 Co-authored-by: Beatrix Becker beatrix.becker@iws.uni-stuttgart.de

---
 dumux/material/fluidstates/2p2c.hh            | 16 +++++++--
 dumux/material/fluidstates/compositional.hh   | 19 ++++++----
 dumux/material/fluidstates/immiscible.hh      | 11 +++++-
 .../fluidstates/isothermalimmiscible.hh       | 11 +++++-
 dumux/material/fluidstates/nonequilibrium.hh  | 11 +++++-
 dumux/material/fluidstates/pseudo1p2c.hh      | 36 ++++++++++++-------
 6 files changed, 80 insertions(+), 24 deletions(-)

diff --git a/dumux/material/fluidstates/2p2c.hh b/dumux/material/fluidstates/2p2c.hh
index 91d746422e..27258cb39f 100644
--- a/dumux/material/fluidstates/2p2c.hh
+++ b/dumux/material/fluidstates/2p2c.hh
@@ -120,8 +120,12 @@ public:
     Scalar density(int phaseIdx) const
     { return density_[phaseIdx]; }
 
-    /*!
-     * \brief The dynamic viscosity \f$\mu_\alpha\f$ of fluid phase \f$\alpha\f$ in \f$\mathrm{[Pa s]}\f$
+    /*! @copydoc CompositionalFluidState::molarDensity()
+     */
+    Scalar molarDensity(int phaseIdx) const
+    { return molarDensity_[phaseIdx]; }
+
+    /*! @copydoc CompositionalFluidState::viscosity()
      */
     Scalar viscosity(int phaseIdx) const
     { return viscosity_[phaseIdx]; }
@@ -246,6 +250,13 @@ public:
      */
     void setDensity(int phaseIdx, Scalar value)
     { density_[phaseIdx] = value; }
+
+    /*!
+     * \brief Set the molar density of a phase \f$\mathrm{[mol / m^3]}\f$
+     */
+    void setMolarDensity(int phaseIdx, Scalar value)
+    { molarDensity_[phaseIdx] = value; }
+
     /*!
      * \brief Sets the saturation of a phase.
      * Internally, only the wetting saturation is stored.
@@ -299,6 +310,7 @@ protected:
     Scalar sw_;
     PhaseVector nu_; //phase mass fraction
     Scalar density_[numPhases];
+    Scalar molarDensity_[numPhases];
     Scalar viscosity_[numPhases];
     Scalar massFraction_[numPhases][numComponents];
     Scalar moleFraction_[numPhases][numComponents];
diff --git a/dumux/material/fluidstates/compositional.hh b/dumux/material/fluidstates/compositional.hh
index f1eed0374f..7706a45e1e 100644
--- a/dumux/material/fluidstates/compositional.hh
+++ b/dumux/material/fluidstates/compositional.hh
@@ -194,15 +194,11 @@ public:
     { return density_[phaseIdx]; }
 
     /*!
-     * \brief The molar density \f$\rho_{mol,\alpha}\f$
-     *   of a fluid phase \f$\alpha\f$ in \f$\mathrm{[mol/m^3]}\f$
-     *
-     * The molar density is defined by the mass density \f$\rho_\alpha\f$ and the mean molar mass \f$\overline M_\alpha\f$:
-     *
-     * \f[\rho_{mol,\alpha} = \frac{\rho_\alpha}{\overline M_\alpha} \;.\f]
+     * \brief The molar density \f$\rho_\alpha\f$ of the fluid phase
+     *  \f$\alpha\f$ in \f$\mathrm{[mol/m^3]}\f$
      */
     Scalar molarDensity(int phaseIdx) const
-    { return density_[phaseIdx]/averageMolarMass(phaseIdx); }
+    { return molarDensity_[phaseIdx]; }
 
     /*!
      * \brief The absolute temperature\f$T_\alpha\f$ of a fluid phase \f$\alpha\f$ in \f$\mathrm{[K]}\f$
@@ -298,6 +294,7 @@ public:
             pressure_[phaseIdx] = fs.pressure(phaseIdx);
             saturation_[phaseIdx] = fs.saturation(phaseIdx);
             density_[phaseIdx] = fs.density(phaseIdx);
+            molarDensity_[phaseIdx] = fs.molarDensity(phaseIdx);
             enthalpy_[phaseIdx] = fs.enthalpy(phaseIdx);
             viscosity_[phaseIdx] = fs.viscosity(phaseIdx);
         }
@@ -427,6 +424,12 @@ public:
     void setDensity(int phaseIdx, Scalar value)
     { density_[phaseIdx] = value; }
 
+    /*!
+     * \brief Set the molar density of a phase \f$\mathrm{[mol / m^3]}\f$
+     */
+    void setMolarDensity(int phaseIdx, Scalar value)
+    { molarDensity_[phaseIdx] = value; }
+
     /*!
      * \brief Set the specific enthalpy of a phase\f$\mathrm{[J/kg]}\f$
      */
@@ -465,6 +468,7 @@ public:
             Valgrind::CheckDefined(pressure_[phaseIdx]);
             Valgrind::CheckDefined(saturation_[phaseIdx]);
             Valgrind::CheckDefined(density_[phaseIdx]);
+            Valgrind::CheckDefined(molarDensity_[phaseIdx]);
             Valgrind::CheckDefined(enthalpy_[phaseIdx]);
             Valgrind::CheckDefined(viscosity_[phaseIdx]);
         }
@@ -482,6 +486,7 @@ protected:
     Scalar pressure_[numPhases];
     Scalar saturation_[numPhases];
     Scalar density_[numPhases];
+    Scalar molarDensity_[numPhases];
     Scalar enthalpy_[numPhases];
     Scalar viscosity_[numPhases];
     Scalar temperature_[numPhases];
diff --git a/dumux/material/fluidstates/immiscible.hh b/dumux/material/fluidstates/immiscible.hh
index 4496472eb1..2cb491eeb0 100644
--- a/dumux/material/fluidstates/immiscible.hh
+++ b/dumux/material/fluidstates/immiscible.hh
@@ -205,7 +205,7 @@ public:
      * \f[\rho_{mol,\alpha} = \frac{\rho_\alpha}{\overline M_\alpha} \;.\f]
      */
     Scalar molarDensity(int phaseIdx) const
-    { return density_[phaseIdx]/averageMolarMass(phaseIdx); }
+    { return molarDensity_[phaseIdx]; }
 
     /*!
      * \brief The absolute temperature\f$T_\alpha\f$ of a fluid phase \f$\alpha\f$ in \f$\mathrm{[K]}\f$
@@ -283,6 +283,7 @@ public:
             pressure_[phaseIdx] = fs.pressure(phaseIdx);
             saturation_[phaseIdx] = fs.saturation(phaseIdx);
             density_[phaseIdx] = fs.density(phaseIdx);
+            molarDensity_[phaseIdx] = fs.molarDensity(phaseIdx);
             enthalpy_[phaseIdx] = fs.enthalpy(phaseIdx);
             viscosity_[phaseIdx] = fs.viscosity(phaseIdx);
             temperature_[phaseIdx] = fs.temperature(0);
@@ -323,6 +324,12 @@ public:
     void setDensity(int phaseIdx, Scalar value)
     { density_[phaseIdx] = value; }
 
+    /*!
+     * \brief Set the molar density of a phase \f$\mathrm{[kg/m^3]}\f$
+     */
+    void setMolarDensity(int phaseIdx, Scalar value)
+    { molarDensity_[phaseIdx] = value; }
+
     /*!
      * \brief Set the specific enthalpy of a phase \f$\mathrm{[J/kg]}\f$
      */
@@ -353,6 +360,7 @@ public:
             Valgrind::CheckDefined(pressure_[i]);
             Valgrind::CheckDefined(saturation_[i]);
             Valgrind::CheckDefined(density_[i]);
+            Valgrind::CheckDefined(molarDensity_[i]);
             //Valgrind::CheckDefined(internalEnergy_[i]);
             Valgrind::CheckDefined(viscosity_[i]);
         }
@@ -365,6 +373,7 @@ protected:
     Scalar pressure_[numPhases];
     Scalar saturation_[numPhases];
     Scalar density_[numPhases];
+    Scalar molarDensity_[numPhases];
     Scalar enthalpy_[numPhases];
     Scalar viscosity_[numPhases];
     Scalar temperature_[numPhases];
diff --git a/dumux/material/fluidstates/isothermalimmiscible.hh b/dumux/material/fluidstates/isothermalimmiscible.hh
index ba5b0e7eae..6f41a63c57 100644
--- a/dumux/material/fluidstates/isothermalimmiscible.hh
+++ b/dumux/material/fluidstates/isothermalimmiscible.hh
@@ -173,7 +173,7 @@ public:
      * \f[\rho_{mol,\alpha} = \frac{\rho_\alpha}{\overline M_\alpha} \;.\f]
      */
     Scalar molarDensity(int phaseIdx) const
-    { return density_[phaseIdx]/averageMolarMass(phaseIdx); }
+    { return molarDensity_[phaseIdx]; }
 
     /*!
      * \brief The temperature of a fluid phase \f$\mathrm{[K]}\f$
@@ -234,6 +234,7 @@ public:
             pressure_[phaseIdx] = fs.pressure(phaseIdx);
             saturation_[phaseIdx] = fs.saturation(phaseIdx);
             density_[phaseIdx] = fs.density(phaseIdx);
+            molarDensity_[phaseIdx] = fs.molarDensity(phaseIdx);
             viscosity_[phaseIdx] = fs.viscosity(phaseIdx);
         }
         temperature_ = fs.temperature(0);
@@ -263,6 +264,12 @@ public:
     void setDensity(int phaseIdx, Scalar value)
     { density_[phaseIdx] = value; }
 
+    /*!
+     * \brief Set the molar density of a phase \f$\mathrm{[kg/m^3]}\f$
+     */
+    void setMolarDensity(int phaseIdx, Scalar value)
+    { molarDensity_[phaseIdx] = value; }
+
     /*!
      * \brief Set the dynamic viscosity of a phase \f$\mathrm{[Pa s]}\f$
      */
@@ -284,6 +291,7 @@ public:
             Valgrind::CheckDefined(pressure_[i]);
             Valgrind::CheckDefined(saturation_[i]);
             Valgrind::CheckDefined(density_[i]);
+            Valgrind::CheckDefined(molarDensity_[i]);
             Valgrind::CheckDefined(viscosity_[i]);
         }
 
@@ -295,6 +303,7 @@ protected:
     Scalar pressure_[numPhases];
     Scalar saturation_[numPhases];
     Scalar density_[numPhases];
+    Scalar molarDensity_[numPhases];
     Scalar viscosity_[numPhases];
     Scalar temperature_;
 };
diff --git a/dumux/material/fluidstates/nonequilibrium.hh b/dumux/material/fluidstates/nonequilibrium.hh
index 0fa62dea73..56e578642f 100644
--- a/dumux/material/fluidstates/nonequilibrium.hh
+++ b/dumux/material/fluidstates/nonequilibrium.hh
@@ -205,7 +205,7 @@ public:
      * \f[\rho_{mol,\alpha} = \frac{\rho_\alpha}{\overline M_\alpha} \;.\f]
      */
     Scalar molarDensity(int phaseIdx) const
-    { return density_[phaseIdx]/averageMolarMass(phaseIdx); }
+    { return molarDensity_[phaseIdx]; }
 
 
     /*!
@@ -374,6 +374,12 @@ public:
     void setDensity(int phaseIdx, Scalar value)
     { density_[phaseIdx] = value; }
 
+    /*!
+     * \brief Set the molar density of a phase \f$\mathrm{[kg / m^3]}\f$
+     */
+    void setMolarDensity(int phaseIdx, Scalar value)
+    { molarDensity_[phaseIdx] = value; }
+
     /*!
      * \brief Set the specific enthalpy of a phase \f$\mathrm{[J/m^3]}\f$
      */
@@ -413,6 +419,7 @@ public:
             pressure_[phaseIdx] = fs.pressure(phaseIdx);
             saturation_[phaseIdx] = fs.saturation(phaseIdx);
             density_[phaseIdx] = fs.density(phaseIdx);
+            molarDensity_[phaseIdx] = fs.molarDensity(phaseIdx);
             enthalpy_[phaseIdx] = fs.enthalpy(phaseIdx);
             viscosity_[phaseIdx] = fs.viscosity(phaseIdx);
             temperature_[phaseIdx] = fs.temperature(phaseIdx);
@@ -439,6 +446,7 @@ public:
             Valgrind::CheckDefined(pressure_[i]);
             Valgrind::CheckDefined(saturation_[i]);
             Valgrind::CheckDefined(density_[i]);
+            Valgrind::CheckDefined(molarDensity_[i]);
             Valgrind::CheckDefined(temperature_[i]);
             Valgrind::CheckDefined(enthalpy_[i]);
             Valgrind::CheckDefined(viscosity_[i]);
@@ -455,6 +463,7 @@ protected:
     Scalar pressure_[numPhases];
     Scalar saturation_[numPhases];
     Scalar density_[numPhases];
+    Scalar molarDensity_[numPhases];
     Scalar enthalpy_[numPhases];
     Scalar viscosity_[numPhases];
     Scalar temperature_[numPhases]; //
diff --git a/dumux/material/fluidstates/pseudo1p2c.hh b/dumux/material/fluidstates/pseudo1p2c.hh
index 6a6b716867..16d76d3dcc 100644
--- a/dumux/material/fluidstates/pseudo1p2c.hh
+++ b/dumux/material/fluidstates/pseudo1p2c.hh
@@ -121,18 +121,18 @@ public:
     }
 
     /*!
-     * \brief Returns the mass fraction \f$X^\kappa_\alpha\f$ of component \f$\kappa\f$ in fluid phase \f$\alpha\f$ in \f$\mathrm{[-]}\f$.
-     *
-     * The mass fraction \f$X^\kappa_\alpha\f$ is defined as the weight of all molecules of a
-     * component divided by the total mass of the fluid phase. It is related with the component's mole fraction by means of the relation
-     *
-     * This is either set to 1 or 0 depending on the phase presence for the
-     * non-wetting phase in general.
-     * It is set to the mass fraction of water or 1-massFractionWater
-     * if the considered component is the main component of the wetting phase.
-     *
-     * \param phaseIdx the index of the phase
-     * \param compIdx the index of the component
+     *  @copydoc CompositionalFluidState::molarDensity()
+     */
+    Scalar molarDensity(int phaseIdx) const
+    {
+        if(phaseIdx == presentPhaseIdx_)
+            return molarDensity_;
+        else
+            return 0.;
+    }
+
+    /*!
+     *  @copydoc CompositionalFluidState::massFraction()
      */
     Scalar massFraction(int phaseIdx, int compIdx) const
     {
@@ -291,6 +291,17 @@ public:
         assert(phaseIdx == presentPhaseIdx_);
         density_ = value;
     }
+    /*!
+     * \brief Set the molar density of a phase \f$\mathrm{[mol / m^3]}\f$
+     *
+     * \param phaseIdx the index of the phase
+     * @param value Value to be stored
+     */
+    void setMolarDensity(int phaseIdx, Scalar value)
+    {
+        assert(phaseIdx == presentPhaseIdx_);
+        molarDensity_ = value;
+    }
     /*!
      * \brief Sets the phase Index that is present in this fluidState.
      * @param phaseIdx the index of the phase
@@ -349,6 +360,7 @@ protected:
     Scalar moleFractionWater_;
     Scalar pressure_[numPhases];
     Scalar density_;
+    Scalar molarDensity_;
     Scalar viscosity_;
     Scalar enthalpy_;
     Scalar temperature_;
-- 
GitLab