Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add support for negative inflow to kinematic_wave #702

Merged
merged 3 commits into from
Aug 8, 2024
Merged

Add support for negative inflow to kinematic_wave #702

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
5a0fc86
Make kinematic_wave more robust
kordejong Aug 6, 2024
e2d8b28
Treat external inflow and internal extraction seperately
kordejong Aug 7, 2024
6ec83a6
Final tweaks, for now
kordejong Aug 8, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
160 changes: 107 additions & 53 deletions source/framework/algorithm/include/lue/framework/algorithm/definition/kinematic_wave.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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>
Expand All @@ -21,53 +23,68 @@ 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},
_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};
// TODO static Float const min_discharge{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});
Expand All @@ -76,40 +93,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

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If new_discharge is updated in each iteration of newton_raphson_iterate, it should be set to non-zero values when it is calculated as negative. The code can be as follows:

auto fq(Float const new_discharge) const -> Float
                {
                   new_discharge = std::max<Float>(discharge_guess, min_discharge);

                    return _time_step_duration_over_channel_length * new_discharge +
                           _alpha * std::pow(new_discharge, _beta) - _known_terms;
                }

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Discharge can't be negative, can it? It does not make sense I think.
The assertion doesn't fire, so in practice, as the code is now, negative discharges never happen.

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

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If new_discharge is updated in each iteration of newton_raphson_iterate, it should be set to non-zero values when it is calculated as negative. The code can be as follows:

auto dfq(Float const new_discharge) const -> Float
                {
                    
                    new_discharge = std::max<Float>(discharge_guess, min_discharge);
                    
                    return _time_step_duration_over_channel_length +
                           _alpha_beta * std::pow(new_discharge, _beta - Float{1});
                }

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

new_discharge does not end up being negative so I don't think this change is necessary.

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;
Expand All @@ -122,53 +147,82 @@ 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)
// 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)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it is not necessary to check. For all cases (whether lateral inflow is greater than 0 or less than 0), the same numerical method can be used while ensuring new_discharge is set to a small non-zero values when it is calculated as negative.

code can be like as follows

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{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.
               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);
                   
                   return std::max<Float>(new_discharge, 0);

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Previously, I tried accepting negative inflows but this made the iteration not converge anymore in many cases. To solve it, I split the procedure in:

  1. Only include positive inflows to newton_raphson_iterate
  2. Use negative inflows (extractions) only to update the new discharge

This way, there is no need to handle extractions while solving for the new discharge. Also, this is what makes the procedure actually work in all cases I tested. Unless it is conceptually wrong, I suggest to keep it like it is for now.

{
// 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()};

// These brackets are crucial for obtaining a good performance
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.
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);

// if (actual_nr_iterations == max_nr_iterations)
// {
// throw std::runtime_error(fmt::format(
// "Iterating to a solution took more iterations than expected (initial guess: {}, "
// "brackets: [{}, {}])",
// discharge_guess,
// min_discharge,
// max_discharge));
// }
}

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;
Comment on lines -138 to +247
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it is not needed. The same numerical method used for positive lateral_inflow can be used for negative lateral_inflow, while new_discharge should be set to a small non-zero value when it is calculated as negative.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That is what I tried initially and I could not get it to work well. In case the extraction is relatively large, the algorithm doesn't converge to a solution. It seems to me that the numerical scheme is not meant to be used like that.

Copy link
Contributor

@saeb-faraji-gargari saeb-faraji-gargari Aug 8, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

-Just as an idea, the Newton-Raphson method may diverge when lateral_inflow is negative. In this case, I would suggest using discharge_guess as new_discharge (new_discharge = discharge_guess;) since discharge_guess is the solution to the linear discretized version of the kinematic-wave equations.

-The following criterion needs to be satisfied for the convergence of the Newton-Raphson method:
| fq * d2fq| < |dfq|^2

Ref: https://math.stackexchange.com/questions/3136446/condition-for-convergence-of-newton-raphson-method

When lateral_inflow < 0, − _known_terms (-C in Chow's book) becomes positive, may lead to a large positive value for fq. This can violate the convergence criteria of the Newton-Raphson method. This issue can be considered in the future to improve the kinematic-wave operator.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK, let's look into that in the near future. I have started initial work on #703, based on Rens' comments . Can you please copy your last comment to that ticket maybe? This PR is closed already.

}

lue_hpx_assert(nr_iterations < max_nr_iterations);
lue_hpx_assert(new_discharge >= Float{0});

return new_discharge;
Expand Down
6 changes: 0 additions & 6 deletions source/framework/algorithm/test/clump_test.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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};
Expand Down Expand Up @@ -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")
Expand Down
18 changes: 11 additions & 7 deletions source/framework/algorithm/test/flow_accumulation.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
);
}


Expand Down
Loading
Loading