The TS-GitLab will have to shut down towards the end of the year — please think about migrating your projects to GitLab.com or GitHub.
(This is still a very early message, meant to keep you informed. There will be more urgent ones in the future.)

Commit 48b2aabe authored by Lukas Riedel's avatar Lukas Riedel
Browse files

Add mass conservation test to testing suite

Test executable is now run with new metaini file.
Limits are set closely to the current mass conservation results.
Moved class definitions into own header file.
Test executable is managed from within the `testing` directory

This fixes #39.
parent 3eec101c
add_executable("test-mass-conservation" test-mass-conservation.cc)
dune_target_link_libraries(test-mass-conservation ${DUNE_LIBS})
\ No newline at end of file
......@@ -3,306 +3,9 @@
#include "config.h"
#endif
#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.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))
{
// evaluate position in element local and global coordinates
const auto p_local = ig.geometryInInside().global(it.position());
const auto p_global = ig.inside().geometry().global(p_local);
// 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;
}
};
} // namespace Dorie
} // namespace Dune
#include "test-mass-conservation.hh"
namespace Dune{
namespace Dorie{
......@@ -314,19 +17,19 @@ namespace Dorie{
template<int dim, int order>
using Simplex = Dune::Dorie::BaseTraits<Dune::UGGrid,
Dune::GeometryType::BasicType::simplex,dim,order,true,false>;
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,true,true>;
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,true,false>;
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,true,true>;
Dune::GeometryType::BasicType::cube,dim,order,false,true>;
int main (int argc, char** argv)
{
......@@ -335,20 +38,54 @@ int main (int argc, char** argv)
// Read ini file
if (argc!=2)
DUNE_THROW(Dune::IOError,"No parameter file specified!");
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
// inifile["adaptivity.useAdaptivity"] = "true";
using TestSim = Dune::Dorie::TestSimulation<Cube<2,1>>;
auto grid = Dune::Dorie::build_grid_cube<Dune::YaspGrid<2>>(inifile, helper);
// 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);
sim.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
#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))
{
// evaluate position in element local and global coordinates
const auto p_local = ig.geometryInInside().global(it.position());
const auto p_global = ig.inside().geometry().global(p_local);
// 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: