[Richards&Transport] Add conservation of mass check optional in inifile

parent ffbfeb84
......@@ -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
......
......@@ -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
......
......@@ -327,21 +327,39 @@ void RichardsSimulation<Traits>::write_data () const
// check relative change of mass
if (Dune::FloatCmp::ne(mass,0.))
{
// change in time
auto dt = time - cache_time;
// change in time
auto dt = time - cache_time;
// change of mass by PDE solution
auto mass_dt = (mass - cache_mass);
// 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);
// 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;
if (Dune::FloatCmp::gt(rel_mass_dt,1E-5)) // TODO: add option on inifile
this->_log->warn("Conservation of mass does not hold:"
" {}% of mass is not conserved at this time-step.",
rel_mass_dt*100);
// 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
......
......@@ -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>
......
......@@ -3,6 +3,7 @@
#include <dune/dorie/common/grid_function_integrator.hh>
#include <dune/dorie/common/statistical_moments.hh>
#include <limits>
#include <cmath>
......@@ -285,10 +286,28 @@ void TransportSimulation<Traits>::write_data () const
// relative mass change
auto rel_mass_dt = mass_dt/mass;
if (Dune::FloatCmp::gt(rel_mass_dt,1E-5)) // TODO: add option on inifile
this->_log->warn("Conservation of mass does not hold:"
" {}% of mass is not conserved at this time-step.",
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
......
......@@ -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>
......
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