Repository: computationalgeography/lue (language: python)

Issue:
Add support for negative inflow to `kinematic_wave`
Lateral inflow represents the addition of water (e.g. precepitation), or removal of water (e.g. infiltration, pumping). Very comparable with `accu_threshold`. The outflow passed in (negative inflow) is a maximum / potential value and may be different from what is actually removed from the cell.

We could add a second return value containing the actual amount of water removed, given the actual amount of water available in each cell. Until then, it is impossible for users to check the water budget when less water is removed than potentially can be removed. See #527 

Lateral inflow can be negative. Make sure the actual outflow (negative inflow) is not larger than the actual amount of water in the cell.


Discussion / hints:

Fix patch:
diff --git a/source/framework/algorithm/include/lue/framework/algorithm/definition/kinematic_wave.hpp b/source/framework/algorithm/include/lue/framework/algorithm/definition/kinematic_wave.hpp
index c287c5a08..2c9677c5d 100644
--- a/source/framework/algorithm/include/lue/framework/algorithm/definition/kinematic_wave.hpp
+++ b/source/framework/algorithm/include/lue/framework/algorithm/definition/kinematic_wave.hpp
@@ -4,6 +4,8 @@
 #include "lue/framework/algorithm/kinematic_wave.hpp"
 #include "lue/framework/algorithm/routing_operation_export.hpp"
 #include "lue/macro.hpp"
+#include <limits>
+// #define BOOST_MATH_INSTRUMENT
 #include <boost/math/tools/roots.hpp>
 #include <fmt/format.h>
 #include <cmath>
@@ -21,7 +23,7 @@ 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,
@@ -29,45 +31,59 @@ namespace lue {
 
                     _upstream_discharge{upstream_discharge},
                     _current_discharge{current_discharge},
-                    _inflow{inflow},
+                    _lateral_inflow{lateral_inflow},
                     _alpha{alpha},
                     _beta{beta},
+                    _alpha_beta{_alpha * _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});
+                    lue_hpx_assert(_time_step_duration > Float{0});
+                    lue_hpx_assert(channel_length > 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
+                auto guess() const -> Float
                 {
                     // Small, but not zero!
-                    Float discharge_guess{1e-30};
+                    static Float const min_discharge{std::numeric_limits<Float>::min()};
+                    Float discharge_guess{min_discharge};
 
                     // 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 *
+                            _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, min_discharge);
                     }
 
                     lue_hpx_assert(discharge_guess > Float{0});
@@ -76,40 +92,48 @@ namespace lue {
                 }
 
 
-                std::pair<Float, Float> operator()(Float const new_discharge) const
+                auto operator()(Float const new_discharge) const -> std::pair<Float, Float>
                 {
                     return std::make_pair(fq(new_discharge), dfq(new_discharge));
                 }
 
 
-                Float fq(Float const new_discharge) const
+                auto fq(Float const new_discharge) const -> Float
                 {
+                    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;
                 }
 
 
-                Float dfq(Float const new_discharge) const
+                auto dfq(Float const new_discharge) const -> Float
                 {
                     lue_hpx_assert(new_discharge > Float{0});  // pow(0, -) is not defined
 
                     return _time_step_duration_over_channel_length +
-                           _alpha * _beta * std::pow(new_discharge, _beta - Float{1});
+                           _alpha_beta * std::pow(new_discharge, _beta - Float{1});
                 }
 
 
             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 _alpha_beta;
+
                 Float _time_step_duration;
 
                 Float _time_step_duration_over_channel_length;
