Commit 06897d0e authored by Timo Koch's avatar Timo Koch
Browse files

Merge branch 'feature/1p-convergence-test' into 'master'

[test][1p] add convergence test

See merge request !1665
parents f9ab560e 0d70cacc
add_subdirectory(convergence)
add_subdirectory(pointsources)
add_subdirectory(incompressible)
add_subdirectory(compressible)
......
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)
// -*- 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);
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;
}
[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
// -*- 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
// -*- 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