Adapt to new material laws

Fixes #21 (closed)

lecture/mm/buckleyleverett/buckleyleverettanalyticsolution.hh

@@ -21,8 +21,7 @@
 
 #include <dumux/common/properties.hh>
 #include <dumux/porousmediumflow/2p/sequential/properties.hh>
-#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
-#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/materiallaw.hh>
 
 /*!
  * \file
@@ -36,33 +35,6 @@ namespace Dumux {
  * the Buckley-Leverett problem
  */
 
-template<typename Scalar, typename Law>
-struct CheckMaterialLaw
-{
-    static bool isLinear()
-    {
-        return false;
-    }
-};
-
-template<typename Scalar>
-struct CheckMaterialLaw<Scalar, LinearMaterial<Scalar> >
-{
-    static bool isLinear()
-    {
-        return true;
-    }
-};
-
-template<typename Scalar>
-struct CheckMaterialLaw<Scalar, EffToAbsLaw< LinearMaterial<Scalar> > >
-{
-    static bool isLinear()
-    {
-        return true;
-    }
-};
-
 template<class TypeTag> class BuckleyLeverettAnalytic
 {
     using Problem = GetPropType<TypeTag, Properties::Problem>;
@@ -70,8 +42,6 @@ template<class TypeTag> class BuckleyLeverettAnalytic
     using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
     using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
-    using MaterialLaw = typename SpatialParams::MaterialLaw;
-    using MaterialLawParams = typename MaterialLaw::Params;
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
     using FluidState = GetPropType<TypeTag, Properties::FluidState>;
     using Indices = GetPropType<TypeTag, Properties::Indices>;
@@ -87,7 +57,6 @@ template<class TypeTag> class BuckleyLeverettAnalytic
     using GlobalPosition = Dune::FieldVector<Scalar, dimworld>;
 
 public:
-
     // functions needed for analytical solution
 
     void initializeAnalytic()
@@ -106,11 +75,11 @@ public:
      */
     void prepareAnalytic()
     {
-        const MaterialLawParams& materialLawParams(problem_.spatialParams().materialLawParams(dummyElement_));
-
-        swr_ = materialLawParams.swr();
-        snr_ = materialLawParams.snr();
-        porosity_ = problem_.spatialParams().porosity(dummyElement_);
+        const auto& dummyElement = *problem_.gridView().template begin<0>();
+        const auto fluidMatrixInteraction = problem_.spatialParams().fluidMatrixInteractionAtPos(dummyElement.geometry().center());
+        swr_ = fluidMatrixInteraction.effToAbsParams().swr();
+        snr_ = fluidMatrixInteraction.effToAbsParams().snr();
+        porosity_ = problem_.spatialParams().porosity(dummyElement);
         time_ = 0;
         satVec_ = swr_;
 
@@ -128,12 +97,11 @@ public:
 
         for (int i = 0; i < pointNum_; i++)
         {
-            fractionalW_[i] = MaterialLaw::krw(materialLawParams, satVec_[i])/viscosityW;
-            fractionalW_[i] /= (fractionalW_[i] + MaterialLaw::krn(materialLawParams, satVec_[i])/viscosityNW);
+            fractionalW_[i] = fluidMatrixInteraction.krw(satVec_[i])/viscosityW;
+            fractionalW_[i] /= (fractionalW_[i] + fluidMatrixInteraction.krn(satVec_[i])/viscosityNW);
         }
 
         dfwdsw_ = 0;
-
         for (int i = 1; i < intervalNum_; i++)
         {
             dfwdsw_[i] = (fractionalW_[i + 1] - fractionalW_[i - 1]) / (satVec_[i + 1] - satVec_[i - 1]);
@@ -177,12 +145,13 @@ public:
         {
             xf_[i] = vTot_ * time_ / porosity_ * dfwdsw_[i];
         }
+
+        // position of maximum xf_
+        int xfmax = 0;
         int xhelp = pointNum_ / 3;
         int xhelpold = 0;
         int xhelpoldold = 0;
-        int xfmax = 0;
-
-        // position of maximum xf_
+        int xhelp2 = 0;
         for (int i = 0; i < pointNum_; i++)
         {
             if (xf_[i] > xf_[i + 1])
@@ -197,13 +166,9 @@ public:
         Scalar A1;
         Scalar A2;
         Scalar b;
-        int xhelp2 = 0;
 
         while (a)
         {
-            if (CheckMaterialLaw<Scalar, MaterialLaw>::isLinear())
-                break;
-
             A1 = 0;
 
             for (int i = 0; i <= xhelp - 1; i++)
@@ -262,7 +227,7 @@ public:
             int index = problem_.variables().index(*eIt);
 
             // account for linear material law
-            if (CheckMaterialLaw<Scalar, MaterialLaw>::isLinear())
+            if constexpr (SpatialParams::pcSwCurveIsLinear())
             {
                 if (globalPos[0] <= xf_[1])
                 {
@@ -313,9 +278,7 @@ public:
     }
 
     BlockVector AnalyticSolution() const
-    {
-        return analyticSolution_;
-    }
+    { return analyticSolution_; }
 
     //Write saturation and pressure into file
     template<class MultiWriter>
@@ -332,9 +295,13 @@ public:
     }
 
     //! Construct an IMPES object.
-    BuckleyLeverettAnalytic(Problem& problem, Scalar totalVelocity = 3e-7) :
-        problem_(problem), analyticSolution_(0), error_(0), elementVolume_(0), size_(problem.gridView().size(0)), vTot_(totalVelocity), dummyElement_(
-                *(problem_.gridView().template begin<0> ()))
+    BuckleyLeverettAnalytic(Problem& problem, Scalar totalVelocity = 3e-7)
+        : problem_(problem)
+        , analyticSolution_(0)
+        , error_(0)
+        , elementVolume_(0)
+        , size_(problem.gridView().size(0))
+        , vTot_(totalVelocity)
     {
         dummyGlobal_ = 0.0;
         dummyGlobal_[0] = 1.0;
@@ -365,7 +332,6 @@ private:
     Dune::FieldVector<Scalar, pointNum_> dfwdsw_;
     Dune::FieldVector<Scalar, pointNum_> xf_;
     int dfwdswmax_;
-    const Element& dummyElement_;
     GlobalPosition dummyGlobal_;
 
 };