@@ -122,53 +146,107 @@ 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)
+            lue_hpx_assert(upstream_discharge >= Float{0});
+            lue_hpx_assert(current_discharge >= Float{0});
+            lue_hpx_assert(alpha >= Float{0});
+            lue_hpx_assert(beta >= Float{0});
+            lue_hpx_assert(time_step_duration > Float{0});
+            lue_hpx_assert(channel_length > Float{0});
+
+            // Lateral inflow can represent two things:
+            // - Actual inflow from an external source (positive value): e.g.: precepitation
+            // - Potential extraction to an internal sink (negative value): e.g.: infiltration, pumping
+            //
+            // Inflow:
+            // Add the amount of water to the discharge computed
+            //
+            // Extraction:
+            // Subtract as much water from the discharge computed as possible. To allow for water balance
+            // checks, we should output the actual amount of water that got extracted from the cell. This is
+            // the difference between the discharge computed and the potential extraction passed in.
+            // https://github.com/computationalgeography/lue/issues/527
+
+            Float new_discharge{0};
+
+            if (upstream_discharge + current_discharge > 0 || lateral_inflow > 0)
             {
-                // It's a dry place. No need to do anything fancy.
-                return Float{0};
+                // The cell receives water, from upstream and/or from an external source
+                Float const inflow = lateral_inflow >= 0 ? lateral_inflow : Float{0};
+
+                NonLinearKinematicWave<Float> kinematic_wave{
+                    upstream_discharge,
+                    current_discharge,
+                    inflow,
+                    alpha,
+                    beta,
+                    time_step_duration,
+                    channel_length};
+
+                Float const discharge_guess{kinematic_wave.guess()};
+
+                // "If the function is 'Really Well Behaved' (is monotonic and has only one root) the bracket
+                // bounds min and max may as well be set to the widest limits"
+                Float const min_discharge{0};
+                Float const max_discharge{std::numeric_limits<Float>::max()};
+                int const digits = static_cast<int>(std::numeric_limits<Float>::digits * 0.6);
+
+                // In general, 2-3 iterations are enough. In rare cases more are needed. The unit tests don't
+                // seem to reach 8, so max 10 should be enough. This max is used in special cases:
+                // - upstream_discharge + current_discharge == 0 and beta <  1
+                std::uintmax_t const max_nr_iterations{10};
+                std::uintmax_t actual_nr_iterations{max_nr_iterations};
+
+                // https://www.boost.org/doc/libs/1_85_0/libs/math/doc/html/math_toolkit/roots_deriv.html
+                // std::cout.precision(std::numeric_limits<Float>::digits10);
+                new_discharge = boost::math::tools::newton_raphson_iterate(
+                    kinematic_wave,
+                    discharge_guess,
+                    min_discharge,
+                    max_discharge,
+                    digits,
+                    actual_nr_iterations);
+
+                // TODO We can't throw an exception as this actually happens once in a while
+                // https://github.com/computationalgeography/lue/issues/703
+                // if (actual_nr_iterations == max_nr_iterations)
+                // {
+                //     // This only seems to happen when upstream_discharge == 0, current_discharge == 0, and
+                //     // lateral_inflow > 0
+                //     throw std::runtime_error(fmt::format(
+                //         "Iterating to a solution took more iterations than expected: "
+                //         "    upstream discharge: {}, "
+                //         "    current discharge: {}, "
+                //         "    lateral inflow: {}, "
+                //         "    alpha: {}, "
+                //         "    beta: {}, "
+                //         "    time step duration: {}, "
+                //         "    channel length: {}, "
+                //         "    initial guess: {}",
+                //         upstream_discharge,
+                //         current_discharge,
+                //         inflow,
+                //         alpha,
+                //         beta,
+                //         time_step_duration,
+                //         channel_length,
+                //         discharge_guess));
+                // }
             }
 
-            NonLinearKinematicWave<Float> kinematic_wave{
-                upstream_discharge,
-                current_discharge,
-                inflow,
-                alpha,
-                beta,
-                time_step_duration,
-                channel_length};
-
-            Float const discharge_guess{kinematic_wave.guess()};
-
-            // These brackets are crucial for obtaining a good performance
-            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};
-
-            Float new_discharge = boost::math::tools::newton_raphson_iterate(
-                kinematic_wave, discharge_guess, min_discharge, max_discharge, digits, nr_iterations);
-
-            if (nr_iterations == max_nr_iterations)
+            if (lateral_inflow < Float{0})
             {
-                throw std::runtime_error(fmt::format(
-                    "Iterating to a solution took more iterations than expected (initial guess: {}, "
-                    "brackets: [{}, {}])",
-                    discharge_guess,
-                    min_discharge,
-                    max_discharge));
+                // Convert units: m³ / m / s → m³ / s
+                Float const extraction{std::min(channel_length * std::abs(lateral_inflow), new_discharge)};
+
+                new_discharge -= extraction;
             }
 
