From 9e8571936b95c20aa6e575fc89376dda6092d733 Mon Sep 17 00:00:00 2001
From: Katharina Heck <katharina.heck@iws.uni-stuttgart.de>
Date: Thu, 8 Nov 2018 11:58:02 +0100
Subject: [PATCH] [mpnc][spatialparams] free spatialparams of dependency on
 fluidsystem, introducse wettingphase in spatialparams as in 2p model

---
 .../fluidmatrixinteractions/mp/2padapter.hh   | 12 ++++++----
 .../porousmediumflow/mpnc/volumevariables.hh  | 14 +++++++----
 .../mpnc/implicit/2p2ccomparison/problem.hh   |  6 ++---
 .../implicit/2p2ccomparison/spatialparams.hh  | 23 ++++++++++++------
 .../mpnc/implicit/kinetic/problem.hh          |  8 +++----
 .../mpnc/implicit/kinetic/spatialparams.hh    | 22 ++++++++++++-----
 .../mpnc/implicit/obstacle/problem.hh         |  6 ++---
 .../mpnc/implicit/obstacle/spatialparams.hh   | 24 ++++++++++++-------
 .../implicit/thermalnonequilibrium/problem.hh |  6 ++---
 .../thermalnonequilibrium/spatialparams.hh    | 22 ++++++++++++-----
 10 files changed, 94 insertions(+), 49 deletions(-)

diff --git a/dumux/material/fluidmatrixinteractions/mp/2padapter.hh b/dumux/material/fluidmatrixinteractions/mp/2padapter.hh
index 08f35e8479..d4f6bdb01e 100644
--- a/dumux/material/fluidmatrixinteractions/mp/2padapter.hh
+++ b/dumux/material/fluidmatrixinteractions/mp/2padapter.hh
@@ -40,11 +40,9 @@ namespace Dumux
  *
  * \sa MpBrookscoreyMaterialParams
  */
-template <int wPhaseIdx, class TwoPLaw >
+template <class TwoPLaw>
 class TwoPAdapter
 {
-    enum { nPhaseIdx = (wPhaseIdx == 0)?1:0 };
-
 public:
     using Params = typename TwoPLaw::Params;
     enum { numPhases = 2 };
@@ -58,8 +56,10 @@ public:
     template <class ContainerT, class FluidState>
     static void capillaryPressures(ContainerT &values,
                                    const Params &params,
-                                   const FluidState &state)
+                                   const FluidState &state,
+                                   const int wPhaseIdx)
     {
+        const int nPhaseIdx = 1 - wPhaseIdx;
         // non-wetting phase gets the capillary pressure added
         values[nPhaseIdx] = 0;
 
@@ -76,8 +76,10 @@ public:
     template <class ContainerT, class FluidState>
     static void relativePermeabilities(ContainerT &values,
                                        const Params &params,
-                                       const FluidState &state)
+                                       const FluidState &state,
+                                       const int wPhaseIdx)
     {
+        const int nPhaseIdx = 1 - wPhaseIdx;
         values[wPhaseIdx] = TwoPLaw::krw(params, state.saturation(wPhaseIdx));
         values[nPhaseIdx] = TwoPLaw::krn(params, state.saturation(wPhaseIdx));
     }
diff --git a/dumux/porousmediumflow/mpnc/volumevariables.hh b/dumux/porousmediumflow/mpnc/volumevariables.hh
index 5d3a4e3459..3f3d3f41a4 100644
--- a/dumux/porousmediumflow/mpnc/volumevariables.hh
+++ b/dumux/porousmediumflow/mpnc/volumevariables.hh
@@ -106,11 +106,13 @@ public:
 
         //calculate the remaining quantities
         const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
+        const int wPhaseIdx = problem.spatialParams().template wettingPhase<FluidSystem>(element, scv, elemSol);
         // relative permeabilities
         using MaterialLaw = typename Problem::SpatialParams::MaterialLaw;
         MaterialLaw::relativePermeabilities(relativePermeability_,
                                             materialParams,
-                                            fluidState_);
+                                            fluidState_,
+                                            wPhaseIdx);
         typename FluidSystem::ParameterCache paramCache;
         paramCache.updateAll(fluidState_);
         if (enableDiffusion)
