Add support for negative inflow to kinematic_wave

[kordejong]

Closes #524

[saeb-faraji-gargari]

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;
                }

[kordejong]

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.

[saeb-faraji-gargari]

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});
                }

[kordejong]

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

[saeb-faraji-gargari]

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);

[kordejong]

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.

[saeb-faraji-gargari]

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.

[kordejong]

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.

[saeb-faraji-gargari]

-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.

[kordejong]

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.

[kordejong]

I propose to merge this PR, as it is an improvement compared to the current version. I will create a new issue (#703) to fine tune the implementation later. OK?

[saeb-faraji-gargari]

I propose to merge this PR, as it is an improvement compared to the current version. I will create a new issue (#703) to fine tune the implementation later. OK?

Ok. It is fine.