From b44a24c9955826554a1c82f4cf15ce85716b7a37 Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Thu, 22 Nov 2018 17:21:59 +0100
Subject: [PATCH] [brine] Remove default for constant salinty

---
 .../material/binarycoefficients/brine_co2.hh  |  4 +--
 dumux/material/components/brine.hh            | 34 ++++++++++---------
 dumux/material/fluidsystems/brineco2.hh       |  8 ++---
 test/material/components/plotproperties.cc    |  4 +++
 .../fluidsystems/test_fluidsystems.cc         |  5 +++
 .../poromechanics/el2p/params.input           |  3 ++
 .../co2/implicit/params.input                 |  2 +-
 .../co2/implicit/paramsni.input               |  2 ++
 8 files changed, 37 insertions(+), 25 deletions(-)

diff --git a/dumux/material/binarycoefficients/brine_co2.hh b/dumux/material/binarycoefficients/brine_co2.hh
index 2ac59a94a7..7d9193a9d2 100644
--- a/dumux/material/binarycoefficients/brine_co2.hh
+++ b/dumux/material/binarycoefficients/brine_co2.hh
@@ -387,7 +387,7 @@ template<class Scalar, class CO2Tables, bool verbose = true>
 class Brine_CO2_Old
 {
     using H2O = Dumux::Components::H2O<Scalar>;
-  using Brine = Dumux::Components::Brine<Scalar,H2O>;
+    using Brine = Dumux::Components::Brine<Scalar,H2O>;
     using CO2 = Dumux::Components::CO2<Scalar, CO2Tables>;
     using IdealGas = Dumux::IdealGas<Scalar>;
 
@@ -419,7 +419,7 @@ public:
         const Scalar Mw = H2O::molarMass(); /* molecular weight of water [kg/mol] */
         const Scalar Ms = 58.8e-3; /* molecular weight of NaCl  [kg/mol] */
 
-        const Scalar X_NaCl = Brine::salinity;
+        const Scalar X_NaCl = Brine::salinity();
         /* salinity: conversion from mass fraction to mole fraction */
         const Scalar x_NaCl = -Mw * X_NaCl / ((Ms - Mw) * X_NaCl - Ms);
 
diff --git a/dumux/material/components/brine.hh b/dumux/material/components/brine.hh
index 25780d2baf..5264cab03f 100644
--- a/dumux/material/components/brine.hh
+++ b/dumux/material/components/brine.hh
@@ -44,6 +44,7 @@ namespace Components {
  * \tparam Scalar The type used for scalar values
  * \tparam H2O Static polymorphism: the Brine class can access all properties of the H2O class
  * \note This is an implementation of brine as a pseudo-component with a constant salinity.
+ * \note the salinity is read from the input file and is a mandatory parameter
  */
 template <class Scalar,
           class H2O_Tabulated = Components::TabulatedComponent<Components::H2O<Scalar>>>
@@ -52,12 +53,10 @@ class Brine
 , public Components::Liquid<Scalar, Brine<Scalar, H2O_Tabulated> >
 , public Components::Gas<Scalar, Brine<Scalar, H2O_Tabulated> >
 {
+    using ThisType = Brine<Scalar, H2O_Tabulated>;
 public:
     using H2O = Components::TabulatedComponent<Dumux::Components::H2O<Scalar>>;
 
-    // The constant salinity
-    static Scalar constantSalinity;
-
     //! The ideal gas constant \f$\mathrm{[J/mol/K]}\f$
     static constexpr Scalar R = Constants<Scalar>::R;
 
@@ -67,15 +66,24 @@ public:
     static std::string name()
     { return "Brine"; }
 
+    /*!
+     * \brief Return the constant salinity
+     */
+    static Scalar salinity()
+    {
+        static const Scalar salinity = getParam<Scalar>("Brine.Salinity");
+        return salinity;
+    }
+
     /*!
      * \brief The molar mass in \f$\mathrm{[kg/mol]}\f$ of brine.
      * This assumes that the salt is pure NaCl.
      */
-   static constexpr Scalar molarMass()
+   static Scalar molarMass()
    {
        const Scalar M1 = H2O::molarMass();
        const Scalar M2 = Components::NaCl<Scalar>::molarMass(); // molar mass of NaCl [kg/mol]
-       return M1*M2/(M2 + constantSalinity*(M1 - M2));
+       return M1*M2/(M2 + ThisType::salinity()*(M1 - M2));
    };
 
     /*!
@@ -113,8 +121,8 @@ public:
         Scalar ps = H2O::vaporPressure(temperature); //Saturation vapor pressure for pure water
         Scalar pi = 0;
         using std::log;
-        if (constantSalinity < 0.26) // here we have hard coded the solubility limit for NaCl
-            pi = (R * temperature * log(1- constantSalinity)); // simplified version of Eq 2.29 in Vishal Jambhekar's Promo
+        if (ThisType::salinity() < 0.26) // here we have hard coded the solubility limit for NaCl
+            pi = (R * temperature * log(1- ThisType::salinity())); // simplified version of Eq 2.29 in Vishal Jambhekar's Promo
         else
             pi = (R * temperature * log(0.74));
         using std::exp;
@@ -164,7 +172,7 @@ public:
         /*Regularization*/
         using std::min;
         using std::max;
-        const Scalar salinity = min(max(constantSalinity,0.0), salSat);
+        const Scalar salinity = min(max(ThisType::salinity(),0.0), salSat);
 
         const Scalar hw = H2O::liquidEnthalpy(T, p)/1E3; /* kJ/kg */
 
@@ -303,7 +311,7 @@ public:
         using std::max;
         const Scalar TempC = temperature - 273.15;
         const Scalar pMPa = pressure/1.0E6;
-        const Scalar salinity = max(0.0, constantSalinity);
+        const Scalar salinity = max(0.0, ThisType::salinity());
 
         const Scalar rhow = H2O::liquidDensity(temperature, pressure);
 
@@ -402,7 +410,7 @@ public:
         // regularisation
         using std::max;
         temperature = max(temperature, 275.0);
-        const Scalar salinity = max(0.0, constantSalinity);
+        const Scalar salinity = max(0.0, ThisType::salinity());
 
         using std::pow;
         using std::exp;
@@ -429,12 +437,6 @@ public:
 template <class Scalar, class H2O>
 struct IsAqueous<Brine<Scalar, H2O>> : public std::true_type {};
 
-/*!
- * \brief Default value for the salinity of the brine (dimensionless).
- */
-template <class Scalar, class H2O>
-Scalar Brine<Scalar, H2O>::constantSalinity = 0.1;
-
 } // end namespace Components
 
 } // end namespace Dumux
diff --git a/dumux/material/fluidsystems/brineco2.hh b/dumux/material/fluidsystems/brineco2.hh
index 17b5e215e3..13e69e40c7 100644
--- a/dumux/material/fluidsystems/brineco2.hh
+++ b/dumux/material/fluidsystems/brineco2.hh
@@ -337,10 +337,6 @@ public:
         std::cout << " - use constant salinity: " << std::boolalpha << Policy::useConstantSalinity() << "\n";
         std::cout << " - use CO2 gas density as gas mixture density: " << std::boolalpha << Policy::useCO2GasDensityAsGasMixtureDensity() << std::endl;
 
-        // maybe set salinity of the constant salinity brine
-        if (useConstantSalinity)
-            ConstantSalinityBrine::constantSalinity = getParam<Scalar>("FluidSystem.Salinity", 0.3);
-
         if (H2O::isTabulated)
             H2O::init(startTemp, endTemp, tempSteps, startPressure, endPressure, pressureSteps);
     }
