Commit 7e9e1bc1 authored by Lukas Riedel's avatar Lukas Riedel 🎧
Browse files

Add simple executable for testing mass conservation

New test implements local operators for accumulating the total water content and Neumann fluxes across the boundaries, respectively. Test is not added to the system tests yet.
parent e1bb918c
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
add_executable("test-mass-conservation" test-mass-conservation.cc)
dune_target_link_libraries(test-mass-conservation ${DUNE_LIBS})
\ 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 <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
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,true,false>;
template<int dim, int order>
using SimplexAdaptive = Dune::Dorie::BaseTraits<Dune::UGGrid,
Dune::GeometryType::BasicType::simplex,dim,order,true,true>;
template<int dim, int order>
using Cube = Dune::Dorie::BaseTraits<Dune::Dorie::YaspGrid,
Dune::GeometryType::BasicType::cube,dim,order,true,false>;
template<int dim, int order>
using CubeAdaptive = Dune::Dorie::BaseTraits<Dune::UGGrid,
Dune::GeometryType::BasicType::cube,dim,order,true,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,"No parameter file specified!");
const std::string inifilename = argv[1];
Dune::ParameterTree inifile;
Dune::ParameterTreeParser ptreeparser;
ptreeparser.readINITree(inifilename,inifile);
// 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();
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