Merge branch 'fix/improve-2p2cni-waterair-test' into 'master'

Fix/improve 2p2cni waterair test

Closes #851

See merge request !1949

test/porousmediumflow/2p2c/implicit/waterair/CMakeLists.txt

@@ -5,31 +5,31 @@ dumux_add_test(NAME test_2p2cni_waterair_box
               LABELS porousmediumflow 2p2c  2p2cni
               SOURCES main.cc
               TIMEOUT 1500
-              COMPILE_DEFINITIONS TYPETAG=WaterAirBox ENABLECACHING=0
+              COMPILE_DEFINITIONS TYPETAG=WaterAirBox
               COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
               CMD_ARGS --script fuzzy
                        --files ${CMAKE_SOURCE_DIR}/test/references/test_2p2cni_waterair_box-reference.vtu
-                               ${CMAKE_CURRENT_BINARY_DIR}/test_2p2cni_waterair_box-00013.vtu
+                               ${CMAKE_CURRENT_BINARY_DIR}/test_2p2cni_waterair_box-00012.vtu
                        --command "${CMAKE_CURRENT_BINARY_DIR}/test_2p2cni_waterair_box params.input -Problem.Name test_2p2cni_waterair_box")
 
 dumux_add_test(NAME test_2p2cni_waterair_buoyancy_box
               LABELS porousmediumflow 2p2c  2p2cni
               TARGET test_2p2cni_waterair_box
               TIMEOUT 1500
-              COMPILE_DEFINITIONS TYPETAG=WaterAirBox ENABLECACHING=0
+              COMPILE_DEFINITIONS TYPETAG=WaterAirBox
               COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
               CMD_ARGS --script fuzzy
                        --files ${CMAKE_SOURCE_DIR}/test/references/test_2p2cni_waterair_buoyancy_box-reference.vtu
-                               ${CMAKE_CURRENT_BINARY_DIR}/test_2p2cni_waterair_buoyancy_box-00015.vtu
+                               ${CMAKE_CURRENT_BINARY_DIR}/test_2p2cni_waterair_buoyancy_box-00013.vtu
                        --command "${CMAKE_CURRENT_BINARY_DIR}/test_2p2cni_waterair_box params.input -Problem.Name test_2p2cni_waterair_buoyancy_box -TimeLoop.DtInitial 100")
 
 dumux_add_test(NAME test_2p2cni_waterair_tpfa
               LABELS porousmediumflow 2p2c  2p2cni
               SOURCES main.cc
               TIMEOUT 1500
-              COMPILE_DEFINITIONS TYPETAG=WaterAirCCTpfa ENABLECACHING=0
+              COMPILE_DEFINITIONS TYPETAG=WaterAirCCTpfa
               COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
               CMD_ARGS --script fuzzy
                        --files ${CMAKE_SOURCE_DIR}/test/references/test_2p2cni_waterair_tpfa-reference.vtu
-                               ${CMAKE_CURRENT_BINARY_DIR}/test_2p2cni_waterair_tpfa-00013.vtu
+                               ${CMAKE_CURRENT_BINARY_DIR}/test_2p2cni_waterair_tpfa-00012.vtu
                        --command "${CMAKE_CURRENT_BINARY_DIR}/test_2p2cni_waterair_tpfa params.input -Problem.Name test_2p2cni_waterair_tpfa")

test/porousmediumflow/2p2c/implicit/waterair/main.cc

@@ -22,18 +22,14 @@
  * \brief Test for the two-phase two-component CC model.
  */
 #include <config.h>
-#include <ctime>
+
 #include <iostream>
 
 #include <dune/common/parallel/mpihelper.hh>
-#include <dune/common/timer.hh>
-#include <dune/grid/io/file/dgfparser/dgfexception.hh>
-#include <dune/grid/io/file/vtk.hh>
-#include <dune/istl/io.hh>
+#include <dune/common/exceptions.hh>
 
 #include <dumux/common/properties.hh>
 #include <dumux/common/parameters.hh>
-#include <dumux/common/valgrind.hh>
 #include <dumux/common/dumuxmessage.hh>
 
 #include <dumux/linear/amgbackend.hh>
