Commit f5c35a65 authored by Lukas Riedel's avatar Lukas Riedel

Merge branch 'feature/test-mass-conservation' into 'master'

Add test for mass conservation

Closes #39

See merge request !28
parents 536610da 3d61ca2f
......@@ -2,6 +2,7 @@ image: dorie/dune-env:2.5.1
variables:
DUNE_CONTROL_PATH: /opt/dune:$CI_PROJECT_DIR
DUNECONTROL: ./dune-common/bin/dunecontrol
before_script:
- cd /opt/dune
......@@ -16,7 +17,8 @@ build:main:
tags:
- demeter
script:
- CMAKE_FLAGS="-DDUNE_PYTHON_VIRTUALENV_SETUP=True -DDUNE_PYTHON_ALLOW_GET_PIP=True" MAKE_FLAGS="-j 2" ./dune-common/bin/dunecontrol --only=dorie all
- CMAKE_FLAGS="-DDUNE_PYTHON_VIRTUALENV_SETUP=True -DDUNE_PYTHON_ALLOW_GET_PIP=True" MAKE_FLAGS="-j 2" $DUNECONTROL --only=dorie all
- $DUNECONTROL --only=dorie make build_tests
artifacts:
name: "$CI_JOB_NAME"
paths:
......@@ -30,9 +32,10 @@ build:update_dune:
only:
- master
script:
- ./dune-common/bin/dunecontrol update
- ./dune-common/bin/dunecontrol exec "rm -rf build-cmake"
- CMAKE_FLAGS="-DDUNE_PYTHON_VIRTUALENV_SETUP=True -DDUNE_PYTHON_ALLOW_GET_PIP=True" MAKE_FLAGS="-j 2" ./dune-common/bin/dunecontrol --module=dorie all
- $DUNECONTROL update
- $DUNECONTROL exec "rm -rf build-cmake"
- CMAKE_FLAGS="-DDUNE_PYTHON_VIRTUALENV_SETUP=True -DDUNE_PYTHON_ALLOW_GET_PIP=True" MAKE_FLAGS="-j 2" $DUNECONTROL --module=dorie all
- $DUNECONTROL --only=dorie make build_tests
artifacts:
name: "$CI_JOB_NAME"
paths:
......@@ -44,7 +47,8 @@ build:debug:
tags:
- demeter
script:
- CMAKE_FLAGS="-DCMAKE_BUILD_TYPE=Debug -DDUNE_PYTHON_VIRTUALENV_SETUP=True -DDUNE_PYTHON_ALLOW_GET_PIP=True" MAKE_FLAGS="-j 2" ./dune-common/bin/dunecontrol --only=dorie all
- CMAKE_FLAGS="-DCMAKE_BUILD_TYPE=Debug -DDUNE_PYTHON_VIRTUALENV_SETUP=True -DDUNE_PYTHON_ALLOW_GET_PIP=True" MAKE_FLAGS="-j 2" $DUNECONTROL --only=dorie all
- $DUNECONTROL --only=dorie make build_tests
test:main:
stage: test
......@@ -54,8 +58,8 @@ test:main:
- build:main
allow_failure: true
script:
- ./dune-common/bin/dunecontrol --only=dorie configure
- ARGS="--output-on-failure -j 2" ./dune-common/bin/dunecontrol --only=dorie make test
- $DUNECONTROL --only=dorie configure
- ARGS="--output-on-failure -j 2" $DUNECONTROL --only=dorie make test
artifacts:
name: "$CI_JOB_NAME"
paths:
......@@ -72,8 +76,8 @@ test:update_dune:
- master
allow_failure: true
script:
- ./dune-common/bin/dunecontrol --only=dorie configure
- ARGS="--output-on-failure -j 2" ./dune-common/bin/dunecontrol --only=dorie make test
- $DUNECONTROL --only=dorie configure
- ARGS="--output-on-failure -j 2" $DUNECONTROL --only=dorie make test
artifacts:
name: "$CI_JOB_NAME"
paths:
......@@ -87,8 +91,8 @@ deploy:build_docs:
dependencies:
- build:main
script:
- ./dune-common/bin/dunecontrol --only=dorie configure
- ./dune-common/bin/dunecontrol --only=dorie make doc
- $DUNECONTROL --only=dorie configure
- $DUNECONTROL --only=dorie make doc
artifacts:
name: "$CI_JOB_NAME"
paths:
......
......@@ -32,9 +32,6 @@ dune_project()
dune_enable_all_packages()
dune_require_cxx_standard(MODULE "dorie" VERSION 14)
# avoid the executables from being 'excluded from all'
set(DUNE_BUILD_TESTS_ON_MAKE_ALL TRUE)
# add subdirectories
add_subdirectory("bin")
add_subdirectory("m4")
......
......@@ -233,8 +233,15 @@ A debugger needs special security privileges usually not provided by the Docker
The debugger `gdb` is pre-installed in the Docker container.
### DORiE is running, but I suspect that something is wrong
You can execute system tests in order to ensure that DORiE is running correctly and producing the expected results:
### Running System Tests
DORiE includes system tests for comparing its results the ones of ODE solvers or
former versions of itself. They ensure that DORiE is running correctly and
producing the expected results. If you installed DORiE, you have to make sure
that all test executables are built by calling
dunecontrol --only=dorie make build_tests
You then execute the system tests with
ARGS="--output-on-failure" dunecontrol --only=dorie make test
......
add_subdirectory(interface)
add_subdirectory(solver)
add_subdirectory(impl)
add_subdirectory(test)
add_executable("dorie" dorie.cc)
dune_target_link_libraries(dorie dorie-impl ${DUNE_LIBS})
\ No newline at end of file
#ifndef TEST_MASS_CONSERVATION_HH
#define TEST_MASS_CONSERVATION_HH
#include <algorithm>
#include <dune/geometry/referenceelements.hh>
#include <dune/pdelab/common/referenceelements.hh>
#include <dune/pdelab/common/quadraturerules.hh>
#include <dune/pdelab/localoperator/pattern.hh>
#include <dune/pdelab/localoperator/flags.hh>
#include <dune/pdelab/localoperator/idefault.hh>
#include <dune/pdelab/localoperator/defaultimp.hh>
#include <dune/pdelab/finiteelement/localbasiscache.hh>
#include <dune/dorie/dorie.hh>
#include <dune/dorie/interface/simulation.cc>
namespace Dune {
namespace Dorie {
/// Local Operator evaluating the water content of a grid cell
template <typename Traits, typename Parameter>
class WaterContentOperator : public Dune::PDELab::FullVolumePattern,
public Dune::PDELab::LocalOperatorDefaultFlags
{
public:
// residual assembly flags
enum { doAlphaVolume = true };
private:
static constexpr int dim = Traits::dim;
using RF = typename Traits::RF;
using Scalar = typename Traits::Scalar;
const Parameter &param; //!< class containing parameterization
const int quadrature_factor; //!< factor to FEM integration order
const int intorderadd; //!< integration order addend
public:
WaterContentOperator(
const Parameter &_param,
const int _quadrature_factor = 2,
const int _intorderadd = 2)
: param(_param),
quadrature_factor(_quadrature_factor),
intorderadd(_intorderadd)
{ }
template <typename EG, typename LFSU, typename X, typename LFSV, typename R>
void alpha_volume(
const EG &eg,
const LFSU &lfsu,
const X &x,
const LFSV &lfsv,
R &r) const
{
const int order = lfsu.finiteElement().localBasis().order();
const int intorder = intorderadd + quadrature_factor * order;
// get element geomentry
auto gt = eg.geometry();
// loop over quadrature points
for (const auto &it : quadratureRule(gt, intorder))
{
// evaluate position in element local and global coordinates
const auto p_local = it.position();
const auto p_global = gt.global(p_local);
// evaluate basis function
std::vector<Scalar> phi(lfsu.size());
lfsu.finiteElement().localBasis().evaluateFunction(p_local, phi);
// evaluate u
RF u = 0.;
for (unsigned int i = 0; i < lfsu.size(); i++)
u += x(lfsu, i) * phi[i];
// calculate water content from matric head
const RF head = param.headMillerToRef(u, p_global);
const RF water_content = param.waterContent(head, p_global);
// integration factor
const RF factor = it.weight() * gt.integrationElement(p_local);
// update residual
r.accumulate(lfsv, 0, water_content * factor);
} // quadrature rule
} // alpha_volume
};
/// Local Operator evaluating the boundary flux for each cell
/** This only includes fluxed through Neumann boundary conditions
*/
template <typename Traits, typename Parameter, typename Boundary>
class FluxOperator : public Dune::PDELab::FullVolumePattern,
public Dune::PDELab::LocalOperatorDefaultFlags,
public Dune::PDELab::InstationaryLocalOperatorDefaultMethods<typename Traits::RangeField>
{
public:
// residual assembly flags
enum { doLambdaBoundary = true };
private:
static constexpr int dim = Traits::dim;
using RF = typename Traits::RF;
using Scalar = typename Traits::Scalar;
const Parameter &param; //!< class containing parameterization
const Boundary &boundary; //!< class containing BC information
const int quadrature_factor; //!< factor to FEM integration order
const int intorderadd; //!< integration order addend
RF time;
public:
FluxOperator(
const Parameter &_param,
const Boundary &_boundary,
const int _quadrature_factor = 2,
const int _intorderadd = 2)
: param(_param),
boundary(_boundary),
quadrature_factor(_quadrature_factor),
intorderadd(_intorderadd)
{ }
void setTime(RF t) { time = t; }
template <typename IG, typename LFSV, typename R>
void lambda_boundary(
const IG &ig,
const LFSV &lfsv_s,
R &r_s) const
{
const int order = lfsv_s.finiteElement().localBasis().order();
const int intorder = intorderadd + quadrature_factor * order;
// get element geomentry
auto gtface = ig.geometry();
// loop over quadrature points
for (const auto &it : quadratureRule(gtface, intorder))
{
// check boundary condition type
const auto bc = boundary.bc(ig.intersection(), it.position(), time);
if (!BoundaryCondition::isNeumann(bc))
continue;
// evaluate flux boundary condition
const RF normal_flux = boundary.j(ig.intersection(), it.position(), time);
// integration factor
const RF factor = it.weight() * ig.geometry().integrationElement(it.position());
// update residual
r_s.accumulate(lfsv_s, 0, - normal_flux * factor);
} // quadrature rule
} // lambda_boundary
};
/// Derived Simulation for testing mass conservation. Assumes
template <typename Traits>
class TestSimulation : public Simulation<Traits>
{
protected:
// fetch typedefs from base class
using Sim = Simulation<Traits>;
using GV = typename Sim::GV;
using RF = typename Sim::RF;
using GFS = typename Sim::GFS;
using MBE = typename Sim::MBE;
using Grid = typename Sim::Grid;
using U = typename Sim::U;
using Parameters = typename Sim::Parameters;
using FlowBoundary = typename Sim::FlowBoundary;
/// Error Function Space Helper
using ERRGFSHelper = Dune::Dorie::LinearSolverGridFunctionSpaceHelper<GV, RF, Traits::GridGeometryType>;
/// Linear solver GFS
using ERRGFS = typename ERRGFSHelper::Type;
/// LSGFS Constraints Container
using ERRCC = typename ERRGFSHelper::CC;
/// Solution vector type
using U0 = Dune::PDELab::Backend::Vector<ERRGFS, RF>;
/// Empty constraints container
using NoTrafo = Dune::PDELab::EmptyTransformation;
/// Grid operator for Error LOP
template <typename LOP>
using ERRGO = Dune::PDELab::GridOperator<GFS, ERRGFS, LOP, MBE, RF, RF, RF, NoTrafo, NoTrafo>;
public:
/// Construct this simulation instance.
/** \todo Give error limits here
*/
TestSimulation(
Dune::MPIHelper &_helper,
std::shared_ptr<Grid> _grid,
Dune::ParameterTree &_inifile)
: Simulation<Traits>(_helper, _grid, _inifile)
{ }
/// Compute the water flux into a certain domain
double water_flux(U &solution, const RF time)
{
using FluxLOP = typename Dune::Dorie::FluxOperator<Traits, Parameters, FlowBoundary>;
using FluxGO = ERRGO<FluxLOP>;
// create GFS and operators
auto errgfs = ERRGFSHelper::create(this->gv);
FluxLOP fluxlop(*this->param, *this->fboundary);
fluxlop.setTime(time);
MBE mbe(0);
FluxGO fluxgo(*this->gfs, errgfs, fluxlop, mbe);
U0 eta(errgfs, 0.0);
// assemble residual and accumulate value
fluxgo.residual(solution, eta);
auto &eta_n = Dune::PDELab::Backend::native(eta);
return std::accumulate(eta_n.begin(), eta_n.end(), 0.0);
}
/// Compute the total water content of a certain solution
double water_content(U &solution)
{
using WCLOP = typename Dune::Dorie::WaterContentOperator<Traits, Parameters>;
using WCGO = ERRGO<WCLOP>;
// create GFS and operators
auto errgfs = ERRGFSHelper::create(this->gv);
WCLOP wclop(*this->param);
MBE mbe(0);
WCGO wcgo(*this->gfs, errgfs, wclop, mbe);
U0 eta(errgfs, 0.0);
// assemble residual and accumulate value
wcgo.residual(solution, eta);
auto &eta_n = Dune::PDELab::Backend::native(eta);
return std::accumulate(eta_n.begin(), eta_n.end(), 0.0);
}
/// Run the simulation
/** \return Error quantity
*/
double run_test()
{
const auto t_start = this->controller->getTime();
this->output->write_vtk_output(t_start);
double acc{0.0}; // accumulated deviation
double acc_square{0.0}; // accumulated squared deviation
while (this->controller->doStep())
{
const auto time = this->controller->getTime();
const auto dt = this->controller->getDT();
// compute current water content
const auto wc_old = water_content(*this->uold);
// compute expected integrated flux (water content change)
const auto flux_int = water_flux(*this->uold, time) * dt;
// compute time step
if (!this->compute_time_step())
{
continue;
}
if (this->adaptivity->adapt_grid(
*this->grid, this->gv, *this->gfs, *this->param, *this->fboundary,
time, *this->uold, *this->unew))
{
// reset operators if grid changes
this->operator_setup();
}
// compute new water content
const auto wc_new = water_content(*this->unew);
const auto deviation = (wc_new - wc_old) - flux_int;
acc += deviation;
acc_square += deviation * deviation;
std::cout << "wc_old: " << wc_old << std::endl;
std::cout << "wc_new: " << wc_new << std::endl;
std::cout << "integrated flux: " << flux_int << std::endl;
std::cout << "deviation: " << deviation << std::endl;
this->output->write_vtk_output(this->controller->getTime());
}
std::cout << "total deviation: " << acc << std::endl;
std::cout << "total normalized deviation: " << std::sqrt(acc_square) << std::endl;
return acc;
}
};
} // namespace Dorie
} // namespace Dune
#endif // TEST_MASS_CONSERVATION_HH
\ No newline at end of file
......@@ -35,6 +35,15 @@ dorie_add_system_test_dependency(dorie_parallel_reference_compare_0001 dorie_par
dorie_add_system_test_dependency(dorie_parallel_reference_compare_0002 dorie_parallel_reference_0001)
dorie_add_system_test_dependency(dorie_parallel_reference_compare_0003 dorie_parallel_reference_0001)
# dorie mass conservation
configure_file(mass_conservation.mini.in ${CMAKE_CURRENT_SOURCE_DIR}/mass_conservation.mini)
dune_add_system_test(
SOURCE test-mass-conservation.cc
BASENAME mass-conservation
CREATED_TARGETS test-mass-conservation
INIFILE mass_conservation.mini)
dune_target_link_libraries(${test-mass-conservation} ${DUNE_LIBS})
# dorie pfg
dorie_add_system_test(dorie-rfg parfield.mini)
dorie_add_system_test(dorie-rfg parfield_muphi.mini)
......@@ -47,3 +56,7 @@ dorie_add_system_test_dependency(dorie_plot dorie_reference_2d_0000)
# dorie create
dorie_add_system_test(dorie create.mini)
# dune excludes test targets from 'make all'; undo that here where applicable
set_property(TARGET dorie PROPERTY EXCLUDE_FROM_ALL 0)
set_property(TARGET dorie-rfg PROPERTY EXCLUDE_FROM_ALL 0)
\ No newline at end of file
spatial_resolution_north 0
spatial_resolution_south 0
spatial_resolution_west -1
spatial_resolution_east -1
number_BC_change_times 1
0 neumann -5e-6 neumann 0
\ No newline at end of file
include ${CMAKE_BINARY_DIR}/doc/default_files/config.ini
__name = mc
_asset_path = "${CMAKE_CURRENT_LIST_DIR}"
_limit_cube_static = 2E-9
_limit_cube_adapt = 2E-4
_limit_simplex_static = 5E-10
_limit_simplex_adapt = 3E-6
_prefix = _simplex
_suffix = _static, _adapt | expand adaptive
limit = {_limit{_prefix}{_suffix}}
output.fileName = mass_conservation | unique
output.outputPath = mass_conservation | unique
output.verbose = 0
boundary.file = "{_asset_path}/bcs/mass_conservation_2d.dat"
parameters.arrayFile = "{_asset_path}/parfields/sand.h5"
time.end = 1E4
time.maxTimestep = 1E4
time.startTimestep = 1E2
adaptivity.useAdaptivity = false, true | expand adaptive
adaptivity.maxLevel = 3
adaptivity.markingStrategy = threshold
adaptivity.refinementFraction = 5E-9
adaptivity.coarseningFraction = 5E-10
adaptivity.threshold = 0.0
grid.FEorder = 1
grid.gridType = simplex
grid.initialLevel = 2
grid.cells = 20 20
grid.gridFile = {_asset_path}/meshes/square.msh
\ No newline at end of file
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <dune/dorie/dorie.hh>
#include <dune/dorie/test/test-mass-conservation.hh>
namespace Dune{
namespace Dorie{
/// Resolve the second (default) template parameter of YaspGrid
template<int dim>
using YaspGrid = Dune::YaspGrid<dim>;
}
}
template<int dim, int order>
using Simplex = Dune::Dorie::BaseTraits<Dune::UGGrid,
Dune::GeometryType::BasicType::simplex,dim,order,false,false>;
template<int dim, int order>
using SimplexAdaptive = Dune::Dorie::BaseTraits<Dune::UGGrid,
Dune::GeometryType::BasicType::simplex,dim,order,false,true>;
template<int dim, int order>
using Cube = Dune::Dorie::BaseTraits<Dune::Dorie::YaspGrid,
Dune::GeometryType::BasicType::cube,dim,order,false,false>;
template<int dim, int order>
using CubeAdaptive = Dune::Dorie::BaseTraits<Dune::UGGrid,
Dune::GeometryType::BasicType::cube,dim,order,false,true>;
int main (int argc, char** argv)
{
// initialize MPI if needed
auto& helper = Dune::MPIHelper::instance(argc, argv);
// Read ini file
if (argc!=2)
DUNE_THROW(Dune::IOError,
"Call this program with arguments <config>");
const std::string inifilename = argv[1];
Dune::ParameterTree inifile;
Dune::ParameterTreeParser ptreeparser;
ptreeparser.readINITree(inifilename,inifile);
const double limit = inifile.get<double>("limit");
// run test
double result {0.0};
const auto adaptive = inifile.get<bool>("adaptivity.useAdaptivity");
const auto grid_type = inifile.get<std::string>("grid.gridType");
if (adaptive) {
if (grid_type == "rectangular") {
using TestSim = Dune::Dorie::TestSimulation<CubeAdaptive<2, 1>>;
auto grid = Dune::Dorie::build_grid_cube<Dune::UGGrid<2>>(inifile, helper);
TestSim sim(helper, grid, inifile);
result = sim.run_test();
}
else if (grid_type == "simplex") {
using TestSim = Dune::Dorie::TestSimulation<SimplexAdaptive<2, 1>>;
auto grid = Dune::Dorie::build_grid_gmsh<Dune::UGGrid<2>>(inifile, helper);
TestSim sim(helper, grid, inifile);
result = sim.run_test();
}
}
else {
if (grid_type == "rectangular") {
using TestSim = Dune::Dorie::TestSimulation<Cube<2, 1>>;
auto grid = Dune::Dorie::build_grid_cube<Dune::YaspGrid<2>>(inifile, helper);
TestSim sim(helper, grid, inifile);
result = sim.run_test();
}
else if (grid_type == "simplex") {
using TestSim = Dune::Dorie::TestSimulation<Simplex<2, 1>>;
auto grid = Dune::Dorie::build_grid_gmsh<Dune::UGGrid<2>>(inifile, helper);
TestSim sim(helper, grid, inifile);
result = sim.run_test();
}
}
if (std::abs(result) > limit) {
std::cerr << "Mass conservation test exceeded limit!" << std::endl;
std::cerr << "Limit: " << limit << std::endl;
std::cerr << "Result: " << result << std::endl;
return 1;
}
return 0;
}
\ No newline at end of file
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment