[remediation] Cleanup implementation. Shorten test.

The input interface got error checking.
The output name is ńow set accordong to the scenario.
Grid file removed and input gridcreator used.
Added EPSILONS (!) for floating point comparisons.
Test time reduced to about 40 seconds.

lecture/mm/remediationscenarios/CMakeLists.txt

@@ -3,9 +3,9 @@ add_input_file_links()
 add_dumux_test(remediationscenariosexercise remediationscenariosexercise remediationscenariosexercise.cc
                python ${dumux_INCLUDE_DIRS}/bin/runtest.py
                  --script fuzzy
-                 --files ${CMAKE_SOURCE_DIR}/lecture/references/kuevette3p3cni-reference.vtu
-                         ${CMAKE_CURRENT_BINARY_DIR}/kuevette3p3cni-00015.vtu
-                 --command "${CMAKE_CURRENT_BINARY_DIR}/remediationscenariosexercise -ParameterFile ${CMAKE_CURRENT_SOURCE_DIR}/remediationscenariosexercise.input -Grid.File ${CMAKE_CURRENT_SOURCE_DIR}/grids/kuevette.dgf -TimeManager.TEnd 6000 -TimeManager.OutputInterval 100")
+                 --files ${CMAKE_SOURCE_DIR}/lecture/references/remediation-steam-reference.vtu
+                         ${CMAKE_CURRENT_BINARY_DIR}/remediation-steam-00006.vtu
+                 --command "${CMAKE_CURRENT_BINARY_DIR}/remediationscenariosexercise -TimeManager.TEnd 300 -TimeManager.OutputInterval 20")
 
 # headers for installation and headercheck
 install(FILES

lecture/mm/remediationscenarios/grids/kuevette.dgf

-DGF
-Interval
-0 0   % first corner
-1.5 0.74   % second corner
-37 19   % 22 cells in x and 12 in y direction
-#
-GridParameter
-overlap 0
-periodic
-closure green
-#GridParameter
-BOUNDARYDOMAIN
-default 1    % all boundaries have id 1
-#BOUNDARYDOMAIN
-# unitcube.dgf

lecture/mm/remediationscenarios/remediationscenariosexercise.input

-
 [TimeManager]
 DtInitial = 1 # [s]
 TEnd = 25920000 # [s] : maybe shorter for some scenarios, can be adapted
+OutputInterval = 25
 
 # Recommended MaxTimeStepSize
 # for steam injection 5 [s]
@@ -9,44 +9,40 @@ TEnd = 25920000 # [s] : maybe shorter for some scenarios, can be adapted
 # for air injection 86400 [s]
 MaxTimeStepSize = 5 # [s]
 
-OutputInterval = 25
-
 [Grid]
-File = ./grids/kuevette.dgf
-
-[Newton]
-UseLineSearch = true
+UpperRight = 1.5 0.74
+Cells = 37 19
 
 [Implicit]
 NumericDifferenceMethod = 1
 
-[]
+[Problem]
 # scenario selection: choose one of the following,
 # all zero means you can use the "free" scenario parameters
-isSteam = 1 # should be default
-isSteamAir = 0
-isAir = 0
+SteamScenario = 1 # should be default
+SteamAirScenario = 0
+AirScenario = 0
 
 # steam scenario, Do not change!
-waterFluxS = -0.3435
-airFluxS = -0.000001
-contaminantFluxS = -0.0000001
-heatFluxS = -15150.0
+WaterFluxS = -0.3435
+AirFluxS = -0.000001
+ContaminantFluxS = -0.0000001
+HeatFluxS = -15150.0
 
 # steamAir scenario, Do not change!
-waterFluxSA = -0.1435
-airFluxSA = -0.2
-contaminantFluxSA = -0.0000001
-heatFluxSA = -6930.0
+WaterFluxSA = -0.1435
+AirFluxSA = -0.2
+ContaminantFluxSA = -0.0000001
+HeatFluxSA = -6930.0
 
 # cool (~20°C) air scenario, Do not change!
-waterFluxA = -0.0082
-airFluxA = -0.3353
-contaminantFluxA = -0.0000001
-heatFluxA = -555.0
+WaterFluxA = -0.0082
+AirFluxA = -0.3353
+ContaminantFluxA = -0.0000001
+HeatFluxA = -555.0
 
 # scenario, free to choose ...
-waterFlux = -0.3325
-airFlux = -0.02
-contaminantFlux = -0.0000001
-heatFlux = -5555.0
+WaterFlux = -0.3325
+AirFlux = -0.02
+ContaminantFlux = -0.0000001
+HeatFlux = -5555.0

lecture/mm/remediationscenarios/remediationscenariosproblem.hh

@@ -26,6 +26,8 @@
 #ifndef DUMUX_REMEDIATIONSCENARIOS_PROBLEM_HH
 #define DUMUX_REMEDIATIONSCENARIOS_PROBLEM_HH
 
+#include <dune/common/exceptions.hh>
+
 #include <dumux/material/fluidsystems/h2oairmesitylene.hh>
 
 #include <dumux/porousmediumflow/3p3c/implicit/model.hh>
@@ -62,7 +64,7 @@ SET_BOOL_PROP(KuevetteProblem, ProblemEnableGravity, true);
 SET_INT_PROP(KuevetteProblem, ImplicitNumericDifferenceMethod, 1);
 
 // Set the maximum time step
-SET_SCALAR_PROP(KuevetteProblem, TimeManagerMaxTimeStepSize, 25.);
+SET_SCALAR_PROP(KuevetteProblem, TimeManagerMaxTimeStepSize, 25.0);
 
 // set newton relative tolerance
 SET_SCALAR_PROP(KuevetteProblem, NewtonMaxRelativeShift, 1e-4);
@@ -155,9 +157,51 @@ public:
      * \param gridView The grid view
      */
     KuevetteProblem(TimeManager &timeManager, const GridView &gridView)
-        : ParentType(timeManager, gridView), eps_(1e-6)
+        : ParentType(timeManager, gridView)
     {
         FluidSystem::init();
+
+        int isSteam = GET_RUNTIME_PARAM(TypeTag, bool, Problem.SteamScenario) ? 1 : 0;
+        int isSteamAir = GET_RUNTIME_PARAM(TypeTag, bool, Problem.SteamAirScenario) ? 1 : 0;
+        int isAir = GET_RUNTIME_PARAM(TypeTag, bool, Problem.AirScenario) ? 1 : 0;
+
+        // error checking
+        if (isSteam + isSteamAir + isAir > 1)
+            DUNE_THROW(Dune::InvalidStateException, "Please choose a single scenario!");
+
+        if (isSteam == 1)
+        {
+            waterFlux_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.WaterFluxS); // 0.3435 [mol/(s m)] in total
+            airFlux_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.AirFluxS);
+            contaminantFlux_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.ContaminantFluxS);
+            heatFlux_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.HeatFluxS);
+            name_ = "remediation-steam";
+        }
+        else if (isSteamAir == 1)
+        {
+            waterFlux_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.WaterFluxSA); // 0.3435 [mol/(s m)] in total
+            airFlux_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.AirFluxSA);
+            contaminantFlux_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.ContaminantFluxSA);
+            heatFlux_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.HeatFluxSA);
+            name_ = "remediation-steamair";
+        }
+        else if (isAir == 1)
+        {
+            waterFlux_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.WaterFluxA); // 0.3435 [mol/(s m)] in total
+            airFlux_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.AirFluxA);
+            contaminantFlux_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.ContaminantFluxA);
+            heatFlux_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.HeatFluxA);
+            name_ = "remediation-air";
+        }
+        else
+        {
+            // free scenario
+            waterFlux_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.WaterFlux); // 0.3435 [mol/(s m)] in total
+            airFlux_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.AirFlux);
+            contaminantFlux_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.ContaminantFlux);
+            heatFlux_ = GET_RUNTIME_PARAM(TypeTag, Scalar, Problem.HeatFlux);
+            name_ = "remediation";
+        }
     }
 
     /*!
@@ -170,8 +214,8 @@ public:
      *
      * This is used as a prefix for files generated by the simulation.
      */