@@ -43,11 +39,8 @@
 #include <dumux/assembly/fvassembler.hh>
 #include <dumux/assembly/diffmethod.hh>
 
-#include <dumux/discretization/method.hh>
-
 #include <dumux/io/vtkoutputmodule.hh>
-#include <dumux/io/grid/gridmanager.hh>
-#include <dumux/io/loadsolution.hh>
+#include <dumux/io/grid/gridmanager_yasp.hh>
 
 // the problem definitions
 #include "problem.hh"
@@ -73,10 +66,6 @@ int main(int argc, char** argv) try
     GridManager<GetPropType<TypeTag, Properties::Grid>> gridManager;
     gridManager.init();
 
-    ////////////////////////////////////////////////////////////
-    // run instationary non-linear problem on this grid
-    ////////////////////////////////////////////////////////////
-
     // we compute on the leaf grid view
     const auto& leafGridView = gridManager.grid().leafGridView();
 
@@ -85,34 +74,18 @@ int main(int argc, char** argv) try
     auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);
     gridGeometry->update();
 
+    ////////////////////////////////////////////////////////////
+    // run instationary non-linear problem on this grid
+    ////////////////////////////////////////////////////////////
+
     // the problem (initial and boundary conditions)
     using Problem = GetPropType<TypeTag, Properties::Problem>;
     auto problem = std::make_shared<Problem>(gridGeometry);
 
-    // get some time loop parameters
-    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
-    const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
-    const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
-    auto dt = getParam<Scalar>("TimeLoop.DtInitial");
-
-    // check if we are about to restart a previously interrupted simulation
-    Scalar restartTime = getParam<Scalar>("Restart.Time", 0);
-
     // the solution vector
     using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
-    SolutionVector x(gridGeometry->numDofs());
-    if (restartTime > 0)
-    {
-        using IOFields = GetPropType<TypeTag, Properties::IOFields>;
-        using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
-        using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
-        using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
-        const auto fileName = getParam<std::string>("Restart.File");
-        const auto pvName = createPVNameFunction<IOFields, PrimaryVariables, ModelTraits, FluidSystem>();
-        loadSolution(x, fileName, pvName, *gridGeometry);
-    }
-    else
-        problem->applyInitialSolution(x);
+    SolutionVector x;
+    problem->applyInitialSolution(x);
     auto xOld = x;
 
     // the grid variables
@@ -121,15 +94,19 @@ int main(int argc, char** argv) try
     gridVariables->init(x);
 
     // intialize the vtk output module
-    using IOFields = GetPropType<TypeTag, Properties::IOFields>;
     VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
     using VelocityOutput = GetPropType<TypeTag, Properties::VelocityOutput>;
     vtkWriter.addVelocityOutput(std::make_shared<VelocityOutput>(*gridVariables));
-    IOFields::initOutputModule(vtkWriter); // Add model specific output fields
-    vtkWriter.write(restartTime);
+    // Add model specific output fields
+    GetPropType<TypeTag, Properties::IOFields>::initOutputModule(vtkWriter);
+    vtkWriter.write(0.0);
 
     // instantiate time loop
-    auto timeLoop = std::make_shared<TimeLoop<Scalar>>(restartTime, dt, tEnd);
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+    const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+    const auto dtInitial = getParam<Scalar>("TimeLoop.DtInitial");
+    auto timeLoop = std::make_shared<TimeLoop<Scalar>>(0.0, dtInitial, tEnd);
     timeLoop->setMaxTimeStepSize(maxDt);
 
     // the assembler with time loop for instationary problem
@@ -160,8 +137,8 @@ int main(int argc, char** argv) try
         // report statistics of this time step
         timeLoop->reportTimeStep();
 
-        // set new dt as suggested by the newton solver
-        timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
+        // set new dt as twice the old time step size
+        timeLoop->setTimeStepSize(2.0*timeLoop->timeStepSize());
 
         // write vtk output
         vtkWriter.write(timeLoop->time());
@@ -183,12 +160,12 @@ int main(int argc, char** argv) try
 
     return 0;
 } // end main