@@ -474,7 +470,7 @@ public:
             // calulate the equilibrium composition for given T & p
             Scalar xlH2O, xgH2O;
             Scalar xlCO2, xgCO2;
-            const Scalar salinity = useConstantSalinity ? ConstantSalinityBrine::constantSalinity
+            const Scalar salinity = useConstantSalinity ? ConstantSalinityBrine::salinity()
                                                         : fluidState.massFraction(liquidPhaseIdx, NaClIdx);
             Brine_CO2::calculateMoleFractions(T, pl, salinity, /*knownGasPhaseIdx=*/-1, xlCO2, xgH2O);
 
@@ -524,7 +520,7 @@ public:
         Scalar xlCO2;
 
         // calulate the equilibrium composition for given T & p
-        const Scalar salinity = useConstantSalinity ? ConstantSalinityBrine::constantSalinity
+        const Scalar salinity = useConstantSalinity ? ConstantSalinityBrine::salinity()
                                                     : fluidState.massFraction(liquidPhaseIdx, NaClIdx);
         Brine_CO2::calculateMoleFractions(T, p, salinity, /*knowgasPhaseIdx=*/-1, xlCO2, xgH2O);
 
diff --git a/test/material/components/plotproperties.cc b/test/material/components/plotproperties.cc
index 4fcbe0ca4b..dd001166e7 100644
--- a/test/material/components/plotproperties.cc
+++ b/test/material/components/plotproperties.cc
@@ -30,6 +30,7 @@
 #include <vector>
 #include <dumux/common/typetraits/isvalid.hh>
 #include <dumux/common/typetraits/typetraits.hh>