-            lue_hpx_assert(nr_iterations < max_nr_iterations);
             lue_hpx_assert(new_discharge >= Float{0});
 
             return new_discharge;


Test patch:
diff --git a/source/framework/algorithm/test/clump_test.cpp b/source/framework/algorithm/test/clump_test.cpp
index 18a1b8960..5dd4bdc75 100644
--- a/source/framework/algorithm/test/clump_test.cpp
+++ b/source/framework/algorithm/test/clump_test.cpp
@@ -1265,9 +1265,6 @@ BOOST_AUTO_TEST_CASE(random_input)
     std::random_device random_device{};
     std::default_random_engine random_number_engine(random_device());
 
-    // Shape const array_shape{10, 10};
-    // Shape const partition_shape{5, 5};
-
     Shape const array_shape = [&]()
     {
         lue::Count const min{100};
@@ -1333,9 +1330,6 @@ BOOST_AUTO_TEST_CASE(random_input)
 
     using namespace lue::value_policies;
 
-    // lue::wait_all(zone_array.partitions());
-    // hpx::cout << zone_array << std::endl;
-
     // Test whether clumping the different partitioned arrays result in the same number of clumps. If so,
     // merging clumps seems to have worked.
     BOOST_TEST_CONTEXT("nondiagonal")
diff --git a/source/framework/algorithm/test/flow_accumulation.cpp b/source/framework/algorithm/test/flow_accumulation.cpp
index ed10adf6b..8ff39272b 100644
--- a/source/framework/algorithm/test/flow_accumulation.cpp
+++ b/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
+        );
     }
 
 
diff --git a/source/framework/algorithm/test/kinematic_wave_test.cpp b/source/framework/algorithm/test/kinematic_wave_test.cpp
index 71c3dd4fc..5a00bf4c0 100644
--- a/source/framework/algorithm/test/kinematic_wave_test.cpp
+++ b/source/framework/algorithm/test/kinematic_wave_test.cpp
@@ -1,13 +1,19 @@
 #define BOOST_TEST_MODULE lue framework algorithm kinematic_wave
 #include "flow_accumulation.hpp"
+#include "lue/framework/algorithm/definition/kinematic_wave.hpp"
 #include "lue/framework/algorithm/value_policies/all.hpp"
 #include "lue/framework/algorithm/value_policies/greater_than.hpp"
 #include "lue/framework/algorithm/value_policies/kinematic_wave.hpp"
+#include "lue/framework/algorithm/value_policies/uniform.hpp"
 #include "lue/framework/test/hpx_unit_test.hpp"
 
 
+namespace tt = boost::test_tools;
+
+
 BOOST_AUTO_TEST_CASE(pcraster_manual_example)
 {
+    // https://pcraster.geo.uu.nl/pcraster/4.4.1/documentation/pcraster_manual/sphinx/op_kinematic.html
     auto flow_direction = lue::test::pcraster_example_flow_direction();
 
     auto const array_shape{flow_direction.shape()};
@@ -15,13 +21,18 @@ BOOST_AUTO_TEST_CASE(pcraster_manual_example)
 
     using Element = double;
 
-    auto const discharge = lue::test::create_partitioned_array<lue::PartitionedArray<Element, 2>>(
+    auto const current_discharge = lue::test::create_partitioned_array<lue::PartitionedArray<Element, 2>>(
         array_shape,
         partition_shape,
+        // clang-format off
         {{
-            10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10,
-            10, 10, 10, 10, 10, 10, 50, 50, 50, 50, 50, 49,
-        }});
+            10, 10, 10, 10, 10,
+            10, 10, 10, 10, 10,
+            10, 10, 10, 10, 10,
+            10, 10, 10, 10, 50,
+            50, 50, 50, 50, 49,
+        }}  // clang-format on
+    );
 
     auto const inflow = lue::create_partitioned_array<Element>(array_shape, partition_shape, 0);
 
