From 087c8e14ac983668a786ed3a786c4714ad7eed7e Mon Sep 17 00:00:00 2001
From: "Dennis.Glaeser" <dennis.glaeser@iws.uni-stuttgart.de>
Date: Thu, 11 Jul 2019 17:32:04 +0200
Subject: [PATCH 1/3] [test][1p] add convergence test

---
 .../1p/implicit/CMakeLists.txt                |   1 +
 .../1p/implicit/convergence/CMakeLists.txt    |  41 ++++
 .../1p/implicit/convergence/main.cc           | 215 ++++++++++++++++++
 .../1p/implicit/convergence/params.input      |  23 ++
 .../1p/implicit/convergence/problem.hh        | 211 +++++++++++++++++
 .../1p/implicit/convergence/solver.hh         | 149 ++++++++++++
 .../1p/implicit/convergence/spatialparams.hh  |  75 ++++++
 7 files changed, 715 insertions(+)
 create mode 100644 test/porousmediumflow/1p/implicit/convergence/CMakeLists.txt
 create mode 100644 test/porousmediumflow/1p/implicit/convergence/main.cc
 create mode 100644 test/porousmediumflow/1p/implicit/convergence/params.input
 create mode 100644 test/porousmediumflow/1p/implicit/convergence/problem.hh
 create mode 100644 test/porousmediumflow/1p/implicit/convergence/solver.hh
 create mode 100644 test/porousmediumflow/1p/implicit/convergence/spatialparams.hh

diff --git a/test/porousmediumflow/1p/implicit/CMakeLists.txt b/test/porousmediumflow/1p/implicit/CMakeLists.txt
index dbbc593c5a..1e457762fa 100644
--- a/test/porousmediumflow/1p/implicit/CMakeLists.txt
+++ b/test/porousmediumflow/1p/implicit/CMakeLists.txt
@@ -1,3 +1,4 @@
+add_subdirectory(convergence)
 add_subdirectory(pointsources)
 add_subdirectory(incompressible)
 add_subdirectory(compressible)
