From e2de77eb1537fbcc5702666ece1d89bff1370d08 Mon Sep 17 00:00:00 2001
From: Gabriele Seitz <gabriele.seitz@iws.uni-stuttgart.de>
Date: Tue, 17 Oct 2017 11:53:26 +0200
Subject: [PATCH] [1pncmin] cleanup and add cc-testproblem

---
 .../1pncmin/implicit/model.hh                 | 10 ++-
 .../1pncmin/implicit/propertydefaults.hh      | 10 +--
 .../1pncmin/implicit/CMakeLists.txt           |  3 +-
 .../1pncmin/implicit/test_box1pncmin.input    | 71 +++++++++++++++++++
 .../1pncmin/implicit/test_cc1pncmin.cc        | 55 ++++++++++++++
 .../1pncmin/implicit/test_cc1pncmin.input     | 71 +++++++++++++++++++
 .../1pncmin/implicit/thermochemproblem.hh     | 17 +++--
 .../implicit/thermochemspatialparams.hh       | 10 +--
 8 files changed, 222 insertions(+), 25 deletions(-)
 create mode 100644 test/porousmediumflow/1pncmin/implicit/test_box1pncmin.input
 create mode 100644 test/porousmediumflow/1pncmin/implicit/test_cc1pncmin.cc
 create mode 100644 test/porousmediumflow/1pncmin/implicit/test_cc1pncmin.input

diff --git a/dumux/porousmediumflow/1pncmin/implicit/model.hh b/dumux/porousmediumflow/1pncmin/implicit/model.hh
index 1b1cbd79fe..c6eedb7f03 100644
--- a/dumux/porousmediumflow/1pncmin/implicit/model.hh
+++ b/dumux/porousmediumflow/1pncmin/implicit/model.hh
@@ -126,7 +126,6 @@ class OnePNCMinModel: public OnePNCModel<TypeTag>
         dim = GridView::dimension,
         dimWorld = GridView::dimensionworld,
 
-
         numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
         numSPhases = GET_PROP_VALUE(TypeTag, NumSPhases),
         numComponents = GET_PROP_VALUE(TypeTag, NumComponents),
@@ -161,6 +160,15 @@ public:
        // register standardized vtk output fields
         auto& vtkOutputModule = problem.vtkOutputModule();
 
+        vtkOutputModule.addSecondaryVariable("Kxx",
+                                             [this](const VolumeVariables& v){ return this->perm_(v.permeability())[0][0]; });
+        if (dim >= 2)
+            vtkOutputModule.addSecondaryVariable("Kyy",
+                                                 [this](const VolumeVariables& v){ return this->perm_(v.permeability())[1][1]; });
+        if (dim >= 3)
+            vtkOutputModule.addSecondaryVariable("Kzz",
+                                                 [this](const VolumeVariables& v){ return this->perm_(v.permeability())[2][2]; });
+
         vtkOutputModule.addSecondaryVariable("permeabilityFactor", [](const VolumeVariables& v)
                                              { return v.permeabilityFactor(); });
 
diff --git a/dumux/porousmediumflow/1pncmin/implicit/propertydefaults.hh b/dumux/porousmediumflow/1pncmin/implicit/propertydefaults.hh
index 47ed91b534..4770a7ea7a 100644
--- a/dumux/porousmediumflow/1pncmin/implicit/propertydefaults.hh
+++ b/dumux/porousmediumflow/1pncmin/implicit/propertydefaults.hh
@@ -33,7 +33,7 @@
 #include "volumevariables.hh"
 #include "properties.hh"
 
-#include <dumux/porousmediumflow/compositional/localresidual.hh>
+// #include <dumux/porousmediumflow/compositional/localresidual.hh>
 #include <dumux/porousmediumflow/nonisothermal/implicit/propertydefaults.hh>
 #include <dumux/porousmediumflow/1pnc/implicit/newtoncontroller.hh>
 #include <dumux/porousmediumflow/implicit/darcyfluxvariables.hh>
@@ -80,9 +80,9 @@ public:
 };
 
 //physical processes to be considered by the isothermal model
