From 775a575f4dfe42fdedc3748fde9618fd117d6b8e Mon Sep 17 00:00:00 2001 From: Takuya Iwanaga Date: Sat, 30 Nov 2024 15:58:49 +1100 Subject: [PATCH] Performance changes for all supported models --- src/GR4JNode.jl | 104 ++++++++++--------- src/HyModNode.jl | 46 +++------ src/IHACRESNode.jl | 244 +++++++++++++++++++-------------------------- src/SYMHYDNode.jl | 112 ++++++++++++--------- src/Streamfall.jl | 15 ++- 5 files changed, 248 insertions(+), 273 deletions(-) diff --git a/src/GR4JNode.jl b/src/GR4JNode.jl index eea236a..39af19f 100644 --- a/src/GR4JNode.jl +++ b/src/GR4JNode.jl @@ -1,5 +1,5 @@ """ -Portions of the GR4J implementation is directly adapted from Python code written by Andrew MacDonald (2014). +Portions of the GR4J implementation is adapted from Python code written by Andrew MacDonald (2014). Per license requirements, the full license conditions included below. @@ -65,8 +65,9 @@ A four-parameter model with two stores. https://github.com/amacd31/gr4j """ -Base.@kwdef mutable struct GR4JNode{P, A<:Real} <: GRNJNode - @network_node +Base.@kwdef mutable struct GR4JNode{P, A<:AbstractFloat} <: GRNJNode + const name::String + const area::A # parameters X1::P = Param(350.0, bounds=(1.0, 1500.0)) @@ -80,21 +81,37 @@ Base.@kwdef mutable struct GR4JNode{P, A<:Real} <: GRNJNode # x4 : time base of unit hydrograph UH1 (days, > 0.5) # stores - p_store::Array{A} = [0.0] - r_store::Array{A} = [0.0] + p_store::Vector{A} = [0.0] + r_store::Vector{A} = [0.0] - UH1::Array{Array{A}} = [] - UH2::Array{Array{A}} = [] + UH1::Vector{Vector{A}} = [] + UH2::Vector{Vector{A}} = [] - outflow::Array{A} = [] + uh1_ordinates::Vector{A} = [] + uh2_ordinates::Vector{A} = [] + + outflow::Vector{A} = [] end -# function prep_state!(node::HyModNode, sim_length::Int64) -# node.UH1 = fill(undef, sim_length) -# node.UH2 = fill(undef, sim_length) -# node.outflow = fill(undef, sim_length) -# end +function prep_state!(node::GR4JNode, timesteps::Int64) + resize!(node.p_store, timesteps+1) + resize!(node.r_store, timesteps+1) + node.p_store[2:end] .= 0.0 + node.r_store[2:end] .= 0.0 + + # Prep cache + X4 = node.X4.val + nUH1 = Int(ceil(X4)) + nUH2 = Int(ceil(2.0*X4)) + cUH1, cUH2 = fill(0.0, nUH1), fill(0.0, nUH2) + node.UH1 = fill(cUH1, timesteps) + node.UH2 = fill(cUH2, timesteps) + node.uh1_ordinates = Vector{Float64}(undef, nUH1) + node.uh2_ordinates = Vector{Float64}(undef, nUH2) + + node.outflow = fill(0.0, timesteps) +end function GR4JNode(name::String, spec::Dict) @@ -129,22 +146,6 @@ function GR4JNode(name::String, spec::Dict) end -""" - run_node!(node::GR4JNode, climate::Climate) - -Run GR4J node for all time steps in given climate sequence. -""" -function run_node!(node::GR4JNode, climate::Climate; inflow=nothing, extraction=nothing, exchange=nothing) - timesteps = sim_length(climate) - # prep_state!(node, sim_length) - for ts in 1:timesteps - run_timestep!(node, climate, ts; inflow=inflow, extraction=extraction, exchange=exchange) - end - - return node.outflow -end - - """ run_timestep!(node::GR4JNode, climate::Climate, timestep::Int; inflow::Float64, extraction::Float64, exchange::Float64) @@ -158,8 +159,10 @@ function run_timestep!(node::GR4JNode, climate::Climate, timestep::Int; inflow=n end -function run_timestep!(node::GR4JNode, rain, et, ts; inflow=nothing, extraction=nothing, exchange=nothing) - res = run_gr4j(rain, et, node.X1, node.X2, node.X3, node.X4, node.area, node.p_store[end], node.r_store[end]) +function run_timestep!(node::GR4JNode, rain::Float64, et::Float64, ts::Int64; inflow=nothing, extraction=nothing, exchange=nothing) + uh1_cache = node.uh1_ordinates + uh2_cache = node.uh2_ordinates + res = run_gr4j(rain, et, node.X1.val, node.X2.val, node.X3.val, node.X4.val, node.area, node.UH1[ts], node.UH2[ts], uh1_cache, uh2_cache, node.p_store[ts], node.r_store[ts]) Q, p_s, r_s, UH1, UH2 = res node_name = node.name @@ -171,7 +174,7 @@ function run_timestep!(node::GR4JNode, rain, et, ts; inflow=nothing, extraction= Q = Q + in_flow + ex - wo end - update_state!(node, p_s, r_s, Q, UH1, UH2) + update_state!(node, ts, p_s, r_s, Q, UH1, UH2) end @@ -182,6 +185,16 @@ function update_state!(node::GR4JNode, ps, rs, q, UH1, UH2) push!(node.UH1, UH1) push!(node.UH2, UH2) end +function update_state!(node::GR4JNode, ts::Int64, ps, rs, q, UH1, UH2)::Nothing + node.p_store[ts+1] = ps + node.r_store[ts+1] = rs + node.outflow[ts] = q + + node.UH1[ts] = UH1 + node.UH2[ts] = UH2 + + return nothing +end function update_params!(node::GR4JNode, X1::Float64, X2::Float64, X3::Float64, X4::Float64)::Nothing @@ -218,12 +231,12 @@ end Determine unit hydrograph ordinates. """ -function s_curve(t, x4; uh2::Bool = false)::Float64 +function s_curve(t::F, x4::F; uh2::Bool = false)::F where {F<:Float64} if t <= 0.0 return 0.0 end - ordinate::Float64 = 0.0 + ordinate::F = 0.0 if t < x4 ordinate = (t/x4)^2.5 else @@ -258,14 +271,10 @@ Generated simulated streamflow for given rainfall and potential evaporation. # Returns - tuple of simulated outflow [ML/day], and intermediate states: p_store, r_store, UH1, UH2 """ -function run_gr4j(P, E, X1, X2, X3, X4, area, p_store=0.0, r_store=0.0)::Tuple - nUH1 = Int(ceil(X4)) - nUH2 = Int(ceil(2.0*X4)) - - uh1_ordinates = Array{Float64}(undef, nUH1) - uh2_ordinates = Array{Float64}(undef, nUH2) - - UH1, UH2 = fill(0.0, nUH1), fill(0.0, nUH2) +function run_gr4j(P::F, E::F, X1::F, X2::F, X3::F, X4::F, area::F, UH1::Vector{Float64}, UH2::Vector{Float64}, uh1_ordinates::Vector{Float64}, uh2_ordinates::Vector{Float64}, p_store=0.0, r_store=0.0)::Tuple where {F<:Float64} + # Main.@infiltrate + nUH1::Int64 = Int(ceil(X4)) + nUH2::Int64 = Int(ceil(2.0*X4)) @inbounds for t in 2:(nUH1+1) t_f = Float64(t) @@ -277,7 +286,7 @@ function run_gr4j(P, E, X1, X2, X3, X4, area, p_store=0.0, r_store=0.0)::Tuple uh2_ordinates[t - 1] = s_curve(t_f, X4, uh2=true) - s_curve(t_f-1.0, X4, uh2=true) end - Q = 0.0 + Q::F = 0.0 if P > E net_evap = 0.0 scaled_net_precip = min((P - E)/X1, 13.0) @@ -303,12 +312,11 @@ function run_gr4j(P, E, X1, X2, X3, X4, area, p_store=0.0, r_store=0.0)::Tuple p_store = p_store - net_evap + reservoir_production - tmp_a = (p_store / 2.25 / X1)^4 - tmp_b = (1 + tmp_a)^0.25 + tmp_a::F = (p_store / 2.25 / X1)^4 + tmp_b::F = (1 + tmp_a)^0.25 percolation = p_store / tmp_b routed_volume = routed_volume + (p_store-percolation) - p_store = percolation @inbounds for i in 1:nUH1-1 UH1[i] = UH1[i+1] + uh1_ordinates[i]*routed_volume @@ -321,7 +329,7 @@ function run_gr4j(P, E, X1, X2, X3, X4, area, p_store=0.0, r_store=0.0)::Tuple UH2[end] = uh2_ordinates[end] * routed_volume tmp_a = (r_store / X3)^3.5 - groundwater_exchange = X2 * tmp_a + groundwater_exchange::F = X2 * tmp_a r_store = max(0.0, r_store + UH1[1] * 0.9 + groundwater_exchange) tmp_a = (r_store / X3)^4 @@ -333,5 +341,5 @@ function run_gr4j(P, E, X1, X2, X3, X4, area, p_store=0.0, r_store=0.0)::Tuple Q = (QR + QD) * area - return Q, p_store, r_store, UH1, UH2 + return Q, percolation, r_store, UH1, UH2 end diff --git a/src/HyModNode.jl b/src/HyModNode.jl index a8b4eeb..c68cbff 100644 --- a/src/HyModNode.jl +++ b/src/HyModNode.jl @@ -17,8 +17,8 @@ Adapted with kind permission from: https://github.com/jdherman/GRA-2020-SALib https://doi.org/10.5194/hess-17-149-2013 """ Base.@kwdef mutable struct SimpleHyModNode{P, A<:AbstractFloat} <: HyModNode - name::String - area::A + const name::String + const area::A # parameters Sm_max::P = Param(250.0, bounds=(1.0, 500.0)) @@ -68,37 +68,23 @@ function SimpleHyModNode(name::String, area::Float64, sm_max::Float64, B::Float6 end -function prep_state!(node::HyModNode, sim_length::Int64) - node.Sm = fill(0.0, sim_length+1) - node.Sf1 = fill(0.0, sim_length+1) - node.Sf2 = fill(0.0, sim_length+1) - node.Sf3 = fill(0.0, sim_length+1) - node.Ss1 = fill(0.0, sim_length+1) - - node.outflow = fill(0.0, sim_length) -end +function prep_state!(node::SimpleHyModNode, sim_length::Int64)::Nothing + resize!(node.Sm, sim_length+1) + resize!(node.Sf1, sim_length+1) + resize!(node.Sf2, sim_length+1) + resize!(node.Sf3, sim_length+1) + resize!(node.Ss1, sim_length+1) + node.Sm[2:end] .= 0.0 + node.Sf1[2:end] .= 0.0 + node.Sf2[2:end] .= 0.0 + node.Sf3[2:end] .= 0.0 + node.Ss1[2:end] .= 0.0 -""" - run_node!(node::HyModNode, climate::Climate; - inflow=nothing, extraction=nothing, exchange=nothing) + resize!(node.outflow, sim_length) + node.outflow .= 0.0 -Run given HyMod node for entire simulation period. -""" -function run_node!( - node::HyModNode, climate::Climate; - inflow=nothing, extraction=nothing, exchange=nothing -) - timesteps = sim_length(climate) - prep_state!(node, timesteps) - - # Identify relevant climate data - P_and_ET = Matrix{Float64}(climate_values(node, climate)) - - for ts in 1:timesteps - run_timestep!(node, P_and_ET[ts, 1], P_and_ET[ts, 2], ts; - inflow=inflow, extraction=extraction, exchange=exchange) - end + return nothing end """ diff --git a/src/IHACRESNode.jl b/src/IHACRESNode.jl index 7f22234..d1cee59 100644 --- a/src/IHACRESNode.jl +++ b/src/IHACRESNode.jl @@ -5,8 +5,9 @@ using ModelParameters abstract type IHACRESNode <: NetworkNode end -Base.@kwdef mutable struct BilinearNode{P, A<:Real} <: IHACRESNode - @network_node +Base.@kwdef mutable struct BilinearNode{P, A<:AbstractFloat} <: IHACRESNode + const name::String + const area::A # https://wiki.ewater.org.au/display/SD41/IHACRES-CMD+-+SRG d::P = Param(200.0, bounds=(10.0, 550.0)) # flow threshold @@ -23,7 +24,7 @@ Base.@kwdef mutable struct BilinearNode{P, A<:Real} <: IHACRESNode storage_coef::P = Param(2.9, bounds=(1e-10, 10.0)) alpha::P = Param(0.1, bounds=(1e-5, 1 - 1/10^9)) - level_params::Array{P, 1} = [ + const level_params::Array{P, 1} = [ Param(-0.01, bounds=(-10.0, -0.01)), # p1 Param(0.8, bounds=(0.0, 1.5)), # p2 Param(4.5, bounds=(0.0, 20.0)), # p3 @@ -50,29 +51,27 @@ Base.@kwdef mutable struct BilinearNode{P, A<:Real} <: IHACRESNode end -# function prep_state!(node::IHACRESNode, sim_length::Int64) -# ini_storage = node.storage[1] -# node.storage = fill(0.0, sim_length+1) -# node.storage[1] = ini_storage +function prep_state!(node::IHACRESNode, timesteps::Int64)::Nothing + resize!(node.storage, timesteps+1) + node.storage[2:end] .= 0.0 -# ini_qs = node.quick_store[1] -# node.quick_store = fill(0.0, sim_length+1) -# node.quick_store[1] = ini_qs + resize!(node.quick_store, timesteps+1) + node.quick_store[2:end] .= 0.0 -# ini_ss = node.slow_store[1] -# node.slow_store = fill(0.0, sim_length+1) -# node.slow_store[1] = ini_ss + resize!(node.slow_store, timesteps+1) + node.slow_store[2:end] .= 0.0 -# node.outflow = fill(0.0, sim_length) -# node.effective_rainfall = fill(0.0, sim_length) -# node.et = fill(0.0, sim_length) -# node.inflow = fill(0.0, sim_length) -# node.level = fill(0.0, sim_length) + node.outflow = fill(0.0, timesteps) + node.effective_rainfall = fill(0.0, timesteps) + node.et = fill(0.0, timesteps) + node.inflow = fill(0.0, timesteps) + node.level = fill(0.0, timesteps) -# ini_gw = node.gw_store[1] -# node.gw_store = fill(0.0, sim_length+1) -# node.gw_store[1] = ini_gw -# end + resize!(node.gw_store, timesteps+1) + node.gw_store[2:end] .= 0.0 + + return nothing +end function BilinearNode(name::String, spec::Dict) @@ -144,9 +143,9 @@ function BilinearNode(name::String, area::Float64, d::Float64, d2::Float64, e::F end -function update_state(s_node::IHACRESNode, storage::Float64, e_rainfall::Float64, et::Float64, +function update_state!(s_node::IHACRESNode, storage::Float64, e_rainfall::Float64, et::Float64, qflow_store::Float64, sflow_store::Float64, outflow::Float64, - level::Float64, gw_store::Float64) + level::Float64, gw_store::Float64)::Nothing push!(s_node.storage, storage) push!(s_node.effective_rainfall, e_rainfall) push!(s_node.et, et) @@ -156,9 +155,36 @@ function update_state(s_node::IHACRESNode, storage::Float64, e_rainfall::Float64 push!(s_node.slow_store, sflow_store) push!(s_node.level, level) push!(s_node.gw_store, gw_store) -end + return nothing +end +function update_state!( + s_node::IHACRESNode, + ts::Int64, + storage::Float64, + e_rainfall::Float64, + et::Float64, + qflow_store::Float64, + sflow_store::Float64, + inflow::Float64, + outflow::Float64, + level::Float64, + gw_store::Float64 +)::Nothing + s_node.storage[ts+1] = storage + s_node.quick_store[ts+1] = qflow_store + s_node.slow_store[ts+1] = sflow_store + s_node.gw_store[ts+1] = gw_store + + s_node.inflow[ts] = inflow + s_node.outflow[ts] = outflow + + s_node.effective_rainfall[ts] = e_rainfall + s_node.et[ts] = et + s_node.level[ts] = level + return nothing +end """ run_node!(node::IHACRESNode, climate::Climate, timestep::Int; @@ -192,103 +218,65 @@ function run_node!(node::IHACRESNode, climate::Climate, ts::Int; return run_timestep!(node, rain, et, ts; inflow=in_flow, extraction=wo, exchange=ex) end - """ - run_step!(s_node::BilinearNode, - rain::Float64, - evap::Float64, - inflow::Float64, - ext::Float64, - gw_exchange::Float64; - current_store::Float64, - quick_store::Float64, - slow_store::Float64, - gw_store::Float64)::Tuple{Float64, Float64} - -Run node with ET data to calculate outflow and update state. + run_timestep!( + node::BilinearNode, climate::Climate, timestep::Int, + inflow::Float64, extraction::Float64, exchange::Float64 + ) -# Arguments -- `s_node::BilinearNode` : IHACRESNode -- rain : rainfall for time step -- evap : evapotranspiration for time step -- inflow : inflow from previous node -- ext : irrigation and other water extractions -- gw_exchange : flux in ML where positive is contribution to stream, negative is infiltration -- current_store : replacement cmd state value (uses last known state if not provided) -- quick_store : replacement quick store state value -- slow_store : replacement slow store state value -- gw_store : replacement groundwater store state value - -# Returns -- float, outflow from node [ML/day], stream level +Run the given BilinearNode for a time step. """ -function run_step!(s_node::IHACRESNode, - rain::Float64, - evap::Float64, - inflow::Float64, - ext::Float64, - gw_exchange::Float64, - current_store::Float64, - quick_store::Float64, - slow_store::Float64, - gw_store::Float64)::Tuple{Float64, Float64} - - s_node.storage[end] = current_store - s_node.quick_store[end] = quick_store - s_node.slow_store[end] = slow_store - s_node.gw_store[end] = gw_store - - return _run_step!(s_node, rain, evap, inflow, ext, gw_exchange, - current_store, quick_store, slow_store, gw_store) -end - - -function run_step!(s_node::IHACRESNode, - rain::Float64, - evap::Float64, - inflow::Float64, - ext::Float64, - gw_exchange::Float64)::Tuple{Float64, Float64} - - current_store = s_node.storage[end] - quick_store = s_node.quick_store[end] - slow_store = s_node.slow_store[end] - gw_store = s_node.gw_store[end] +function run_timestep!( + node::IHACRESNode, climate::Climate, timestep::Int; + inflow=nothing, extraction=nothing, exchange=nothing +)::Float64 + ts = timestep + P, ET = climate_values(node, climate, ts) - return _run_step!(s_node, rain, evap, inflow, ext, gw_exchange, - current_store, quick_store, slow_store, gw_store) + return run_timestep!(node, ts, P, ET, in_flow, ext, flux) end +function run_timestep!(node, P, ET, ts; inflow=nothing, extraction=nothing, exchange=nothing) + in_flow = timestep_value(ts, node.name, "inflow", inflow) + ext = timestep_value(ts, node.name, "extraction", extraction) + flux = timestep_value(ts, node.name, "exchange", exchange) + return run_ihacres!(node, P, ET, ts, in_flow, ext, flux) +end -function _run_step!(s_node::IHACRESNode, - rain::Float64, - evap::Float64, - inflow::Float64, - ext::Float64, - gw_exchange::Float64, - current_store::Float64, - quick_store::Float64, - slow_store::Float64, - gw_store::Float64)::Tuple{Float64, Float64} +function run_ihacres!( + s_node::IHACRESNode, + rain::F, + evap::F, + ts::Int64, + inflow::F, + ext::F, + gw_exchange::F +)::Tuple{F, F} where {F<:Float64} + current_store::F = s_node.storage[ts] + quick_store::F = s_node.quick_store[ts] + slow_store::F = s_node.slow_store[ts] + gw_store::F = s_node.gw_store[ts] cache_res = s_node.cache.i_cache - @ccall IHACRES.calc_ft_interim_cmd(cache_res::Ptr{Cdouble}, - current_store::Cdouble, - rain::Cdouble, - s_node.d::Cdouble, - s_node.d2::Cdouble, - s_node.alpha::Cdouble)::Cvoid + @ccall IHACRES.calc_ft_interim_cmd( + cache_res::Ptr{Cdouble}, + current_store::Cdouble, + rain::Cdouble, + s_node.d.val::Cdouble, + s_node.d2.val::Cdouble, + s_node.alpha.val::Cdouble + )::Cvoid (mf, e_rainfall, recharge) = cache_res - et::Float64 = @ccall IHACRES.calc_ET( - s_node.e::Cdouble, + et::F = @ccall IHACRES.calc_ET( + s_node.e.val::Cdouble, evap::Cdouble, mf::Cdouble, - s_node.f::Cdouble, - s_node.d::Cdouble + s_node.f.val::Cdouble, + s_node.d.val::Cdouble )::Cdouble - cmd::Float64 = @ccall IHACRES.calc_cmd( + cmd::F = @ccall IHACRES.calc_cmd( current_store::Cdouble, rain::Cdouble, et::Cdouble, @@ -303,8 +291,8 @@ function _run_step!(s_node::IHACRESNode, e_rainfall::Cdouble, recharge::Cdouble, s_node.area::Cdouble, - s_node.a::Cdouble, - s_node.b::Cdouble + s_node.a.val::Cdouble, + s_node.b.val::Cdouble )::Cvoid (nq_store, ns_store, outflow) = cache_res @@ -312,7 +300,7 @@ function _run_step!(s_node::IHACRESNode, @ccall IHACRES.routing( routing_res::Ptr{Cdouble}, gw_store::Cdouble, - s_node.storage_coef::Cdouble, + s_node.storage_coef.val::Cdouble, inflow::Cdouble, outflow::Cdouble, ext::Cdouble, @@ -321,50 +309,18 @@ function _run_step!(s_node::IHACRESNode, (gw_store, outflow) = routing_res - level_params = Array{Float64}(s_node.level_params) + level_params = Vector{Float64}(s_node.level_params) level::Float64 = @ccall IHACRES.calc_ft_level( outflow::Cdouble, level_params::Ptr{Cdouble} )::Cdouble - push!(s_node.inflow, inflow) - update_state(s_node, cmd, e_rainfall, et, nq_store, ns_store, outflow, level, gw_store) + update_state!(s_node, ts, cmd, e_rainfall, et, nq_store, ns_store, inflow, outflow, level, gw_store) return outflow, level end -""" - run_timestep!(s_node::IHACRESNode, - rain::Float64, - evap::Float64, - timestep::Union{Int64, Nothing}; - inflow::Float64, - extraction::Float64, - exchange::Float64)::Tuple{Float64, Float64} - -Run node for a given time step. -""" -function run_timestep!(s_node::IHACRESNode, - rain::Float64, - evap::Float64, - ts::Union{Int64, Nothing}; - inflow::Float64=0.0, - extraction::Float64=0.0, - exchange::Float64=0.0)::Tuple{Float64, Float64} - if !isnothing(ts) - current_store = s_node.storage[ts] - quick_store = s_node.quick_store[ts] - slow_store = s_node.slow_store[ts] - gw_store = s_node.gw_store[ts] - return run_step!(s_node, rain, evap, inflow, extraction, exchange, - current_store, quick_store, slow_store, gw_store) - end - - return run_step!(s_node, rain, evap, inflow, extraction, exchange) -end - - """ param_info(node::IHACRESNode; with_level::Bool = true)::Tuple diff --git a/src/SYMHYDNode.jl b/src/SYMHYDNode.jl index bcf23b5..15b7357 100644 --- a/src/SYMHYDNode.jl +++ b/src/SYMHYDNode.jl @@ -7,8 +7,9 @@ const SYMHYD_SOIL_ET_CONST = 10.0 """ """ -Base.@kwdef mutable struct SYMHYDNode{P, A<:Real} <: NetworkNode - @network_node +Base.@kwdef mutable struct SYMHYDNode{P, A<:AbstractFloat} <: NetworkNode + name::String + area::A # parameters baseflow_coef::P = Param(0.5, bounds=(0.0, 1.0)) @@ -69,18 +70,19 @@ function SYMHYDNode(name::String, spec::Dict) return n end +function prep_state!(node::SYMHYDNode, sim_length::Int64)::Nothing + resize!(node.sm_store, sim_length+1) + resize!(node.gw_store, sim_length+1) + resize!(node.total_store, sim_length+1) + node.sm_store[2:end] .= 0.0 + node.gw_store[2:end] .= 0.0 + node.total_store[2:end] .= 0.0 -""" - -Run SYMHYD for all time steps within a climate scenario. -""" -function run_node!(node::SYMHYDNode, climate::Climate; inflow=nothing, extraction=nothing, exchange=nothing) - timesteps = length(climate) - for ts in 1:timesteps - run_timestep!(node, climate, ts; inflow=inflow, extraction=extraction, exchange=exchange) - end + node.outflow = fill(0.0, sim_length) + node.baseflow = fill(0.0, sim_length) + node.quickflow = fill(0.0, sim_length) - return node.outflow + return nothing end @@ -88,16 +90,24 @@ end Run SYMHYD for a given time step """ -function run_timestep!(node::SYMHYDNode, climate::Climate, timestep::Int; inflow=nothing, extraction=extraction, exchange=nothing) - P, E = climate_values(node, climate, timestep) +function run_timestep!(node::SYMHYDNode, climate::Climate, ts::Int; inflow=nothing, extraction=extraction, exchange=nothing)::Nothing + P, E = climate_values(node, climate, ts) - run_timestep!(node, P, E, timestep; inflow=inflow, extraction=extraction, exchange=exchange) -end + run_timestep!(node, P, E, ts; inflow=inflow, extraction=extraction, exchange=exchange) + return nothing +end -function run_timestep!(node::SYMHYDNode, rain, et, ts; inflow=nothing, extraction=nothing, exchange=nothing) - res = run_symhyd(node, rain, et) - sm_store, gw_store, total_store, total_runoff, baseflow, event_runoff = res +function run_timestep!( + node::SYMHYDNode, + rain::F, + et::F, + ts::Int; + inflow=nothing, + extraction=nothing, + exchange=nothing +)::Nothing where {F<:AbstractFloat} + sm_store, gw_store, total_store, total_runoff, baseflow, event_runoff = run_symhyd(node, rain, et, ts) node_name = node.name wo = timestep_value(ts, node_name, "releases", extraction) @@ -108,8 +118,10 @@ function run_timestep!(node::SYMHYDNode, rain, et, ts; inflow=nothing, extractio total_runoff = total_runoff + in_flow + ex - wo end - area = node.area - update_state!(node, sm_store, gw_store, total_store, total_runoff * area, baseflow * area, event_runoff * area) + area::F = node.area + update_state!(node, ts, sm_store, gw_store, total_store, total_runoff * area, baseflow * area, event_runoff * area) + + return nothing end @@ -121,6 +133,14 @@ function update_state!(node::SYMHYDNode, sm_store, gw_store, total_store, outflo push!(node.baseflow, baseflow) push!(node.quickflow, quickflow) end +function update_state!(node::SYMHYDNode, ts::Int64, sm_store, gw_store, total_store, outflow, baseflow, quickflow) + node.sm_store[ts+1] = sm_store + node.gw_store[ts+1] = gw_store + node.total_store[ts+1] = total_store + node.outflow[ts] = outflow + node.baseflow[ts] = baseflow + node.quickflow[ts] = quickflow +end function update_params!(node::SYMHYDNode, baseflow_coef::Float64, impervious_threshold::Float64, @@ -169,33 +189,33 @@ end Run SYMHYD for a single time step with given inputs and state variables. """ -function run_symhyd(node::SYMHYDNode, P, ET) +function run_symhyd(node::SYMHYDNode, P::F, ET::F, ts::Int64)::NTuple{6,F} where {F<:Float64} - sm_store = node.sm_store[end] - gw_store = node.gw_store[end] - total_store = node.total_store[end] + sm_store::F = node.sm_store[ts] + gw_store::F = node.gw_store[ts] + total_store::F = node.total_store[ts] - pervious_incident = P - impervious_incident = P + pervious_incident::F = P + impervious_incident::F = P - impervious_ET = min(node.impervious_threshold, impervious_incident) - impervious_runoff = impervious_incident - impervious_ET + impervious_ET::F = min(node.impervious_threshold, impervious_incident) + impervious_runoff::F = impervious_incident - impervious_ET - interception_ET = min(pervious_incident, min(ET, node.risc)) - throughfall = pervious_incident - impervious_ET + interception_ET::F = min(pervious_incident, min(ET, node.risc)) + throughfall::F = pervious_incident - impervious_ET - smsc = node.smsc - sm_fraction = sm_store / smsc + smsc::F = node.smsc.val + sm_fraction::F = sm_store / smsc - infiltration_capacity = node.infiltration_coef * exp(-node.infiltration_shape*sm_fraction) - infiltration = min(throughfall, infiltration_capacity) - infiltration_Xs_runoff = throughfall - infiltration + infiltration_capacity::F = node.infiltration_coef.val::F * exp(-node.infiltration_shape.val::F * sm_fraction)::F + infiltration::F = min(throughfall, infiltration_capacity) + infiltration_Xs_runoff::F = throughfall - infiltration - interflow_runoff = node.interflow_coef * sm_fraction * infiltration - infiltration_after_interflow = infiltration - interflow_runoff + interflow_runoff::F = node.interflow_coef.val::F * sm_fraction * infiltration + infiltration_after_interflow::F = infiltration - interflow_runoff - recharge = node.recharge_coef * sm_fraction * infiltration_after_interflow - soil_input = infiltration_after_interflow - recharge + recharge::F = node.recharge_coef.val::F * sm_fraction * infiltration_after_interflow + soil_input::F = infiltration_after_interflow - recharge # Update states... sm_store += soil_input @@ -208,18 +228,18 @@ function run_symhyd(node::SYMHYDNode, P, ET) sm_fraction = 1.0 end - baseflow_runoff = node.baseflow_coef * gw_store + baseflow_runoff::F = node.baseflow_coef.val::F * gw_store gw_store -= baseflow_runoff - soil_ET = min(sm_store, min(ET - interception_ET, sm_fraction*SYMHYD_SOIL_ET_CONST)) + soil_ET::F = min(sm_store, min(ET - interception_ET, sm_fraction*SYMHYD_SOIL_ET_CONST)) sm_store -= soil_ET - pervious_frac = node.pervious_fraction + pervious_frac::F = node.pervious_fraction.val::F total_store = (sm_store + gw_store) * pervious_frac - event_runoff = (1.0 - pervious_frac) * impervious_runoff + pervious_frac * (infiltration_Xs_runoff + interflow_runoff) - total_runoff = event_runoff + pervious_frac * baseflow_runoff - baseflow = baseflow_runoff * pervious_frac + event_runoff::F = (1.0 - pervious_frac) * impervious_runoff + pervious_frac * (infiltration_Xs_runoff + interflow_runoff) + total_runoff::F = event_runoff + pervious_frac * baseflow_runoff + baseflow::F = baseflow_runoff * pervious_frac # values for time step... return sm_store, gw_store, total_store, total_runoff, baseflow, event_runoff diff --git a/src/Streamfall.jl b/src/Streamfall.jl index ed1a174..f5465b1 100644 --- a/src/Streamfall.jl +++ b/src/Streamfall.jl @@ -45,6 +45,7 @@ end include("Network.jl") include("Node.jl") include("Climate.jl") +# include("Operation.jl") include("IHACRESNode.jl") include("IHACRESExpuhNode.jl") include("GR4JNode.jl") @@ -222,22 +223,26 @@ Run a specific node, and only that node, for all time steps. - `extraction::Union{DataFrame, Vector, Number}` : Extractions from this subcatchment - `exchange::Union{DataFrame, Vector, Number}` : Groundwater flux """ -function run_node!(node::NetworkNode, climate::Climate; inflow=nothing, extraction=nothing, exchange=nothing) +function run_node!(node::NetworkNode, climate::Climate; inflow=nothing, extraction=nothing, exchange=nothing)::Nothing timesteps = sim_length(climate) + prep_state!(node, timesteps) + + P_and_ET = climate_values(node, climate) node_name = node.name @inbounds for ts in 1:timesteps - rain, et = climate_values(node, climate, ts) - # Check if values are provided, otherwise default to 0. ext = timestep_value(ts, node_name, "extraction", extraction) flux = timestep_value(ts, node_name, "exchange", exchange) in_flow = timestep_value(ts, node_name, "inflow", inflow) - run_timestep!(node, rain, et, ts; inflow=in_flow, extraction=ext, exchange=flux) + run_timestep!( + node, P_and_ET[ts, 1], P_and_ET[ts, 2], ts; + inflow=in_flow, extraction=ext, exchange=flux + ) end - return node.outflow, node.level + return nothing end