+#include <dumux/common/parameters.hh>
 #include <dumux/io/gnuplotinterface.hh>
 #include <dumux/material/components/air.hh>
 #include <dumux/material/components/benzene.hh>
@@ -339,7 +340,10 @@ int main(int argc, char *argv[])
     else if (compName == "Benzene")
         plotStuff< Components::Benzene<double> >(openPlotWindow);
     else if (compName == "Brine")
+    {
+        Parameters::init([](auto& params){ params["Brine.Salinity"] = "0.1"; });
         plotStuff< Components::Brine<double> >(openPlotWindow);
+    }
     else if (compName == "Calcite")
         plotStuff< Components::Calcite<double> >(openPlotWindow);
     else if (compName == "CalciumIon")
diff --git a/test/material/fluidsystems/test_fluidsystems.cc b/test/material/fluidsystems/test_fluidsystems.cc
index 2f7f8481d7..072ca26010 100644
--- a/test/material/fluidsystems/test_fluidsystems.cc
+++ b/test/material/fluidsystems/test_fluidsystems.cc
@@ -25,6 +25,7 @@
 #include <config.h>
 
 #include "checkfluidsystem.hh"
+#include <dumux/common/parameters.hh>
 
 // include all fluid systems in dumux-stable
 #include <dumux/material/fluidsystems/2pimmiscible.hh>