-    const char *name() const
-    { return "kuevette3p3cni"; }
+    const std::string name() const
+    { return name_; }
 
     void sourceAtPos(PrimaryVariables &values,
                      const GlobalPosition &globalPos) const
@@ -245,45 +289,19 @@ public:
         const GlobalPosition &globalPos = element.geometry().corner(scvIdx);
 
         // negative values for injection
-
-        if (globalPos[0] > 0.01)
+        if (globalPos[0] > 0.01 + eps_)
         {
-            values[Indices::contiWEqIdx] = 0.;
-            values[Indices::contiGEqIdx] = 0.;
-            values[Indices::contiNEqIdx] = 0.;
-            values[Indices::energyEqIdx] = 0.;
+            values[Indices::contiWEqIdx] = 0.0;
+            values[Indices::contiGEqIdx] = 0.0;
+            values[Indices::contiNEqIdx] = 0.0;
+            values[Indices::energyEqIdx] = 0.0;
         }
         else
         {
-            bool isSteam, isSteamAir, isAir;
-            isSteam = GET_RUNTIME_PARAM(TypeTag, Scalar, isSteam);
-            isSteamAir = GET_RUNTIME_PARAM(TypeTag, Scalar, isSteamAir);
-            isAir = GET_RUNTIME_PARAM(TypeTag, Scalar, isAir);
-
-            if (isSteam){
-            values[Indices::contiWEqIdx] = GET_RUNTIME_PARAM(TypeTag, Scalar, waterFluxS); // 0.3435 [mol/(s m)] in total
-            values[Indices::contiGEqIdx] = GET_RUNTIME_PARAM(TypeTag, Scalar, airFluxS);
-            values[Indices::contiNEqIdx] = GET_RUNTIME_PARAM(TypeTag, Scalar, contaminantFluxS);
-            values[Indices::energyEqIdx] = GET_RUNTIME_PARAM(TypeTag, Scalar, heatFluxS);
-            } else
-            if (isSteamAir){
-            values[Indices::contiWEqIdx] = GET_RUNTIME_PARAM(TypeTag, Scalar, waterFluxSA); // 0.3435 [mol/(s m)] in total
-            values[Indices::contiGEqIdx] = GET_RUNTIME_PARAM(TypeTag, Scalar, airFluxSA);
-            values[Indices::contiNEqIdx] = GET_RUNTIME_PARAM(TypeTag, Scalar, contaminantFluxSA);
-            values[Indices::energyEqIdx] = GET_RUNTIME_PARAM(TypeTag, Scalar, heatFluxSA);
-            } else
-            if (isAir){
-            values[Indices::contiWEqIdx] = GET_RUNTIME_PARAM(TypeTag, Scalar, waterFluxA); // 0.3435 [mol/(s m)] in total
-            values[Indices::contiGEqIdx] = GET_RUNTIME_PARAM(TypeTag, Scalar, airFluxA);
-            values[Indices::contiNEqIdx] = GET_RUNTIME_PARAM(TypeTag, Scalar, contaminantFluxA);
-            values[Indices::energyEqIdx] = GET_RUNTIME_PARAM(TypeTag, Scalar, heatFluxA);
-            } else {  // free scenario
-            values[Indices::contiWEqIdx] = GET_RUNTIME_PARAM(TypeTag, Scalar, waterFlux); // 0.3435 [mol/(s m)] in total
-            values[Indices::contiGEqIdx] = GET_RUNTIME_PARAM(TypeTag, Scalar, airFlux);
-            values[Indices::contiNEqIdx] = GET_RUNTIME_PARAM(TypeTag, Scalar, contaminantFlux);
-            values[Indices::energyEqIdx] = GET_RUNTIME_PARAM(TypeTag, Scalar, heatFlux);
-            }
-
+            values[Indices::contiWEqIdx] = waterFlux_;
+            values[Indices::contiGEqIdx] = airFlux_;
+            values[Indices::contiNEqIdx] = contaminantFlux_;
+            values[Indices::energyEqIdx] = heatFlux_;
         }
     }
 
