Check if all states are covered (#359)

* add glacierstore to statevars * added state checks removed statevars() functions; added all state-related function to states.jl * run JuliaFormatter * fix function for other model configurations also fix sbm+1d floodplain test to ensure required states are provided * put snow and glacier flags back * improve floodplain handling * added tests; update changelog * fix test_throws * fix review comments Co-Authored-By: Willem van Verseveld <[email protected]> * run JuliaFormatter --------- Co-authored-by: Willem van Verseveld <[email protected]>
Deltares · Feb 20, 2024 · 9037c05 · 9037c05
1 parent bdc662b
commit 9037c05
Show file tree

Hide file tree

Showing 24 changed files with 384 additions and 156 deletions.
diff --git a/build/create_binaries/add_metadata.jl b/build/create_binaries/add_metadata.jl
@@ -20,18 +20,18 @@ function add_metadata(project_dir, license_file, output_dir, git_repo)
     cp(
         normpath(project_dir, "Project.toml"),
         normpath(output_dir, "share/julia/Project.toml");
-        force=true
+        force = true,
     )
     # the Manifest.toml always gives the exact version of Wflow that was built
     cp(
         normpath(project_dir, "Manifest.toml"),
         normpath(output_dir, "share/julia/Manifest.toml");
-        force=true
+        force = true,
     )
 
     # put the LICENSE in the top level directory
-    cp(license_file, normpath(output_dir, "LICENSE"); force=true)
-    cp(normpath(project_dir, "README.md"), normpath(output_dir, "README.md"); force=true)
+    cp(license_file, normpath(output_dir, "LICENSE"); force = true)
+    cp(normpath(project_dir, "README.md"), normpath(output_dir, "README.md"); force = true)
     open(normpath(output_dir, "README.md"), "a") do io
         # since the exact Wflow version may be hard to find in the Manifest.toml file
         # we can also extract that information, and add it to the README.md

diff --git a/build/create_binaries/create_app.jl b/build/create_binaries/create_app.jl
@@ -15,10 +15,10 @@ create_app(
     project_dir,
     output_dir;
     # map from binary name to julia function name
-    executables=["wflow_cli" => "julia_main"],
-    precompile_execution_file="precompile.jl",
-    filter_stdlibs=false,
-    force=true,
+    executables = ["wflow_cli" => "julia_main"],
+    precompile_execution_file = "precompile.jl",
+    filter_stdlibs = false,
+    force = true,
 )
 
 include("add_metadata.jl")

diff --git a/docs/src/changelog.md b/docs/src/changelog.md
@@ -29,6 +29,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Total water storage as an export variable for `SBM` concept. This is the total water stored
   per grid cell in millimeters. Excluded from this variable are the floodplain, lakes and
   reservoirs.
+- Checks to see if all states are covered in the .toml file. If not all states are covered,
+  an error is thrown. If there are more states specified than required, these states are
+  ignored (with a warning in the logging), and the simulation will continue.
 
 ## v0.7.3 - 2024-01-12
 
@@ -47,7 +50,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   terms of the quadratic equation (and solution) were fixed.
 - Use `kvfrac` for the computation of vertical saturated hydraulic conductivity at the
   bottom of the soil layer, since `kvfrac` is also used for the computation of vertical
-  unsaturated flow. 
+  unsaturated flow.
 
 ### Changed
 - For cyclic parameters different cyclic time inputs are supported (only one common cyclic

diff --git a/src/Wflow.jl b/src/Wflow.jl
@@ -52,15 +52,15 @@ end
 
 function Clock(config, reader)
     nctimes = reader.dataset["time"][:]
-    
+
     # if the timestep is not given, use the difference between NetCDF time 1 and 2
     timestepsecs = get(config, "timestepsecs", nothing)
     if timestepsecs === nothing
         timestepsecs = Dates.value(Second(nctimes[2] - nctimes[1]))
         config.timestepsecs = timestepsecs
     end
     Δt = Second(timestepsecs)
-    
+
     # if the config file does not have a start or endtime, follow the NetCDF times
     # and add them to the config
     starttime = get(config, "starttime", nothing)
@@ -80,7 +80,7 @@ function Clock(config, reader)
         config.starttime = starttime + Δt
     end
     starttime = cftime(config.starttime, calendar)
-    
+
     Clock(starttime, 0, Δt)
 end
 
@@ -133,6 +133,7 @@ include("utils.jl")
 include("bmi.jl")
 include("subdomains.jl")
 include("logging.jl")
+include("states.jl")
 
 """
     run(tomlpath::AbstractString; silent=false)
@@ -218,8 +219,8 @@ function run(model::Model; close_files = true)
     # determine timesteps to run
     calendar = get(config, "calendar", "standard")::String
     @warn string(
-        "The definition of `starttime` has changed (equal to model state time).\n Please", 
-        " update your settings TOML file by subtracting one model timestep Δt from the", 
+        "The definition of `starttime` has changed (equal to model state time).\n Please",
+        " update your settings TOML file by subtracting one model timestep Δt from the",
         " `starttime`, if it was used with a Wflow version up to v0.6.3.",
     )
     starttime = clock.time

diff --git a/src/flextopo.jl b/src/flextopo.jl
@@ -6,17 +6,17 @@
     # Number of cells
     n::Int | "-" | 0 | "none" | "none"
     #dictionary with all possible functions for each store
-    dic_function::Dict |"-"| 0 | "none" | "none"
+    dic_function::Dict | "-" | 0 | "none" | "none"
     #current class
-    kclass::Vector{Int64} |"-"| 0 | "none" | "none"
-    classes::Vector{String} |"-"| 0 | "none" | "none"
-    select_snow::Vector{String} |"-"| 0 | "none" | "none"
-    select_interception::Vector{String} |"-"| 0 | "none" | "none"
-    select_hortonponding::Vector{String} |"-"| 0 | "none" | "none"
-    select_hortonrunoff::Vector{String} |"-"| 0 | "none" | "none"
-    select_rootzone::Vector{String} |"-"| 0 | "none" | "none"
-    select_fast::Vector{String} |"-"| 0 | "none" | "none"
-    select_slow::Vector{String} |"-"| 0 | "none" | "none"
+    kclass::Vector{Int64} | "-" | 0 | "none" | "none"
+    classes::Vector{String} | "-" | 0 | "none" | "none"
+    select_snow::Vector{String} | "-" | 0 | "none" | "none"
+    select_interception::Vector{String} | "-" | 0 | "none" | "none"
+    select_hortonponding::Vector{String} | "-" | 0 | "none" | "none"
+    select_hortonrunoff::Vector{String} | "-" | 0 | "none" | "none"
+    select_rootzone::Vector{String} | "-" | 0 | "none" | "none"
+    select_fast::Vector{String} | "-" | 0 | "none" | "none"
+    select_slow::Vector{String} | "-" | 0 | "none" | "none"
 
     #fraction of each class
     hrufrac::Vector{SVector{N,T}} | "-"
@@ -223,16 +223,6 @@
 
 end
 
-statevars(::FLEXTOPO) = (
-    :snow,
-    :snowwater,
-    :interceptionstorage,
-    :hortonpondingstorage,
-    :hortonrunoffstorage,
-    :rootzonestorage,
-    :faststorage,
-    :slowstorage,
-)
 
 function common_snow_hbv(flextopo::FLEXTOPO)
     for i = 1:flextopo.n

diff --git a/src/flextopo_model.jl b/src/flextopo_model.jl
@@ -722,8 +722,7 @@ function set_states(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:FlextopoModel
     # read and set states in model object if reinit=true
     if reinit == false
         instate_path = input_path(config, config.state.path_input)
-        state_ncnames = ncnames(config.state)
-        set_states(instate_path, model, state_ncnames; type = Float, dimname = :classes)
+        set_states(instate_path, model; type = Float, dimname = :classes)
 
         # update kinematic wave volume for river and land domain
         lateral.land.volume .= lateral.land.h .* lateral.land.width .* lateral.land.dl
@@ -738,4 +737,4 @@ function set_states(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:FlextopoModel
         end
     end
     return model
-end
+end
diff --git a/src/flow.jl b/src/flow.jl
@@ -1,7 +1,8 @@
 
 abstract type SurfaceFlow end
 
-@get_units @exchange @grid_type @grid_location @with_kw struct SurfaceFlowRiver{T,R,L} <: SurfaceFlow
+@get_units @exchange @grid_type @grid_location @with_kw struct SurfaceFlowRiver{T,R,L} <:
+                                                               SurfaceFlow
     β::T | "-" | 0 | "scalar"                    # constant in Manning's equation
     sl::Vector{T} | "m m-1"                      # Slope [m m⁻¹]
     n::Vector{T} | "s m-1/3"                     # Manning's roughness [s m⁻⅓]
@@ -36,8 +37,9 @@ abstract type SurfaceFlow end
     # end
 end
 
-@get_units @exchange @grid_type @grid_location @with_kw struct SurfaceFlowLand{T} <: SurfaceFlow
-    β::T | "-"  | 0 | "scalar"                      # constant in Manning's equation
+@get_units @exchange @grid_type @grid_location @with_kw struct SurfaceFlowLand{T} <:
+                                                               SurfaceFlow
+    β::T | "-" | 0 | "scalar"                       # constant in Manning's equation
     sl::Vector{T} | "m m-1"                         # Slope [m m⁻¹]
     n::Vector{T} | "s m-1/3"                        # Manning's roughness [s m⁻⅓]
     dl::Vector{T} | "m"                             # Drain length [m]
@@ -57,7 +59,7 @@ end
     α::Vector{T} | "s3/5 m1/5"                      # Constant in momentum equation A = αQᵝ, based on Manning's equation
     cel::Vector{T} | "m s-1"                        # Celerity of the kinematic wave
     to_river::Vector{T} | "m3 s-1"                  # Part of overland flow [m³ s⁻¹] that flows to the river
-    kinwave_it::Bool | "-" | 0 |  "none" | "none"   # Boolean for iterations kinematic wave
+    kinwave_it::Bool | "-" | 0 | "none" | "none"    # Boolean for iterations kinematic wave
 end
 
 function initialize_surfaceflow_land(nc, config, inds; sl, dl, width, iterate, tstep, Δt)
@@ -177,8 +179,6 @@ function initialize_surfaceflow_river(
     return sf_river
 end
 
-statevars(::SurfaceFlowRiver) = (:q, :h, :h_av)
-statevars(::SurfaceFlowLand) = (:q, :h, :h_av)
 
 function update(sf::SurfaceFlowLand, network, frac_toriver)
     @unpack graph, subdomain_order, topo_subdomain, indices_subdomain, upstream_nodes =
@@ -388,7 +388,7 @@ function stable_timestep(sf::S) where {S<:SurfaceFlow}
     return Δt, its
 end
 
-@get_units @exchange @grid_type  @grid_location @with_kw struct LateralSSF{T}
+@get_units @exchange @grid_type @grid_location @with_kw struct LateralSSF{T}
     kh₀::Vector{T} | "m d-1"               # Horizontal hydraulic conductivity at soil surface [m d⁻¹]
     f::Vector{T} | "m-1"                   # A scaling parameter [m⁻¹] (controls exponential decline of kh₀)
     soilthickness::Vector{T} | "m"         # Soil thickness [m]
@@ -412,7 +412,6 @@ end
     end
 end
 
-statevars(::LateralSSF) = (:ssf,)
 
 function update(ssf::LateralSSF, network, frac_toriver)
     @unpack subdomain_order, topo_subdomain, indices_subdomain, upstream_nodes = network
@@ -475,7 +474,7 @@ end
     h_thresh::T | "m" | 0 | "scalar"                        # depth threshold for calculating flow
     Δt::T | "s" | 0 | "none" | "none"                       # model time step [s]
     q::Vector{T} | "m3 s-1" | _ | "edge"                    # river discharge (subgrid channel)
-    q0::Vector{T} | "m3 s-1"| _ | "edge"                    # river discharge (subgrid channel) at previous time step
+    q0::Vector{T} | "m3 s-1" | _ | "edge"                   # river discharge (subgrid channel) at previous time step
     q_av::Vector{T} | "m3 s-1" | _ | "edge"                 # average river channel (+ floodplain) discharge [m³ s⁻¹]
     q_channel_av::Vector{T} | "m3 s-1"                      # average river channel discharge [m³ s⁻¹]
     zb_max::Vector{T} | "m"                                 # maximum channel bed elevation
@@ -489,7 +488,7 @@ end
     h_av::Vector{T} | "m"                                   # average water depth
     dl::Vector{T} | "m"                                     # river length
     dl_at_link::Vector{T} | "m" | _ | "edge"                # river length at edge/link
-    width::Vector{T} | "m"                      	        # river width
+    width::Vector{T} | "m"                                  # river width
     width_at_link::Vector{T} | "m" | _ | "edge"             # river width at edge/link
     a::Vector{T} | "m2" | _ | "edge"                        # flow area at edge/link
     r::Vector{T} | "m" | _ | "edge"                         # wetted perimeter at edge/link
@@ -972,7 +971,7 @@ const dirs = (:yd, :xd, :xu, :yu)
     runoff::Vector{T} | "m3 s-1"                            # runoff from hydrological model
     h::Vector{T} | "m"                                      # water depth of cell (for river cells the reference is the river bed elevation `zb`)
     z::Vector{T} | "m"                                      # elevation of cell
-    froude_limit::Bool | "-" | 0  |"none" | "none"          # if true a check is performed if froude number > 1.0 (algorithm is modified)
+    froude_limit::Bool | "-" | 0 | "none" | "none"          # if true a check is performed if froude number > 1.0 (algorithm is modified)
     rivercells::Vector{Bool} | "-"                          # river cells
     h_av::Vector{T} | "m"                                   # average water depth (for river cells the reference is the river bed elevation `zb`)
 end
@@ -1319,7 +1318,7 @@ function update(sw::ShallowWaterLand{T}, swr::ShallowWaterRiver{T}, network, Δt
     end
 end
 
-""" 
+"""
     FloodPlainProfile
 
 Floodplain `volume` is a function of `depth` (flood depth intervals). Based on the

diff --git a/src/groundwater/boundary_conditions.jl b/src/groundwater/boundary_conditions.jl
@@ -41,7 +41,8 @@ function flux!(Q, river::River, aquifer)
 end
 
 
-@get_units @exchange @grid_type @grid_location struct Drainage{T} <: AquiferBoundaryCondition
+@get_units @exchange @grid_type @grid_location struct Drainage{T} <:
+                                                      AquiferBoundaryCondition
     elevation::Vector{T} | "m"
     conductance::Vector{T} | "m2 d-1"
     flux::Vector{T} | "m3 d-1"
@@ -79,7 +80,8 @@ function flux!(Q, headboundary::HeadBoundary, aquifer)
 end
 
 
-@get_units @exchange @grid_type @grid_location struct Recharge{T} <: AquiferBoundaryCondition
+@get_units @exchange @grid_type @grid_location struct Recharge{T} <:
+                                                      AquiferBoundaryCondition
     rate::Vector{T} | "m d-1"
     flux::Vector{T} | "m3 d-1"
     index::Vector{Int} | "-"

diff --git a/src/hbv.jl b/src/hbv.jl
@@ -57,14 +57,6 @@
     runoff::Vector{T}                   # Total specific runoff per cell [mm Δt⁻¹]
 end
 
-statevars(::HBV) = (
-    :soilmoisture,
-    :snow,
-    :snowwater,
-    :upperzonestorage,
-    :lowerzonestorage,
-    :interceptionstorage,
-)
 
 function update_until_snow(hbv::HBV, config)
 

diff --git a/src/hbv_model.jl b/src/hbv_model.jl
@@ -421,8 +421,7 @@ function set_states(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:HbvModel}
     if reinit == false
         instate_path = input_path(config, config.state.path_input)
         @info "Set initial conditions from state file `$instate_path`."
-        state_ncnames = ncnames(config.state)
-        set_states(instate_path, model, state_ncnames; type = Float)
+        set_states(instate_path, model; type = Float)
         # update kinematic wave volume for river and land domain
         @unpack lateral = model
         # makes sure land cells with zero flow width are set to zero q and h
@@ -447,4 +446,4 @@ function set_states(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:HbvModel}
     end
 
     return model
-end
+end
diff --git a/src/io.jl b/src/io.jl
@@ -955,7 +955,7 @@ function prepare_writer(
     # create a separate state output NetCDF that will hold the last timestep of all states
     # but only if config.state.path_output has been set
     if haskey(config, "state") && haskey(config.state, "path_output")
-        state_ncnames = ncnames(config.state)
+        state_ncnames = check_states(config)
         state_map = out_map(state_ncnames, modelmap)
         nc_state_path = output_path(config, config.state.path_output)
         @info "Create a state output NetCDF file `$nc_state_path`."

diff --git a/src/reservoir_lake.jl b/src/reservoir_lake.jl
@@ -22,7 +22,6 @@
     end
 end
 
-statevars(::SimpleReservoir) = (:volume,)
 
 function initialize_simple_reservoir(config, nc, inds_riv, nriv, pits, Δt)
     # read only reservoir data if reservoirs true
@@ -461,18 +460,6 @@ function maximum_storage(storfunc, outflowfunc, area, sh, hq)
     return maxstorage
 end
 
-statevars(::Lake) = (:waterlevel,)
-
-"""
-    statevars(::SBM; snow=false)
-    statevars(::LateralSSF)
-    statevars(::SurfaceFlow)
-    statevars(::NaturalLake)
-    statevars(::SimpleReservoir)
-
-Get a tuple of symbols representing the fields that are model states.
-"""
-function statevars end
 
 function interpolate_linear(x, xp, fp)
     if x <= minimum(xp)