@@ -146,6 +147,7 @@ int main()
     {   using H2OType = Components::SimpleH2O<Scalar>;
         using FluidSystem = FluidSystems::BrineCO2< Scalar, HeterogeneousCO2Tables::CO2Tables,
                                                     H2OType, FluidSystems::BrineCO2DefaultPolicy</*useConstantSalinity=*/true> >;
+        Parameters::init([](auto& params){ params["Brine.Salinity"] = "0.3"; });
         success += checkFluidSystem<Scalar, FluidSystem>( false ); }
     {   using H2OType = Components::SimpleH2O<Scalar>;
         using FluidSystem = FluidSystems::BrineCO2< Scalar, HeterogeneousCO2Tables::CO2Tables,
@@ -154,6 +156,7 @@ int main()
     {   using H2OType = Components::H2O<Scalar>;
         using FluidSystem = FluidSystems::BrineCO2< Scalar, HeterogeneousCO2Tables::CO2Tables,
                                                     H2OType, FluidSystems::BrineCO2DefaultPolicy</*useConstantSalinity=*/true> >;
+        Parameters::init([](auto& params){ params["Brine.Salinity"] = "0.3"; });
         success += checkFluidSystem<Scalar, FluidSystem>( false ); }
     {   using H2OType = Components::H2O<Scalar>;
         using FluidSystem = FluidSystems::BrineCO2< Scalar, HeterogeneousCO2Tables::CO2Tables,
@@ -162,6 +165,7 @@ int main()
     {   using H2OType = Components::TabulatedComponent<Components::H2O<Scalar>>;
         using FluidSystem = FluidSystems::BrineCO2< Scalar, HeterogeneousCO2Tables::CO2Tables,
                                                     H2OType, FluidSystems::BrineCO2DefaultPolicy</*useConstantSalinity=*/true> >;
+        Parameters::init([](auto& params){ params["Brine.Salinity"] = "0.3"; });
         success += checkFluidSystem<Scalar, FluidSystem>( false ); }
     {   using H2OType = Components::TabulatedComponent<Components::H2O<Scalar>>;
         using FluidSystem = FluidSystems::BrineCO2< Scalar, HeterogeneousCO2Tables::CO2Tables,
@@ -171,6 +175,7 @@ int main()
         using FluidSystem = FluidSystems::BrineCO2< Scalar, HeterogeneousCO2Tables::CO2Tables,
                                                     H2OType, FluidSystems::BrineCO2DefaultPolicy</*useConstantSalinity=*/true,
                                                     /*fastButSimplifiedRelations*/true> >;
+        Parameters::init([](auto& params){ params["Brine.Salinity"] = "0.3"; });
         success += checkFluidSystem<Scalar, FluidSystem>( false ); }
     {   using H2OType = Components::TabulatedComponent<Components::H2O<Scalar>>;
         using FluidSystem = FluidSystems::BrineCO2< Scalar, HeterogeneousCO2Tables::CO2Tables,
diff --git a/test/multidomain/poromechanics/el2p/params.input b/test/multidomain/poromechanics/el2p/params.input
index c6c130e727..b50e7d2159 100644
--- a/test/multidomain/poromechanics/el2p/params.input
+++ b/test/multidomain/poromechanics/el2p/params.input
@@ -29,6 +29,9 @@ MaxIterations = 2000
 [Component]
 SolidDensity = 2650
 
+[Brine]
+Salinity = 0.1
+
 [Newton]
 TargetSteps = 15
 MaxRelativeShift = 1e-6
diff --git a/test/porousmediumflow/co2/implicit/params.input b/test/porousmediumflow/co2/implicit/params.input
index 2360bbfcf4..02d69f446b 100644
--- a/test/porousmediumflow/co2/implicit/params.input
+++ b/test/porousmediumflow/co2/implicit/params.input
@@ -19,7 +19,7 @@ EnableGravity = true
 DepthBOR = 1200# [m] depth below ground surface
 InjectionRate = 1e-4 # [kg/sq/s]
 
-[FluidSystem]
+[Brine]
 Salinity = 1e-1
 
 [LinearSolver]
diff --git a/test/porousmediumflow/co2/implicit/paramsni.input b/test/porousmediumflow/co2/implicit/paramsni.input
index cad66cb663..987cda134e 100644
--- a/test/porousmediumflow/co2/implicit/paramsni.input
+++ b/test/porousmediumflow/co2/implicit/paramsni.input
@@ -12,6 +12,8 @@ PressureLow = 1e5 # [Pa] low end for tabularization of fluid properties
 PressureHigh = 3e7 # [Pa] high end for tabularization of fluid properties
 TemperatureLow = 290.15 # [Pa] low end for tabularization of fluid properties
 TemperatureHigh = 330.15 # [Pa] high end for tabularization of fluid properties
+
+[Brine]
 Salinity = 0.1 # [-] salinity of brine
 
 [Problem]
-- 
GitLab