[tests][1p] Add draft for anisotropic convergence test case

test/porousmediumflow/1p/implicit/convergence/CMakeLists.txt

 add_subdirectory(analyticsolution)
 add_subdirectory(discretesolution)
+add_subdirectory(unstructured)

test/porousmediumflow/1p/implicit/convergence/unstructured/CMakeLists.txt

+add_input_file_links()
+dune_symlink_to_source_files(FILES "convergencetest.py")
+
+dumux_add_test(NAME test_1p_convergence_structured
+              LABELS porousmediumflow 1p
+              SOURCES main.cc
+              COMPILE_DEFINITIONS GRIDTYPE=Dune::YaspGrid<2>
+              COMPILE_DEFINITIONS TYPETAG=OnePConvergence
+              TIMEOUT 1000
+              CMAKE_GUARD HAVE_UMFPACK
+              COMMAND ./convergencetest.py
+              CMD_ARGS test_1p_convergence_structured params.input
+                       -Problem.Name test_1p_convergence_structured)
+
+dumux_add_test(NAME test_1p_convergence_unstructured
+              LABELS porousmediumflow 1p
+              SOURCES main.cc
+              COMPILE_DEFINITIONS GRIDTYPE=Dune::UGGrid<2>
+              COMPILE_DEFINITIONS TYPETAG=OnePConvergence
+              TIMEOUT 1000
+              CMAKE_GUARD HAVE_UMFPACK
+              COMMAND ./convergencetest.py
+              CMD_ARGS test_1p_convergence_unstructured params.input
+                       -Problem.Name test_1p_convergence_unstructured
+                       -Grid.File ../../incompressible/grids/randomlydistorted.dgf)

test/porousmediumflow/1p/implicit/convergence/unstructured/convergencetest.py

+#!/usr/bin/env python3
+
+from math import *
+import subprocess
+import sys
+
+if len(sys.argv) < 2:
+    sys.stderr.write('Please provide a single argument <testname> to the script\n')
+    sys.exit(1)
+
+executableName = str(sys.argv[1])
+testargs = [str(i) for i in sys.argv][2:]
+testname = testargs[2]
+
+# remove the old log files
+subprocess.call(['rm', testname + '.log'])
+print("Removed old log file ({})!".format(testname + '.log'))
+
+# do the runs with different refinement
+for i in [0, 1, 2, 3]:
+    subprocess.call(['./' + executableName] + testargs + ['-Grid.Refinement', str(i)])
+
+def checkRates():
+    # check the rates and append them to the log file
+    logfile = open(testname + '.log', "r+")
+
+    errorP = []
+    for line in logfile:
+        line = line.strip("\n")
+        line = line.strip("\[ConvergenceTest\]")
+        line = line.split()
+        errorP.append(float(line[2]))
+
+    resultsP = []
+    logfile.truncate(0)
+    logfile.write("n\terrorP\t\trateP\n")
+    logfile.write("-"*50 + "\n")
+    for i in range(len(errorP)-1):
+        if isnan(errorP[i]) or isinf(errorP[i]):
+            continue
+        if not ((errorP[i] < 1e-12 or errorP[i+1] < 1e-12)):
+            rateP = (log(errorP[i])-log(errorP[i+1]))/log(2)
+            message = "{}\t{:0.4e}\t{:0.4e}\n".format(i, errorP[i], rateP)
+            logfile.write(message)
+            resultsP.append(rateP)
+        else:
+            logfile.write("error: exact solution!?")
+    i = len(errorP)-1
+    message = "{}\t{:0.4e}\n".format(i, errorP[i], "")
+    logfile.write(message)
+
+    logfile.close()
+    print("\nComputed the following convergence rates for {}:\n".format(testname))
+
+    subprocess.call(['cat', testname + '_darcy.log'])
+
+    return {"p" : resultsP}
+
+def checkConvRates():
+    rates = checkRates()
+
+    def mean(numbers):
+        return float(sum(numbers)) / len(numbers)
+
+    # check the rates, we expect rates around 2
+    if mean(rates["p"]) < 2.05 and mean(rates["p"]) < 1.95:
+        sys.stderr.write("*"*70 + "\n" + "The convergence rates for pressure were not close enough to 2! Test failed.\n" + "*"*70 + "\n")
+        sys.exit(1)
+
+
+checkConvRates()

test/porousmediumflow/1p/implicit/convergence/unstructured/main.cc