-SET_BOOL_PROP(OnePNCMin, EnableAdvection, true);
-SET_BOOL_PROP(OnePNCMin, EnableMolecularDiffusion, true);
-SET_BOOL_PROP(OnePNCMin, EnableEnergyBalance, false);
+// SET_BOOL_PROP(OnePNCMin, EnableAdvection, true);
+// SET_BOOL_PROP(OnePNCMin, EnableMolecularDiffusion, true);
+// SET_BOOL_PROP(OnePNCMin, EnableEnergyBalance, false);
 
 /*!
  * \brief Set the property for the number of equations.
@@ -135,7 +135,7 @@ SET_TYPE_PROP(OnePNCMin, Indices, OnePNCMinIndices <TypeTag, /*PVOffset=*/0>);
 //        Actually the Forchheimer coefficient is also a function of the dimensions of the
 //        porous medium. Taking it as a constant is only a first approximation
 //        (Nield, Bejan, Convection in porous media, 2006, p. 10)
-SET_SCALAR_PROP(OnePNCMin, SpatialParamsForchCoeff, 0.55);
+// SET_SCALAR_PROP(OnePNCMin, SpatialParamsForchCoeff, 0.55);
 
 //set isothermal VolumeVariables
 SET_TYPE_PROP(OnePNCMin, IsothermalVolumeVariables, OnePNCMinVolumeVariables<TypeTag>);
diff --git a/test/porousmediumflow/1pncmin/implicit/CMakeLists.txt b/test/porousmediumflow/1pncmin/implicit/CMakeLists.txt
index d4c9a5d1c9..91d858c04b 100644
--- a/test/porousmediumflow/1pncmin/implicit/CMakeLists.txt
+++ b/test/porousmediumflow/1pncmin/implicit/CMakeLists.txt
@@ -3,8 +3,9 @@
 add_input_file_links()
 
 add_executable(test_box1pncmin test_box1pncmin.cc)
+add_executable(test_cc1pncmin test_cc1pncmin.cc)
 
-#add_executable_all(test_box1p2c test_box1p2c.cc)
+# add_executable_all(test_box1p2c test_box1p2c.cc)
 #
 # add_executable_all(test_cc1p2c test_cc1p2c.cc)
 