@@ -180,12 +182,13 @@ public:
         // set the phase pressures
         /////////////
         // capillary pressure parameters
+        const int wPhaseIdx = problem.spatialParams().template wettingPhase<FluidSystem>(element, scv, elemSol);
         const auto& materialParams =
             problem.spatialParams().materialLawParams(element, scv, elemSol);
         // capillary pressures
         Scalar capPress[numPhases()];
         using MaterialLaw = typename Problem::SpatialParams::MaterialLaw;
-        MaterialLaw::capillaryPressures(capPress, materialParams, fluidState);
+        MaterialLaw::capillaryPressures(capPress, materialParams, fluidState, wPhaseIdx);
         // add to the pressure of the first fluid phase
 
         // depending on which pressure is stored in the primary variables
@@ -558,10 +561,12 @@ public:
         const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
 
         // relative permeabilities
+        const int wPhaseIdx = problem.spatialParams().template wettingPhase<FluidSystem>(element, scv, elemSol);
         using MaterialLaw = typename Problem::SpatialParams::MaterialLaw;
         MaterialLaw::relativePermeabilities(relativePermeability_,
                                             materialParams,
-                                            fluidState_);
+                                            fluidState_,
+                                            wPhaseIdx);
         typename FluidSystem::ParameterCache paramCache;
         paramCache.updateAll(fluidState_);
         if (enableDiffusion)
@@ -632,10 +637,11 @@ public:
         // capillary pressure parameters
         const auto& materialParams =
             problem.spatialParams().materialLawParams(element, scv, elemSol);
+        const int wPhaseIdx = problem.spatialParams().template wettingPhase<FluidSystem>(element, scv, elemSol);
         // capillary pressures
         Scalar capPress[numPhases()];
         using MaterialLaw = typename Problem::SpatialParams::MaterialLaw;
-        MaterialLaw::capillaryPressures(capPress, materialParams, fluidState);
+        MaterialLaw::capillaryPressures(capPress, materialParams, fluidState, wPhaseIdx);
         // add to the pressure of the first fluid phase
 
         // depending on which pressure is stored in the primary variables
diff --git a/test/porousmediumflow/mpnc/implicit/2p2ccomparison/problem.hh b/test/porousmediumflow/mpnc/implicit/2p2ccomparison/problem.hh
index 43fd63f3ac..924c7def94 100644
--- a/test/porousmediumflow/mpnc/implicit/2p2ccomparison/problem.hh
+++ b/test/porousmediumflow/mpnc/implicit/2p2ccomparison/problem.hh
@@ -71,8 +71,7 @@ struct SpatialParams<TypeTag, TTag::MPNCComparison>
 {
     using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
-    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
-    using type = MPNCComparisonSpatialParams<FVGridGeometry, Scalar, FluidSystem>;
+    using type = MPNCComparisonSpatialParams<FVGridGeometry, Scalar>;
 };
 
 // Set fluid configuration
@@ -283,8 +282,9 @@ private:
         const auto& matParams =
             this->spatialParams().materialLawParamsAtPos(globalPos);
         PhaseVector pc;
+        const int wPhaseIdx = this->spatialParams().template wettingPhaseAtPos<FluidSystem>(globalPos);
         using MaterialLaw = typename ParentType::SpatialParams::MaterialLaw;
-        MaterialLaw::capillaryPressures(pc, matParams, fs);
+        MaterialLaw::capillaryPressures(pc, matParams, fs, wPhaseIdx);
         fs.setPressure(liquidPhaseIdx,
                        fs.pressure(gasPhaseIdx) + pc[liquidPhaseIdx] - pc[gasPhaseIdx]);
 