+// -*- 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 Test for the one-phase CC model
+ */
+
+#include <config.h>
+
+#include <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/fvector.hh>
+
+#include <dumux/nonlinear/newtonsolver.hh>
+#include <dumux/linear/seqsolverbackend.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 "problem.hh"
+
+/*!
+* \brief Creates analytical solution.
+* Returns a tuple of the analytical solution for the pressure, the velocity and the velocity at the faces
+* \param problem the problem for which to evaluate the analytical solution
+*/
+template<class Scalar, class Problem>
+auto createAnalyticalSolution(const Problem& problem)
+{
+    const auto& gridGeometry = problem.gridGeometry();
+    using GridView = typename std::decay_t<decltype(gridGeometry)>::GridView;
+
+    static constexpr auto dim = GridView::dimension;
+    static constexpr auto dimWorld = GridView::dimensionworld;
+
+    using VelocityVector = Dune::FieldVector<Scalar, dimWorld>;
+
+    std::vector<Scalar> analyticalPressure;
+    std::vector<VelocityVector> analyticalVelocity;
+
+    analyticalPressure.resize(gridGeometry.numDofs());
+    analyticalVelocity.resize(gridGeometry.numDofs());
+
+    for (const auto& element : elements(gridGeometry.gridView()))
+    {
+        auto fvGeometry = localView(gridGeometry);
+        fvGeometry.bindElement(element);
+        for (auto&& scv : scvs(fvGeometry))
+        {
+            const auto ccDofIdx = scv.dofIndex();
+            const auto ccDofPosition = scv.dofPosition();
+            const auto analyticalSolutionAtCc = problem.analyticalSolution(ccDofPosition);
+            analyticalPressure[ccDofIdx] = analyticalSolutionAtCc[dim];
+
+            for (int dirIdx = 0; dirIdx < dim; ++dirIdx)
+                analyticalVelocity[ccDofIdx][dirIdx] = analyticalSolutionAtCc[dirIdx];
+        }
+    }
+
+    return std::make_tuple(analyticalPressure, analyticalVelocity);
+}
+
+template<class Problem, class SolutionVector>
+void printL2Error(const Problem& problem, const SolutionVector& x)
+{
+    using namespace Dumux;
+    using Scalar = double;
+
+    Scalar l2error = 0.0;
+
+    for (const auto& element : elements(problem.gridGeometry().gridView()))
+    {
+        auto fvGeometry = localView(problem.gridGeometry());
+        fvGeometry.bindElement(element);
+
+        for (auto&& scv : scvs(fvGeometry))
+        {
+            const auto dofIdx = scv.dofIndex();
+            const Scalar delta = x[dofIdx] - problem.analyticalSolution(scv.center())[2/*pressureIdx*/];
+            l2error += scv.volume()*(delta*delta);
+        }
+    }
+    using std::sqrt;
+    l2error = sqrt(l2error);
+
+    const auto numDofs = problem.gridGeometry().numDofs();
+    std::ostream tmp(std::cout.rdbuf());
+    tmp << std::setprecision(8) << "** L2 error (abs) for "
+            << std::setw(6) << numDofs << " cc dofs "
+            << std::scientific
+            << "L2 error = " << l2error
+            << std::endl;
+
+    // write the norm into a log file
+    std::ofstream logFile;
+    logFile.open(problem.name() + ".log", std::ios::app);
+    logFile << "[ConvergenceTest] L2(p) = " << l2error << std::endl;
+    logFile.close();
+}
+
+int main(int argc, char** argv)
+{
+    using namespace Dumux;
+
+    using TypeTag = Properties::TTag::TYPETAG;
+
+    // initialize MPI, finalize is done automatically on exit
+    const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+    ////////////////////////////////////////////////////////////
+    // parse the command line arguments and input file
+    ////////////////////////////////////////////////////////////
+
+    // parse command line arguments
+    Parameters::init(argc, argv);
+
+    //////////////////////////////////////////////////////////////////////
+    // try to create a grid (from the given grid file or the input file)
+    /////////////////////////////////////////////////////////////////////
+    using Grid = GetPropType<TypeTag, Properties::Grid>;
+    GridManager<Grid> gridManager;
+    gridManager.init();
+
+    // we compute on the leaf grid view
+    const auto& leafGridView = gridManager.grid().leafGridView();
+
+    // create the finite volume grid geometry
+    using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+    auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);
+    gridGeometry->update();
+
+    // the problem (boundary conditions)
+    using Problem = GetPropType<TypeTag, Properties::Problem>;
+    auto problem = std::make_shared<Problem>(gridGeometry);
+
+    // the solution vector
+    using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
+    SolutionVector x(gridGeometry->numDofs());
+
+    // the grid variables
+    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
+    auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
+    gridVariables->init(x);
+
+    // intialize the vtk output module
+    VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
+    using VelocityOutput = GetPropType<TypeTag, Properties::VelocityOutput>;
+    vtkWriter.addVelocityOutput(std::make_shared<VelocityOutput>(*gridVariables));
+    using IOFields = GetPropType<TypeTag, Properties::IOFields>;
+    IOFields::initOutputModule(vtkWriter); // Add model specific output fields
+    vtkWriter.write(0.0);
+
+    // create assembler & linear solver
+    using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
+    auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables);
+
+    using LinearSolver = UMFPackBackend;
+    auto linearSolver = std::make_shared<LinearSolver>();
+
+    // the non-linear solver
+    using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
+    NewtonSolver nonLinearSolver(assembler, linearSolver);
+
+    // linearize & solve
+    nonLinearSolver.solve(x);
+
+    // output result to vtk
+    vtkWriter.write(1.0);
+
+    printL2Error(*problem, x);
+
+    if (mpiHelper.rank() == 0)
+        Parameters::print();
+
+    return 0;
+
+}