diff --git a/test/porousmediumflow/1p/implicit/convergence/CMakeLists.txt b/test/porousmediumflow/1p/implicit/convergence/CMakeLists.txt
new file mode 100644
index 0000000000..620355cb05
--- /dev/null
+++ b/test/porousmediumflow/1p/implicit/convergence/CMakeLists.txt
@@ -0,0 +1,41 @@
+dune_symlink_to_source_files(FILES "params.input" "ploterrors.py")
+
+# executable for tpfa tests
+add_executable(test_1p_convergence_tpfa EXCLUDE_FROM_ALL main.cc)
+target_compile_definitions(test_1p_convergence_tpfa PUBLIC "TYPETAG=OnePIncompressibleTpfa")
+
+# executable for box tests
+add_executable(test_1p_convergence_box EXCLUDE_FROM_ALL main.cc)
+target_compile_definitions(test_1p_convergence_box PUBLIC "TYPETAG=OnePIncompressibleBox")
+
+# using tpfa and conforming refinement
+dune_add_test(NAME test_1p_convergence_tpfa_conforming
+              TARGET test_1p_convergence_tpfa
+              CMAKE_GUARD "( dune-functions_FOUND )"
+              LABELS porousmediumflow 1p
+              COMMAND ./test_1p_convergence_tpfa
+              CMD_ARGS params.input -Problem.Name test_1p_convergence_tpfa_conforming)
+
+# using tpfa and nonconforming refinement
+dune_add_test(NAME test_1p_convergence_tpfa_nonconforming
+              TARGET test_1p_convergence_tpfa
+              CMAKE_GUARD "( dune-functions_FOUND )"
+              LABELS porousmediumflow 1p
+              COMMAND ./test_1p_convergence_tpfa
+              CMD_ARGS params.input -Problem.Name test_1p_convergence_tpfa_conforming)
+
+# using box and conforming refinement
+dune_add_test(NAME test_1p_convergence_box_conforming
+              TARGET test_1p_convergence_box
+              CMAKE_GUARD "( dune-functions_FOUND )"
+              LABELS porousmediumflow 1p
+              COMMAND ./test_1p_convergence_box
+              CMD_ARGS params.input -Problem.Name test_1p_convergence_box_conforming)
+
+# using box and nonconforming refinement
+dune_add_test(NAME test_1p_convergence_box_nonconforming
+              TARGET test_1p_convergence_box
+              CMAKE_GUARD "( dune-functions_FOUND )"
+              LABELS porousmediumflow 1p
+              COMMAND ./test_1p_convergence_box
+              CMD_ARGS params.input -Problem.Name test_1p_convergence_box_conforming)
diff --git a/test/porousmediumflow/1p/implicit/convergence/main.cc b/test/porousmediumflow/1p/implicit/convergence/main.cc
new file mode 100644
index 0000000000..a0b0890b8e
--- /dev/null
+++ b/test/porousmediumflow/1p/implicit/convergence/main.cc
@@ -0,0 +1,215 @@
+// -*- 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 Convergence test for single-phase flow
+ */
+#include <config.h>
+
+#include <iostream>
+#include <fstream>
+#include <utility>
+#include <algorithm>
+#include <type_traits>
+
+#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+#include <dune/functions/gridfunctions/analyticgridviewfunction.hh>
+#include <dune/common/parallel/mpihelper.hh>
+
+#include <dumux/common/parameters.hh>
+#include <dumux/common/properties.hh>
+#include <dumux/common/dumuxmessage.hh>
+
+#include "solver.hh"
+
+#include <dumux/common/integrate.hh>
+#include <dumux/multidomain/glue.hh>
+
+#include <dumux/discretization/functionspacebasis.hh>
+#include <dumux/discretization/projection/projector.hh>
+
+// main program
+int main(int argc, char** argv) try
+{
+    using namespace Dumux;
+    using TypeTag = Properties::TTag::TYPETAG;
+    static constexpr auto dm = GetPropType<TypeTag, Properties::FVGridGeometry>::discMethod;
+    static constexpr bool isBox = dm == DiscretizationMethod::box;
+
+    // initialize MPI, finalize is done automatically on exit
+    const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+    // print dumux start message
+    if (mpiHelper.rank() == 0)
+        DumuxMessage::print(/*firstCall=*/true);
+
+    // initialize parameter tree
+    Parameters::init(argc, argv);
+
+    // obtain the sequence of discretizations
+    auto cellsVector = getParam< std::vector<int> >("Grid.CellsSequence");
+
+    std::sort(cellsVector.begin(), cellsVector.end());
+    cellsVector.erase(std::unique(cellsVector.begin(), cellsVector.end()), cellsVector.end());
+
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    std::vector<Scalar> dxVector(cellsVector.size());
+    std::transform(cellsVector.begin(),
+                   cellsVector.end(),
+                   dxVector.begin(),
+                   [] (auto cells) { return 1.0/cells; });
+
+    // run the finest-resolved simulation first
+    std::cout << "\n --- Solving finest solution (dx = " << 1.0/cellsVector.back() << ") --- \n\n";
+    const auto finestStorage = solveRefinementLevel<TypeTag>(cellsVector.back());
+    const auto& finestGridGeometry = *finestStorage.gridGeometry;
+    const auto& finestSolution = *finestStorage.solution;
+    const auto& finestBasis = getFunctionSpaceBasis(finestGridGeometry);
+
+    // grid function with the discrete solution
+    using namespace Dune::Functions;
+    using BlockType = typename std::decay_t<decltype(finestSolution)>::block_type;
+    const auto gfFine = makeDiscreteGlobalBasisFunction<BlockType>(finestBasis, finestSolution);
+
+    // grid function with the exact solution
+    using Problem = GetPropType<TypeTag, Properties::Problem>;
+    auto analyticSol = [] (const auto x) { return Problem::exact(x); };
+    const auto gfAnalyticFine = makeAnalyticGridViewFunction(analyticSol, finestBasis.gridView());
+
+    // run simulations with coarser grids and compute errors
+    std::vector<Scalar> discreteErrors; discreteErrors.reserve(cellsVector.size()-1);
+    std::vector<Scalar> analyticErrors; analyticErrors.reserve(cellsVector.size());
+
+    cellsVector.pop_back();
+    const auto intOrder = getParam<int>("L2Norm.IntegrationOrder");
+    for (auto dx : cellsVector)
+    {
+        std::cout << "\n --- Solving with dx = " << dx << " --- \n\n";
+        const auto storage = solveRefinementLevel<TypeTag>(dx);
+        const auto& gridGeometry = *storage.gridGeometry;
+        const auto& solution = *storage.solution;
+        const auto& basis = getFunctionSpaceBasis(gridGeometry);
+
+        const auto gf = makeDiscreteGlobalBasisFunction<BlockType>(basis, solution);
+        const auto gfAnalytic = makeAnalyticGridViewFunction(analyticSol, basis.gridView());
+
+        std::cout << "Projecting solution into reference space" << std::endl;
+        const auto glue = makeGlue(gridGeometry, finestGridGeometry);
+        const auto projector = makeProjector(basis, finestBasis, glue);
+
+        auto params = projector.defaultParams();
+        params.residualReduction = 1e-16;
+
+        const auto projectedSolution = projector.project(solution, params);
+        const auto gfProjected = makeDiscreteGlobalBasisFunction<BlockType>(finestBasis, projectedSolution);
+
+        std::cout << "Computing error norms" << std::endl;
+        discreteErrors.push_back( integrateL2Error(finestBasis.gridView(), gfFine, gfProjected, intOrder) );
+        analyticErrors.push_back( integrateL2Error(basis.gridView(), gf, gfAnalytic, intOrder) );
+    }
+
+    // compute analytical error for finest solution
+    std::cout << "Computing error norm for finest grid w.r.t. analytical solution" << std::endl;
+    analyticErrors.push_back( integrateL2Error(finestBasis.gridView(), gfFine, gfAnalyticFine, intOrder) );
+
+    using std::log;
+    std::cout << "\nComputed the following errors/rates w.r.t. the analytical solution:\n";
+    for (unsigned int i = 0; i < analyticErrors.size(); ++i)
+    {
+        std::cout << analyticErrors[i];
+        if (i == 0)
+            std::cout << "\n";
+        else
+        {
+            const auto rate = (log(analyticErrors[i]) - log(analyticErrors[i-1]))
+                              /(log(dxVector[i]) - log(dxVector[i-1]));
+
+            std::cout << " \t--->\t " << rate << std::endl;
+            if (i == analyticErrors.size()-1)
+                if ( (isBox && rate < 1.95) || (!isBox && rate < 0.95) )
+                    DUNE_THROW(Dune::InvalidStateException, "Computed rate is below expected value: " << rate);
+        }
+    }
+
+    std::cout << "\nComputed the following errors/rates w.r.t. the discrete solution:\n";
+    for (unsigned int i = 0; i < discreteErrors.size(); ++i)
+    {
+        std::cout << discreteErrors[i];
+        if (i == 0)
+            std::cout << "\n";
+        else
+        {
+            const auto rate = (log(discreteErrors[i]) - log(discreteErrors[i-1]))
+                              /(log(dxVector[i]) - log(dxVector[i-1]));
+
+            std::cout << " \t--->\t " << rate << std::endl;
+            if (i == discreteErrors.size()-1)
+                if ( (isBox && rate < 1.95) || (!isBox && rate < 0.95) )
+                    DUNE_THROW(Dune::InvalidStateException, "Computed rate is below expected value: " << rate);
+        }
+    }
+
+    // maybe write the errors to disk
+    if (getParam<bool>("L2Norm.WriteLogFiles", false))
+    {
+        std::ofstream analyticFile(getParam<std::string>("L2Norm.AnalyticErrorsFile"));
+        std::ofstream discreteFile(getParam<std::string>("L2Norm.DiscreteErrorsFile"));
+
+        for (unsigned int i = 0; i < analyticErrors.size(); ++i)
+            analyticFile << dxVector[i] << "," << analyticErrors[i] << "\n";
+        for (unsigned int i = 0; i < discreteErrors.size(); ++i)
+            discreteFile << dxVector[i] << "," << discreteErrors[i] << "\n";
+
+        // show parameters used
+        Parameters::print();
+
+        // print dumux good bye message
+        if (mpiHelper.rank() == 0)
+            DumuxMessage::print(/*firstCall=*/false);
+
+        analyticFile.close();
+        discreteFile.close();
+    }
+
+    return 0;
+}
+catch (Dumux::ParameterException &e)
+{
+    std::cerr << std::endl << e << " ---> Abort!" << std::endl;
+    return 1;
+}
+catch (Dune::DGFException & e)
+{
+    std::cerr << "DGF exception thrown (" << e <<
+                 "). Most likely, the DGF file name is wrong "
+                 "or the DGF file is corrupted, "
+                 "e.g. missing hash at end of file or wrong number (dimensions) of entries."
+                 << " ---> Abort!" << std::endl;
+    return 2;
+}
+catch (Dune::Exception &e)
+{
+    std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
+    return 3;
+}
+catch (...)
+{
+    std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
+    return 4;
+}
diff --git a/test/porousmediumflow/1p/implicit/convergence/params.input b/test/porousmediumflow/1p/implicit/convergence/params.input
new file mode 100644
index 0000000000..d4f8a26840
--- /dev/null
+++ b/test/porousmediumflow/1p/implicit/convergence/params.input
@@ -0,0 +1,23 @@
+[Problem]
+Name = 1p_convergence
+SourceIntegrationOrder = 3
+ExactSolPeriodLength = 6.283 # ~ 2*pi
+EnableGravity = false
+
+[Grid]
+LowerLeft = 0 0
+UpperRight = 1 1
+CellsSequence = 10 20 40 80
+
+[IO]
+WriteVTK = true
+
+[L2Norm]
+IntegrationOrder = 3
+WriteLogFiles = false
+AnalyticErrorsFile = conforming_analytic.log
+DiscreteErrorsFile = conforming_discrete.log
+
+[Component]
+LiquidDensity = 1
+LiquidKinematicViscosity = 1
diff --git a/test/porousmediumflow/1p/implicit/convergence/problem.hh b/test/porousmediumflow/1p/implicit/convergence/problem.hh
new file mode 100644
index 0000000000..f0bfd880a9
--- /dev/null
+++ b/test/porousmediumflow/1p/implicit/convergence/problem.hh
@@ -0,0 +1,211 @@
+// -*- 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 3 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
+ * \ingroup OnePTests
+ * \brief The properties & problem setup for the convergence test
+ */
+#ifndef DUMUX_INCOMPRESSIBLE_ONEP_CONVERGENCETEST_PROBLEM_HH
+#define DUMUX_INCOMPRESSIBLE_ONEP_CONVERGENCETEST_PROBLEM_HH
+
+#include <cmath>
+#include <dune/grid/yaspgrid.hh>
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dumux/discretization/cctpfa.hh>
+#include <dumux/discretization/ccmpfa.hh>
+#include <dumux/discretization/box.hh>
+
+#include <dumux/porousmediumflow/problem.hh>
+#include <dumux/porousmediumflow/1p/model.hh>
+#include <dumux/porousmediumflow/1p/incompressiblelocalresidual.hh>
+
+#include <dumux/material/components/constant.hh>
+#include <dumux/material/fluidsystems/1pliquid.hh>
+
+#include "spatialparams.hh"
+
+namespace Dumux {
+// forward declarations
+template<class TypeTag> class OnePTestProblem;
+
+namespace Properties {
+
+// create the type tag nodes
+namespace TTag {
+struct OnePIncompressible { using InheritsFrom = std::tuple<OneP>; };
+struct OnePIncompressibleTpfa { using InheritsFrom = std::tuple<OnePIncompressible, CCTpfaModel>; };
+struct OnePIncompressibleMpfa { using InheritsFrom = std::tuple<OnePIncompressible, CCMpfaModel>; };
+struct OnePIncompressibleBox { using InheritsFrom = std::tuple<OnePIncompressible, BoxModel>; };
+} // end namespace TTag
+
+// Set the grid type
+template<class TypeTag>
+struct Grid<TypeTag, TTag::OnePIncompressible> { using type = Dune::YaspGrid<2>; };
+
+// Set the problem type
+template<class TypeTag>
+struct Problem<TypeTag, TTag::OnePIncompressible> { using type = OnePTestProblem<TypeTag>; };
+
+// set the spatial params
+template<class TypeTag>
+struct SpatialParams<TypeTag, TTag::OnePIncompressible>
+{
+    using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using type = OnePTestSpatialParams<FVGridGeometry, Scalar>;
+};
+
+// use the incompressible local residual (provides analytic jacobian)
+template<class TypeTag>
+struct LocalResidual<TypeTag, TTag::OnePIncompressible> { using type = OnePIncompressibleLocalResidual<TypeTag>; };
+
+// the fluid system
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::OnePIncompressible>
+{
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using type = FluidSystems::OnePLiquid<Scalar, Components::Constant<0, Scalar> >;
+};
+
+} // end namespace Properties
+
+/*!
+ * \ingroup OnePTests
+ * \brief problem setup for the convergence test
+ */
+template<class TypeTag>
+class OnePTestProblem : public PorousMediumFlowProblem<TypeTag>
+{
+    using ParentType = PorousMediumFlowProblem<TypeTag>;
+
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
+    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+    using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
+
+    using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
+    using FVElementGeometry = typename FVGridGeometry::LocalView;
+    using SubControlVolume = typename FVGridGeometry::SubControlVolume;
+    using GridView = typename FVGridGeometry::GridView;
+    using Element = typename GridView::template Codim<0>::Entity;
+    using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+public:
+    /*!
+     * \brief The constructor.
+     * \param fvGridGeometry The finite-volume grid geometry
+     */
+    OnePTestProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+    : ParentType(fvGridGeometry)
+    {}
+
+    /*!
+     * \brief Specifies which kind of boundary condition should be
+     *        used for which equation on a given boundary control volume.
+     * \param globalPos The position of the center of the finite volume
+     */
+    BoundaryTypes boundaryTypesAtPos(const GlobalPosition& globalPos) const
+    {
+        BoundaryTypes values;
+        values.setAllDirichlet();
+        return values;
+    }
+
+    /*!
+     * \brief Evaluates the boundary conditions for a Dirichlet control volume.
+     * \param globalPos The center of the finite volume for which it is to be set.
+     */
+    PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
+    { return exact(globalPos); }
+
+    /*!
+     * \brief Evaluates the source term within a sub-control volume.
+     * \param element The finite element
+     * \param fvGeometry The element finite-volume geometry
+     * \param elemVolVars The element volume variables
+     * \param scv The sub-control volume for which the source term is evaluated
+     */
+    template<class ElementVolumeVariables>
+    NumEqVector source(const Element& element,
+                       const FVElementGeometry& fvGeometry,
+                       const ElementVolumeVariables& elemVolVars,
+                       const SubControlVolume& scv) const
+    {
+        static const auto order = getParam<Scalar>("Problem.SourceIntegrationOrder");
+        static const auto periodLength = getParam<Scalar>("Problem.ExactSolPeriodLength");
+        const auto& k = this->spatialParams().permeabilityAtPos(scv.center());
+
+        using std::sin;
+        using std::cos;
+
+        const auto eg = element.geometry();
+        const auto rule = Dune::QuadratureRules<Scalar, GridView::dimension>::rule(eg.type(), order);
+
+        Scalar source = 0.0;
+        for (auto qp : rule)
+        {
+            const auto p = eg.global(qp.position());
+            const auto x = p[0];
+            const auto y = p[1];
+
+            const auto preFactor = -1.0*periodLength*periodLength*M_PI*M_PI;
+            const auto preFactorArg = periodLength*M_PI;
+            const auto sineTerm = sin(preFactorArg*x);
+            const auto cosTerm = cos(preFactorArg*y);
+            const auto secondDeriv = preFactor*sineTerm*cosTerm;
+
+            // derivative in x and y are identical
+            source -= 2.0*k*secondDeriv*qp.weight()*eg.integrationElement(qp.position());
+        }
+
+        source /= eg.volume();
+        return NumEqVector(source);
+    }
+
+    /*!
+     * \brief Returns the temperature \f$\mathrm{[K]}\f$ for an isothermal problem.
+     * \param globalPos The center of the finite volume for which it is to be set.
+     */
+    Scalar temperature() const
+    { return 283.15; }
+
+    /*!
+     * \brief Returns the exact solution at a position.
+     * \param globalPos The center of the finite volume for which it is to be set.
+     */
+    static PrimaryVariables exact(const GlobalPosition& globalPos)
+    {
+        const auto x = globalPos[0];
+        const auto y = globalPos[1];
+
+        using std::sin;
+        using std::cos;
+
+        static const auto periodLength = getParam<Scalar>("Problem.ExactSolPeriodLength");
+        const auto preFactorArg = periodLength*M_PI;
+        const auto u = sin(preFactorArg*x)*cos(preFactorArg*y);
+
+        return PrimaryVariables(u);
+    }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/test/porousmediumflow/1p/implicit/convergence/solver.hh b/test/porousmediumflow/1p/implicit/convergence/solver.hh
new file mode 100644
index 0000000000..93eef93e25
--- /dev/null
+++ b/test/porousmediumflow/1p/implicit/convergence/solver.hh
@@ -0,0 +1,149 @@
+// -*- 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 3 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
+ * \ingroup OnePTests
+ * \brief The solver of the single-phase convergence test
+ */
+#ifndef DUMUX_INCOMPRESSIBLE_ONEP_CONVERGENCETEST_SOLVER_HH
+#define DUMUX_INCOMPRESSIBLE_ONEP_CONVERGENCETEST_SOLVER_HH
+
+#include <iostream>
+#include <string>
+
+#include <dune/common/timer.hh>
+#include <dune/grid/io/file/dgfparser/dgfexception.hh>
+#include <dune/grid/io/file/vtk.hh>
+
+#include <dumux/linear/seqsolverbackend.hh>
+#include <dumux/linear/pdesolver.hh>
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/grid/gridmanager.hh>
+
+#include <dumux/assembly/fvassembler.hh>
+#include <dumux/assembly/diffmethod.hh>
+
+#include "problem.hh"
+
+// return type of solveRefinementLevel()
+// stores the grid geometry and the produced solution
+template<class TypeTag>
+struct SolutionStorage
+{
+    using Grid = Dumux::GetPropType<TypeTag, Dumux::Properties::Grid>;
+    using GridGeometry = Dumux::GetPropType<TypeTag, Dumux::Properties::FVGridGeometry>;
+    using SolutionVector = Dumux::GetPropType<TypeTag, Dumux::Properties::SolutionVector>;
+
+public:
+    Dumux::GridManager<Grid> gridManager;
+    std::shared_ptr<GridGeometry> gridGeometry;
+    std::shared_ptr<SolutionVector> solution;
+};
+
+/*!
+ * \brief Solves the problem for a given number of cells per direction.
+ * \param numCells the number of cells per direction to be used
+ * \return returns an object of SolutionStorage, which carries
+ *         the grid and the solution after solving the problem
+ */
+template<class TypeTag>
+SolutionStorage<TypeTag> solveRefinementLevel(int numCells)
+{
+    using namespace Dumux;
+
+    // adapt the parameter tree to carry given number of cells
+    const auto c = std::to_string(numCells);
+    Parameters::init([&c] (auto& tree) { tree["Grid.Cells"] = c + " " + c; });
+
+    // the returned object of this function
+    // we create the grid in there directly
+    SolutionStorage<TypeTag> storage;
+    auto& gridManager = storage.gridManager;
+    gridManager.init();
+
+    // we compute on the leaf grid view
+    const auto& leafGridView = gridManager.grid().leafGridView();
+
+    // start timer
+    Dune::Timer timer;
+
+    // create the finite volume grid geometry
+    using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
+    auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
+    fvGridGeometry->update();
+
+    // the problem (boundary conditions)
+    using Problem = GetPropType<TypeTag, Properties::Problem>;
+    auto problem = std::make_shared<Problem>(fvGridGeometry);
+
+    // the solution vector
+    using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
+    auto x = std::make_shared<SolutionVector>(fvGridGeometry->numDofs());
+
+    // the grid variables
+    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
+    auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
+    gridVariables->init(*x);
+
+    // create assembler & linear solver
+    using Assembler = FVAssembler<TypeTag, DiffMethod::analytic>;
+    auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables);
+
+    using LinearSolver = ILU0BiCGSTABBackend;
+    auto linearSolver = std::make_shared<LinearSolver>();
+
+    // solver the linear problem
+    LinearPDESolver<Assembler, LinearSolver> solver(assembler,  linearSolver);
+    solver.solve(*x);
+
+    // maybe output result to vtk
+    if (getParam<bool>("IO.WriteVTK", false))
+    {
+        VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, *x, problem->name());
+        using IOFields = GetPropType<TypeTag, Properties::IOFields>;
+        IOFields::initOutputModule(vtkWriter); // Add model specific output fields
+
+        // add exact solution
+        using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+        std::vector<Scalar> exact(fvGridGeometry->numDofs());
+
+        for (const auto& element : elements(fvGridGeometry->gridView()))
+        {
+            auto fvGeometry = localView(*fvGridGeometry);
+            fvGeometry.bindElement(element);
+
+            for (const auto& scv : scvs(fvGeometry))
+                exact[scv.dofIndex()] = problem->exact(scv.dofPosition());
+        }
+
+        vtkWriter.addField(exact, "p_exact");
+        vtkWriter.write(1.0);
+    }
+
+    // fill storage and return
+    storage.gridGeometry = fvGridGeometry;
+    storage.solution = x;
+    return storage;
+}
+
+#endif
diff --git a/test/porousmediumflow/1p/implicit/convergence/spatialparams.hh b/test/porousmediumflow/1p/implicit/convergence/spatialparams.hh
new file mode 100644
index 0000000000..9c245763e3
--- /dev/null
+++ b/test/porousmediumflow/1p/implicit/convergence/spatialparams.hh
@@ -0,0 +1,75 @@
+// -*- 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 3 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
+ * \ingroup OnePTests
+ * \brief The spatial params of the incompressible single-phase convergence test
+ */
+#ifndef DUMUX_INCOMPRESSIBLE_ONEP_CONVERGENCETEST_SPATIAL_PARAMS_HH
+#define DUMUX_INCOMPRESSIBLE_ONEP_CONVERGENCETEST_SPATIAL_PARAMS_HH
+
+#include <dumux/porousmediumflow/properties.hh>
+#include <dumux/material/spatialparams/fv1p.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup OnePTests
+ * \brief The spatial params of the incompressible single-phase convergence test
+ */
+template<class FVGridGeometry, class Scalar>
+class OnePTestSpatialParams
+: public FVSpatialParamsOneP<FVGridGeometry, Scalar,
+                             OnePTestSpatialParams<FVGridGeometry, Scalar>>
+{
+    using ThisType = OnePTestSpatialParams<FVGridGeometry, Scalar>;
+    using ParentType = FVSpatialParamsOneP<FVGridGeometry, Scalar, ThisType>;
+
+    using Element = typename FVGridGeometry::GridView::template Codim<0>::Entity;
+    using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+public:
+    using PermeabilityType = Scalar;
+
+    /*!
+     * \brief The constructor.
+     * \param fvGridGeometry The finite-volume grid geometry
+     */
+    OnePTestSpatialParams(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+    : ParentType(fvGridGeometry)
+    {}
+
+    /*!
+     * \brief Function for defining the (intrinsic) permeability \f$[m^2]\f$.
+     * \param globalPos The global position
+     */
+    PermeabilityType permeabilityAtPos(const GlobalPosition& globalPos) const
+    { return 1.0; }
+
+    /*!
+     * \brief Defines the porosity \f$\mathrm{[-]}\f$.
+     * \param globalPos The global position
+     */
+    Scalar porosityAtPos(const GlobalPosition& globalPos) const
+    { return 1.0; }
+};
+
+} // end namespace Dumux
+
+#endif
-- 
GitLab


