From c217b2b91e096e4a381735981247cd5842eb0bcd Mon Sep 17 00:00:00 2001
From: Katharina Heck <katharina.heck@iws.uni-stuttgart.de>
Date: Thu, 2 Nov 2017 14:36:43 +0100
Subject: [PATCH] [fix][cleanup] use residualreducition in cc test and cleanup

---
 .../3p/implicit/3pniconductionproblem.hh      |  2 +-
 .../3p/implicit/infiltration3pproblem.hh      | 40 +++----------------
 .../implicit/infiltration3pspatialparams.hh   | 27 +++++++++----
 .../3p/implicit/test_box3p.cc                 |  4 +-
 .../3p/implicit/test_box3p.input              |  5 ++-
 .../porousmediumflow/3p/implicit/test_cc3p.cc |  5 ++-
 .../3p/implicit/test_cc3p.input               |  5 +--
 .../3p/implicit/test_cc3pniconduction.input   |  3 ++
 8 files changed, 40 insertions(+), 51 deletions(-)

diff --git a/test/porousmediumflow/3p/implicit/3pniconductionproblem.hh b/test/porousmediumflow/3p/implicit/3pniconductionproblem.hh
index fd1df77e31..4b7862c146 100644
--- a/test/porousmediumflow/3p/implicit/3pniconductionproblem.hh
+++ b/test/porousmediumflow/3p/implicit/3pniconductionproblem.hh
@@ -71,7 +71,7 @@ SET_TYPE_PROP(ThreePNIConductionProblem,
               SpatialParams,
               ThreePNISpatialParams<TypeTag>);
 
-}
+}// end namespace Properties
 
 
 /*!
diff --git a/test/porousmediumflow/3p/implicit/infiltration3pproblem.hh b/test/porousmediumflow/3p/implicit/infiltration3pproblem.hh
index dc2fc96c02..0d41288f80 100644
--- a/test/porousmediumflow/3p/implicit/infiltration3pproblem.hh
+++ b/test/porousmediumflow/3p/implicit/infiltration3pproblem.hh
@@ -55,10 +55,7 @@ SET_TYPE_PROP(InfiltrationThreePProblem,
               FluidSystem,
               FluidSystems::H2OAirMesitylene<typename GET_PROP_TYPE(TypeTag, Scalar)>);
 
-
-// Maximum tolerated relative error in the Newton method
-SET_SCALAR_PROP(InfiltrationThreePProblem, NewtonMaxRelativeShift, 1e-4);
-}
+}// end namespace Properties
 
 /*!
  * \ingroup ThreePModel
@@ -204,9 +201,7 @@ public:
     {
         BoundaryTypes values;
 
-        if(globalPos[0] > 500. - eps_)
-            values.setAllDirichlet();
-        else if(globalPos[0] < eps_)
+        if(globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_ || globalPos[0] < this->fvGridGeometry().bBoxMin()[0] +eps_)
             values.setAllDirichlet();
         else
             values.setAllNeumann();
@@ -226,30 +221,7 @@ public:
     PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
     {
         PrimaryVariables values(0.0);
-
-        Scalar y = globalPos[1];
-        Scalar x = globalPos[0];
-        Scalar sw, swr=0.12, sgr=0.03;
-
-        if(y > (-1.E-3*x+5) - eps_)
-        {
-            Scalar pc = 9.81 * 1000.0 * (y - (-5E-4*x+5));
-            if (pc < 0.0) pc = 0.0;
-
-            sw = invertPcgw_(pc,
-                             this->spatialParams().materialLawParamsAtPos(globalPos));
-            if (sw < swr) sw = swr;
-            if (sw > 1.-sgr) sw = 1.-sgr;
-
-            values[pressureIdx] = 1e5 ;
-            values[swIdx] = sw;
-            values[snIdx] = 0.;
-        }else {
-            values[pressureIdx] = 1e5 + 9.81 * 1000.0 * ((-5E-4*x+5) - y);
-            values[swIdx] = 1.-sgr;
-            values[snIdx] = 0.;
-        }
-
+        initial_(values, globalPos);
         return values;
     }
  /*!
@@ -265,15 +237,15 @@ public:
      */
     NeumannFluxes neumannAtPos(const GlobalPosition &globalPos) const
     {
-        PrimaryVariables values(0.0);
+        NeumannFluxes values(0.0);
 
         // negative values for injection
         if (time_<2592000.)
         {
-            if ((globalPos[0] <= 175.+eps_) && (globalPos[0] >= 155.-eps_) && (globalPos[1] >= 10.-eps_))
+            if ((globalPos[0] < 175.0+eps_) && (globalPos[0] > 155.0-eps_) && (globalPos[1] > this->fvGridGeometry().bBoxMax()[1]-eps_))
             {
                 values[Indices::contiWEqIdx] = -0.0;
-                values[Indices::contiNEqIdx] = -0.001, // /*Molfluss, umr. über M(Mesit.)=0,120 kg/mol --> 1.2e-4  kg/(sm)
+                values[Indices::contiNEqIdx] = -0.001; // /*Molfluss, umr. über M(Mesit.)=0,120 kg/mol --> 1.2e-4  kg/(sm)
                 values[Indices::contiGEqIdx] = -0.0;
             }
         }
diff --git a/test/porousmediumflow/3p/implicit/infiltration3pspatialparams.hh b/test/porousmediumflow/3p/implicit/infiltration3pspatialparams.hh
index d83cf0b9c8..b2f00a114c 100644
--- a/test/porousmediumflow/3p/implicit/infiltration3pspatialparams.hh
+++ b/test/porousmediumflow/3p/implicit/infiltration3pspatialparams.hh
@@ -74,6 +74,9 @@ class InfiltrationThreePSpatialParams : public ImplicitSpatialParams<TypeTag>
     using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
     using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using Element = typename GridView::template Codim<0>::Entity;
+    using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+    using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag,ElementSolutionVector);
     enum {
         dimWorld=GridView::dimensionworld
     };
