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 +// #define BOOST_MATH_INSTRUMENT #include #include #include @@ -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::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(discharge_guess, min_discharge); } lue_hpx_assert(discharge_guess > Float{0}); @@ -76,40 +92,48 @@ namespace lue { } - std::pair operator()(Float const new_discharge) const + auto operator()(Float const new_discharge) const -> std::pair { 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 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::max()}; + int const digits = static_cast(std::numeric_limits::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::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 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(std::numeric_limits::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; 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; - 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( 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>( + auto const current_discharge = lue::test::create_partitioned_array>( 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(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>( 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(array_shape, partition_shape, 0); + auto const current_discharge = lue::create_partitioned_array(array_shape, partition_shape, 0); auto const inflow = lue::create_partitioned_array(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(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(array_shape, partition_shape, 1); + auto const current_discharge = lue::create_partitioned_array(array_shape, partition_shape, 1); auto const inflow = lue::create_partitioned_array(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(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(array_shape, partition_shape, 0); + auto const current_discharge = lue::create_partitioned_array(array_shape, partition_shape, 0); auto const inflow = lue::create_partitioned_array(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(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( + 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( + 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( + 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( + 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( + 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( +// 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( + 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::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 discharge_distribution{0, 1000}; + std::uniform_real_distribution lateral_inflow_distribution{-1000, 1000}; + std::uniform_real_distribution alpha_distribution{0.5, 6.0}; + std::uniform_real_distribution beta_distribution{0.5, 2.0}; + std::uniform_real_distribution time_step_duration_distribution{1, 100}; + std::uniform_real_distribution 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( + upstream_discharge, + current_discharge, + lateral_inflow, + alpha, + beta, + time_step_duration, + channel_length); + + BOOST_CHECK(new_discharge >= 0); + } +}