From 6b7cad65d53cc7360bb139d1251030cddc7f6ff7 Mon Sep 17 00:00:00 2001
From: Johannes Hommel <johannes.hommel@iws.uni-stuttgart.de>
Date: Wed, 31 Jan 2018 12:04:01 +0100
Subject: [PATCH] [2pncmin][poro-permLaws] renamed initial to reference values

---
 .../permeabilitykozenycarman.hh               |  6 +-
 .../porosityprecipitation.hh                  |  2 +-
 .../2pncmin/implicit/dissolutionproblem.hh    | 12 ++--
 .../implicit/dissolutionspatialparams.hh      | 72 ++++++++++---------
 .../2pncmin/implicit/test_2pncmin.input       |  7 +-
 5 files changed, 50 insertions(+), 49 deletions(-)

diff --git a/dumux/material/fluidmatrixinteractions/permeabilitykozenycarman.hh b/dumux/material/fluidmatrixinteractions/permeabilitykozenycarman.hh
index 2e6842c304..3650c15766 100644
--- a/dumux/material/fluidmatrixinteractions/permeabilitykozenycarman.hh
+++ b/dumux/material/fluidmatrixinteractions/permeabilitykozenycarman.hh
@@ -76,12 +76,12 @@ public:
                                   const SubControlVolume& scv,
                                   const ElementSolution& elemSol) const
     {
-        auto initPoro = spatialParams().initialPorosity(element, scv);
+        auto refPoro = spatialParams().referencePorosity(element, scv);
         auto poro = spatialParams().porosity(element, scv, elemSol);
 
         using std::pow;
-        auto factor = pow((1.0 - initPoro)/(1.0 - poro), 2) * pow(poro/initPoro, 3);
-        return applyFactorToPermeability_(spatialParams().initialPermeability(element, scv), factor);
+        auto factor = pow((1.0 - refPoro)/(1.0 - poro), 2) * pow(poro/refPoro, 3);
+        return applyFactorToPermeability_(spatialParams().referencePermeability(element, scv), factor);
     }
 
 private:
diff --git a/dumux/material/fluidmatrixinteractions/porosityprecipitation.hh b/dumux/material/fluidmatrixinteractions/porosityprecipitation.hh
index 27edea92b9..4440d126d4 100644
--- a/dumux/material/fluidmatrixinteractions/porosityprecipitation.hh
+++ b/dumux/material/fluidmatrixinteractions/porosityprecipitation.hh
@@ -76,7 +76,7 @@ public:
         auto minPoro = spatialParams_().minPorosity(element, scv);
 
         using std::max;
-        return max(minPoro, spatialParams_().initialPorosity(element, scv) - sumPrecipitates);
+        return max(minPoro, spatialParams_().referencePorosity(element, scv) - sumPrecipitates);
     }
 
 private:
diff --git a/test/porousmediumflow/2pncmin/implicit/dissolutionproblem.hh b/test/porousmediumflow/2pncmin/implicit/dissolutionproblem.hh
index 2779900a94..376abd22eb 100644
--- a/test/porousmediumflow/2pncmin/implicit/dissolutionproblem.hh
+++ b/test/porousmediumflow/2pncmin/implicit/dissolutionproblem.hh
@@ -144,8 +144,7 @@ public:
         temperature_            = getParam<Scalar>("Problem.Temperature");
         reservoirPressure_      = getParam<Scalar>("Problem.ReservoirPressure");
         initLiqSaturation_      = getParam<Scalar>("Problem.LiquidSaturation");
-        initPrecipitatedSalt1_  = getParam<Scalar>("Problem.InitPrecipitatedSalt1");
-        initPrecipitatedSalt2_  = getParam<Scalar>("Problem.InitPrecipitatedSalt2");
+        initPrecipitatedSaltBlock_  = getParam<Scalar>("Problem.InitPrecipitatedSaltBlock");
 
         outerLiqSaturation_     = getParam<Scalar>("Problem.OuterLiqSaturation");
         innerLiqSaturation_     = getParam<Scalar>("Problem.InnerLiqSaturation");
@@ -284,9 +283,9 @@ public:
         priVars[switchIdx]   = initLiqSaturation_;                 // Sw primary variable
         priVars[xwNaClIdx]   = massToMoleFrac_(outerSalinity_);     // mole fraction
         if(globalPos[0] > 5.0 - eps_ && globalPos[0] < 19.0 + eps_)
-            priVars[precipNaClIdx] = initPrecipitatedSalt2_; // [kg/m^3]
+            priVars[precipNaClIdx] = initPrecipitatedSaltBlock_; // [kg/m^3]
         else
-            priVars[precipNaClIdx] = initPrecipitatedSalt1_; // [kg/m^3]
+            priVars[precipNaClIdx] = 0.0; // [kg/m^3]
 
         return priVars;
     }
@@ -350,7 +349,7 @@ public:
         if (precipSalt*timeStepSize_ + volVars.precipitateVolumeFraction(sPhaseIdx)* volVars.molarDensity(sPhaseIdx)< 0)
             precipSalt = -volVars.precipitateVolumeFraction(sPhaseIdx)* volVars.molarDensity(sPhaseIdx)/timeStepSize_;
 
-        if (volVars.precipitateVolumeFraction(sPhaseIdx) >= this->spatialParams().initialPorosity(element, scv) - saltPorosity  && precipSalt > 0)
+        if (volVars.precipitateVolumeFraction(sPhaseIdx) >= this->spatialParams().referencePorosity(element, scv) - saltPorosity  && precipSalt > 0)
             precipSalt = 0;
 
         source[conti0EqIdx + NaClIdx] += -precipSalt;
@@ -425,8 +424,7 @@ private:
     Scalar initLiqSaturation_;
     Scalar outerLiqSaturation_;
     Scalar innerLiqSaturation_;
-    Scalar initPrecipitatedSalt1_;
-    Scalar initPrecipitatedSalt2_;
+    Scalar initPrecipitatedSaltBlock_;
     Scalar innerSalinity_;
     Scalar time_ = 0.0;
     Scalar timeStepSize_ = 0.0;
diff --git a/test/porousmediumflow/2pncmin/implicit/dissolutionspatialparams.hh b/test/porousmediumflow/2pncmin/implicit/dissolutionspatialparams.hh
index 30585172c1..a8ff76e143 100644
--- a/test/porousmediumflow/2pncmin/implicit/dissolutionspatialparams.hh
+++ b/test/porousmediumflow/2pncmin/implicit/dissolutionspatialparams.hh
@@ -72,12 +72,14 @@ class DissolutionSpatialparams : public FVSpatialParams<TypeTag>
     using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
     using MaterialLawParams = typename GET_PROP_TYPE(TypeTag, MaterialLaw)::Params;
+    using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
     using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
     using CoordScalar = typename GridView::ctype;
     enum { dimWorld=GridView::dimensionworld };
 
     using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
     using Tensor = Dune::FieldMatrix<CoordScalar, dimWorld, dimWorld>;
+    using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
     using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry)::LocalView;
     using SubControlVolume = typename FVElementGeometry::SubControlVolume;
     using Element = typename GridView::template Codim<0>::Entity;
@@ -93,8 +95,8 @@ public:
     : ParentType(problem)
     {
         solubilityLimit_     = getParam<Scalar>("SpatialParams.SolubilityLimit", 0.26);
-        initialPorosity_     = getParam<Scalar>("SpatialParams.Porosity", 0.11);
-        initialPermeability_ = getParam<Scalar>("SpatialParams.Permeability", 2.23e-14);
+        referencePorosity_     = getParam<Scalar>("SpatialParams.referencePorosity", 0.11);
+        refPermeability_ = getParam<Scalar>("SpatialParams.referencePermeability", 2.23e-14);
         irreducibleLiqSat_   = getParam<Scalar>("SpatialParams.IrreducibleLiqSat", 0.2);
         irreducibleGasSat_   = getParam<Scalar>("SpatialParams.IrreducibleGasSat", 1e-3);
         pEntry1_             = getParam<Scalar>("SpatialParams.Pentry1", 500);
@@ -108,28 +110,13 @@ public:
         materialParams_.setPe(pEntry1_);
         materialParams_.setLambda(bcLambda1_);
 
-        // set main diagonal entries of the permeability tensor to a value
-        // setting to one value means: isotropic, homogeneous
-        for (int i = 0; i < dimWorld; i++)  //TODO make this nice and dependend on PermeabilityType!
-            initK_[i][i] = initialPermeability_;
-
         //! Intitialize the parameter laws
         poroLaw_.init(*this);
         permLaw_.init(*this);
-    }
 
-    /*! Intrinsic permeability tensor K \f$[m^2]\f$ depending
-     *  on the position in the domain
-     *
-     *  \param element The finite volume element
-     *  \param scv The sub-control volume
-     *
-     *  Solution dependent permeability function
-     */
-    PermeabilityType permeability(const Element& element,
-                        const SubControlVolume& scv,
-                        const ElementSolutionVector& elemSol) const
-    { return permLaw_.evaluatePermeability(element, scv, elemSol); }
+        for (int i = 0; i < dimWorld; i++)
+            referencePermeability_[i][i] = refPermeability_;
+    }
 
     /*!
      *  \brief Define the minimum porosity \f$[-]\f$ distribution
@@ -141,34 +128,50 @@ public:
     { return 1e-5; }
 
     /*!
-     *  \brief Define the initial porosity \f$[-]\f$ distribution
+     *  \brief Define the reference porosity \f$[-]\f$ distribution.
+     *  This is the porosity of the porous medium without any of the
+     *  considered solid phases.
      *
      *  \param element The finite element
      *  \param scv The sub-control volume
      */