diff --git a/test/porousmediumflow/1pncmin/implicit/test_box1pncmin.input b/test/porousmediumflow/1pncmin/implicit/test_box1pncmin.input
new file mode 100644
index 0000000000..2f267580bd
--- /dev/null
+++ b/test/porousmediumflow/1pncmin/implicit/test_box1pncmin.input
@@ -0,0 +1,71 @@
+# Parameter file for test case 1pncmin.
+# Everything behind a '#' is a comment.
+# Type "./test_1pncmin --help" for more information.
+
+[TimeManager]
+DtInitial = 1                      # [s] initial time step size
+MaxTimeStepSize = 50               # [s] maximum time step size
+TEnd= 5000                        # [s] duration of the simulation
+FreqOutput          = 10           # frequency of VTK output
+WriteRestartFile    = 1            # Boolean. Should restart files be written? (1) Yes (0) No
+
+[Grid]
+UpperRight = 0.08 0.01 #20 10 #     # [m] upper right corner coordinates
+Cells =  40 2 # 21 6                # [-] number of cells in x,y-direction
+
+[FluidSystem]
+NTemperature = 10                   # [-] number of tabularization entries
+NPressure = 100                     # [-] number of tabularization entries
+PressureLow = 1E5                   # [Pa]low end for tabularization of fluid properties
+PressureHigh = 1E6                  # [Pa]high end for tabularization of fluid properties
+TemperatureLow = 373.15             # [Pa]low end for tabularization of fluid properties
+TemperatureHigh = 873.15            # [Pa]high end for tabularization of fluid properties
+
+[Problem]
+Name =  box1pncmin
+IsCharge = 1                       # Bool: 1: charge; 0: discharge
+
+[Charge]
+PressureInitial = 1E5              # [Pa] Initial reservoir pressure
+TemperatureInitial = 773.15        # [K]  reservoir temperature
+VaporInitial = 0.01                # [-]  initial mole fraction of water
+CaOInitial = 0.0                   # [-]  molefraction in the solid phase;  0 dehydration/charge, 1  hydration/discharge
+CaO2H2Initial = 0.3960             # [-]  molefraction in the solid phase;  0 dehydration/charge, 1  hydration/discharge
+PressureIn = 1.02e5                # [Pa] Inlet pressure; charge
+PressureOut = 1e5                  # [Pa] outlet pressure
+TemperatureIn = 773.15             # [K] inlet temperature: // Shao 500 °C
+TemperatureOut = 673.15            # [K] outlet temperature: // Shao noflow
+InFlow = 5                         # [mol/s] Inflow of HTF // Shao 0.309 g/s (Area (0.015)^2pi m^2) --> here 1.55 mol/s
+VaporIn = 0.01                     # [] molefraction // Shao 0.01 [g/g]
+
+[Discharge]
+PressureInitial = 1E5              # [Pa] Initial reservoir pressure
+TemperatureInitial = 573.15        # [K]  reservoir temperature
+VaporInitial = 0.35                # [-]  initial mole fraction of water
+CaOInitial = 0.2                   # [-]  molefraction in the solid phase;  0 dehydration/charge, 1  hydration/discharge
+CaO2H2Initial = 0.0
+PressureIn = 1.05e5                # [Pa] Inlet pressure; charge
+PressureOut = 1e5                  # [Pa] outlet pressure
+TemperatureIn = 573.15             # [K] inlet temperature: charge: 873 K ; discharge: 473K
+TemperatureOut = 573.15            # [K] outlet temperature: charge:  573; discharge: outflow
+InFlow = 5                         # 0.277[mol/s] Inflow of HTF
+VaporIn = 0.36                     # [] molefraction
+
+[Vtk]
+AddVelocity         = 0            # Add extra information
+VtuWritingFreq      = 1            # 1: write a vtu file at every timestep, 2: write a vtu file every second timestep ...
+
+[LinearSolver]
+ResidualReduction = 1e-6
+
+[Newton]
+MaxRelativeShift = 1e-6
+2WriteConvergence = 1
+
+[Output]
+#Frequency of restart file, flux and VTK output
+FreqRestart = 1000                 # how often restart files are written out
+FreqOutput = 50                    # frequency of VTK output
+FreqMassOutput = 2                 # frequency of mass and evaporation rate output (Darcy)
+FreqFluxOutput = 1000              # frequency of detailed flux output
+FreqVaporFluxOutput = 2            # frequency of summarized flux output
diff --git a/test/porousmediumflow/1pncmin/implicit/test_cc1pncmin.cc b/test/porousmediumflow/1pncmin/implicit/test_cc1pncmin.cc
new file mode 100644
index 0000000000..953a59ddf4
--- /dev/null
+++ b/test/porousmediumflow/1pncmin/implicit/test_cc1pncmin.cc
@@ -0,0 +1,55 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief Test for the 2pnc box model used for water management in PEM fuel cells.
+ */
+#include <config.h>
+#include "thermochemproblem.hh"
+#include <dumux/common/start.hh>
+
+/*!
+ * \brief Provides an interface for customizing error messages associated with
+ *        reading in parameters.
+ *
+ * \param progName  The name of the program, that was tried to be started.
+ * \param errorMsg  The error message that was issued by the start function.
+ *                  Comprises the thing that went wrong and a general help message.
+ */
+void usage(const char *progName, const std::string &errorMsg)
+{
+    if (errorMsg.size() > 0) {
+        std::string errorMessageOut = "\nUsage: ";
+                    errorMessageOut += progName;
+                    errorMessageOut += " [options]\n";
+                    errorMessageOut += errorMsg;
+                    errorMessageOut += "\n\nThe list of mandatory options for this program is:\n"
+                                        "\t-ParameterFile Parameter file (Input file) \n";
+
+        std::cout << errorMessageOut
+                  << "\n";
+    }
+}
+
+int main(int argc, char** argv)
+{
+    typedef TTAG(ThermoChemCCProblem) ProblemTypeTag;
+    return Dumux::start<ProblemTypeTag>(argc, argv, usage);
+}
diff --git a/test/porousmediumflow/1pncmin/implicit/test_cc1pncmin.input b/test/porousmediumflow/1pncmin/implicit/test_cc1pncmin.input
new file mode 100644
index 0000000000..c9e3985775
--- /dev/null
+++ b/test/porousmediumflow/1pncmin/implicit/test_cc1pncmin.input
@@ -0,0 +1,71 @@
+# Parameter file for test case 1pncmin.
+# Everything behind a '#' is a comment.
+# Type "./test_1pncmin --help" for more information.
+
+[TimeManager]
+DtInitial = 1                      # [s] initial time step size
+MaxTimeStepSize = 50               # [s] maximum time step size
+TEnd= 5000                        # [s] duration of the simulation
+FreqOutput          = 10           # frequency of VTK output
+WriteRestartFile    = 1            # Boolean. Should restart files be written? (1) Yes (0) No
+
+[Grid]
+UpperRight = 0.08 0.01 #20 10 #     # [m] upper right corner coordinates
+Cells =  40 2 # 21 6                # [-] number of cells in x,y-direction
+
+[FluidSystem]
+NTemperature = 10                   # [-] number of tabularization entries
+NPressure = 100                     # [-] number of tabularization entries
+PressureLow = 1E5                   # [Pa]low end for tabularization of fluid properties
+PressureHigh = 1E6                  # [Pa]high end for tabularization of fluid properties
+TemperatureLow = 373.15             # [Pa]low end for tabularization of fluid properties
+TemperatureHigh = 873.15            # [Pa]high end for tabularization of fluid properties
+
+[Problem]
+Name = cc1pncmin
+IsCharge = 1                       # Bool: 1: charge; 0: discharge
+
+[Charge]
+PressureInitial = 1E5              # [Pa] Initial reservoir pressure
+TemperatureInitial = 773.15        # [K]  reservoir temperature
+VaporInitial = 0.01                # [-]  initial mole fraction of water
+CaOInitial = 0.0                   # [-]  molefraction in the solid phase;  0 dehydration/charge, 1  hydration/discharge
+CaO2H2Initial = 0.3960             # [-]  molefraction in the solid phase;  0 dehydration/charge, 1  hydration/discharge
+PressureIn = 1.02e5                # [Pa] Inlet pressure; charge
+PressureOut = 1e5                  # [Pa] outlet pressure
+TemperatureIn = 773.15             # [K] inlet temperature: // Shao 500 °C
+TemperatureOut = 673.15            # [K] outlet temperature: // Shao noflow
+InFlow = 5                         # [mol/s] Inflow of HTF // Shao 0.309 g/s (Area (0.015)^2pi m^2) --> here 1.55 mol/s
+VaporIn = 0.01                     # [] molefraction // Shao 0.01 [g/g]
+
+[Discharge]
+PressureInitial = 1E5              # [Pa] Initial reservoir pressure
+TemperatureInitial = 573.15        # [K]  reservoir temperature
+VaporInitial = 0.35                # [-]  initial mole fraction of water
+CaOInitial = 0.2                   # [-]  molefraction in the solid phase;  0 dehydration/charge, 1  hydration/discharge
+CaO2H2Initial = 0.0
+PressureIn = 1.05e5                # [Pa] Inlet pressure; charge
+PressureOut = 1e5                  # [Pa] outlet pressure
+TemperatureIn = 573.15             # [K] inlet temperature: charge: 873 K ; discharge: 473K
+TemperatureOut = 573.15            # [K] outlet temperature: charge:  573; discharge: outflow
+InFlow = 5                         # 0.277[mol/s] Inflow of HTF
+VaporIn = 0.36                     # [] molefraction
+
+[Vtk]
+AddVelocity         = 0            # Add extra information
+VtuWritingFreq      = 1            # 1: write a vtu file at every timestep, 2: write a vtu file every second timestep ...
+
+[LinearSolver]
+ResidualReduction = 1e-6
+
+[Newton]
+MaxRelativeShift = 1e-6
+2WriteConvergence = 1
+
+[Output]
+#Frequency of restart file, flux and VTK output
+FreqRestart = 1000                 # how often restart files are written out
+FreqOutput = 50                    # frequency of VTK output
+FreqMassOutput = 2                 # frequency of mass and evaporation rate output (Darcy)
+FreqFluxOutput = 1000              # frequency of detailed flux output
+FreqVaporFluxOutput = 2            # frequency of summarized flux output
diff --git a/test/porousmediumflow/1pncmin/implicit/thermochemproblem.hh b/test/porousmediumflow/1pncmin/implicit/thermochemproblem.hh
index 366843be44..21847e375f 100644
--- a/test/porousmediumflow/1pncmin/implicit/thermochemproblem.hh
+++ b/test/porousmediumflow/1pncmin/implicit/thermochemproblem.hh
@@ -25,6 +25,9 @@
 #define DUMUX_THERMOCHEM_PROBLEM_HH
 
 #include <dumux/porousmediumflow/1pncmin/implicit/model.hh>
