Commit c2c3f958 authored by Martin Schneider's avatar Martin Schneider
Browse files

[tests][wmpfa] Add wmpfa to analytic convergence test and delete obsolete old test

parent 7de9c413
add_subdirectory(analyticsolution)
add_subdirectory(discretesolution)
add_subdirectory(unstructured)
......@@ -6,15 +6,18 @@ dune_symlink_to_source_files(FILES "convergencetest.py")
# the tests on unstructured grids to be skipped
add_executable(test_1p_convergence_analytic_tpfa EXCLUDE_FROM_ALL main.cc)
add_executable(test_1p_convergence_analytic_mpfa EXCLUDE_FROM_ALL main.cc)
add_executable(test_1p_convergence_analytic_wmpfa EXCLUDE_FROM_ALL main.cc)
add_executable(test_1p_convergence_analytic_box EXCLUDE_FROM_ALL main.cc)
if (dune-uggrid_FOUND)
target_compile_definitions(test_1p_convergence_analytic_tpfa PUBLIC "TYPETAG=OnePConvergenceTpfa" PUBLIC "GRIDTYPE=Dune::UGGrid<2>")
target_compile_definitions(test_1p_convergence_analytic_mpfa PUBLIC "TYPETAG=OnePConvergenceMpfa" PUBLIC "GRIDTYPE=Dune::UGGrid<2>")
target_compile_definitions(test_1p_convergence_analytic_wmpfa PUBLIC "TYPETAG=OnePConvergenceWMpfa" PUBLIC "GRIDTYPE=Dune::UGGrid<2>")
target_compile_definitions(test_1p_convergence_analytic_box PUBLIC "TYPETAG=OnePConvergenceBox" PUBLIC "GRIDTYPE=Dune::UGGrid<2>")
else()
target_compile_definitions(test_1p_convergence_analytic_tpfa PUBLIC "TYPETAG=OnePConvergenceTpfa" PUBLIC "GRIDTYPE=Dune::YaspGrid<2>")
target_compile_definitions(test_1p_convergence_analytic_mpfa PUBLIC "TYPETAG=OnePConvergenceMpfa" PUBLIC "GRIDTYPE=Dune::YaspGrid<2>")
target_compile_definitions(test_1p_convergence_analytic_wmpfa PUBLIC "TYPETAG=OnePConvergenceWMpfa" PUBLIC "GRIDTYPE=Dune::YaspGrid<2>")
target_compile_definitions(test_1p_convergence_analytic_box PUBLIC "TYPETAG=OnePConvergenceBox" PUBLIC "GRIDTYPE=Dune::YaspGrid<2>")
endif()
......@@ -33,6 +36,13 @@ dumux_add_test(NAME test_1p_convergence_analytic_mpfa_structured
COMMAND ./convergencetest.py
CMD_ARGS test_1p_convergence_analytic_mpfa mpfa_structured)
# using wmpfa and structured grid
dumux_add_test(NAME test_1p_convergence_analytic_wmpfa_structured
TARGET test_1p_convergence_analytic_wmpfa
LABELS porousmediumflow 1p
COMMAND ./convergencetest.py
CMD_ARGS test_1p_convergence_analytic_wmpfa wmpfa_structured)
# using box and structured grid
dumux_add_test(NAME test_1p_convergence_analytic_box_structured
TARGET test_1p_convergence_analytic_box
......@@ -50,6 +60,14 @@ if (dune-uggrid_FOUND)
CMD_ARGS test_1p_convergence_analytic_mpfa mpfa_unstructured
params.input -Grid.File ../../incompressible/grids/randomlydistorted.dgf)
# using wmpfa and unstructured grid
dumux_add_test(NAME test_1p_convergence_analytic_wmpfa_unstructured
TARGET test_1p_convergence_analytic_wmpfa
LABELS porousmediumflow 1p
COMMAND ./convergencetest.py
CMD_ARGS test_1p_convergence_analytic_wmpfa wmpfa_unstructured
params.input -Grid.File ../../incompressible/grids/randomlydistorted.dgf)
# using box and unstructured grid
dumux_add_test(NAME test_1p_convergence_analytic_box_unstructured
TARGET test_1p_convergence_analytic_box
......@@ -65,6 +83,11 @@ else()
CMAKE_GUARD dune-uggrid_FOUND
LABELS porousmediumflow 1p)
dumux_add_test(NAME test_1p_convergence_analytic_wmpfa_unstructured
SOURCES main.cc
CMAKE_GUARD dune-uggrid_FOUND
LABELS porousmediumflow 1p)
dumux_add_test(NAME test_1p_convergence_analytic_box_unstructured
SOURCES main.cc
CMAKE_GUARD dune-uggrid_FOUND
......
......@@ -69,6 +69,6 @@ def mean(numbers):
return float(sum(numbers)) / len(numbers)
# check the rates, we expect rates around 2
if mean(rates["p"]) < 1.8:
if mean(rates["p"]) < 1.75:
sys.stderr.write("*"*70 + "\n" + "The convergence rates for pressure were not close enough to 2! Test failed.\n" + "*"*70 + "\n")
sys.exit(1)
......@@ -31,6 +31,8 @@
#include <dumux/discretization/cctpfa.hh>
#include <dumux/discretization/ccmpfa.hh>
#include <dumux/discretization/ccwmpfa.hh>
#include <dumux/discretization/cellcentered/wmpfa/methods.hh>
#include <dumux/discretization/box.hh>
#include <dumux/material/components/constant.hh>
......@@ -47,6 +49,7 @@ namespace TTag {
struct OnePConvergence { using InheritsFrom = std::tuple<OneP>; };
struct OnePConvergenceTpfa { using InheritsFrom = std::tuple<OnePConvergence, CCTpfaModel>; };
struct OnePConvergenceMpfa { using InheritsFrom = std::tuple<OnePConvergence, CCMpfaModel>; };
struct OnePConvergenceWMpfa { using InheritsFrom = std::tuple<OnePConvergence, CCWMpfaModel>; };
struct OnePConvergenceBox { using InheritsFrom = std::tuple<OnePConvergence, BoxModel>; };
} // end namespace TTag
......@@ -82,6 +85,9 @@ struct EnableGridFluxVariablesCache<TypeTag, TTag::OnePConvergence> { static con
template<class TypeTag>
struct EnableGridGeometryCache<TypeTag, TTag::OnePConvergence> { static constexpr bool value = true; };
template<class TypeTag>
struct DiscretizationSubmethod<TypeTag, TTag::OnePConvergenceWMpfa> { static constexpr WMpfaMethod value = WMpfaMethod::nltpfa; };
} // end namespace Dumux::Properties
#endif
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)
#!/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 + '.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()
// -*- 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;
}
[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
// -*- 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
*