Commit 4c8f8446 authored by Lukas Riedel's avatar Lukas Riedel 📝 Committed by Santiago Ospina De Los Ríos

Calculate velocity instead of flux in CFL condition

* Pass water content grid function to CFL condition and calculate
  actual velocity together with flux.
* Read courant number once and throw error if it is negative.
* Print CFL condition result into trace log.
parent f894a0ff
......@@ -71,6 +71,7 @@
* Transport model option `dirichletMode` was not working properly !121
* Use unsafe loader of PyYAML v5.2 for loading parameter scraper data !177
* CFL condition in explicit Transport model serves as time step upper limit !179
* Use apparent solute velocity instead of water flux in CFL condition !180
### Deprecated
......
......@@ -223,15 +223,14 @@ A full list is found at https://www.w3.org/TR/REC-html40/sgml/entities.html
</parameter>
<parameter name="courant">
<definition> Courant number for explicit methods. It is a scale for the
maximum stable time step according to the CFL condition. The CFL
condition penalizes high velocities, low dispersion coefficients, and
small grid elements. Hence, Courant numbers near to 1 tend to be more
instable while numbers near to 0 tend to considerably limit the
maximum time step.
<definition> Courant number for explicit methods. Serves as multiplier
of the CFL condition which estimates the maximum time step based on
effective velocity and grid element size. Reduce this value if
the solution is unstable. Notice that doing so also limits the maximum
time step and hence increases the overall simulation time.
</definition>
<values> 0 &lt; float &lt; 1 </values>
<suggestion> 0.8 </suggestion>
<suggestion> 0.5 </suggestion>
</parameter>
<parameter name="FEorder">
......
......@@ -39,14 +39,12 @@ namespace Dorie{
*
* @return The CFL-condition (\f$\mathcal{CFL}\f$).
*/
template<class GF, class RF = double>
RF cfl_condition(const GF& gf, unsigned int intorder = 1)
template<class VectorGF, class ScalarGF, class RF = double>
RF cfl_condition(const VectorGF& water_flux,
const ScalarGF& water_content,
const unsigned int intorder = 1)
{
using Traits = typename GF::Traits;
using DomainField = typename Traits::DomainFieldType;
using RangeField = typename Traits::RangeFieldType;
const auto& gv = gf.getGridView();
const auto& gv = water_flux.getGridView();
RF suggested_dt = std::numeric_limits<RF>::max();
......@@ -57,30 +55,35 @@ namespace Dorie{
// Check that cell is not a vertex
assert(geo.corners()>1);
DomainField h_T = 0.;
// Cell diameter: maximum distance between all vertices
using DF = typename VectorGF::Traits::DomainFieldType;
DF h_T = std::numeric_limits<DF>::max();
// Compute shortest axis of cell
for (int i = 1; i < geo.corners(); ++i)
for (int j = 0; j < i-1; ++j)
h_T = std::max(h_T,(geo.corner(i)-geo.corner(j)).two_norm());
h_T = std::min(h_T, (geo.corner(i)-geo.corner(j)).two_norm());
// Check that h_T is not zero in debug mode, who knows.
assert(Dune::FloatCmp::gt(h_T,0.));
RangeField max_water_flux = 0.;
typename VectorGF::Traits::RangeFieldType max_velocity = 0.;
// Find the maximum velocity in the cell
for (const auto& point : quadratureRule(geo, intorder))
{
typename Traits::RangeType water_flux;
typename VectorGF::Traits::RangeType velocity;
typename ScalarGF::Traits::RangeType scale;
const auto& x = point.position();
gf.evaluate(cell,x,water_flux);
max_water_flux = std::max(max_water_flux,water_flux.two_norm());
water_flux.evaluate(cell, x, velocity);
water_content.evaluate(cell, x, scale);
velocity /= scale;
max_velocity = std::max(max_velocity, velocity.two_norm());
}
// Calculate the cfl for this cell
if (Dune::FloatCmp::gt(max_water_flux,0.))
suggested_dt = std::min(suggested_dt,h_T/max_water_flux);
if (Dune::FloatCmp::gt(max_velocity, 0.))
suggested_dt = std::min(suggested_dt, h_T/max_velocity);
}
// Communicate the lowest cfl number to all procesors
......@@ -90,4 +93,4 @@ namespace Dorie{
}
}
#endif
\ No newline at end of file
#endif
......@@ -21,6 +21,7 @@ ModelTransport<Traits>::ModelTransport(
, mbe_slop(estimate_mbe_entries<typename MBE::size_type>(Traits::dim,Traits::GridGeometryType))
, mbe_tlop(1)
, time_steps(0)
, _courant_number(inifile.get<RF>("numerics.courant"))
, enable_fluxrc(_inifile.get<bool>("fluxReconstruction.enable"))
{
// --- Output Policy ---
......@@ -76,6 +77,18 @@ ModelTransport<Traits>::ModelTransport(
else
DUNE_THROW(Dune::NotImplemented,"Time method not supported!");
if (not ts_param->implicit())
{
if (FloatCmp::le(_courant_number, 0.0)) {
this->_log->error("Courant number must be larger than zero");
DUNE_THROW(IOError, "Invalid courant number");
}
else if (FloatCmp::gt(_courant_number, 1.0)) {
this->_log->warn("Courant number larger than 1.0! Solver might not "
"converge, leading to instabilities in the solution.");
}
}
// --- Writer Setup --- //
if (output_policy() != OutputPolicy::None)
{
......@@ -159,12 +172,14 @@ void ModelTransport<Traits>::step()
}
// Suggest times step for explicit methods.
// This assumes that timestep is never increased on failure!
if (not ts_param->implicit()) {
double courant = inifile.get<double>("numerics.courant");
assert(courant>0.);
const auto cfl = Dune::Dorie::cfl_condition(*water_flux_gf);
time_step_ctrl.suggest_time_step(courant * cfl);
// NOTE: Do not futher increase time step on failure!
if (not ts_param->implicit())
{
const auto cfl = Dune::Dorie::cfl_condition(*water_flux_gf,
*water_content_gf,
order);
this->_log->trace("CFL condition: {}s", cfl);
time_step_ctrl.suggest_time_step(_courant_number * cfl);
}
// get time step suggestion from boundary conditions
......
......@@ -325,6 +325,8 @@ protected:
unsigned int time_steps;
const RF _courant_number;
const bool enable_fluxrc;
public:
......
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