+#include <dumux/implicit/box/properties.hh>
+#include <dumux/implicit/cellcentered/tpfa/properties.hh>
+#include <dumux/implicit/cellcentered/mpfa/properties.hh>
 #include <dumux/porousmediumflow/implicit/problem.hh>
 #include <dumux/material/fluidsystems/simplesteamaircao2h2.hh>
 #include <dumux/material/fluidmatrixinteractions/1p/thermalconductivityaverage.hh>
@@ -45,11 +48,11 @@ namespace Properties
 #if NONISOTHERMAL
 NEW_TYPE_TAG(ThermoChemProblem, INHERITS_FROM(OnePNCMinNI, ThermoChemSpatialParams));
 NEW_TYPE_TAG(ThermoChemBoxProblem, INHERITS_FROM(BoxModel, ThermoChemProblem));
-NEW_TYPE_TAG(ThermoChemCCProblem, INHERITS_FROM(CCModel, ThermoChemProblem));
+NEW_TYPE_TAG(ThermoChemCCProblem, INHERITS_FROM(CCTpfaModel, ThermoChemProblem));
 #else
 NEW_TYPE_TAG(ThermoChemProblem, INHERITS_FROM(OnePNCMin, ThermoChemSpatialParams));
 NEW_TYPE_TAG(ThermoChemBoxProblem, INHERITS_FROM(BoxModel, ThermoChemProblem));
