From 9f5e224510a200b3f7df16b093e6d5311ed26288 Mon Sep 17 00:00:00 2001 From: Willem van Verseveld Date: Tue, 3 Dec 2024 08:23:30 +0100 Subject: [PATCH] Stable timestep kinematic wave and local inertial Fix computing `dt_min` for kinematic wave: value of NaN could occur when only one cell q > 0.0 (quantile!). Increased default `dt_min` for local inertial routing to 60.0 s. --- src/routing/surface_kinwave.jl | 12 ++++++++---- src/routing/surface_local_inertial.jl | 4 ++-- test/run_sbm.jl | 28 +++++++++++++-------------- test/run_sbm_gwf.jl | 2 +- 4 files changed, 25 insertions(+), 21 deletions(-) diff --git a/src/routing/surface_kinwave.jl b/src/routing/surface_kinwave.jl index dd48e014b..7160117aa 100644 --- a/src/routing/surface_kinwave.jl +++ b/src/routing/surface_kinwave.jl @@ -538,12 +538,16 @@ function stable_timestep( stable_timesteps[k] = (flow_length[i] / c) end end - if k > 0 - dt_s = quantile!(@view(stable_timesteps[1:k]), p) + + dt_min = if k == 1 + stable_timesteps[k] + elseif k > 0 + quantile!(@view(stable_timesteps[1:k]), p) else - dt_s = 600.0 + 600.0 end - return dt_s + + return dt_min end """ diff --git a/src/routing/surface_local_inertial.jl b/src/routing/surface_local_inertial.jl index 4987afd21..3050a9a11 100644 --- a/src/routing/surface_local_inertial.jl +++ b/src/routing/surface_local_inertial.jl @@ -839,7 +839,7 @@ function stable_timestep(model::LocalInertialRiverFlow{T})::T where {T} @fastmath @inbounds dt = cfl * flow_length[i] / sqrt(g * h[i]) dt_min = min(dt, dt_min) end - dt_min = isinf(dt_min) ? T(10.0) : dt_min + dt_min = isinf(dt_min) ? T(60.0) : dt_min return dt_min end @@ -856,7 +856,7 @@ function stable_timestep(model::LocalInertialOverlandFlow{T})::T where {T} end dt_min = min(dt, dt_min) end - dt_min = isinf(dt_min) ? T(10.0) : dt_min + dt_min = isinf(dt_min) ? T(60.0) : dt_min return dt_min end diff --git a/test/run_sbm.jl b/test/run_sbm.jl index 745dbddb6..282c5bc2d 100644 --- a/test/run_sbm.jl +++ b/test/run_sbm.jl @@ -266,11 +266,11 @@ Wflow.run_timestep!(model) @testset "river flow and depth (local inertial)" begin q = model.lateral.river.variables.q_av @test sum(q) ≈ 3854.717278465037f0 - @test q[1622] ≈ 7.311274274855274f-5 + @test q[1622] ≈ 7.296767063082754f-5 @test q[43] ≈ 11.716766734364437f0 @test q[501] ≈ 3.4819773071884716f0 h = model.lateral.river.variables.h_av - @test h[1622] ≈ 0.001987887644883981f0 + @test h[1622] ≈ 0.001986669483044286f0 @test h[43] ≈ 0.43311924038778815f0 @test h[501] ≈ 0.05635210581824346f0 q_channel = model.lateral.river.variables.q_channel_av @@ -289,17 +289,17 @@ Wflow.run_timestep!(model) @testset "river and overland flow and depth (local inertial)" begin q = model.lateral.river.variables.q_av @test sum(q) ≈ 2380.64389229669f0 - @test q[1622] ≈ 7.328535246760549f-5 + @test q[1622] ≈ 7.30561606758937f-5 @test q[43] ≈ 5.3566292152594155f0 - @test q[501] ≈ 1.6042388573126602f0 + @test q[501] ≈ 1.602564408503896f0 h = model.lateral.river.variables.h_av - @test h[1622] ≈ 0.0019891342000364796f0 + @test h[1622] ≈ 0.001987528017923597f0 @test h[43] ≈ 0.30026439683630496f0 - @test h[501] ≈ 0.03195324587192846f0 + @test h[501] ≈ 0.031933708617123746f0 qx = model.lateral.land.variables.qx qy = model.lateral.land.variables.qy - @test qx[[26, 35, 631]] ≈ [0.1939736998417174f0, 0.026579954465883678f0, 0.0f0] - @test qy[[26, 35, 631]] ≈ [0.12906530420401777f0, 1.7225115950614904f0, 0.0f0] + @test qx[[26, 35, 631]] ≈ [0.18343478752498582f0, 0.000553471702071059f0, 0.0f0] + @test qy[[26, 35, 631]] ≈ [0.12607229901243375f0, 0.019605967561619194f0, 0.0f0] h = model.lateral.land.variables.h @test h[[26, 35, 631]] ≈ [0.07367301172613304f0, 0.009139882310161706f0, 0.0007482998926237368f0] @@ -420,12 +420,12 @@ Wflow.run_timestep!(model) @testset "river flow (local inertial) with floodplain schematization simulation" begin q = model.lateral.river.variables.q_av @test sum(q) ≈ 3843.944494991296f0 - @test q[1622] ≈ 7.311274278896417f-5 + @test q[1622] ≈ 7.296767071929629f-5 @test q[43] ≈ 11.716766734364418f0 @test q[501] ≈ 3.424364314225289f0 - @test q[5808] ≈ 0.0022274679501657693f0 + @test q[5808] ≈ 0.002228981516146531f0 h = model.lateral.river.variables.h_av - @test h[1622] ≈ 0.001987887580593841f0 + @test h[1622] ≈ 0.0019866694251020806f0 @test h[43] ≈ 0.433119240388070f0 @test h[501] ≈ 0.055832770820860404f0 @test h[5808] ≈ 0.005935591961908253f0 @@ -441,12 +441,12 @@ Wflow.run_timestep!(model) @testset "change boundary condition for local inertial routing (including floodplain)" begin q = model.lateral.river.variables.q_av @test sum(q) ≈ 3844.1544889903134f0 - @test q[1622] ≈ 7.311151138631112f-5 + @test q[1622] ≈ 7.296767071929629f-5 @test q[43] ≈ 11.716766734416717f0 @test q[501] ≈ 3.424329413571391f0 - @test q[5808] ≈ 0.060518234525259465f0 + @test q[5808] ≈ 0.055269620065756024f0 h = model.lateral.river.variables.h_av - @test h[1622] ≈ 0.0019878952928530183f0 + @test h[1622] ≈ 0.0019866694251020806f0 @test h[43] ≈ 0.4331192403230577f0 @test h[501] ≈ 0.0558281185092927f0 @test h[5808] ≈ 2.0000006940603936f0 diff --git a/test/run_sbm_gwf.jl b/test/run_sbm_gwf.jl index 7fd9eb4d3..81038ed5a 100644 --- a/test/run_sbm_gwf.jl +++ b/test/run_sbm_gwf.jl @@ -96,7 +96,7 @@ Wflow.run_timestep!(model) @test sum(q) ≈ 0.02727911500112358f0 @test q[6] ≈ 0.006111263175002127f0 @test river.variables.volume[6] ≈ 7.6120096530771075f0 - @test river.boundary_conditions.inwater[6] ≈ 0.00022087679662860144f0 + @test river.boundary_conditions.inwater[6] ≈ 0.0002210785332342944f0 @test q[13] ≈ 0.0004638698607639214f0 @test q[5] ≈ 0.0064668491697542786f0 end