Make kinematic_wave more robust

source/framework/algorithm/include/lue/framework/algorithm/definition/kinematic_wave.hpp

@@ -9,6 +9,82 @@
 #include <cmath>
 
 
+// /* Using Newton-Raphson Method
+// */
+// typedef long double REAL;
+// REAL Qk1; /* Q at loop k+1 for i+1, j+1 */
+// REAL ab_pQ, deltaTX, C;
+// int count;
+//
+// REAL Qkx;
+// REAL fQkx;
+// REAL dfQkx;
+// POSTCOND(sizeof(REAL) >= 8);
+//
+// /* if no input then output = 0 */
+// if ((Qin+Qold+q) == 0) /* +q CW NEW! */
+// return(0);
+//
+// /* common terms */
+// ab_pQ = alpha*beta*pow(((Qold+Qin)/2),beta-1);
+// deltaTX = deltaT/deltaX;
+// C = deltaTX*Qin + alpha*pow(Qold,beta) + deltaT*q;
+//
+// /* 1. Initial guess Qk1. */
+// /* 2. Evaluate function f at Qkx. */
+// /* 3. Evaluate derivative df at Qkx. */
+// /* 4. Check convergence. */
+//
+// /*
+// * There's a problem with the first guess of Qkx. fQkx is only defined
+// * for Qkx's > 0. Sometimes the first guess results in a Qkx+1 which is
+// * negative or 0. In that case we change Qkx+1 to 1e-30. This keeps the
+// * convergence loop healthy.
+// */
+//
+// Qkx = (deltaTX * Qin + Qold * ab_pQ + deltaT * q) / (deltaTX + ab_pQ);
+// Qkx = MAX(Qkx, 1e-30); /* added test-case calc::KinematicTest::iterate1 */
+// fQkx = deltaTX * Qkx + alpha * pow(Qkx, beta) - C; /* Current k */
+// dfQkx = deltaTX + alpha * beta * pow(Qkx, beta - 1); /* Current k */
+//
+//
+//
+//
+//
+// Qkx -= fQkx / dfQkx; /* Next k */
+// Qkx = MAX(Qkx, 1e-30);
+// count = 0;
+// do {
+// fQkx = deltaTX * Qkx + alpha * pow(Qkx, beta) - C; /* Current k */
+// dfQkx = deltaTX + alpha * beta * pow(Qkx, beta - 1); /* Current k */
+// Qkx -= fQkx / dfQkx; /* Next k */
+// Qkx = MAX(Qkx, 1e-30);
+// count++;
+// } while(fabsl(fQkx) > epsilon && count < MAX_ITERS);
+//
+// #ifdef DEBUG_DEVELOP
+// /* Our loop should converge in around 2 steps, otherwise something's
+// * wrong.
+// */
+// /* test-case calc::KinematicTest::iterate2
+// * is such a case, but values are very low
+// * 1e-30 is returned
+// */
+// if (0 && count == MAX_ITERS) {
+// printf("\nfQkx %g Qkx %g\n",(double)fQkx, (double)Qkx);
+// printf("Qin %g \n", Qin);
+// printf("Qold %g \n", Qold);
+// printf("q %g \n", q);
+// printf("alpha %g \n",alpha);
+// printf("beta %g \n", beta);
+// printf("deltaT %g \n", deltaT);
+// printf("deltaX %g \n", deltaX);
+// }
+// #endif
+// Qk1 = Qkx;
+// return(MAX(Qk1,0));
+
+
 namespace lue {
  namespace detail {
 
@@ -21,53 +97,67 @@ namespace lue {
  NonLinearKinematicWave(
  Float const upstream_discharge,
  Float const current_discharge,
- Float const inflow,
+ Float const lateral_inflow,
  Float const alpha,
  Float const beta,
  Float const time_step_duration,
  Float const channel_length):
 
  _upstream_discharge{upstream_discharge},
  _current_discharge{current_discharge},
- _inflow{inflow},
+ _lateral_inflow{lateral_inflow},
  _alpha{alpha},
  _beta{beta},
  _time_step_duration{time_step_duration}
 
  {
+ lue_hpx_assert(_upstream_discharge >= Float{0});
+ lue_hpx_assert(_current_discharge >= Float{0});
+ lue_hpx_assert(_upstream_discharge + _current_discharge + _lateral_inflow > Float{0});
+ lue_hpx_assert(_alpha > Float{0});
+ lue_hpx_assert(_beta > Float{0}); // TODO > 1?
+ lue_hpx_assert(_time_step_duration > Float{0});
+
  // Known terms, independent of new discharge
  _time_step_duration_over_channel_length = _time_step_duration / channel_length;
 
  lue_hpx_assert(_time_step_duration_over_channel_length > Float{0});
 
  _known_terms = _time_step_duration_over_channel_length * _upstream_discharge +
  _alpha * std::pow(_current_discharge, _beta) +
- _time_step_duration * _inflow;
+ _time_step_duration * _lateral_inflow;
  }
 
 
  /*!
  @brief Return a valid initial guess for the new discharge
+
+ Note that fq is only defined for discharges larger than zero. In case the initial guess
+ ends up being zero or negative, a small positive value is returned.
+
  */
  Float guess() const
  {
  // Small, but not zero!
  Float discharge_guess{1e-30};
 
  // pow(0, -) is not defined
- if (_current_discharge + _upstream_discharge > Float{0})
+ if ((_current_discharge + _upstream_discharge != Float{0}) || _beta >= Float{1})
  {
  Float const a_b_pq =
  _alpha * _beta *
  std::pow((_current_discharge + _upstream_discharge) / Float{2}, _beta - Float{1});
 
  lue_hpx_assert(!std::isnan(a_b_pq));
 
- discharge_guess = (_time_step_duration_over_channel_length * _upstream_discharge +
- a_b_pq * _current_discharge + _time_step_duration * _inflow) /
- (_time_step_duration_over_channel_length + a_b_pq);
+ discharge_guess =
+ (_time_step_duration_over_channel_length * _upstream_discharge +
+ a_b_pq * _current_discharge + _time_step_duration * _lateral_inflow) /
+ (_time_step_duration_over_channel_length + a_b_pq);
 
  lue_hpx_assert(!std::isnan(discharge_guess));
+
+ discharge_guess = std::max<Float>(discharge_guess, 1e-30);
  }
 
  lue_hpx_assert(discharge_guess > Float{0});
@@ -84,6 +174,8 @@ namespace lue {
 
  Float fq(Float const new_discharge) const
  {
+ lue_hpx_assert(new_discharge > Float{0}); // pow(0, -) is not defined
+
  return _time_step_duration_over_channel_length * new_discharge +
  _alpha * std::pow(new_discharge, _beta) - _known_terms;
  }
@@ -100,14 +192,18 @@ namespace lue {
 
  private:
 
+ //! Updated / new discharge in the upstream cell
  Float _upstream_discharge;
 
+ //! Current / previous discharge in the current cell
  Float _current_discharge;
 
- Float _inflow;
+ //! Lateral inflow
+ Float _lateral_inflow;
 
  Float _alpha;
 
+ //! Momentum coefficient / Boussinesq coefficient [1.01, 1.33] (Chow, p278)
  Float _beta;
 
  Float _time_step_duration;
@@ -122,13 +218,13 @@ namespace lue {
  Float iterate_to_new_discharge(
  Float const upstream_discharge, // Summed discharge for cells draining into current cell
  Float const current_discharge,
- Float const inflow,
+ Float const lateral_inflow,
  Float const alpha,
  Float const beta,
  Float const time_step_duration,
  Float const channel_length)
  {
- if (upstream_discharge + current_discharge + inflow == 0)
+ if (upstream_discharge + current_discharge + lateral_inflow <= 0)
  {
  // It's a dry place. No need to do anything fancy.
  return Float{0};
@@ -137,7 +233,7 @@ namespace lue {
  NonLinearKinematicWave<Float> kinematic_wave{
  upstream_discharge,
  current_discharge,
- inflow,
+ lateral_inflow,
  alpha,
  beta,
  time_step_duration,
@@ -146,19 +242,21 @@ namespace lue {
  Float const discharge_guess{kinematic_wave.guess()};
 
  // These brackets are crucial for obtaining a good performance
+ // Float const min_discharge{0};
+ // Float const max_discharge{std::numeric_limits<Float>::max()};
  Float const min_discharge{discharge_guess < Float{1} ? Float{0} : discharge_guess / Float{1000}};
  Float const max_discharge{
  discharge_guess < Float{1} ? Float{1000} : Float{1000} * discharge_guess};
 
  int const digits = static_cast<int>(std::numeric_limits<Float>::digits * 0.6);
 
  std::uintmax_t const max_nr_iterations{20};
- std::uintmax_t nr_iterations{max_nr_iterations};
+ std::uintmax_t actual_nr_iterations{max_nr_iterations};
 
  Float new_discharge = boost::math::tools::newton_raphson_iterate(
- kinematic_wave, discharge_guess, min_discharge, max_discharge, digits, nr_iterations);
+ kinematic_wave, discharge_guess, min_discharge, max_discharge, digits, actual_nr_iterations);
 
- if (nr_iterations == max_nr_iterations)
+ if (actual_nr_iterations == max_nr_iterations)
  {
  throw std::runtime_error(fmt::format(
  "Iterating to a solution took more iterations than expected (initial guess: {}, "
@@ -168,7 +266,7 @@ namespace lue {
  max_discharge));
  }
 
- lue_hpx_assert(nr_iterations < max_nr_iterations);
+ lue_hpx_assert(actual_nr_iterations < max_nr_iterations);
  lue_hpx_assert(new_discharge >= Float{0});
 
  return new_discharge;

source/framework/algorithm/test/flow_accumulation.cpp

@@ -126,17 +126,21 @@ namespace lue::test {
  {
  using Shape = ShapeT<FlowDirectionArray>;
 
- Shape array_shape{{5, 5}};
- Shape partition_shape{{5, 5}};
+ Shape const array_shape{{5, 5}};
+ Shape const partition_shape{{5, 5}};
 
  return create_partitioned_array<FlowDirectionArray>(
  array_shape,
  partition_shape,
- {
- {
- s, s, s, sw, sw, s, s, sw, sw, sw, se, s, sw, w, sw, se, s, sw, w, w, e, p, w, w, w,
- },
- });
+ // clang-format off
+ {{
+ s, s, s, sw, sw,
+ s, s, sw, sw, sw,
+ se, s, sw, w, sw,
+ se, s, sw, w, w,
+ e, p, w, w, w,
+ }} // clang-format on
+ );
  }