-    Scalar initialPorosity(const Element& element, const SubControlVolume &scv) const
-    { return initialPorosity_; }
+    Scalar referencePorosity(const Element& element, const SubControlVolume &scv) const
+    { return referencePorosity_; }
 
     /*!
-     *  \brief Define the initial permeability \f$[m^2]\f$ distribution
+     *  \brief Return the actual recent porosity \f$[-]\f$ accounting for
+     *  clogging caused by mineralization
      *
      *  \param element The finite element
      *  \param scv The sub-control volume
      */
-    PermeabilityType initialPermeability(const Element& element, const SubControlVolume &scv) const
-    { return initK_; }
+    Scalar porosity(const Element& element,
+                    const SubControlVolume& scv,
+                    const ElementSolutionVector& elemSol) const
+    { return poroLaw_.evaluatePorosity(element, scv, elemSol); }
 
     /*!
-     *  \brief Define the minimum porosity \f$[-]\f$ after clogging caused by mineralization
+     *  \brief Define the reference permeability \f$[m^2]\f$ distribution
      *
      *  \param element The finite element
      *  \param scv The sub-control volume
      */
-    Scalar porosity(const Element& element,
-                    const SubControlVolume& scv,
-                    const ElementSolutionVector& elemSol) const
-    { return poroLaw_.evaluatePorosity(element, scv, elemSol); }
+    PermeabilityType referencePermeability(const Element& element, const SubControlVolume &scv) const
+    { return referencePermeability_; }
+
 
+    /*! Intrinsic permeability tensor K \f$[m^2]\f$ depending
+     *  on the position in the domain
+     *
+     *  \param element The finite volume element
+     *  \param scv The sub-control volume
+     *
+     *  Solution dependent permeability function
+     */
+    PermeabilityType permeability(const Element& element,
+                        const SubControlVolume& scv,
+                        const ElementSolutionVector& elemSol) const
+    { return permLaw_.evaluatePermeability(element, scv, elemSol); }
 
     Scalar solidity(const SubControlVolume &scv) const
     { return 1.0 - porosityAtPos(scv.center()); }
@@ -184,14 +187,15 @@ public:
     { return materialParams_; }
 
 private:
+
     MaterialLawParams materialParams_;
 
     PorosityLaw poroLaw_;
     PermeabilityLaw permLaw_;
     Scalar solubilityLimit_;
-    Scalar initialPorosity_;
-    Scalar initialPermeability_;
-    PermeabilityType initK_= 0.0;
+    Scalar referencePorosity_;
+    Scalar refPermeability_;
+    PermeabilityType referencePermeability_ = 0.0;
     Scalar irreducibleLiqSat_;
     Scalar irreducibleGasSat_;
     Scalar pEntry1_;
diff --git a/test/porousmediumflow/2pncmin/implicit/test_2pncmin.input b/test/porousmediumflow/2pncmin/implicit/test_2pncmin.input
index 04e398d6a5..0fd2372265 100644
--- a/test/porousmediumflow/2pncmin/implicit/test_2pncmin.input
+++ b/test/porousmediumflow/2pncmin/implicit/test_2pncmin.input
@@ -31,13 +31,12 @@ OuterPressure       = 11e6         # [Pa] reservoir boundary pressure
 OuterLiqSaturation  = 0.2          # [-]  liquid saturation at outer boundary
 OuterSalinity       = 0.4          # [-]  Initial salinity
 LiquidSaturation    = 0.4          # [-]  initial liquid saturation
-InitPrecipitatedSalt1    = 0.0     # [-]  initial precipitated salt
-InitPrecipitatedSalt2    = 0.05    # [-]  initial precipitated salt in the blocked part
+InitPrecipitatedSaltBlock    = 0.05    # [-]  initial precipitated salt in the blocked part
 
 [SpatialParams]
 SolubilityLimit     = 0.26 	       # [-]  solubility limit of salt in brine
-Porosity            = 0.11         # [-]  initial porosity
-Permeability        = 2.23e-14
+referencePorosity   = 0.11         # [-]  initial porosity
+referencePermeability   = 2.23e-14
 IrreducibleLiqSat   = 0.2          # [-]  irreducible liquid saturation
 IrreducibleGasSat   = 0.001        # [-]  irreducible gas saturation
 Pentry1             = 500          # [Pa]
-- 
GitLab