-catch (Dumux::ParameterException &e)
+catch (const Dumux::ParameterException &e)
 {
     std::cerr << std::endl << e << " ---> Abort!" << std::endl;
     return 1;
 }
-catch (Dune::DGFException & e)
+catch (const Dune::DGFException & e)
 {
     std::cerr << "DGF exception thrown (" << e <<
                  "). Most likely, the DGF file name is wrong "
@@ -197,7 +174,7 @@ catch (Dune::DGFException & e)
                  << " ---> Abort!" << std::endl;
     return 2;
 }
-catch (Dune::Exception &e)
+catch (const Dune::Exception &e)
 {
     std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
     return 3;

test/porousmediumflow/2p2c/implicit/waterair/problem.hh

@@ -30,7 +30,6 @@
 
 #include <dumux/discretization/cctpfa.hh>
 #include <dumux/discretization/box.hh>
-#include <dumux/discretization/method.hh>
 
 #include <dumux/material/components/n2.hh>
 #include <dumux/material/fluidsystems/h2on2.hh>
@@ -76,6 +75,14 @@ struct SpatialParams<TypeTag, TTag::WaterAir>
 // Define whether mole(true) or mass (false) fractions are used
 template<class TypeTag>
 struct UseMoles<TypeTag, TTag::WaterAir> { static constexpr bool value = true; };
+
+// Enable caching
+template<class TypeTag>
+struct EnableGridGeometryCache<TypeTag, TTag::WaterAir> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridVolumeVariablesCache<TypeTag, TTag::WaterAir> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridFluxVariablesCache<TypeTag, TTag::WaterAir> { static constexpr bool value = true; };
 } // end namespace Dumux
 
 /*!
@@ -118,7 +125,8 @@ class WaterAirProblem : public PorousMediumFlowProblem<TypeTag>
     using ParentType = PorousMediumFlowProblem<TypeTag>;
 
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
-    using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
+    using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+    using GridView = typename GridGeometry::GridView;
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
     using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>;
     using Indices = typename ModelTraits::Indices;
@@ -149,10 +157,6 @@ class WaterAirProblem : public PorousMediumFlowProblem<TypeTag>
     using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
     using Element = typename GridView::template Codim<0>::Entity;
     using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
-    using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
-    using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
-    using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
-    using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
 
     //! Property that defines whether mole or mass fractions are used
     static constexpr bool useMoles = ModelTraits::useMoles();

test/porousmediumflow/2p2c/implicit/waterair/spatialparams.hh

@@ -44,27 +44,23 @@ class WaterAirSpatialParams
 : public FVSpatialParams<GridGeometry, Scalar,
                          WaterAirSpatialParams<GridGeometry, Scalar>>
 {
+    using ThisType = WaterAirSpatialParams<GridGeometry, Scalar>;
+    using ParentType = FVSpatialParams<GridGeometry, Scalar, ThisType>;
     using GridView = typename GridGeometry::GridView;
     using FVElementGeometry = typename GridGeometry::LocalView;
-    using SubControlVolume = typename FVElementGeometry::SubControlVolume;
     using Element = typename GridView::template Codim<0>::Entity;
-    using ParentType = FVSpatialParams<GridGeometry, Scalar,
-                                       WaterAirSpatialParams<GridGeometry, Scalar>>;
-
     static constexpr int dimWorld = GridView::dimensionworld;
-
-    using EffectiveLaw = RegularizedBrooksCorey<Scalar>;
-
     using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
 
 public:
     //! Export the type used for the permeability
     using PermeabilityType = Scalar;
     //! Export the type used for the material law
-    using MaterialLaw = EffToAbsLaw<EffectiveLaw>;
+    using MaterialLaw = EffToAbsLaw<RegularizedBrooksCorey<Scalar>>;
     using MaterialLawParams = typename MaterialLaw::Params;
 
-    WaterAirSpatialParams(std::shared_ptr<const GridGeometry> gridGeometry) : ParentType(gridGeometry)
+    WaterAirSpatialParams(std::shared_ptr<const GridGeometry> gridGeometry)
+    : ParentType(gridGeometry)
     {
         layerBottom_ = 22.0;