@@ -34,16 +45,20 @@ BOOST_AUTO_TEST_CASE(pcraster_manual_example)
     auto const discharge_we_want = lue::test::create_partitioned_array<lue::PartitionedArray<Element, 2>>(
         array_shape,
         partition_shape,
+        // clang-format off
         {{
-            2.32293, 2.32293, 2.32293, 2.32293,  2.32293,  4.00491,  4.00491, 5.76591,  4.00491,
-            2.32293, 5.27401, 9.81514, 8.40418,  4.00491,  2.32293,  2.32293, 21.22170, 9.68365,
-            9.60760, 7.19019, 7.19019, 57.27763, 17.74742, 12.90772, 7.09124,
-        }});
-
-    auto const kinematic_wave = lue::value_policies::kinematic_wave(
-        flow_direction, discharge, inflow, alpha, beta, time_step_duration, channel_length);
-
-    lue::test::check_arrays_are_close(kinematic_wave, discharge_we_want);
+            2.32293,  2.32293,  2.32293,  2.32293, 2.32293,
+            4.00491,  4.00491,  5.76591,  4.00491, 2.32293,
+            5.27401,  9.81514,  8.40418,  4.00491, 2.32293,
+            2.32293, 21.22170,  9.68365,  9.60760, 7.19019,
+            7.19019, 57.27763, 17.74742, 12.90772, 7.09124,
+        }}  // clang-format on
+    );
+
+    auto const new_discharge = lue::value_policies::kinematic_wave(
+        flow_direction, current_discharge, inflow, alpha, beta, time_step_duration, channel_length);
+
+    lue::test::check_arrays_are_close(new_discharge, discharge_we_want);
 }
 
 
@@ -56,7 +71,7 @@ BOOST_AUTO_TEST_CASE(zero_discharge_and_inflow)
 
     using Element = double;
 
-    auto const discharge = lue::create_partitioned_array<Element>(array_shape, partition_shape, 0);
+    auto const current_discharge = lue::create_partitioned_array<Element>(array_shape, partition_shape, 0);
     auto const inflow = lue::create_partitioned_array<Element>(array_shape, partition_shape, 0);
 
     Element const alpha = 1.5;
@@ -67,10 +82,10 @@ BOOST_AUTO_TEST_CASE(zero_discharge_and_inflow)
 
     auto const discharge_we_want = lue::create_partitioned_array<Element>(array_shape, partition_shape, 0);
 
-    auto const kinematic_wave = lue::value_policies::kinematic_wave(
-        flow_direction, discharge, inflow, alpha, beta, time_step_duration, channel_length);
+    auto const new_discharge = lue::value_policies::kinematic_wave(
+        flow_direction, current_discharge, inflow, alpha, beta, time_step_duration, channel_length);
 
-    lue::test::check_arrays_are_close(kinematic_wave, discharge_we_want);
+    lue::test::check_arrays_are_close(new_discharge, discharge_we_want);
 }
 
 
@@ -83,7 +98,7 @@ BOOST_AUTO_TEST_CASE(non_zero_discharge_and_zero_inflow)
 
     using Element = double;
 
-    auto const discharge = lue::create_partitioned_array<Element>(array_shape, partition_shape, 1);
+    auto const current_discharge = lue::create_partitioned_array<Element>(array_shape, partition_shape, 1);
     auto const inflow = lue::create_partitioned_array<Element>(array_shape, partition_shape, 0);
 
     Element const alpha = 1.5;
@@ -92,12 +107,12 @@ BOOST_AUTO_TEST_CASE(non_zero_discharge_and_zero_inflow)
 
     auto const channel_length = lue::create_partitioned_array<Element>(array_shape, partition_shape, 10);
 