@@ -102,15 +105,16 @@ public:
 
         // porosities
         porosity_ = getParam<Scalar>("SpatialParams.porosity");
-
+        vGAlpha_ = getParam<Scalar>("SpatialParams.vanGenuchtenAlpha");
+        vGN_ = getParam<Scalar>("SpatialParams.vanGenuchtenN");
         // residual saturations
         materialParams_.setSwr(0.12);
         materialParams_.setSnr(0.07);
         materialParams_.setSgr(0.03);
 
         // parameters for the 3phase van Genuchten law
-        materialParams_.setVgAlpha(getParam<Scalar>("SpatialParams.vanGenuchtenAlpha"));
-        materialParams_.setVgn(getParam<Scalar>("SpatialParams.vanGenuchtenN"));
+        materialParams_.setVgAlpha(vGAlpha_);
+        materialParams_.setVgn(vGN_);
         materialParams_.setKrRegardsSnr(false);
 
         // parameters for adsorption
@@ -135,14 +139,19 @@ public:
         plotMaterialLaw.plotkr(materialParams_);
     }
 
-    /*!
-     * \brief Returns the scalar intrinsic permeability \f$[m^2]\f$
+      /*!
+     * \brief Function for defining the (intrinsic) permeability \f$[m^2]\f$.
      *
-     * \param globalPos The global position
+     * \param element The element
+     * \param scv The sub control volume
+     * \param elemSol The element solution vector
+     * \return the intrinsic permeability
      */
-    Scalar permeabilityAtPos(const GlobalPosition& globalPos) const
+    PermeabilityType permeability(const Element& element,
+                                  const SubControlVolume& scv,
+                                  const ElementSolutionVector& elemSol) const
     {
-        if (isFineMaterial_(globalPos))
+        if (isFineMaterial_(scv.dofPosition()))
             return fineK_;
         return coarseK_;
     }
@@ -177,6 +186,8 @@ private:
     Scalar coarseK_;
 
     Scalar porosity_;
