Commit 52fa9ee5 authored by Santiago Ospina De Los Ríos's avatar Santiago Ospina De Los Ríos

Merge branch...

Merge branch '171-negative-solute-concentrations-when-pulses-pass-through-soil-layer-interfaces' into 'master'

Resolve "Negative solute concentrations when pulses pass through soil layer interfaces"

Closes #171

See merge request !181
parents 35e9cd15 37da8371
......@@ -34,6 +34,7 @@
* CMake option to enable code coverage flags on all targets !172
* Steady state initial condition in Richards model !176
* Changes to config file parameters listed per version in user docs !175
* Control negative transport solution by a check policy !181
### Changed
* Data structures for storing and accessing parameter information !55
......
......@@ -268,6 +268,38 @@ A full list is found at https://www.w3.org/TR/REC-html40/sgml/entities.html
</parameter>
</category>
<category name="solutionCheck">
<parameter name="negativeConcentration">
<definition> Check for negative solute concentrations in the computed
solution.
* ``pass``: Passes silently when negative values are encountered.
* ``warn``: Issues a warning when finding negative values, but
continues the computation.
* ``retry``: Discard the computed solution, and repeat it with a
reduced time step size according to the ``transport.time`` settings.
</definition>
<suggestion> retry </suggestion>
<values> pass, warn, retry </values>
<comment> Choose handling: pass, warn, retry </comment>
<versionadded version="unreleased"> </versionadded>
</parameter>
<parameter name="negativeTol">
<definition> Relative tolerance when checking negative values. This
factor is multiplied with the *largest* value of the current solute
concentration.
</definition>
<values> float &gt; 0 </values>
<suggestion> 1E-6 </suggestion>
<versionadded version="unreleased"> </versionadded>
</parameter>
</category>
<category name="solverParameters">
<parameter name="reduction">
<definition> Required relative defect reduction.
......
......@@ -37,8 +37,23 @@ private:
public:
static_assert(ModelRichardsTraits::enable_rt_engine,
"Flux reconstruction is not available in the Richards model!");
using GFWaterContentContainer = GridFunctionContainer<const typename ModelRichardsTraits::GFWaterContent>;
/// Container for the water content grid function
/** \note Interpolates linearly between time steps.
* \see GridFunctionContainerPolicy
*/
using GFWaterContentContainer = GridFunctionContainer<
const typename ModelRichardsTraits::GFWaterContent,
double,
GridFunctionContainerPolicy::LinearInterpolation
>;
/// Container for the water flux grid function
/** \note Returns the next future value
* \see GridFunctionContainerPolicy
*/
using GFWaterFluxContainer = GridFunctionContainer<const typename ModelRichardsTraits::GFFluxReconstruction>;
private:
using ModelTransportTraits = Dune::Dorie::ModelTransportTraits<
BaseTraits,
......
......@@ -39,6 +39,24 @@ ModelTransport<Traits>::ModelTransport(
DUNE_THROW(NotImplemented,"Invalid output policy: " << output_policy_str);
}
// --- Solution Check Policy --- //
const auto check_negative_policy_str
= inifile.get<std::string>("solutionCheck.negativeConcentration");
if (check_negative_policy_str == "pass") {
check_negative_policy = CheckNegativePolicy::Pass;
}
else if (check_negative_policy_str == "warn") {
check_negative_policy = CheckNegativePolicy::Warn;
}
else if (check_negative_policy_str == "retry") {
check_negative_policy = CheckNegativePolicy::Retry;
}
else {
this->_log->error("Invalid policy for negative solute concentrations: {}",
check_negative_policy_str);
DUNE_THROW(IOError, "Invalid policy for negative concentrations");
}
// initialize time
this->_time = time_step_ctrl.get_time_begin();
......@@ -185,54 +203,53 @@ void ModelTransport<Traits>::step()
// get time step suggestion from boundary conditions
sboundary->set_time(this->_time);
const auto dt_suggestion = sboundary->max_time_step(this->_time);
time_step_ctrl.suggest_time_step(dt_suggestion);
// Allocate DOF vector for next time step
auto unext = std::make_shared<U>(*u);
// Obtain time variables
const auto time = this->_time;
auto dt = time_step_ctrl.get_time_step(time);
// Solver loop
// NOTE: Runs infinitely until solution is computed and valid (causes break)
// or exception is thrown.
while (true)
{
// Obtain time variables
time_step_ctrl.suggest_time_step(dt_suggestion);
const auto time = this->_time;
const auto dt = time_step_ctrl.get_time_step(time);
try {
// long time step log for level 'debug'
if (this->_log->should_log(spdlog::level::debug)) {
this->_log->debug(" Time {:.2e} + {:.2e} -> {:.2e}",
time, dt, time+dt);
}
std::shared_ptr<U> unext = std::make_shared<U>(*u);
// create a DOF vector for next step
// Long time step log for level 'debug'
if (this->_log->should_log(spdlog::level::debug)) {
this->_log->debug(" Time {:.2e} + {:.2e} -> {:.2e}",
time, dt, time+dt);
}
try
{
// Apply the solver
dwarn.push(false);
osm->apply(time, dt, *u, *unext);
dwarn.pop();
// short time step log for level 'info' (after success)
if (not this->_log->should_log(spdlog::level::debug)) {
this->_log->info("Time Step {}: {:.2e} + {:.2e} -> {:.2e}",
time_steps,
time,
dt,
time+dt);
// Check the solution
if (auto [valid, msg] = check_solution(*u, *unext);
valid)
{
// Solution computed and valid: leave loop
break;
}
else
{
this->_log->debug(" Invalid solution");
this->_log->trace(" Error during solution check: {}", msg);
}
u = unext;
// update time and time step
this->_time += dt;
time_step_ctrl.increase_time_step();
// leave loop
break;
}
catch (Dune::ISTLError &e){
this->_log->debug(" Solver not converged");
this->_log->trace(" Error in linear solver: {}", e.what());
}
catch (Dune::Exception &e){
this->_log->debug(" Solver not converged");
this->_log->error(" Unexpected error in solver: {}", e.what());
DUNE_THROW(Dune::Exception, "Critical error in linear solver");
DUNE_THROW(Dune::Exception, "Critical error in solver");
}
// rethrow anything else
catch (...) {
......@@ -240,18 +257,36 @@ void ModelTransport<Traits>::step()
}
// abort if minimal time step is already reached
if (Dune::FloatCmp::eq(dt, time_step_ctrl.get_time_step_min())) {
if (Dune::FloatCmp::le(dt, time_step_ctrl.get_time_step_min())) {
this->_log->error("Computing a solution for the minimal time "
"step failed");
"step failed");
DUNE_THROW(Exception, "No solution for minimal time step");
}
// Reduce time step and try again
time_step_ctrl.decrease_time_step();
dt = time_step_ctrl.get_time_step(time);
} // while (true)
cache_solute_flux_gf_rt.reset();
// Short time step log for level 'info' (after success)
if (not this->_log->should_log(spdlog::level::debug)) {
this->_log->info("Time Step {}: {:.2e} + {:.2e} -> {:.2e}",
time_steps,
time,
dt,
time+dt);
}
// Overwrite last solution with new solution
u = unext;
// Update time and time step
this->_time += dt;
time_step_ctrl.increase_time_step();
time_steps++;
cache_solute_flux_gf_rt.reset();
if (this->output_policy() == OutputPolicy::EndOfTransportStep)
write_data();
}
......@@ -311,5 +346,48 @@ void ModelTransport<Traits>::write_data () const
}
}
template<typename Traits>
std::pair<bool, std::string> ModelTransport<Traits>::check_solution(
const U& u_before,
const U& u_after
) const
{
// Check for negative values
if (check_negative_policy != CheckNegativePolicy::Pass)
{
// Compute tolerance from largest solute concentration and user input
const auto max = *std::max_element(u_after.begin(), u_after.end());
const auto eps = max * inifile.get<RF>("solutionCheck.negativeTol");
// Check if values are negative and record minimum value encountered
auto min = std::numeric_limits<RF>::max();
auto check_negative = [eps, &min](const auto val)
{
min = std::min(min, val);
return val < -eps;
};
// Perform the check
if (std::any_of(u_after.begin(), u_after.end(), check_negative))
{
// NOTE: fmt library included via spdlog
const std::string msg = fmt::format("Negative value in solution: {:.2e}",
min);
// Issue warning, but pass check
if (check_negative_policy == CheckNegativePolicy::Warn) {
this->_log->warn(msg);
}
// Do not pass check
else {
return std::make_pair(false, msg);
}
}
}
// Default: Pass without message
return std::make_pair(true, "");
}
} // namespace Dorie
} // namespace Dune
......@@ -5,9 +5,12 @@
#include <string>
#include <iostream>
#include <limits>
#include <algorithm>
#include <utility>
#include <dune/common/parametertree.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/float_cmp.hh>
#include <dune/pdelab/backend/istl.hh>
#include <dune/pdelab/gridoperator/gridoperator.hh>
......@@ -24,6 +27,7 @@
#include <dune/dorie/common/cfl_condition.hh>
#include <dune/dorie/common/flux_reconstruction/raviart_thomas.hh>
#include <dune/dorie/common/boundary_condition/manager.hh>
#include <dune/dorie/common/logging.hh>
#include <dune/dorie/model/base.hh>
#include <dune/dorie/model/transport/adapters/total_solute.hh>
......@@ -175,6 +179,14 @@ struct ModelTransportTraits : public BaseTraits
// -- Utility Class Definitions -- //
/// VTK Output writer base class
using Writer = Dune::Dorie::GridFunctionVTKSequenceWriter<GV>;
/// Policy when checking for negative values in the solution
enum class CheckNegativePolicy
{
Pass, //!< Always pass the check
Warn, //!< Pass the check but issue warnings for negative values
Retry //!< Reject the solution for negative values
};
};
......@@ -283,6 +295,9 @@ protected:
/// Time stepping scheme
using TimeSteppingParameter = Dune::PDELab::TimeSteppingParameterInterface<RF>;
/// Policy for negative values in solution
using CheckNegativePolicy = typename Traits::CheckNegativePolicy;
// -- Member variables -- //
Dune::MPIHelper& helper;
......@@ -329,6 +344,9 @@ protected:
const bool enable_fluxrc;
/// Stored policy for negative values in solution
CheckNegativePolicy check_negative_policy;
public:
/*-----------------------------------------------------------------------*//**
......@@ -558,6 +576,23 @@ protected:
*/
void operator_setup ();
/// Check the solution of a time step
/** This function currently performs the following checks:
*
* * No negative values in the resulting solution \p u_after
*
* The method returns a pair of `bool` and `string`. The boolean indicates
* if the check was successful, and the string contains a check error
* message if it was not.
*
* \param u_before Solution at the start of the time step
* \param u_after Computed solution at the end of the time step
* \return Pair of `bool` (check passed) and `string` (check fail
* message, if any)
*/
std::pair<bool, std::string> check_solution (const U& u_before,
const U& u_after) const;
private:
#ifdef GTEST
......
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