-NEW_TYPE_TAG(ThermoChemCCProblem, INHERITS_FROM(CCModel, ThermoChemProblem));
+NEW_TYPE_TAG(ThermoChemCCProblem, INHERITS_FROM(CCTpfaModel, ThermoChemProblem));
 #endif
 // Set the grid type
 SET_TYPE_PROP(ThermoChemProblem, Grid, Dune::YaspGrid<2>);
@@ -72,7 +75,7 @@ SET_TYPE_PROP(ThermoChemProblem, LinearSolver, UMFPackBackend<TypeTag>);
 
 
 // Set the spatial parameters
-SET_TYPE_PROP(ThermoChemProblem, SpatialParams, Dumux::ThermoChemSpatialParams<TypeTag>);
+SET_TYPE_PROP(ThermoChemProblem, SpatialParams, ThermoChemSpatialParams<TypeTag>);
 }
 
 /*!
@@ -137,13 +140,9 @@ class ThermoChemProblem : public ImplicitPorousMediaProblem<TypeTag>
     using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
     using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
 
-    using GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator);
-    using DimVector = Dune::FieldVector<Scalar, dim> ;
 
-//     typedef Dumux::Constants<Scalar> Constant;
-
-    enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
-    enum { dofCodim = isBox ? dim : 0 };
+//     enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+//     enum { dofCodim = isBox ? dim : 0 };
 
 public:
     /*!
diff --git a/test/porousmediumflow/1pncmin/implicit/thermochemspatialparams.hh b/test/porousmediumflow/1pncmin/implicit/thermochemspatialparams.hh
index b11543f55a..36ee0d5715 100644
--- a/test/porousmediumflow/1pncmin/implicit/thermochemspatialparams.hh
+++ b/test/porousmediumflow/1pncmin/implicit/thermochemspatialparams.hh
@@ -28,10 +28,6 @@
 
 #include <dumux/material/spatialparams/implicit1p.hh>
 #include <dumux/porousmediumflow/1pncmin/implicit/indices.hh>
-// #include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
-// #include <dumux/material/fluidmatrixinteractions/2p/regularizedbrookscorey.hh>
-// #include <dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh>
-// #include <dumux/material/fluidmatrixinteractions/2p/philtophoblaw.hh>
 #include <dumux/material/fluidmatrixinteractions/porosityreactivebed.hh>
 #include <dumux/material/fluidmatrixinteractions/permeabilitykozenycarman.hh>
 
@@ -63,7 +59,6 @@ class ThermoChemSpatialParams : public ImplicitSpatialParamsOneP<TypeTag>
 {
     using ThisType = ThermoChemSpatialParams<TypeTag>;
     using ParentType = ImplicitSpatialParamsOneP<TypeTag>;
-//     typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
     using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
     using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
@@ -80,11 +75,8 @@ class ThermoChemSpatialParams : public ImplicitSpatialParamsOneP<TypeTag>
     };
 
     using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
-//     typedef Dune::FieldVector<CoordScalar,dim> DimVector;
-//     typedef Dune::FieldMatrix<CoordScalar,dim,dim> DimMatrix;
     using Tensor = Dune::FieldMatrix<CoordScalar, dimWorld, dimWorld>;
-    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
-//     typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+//     using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
     using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
     using Element = typename GridView::template Codim<0>::Entity;
 
-- 
GitLab