|
53 | 53 | Compute the area and level of a basin given its storage. |
54 | 54 | """ |
55 | 55 | function get_area_and_level(basin::Basin, state_idx::Int, storage::T)::Tuple{T, T} where {T} |
56 | | - level = basin.storage_to_level[state_idx](max(storage, 0.0)) |
57 | | - area = basin.level_to_area[state_idx](level) |
58 | | - |
| 56 | + storage_to_level = basin.storage_to_level[state_idx] |
| 57 | + level_to_area = basin.level_to_area[state_idx] |
| 58 | + if storage >= 0 |
| 59 | + level = storage_to_level(storage) |
| 60 | + else |
| 61 | + # Negative storage is not feasible and this yields a level |
| 62 | + # below the basin bottom, but this does yield usable gradients |
| 63 | + # for the non-linear solver |
| 64 | + bottom = first(level_to_area.t) |
| 65 | + level = bottom + derivative(storage_to_level, 0.0) * storage |
| 66 | + end |
| 67 | + area = level_to_area(level) |
59 | 68 | return area, level |
60 | 69 | end |
61 | 70 |
|
|
887 | 896 | (A::AbstractInterpolation)(t::GradientTracer) = t |
888 | 897 | reduction_factor(x::GradientTracer, threshold::Real) = x |
889 | 898 | relaxed_root(x::GradientTracer, threshold::Real) = x |
| 899 | +get_area_and_level(basin::Basin, state_idx::Int, storage::GradientTracer) = storage, storage |
| 900 | +stop_declining_negative_storage!(du, u::ComponentVector{<:GradientTracer}) = nothing |
| 901 | + |
| 902 | +@kwdef struct MonitoredBackTracking{B, V} |
| 903 | + linesearch::B = BackTracking() |
| 904 | + dz_tmp::V = [] |
| 905 | + z_tmp::V = [] |
| 906 | +end |
| 907 | + |
| 908 | +""" |
| 909 | +Compute the residual of the non-linear solver, i.e. a measure of the |
| 910 | +error in the solution to the implicit equation defined by the solver algorithm |
| 911 | +""" |
| 912 | +function residual(z, integrator, nlsolver, f) |
| 913 | + (; uprev, t, p, dt, opts, isdae) = integrator |
| 914 | + (; tmp, ztmp, γ, α, cache, method) = nlsolver |
| 915 | + (; ustep, atmp, tstep, k, invγdt, tstep, k, invγdt) = cache |
| 916 | + if isdae |
| 917 | + _uprev = get_dae_uprev(integrator, uprev) |
| 918 | + b, ustep2 = |
| 919 | + _compute_rhs!(tmp, ztmp, ustep, α, tstep, k, invγdt, p, _uprev, f::TF, z) |
| 920 | + else |
| 921 | + b, ustep2 = |
| 922 | + _compute_rhs!(tmp, ztmp, ustep, γ, α, tstep, k, invγdt, method, p, dt, f, z) |
| 923 | + end |
| 924 | + calculate_residuals!( |
| 925 | + atmp, |
| 926 | + b, |
| 927 | + uprev, |
| 928 | + ustep2, |
| 929 | + opts.abstol, |
| 930 | + opts.reltol, |
| 931 | + opts.internalnorm, |
| 932 | + t, |
| 933 | + ) |
| 934 | + ndz = opts.internalnorm(atmp, t) |
| 935 | + return ndz |
| 936 | +end |
| 937 | + |
| 938 | +""" |
| 939 | +MonitoredBackTracing is a thin wrapper of BackTracking, making sure that |
| 940 | +the BackTracking relaxation is rejected if it results in a residual increase |
| 941 | +""" |
| 942 | +function OrdinaryDiffEq.relax!( |
| 943 | + dz, |
| 944 | + nlsolver::AbstractNLSolver, |
| 945 | + integrator::DEIntegrator, |
| 946 | + f, |
| 947 | + linesearch::MonitoredBackTracking, |
| 948 | +) |
| 949 | + (; linesearch, dz_tmp, z_tmp) = linesearch |
| 950 | + |
| 951 | + # Store step before relaxation |
| 952 | + @. dz_tmp = dz |
| 953 | + |
| 954 | + # Apply relaxation and measure the residual change |
| 955 | + @. z_tmp = nlsolver.z + dz |
| 956 | + resid_before = residual(z_tmp, integrator, nlsolver, f) |
| 957 | + relax!(dz, nlsolver, integrator, f, linesearch) |
| 958 | + @. z_tmp = nlsolver.z + dz |
| 959 | + resid_after = residual(z_tmp, integrator, nlsolver, f) |
| 960 | + |
| 961 | + # If the residual increased due to the relaxation, reject it |
| 962 | + if resid_after > resid_before |
| 963 | + @. dz = dz_tmp |
| 964 | + end |
| 965 | +end |
0 commit comments