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