From 089b2c2b60d944b2704e0b56d27173c796366fde Mon Sep 17 00:00:00 2001 From: Willem van Verseveld Date: Tue, 12 Nov 2024 14:54:47 +0100 Subject: [PATCH] Refactor `SimpleReservoir` and `Lake` structs --- src/flow.jl | 44 +++--- src/io.jl | 8 +- src/reservoir_lake.jl | 350 +++++++++++++++++++++++++---------------- src/sbm_model.jl | 5 +- src/states.jl | 4 +- test/bmi.jl | 6 +- test/io.jl | 11 +- test/reservoir_lake.jl | 147 ++++++++++------- test/run_sbm.jl | 26 +-- test/sbm_config.toml | 10 +- 10 files changed, 369 insertions(+), 242 deletions(-) diff --git a/src/flow.jl b/src/flow.jl index f61176e3c..45309eb85 100644 --- a/src/flow.jl +++ b/src/flow.jl @@ -326,14 +326,14 @@ function update!(sf::SurfaceFlowRiver, network, doy) # because of possible iterations set reservoir and lake inflow and total outflow at # start to zero, the total sum of inflow and outflow at each sub time step is calculated if !isnothing(reservoir) - reservoir.inflow .= 0.0 - reservoir.totaloutflow .= 0.0 - reservoir.actevap .= 0.0 + reservoir.boundary_conditions.inflow .= 0.0 + reservoir.variables.totaloutflow .= 0.0 + reservoir.variables.actevap .= 0.0 end if !isnothing(lake) - lake.inflow .= 0.0 - lake.totaloutflow .= 0.0 - lake.actevap .= 0.0 + lake.boundary_conditions.inflow .= 0.0 + lake.variables.totaloutflow .= 0.0 + lake.variables.actevap .= 0.0 end dt, its = stable_timestep(sf) @@ -376,7 +376,7 @@ function update!(sf::SurfaceFlowRiver, network, doy) n_downstream = length(downstream_nodes) if n_downstream == 1 j = only(downstream_nodes) - qin[j] = reservoir.outflow[i] + qin[j] = reservoir.variables.outflow[i] elseif n_downstream == 0 error( """A reservoir without a downstream river node is not supported. @@ -397,7 +397,7 @@ function update!(sf::SurfaceFlowRiver, network, doy) n_downstream = length(downstream_nodes) if n_downstream == 1 j = only(downstream_nodes) - qin[j] = lake.outflow[i] + qin[j] = lake.variables.outflow[i] elseif n_downstream == 0 error( """A lake without a downstream river node is not supported. @@ -1108,7 +1108,7 @@ function shallowwater_river_update!(sw::ShallowWaterRiver, network, dt, doy, upd q_in = get_inflow_waterbody(sw, links_at_node.src[i]) update!(reservoir, v, q_in + inflow_wb[i], dt) - river_v.q[i] = reservoir.outflow[v] + river_v.q[i] = reservoir.variables.outflow[v] river_v.q_av[i] += river_v.q[i] * dt end (; lake, lake_index, inflow_wb) = sw.boundary_conditions @@ -1117,7 +1117,7 @@ function shallowwater_river_update!(sw::ShallowWaterRiver, network, dt, doy, upd q_in = get_inflow_waterbody(sw, links_at_node.src[i]) update!(lake, v, q_in + inflow_wb[i], doy, dt) - river_v.q[i] = lake.outflow[v] + river_v.q[i] = lake.variables.outflow[v] river_v.q_av[i] += river_v.q[i] * dt end if update_h @@ -1173,14 +1173,14 @@ end function update!(sw::ShallowWaterRiver{T}, network, doy; update_h = true) where {T} (; reservoir, lake) = sw.boundary_conditions if !isnothing(reservoir) - reservoir.inflow .= 0.0 - reservoir.totaloutflow .= 0.0 - reservoir.actevap .= 0.0 + reservoir.boundary_conditions.inflow .= 0.0 + reservoir.variables.totaloutflow .= 0.0 + reservoir.variables.actevap .= 0.0 end if !isnothing(lake) - lake.inflow .= 0.0 - lake.totaloutflow .= 0.0 - lake.actevap .= 0.0 + lake.boundary_conditions.inflow .= 0.0 + lake.variables.totaloutflow .= 0.0 + lake.variables.actevap .= 0.0 end if !isnothing(sw.floodplain) sw.floodplain.variables.q_av .= 0.0 @@ -1485,14 +1485,14 @@ function update!( (; reservoir, lake) = swr.boundary_conditions if !isnothing(reservoir) - reservoir.inflow .= 0.0 - reservoir.totaloutflow .= 0.0 - reservoir.actevap .= 0.0 + reservoir.boundary_conditions.inflow .= 0.0 + reservoir.variables.totaloutflow .= 0.0 + reservoir.variables.actevap .= 0.0 end if !isnothing(lake) - lake.inflow .= 0.0 - lake.totaloutflow .= 0.0 - lake.actevap .= 0.0 + lake.boundary_conditions.inflow .= 0.0 + lake.variables.totaloutflow .= 0.0 + lake.variables.actevap .= 0.0 end swr.variables.q_av .= 0.0 swr.variables.h_av .= 0.0 diff --git a/src/io.jl b/src/io.jl index 8c4d5be28..90d2da595 100644 --- a/src/io.jl +++ b/src/io.jl @@ -184,18 +184,18 @@ end function get_param_res(model) return Dict( symbols"vertical.atmospheric_forcing.precipitation" => - model.lateral.river.boundary_conditions.reservoir.precipitation, + model.lateral.river.boundary_conditions.reservoir.boundary_conditions.precipitation, symbols"vertical.atmospheric_forcing.potential_evaporation" => - model.lateral.river.boundary_conditions.reservoir.evaporation, + model.lateral.river.boundary_conditions.reservoir.boundary_conditions.evaporation, ) end function get_param_lake(model) return Dict( symbols"vertical.atmospheric_forcing.precipitation" => - model.lateral.river.boundary_conditions.lake.precipitation, + model.lateral.river.boundary_conditions.lake.boundary_conditions.precipitation, symbols"vertical.atmospheric_forcing.potential_evaporation" => - model.lateral.river.boundary_conditions.lake.evaporation, + model.lateral.river.boundary_conditions.lake.boundary_conditions.evaporation, ) end diff --git a/src/reservoir_lake.jl b/src/reservoir_lake.jl index edc55a91f..e7d8a6577 100644 --- a/src/reservoir_lake.jl +++ b/src/reservoir_lake.jl @@ -1,28 +1,15 @@ -@get_units @grid_loc @with_kw struct SimpleReservoir{T} - dt::T # Model time step [s] - maxvolume::Vector{T} | "m3" # maximum storage (above which water is spilled) [m³] - area::Vector{T} | "m2" # reservoir area [m²] - maxrelease::Vector{T} | "m3 s-1" # maximum amount that can be released if below spillway [m³ s⁻¹] - demand::Vector{T} | "m3 s-1" # minimum (environmental) flow requirement downstream of the reservoir [m³ s⁻¹] - targetminfrac::Vector{T} | "-" # target minimum full fraction (of max storage) [-] - targetfullfrac::Vector{T} | "-" # target fraction full (of max storage) [-] - volume::Vector{T} | "m3" # reservoir volume [m³] - inflow::Vector{T} | "m3" # total inflow into reservoir [m³] - outflow::Vector{T} | "m3 s-1" # outflow from reservoir [m³ s⁻¹] - totaloutflow::Vector{T} | "m3" # total outflow from reservoir [m³] - percfull::Vector{T} | "-" # fraction full (of max storage) [-] - demandrelease::Vector{T} | "m3 s-1" # minimum (environmental) flow released from reservoir [m³ s⁻¹] - precipitation::Vector{T} # average precipitation for reservoir area [mm Δt⁻¹] - evaporation::Vector{T} # average potential evaporation for reservoir area [mm Δt⁻¹] - actevap::Vector{T} # average actual evaporation for reservoir area [mm Δt⁻¹] - - function SimpleReservoir{T}(args...) where {T} - equal_size_vectors(args) - return new(args...) - end + +@get_units @grid_loc @with_kw struct ReservoirParameters{T} + dt::T # Model time step [s] + maxvolume::Vector{T} | "m3" # maximum storage (above which water is spilled) [m³] + area::Vector{T} | "m2" # reservoir area [m²] + maxrelease::Vector{T} | "m3 s-1" # maximum amount that can be released if below spillway [m³ s⁻¹] + demand::Vector{T} | "m3 s-1" # minimum (environmental) flow requirement downstream of the reservoir [m³ s⁻¹] + targetminfrac::Vector{T} | "-" # target minimum full fraction (of max storage) [-] + targetfullfrac::Vector{T} | "-" # target fraction full (of max storage end -function initialize_simple_reservoir(config, nc, inds_riv, nriv, pits, dt) +function ReservoirParameters(config, nc, inds_riv, nriv, pits, dt) # read only reservoir data if reservoirs true # allow reservoirs only in river cells # note that these locations are only the reservoir outlet pixels @@ -129,9 +116,13 @@ function initialize_simple_reservoir(config, nc, inds_riv, nriv, pits, dt) # all upstream flow goes to the river and flows into the reservoir pits[inds_res] .= true - n = length(resarea) - @info "Read `$n` reservoir locations." - reservoirs = SimpleReservoir{Float}(; + reservoir_indices = ( + indices_outlet = inds_res, + indices_coverage = inds_res_cov, + reverse_indices = rev_inds_reservoir, + ) + + parameters = ReservoirParameters{Float}(; dt = dt, demand = resdemand, maxrelease = resmaxrelease, @@ -139,25 +130,66 @@ function initialize_simple_reservoir(config, nc, inds_riv, nriv, pits, dt) area = resarea, targetfullfrac = res_targetfullfrac, targetminfrac = res_targetminfrac, - volume = res_targetfullfrac .* resmaxvolume, - inflow = fill(mv, n), + ) + + return parameters, reservoir_indices, resindex, pits +end + +@get_units @grid_loc @with_kw struct ReservoirVariables{T} + volume::Vector{T} | "m3" # reservoir volume [m³] + outflow::Vector{T} | "m3 s-1" # outflow from reservoir [m³ s⁻¹] + totaloutflow::Vector{T} | "m3" # total outflow from reservoir [m³] + percfull::Vector{T} | "-" # fraction full (of max storage) [-] + demandrelease::Vector{T} | "m3 s-1" # minimum (environmental) flow released from reservoir [m³ s⁻¹] + actevap::Vector{T} # average actual evaporation for reservoir area [mm Δt⁻¹] +end + +function ReservoirVariables(n, parameters) + (; targetfullfrac, maxvolume) = parameters + variables = ReservoirVariables{Float}(; + volume = targetfullfrac .* maxvolume, outflow = fill(mv, n), totaloutflow = fill(mv, n), percfull = fill(mv, n), demandrelease = fill(mv, n), + actevap = fill(mv, n), + ) + return variables +end + +@get_units @grid_loc @with_kw struct ReservoirBC{T} + inflow::Vector{T} | "m3" # total inflow into reservoir [m³] + precipitation::Vector{T} # average precipitation for reservoir area [mm Δt⁻¹] + evaporation::Vector{T} # average potential evaporation for reservoir area [mm Δt⁻¹] +end + +function ReservoirBC(n) + bc = ReservoirBC{Float}(; + inflow = fill(mv, n), precipitation = fill(mv, n), evaporation = fill(mv, n), - actevap = fill(mv, n), ) + return bc +end - return reservoirs, - resindex, - ( - indices_outlet = inds_res, - indices_coverage = inds_res_cov, - reverse_indices = rev_inds_reservoir, - ), - pits +@with_kw struct SimpleReservoir{T} + boundary_conditions::ReservoirBC{T} + parameters::ReservoirParameters{T} + variables::ReservoirVariables{T} +end + +function SimpleReservoir(config, nc, inds_riv, nriv, pits, dt) + parameters, reservoir_indices, resindex, pits = + ReservoirParameters(config, nc, inds_riv, nriv, pits, dt) + + n = length(parameters.area) + @info "Read `$n` reservoir locations." + + variables = ReservoirVariables(n, parameters) + boundary_conditions = ReservoirBC(n) + reservoirs = SimpleReservoir{Float}(; boundary_conditions, parameters, variables) + + return reservoirs, resindex, reservoir_indices, pits end """ @@ -167,70 +199,62 @@ This is called from within the kinematic wave loop, therefore updating only for element rather than all at once. """ function update!(res::SimpleReservoir, i, inflow, timestepsecs) + res_bc = res.boundary_conditions + res_p = res.parameters + res_v = res.variables # limit lake evaporation based on total available volume [m³] - precipitation = 0.001 * res.precipitation[i] * (timestepsecs / res.dt) * res.area[i] - available_volume = res.volume[i] + inflow * timestepsecs + precipitation - evap = 0.001 * res.evaporation[i] * (timestepsecs / res.dt) * res.area[i] + precipitation = + 0.001 * res_bc.precipitation[i] * (timestepsecs / res_p.dt) * res_p.area[i] + available_volume = res_v.volume[i] + inflow * timestepsecs + precipitation + evap = 0.001 * res_bc.evaporation[i] * (timestepsecs / res_p.dt) * res_p.area[i] actevap = min(available_volume, evap) # [m³/timestepsecs] - vol = res.volume[i] + (inflow * timestepsecs) + precipitation - actevap + vol = res_v.volume[i] + (inflow * timestepsecs) + precipitation - actevap vol = max(vol, 0.0) - percfull = vol / res.maxvolume[i] + percfull = vol / res_p.maxvolume[i] # first determine minimum (environmental) flow using a simple sigmoid curve to scale for target level - fac = scurve(percfull, res.targetminfrac[i], Float(1.0), Float(30.0)) - demandrelease = min(fac * res.demand[i] * timestepsecs, vol) + fac = scurve(percfull, res_p.targetminfrac[i], Float(1.0), Float(30.0)) + demandrelease = min(fac * res_p.demand[i] * timestepsecs, vol) vol = vol - demandrelease - wantrel = max(0.0, vol - (res.maxvolume[i] * res.targetfullfrac[i])) + wantrel = max(0.0, vol - (res_p.maxvolume[i] * res_p.targetfullfrac[i])) # Assume extra maximum Q if spilling - overflow_q = max((vol - res.maxvolume[i]), 0.0) - torelease = min(wantrel, overflow_q + res.maxrelease[i] * timestepsecs - demandrelease) + overflow_q = max((vol - res_p.maxvolume[i]), 0.0) + torelease = + min(wantrel, overflow_q + res_p.maxrelease[i] * timestepsecs - demandrelease) vol = vol - torelease outflow = torelease + demandrelease - percfull = vol / res.maxvolume[i] + percfull = vol / res_p.maxvolume[i] # update values in place - res.outflow[i] = outflow / timestepsecs - res.inflow[i] += inflow * timestepsecs - res.totaloutflow[i] += outflow - res.demandrelease[i] = demandrelease / timestepsecs - res.percfull[i] = percfull - res.volume[i] = vol - res.actevap[i] += 1000.0 * (actevap / res.area[i]) + res_v.outflow[i] = outflow / timestepsecs + res_bc.inflow[i] += inflow * timestepsecs + res_v.totaloutflow[i] += outflow + res_v.demandrelease[i] = demandrelease / timestepsecs + res_v.percfull[i] = percfull + res_v.volume[i] = vol + res_v.actevap[i] += 1000.0 * (actevap / res_p.area[i]) return nothing end -@get_units @grid_loc @with_kw struct Lake{T} - dt::T # Model time step [s] - lowerlake_ind::Vector{Int} | "-" # Index of lower lake (linked lakes) - area::Vector{T} | "m2" # lake area [m²] - maxstorage::Vector{Union{T, Missing}} | "m3" # lake maximum storage from rating curve 1 [m³] - threshold::Vector{T} | "m" # water level threshold H₀ [m] below that level outflow is zero - storfunc::Vector{Int} | "-" # type of lake storage curve, 1: S = AH, 2: S = f(H) from lake data and interpolation - outflowfunc::Vector{Int} | "-" # type of lake rating curve, 1: Q = f(H) from lake data and interpolation, 2: General Q = b(H - H₀)ᵉ, 3: Case of Puls Approach Q = b(H - H₀)² - b::Vector{T} | "m3/2 s-1 (if e=3/2)" # rating curve coefficient - e::Vector{T} | "-" # rating curve exponent - sh::Vector{Union{SH, Missing}} # data for storage curve - hq::Vector{Union{HQ, Missing}} # data for rating curve - waterlevel::Vector{T} | "m" # waterlevel H [m] of lake - inflow::Vector{T} | "m3" # inflow to the lake [m³] - storage::Vector{T} | "m3" # storage lake [m³] - outflow::Vector{T} | "m3 s-1" # outflow lake [m³ s⁻¹] - totaloutflow::Vector{T} | "m3" # total outflow lake [m³] - precipitation::Vector{T} # average precipitation for lake area [mm Δt⁻¹] - evaporation::Vector{T} # average potential evaporation for lake area [mm Δt⁻¹] - actevap::Vector{T} # average actual evapotranspiration for lake area [mm Δt⁻¹] - - function Lake{T}(args...) where {T} - equal_size_vectors(args) - return new(args...) - end +@get_units @grid_loc @with_kw struct LakeParameters{T} + dt::T # Model time step [s] + lowerlake_ind::Vector{Int} | "-" # Index of lower lake (linked lakes) + area::Vector{T} | "m2" # lake area [m²] + maxstorage::Vector{Union{T, Missing}} | "m3" # lake maximum storage from rating curve 1 [m³] + threshold::Vector{T} | "m" # water level threshold H₀ [m] below that level outflow is zero + storfunc::Vector{Int} | "-" # type of lake storage curve, 1: S = AH, 2: S = f(H) from lake data and interpolation + outflowfunc::Vector{Int} | "-" # type of lake rating curve, 1: Q = f(H) from lake data and interpolation, 2: General Q = b(H - H₀)ᵉ, 3: Case of Puls Approach Q = b(H - H₀)² + b::Vector{T} | "m3/2 s-1 (if e=3/2)" # rating curve coefficient + e::Vector{T} | "-" # rating curve exponent + sh::Vector{Union{SH, Missing}} # data for storage curve + hq::Vector{Union{HQ, Missing}} # data for rating curve end -function initialize_lake(config, nc, inds_riv, nriv, pits, dt) +function LakeParameters(config, nc, inds_riv, nriv, pits, dt) # read only lake data if lakes true # allow lakes only in river cells # note that these locations are only the lake outlet pixels @@ -396,8 +420,7 @@ function initialize_lake(config, nc, inds_riv, nriv, pits, dt) ) end end - n = length(lakearea) - lakes = Lake{Float}(; + parameters = LakeParameters{Float}(; dt = dt, lowerlake_ind = lowerlake_ind, area = lakearea, @@ -407,26 +430,69 @@ function initialize_lake(config, nc, inds_riv, nriv, pits, dt) outflowfunc = lake_outflowfunc, b = lake_b, e = lake_e, - waterlevel = lake_waterlevel, sh = sh, hq = hq, + ) + lake_indices = ( + indices_outlet = inds_lake, + indices_coverage = inds_lake_cov, + reverse_indices = rev_inds_lake, + ) + return parameters, lake_indices, lakeindex, lake_Waterlevel, pits +end + +@get_units @grid_loc @with_kw struct LakeVariables{T} + waterlevel::Vector{T} | "m" # waterlevel H [m] of lake + storage::Vector{T} | "m3" # storage lake [m³] + outflow::Vector{T} | "m3 s-1" # outflow lake [m³ s⁻¹] + totaloutflow::Vector{T} | "m3" # total outflow lake [m³] + actevap::Vector{T} # average actual evapotranspiration for lake area [mm Δt⁻¹] +end + +function LakeVariables(n, lake_waterlevel) + variables = LakeVariables{Float}(; + waterlevel = lake_waterlevel, inflow = fill(mv, n), storage = initialize_storage(lake_storfunc, lakearea, lake_waterlevel, sh), outflow = fill(mv, n), totaloutflow = fill(mv, n), + actevap = fill(mv, n), + ) + return variables +end + +@get_units @grid_loc @with_kw struct LakeBC{T} + inflow::Vector{T} | "m3" # inflow to the lake [m³] + precipitation::Vector{T} # average precipitation for lake area [mm Δt⁻¹] + evaporation::Vector{T} # average potential evaporation for lake area [mm Δt⁻¹] +end + +function LakeBC(n) + bc = LakeBC{Float}(; + inflow = fill(mv, n), precipitation = fill(mv, n), evaporation = fill(mv, n), - actevap = fill(mv, n), ) + return bc +end - return lakes, - lakeindex, - ( - indices_outlet = inds_lake, - indices_coverage = inds_lake_cov, - reverse_indices = rev_inds_lake, - ), - pits +@with_kw struct Lake{T} + boundary_conditions::LakeBC{T} + parameters::LakeParameters{T} + variables::LakeVariables{T} +end + +function Lake(config, nc, inds_riv, nriv, pits, dt) + parameters, lake_indices, lakeindex, lake_waterlevel, pits = + LakeParameters(config, nc, inds_riv, nriv, pits, dt) + + n = length(parameters.area) + variables = LakeVariables(n, lake_waterlevel) + boundary_conditions = LakeBC(n) + + lakes = Lake{Float}(; boundary_conditions, parameters, variables) + + return lakes, lakeindex, lake_indices, pits end "Determine the initial storage depending on the storage function" @@ -492,23 +558,28 @@ This is called from within the kinematic wave loop, therefore updating only for element rather than all at once. """ function update!(lake::Lake, i, inflow, doy, timestepsecs) - lo = lake.lowerlake_ind[i] + lake_bc = lake.boundary_conditions + lake_p = lake.parameters + lake_v = lake.variables + + lo = lake_p.lowerlake_ind[i] has_lowerlake = lo != 0 # limit lake evaporation based on total available volume [m³] - precipitation = 0.001 * lake.precipitation[i] * (timestepsecs / lake.dt) * lake.area[i] - available_volume = lake.storage[i] + inflow * timestepsecs + precipitation - evap = 0.001 * lake.evaporation[i] * (timestepsecs / lake.dt) * lake.area[i] + precipitation = + 0.001 * lake_bc.precipitation[i] * (timestepsecs / lake_p.dt) * lake_p.area[i] + available_volume = lake_v.storage[i] + inflow * timestepsecs + precipitation + evap = 0.001 * lake_bc.evaporation[i] * (timestepsecs / lake_p.dt) * lake_p.area[i] actevap = min(available_volume, evap) # [m³/timestepsecs] ### Modified Puls Approach (Burek et al., 2013, LISFLOOD) ### # outflowfunc = 3 # Calculate lake factor and SI parameter - if lake.outflowfunc[i] == 3 - lakefactor = lake.area[i] / (timestepsecs * pow(lake.b[i], 0.5)) - si_factor = (lake.storage[i] + precipitation - actevap) / timestepsecs + inflow + if lake_p.outflowfunc[i] == 3 + lakefactor = lake_p.area[i] / (timestepsecs * pow(lake_p.b[i], 0.5)) + si_factor = (lake_v.storage[i] + precipitation - actevap) / timestepsecs + inflow # Adjust SIFactor for ResThreshold != 0 - si_factor_adj = si_factor - lake.area[i] * lake.threshold[i] / timestepsecs + si_factor_adj = si_factor - lake_p.area[i] * lake_p.threshold[i] / timestepsecs # Calculate the new lake outflow/waterlevel/storage if si_factor_adj > 0.0 quadratic_sol_term = @@ -523,33 +594,37 @@ function update!(lake::Lake, i, inflow, doy, timestepsecs) end outflow = min(outflow, si_factor) storage = (si_factor - outflow) * timestepsecs - waterlevel = storage / lake.area[i] + waterlevel = storage / lake_p.area[i] end ### Linearisation for specific storage/rating curves ### - if lake.outflowfunc[i] == 1 || lake.outflowfunc[i] == 2 - diff_wl = has_lowerlake ? lake.waterlevel[i] - lake.waterlevel[lo] : 0.0 - - storage_input = (lake.storage[i] + precipitation - actevap) / timestepsecs + inflow - if lake.outflowfunc[i] == 1 - outflow = - interpolate_linear(lake.waterlevel[i], lake.hq[i].H, lake.hq[i].Q[:, doy]) + if lake_p.outflowfunc[i] == 1 || lake_p.outflowfunc[i] == 2 + diff_wl = has_lowerlake ? lake_v.waterlevel[i] - lake_v.waterlevel[lo] : 0.0 + + storage_input = + (lake_v.storage[i] + precipitation - actevap) / timestepsecs + inflow + if lake_p.outflowfunc[i] == 1 + outflow = interpolate_linear( + lake_v.waterlevel[i], + lake_p.hq[i].H, + lake_p.hq[i].Q[:, doy], + ) outflow = min(outflow, storage_input) else if diff_wl >= 0.0 - if lake.waterlevel[i] > lake.threshold[i] - dh = lake.waterlevel[i] - lake.threshold[i] - outflow = lake.b[i] * pow(dh, lake.e[i]) - maxflow = (dh * lake.area[i]) / timestepsecs + if lake_v.waterlevel[i] > lake_p.threshold[i] + dh = lake_v.waterlevel[i] - lake_p.threshold[i] + outflow = lake_p.b[i] * pow(dh, lake_p.e[i]) + maxflow = (dh * lake_p.area[i]) / timestepsecs outflow = min(outflow, maxflow) else outflow = Float(0) end else - if lake.waterlevel[lo] > lake.threshold[i] - dh = lake.waterlevel[lo] - lake.threshold[i] - outflow = -1.0 * lake.b[i] * pow(dh, lake.e[i]) - maxflow = (dh * lake.area[lo]) / timestepsecs + if lake_v.waterlevel[lo] > lake_p.threshold[i] + dh = lake_v.waterlevel[lo] - lake_p.threshold[i] + outflow = -1.0 * lake_p.b[i] * pow(dh, lake_p.e[i]) + maxflow = (dh * lake_p.area[lo]) / timestepsecs outflow = max(outflow, -maxflow) else outflow = Float(0) @@ -559,43 +634,44 @@ function update!(lake::Lake, i, inflow, doy, timestepsecs) storage = (storage_input - outflow) * timestepsecs # update storage and outflow for lake with rating curve of type 1. - if lake.outflowfunc[i] == 1 - overflow = max(0.0, (storage - lake.maxstorage[i]) / timestepsecs) - storage = min(storage, lake.maxstorage[i]) + if lake_p.outflowfunc[i] == 1 + overflow = max(0.0, (storage - lake_p.maxstorage[i]) / timestepsecs) + storage = min(storage, lake_p.maxstorage[i]) outflow = outflow + overflow end - waterlevel = if lake.storfunc[i] == 1 - lake.waterlevel[i] + (storage - lake.storage[i]) / lake.area[i] + waterlevel = if lake_p.storfunc[i] == 1 + lake_v.waterlevel[i] + (storage - lake_v.storage[i]) / lake_p.area[i] else - interpolate_linear(storage, lake.sh[i].S, lake.sh[i].H) + interpolate_linear(storage, lake_p.sh[i].S, lake_p.sh[i].H) end # update lower lake (linked lakes) in case flow from lower lake to upper lake occurs if diff_wl < 0.0 - lowerlake_storage = lake.storage[lo] + outflow * timestepsecs + lowerlake_storage = lake_v.storage[lo] + outflow * timestepsecs - lowerlake_waterlevel = if lake.storfunc[lo] == 1 - lake.waterlevel[lo] + (lowerlake_storage - lake.storage[lo]) / lake.area[lo] + lowerlake_waterlevel = if lake_p.storfunc[lo] == 1 + lake_v.waterlevel[lo] + + (lowerlake_storage - lake_v.storage[lo]) / lake_p.area[lo] else - interpolate_linear(lowerlake_storage, lake.sh[lo].S, lake.sh[lo].H) + interpolate_linear(lowerlake_storage, lake_p.sh[lo].S, lake_p.sh[lo].H) end # update values for the lower lake in place - lake.outflow[lo] = -outflow - lake.totaloutflow[lo] += -outflow * timestepsecs - lake.storage[lo] = lowerlake_storage - lake.waterlevel[lo] = lowerlake_waterlevel + lake_v.outflow[lo] = -outflow + lake_v.totaloutflow[lo] += -outflow * timestepsecs + lake_v.storage[lo] = lowerlake_storage + lake_v.waterlevel[lo] = lowerlake_waterlevel end end # update values in place - lake.outflow[i] = max(outflow, 0.0) # for a linked lake flow can be negative - lake.waterlevel[i] = waterlevel - lake.inflow[i] += inflow * timestepsecs - lake.totaloutflow[i] += outflow * timestepsecs - lake.storage[i] = storage - lake.actevap[i] += 1000.0 * (actevap / lake.area[i]) + lake_v.outflow[i] = max(outflow, 0.0) # for a linked lake flow can be negative + lake_v.waterlevel[i] = waterlevel + lake_bc.inflow[i] += inflow * timestepsecs + lake_v.totaloutflow[i] += outflow * timestepsecs + lake_v.storage[i] = storage + lake_v.actevap[i] += 1000.0 * (actevap / lake_p.area[i]) return nothing end diff --git a/src/sbm_model.jl b/src/sbm_model.jl index f06be6234..404122042 100644 --- a/src/sbm_model.jl +++ b/src/sbm_model.jl @@ -77,7 +77,7 @@ function initialize_sbm_model(config::Config) pits = zeros(Bool, modelsize_2d) if do_reservoirs reservoirs, resindex, reservoir, pits = - initialize_simple_reservoir(config, nc, inds_riv, nriv, pits, tosecond(dt)) + SimpleReservoir(config, nc, inds_riv, nriv, pits, tosecond(dt)) else reservoir = () reservoirs = nothing @@ -86,8 +86,7 @@ function initialize_sbm_model(config::Config) # lakes if do_lakes - lakes, lakeindex, lake, pits = - initialize_lake(config, nc, inds_riv, nriv, pits, tosecond(dt)) + lakes, lakeindex, lake, pits = Lake(config, nc, inds_riv, nriv, pits, tosecond(dt)) else lake = () lakes = nothing diff --git a/src/states.jl b/src/states.jl index 72c8f5886..13e598e90 100644 --- a/src/states.jl +++ b/src/states.jl @@ -241,13 +241,13 @@ function extract_required_states(config::Config) # Add lake states to dict required_states = add_to_required_states( required_states, - (:lateral, :river, :boundary_conditions, :lake), + (:lateral, :river, :boundary_conditions, :lake, :variables), lake_states, ) # Add reservoir states to dict required_states = add_to_required_states( required_states, - (:lateral, :river, :boundary_conditions, :reservoir), + (:lateral, :river, :boundary_conditions, :reservoir, :variables), reservoir_states, ) # Add paddy states to dict diff --git a/test/bmi.jl b/test/bmi.jl index 01adf09b7..aa7719c77 100644 --- a/test/bmi.jl +++ b/test/bmi.jl @@ -23,7 +23,7 @@ tomlpath = joinpath(@__DIR__, "sbm_config.toml") "vertical.soil.parameters.nlayers", "vertical.soil.parameters.theta_r", "lateral.river.variables.q", - "lateral.river.boundary_conditions.reservoir.outflow", + "lateral.river.boundary_conditions.reservoir.variables.outflow", ] retrieved_vars = BMI.get_input_var_names(model) @test all(x -> x in retrieved_vars, to_check) @@ -36,7 +36,7 @@ tomlpath = joinpath(@__DIR__, "sbm_config.toml") @test BMI.get_var_grid(model, "lateral.river.variables.h") == 3 @test BMI.get_var_grid( model, - "lateral.river.boundary_conditions.reservoir.inflow", + "lateral.river.boundary_conditions.reservoir.boundary_conditions.inflow", ) == 0 @test_throws ErrorException BMI.get_var_grid( model, @@ -44,7 +44,7 @@ tomlpath = joinpath(@__DIR__, "sbm_config.toml") ) @test BMI.get_var_type( model, - "lateral.river.boundary_conditions.reservoir.inflow", + "lateral.river.boundary_conditions.reservoir.boundary_conditions.inflow", ) == "$Float" @test BMI.get_var_units(model, "vertical.soil.parameters.theta_s") == "-" @test BMI.get_var_itemsize(model, "lateral.subsurface.variables.ssf") == diff --git a/test/io.jl b/test/io.jl index ed13d4047..aa9b0e37d 100644 --- a/test/io.jl +++ b/test/io.jl @@ -213,8 +213,10 @@ end @test Wflow.param(model, "lateral.doesnt_exist", -1) == -1 @testset "warm states" begin - @test Wflow.param(model, "lateral.river.boundary_conditions.reservoir.volume")[1] ≈ - 3.2807224993363418e7 + @test Wflow.param( + model, + "lateral.river.boundary_conditions.reservoir.variables.volume", + )[1] ≈ 3.2807224993363418e7 @test Wflow.param(model, "vertical.soil.variables.satwaterdepth")[9115] ≈ 477.13548089422125 @test Wflow.param(model, "vertical.snow.variables.snow_storage")[5] ≈ 11.019233179897599 @@ -458,7 +460,10 @@ end @test (:lateral, :river, :variables, :q) in required_states @test (:lateral, :river, :variables, :h_av) in required_states @test (:lateral, :land, :variables, :h_av) in required_states - @test !((:lateral, :river, :boundary_conditions, :lake, :waterlevel) in required_states) + @test !( + (:lateral, :river, :boundary_conditions, :lake, :variables, :waterlevel) in + required_states + ) # Adding an unused state the see if the right warning message is thrown config.state.vertical.soil.variables.additional_state = "additional_state" diff --git a/test/reservoir_lake.jl b/test/reservoir_lake.jl index 261021927..02d569341 100644 --- a/test/reservoir_lake.jl +++ b/test/reservoir_lake.jl @@ -1,34 +1,43 @@ -res = Wflow.SimpleReservoir{Float64}(; + +res_bc = + Wflow.ReservoirBC{Float}(; inflow = [0.0], precipitation = [4.2], evaporation = [1.5]) +res_params = Wflow.ReservoirParameters{Float64}(; dt = 86400.0, demand = [52.523], maxrelease = [420.184], maxvolume = [25_000_000.0], - volume = [1.925e7], - totaloutflow = [0.0], - inflow = [0.0], area = [1885665.353626924], targetfullfrac = [0.8], targetminfrac = [0.2425554726620697], - precipitation = [4.2], - evaporation = [1.5], +) +res_vars = Wflow.ReservoirVariables{Float64}(; + volume = [1.925e7], + totaloutflow = [0.0], actevap = [0.0], outflow = [NaN], percfull = [NaN], demandrelease = [NaN], ) + +res = Wflow.SimpleReservoir{Float64}(; + boundary_conditions = res_bc, + parameters = res_params, + variables = res_vars, +) @testset "Update reservoir simple" begin Wflow.update!(res, 1, 100.0, 86400.0) - @test res.outflow[1] ≈ 91.3783714867453 - @test res.totaloutflow[1] ≈ 7.895091296454794e6 - @test res.volume[1] ≈ 2.0e7 - @test res.percfull[1] ≈ 0.80 - @test res.demandrelease[1] ≈ 52.5229994727611 - @test res.precipitation[1] ≈ 4.2 - @test res.evaporation[1] ≈ 1.5 - @test res.actevap[1] ≈ 1.5 + @test res.variables.outflow[1] ≈ 91.3783714867453 + @test res.variables.totaloutflow[1] ≈ 7.895091296454794e6 + @test res.variables.volume[1] ≈ 2.0e7 + @test res.variables.percfull[1] ≈ 0.80 + @test res.variables.demandrelease[1] ≈ 52.5229994727611 + @test res.boundary_conditions.precipitation[1] ≈ 4.2 + @test res.boundary_conditions.evaporation[1] ≈ 1.5 + @test res.variables.actevap[1] ≈ 1.5 end -lake = Wflow.Lake{Float64}(; +lake_bc = Wflow.LakeBC{Float}(; inflow = [0.0], precipitation = [20.0], evaporation = [3.2]) +lake_params = Wflow.LakeParameters{Float}(; dt = 86400.0, lowerlake_ind = [0], area = [180510409.0], @@ -36,30 +45,38 @@ lake = Wflow.Lake{Float64}(; threshold = [0.0], storfunc = [1], outflowfunc = [3], - totaloutflow = [0.0], - inflow = [0.0], b = [0.22], e = [2.0], sh = [missing], hq = [missing], +) +lake_vars = Wflow.LakeVariables{Float}(; + totaloutflow = [0.0], storage = Wflow.initialize_storage([1], [180510409.0], [18.5], [missing]), waterlevel = [18.5], - precipitation = [20.0], - evaporation = [3.2], actevap = [0.0], outflow = [NaN], ) + +lake = Wflow.Lake{Float64}(; + boundary_conditions = lake_bc, + parameters = lake_params, + variables = lake_vars, +) @testset "Update lake" begin + lake_p = lake.parameters + lake_v = lake.variables + lake_bc = lake.boundary_conditions Wflow.update!(lake, 1, 2500.0, 181, 86400.0) - @test Wflow.waterlevel(lake.storfunc, lake.area, lake.storage, lake.sh)[1] ≈ + @test Wflow.waterlevel(lake_p.storfunc, lake_p.area, lake_v.storage, lake_p.sh)[1] ≈ 19.672653848925634 - @test lake.outflow[1] ≈ 85.14292808113598 - @test lake.totaloutflow[1] ≈ 7.356348986210149e6 - @test lake.storage[1] ≈ 3.55111879238499e9 - @test lake.waterlevel[1] ≈ 19.672653848925634 - @test lake.precipitation[1] ≈ 20.0 - @test lake.evaporation[1] ≈ 3.2 - @test lake.actevap[1] ≈ 3.2 + @test lake_v.outflow[1] ≈ 85.14292808113598 + @test lake_v.totaloutflow[1] ≈ 7.356348986210149e6 + @test lake_v.storage[1] ≈ 3.55111879238499e9 + @test lake_v.waterlevel[1] ≈ 19.672653848925634 + @test lake_bc.precipitation[1] ≈ 20.0 + @test lake_bc.evaporation[1] ≈ 3.2 + @test lake_v.actevap[1] ≈ 3.2 end datadir = joinpath(@__DIR__, "data") @@ -71,7 +88,7 @@ sh = [ @test keys(sh[1]) == (:H, :S) @test typeof(values(sh[1])) == Tuple{Vector{Float}, Vector{Float}} - lake = Wflow.Lake{Float}(; + lake_params = Wflow.LakeParameters{Float}(; dt = 86400.0, lowerlake_ind = [2, 0], area = [472461536.0, 60851088.0], @@ -84,16 +101,15 @@ sh = [ ), threshold = [393.7, 0.0], storfunc = [2, 2], - inflow = [0.0, 0.0], - totaloutflow = [0.0, 0.0], outflowfunc = [2, 1], b = [140.0, 0.0], e = [1.5, 1.5], sh = sh, hq = [missing, Wflow.read_hq_csv(joinpath(datadir, "input", "lake_hq_2.csv"))], + ) + lake_vars = Wflow.LakeVariables{Float}(; + totaloutflow = [0.0, 0.0], waterlevel = [395.03027, 394.87833], - precipitation = [10.0, 10.0], - evaporation = [2.0, 2.0], actevap = [0.0, 0.0], outflow = [NaN, NaN], storage = Wflow.initialize_storage( @@ -103,27 +119,42 @@ sh = [ sh, ), ) + lake_bc = Wflow.LakeBC{Float}(; + inflow = [0.0, 0.0], + precipitation = [10.0, 10.0], + evaporation = [2.0, 2.0], + ) + + lake = Wflow.Lake{Float}(; + boundary_conditions = lake_bc, + parameters = lake_params, + variables = lake_vars, + ) Wflow.update!(lake, 1, 500.0, 15, 86400.0) Wflow.update!(lake, 2, 500.0, 15, 86400.0) - @test lake.outflow ≈ [214.80170846121263, 236.83281600000214] atol = 1e-2 - @test lake.totaloutflow ≈ [1.855886761104877e7, 2.0462355302400187e7] atol = 1e3 - @test lake.storage ≈ [1.2737435094769483e9, 2.6019755340159863e8] atol = 1e4 - @test lake.waterlevel ≈ [395.0912274997361, 395.2101079057371] atol = 1e-2 - lake.actevap .= 0.0 - lake.totaloutflow .= 0.0 - lake.inflow .= 0.0 + lake_v = lake.variables + lake_bc = lake.boundary_conditions + @test lake_v.outflow ≈ [214.80170846121263, 236.83281600000214] atol = 1e-2 + @test lake_v.totaloutflow ≈ [1.855886761104877e7, 2.0462355302400187e7] atol = 1e3 + @test lake_v.storage ≈ [1.2737435094769483e9, 2.6019755340159863e8] atol = 1e4 + @test lake_v.waterlevel ≈ [395.0912274997361, 395.2101079057371] atol = 1e-2 + lake_v.actevap .= 0.0 + lake_v.totaloutflow .= 0.0 + lake_bc.inflow .= 0.0 Wflow.update!(lake, 1, 500.0, 15, 86400.0) Wflow.update!(lake, 2, 500.0, 15, 86400.0) - @test lake.outflow ≈ [0.0, 239.66710359986183] atol = 1e-2 - @test lake.totaloutflow ≈ [-2.2446764487487033e7, 4.3154002238515094e7] atol = 1e3 - @test lake.storage ≈ [1.3431699662524352e9, 2.6073035986708355e8] atol = 1e4 - @test lake.waterlevel ≈ [395.239782021054, 395.21771942667266] atol = 1e-2 - @test lake.actevap ≈ [2.0, 2.0] + @test lake_v.outflow ≈ [0.0, 239.66710359986183] atol = 1e-2 + @test lake_v.totaloutflow ≈ [-2.2446764487487033e7, 4.3154002238515094e7] atol = 1e3 + @test lake_v.storage ≈ [1.3431699662524352e9, 2.6073035986708355e8] atol = 1e4 + @test lake_v.waterlevel ≈ [395.239782021054, 395.21771942667266] atol = 1e-2 + @test lake_v.actevap ≈ [2.0, 2.0] end @testset "overflowing lake with sh and hq" begin - lake = Wflow.Lake{Float}(; + lake_bc = + Wflow.LakeBC{Float}(; inflow = [0.00], precipitation = [10.0], evaporation = [2.0]) + lake_params = Wflow.LakeParameters{Float}(; dt = 86400.0, lowerlake_ind = [0], area = [200_000_000], @@ -136,26 +167,32 @@ end ), threshold = [0.0], storfunc = [2], - inflow = [0.00], - totaloutflow = [0.0], outflowfunc = [1], b = [0.0], e = [0.0], sh = [Wflow.read_sh_csv(joinpath(datadir, "input", "lake_sh_2.csv"))], hq = [Wflow.read_hq_csv(joinpath(datadir, "input", "lake_hq_2.csv"))], + ) + lake_vars = Wflow.LakeVariables{Float}(; + totaloutflow = [0.0], waterlevel = [397.75], - precipitation = [10.0], - evaporation = [2.0], actevap = [0.0], outflow = [NaN], storage = [410_760_000], ) + lake = Wflow.Lake{Float}(; + boundary_conditions = lake_bc, + parameters = lake_params, + variables = lake_vars, + ) Wflow.update!(lake, 1, 1500.0, 15, 86400.0) - @test Wflow.waterlevel(lake.storfunc, lake.area, lake.storage, lake.sh) ≈ [398.0] atol = - 1e-2 - @test lake.outflow ≈ [1303.67476852] atol = 1e-2 - @test lake.totaloutflow ≈ [11.26375000e7] atol = 1e3 - @test lake.storage ≈ [4.293225e8] atol = 1e4 - @test lake.waterlevel ≈ [398.000000] atol = 1e-2 + lake_p = lake.parameters + lake_v = lake.variables + @test Wflow.waterlevel(lake_p.storfunc, lake_p.area, lake_v.storage, lake_p.sh) ≈ + [398.0] atol = 1e-2 + @test lake_v.outflow ≈ [1303.67476852] atol = 1e-2 + @test lake_v.totaloutflow ≈ [11.26375000e7] atol = 1e3 + @test lake_v.storage ≈ [4.293225e8] atol = 1e4 + @test lake_v.waterlevel ≈ [398.000000] atol = 1e-2 end diff --git a/test/run_sbm.jl b/test/run_sbm.jl index b882597a2..8a1df28b2 100644 --- a/test/run_sbm.jl +++ b/test/run_sbm.jl @@ -132,18 +132,19 @@ end @testset "reservoir simple" begin res = model.lateral.river.boundary_conditions.reservoir - @test res.outflow[1] ≈ 0.21750000119148086f0 - @test res.inflow[1] ≈ 43.18479982574888f0 - @test res.volume[1] ≈ 2.751299001489657f7 - @test res.precipitation[1] ≈ 0.17999997735023499f0 - @test res.evaporation[1] ≈ 0.5400000810623169f0 + @test res.variables.outflow[1] ≈ 0.21750000119148086f0 + @test res.boundary_conditions.inflow[1] ≈ 43.18479982574888f0 + @test res.variables.volume[1] ≈ 2.751299001489657f7 + @test res.boundary_conditions.precipitation[1] ≈ 0.17999997735023499f0 + @test res.boundary_conditions.evaporation[1] ≈ 0.5400000810623169f0 end # set these variables for comparison in "changed dynamic parameters" precip = copy(model.vertical.atmospheric_forcing.precipitation) evap = copy(model.vertical.atmospheric_forcing.potential_evaporation) lai = copy(model.vertical.vegetation_parameter_set.leaf_area_index) -res_evap = copy(model.lateral.river.boundary_conditions.reservoir.evaporation) +res_evap = + copy(model.lateral.river.boundary_conditions.reservoir.boundary_conditions.evaporation) Wflow.close_files(model; delete_output = false) @@ -198,7 +199,8 @@ Wflow.run_timestep!(model) @test (vertical.atmospheric_forcing.potential_evaporation[100] - 1.50) / evap[100] ≈ 3.0f0 @test vertical.vegetation_parameter_set.leaf_area_index[100] / lai[100] ≈ 1.6f0 - @test (res.evaporation[2] - 1.50) / res_evap[2] ≈ 3.0000012203408635f0 + @test (res.boundary_conditions.evaporation[2] - 1.50) / res_evap[2] ≈ + 3.0000012203408635f0 end # test cyclic river inflow @@ -230,7 +232,10 @@ Wflow.load_fixed_forcing!(model) @test maximum(model.vertical.atmospheric_forcing.precipitation) ≈ 2.5 @test minimum(model.vertical.atmospheric_forcing.precipitation) ≈ 0.0 @test all( - isapprox.(model.lateral.river.boundary_conditions.reservoir.precipitation, 2.5), + isapprox.( + model.lateral.river.boundary_conditions.reservoir.boundary_conditions.precipitation, + 2.5, + ), ) end @@ -240,7 +245,10 @@ Wflow.run_timestep!(model) @test maximum(model.vertical.atmospheric_forcing.precipitation) ≈ 2.5 @test minimum(model.vertical.atmospheric_forcing.precipitation) ≈ 0.0 @test all( - isapprox.(model.lateral.river.boundary_conditions.reservoir.precipitation, 2.5), + isapprox.( + model.lateral.river.boundary_conditions.reservoir.boundary_conditions.precipitation, + 2.5, + ), ) end diff --git a/test/sbm_config.toml b/test/sbm_config.toml index d28b87705..9e805a61f 100644 --- a/test/sbm_config.toml +++ b/test/sbm_config.toml @@ -36,7 +36,7 @@ h = "h_river" h_av = "h_av_river" q = "q_river" -[state.lateral.river.boundary_conditions.reservoir] +[state.lateral.river.boundary_conditions.reservoir.variables] volume = "volume_reservoir" [state.lateral.subsurface.variables] @@ -165,7 +165,7 @@ snow_water = "snowwater" h = "h_river" q = "q_river" -[output.lateral.river.boundary_conditions.reservoir] +[output.lateral.river.boundary_conditions.reservoir.variables] volume = "volume_reservoir" [output.lateral.subsurface.variables] @@ -208,7 +208,7 @@ reducer = "maximum" [[csv.column]] header = "volume" index = 1 -parameter = "lateral.river.boundary_conditions.reservoir.volume" +parameter = "lateral.river.boundary_conditions.reservoir.variables.volume" [[csv.column]] coordinate.x = 6.255 @@ -267,5 +267,7 @@ components = [ "lateral.river.parameters", "lateral.river.parameters.flow", "lateral.river.boundary_conditions", - "lateral.river.boundary_conditions.reservoir", + "lateral.river.boundary_conditions.reservoir.boundary_conditions", + "lateral.river.boundary_conditions.reservoir.parameters", + "lateral.river.boundary_conditions.reservoir.variables", ]