...
  View open merge request
Commits (13)
......@@ -268,7 +268,7 @@ adding an empty line, make text **bold** or ``monospaced``.
</parameter>
<parameter name="checkTol">
<definition> Whenever ``checkJumps`` is activated, it check that flux
<definition> Whenever ``checkJumps`` is activated, it checks that flux
reconstruction engine is creating conforming normal fluxes up to ``jumpTol``.
</definition>
<values> float &gt; 0 </values>
......@@ -276,6 +276,24 @@ adding an empty line, make text **bold** or ``monospaced``.
</parameter>
</category>
<category name="conservationOfMass">
<parameter name="check">
<definition> Check that mass is conserved in the domain with respect to
boundary fluxes.
</definition>
<suggestion> none </suggestion>
<values> none, warn, error </values>
</parameter>
<parameter name="checkTol">
<definition> Whenever ``checkJumps`` is activated, it check that mass is
conserved in the domain with respect to boundary fluxes. ``checkTol``
is the persentage of mass that can be lost without a warning/error.
</definition>
<values> float &gt; 0 </values>
<suggestion> 1E-10 </suggestion>
</parameter>
</category>
<category name="NewtonParameters">
<parameter name="ForceIteration" hidden="true">
......
......@@ -280,7 +280,7 @@ adding an empty line, make text **bold** or ``monospaced``.
</parameter>
<parameter name="checkTol">
<definition> Whenever ``checkJumps`` is activated, it check that flux
<definition> Whenever ``checkJumps`` is activated, it checks that flux
reconstruction engine is creating conforming normal fluxes up to ``jumpTol``.
</definition>
<values> float &gt; 0 </values>
......@@ -288,6 +288,25 @@ adding an empty line, make text **bold** or ``monospaced``.
</parameter>
</category>
<category name="conservationOfMass">
<parameter name="check">
<definition> Check that mass is conserved in the domain with respect to
boundary fluxes.
</definition>
<suggestion> none </suggestion>
<values> none, warn, error </values>
</parameter>
<parameter name="checkTol">
<definition> Whenever ``checkJumps`` is activated, it check that mass is
conserved in the domain with respect to boundary fluxes. ``checkTol``
is the persentage of mass that can be lost without a warning/error.
</definition>
<values> float &gt; 0 </values>
<suggestion> 1E-10 </suggestion>
</parameter>
</category>
<category name="solverParameters">
<parameter name="reduction">
<definition> Required relative defect reduction.
......
......@@ -283,9 +283,6 @@ public:
}
private:
/// Controls how to deal with jumps occurring in the reconstructed flux
enum class CheckLevel { none, warn, error };
std::shared_ptr<spdlog::logger> _log;
GV _gv;
FEMU _femu;
......@@ -297,25 +294,6 @@ private:
/// The jump check level
const CheckLevel _check_lvl;
const RF _check_tol;
/// Return a check level for a given string
CheckLevel to_check_lvl (const std::string check_lvl_str)
{
if (to_lower(check_lvl_str) == "none") {
return CheckLevel::none;
}
else if (check_lvl_str == "warn") {
return CheckLevel::warn;
}
else if (check_lvl_str == "error") {
return CheckLevel::error;
}
else {
_log->error("Unknown check for jumps in reconstructed flux: {}",
check_lvl_str);
DUNE_THROW(IOError, "Unknown check for flux reconstruction jumps");
}
}
};
} // namespace Dorie
......
#ifndef DUNE_DORIE_GRID_FUNCTION_INTEGRATOR_HH
#define DUNE_DORIE_GRID_FUNCTION_INTEGRATOR_HH
namespace Dune {
namespace Dorie {
/*-------------------------------------------------------------------------*//**
* @brief Integrate grid function
*
* @author Santiago Ospina De Los Ríos
* @date 2019
*
* @param[in] gf The grid function
* @param[in] int_order The integration order
*
* @tparam GF The grid function type
*
* @return The integral over the whole domain
*/
template<class GF>
typename GF::Traits::RangeType
gf_integrator(const GF& gf, unsigned int int_order)
{
using Traits = typename GF::Traits;
using Range = typename Traits::RangeType;
Range integral(0.);
const auto& gv = gf.getGridView();
// Loop over the grid
for (const auto& cell : Dune::elements(gv,Dune::Partitions::interior))
{
// Get geometry
auto geo = cell.geometry();
// Integrate on the cell
for (const auto& point : quadratureRule(geo, int_order))
{
// Evaluate function
Range y;
const auto& x = point.position();
gf.evaluate(cell,x,y);
auto factor = point.weight() * geo.integrationElement(x);
// Accumulate
integral += y*factor;
}
}
// Accumulate over processors
return gv.comm().sum(integral);
}
/*-------------------------------------------------------------------------*//**
* @brief Integrate grid function on the boundary of the domain
* @details If the grid function is a vector valued function, it will be
* integrated in the normal direction of the boundary.
*
* @author Santiago Ospina De Los Ríos
* @date 2019
*
* @param[in] gf The grid function
* @param[in] int_order The integration order
*
* @tparam GF The grid function type
*
* @return The integral over the whole domain boundary
*/
template<class GF>
Dune::FieldVector<typename GF::Traits::RangeFieldType,1>
boundary_gf_integrator(const GF& gf, unsigned int int_order)
{
using Traits = typename GF::Traits;
using Range = typename Traits::RangeType;
constexpr int dim_range = Traits::dimRange;
constexpr int dim_domain = Traits::dimDomain;
static_assert(dim_range == 1 || dim_range == dim_domain);
Dune::FieldVector<typename GF::Traits::RangeFieldType,1> integral(0.);
const auto& gv = gf.getGridView();
// Loop over the grid
for (const auto& cell : Dune::elements(gv,Dune::Partitions::interior))
{
// Loop over the intersections
for (const auto& face : Dune::intersections(gv,cell))
{
// Only integrate on boundary
if (face.boundary() and not face.neighbor())
{
// Get geometry
auto geo_f = face.geometry();
auto geo_in_i = face.geometryInInside();
// Integrate on the face
for (const auto& point : quadratureRule(geo_f, int_order))
{
// Relative positions
const auto& position_f = point.position();
const auto position_i = geo_in_i.global(position_f);
// Face normal vector
const auto normal = face.unitOuterNormal(position_f);
// Evaluate function
Range y;
gf.evaluate(cell,position_i,y);
const auto factor = point.weight() * geo_f.integrationElement(position_f);
// Accumulate
if constexpr (dim_range == dim_domain)
integral += y * normal * factor;
else
integral += y * factor;
}
}
}
}
// Accumulate over processors
return gv.comm().sum(integral);
}
}
}
#endif
#ifndef DUNE_DORIE_STATISTICAL_MOMENTS_HH
#define DUNE_DORIE_STATISTICAL_MOMENTS_HH
namespace Dune {
namespace Dorie {
template<class GF>
std::vector<typename GF::Traits::DomainType>
statistical_moments(
const GF& gf,
unsigned int int_order,
unsigned int max_order = 2)
{
using Traits = typename GF::Traits;
using Domain = typename Traits::DomainType;
constexpr unsigned int dim = Traits::dimDomain;
assert(max_order>0);
std::vector<Domain> moment(max_order);
const auto& gv = gf.getGridView();
// Loop over the grid
for (const auto& cell : Dune::elements(gv,Dune::Partitions::interior))
{
auto geo = cell.geometry();
for (const auto& point : quadratureRule(geo, int_order))
{
const auto& x = point.position();
auto x_gl = geo.global(x);
typename Traits::RangeType y;
gf.evaluate(cell,x,y);
auto factor = point.weight() * geo.integrationElement(x);
Domain weight(1.);
for (unsigned int i = 0; i < max_order; ++i) {
for (unsigned int j = 0; j < dim; ++j)
{
weight[j] *= x_gl[j];
moment[i][j] += weight[j]*y*factor;
}
}
}
}
// Accumulate over processors
for (unsigned int i = 0; i < max_order; ++i)
moment[i] = gv.comm().sum(moment[i]);
return moment;
}
}
}
#endif
......@@ -31,6 +31,10 @@ enum class OutputPolicy {None,EndOfRichardsStep,EndOfTransportStep};
*/
enum class AdaptivityPolicy {None,WaterFlux,SoluteFlux};
/// Controls how to deal with run-time checks
enum class CheckLevel { none, warn, error };
/// Return the estimation of entries per matrix row for the spatial GridOperator
/** This supposedly decreases matrix assembly time.
* The values specify the *blocks* per row. DG assembles one block for the
......
......@@ -11,6 +11,7 @@
#include <dune/common/fvector.hh>
#include <dune/grid/common/gridview.hh>
#include <dune/dorie/common/util.hh>
namespace Dune {
namespace Dorie {
......@@ -113,6 +114,23 @@ inline bool is_none (const std::string& in)
return false;
}
/// Return a check level for a given string
CheckLevel to_check_lvl (const std::string check_lvl_str)
{
if (to_lower(check_lvl_str) == "none") {
return CheckLevel::none;
}
else if (check_lvl_str == "warn") {
return CheckLevel::warn;
}
else if (check_lvl_str == "error") {
return CheckLevel::error;
}
else {
DUNE_THROW(IOError, "Unknown check level");
}
}
} // namespace Dorie
} // namespace Dune
......
#include <dune/pdelab/function/scalarscaled.hh>
#include <dune/dorie/model/richards/richards.hh>
#include <dune/dorie/common/grid_function_integrator.hh>
#include <dune/dorie/common/statistical_moments.hh>
#include <limits>
#include <cmath>
namespace Dune{
namespace Dorie{
......@@ -94,6 +100,16 @@ RichardsSimulation<Traits>::RichardsSimulation (
inifile.get<std::string>("output.fileName"),
inifile.get<std::string>("output.outputPath"),
"./");
// Configure statistics file
std::ofstream stats_file;
stats_file.open(inifile.get<std::string>("output.fileName") + ".dat", std::ofstream::trunc);
stats_file << "# Statistics of the Richards Equation generated by DORiE.\n";
stats_file << "#\n";
stats_file << "# run_time\t time\t mass\t boundary_flux\t";
stats_file << "mean_cw_x\t mean_cw_y\t mean_cw_z\t";
stats_file << "variance_cw_x\t variance_cw_y\t variance_cw_z\n";
stats_file.close();
}
// --- Adaptivity Setup --- //
......@@ -229,6 +245,8 @@ void RichardsSimulation<Traits>::step()
// invalidate cache for water flux reconstruction
cache_water_flux_gf_rt.reset();
time_steps++;
if (this->output_policy() == OutputPolicy::EndOfRichardsStep)
write_data();
}
......@@ -263,15 +281,102 @@ void RichardsSimulation<Traits>::post_adapt_grid()
template<typename Traits>
void RichardsSimulation<Traits>::write_data () const
{
if (vtkwriter)
const auto time = controller->getTime();
// create adapters
auto head = get_matric_head();
auto wflux = std::make_shared<GFWaterFlux>(gfs, u, fparam);
auto cond = std::make_shared<GFConductivity>(gv, fparam);
auto wc = std::make_shared<GFWaterContent>(head, gv, fparam);
auto sat = std::make_shared<GFSaturation>(head, gv, fparam);
// Compute mass of WC
RF mass = gf_integrator(*wc,2*order+2);
// Compute statistical moments of WC
using PDF = Dune::PDELab::ScalarScaledGridFunctionAdapter<GFWaterContent>;
PDF pdf(1./mass,*wc);
auto moments = statistical_moments(pdf,order*4,2);
auto moments_squared = moments[0];
for (int i = 0; i < dim; ++i) moments_squared[i]*=moments_squared[i];
const auto variance = moments[1] - moments_squared;
Dune::FieldVector<double,3> moments0_3d(0.), variance_3d(0.);
for (int i = 0; i < dim; ++i)
{
// create adapters
auto head = get_matric_head();
auto wflux = std::make_shared<GFWaterFlux>(gfs, u, fparam);
auto cond = std::make_shared<GFConductivity>(gv, fparam);
auto wc = std::make_shared<GFWaterContent>(head, gv, fparam);
auto sat = std::make_shared<GFSaturation>(head, gv, fparam);
moments0_3d[i] = moments[0][i];
variance_3d[i] = variance[i];
}
RF boundary_flux = std::numeric_limits<RF>::quiet_NaN();;
if constexpr (enable_rt_engine) {
if (enable_fluxrc) {
auto wfluxr = get_water_flux_reconstructed();
boundary_flux = boundary_gf_integrator(*wfluxr,2*flux_order+2);
}
}
if (time_steps==0) cache_time = time;
if (time_steps==0) cache_mass = mass;
if (time_steps==0) cache_boundary_flux = boundary_flux;
// hard check on mass
if (std::isnan(mass)) {
this->_log->error("Mass of water is NaN: PDE solver has diverged.");
DUNE_THROW(MathError,"PDE solver has diverged!");
}
// check relative change of mass
if (Dune::FloatCmp::ne(mass,0.))
{
// change in time
auto dt = time - cache_time;
// change of mass by PDE solution
auto mass_dt = (mass - cache_mass);
// change of mass by boundary conditions
mass_dt += 0.5*dt*(boundary_flux + cache_boundary_flux);
// relative mass change
auto rel_mass_dt = mass_dt/mass;
auto mass_tol = inifile.get<double>("conservationOfMass.checkTol");
auto check_mass_str = inifile.get<std::string>("conservationOfMass.check");
auto check_mass_lvl = Dorie::to_check_lvl(check_mass_str);
const bool passed = rel_mass_dt < mass_tol;
if (passed) {
this->_log->trace(" Conservation of mass does not hold for"
" {}% of mass",
rel_mass_dt);
}
else if (not passed && check_mass_lvl == CheckLevel::warn) {
this->_log->warn("Conservation of mass above tolerance:"
" {}% of mass is not conserved.",
rel_mass_dt);
}
else if (not passed && check_mass_lvl == CheckLevel::error) {
this->_log->error("Conservation of mass above tolerance:"
" {}% of mass is not conserved.",
rel_mass_dt);
DUNE_THROW(MathError, "Conservation of mass does not obey tolerance");
}
}
// update cache
cache_time = time;
cache_mass = mass;
cache_boundary_flux = boundary_flux;
// Configure statistics file
std::ofstream stats_file;
stats_file.open(inifile.get<std::string>("output.fileName") + ".dat", std::ofstream::app);
stats_file << timer.elapsed() << "\t" << time << "\t" << mass << "\t";
stats_file << boundary_flux << "\t" << moments0_3d << "\t" << variance_3d;
stats_file << std::endl;
stats_file.close();
if (vtkwriter)
{
if (inifile.get<bool>("output.vertexData")) {
vtkwriter->template addVertexData<GFMatricHead>(head,"head");
vtkwriter->template addVertexData<GFWaterFlux>(wflux,"flux");
......@@ -299,7 +404,6 @@ void RichardsSimulation<Traits>::write_data () const
}
try{
const auto time = controller->getTime();
this->_log->trace("Writing solution at time {:.2e}", time);
vtkwriter->write(time, output_type);
}
......
......@@ -15,6 +15,7 @@
#include <dune/dorie/common/simulation_base.hh>
#include <dune/dorie/common/util.hh>
#include <dune/dorie/common/utility.hh>
#include <dune/dorie/common/interpolator.hh>
#include <dune/dorie/common/time_controller.hh>
#include <dune/dorie/common/grid_creator.hh>
......@@ -190,6 +191,8 @@ private:
static constexpr int order = Traits::order;
static constexpr int flux_order = Traits::flux_order;
static constexpr bool enable_rt_engine = Traits::enable_rt_engine;
using TimeField = typename Traits::TimeField;
using RF = typename Traits::RF;
using Grid = typename Traits::Grid;
using GV = typename Traits::GV;
......@@ -316,12 +319,18 @@ protected:
std::unique_ptr<Writer> vtkwriter;
std::unique_ptr<AdaptivityHandler> adaptivity;
mutable std::shared_ptr<GFFluxReconstruction> cache_water_flux_gf_rt;
Dune::Timer timer;
RF time_before;
RF dt_before;
const bool enable_fluxrc;
unsigned int time_steps;
mutable std::shared_ptr<GFFluxReconstruction> cache_water_flux_gf_rt;
mutable RF cache_mass;
mutable RF cache_boundary_flux;
mutable TimeField cache_time;
public:
/*------------------------------------------------------------------------*//**
......
#include <dune/pdelab/function/scalarscaled.hh>
#include <dune/dorie/model/transport/transport.hh>
#include <dune/dorie/common/grid_function_integrator.hh>
#include <dune/dorie/common/statistical_moments.hh>
#include <limits>
#include <cmath>
namespace Dune{
namespace Dorie{
......@@ -81,6 +88,16 @@ TransportSimulation<Traits>::TransportSimulation(
inifile.get<std::string>("output.fileName"),
inifile.get<std::string>("output.outputPath"),
"./");
// Configure statistics file
std::ofstream stats_file;
stats_file.open(inifile.get<std::string>("output.fileName") + ".dat", std::ofstream::trunc);
stats_file << "# Statistics of the solute transport equation generated by DORiE.\n";
stats_file << "#\n";
stats_file << "# run_time\t time\t mass\t boundary_flux\t";
stats_file << "mean_cw_x\t mean_cw_y\t mean_cw_z\t";
stats_file << "variance_cw_x\t variance_cw_y\t variance_cw_z\n";
stats_file.close();
}
// --- Operator Setup ---
......@@ -217,30 +234,118 @@ void TransportSimulation<Traits>::step()
template<typename Traits>
void TransportSimulation<Traits>::write_data () const
{
const auto time = controller->getTime();
auto solute_gf = std::make_shared<GFSolute>(gfs,u);
auto solute_t_gf = std::make_shared<GFTotalSolute>(solute_gf, water_content_gf);
// Compute statistical moments of Cw
using PDF = Dune::PDELab::ScalarScaledGridFunctionAdapter<GFSolute>;
PDF pdf(1./gf_integrator(*solute_gf,2*order+2),*solute_gf);
auto moments = statistical_moments(pdf,order*4,2);
auto moments_squared = moments[0];
for (int i = 0; i < dim; ++i) moments_squared[i]*=moments_squared[i];
const auto variance = moments[1] - moments_squared;
Dune::FieldVector<double,3> moments0_3d(0.), variance_3d(0.);
for (int i = 0; i < dim; ++i)
{
moments0_3d[i] = moments[0][i];
variance_3d[i] = variance[i];
}
RF mass = gf_integrator(*solute_t_gf,2*order+2);
RF boundary_flux = std::numeric_limits<RF>::quiet_NaN();
if constexpr (enable_rt_engine) {
if (enable_fluxrc) {
auto cwfluxr = get_solute_flux_reconstructed();
boundary_flux = boundary_gf_integrator(*cwfluxr,2*flux_order+2);
}
}
if (time_steps==0) cache_time = time;
if (time_steps==0) cache_mass = mass;
if (time_steps==0) cache_boundary_flux = boundary_flux;
// hard check on mass
if (std::isnan(mass)) {
this->_log->error("Mass of solute is NaN: PDE solver has diverged.");
DUNE_THROW(MathError,"PDE solver has diverged!");
}
// check relative change of mass
if (Dune::FloatCmp::ne(mass,0.))
{
// change in time
auto dt = time - cache_time;
// change of mass by PDE solution
auto mass_dt = (mass - cache_mass);
// change of mass by boundary conditions
mass_dt += 0.5*dt*(boundary_flux + cache_boundary_flux);
// relative mass change
auto rel_mass_dt = mass_dt/mass;
auto mass_tol = inifile.get<double>("conservationOfMass.checkTol");
auto check_mass_str = inifile.get<std::string>("conservationOfMass.check");
auto check_mass_lvl = Dorie::to_check_lvl(check_mass_str);
const bool passed = rel_mass_dt < mass_tol;
if (passed) {
this->_log->trace(" Conservation of mass does not hold for"
" {}% of mass",
rel_mass_dt);
}
else if (not passed && check_mass_lvl == CheckLevel::warn) {
this->_log->warn("Conservation of mass above tolerance:"
" {}% of mass is not conserved.",
rel_mass_dt);
}
else if (not passed && check_mass_lvl == CheckLevel::error) {
this->_log->error("Conservation of mass above tolerance:"
" {}% of mass is not conserved.",
rel_mass_dt);
DUNE_THROW(MathError, "Conservation of mass does not obey tolerance");
}
}
// update cache
cache_time = time;
cache_mass = mass;
cache_boundary_flux = boundary_flux;
// Configure statistics file
std::ofstream stats_file;
stats_file.open(inifile.get<std::string>("output.fileName") + ".dat", std::ofstream::app);
stats_file << timer.elapsed() << "\t" << time << "\t" << mass << "\t";
stats_file << boundary_flux << "\t" << moments0_3d << "\t" << variance_3d;
stats_file << std::endl;
stats_file.close();
if (vtkwriter)
{
if (inifile.get<bool>("output.vertexData")) {
vtkwriter->template addVertexData<GFSolute>(get_solute(),"solute");
vtkwriter->template addVertexData<GFTotalSolute>(get_total_solute(),"total_solute");
vtkwriter->template addVertexData<GFSolute>(solute_gf,"solute");
vtkwriter->template addVertexData<GFTotalSolute>(solute_t_gf,"total_solute");
if constexpr (enable_rt_engine)
if (enable_fluxrc) {
auto cwfluxr = get_solute_flux_reconstructed();
auto RT_name = "flux_RT" + std::to_string(flux_order);
vtkwriter->template addVertexData<GFFluxReconstruction>
(get_solute_flux_reconstructed(),RT_name);
vtkwriter->template addVertexData<GFFluxReconstruction>(cwfluxr,RT_name);
}
} else {
vtkwriter->template addCellData<GFSolute>(get_solute(),"solute");
vtkwriter->template addCellData<GFTotalSolute>(get_total_solute(),"total_solute");
vtkwriter->template addCellData<GFSolute>(solute_gf,"solute");
vtkwriter->template addCellData<GFTotalSolute>(solute_t_gf,"total_solute");
if constexpr (enable_rt_engine)
if (enable_fluxrc) {
auto cwfluxr = get_solute_flux_reconstructed();
auto RT_name = "flux_RT" + std::to_string(flux_order);
vtkwriter->template addCellData<GFFluxReconstruction>
(get_solute_flux_reconstructed(),RT_name);
vtkwriter->template addCellData<GFFluxReconstruction>(cwfluxr,RT_name);
}
}
try{
const auto time = controller->getTime();
this->_log->trace("Writing solution at time {:.2e}", time);
vtkwriter->write(time, output_type);
}
......
......@@ -16,6 +16,7 @@
#include <dune/dorie/common/simulation_base.hh>
#include <dune/dorie/common/util.hh>
#include <dune/dorie/common/utility.hh>
#include <dune/dorie/common/interpolator.hh>
#include <dune/dorie/common/time_controller.hh>
#include <dune/dorie/common/grid_creator.hh>
......@@ -296,12 +297,15 @@ protected:
std::unique_ptr<Writer> vtkwriter;
mutable std::shared_ptr<GFFluxReconstruction> cache_solute_flux_gf_rt;
unsigned int time_steps;
Dune::Timer timer;
const bool enable_fluxrc;
unsigned int time_steps;
mutable std::shared_ptr<GFFluxReconstruction> cache_solute_flux_gf_rt;
mutable RF cache_mass;
mutable RF cache_boundary_flux;
mutable TimeField cache_time;
public:
/*-----------------------------------------------------------------------*//**
......