From aa09e22458aff0773890c33dcdaa305a14cc640a Mon Sep 17 00:00:00 2001
From: "Dennis.Glaeser" <dennis.glaeser@iws.uni-stuttgart.de>
Date: Mon, 15 Jul 2019 08:45:30 +0200
Subject: [PATCH 2/3] [test][1p][convergence] always print parameters and
 goodbye message

---
 .../1p/implicit/convergence/main.cc                | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/test/porousmediumflow/1p/implicit/convergence/main.cc b/test/porousmediumflow/1p/implicit/convergence/main.cc
index a0b0890b8e..3f93cbca75 100644
--- a/test/porousmediumflow/1p/implicit/convergence/main.cc
+++ b/test/porousmediumflow/1p/implicit/convergence/main.cc
@@ -176,17 +176,17 @@ int main(int argc, char** argv) try
         for (unsigned int i = 0; i < discreteErrors.size(); ++i)
             discreteFile << dxVector[i] << "," << discreteErrors[i] << "\n";
 
-        // show parameters used
-        Parameters::print();
-
-        // print dumux good bye message
-        if (mpiHelper.rank() == 0)
-            DumuxMessage::print(/*firstCall=*/false);
-
         analyticFile.close();
         discreteFile.close();
     }
 
+    // show parameters used
+    Parameters::print();
+
+    // print dumux good bye message
+    if (mpiHelper.rank() == 0)
+        DumuxMessage::print(/*firstCall=*/false);
+
     return 0;
 }
 catch (Dumux::ParameterException &e)
-- 
GitLab


From 0d70cacc2c886b95cc0cb60bdec59d5b05994dac Mon Sep 17 00:00:00 2001
From: "Dennis.Glaeser" <dennis.glaeser@iws.uni-stuttgart.de>
Date: Mon, 15 Jul 2019 08:47:40 +0200
Subject: [PATCH 3/3] [test][1p][convergence] remove obsolete file closure
 statements

---
 test/porousmediumflow/1p/implicit/convergence/main.cc | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/test/porousmediumflow/1p/implicit/convergence/main.cc b/test/porousmediumflow/1p/implicit/convergence/main.cc
index 3f93cbca75..ab4ddf5eac 100644
--- a/test/porousmediumflow/1p/implicit/convergence/main.cc
+++ b/test/porousmediumflow/1p/implicit/convergence/main.cc
@@ -175,9 +175,6 @@ int main(int argc, char** argv) try
             analyticFile << dxVector[i] << "," << analyticErrors[i] << "\n";
         for (unsigned int i = 0; i < discreteErrors.size(); ++i)
             discreteFile << dxVector[i] << "," << discreteErrors[i] << "\n";
-
-        analyticFile.close();
-        discreteFile.close();
     }
 
     // show parameters used
-- 
GitLab