@@ -323,24 +341,41 @@ public:
      * \param globalIdx The index of the global vertex
      * \param globalPos The global position
      */
-    int initialPhasePresence(const Vertex &vert,
+    int initialPhasePresence(const Vertex &vertex,
                              const int &globalIdx,
                              const GlobalPosition &globalPos) const
     {
-        if((globalPos[0] >= 0.20) && (globalPos[0] <= 0.80) && (globalPos[1] >= 0.4) && (globalPos[1] <= 0.65))
+        if((globalPos[0] > 0.20 - eps_) && (globalPos[0] < 0.80 + eps_) &&
+           (globalPos[1] > 0.4 - eps_) && (globalPos[1] < 0.65 + eps_))
         {
             return threePhases;
         }
-        else return wgPhaseOnly;
+        else
+        {
+            return wgPhaseOnly;
+        }
     }
 
     bool shouldWriteOutput() const
     {
 
-    int outputInterval = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, TimeManager, OutputInterval);
+        int outputInterval = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, TimeManager, OutputInterval);
+
+        return (this->timeManager().timeStepIndex() % outputInterval == 0 ||
+                this->timeManager().willBeFinished());
+    }
 
-    return (this->timeManager().timeStepIndex() % outputInterval == 0 ||
-            this->timeManager().willBeFinished());
+    /*!
+     * \brief Returns true if a restart file should be written to
+     *        disk.
+     *
+     * The default behavior is to write one restart file every 5 time
+     * steps. This file is intended to be overwritten by the
+     * implementation.
+     */
+    bool shouldWriteRestartFile() const
+    {
+        return false;
     }
 
 