+    Scalar vGN_;
+    Scalar vGAlpha_;
 
     MaterialLawParams materialParams_;
 
diff --git a/test/porousmediumflow/3p/implicit/test_box3p.cc b/test/porousmediumflow/3p/implicit/test_box3p.cc
index 3185e7d849..311bb5218f 100644
--- a/test/porousmediumflow/3p/implicit/test_box3p.cc
+++ b/test/porousmediumflow/3p/implicit/test_box3p.cc
@@ -157,8 +157,8 @@ int main(int argc, char** argv) try
     auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop);
 
     // the linear solver
-    using LinearSolver = Dumux::UMFPackBackend<TypeTag>;
-    auto linearSolver = std::make_shared<LinearSolver>();
+    using LinearSolver = Dumux::AMGBackend<TypeTag>;
+    auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->elementMapper());
 
     // the non-linear solver
     using NewtonController = typename GET_PROP_TYPE(TypeTag, NewtonController);
diff --git a/test/porousmediumflow/3p/implicit/test_box3p.input b/test/porousmediumflow/3p/implicit/test_box3p.input
index 51b596315f..ec23bdcabc 100644
--- a/test/porousmediumflow/3p/implicit/test_box3p.input
+++ b/test/porousmediumflow/3p/implicit/test_box3p.input
@@ -19,4 +19,7 @@ vanGenuchtenN = 4.0
 PlotFluidMatrixInteractions = false
 
 [LinearSolver]
-ResidualReduction = 1e-6
+ResidualReduction = 1e-13
+
+[Newton]
+MaxRelativeShift = 1e-4
diff --git a/test/porousmediumflow/3p/implicit/test_cc3p.cc b/test/porousmediumflow/3p/implicit/test_cc3p.cc
index d0f712f313..d7215ee9e4 100644
--- a/test/porousmediumflow/3p/implicit/test_cc3p.cc
+++ b/test/porousmediumflow/3p/implicit/test_cc3p.cc
@@ -158,8 +158,8 @@ int main(int argc, char** argv) try
     auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop);
 
     // the linear solver
-    using LinearSolver = Dumux::AMGBackend<TypeTag>;
-    auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->elementMapper());
+    using LinearSolver = UMFPackBackend<TypeTag>;
+    auto linearSolver = std::make_shared<LinearSolver>();
 
     // the non-linear solver
     using NewtonController = typename GET_PROP_TYPE(TypeTag, NewtonController);
@@ -207,6 +207,7 @@ int main(int argc, char** argv) try
 
         // set new dt as suggested by newton controller
         timeLoop->setTimeStepSize(newtonController->suggestTimeStepSize(timeLoop->timeStepSize()));
+         problem->setTime(timeLoop->time()+timeLoop->timeStepSize());
 
     } while (!timeLoop->finished());
 
diff --git a/test/porousmediumflow/3p/implicit/test_cc3p.input b/test/porousmediumflow/3p/implicit/test_cc3p.input
index 9eefdbfd11..5463ca72d5 100644
--- a/test/porousmediumflow/3p/implicit/test_cc3p.input
+++ b/test/porousmediumflow/3p/implicit/test_cc3p.input
@@ -18,6 +18,5 @@ vanGenuchtenN = 4.0
 [Output]
 PlotFluidMatrixInteractions = false
 
-[LinearSolver]
-ResidualReduction = 1e-13
-
+[Newton]
+EnableResidualCriterion = true
diff --git a/test/porousmediumflow/3p/implicit/test_cc3pniconduction.input b/test/porousmediumflow/3p/implicit/test_cc3pniconduction.input
index 11dff4a2c5..9e83ca9efa 100644
--- a/test/porousmediumflow/3p/implicit/test_cc3pniconduction.input
+++ b/test/porousmediumflow/3p/implicit/test_cc3pniconduction.input
@@ -14,3 +14,6 @@ EnableGravity = 0 # disable gravity
 
 [Vtk]
 AddVelocity = 1 #Enable velocity output
+
+[Newton]
+EnableResidualCriterion = true
-- 
GitLab