test/porousmediumflow/1p/implicit/convergence/unstructured/params.input

+[Grid]
+UpperRight = 1 1
+Cells = 40 40
+Refinement = 0
+
+[Problem]
+Name = convergence
+
+[Problem]
+EnableGravity = false
+
+[Component]
+LiquidDensity = 1.0
+LiquidKinematicViscosity = 1.0
+
+[Assembly]
+NumericDifference.BaseEpsilon = 1e-4

test/porousmediumflow/1p/implicit/convergence/unstructured/problem.hh

+// -*- 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 BoundaryTests
+ * \brief The convergence test
+ */
+
+#ifndef DUMUX_CONVERGENCE_TEST_ONEP_PROBLEM_HH
+#define DUMUX_CONVERGENCE_TEST_ONEP_PROBLEM_HH
+
+#include <dune/common/fvector.hh>
+#include <dune/grid/yaspgrid.hh>
+
+#include <dumux/discretization/cctpfa.hh>
+#include <dumux/discretization/ccmpfa.hh>
+#include <dumux/discretization/ccwmpfa.hh>
+#include <dumux/discretization/cellcentered/wmpfa/methods.hh>
+#include <dumux/common/boundarytypes.hh>
+
+#include <dumux/material/components/constant.hh>
+#include <dumux/material/fluidsystems/1pliquid.hh>
+
+#include <dumux/porousmediumflow/1p/model.hh>
+#include <dumux/porousmediumflow/problem.hh>
+
+#include "spatialparams.hh"
+
+namespace Dumux {
+template <class TypeTag>
+class ConvergenceProblem;
+
+namespace Properties {
+// Create new type tags
+namespace TTag {
+struct OnePConvergence { using InheritsFrom = std::tuple<OneP, CCWMpfaModel>; };
+} // end namespace TTag
+
+// Set the problem property
+template<class TypeTag>
+struct Problem<TypeTag, TTag::OnePConvergence> { using type = Dumux::ConvergenceProblem<TypeTag>; };
+
+// the fluid system
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::OnePConvergence>
+{
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using type = FluidSystems::OnePLiquid<Scalar, Dumux::Components::Constant<1, Scalar> > ;
+};
+
+// Set the grid type
+template<class TypeTag>
+struct Grid<TypeTag, TTag::OnePConvergence> { using type = GRIDTYPE; };
+
+template<class TypeTag>
+struct SpatialParams<TypeTag, TTag::OnePConvergence>
+{
+    using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using type = ConvergenceTestSpatialParams<GridGeometry, Scalar>;
+};
+
+template<class TypeTag>
+struct EnableGridVolumeVariablesCache<TypeTag, TTag::OnePConvergence> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridFluxVariablesCache<TypeTag, TTag::OnePConvergence> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridGeometryCache<TypeTag, TTag::OnePConvergence> { static constexpr bool value = true; };
+
+template<class TypeTag>
+struct DiscretizationSubmethod<TypeTag, TTag::OnePConvergence> { static constexpr WMpfaMethod value = WMpfaMethod::nltpfa; };
+
+} // end namespace Properties
+
+/*!
+ * \ingroup BoundaryTests
+ * \brief The convergence test
+ */
+template <class TypeTag>
+class ConvergenceProblem : public PorousMediumFlowProblem<TypeTag>
+{
+    using ParentType = PorousMediumFlowProblem<TypeTag>;
+    using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
+    using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
+    using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
+    using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
+    using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+    using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+    using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+    using Element = typename GridView::template Codim<0>::Entity;
+    using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+    static constexpr auto velocityXIdx = 0;
+    static constexpr auto velocityYIdx = 1;
+    static constexpr auto pressureIdx = 2;
+
+public:
+    //! export the Indices
+    using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+
+    ConvergenceProblem(std::shared_ptr<const GridGeometry> gridGeometry)
+    : ParentType(gridGeometry)
+    {}
+
+    /*!
+     * \name Problem parameters
+     */
+    // \{
+
+    /*!
+     * \brief Returns the temperature within the domain in [K].
+     *
+     */
+    Scalar temperature() const
+    { return 273.15 + 10; } // 10°C
+    // \}
+
+    /*!
+     * \name Boundary conditions
+     */
+    // \{
+
+    /*!
+      * \brief Specifies which kind of boundary condition should be
+      *        used for which equation on a given boundary control volume.
+      *
+      * \param element The element
+      * \param scvf The boundary sub control volume face
+      */
+    BoundaryTypes boundaryTypes(const Element &element, const SubControlVolumeFace &scvf) const
+    {
+        BoundaryTypes values;
+
+        values.setAllDirichlet();
+
+        return values;
+    }
+
+        /*!
+     * \brief Evaluates the boundary conditions for a Dirichlet control volume.
+     *
+     * \param element The element for which the Dirichlet boundary condition is set
+     * \param scvf The boundary subcontrolvolumeface
+     *
+     * For this method, the \a values parameter stores primary variables.
+     */
+    PrimaryVariables dirichlet(const Element &element, const SubControlVolumeFace &scvf) const
+    {
+        const auto p = analyticalSolution(scvf.center())[pressureIdx];
+        return PrimaryVariables(p);
+    }
+
+    // \}
+
+    /*!
+     * \name Volume terms
+     */
+    // \{
+    /*!
+     * \brief Evaluates the source term for all phases within a given
+     *        sub control volume.
+     *
+     * \param globalPos The global position
+     */
+    NumEqVector sourceAtPos(const GlobalPosition& globalPos) const
+    {
+        const Scalar x = globalPos[0];
+        const Scalar y = globalPos[1];
+        using std::exp; using std::sin; using std::cos;
+        static constexpr Scalar omega = M_PI;
+        static constexpr Scalar c = 0.9;
+        const Scalar cosOmegaX = cos(omega*x);
+        static const Scalar expTwo = exp(2);
+        const Scalar expYPlusOne = exp(y+1);
+
+        const Scalar result = (-(c*cosOmegaX + 1)*exp(y - 1)
+                                             + 1.5*c*expYPlusOne*cosOmegaX
+                                             + omega*omega*(expYPlusOne - expTwo + 2))
+                                             *sin(omega*x);
+        return NumEqVector(result);
+    }
+
+    // \}
+
+    /*!
+     * \brief Evaluates the initial value for a control volume.
+     *
+     * \param element The element
+     *
+     * For this method, the \a priVars parameter stores primary
+     * variables.
+     */
+    PrimaryVariables initial(const Element &element) const
+    {
+        return PrimaryVariables(0.0);
+    }
+
+    /*!
+     * \brief Returns the analytical solution of the problem at a given position.
+     *
+     * \param globalPos The global position
+     */
+    auto analyticalSolution(const GlobalPosition& globalPos) const
+    {
+        Dune::FieldVector<Scalar, 3> sol(0.0);
+        const Scalar x = globalPos[0];
+        const Scalar y = globalPos[1];
+        static constexpr Scalar omega = M_PI;
+        static constexpr Scalar c = 0.9;
+        using std::exp; using std::sin; using std::cos;
+        const Scalar sinOmegaX = sin(omega*x);
+        const Scalar cosOmegaX = cos(omega*x);
+        static const Scalar expTwo = exp(2);
+        const Scalar expYPlusOne = exp(y+1);
+
+        sol[pressureIdx] = (expYPlusOne + 2 - expTwo)*sinOmegaX + 10.0;
+        sol[velocityXIdx] = c/(2*omega)*expYPlusOne*sinOmegaX*sinOmegaX
+                            -omega*(expYPlusOne + 2 - expTwo)*cosOmegaX;
+        sol[velocityYIdx] = (0.5*c*(expYPlusOne + 2 - expTwo)*cosOmegaX
+                            -(c*cosOmegaX + 1)*exp(y-1))*sinOmegaX;
+
+        return sol;
+    }
+
+    // \}
+
+private:
+    static constexpr Scalar eps_ = 1e-7;
+};
+} // end namespace Dumux
+
+#endif //DUMUX_DARCY_SUBPROBLEM_HH

test/porousmediumflow/1p/implicit/convergence/unstructured/spatialparams.hh

+// -*- 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.                      *