diff --git a/test/porousmediumflow/mpnc/implicit/2p2ccomparison/spatialparams.hh b/test/porousmediumflow/mpnc/implicit/2p2ccomparison/spatialparams.hh
index c52f0452d8..d23513d008 100644
--- a/test/porousmediumflow/mpnc/implicit/2p2ccomparison/spatialparams.hh
+++ b/test/porousmediumflow/mpnc/implicit/2p2ccomparison/spatialparams.hh
@@ -42,10 +42,10 @@ namespace Dumux {
  * \brief Definition of the spatial params properties for the obstacle problem
  *
  */
-template<class FVGridGeometry, class Scalar, class FluidSystem>
+template<class FVGridGeometry, class Scalar>
 class MPNCComparisonSpatialParams
 : public FVSpatialParams<FVGridGeometry, Scalar,
-                         MPNCComparisonSpatialParams<FVGridGeometry, Scalar, FluidSystem>>
+                         MPNCComparisonSpatialParams<FVGridGeometry, Scalar>>
 {
     using GridView = typename FVGridGeometry::GridView;
     using FVElementGeometry = typename FVGridGeometry::LocalView;
@@ -53,17 +53,14 @@ class MPNCComparisonSpatialParams
 
     using Element = typename GridView::template Codim<0>::Entity;
     using ParentType = FVSpatialParams<FVGridGeometry, Scalar,
-                                       MPNCComparisonSpatialParams<FVGridGeometry, Scalar, FluidSystem>>;
+                                       MPNCComparisonSpatialParams<FVGridGeometry, Scalar>>;
 
     using GlobalPosition = typename SubControlVolume::GlobalPosition;
-
-    enum {liquidPhaseIdx = FluidSystem::liquidPhaseIdx};
-
     using EffectiveLaw = RegularizedBrooksCorey<Scalar>;
 
 public:
     using PermeabilityType = Scalar;
-    using MaterialLaw = TwoPAdapter<liquidPhaseIdx, EffToAbsLaw<EffectiveLaw>>;
+    using MaterialLaw = TwoPAdapter<EffToAbsLaw<EffectiveLaw>>;
     using MaterialLawParams = typename MaterialLaw::Params;
 
     //! The constructor
@@ -127,6 +124,18 @@ public:
             return coarseMaterialParams_;
     }
 
+    /*!
+     * \brief Function for defining which phase is to be considered as the wetting phase.
+     *
+     * \return the wetting phase index
+     * \param globalPos The global position
+     */
+    template<class FluidSystem>
+    int wettingPhaseAtPos(const GlobalPosition& globalPos) const
+    {
+        return FluidSystem::phase0Idx;
+    }
+
 private:
     /*!
      * \brief Returns whether a given global position is in the
diff --git a/test/porousmediumflow/mpnc/implicit/kinetic/problem.hh b/test/porousmediumflow/mpnc/implicit/kinetic/problem.hh
index 35e0ffa78b..a063879c14 100644
--- a/test/porousmediumflow/mpnc/implicit/kinetic/problem.hh
+++ b/test/porousmediumflow/mpnc/implicit/kinetic/problem.hh
@@ -122,8 +122,7 @@ struct SpatialParams<TypeTag, TTag::EvaporationAtmosphere>
 {
     using FVGridGeometry = GetPropType<TypeTag, FVGridGeometry>;
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
-    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
-    using type = EvaporationAtmosphereSpatialParams<FVGridGeometry, Scalar, FluidSystem>;
+    using type = EvaporationAtmosphereSpatialParams<FVGridGeometry, Scalar>;
 };
 
 // Set the interfacial area relation: wetting -- non-wetting
@@ -458,10 +457,11 @@ private:
         const auto &materialParams =
             this->spatialParams().materialLawParamsAtPos(globalPos);
         Scalar capPress[numPhases];
-
+        //get the index for the wettingphase
+        const int wPhaseIdx = this->spatialParams().template wettingPhaseAtPos<FluidSystem>(globalPos);
         //obtain pc according to saturation
         using MaterialLaw = typename ParentType::SpatialParams::MaterialLaw;
-        MaterialLaw::capillaryPressures(capPress, materialParams, equilibriumFluidState);
+        MaterialLaw::capillaryPressures(capPress, materialParams, equilibriumFluidState, wPhaseIdx);
 
         Scalar p[numPhases];
         if (this->spatialParams().inPM_(globalPos)){
diff --git a/test/porousmediumflow/mpnc/implicit/kinetic/spatialparams.hh b/test/porousmediumflow/mpnc/implicit/kinetic/spatialparams.hh
index 0d65390fc2..075a567500 100644
--- a/test/porousmediumflow/mpnc/implicit/kinetic/spatialparams.hh
+++ b/test/porousmediumflow/mpnc/implicit/kinetic/spatialparams.hh
@@ -47,28 +47,26 @@ namespace Dumux {
 /**
  * \brief Definition of the spatial parameters for the evaporation atmosphere Problem (using a "poor man's coupling")
  */
-template<class FVGridGeometry, class Scalar, class FluidSystem>
+template<class FVGridGeometry, class Scalar>
 class EvaporationAtmosphereSpatialParams
 : public FVNonEquilibriumSpatialParams<FVGridGeometry, Scalar,
-                                       EvaporationAtmosphereSpatialParams<FVGridGeometry, Scalar, FluidSystem>>
+                                       EvaporationAtmosphereSpatialParams<FVGridGeometry, Scalar>>
 {
     using GridView = typename FVGridGeometry::GridView;
     using FVElementGeometry = typename FVGridGeometry::LocalView;
     using SubControlVolume = typename FVElementGeometry::SubControlVolume;
     using Element = typename GridView::template Codim<0>::Entity;
-    using ThisType = EvaporationAtmosphereSpatialParams<FVGridGeometry, Scalar, FluidSystem>;
+    using ThisType = EvaporationAtmosphereSpatialParams<FVGridGeometry, Scalar>;
     using ParentType = FVNonEquilibriumSpatialParams<FVGridGeometry, Scalar, ThisType>;
 
     using GlobalPosition = Dune::FieldVector<Scalar, GridView::dimension>;
 
     enum { dimWorld = GridView::dimensionworld };
-    enum { numPhases = FluidSystem::numPhases };
-    enum { liquidPhaseIdx   = FluidSystem::liquidPhaseIdx };
 public:
     //! export the type used for the permeability
     using PermeabilityType = Scalar;
     //! export the material law type used
-    using MaterialLaw = TwoPAdapter<liquidPhaseIdx, EffToAbsLaw<RegularizedBrooksCorey<Scalar>>>;
+    using MaterialLaw = TwoPAdapter<EffToAbsLaw<RegularizedBrooksCorey<Scalar>>>;
     //! convenience aliases of the law parameters
     using MaterialLawParams = typename MaterialLaw::Params;
 
@@ -303,6 +301,18 @@ public:
         else DUNE_THROW(Dune::InvalidStateException, "You should not be here: x=" << globalPos[0] << " y= "<< globalPos[dimWorld-1]);
     }
 
+    /*!
+     * \brief Function for defining which phase is to be considered as the wetting phase.
+     *
+     * \return the wetting phase index
+     * \param globalPos The global position
+     */
+    template<class FluidSystem>
+    int wettingPhaseAtPos(const GlobalPosition& globalPos) const
+    {
+        return FluidSystem::phase0Idx;
+    }
+
     /*!\brief Give back whether the tested position (input) is a specific region (porous medium part) in the domain
      *
      * This setting ensures, that the boundary between the two domains has porous medium properties.
diff --git a/test/porousmediumflow/mpnc/implicit/obstacle/problem.hh b/test/porousmediumflow/mpnc/implicit/obstacle/problem.hh
index 1ebb273a5d..1ce66e19b9 100644
--- a/test/porousmediumflow/mpnc/implicit/obstacle/problem.hh
+++ b/test/porousmediumflow/mpnc/implicit/obstacle/problem.hh
@@ -74,8 +74,7 @@ struct SpatialParams<TypeTag, TTag::Obstacle>
 {
     using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
-    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
-    using type = ObstacleSpatialParams<FVGridGeometry, Scalar, FluidSystem>;
+    using type = ObstacleSpatialParams<FVGridGeometry, Scalar>;
 };
 
 // Set fluid configuration
@@ -361,8 +360,9 @@ private:
         // calulate the capillary pressure
         const auto& matParams = this->spatialParams().materialLawParamsAtPos(globalPos);
         PhaseVector pc;
+        const int wPhaseIdx = this->spatialParams().template wettingPhaseAtPos<FluidSystem>(globalPos);
         using MaterialLaw = typename ParentType::SpatialParams::MaterialLaw;
-        MaterialLaw::capillaryPressures(pc, matParams, fs);
+        MaterialLaw::capillaryPressures(pc, matParams, fs, wPhaseIdx);
         fs.setPressure(otherPhaseIdx,
                        fs.pressure(refPhaseIdx)
                        + (pc[otherPhaseIdx] - pc[refPhaseIdx]));
diff --git a/test/porousmediumflow/mpnc/implicit/obstacle/spatialparams.hh b/test/porousmediumflow/mpnc/implicit/obstacle/spatialparams.hh
index e8a5b7831c..772af32182 100644
--- a/test/porousmediumflow/mpnc/implicit/obstacle/spatialparams.hh
+++ b/test/porousmediumflow/mpnc/implicit/obstacle/spatialparams.hh
@@ -30,7 +30,6 @@
 #include <dumux/material/fluidmatrixinteractions/2p/regularizedlinearmaterial.hh>
 #include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
 
-#include <dumux/material/fluidmatrixinteractions/mp/mplinearmaterial.hh>
 #include <dumux/material/fluidmatrixinteractions/mp/2padapter.hh>
 
 namespace Dumux {
@@ -41,30 +40,27 @@ namespace Dumux {
  * \brief Definition of the spatial params properties for the obstacle problem
  *
  */
-template<class FVGridGeometry, class Scalar, class FluidSystem>
+template<class FVGridGeometry, class Scalar>
 class ObstacleSpatialParams
 : public FVSpatialParams<FVGridGeometry, Scalar,
-                         ObstacleSpatialParams<FVGridGeometry, Scalar, FluidSystem>>
+                         ObstacleSpatialParams<FVGridGeometry, Scalar>>
 {
     using GridView = typename FVGridGeometry::GridView;
     using FVElementGeometry = typename FVGridGeometry::LocalView;
     using SubControlVolume = typename FVElementGeometry::SubControlVolume;
     using Element = typename GridView::template Codim<0>::Entity;
     using ParentType = FVSpatialParams<FVGridGeometry, Scalar,
-                                       ObstacleSpatialParams<FVGridGeometry, Scalar, FluidSystem>>;
+                                       ObstacleSpatialParams<FVGridGeometry, Scalar>>;
 
     enum {dimWorld=GridView::dimensionworld};
     using GlobalPosition = typename SubControlVolume::GlobalPosition;
-
-    enum {liquidPhaseIdx = FluidSystem::liquidPhaseIdx};
-
     using EffectiveLaw = RegularizedLinearMaterial<Scalar>;
 
 public:
     //! export the type used for the permeability
     using PermeabilityType = Scalar;
     //! export the material law type used
-    using MaterialLaw = TwoPAdapter<liquidPhaseIdx, EffToAbsLaw<EffectiveLaw>>;
+    using MaterialLaw = TwoPAdapter<EffToAbsLaw<EffectiveLaw>>;
     using MaterialLawParams = typename MaterialLaw::Params;
 
     //! the constructor
@@ -127,6 +123,18 @@ public:
             return coarseMaterialParams_;
     }
 
+    /*!
+     * \brief Function for defining which phase is to be considered as the wetting phase.
+     *
+     * \return the wetting phase index
+     * \param globalPos The global position
+     */
+    template<class FluidSystem>
+    int wettingPhaseAtPos(const GlobalPosition& globalPos) const
+    {
+        return FluidSystem::phase0Idx;
+    }
+
 private:
     /*!
      * \brief Returns whether a given global position is in the
diff --git a/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/problem.hh b/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/problem.hh
index 2ef8a197bd..cb424c04fc 100644
--- a/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/problem.hh
+++ b/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/problem.hh
@@ -80,8 +80,7 @@ struct SpatialParams<TypeTag, TTag::CombustionOneComponent>
 {
     using FVGridGeometry = GetPropType<TypeTag, FVGridGeometry>;
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
-        using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
-    using type = CombustionSpatialParams<FVGridGeometry, Scalar, FluidSystem>;
+    using type = CombustionSpatialParams<FVGridGeometry, Scalar>;
 };
 
 template<class TypeTag>
@@ -455,8 +454,9 @@ private:
         //obtain pc according to saturation
         const auto &materialParams =
         this->spatialParams().materialLawParamsAtPos(globalPos);
+        const int wPhaseIdx = this->spatialParams().template wettingPhaseAtPos<FluidSystem>(globalPos);
         using MaterialLaw = typename ParentType::SpatialParams::MaterialLaw;
-         MaterialLaw::capillaryPressures(capPress, materialParams, fluidState);
+        MaterialLaw::capillaryPressures(capPress, materialParams, fluidState, wPhaseIdx);
 
         Scalar p[numPhases];
 
diff --git a/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/spatialparams.hh b/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/spatialparams.hh
index a0427f6411..27cbc72f14 100644
--- a/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/spatialparams.hh
+++ b/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/spatialparams.hh
@@ -41,16 +41,16 @@ namespace Dumux {
  * \brief Definition of the spatial parameters for the one component combustion problem
  *
  */
-template<class FVGridGeometry, class Scalar, class FluidSystem>
+template<class FVGridGeometry, class Scalar>
 class CombustionSpatialParams
 : public FVNonEquilibriumSpatialParams<FVGridGeometry, Scalar,
-                                       CombustionSpatialParams<FVGridGeometry, Scalar, FluidSystem>>
+                                       CombustionSpatialParams<FVGridGeometry, Scalar>>
 {
     using GridView = typename FVGridGeometry::GridView;
     using FVElementGeometry = typename FVGridGeometry::LocalView;
     using SubControlVolume = typename FVElementGeometry::SubControlVolume;
     using Element = typename GridView::template Codim<0>::Entity;
-    using ThisType = CombustionSpatialParams<FVGridGeometry, Scalar, FluidSystem>;
+    using ThisType = CombustionSpatialParams<FVGridGeometry, Scalar>;
     using ParentType = FVNonEquilibriumSpatialParams<FVGridGeometry, Scalar, ThisType>;
 
     enum {dimWorld = GridView::dimensionworld};
@@ -58,13 +58,11 @@ class CombustionSpatialParams
 
     using EffectiveLaw = HeatPipeLaw<Scalar>;
 
-    enum {wPhaseIdx = FluidSystem::wPhaseIdx};
-
 public:
     //! export the type used for the permeability
     using PermeabilityType = Scalar;
     //! export the material law type used
-    using MaterialLaw = TwoPAdapter<wPhaseIdx, EffToAbsLaw<EffectiveLaw>>;
+    using MaterialLaw = TwoPAdapter<EffToAbsLaw<EffectiveLaw>>;
     using MaterialLawParams = typename MaterialLaw::Params;
     using FluidSolidInterfacialAreaFormulation = FluidSolidInterfacialAreaShiWang<Scalar>;
 
@@ -159,6 +157,18 @@ public:
         }
     }
 
+    /*!
+     * \brief Function for defining which phase is to be considered as the wetting phase.
+     *
+     * \return the wetting phase index
+     * \param globalPos The global position
+     */
+    template<class FluidSystem>
+    int wettingPhaseAtPos(const GlobalPosition& globalPos) const
+    {
+        return FluidSystem::phase0Idx;
+    }
+
     /*!
      * \brief Return a reference to the material parameters of the material law.
      * \param globalPos The position in global coordinates. */
-- 
GitLab