@@ -355,14 +390,18 @@ private:
         values[switch2Idx] = 1.e-6;
         values[temperatureIdx] = 293.0;
 
-        if((globalPos[0] >= 0.20) && (globalPos[0] <= 0.80) && (globalPos[1] >= 0.4) && (globalPos[1] <= 0.65))
+        if((globalPos[0] > 0.20 - eps_) && (globalPos[0] < 0.80 + eps_) &&
+           (globalPos[1] > 0.4 - eps_) && (globalPos[1] < 0.65 + eps_))
         {
             values[switch2Idx] = 0.07;
         }
     }
 
-    const Scalar eps_;
+    static constexpr Scalar eps_ = 1.5e-7;
+    Scalar waterFlux_, airFlux_, contaminantFlux_, heatFlux_;
+    std::string name_;
 };
+
 } //end namespace
 
 #endif

lecture/mm/remediationscenarios/remediationscenariosspatialparams.hh

@@ -136,26 +136,12 @@ public:
         fineMaterialParams_.setKrRegardsSnr(true);
 
         // parameters for adsorption
-        coarseMaterialParams_.setKdNAPL(0.);
-        coarseMaterialParams_.setRhoBulk(1500.);
-        fineMaterialParams_.setKdNAPL(0.);
-        fineMaterialParams_.setRhoBulk(1500.);
+        coarseMaterialParams_.setKdNAPL(0.0);
+        coarseMaterialParams_.setRhoBulk(1500);
+        fineMaterialParams_.setKdNAPL(0.0);
+        fineMaterialParams_.setRhoBulk(1500);
     }
 
-    ~KuevetteSpatialParams()
-    {}
-
-
-    /*!
-     * \brief Update the spatial parameters with the flow solution
-     *        after a timestep.
-     *
-     * \param globalSolution The global solution vector
-     */
-    void update(const SolutionVector &globalSolution)
-    {
-    };
-
     /*!
      * \brief Apply the intrinsic permeability tensor to a pressure
      *        potential gradient.
@@ -262,9 +248,15 @@ public:
 private:
     bool isFineMaterial_(const GlobalPosition &pos) const
     {
-        if (0.13 <= pos[0] && 1.20 >= pos[0] && 0.32 <= pos[1] && pos[1] <= 0.57)
+        if (0.13 < pos[0] + eps_ && 1.20 > pos[0] - eps_ &&
+            0.32 < pos[1] + eps_ && pos[1] < 0.57 - eps_)
+        {
             return true;
-        else return false;
+        }
+        else
+        {
+            return false;
+        }
     };
 
     Scalar fineK_;
@@ -280,6 +272,8 @@ private:
     MaterialLawParams coarseMaterialParams_;
 
     Scalar lambdaSolid_;
+
+    static constexpr Scalar eps_ = 1.5e-7;
 };
 
 }