-    auto const kinematic_wave = lue::value_policies::kinematic_wave(
-        flow_direction, discharge, inflow, alpha, beta, time_step_duration, channel_length);
+    auto const new_discharge = lue::value_policies::kinematic_wave(
+        flow_direction, current_discharge, inflow, alpha, beta, time_step_duration, channel_length);
 
     using namespace lue::value_policies;
 
-    BOOST_CHECK(lue::value_policies::all(kinematic_wave > Element{0}).get());
+    BOOST_CHECK(lue::value_policies::all(new_discharge > Element{0}).get());
 }
 
 
@@ -110,7 +125,7 @@ BOOST_AUTO_TEST_CASE(zero_discharge_and_non_zero_inflow)
 
     using Element = double;
 
-    auto const discharge = lue::create_partitioned_array<Element>(array_shape, partition_shape, 0);
+    auto const current_discharge = lue::create_partitioned_array<Element>(array_shape, partition_shape, 0);
     auto const inflow = lue::create_partitioned_array<Element>(array_shape, partition_shape, 1);
 
     Element const alpha = 1.5;
@@ -119,14 +134,180 @@ BOOST_AUTO_TEST_CASE(zero_discharge_and_non_zero_inflow)
 
     auto const channel_length = lue::create_partitioned_array<Element>(array_shape, partition_shape, 10);
 
-    auto const kinematic_wave = lue::value_policies::kinematic_wave(
-        flow_direction, discharge, inflow, alpha, beta, time_step_duration, channel_length);
+    auto const new_discharge = lue::value_policies::kinematic_wave(
+        flow_direction, current_discharge, inflow, alpha, beta, time_step_duration, channel_length);
 
     using namespace lue::value_policies;
 
-    BOOST_CHECK(lue::value_policies::all(kinematic_wave > Element{0}).get());
+    BOOST_CHECK(lue::value_policies::all(new_discharge > Element{0}).get());
 }
 
 
 // TODO Add tests for extreme inputs:
 // - What about "wrong" values for alpha, beta, time_step_duration?
+
+
+BOOST_AUTO_TEST_CASE(dry_cell)
+{
+    {
+        double const new_discharge{lue::detail::iterate_to_new_discharge<double>(
+            0,     // upstream_discharge
+            0,     // current_discharge
+            0,     // lateral_inflow
+            1.5,   // alpha
+            0.6,   // beta
+            15,    // time_step_duration
+            10)};  // channel_length
+        double const discharge_we_want{0};
+
+        BOOST_TEST(new_discharge == discharge_we_want, tt::tolerance(1e-6));
+    }
+
+    {
+        double const new_discharge{lue::detail::iterate_to_new_discharge<double>(
+            1,     // upstream_discharge
+            0,     // current_discharge
+            -1,    // lateral_inflow
+            1.5,   // alpha
+            0.6,   // beta
+            15,    // time_step_duration
+            10)};  // channel_length
+        double const discharge_we_want{0};
+
+        BOOST_TEST(new_discharge == discharge_we_want, tt::tolerance(1e-6));
+    }
+
+    {
+        double const new_discharge{lue::detail::iterate_to_new_discharge<double>(
+            1,     // upstream_discharge
+            1,     // current_discharge
+            -2,    // lateral_inflow
+            1.5,   // alpha
+            0.6,   // beta
+            15,    // time_step_duration
+            10)};  // channel_length
+        double const discharge_we_want{0};
+
+        BOOST_TEST(new_discharge == discharge_we_want, tt::tolerance(1e-6));
+    }
+
+    {
+        double const new_discharge{lue::detail::iterate_to_new_discharge<double>(
+            0,     // upstream_discharge
+            1,     // current_discharge
+            -1,    // lateral_inflow
+            1.5,   // alpha
+            0.6,   // beta
+            15,    // time_step_duration
+            10)};  // channel_length
+        double const discharge_we_want{0};
+
+        BOOST_TEST(new_discharge == discharge_we_want, tt::tolerance(1e-6));
+    }
+
+    {
+        double const new_discharge{lue::detail::iterate_to_new_discharge<double>(
+            0,     // upstream_discharge
+            0,     // current_discharge
+            -1,    // lateral_inflow
+            1.5,   // alpha
+            0.6,   // beta
+            15,    // time_step_duration
+            10)};  // channel_length
+        double const discharge_we_want{0};
+
+        BOOST_TEST(new_discharge == discharge_we_want, tt::tolerance(1e-6));
+    }
+}
+
+
+// BOOST_AUTO_TEST_CASE(crashed_in_pcraster1)
+// {
+//     double const new_discharge{lue::detail::iterate_to_new_discharge<double>(
+//         0.000201343,   // upstream_discharge
+//         0.000115866,   // current_discharge
+//         -0.000290263,  // lateral_inflow
+//         1.73684,       // alpha
+//         0.6,           // beta
+//         15,            // time_step_duration
+//         10)};          // channel_length
+//     // /* epsilon */ 1E-12);  // epsilon
+//     double const discharge_we_want{0.000031450866300937};
+//
+//     BOOST_TEST(new_discharge == discharge_we_want, tt::tolerance(1e-6));
+// }
+
+
+BOOST_AUTO_TEST_CASE(crashed_in_pcraster2)
+{
+    double const new_discharge{lue::detail::iterate_to_new_discharge<double>(
+        0,             // upstream_discharge
+        1.11659e-07,   // current_discharge
+        -1.32678e-05,  // lateral_inflow
+        1.6808,        // alpha
+        0.6,           // beta
+        15,            // time_step_duration
+        10)};          // channel_length
+
+    double const discharge_we_want{std::numeric_limits<double>::min()};
+
+    BOOST_TEST(new_discharge == discharge_we_want, tt::tolerance(1e-6));
+}
+
+
+BOOST_AUTO_TEST_CASE(random_input)
+{
+    using FloatElement = double;
+
+    std::random_device random_device{};
+    std::default_random_engine random_number_engine(random_device());
+
+    std::uniform_real_distribution<FloatElement> discharge_distribution{0, 1000};
+    std::uniform_real_distribution<FloatElement> lateral_inflow_distribution{-1000, 1000};
+    std::uniform_real_distribution<FloatElement> alpha_distribution{0.5, 6.0};
+    std::uniform_real_distribution<FloatElement> beta_distribution{0.5, 2.0};
+    std::uniform_real_distribution<FloatElement> time_step_duration_distribution{1, 100};
+    std::uniform_real_distribution<FloatElement> channel_length_distribution{1, 100};
+
+    for (std::size_t i = 0; i < 10000; ++i)
+    {
+        FloatElement const upstream_discharge{discharge_distribution(random_number_engine)};
+        FloatElement const current_discharge{discharge_distribution(random_number_engine)};
+        FloatElement const alpha{alpha_distribution(random_number_engine)};
+        FloatElement const beta{beta_distribution(random_number_engine)};
+        FloatElement const time_step_duration{time_step_duration_distribution(random_number_engine)};
+        FloatElement const channel_length{channel_length_distribution(random_number_engine)};
+
+        FloatElement const lateral_inflow{lateral_inflow_distribution(random_number_engine) / channel_length};
+
+        FloatElement new_discharge{-1};
+
+        BOOST_TEST_INFO(fmt::format(
+            "upstream_discharge: {}\n"
+            "current_discharge: {}\n"
+            "lateral_inflow: {}\n"
+            "alpha: {}\n"
+            "beta: {}\n"
+            "time_step_duration: {}\n"
+            "channel_length: {}\n",
+            upstream_discharge,
+            current_discharge,
+            lateral_inflow,
+            alpha,
+            beta,
+            time_step_duration,
+            channel_length));
+
+        // This call should not throw an exception
+        new_discharge = lue::detail::iterate_to_new_discharge<FloatElement>(
+            upstream_discharge,
+            current_discharge,
+            lateral_inflow,
+            alpha,
+            beta,
+            time_step_duration,
+            channel_length);
+
+        BOOST_CHECK(new_discharge >= 0);
+    }
+}