Add Basin / subgrid_time table

core/src/Ribasim.jl

@@ -65,6 +65,7 @@ using PreallocationTools: LazyBufferCache
 # basin profiles and TabulatedRatingCurve. See also the node
 # references in the docs.
 using DataInterpolations:
+    ConstantInterpolation,
     LinearInterpolation,
     LinearInterpolationIntInv,
     invert_integral,

core/src/callback.jl

@@ -671,12 +671,28 @@ function apply_parameter_update!(parameter_update)::Nothing
 end
 
 function update_subgrid_level!(integrator)::Nothing
-    (; p) = integrator
+    (; p, t) = integrator
     du = get_du(integrator)
     basin_level = p.basin.current_properties.current_level[parent(du)]
     subgrid = integrator.p.subgrid
-    for (i, (index, interp)) in enumerate(zip(subgrid.basin_index, subgrid.interpolations))
-        subgrid.level[i] = interp(basin_level[index])
+
+    # First update the all the subgrids with static h(h) relations
+    for (level_index, basin_index, hh_itp) in zip(
+        subgrid.level_index_static,
+        subgrid.basin_index_static,
+        subgrid.interpolations_static,
+    )
+        subgrid.level[level_index] = hh_itp(basin_level[basin_index])
+    end
+    # Then update the subgrids with dynamic h(h) relations
+    for (level_index, basin_index, lookup) in zip(
+        subgrid.level_index_time,
+        subgrid.basin_index_time,
+        subgrid.current_interpolation_index,
+    )
+        itp_index = lookup(t)
+        hh_itp = subgrid.interpolations_time[itp_index]
+        subgrid.level[level_index] = hh_itp(basin_level[basin_index])
     end
 end
 

core/src/parameter.jl

@@ -105,6 +105,7 @@ end
 
 Base.to_index(id::NodeID) = Int(id.value)
 
+"LinearInterpolation from a Float64 to a Float64"
 const ScalarInterpolation = LinearInterpolation{
     Vector{Float64},
     Vector{Float64},
@@ -114,6 +115,10 @@ const ScalarInterpolation = LinearInterpolation{
     (1,),
 }
 
+"ConstantInterpolation from a Float64 to an Int, used to look up indices"
+const IndexLookup =
+    ConstantInterpolation{Vector{Int64}, Vector{Float64}, Vector{Float64}, Int64, (1,)}
+
 set_zero!(v) = v .= zero(eltype(v))
 const Cache = LazyBufferCache{Returns{Int}, typeof(set_zero!)}
 
@@ -876,10 +881,26 @@ end
 
 "Subgrid linearly interpolates basin levels."
 @kwdef struct Subgrid
-    subgrid_id::Vector{Int32}
-    basin_index::Vector{Int32}
-    interpolations::Vector{ScalarInterpolation}
+    # level of each subgrid (static and dynamic) ordered by subgrid_id
     level::Vector{Float64}
+    # static
+    subgrid_id_static::Vector{Int32}
+    # index into the basin.current_level vector for each static subgrid_id
+    basin_index_static::Vector{Int}
+    # index into the subgrid.level vector for each static subgrid_id
+    level_index_static::Vector{Int}
+    # per subgrid one relation
+    interpolations_static::Vector{ScalarInterpolation}
+    # dynamic
+    subgrid_id_time::Vector{Int32}
+    # index into the basin.current_level vector for each dynamic subgrid_id
+    basin_index_time::Vector{Int}
+    # index into the subgrid.level vector for each dynamic subgrid_id
+    level_index_time::Vector{Int}
+    # per subgrid n relations, n being the number of timesteps for that subgrid
+    interpolations_time::Vector{ScalarInterpolation}
+    # per subgrid 1 lookup from t to an index in interpolations_time
+    current_interpolation_index::Vector{IndexLookup}
 end
 
 """

core/src/read.jl

@@ -187,11 +187,18 @@
     return out, !errors
 end
 
+"""
+Validate the split of node IDs between static and time tables.
+
+For node types that can have a part of the parameters defined statically and a part dynamically,
+this checks if each ID is defined exactly once in either table.
+"""
 function static_and_time_node_ids(
     db::DB,
     static::StructVector,
     time::StructVector,
-    node_type::String,
+    node_type::String;
+    is_complete::Bool = true,
 )::Tuple{Set{NodeID}, Set{NodeID}, Vector{NodeID}, Bool}
     ids = get_ids(db, node_type)
     idx = searchsortedfirst.(Ref(ids), static.node_id)
@@ -205,7 +212,7 @@
         errors = true
         @error "$node_type cannot be in both static and time tables, found these node IDs in both: $doubles."
     end
-    if !issetequal(node_ids, union(static_node_ids, time_node_ids))
+    if is_complete && !issetequal(node_ids, union(static_node_ids, time_node_ids))
         errors = true
         @error "$node_type node IDs don't match."
     end
@@ -1226,17 +1233,42 @@
     )
 end
 
+function push_lookup!(
+    current_interpolation_index::Vector{IndexLookup},
+    lookup_index::Vector{Int},
+    lookup_time::Vector{Float64},
+)
+    index_lookup = ConstantInterpolation(
+        lookup_index,
+        lookup_time;
+        extrapolate = true,
+        cache_parameters = true,
+    )
+    push!(current_interpolation_index, index_lookup)
+end
+
 function Subgrid(db::DB, config::Config, basin::Basin)::Subgrid
-    node_to_basin = Dict(node_id => index for (index, node_id) in enumerate(basin.node_id))
-    tables = load_structvector(db, config, BasinSubgridV1)
+    time = load_structvector(db, config, BasinSubgridTimeV1)
+    static = load_structvector(db, config, BasinSubgridV1)
 
-    subgrid_ids = Int32[]
-    basin_index = Int32[]
-    interpolations = ScalarInterpolation[]
+    # Since not all Basins need to have subgrids, don't enforce completeness.
+    _, _, _, valid =
+        static_and_time_node_ids(db, static, time, "Basin"; is_complete = false)
+    if !valid
+        error("Problems encountered when parsing Subgrid static and time node IDs.")
+    end
+
+    node_to_basin = Dict{Int32, Int}(
+        Int32(node_id) => index for (index, node_id) in enumerate(basin.node_id)
+    )
+    subgrid_id_static = Int32[]
+    basin_index_static = Int[]
+    interpolations_static = ScalarInterpolation[]
     has_error = false
-    for group in IterTools.groupby(row -> row.subgrid_id, tables)
+
+    for group in IterTools.groupby(row -> row.subgrid_id, static)
         subgrid_id = first(getproperty.(group, :subgrid_id))
-        node_id = NodeID(NodeType.Basin, first(getproperty.(group, :node_id)), db)
+        node_id = first(getproperty.(group, :node_id))
         basin_level = getproperty.(group, :basin_level)
         subgrid_level = getproperty.(group, :subgrid_level)
 
@@ -1247,24 +1279,108 @@
             # Ensure it doesn't extrapolate before the first value.
             pushfirst!(subgrid_level, first(subgrid_level))
             pushfirst!(basin_level, nextfloat(-Inf))
-            new_interp = LinearInterpolation(
+            hh_itp = LinearInterpolation(
                 subgrid_level,
                 basin_level;
                 extrapolate = true,
                 cache_parameters = true,
             )
-            push!(subgrid_ids, subgrid_id)
-            push!(basin_index, node_to_basin[node_id])
-            push!(interpolations, new_interp)
+            push!(subgrid_id_static, subgrid_id)
+            push!(basin_index_static, node_to_basin[node_id])
+            push!(interpolations_static, hh_itp)
         else
             has_error = true
         end
     end
 
     has_error && error("Invalid Basin / subgrid table.")
-    level = fill(NaN, length(subgrid_ids))
+    subgrid_id_time = Int32[]
+    basin_index_time = Int[]
+    interpolations_time = ScalarInterpolation[]
+    current_interpolation_index = IndexLookup[]
+
+    # Push the first subgrid_id and basin_index
+    if length(time) > 0
+        push!(subgrid_id_time, first(time.subgrid_id))
+        push!(basin_index_time, node_to_basin[first(time.node_id)])
+    end
+
+    # Initialize index_lookup contents
+    lookup_time = Float64[]
+    lookup_index = Int[]
+
+    interpolation_index = 0
+    for group in IterTools.groupby(row -> (row.subgrid_id, row.time), time)
+        interpolation_index += 1
+        subgrid_id = first(getproperty.(group, :subgrid_id))
+        time_group = seconds_since(first(getproperty.(group, :time)), config.starttime)
+        node_id = first(getproperty.(group, :node_id))
+        basin_level = getproperty.(group, :basin_level)
+        subgrid_level = getproperty.(group, :subgrid_level)
+
+        is_valid =
+            valid_subgrid(subgrid_id, node_id, node_to_basin, basin_level, subgrid_level)
 
-    return Subgrid(; subgrid_id = subgrid_ids, basin_index, interpolations, level)
+        if is_valid
+            # Ensure it doesn't extrapolate before the first value.
+            pushfirst!(subgrid_level, first(subgrid_level))
+            pushfirst!(basin_level, nextfloat(-Inf))
+            hh_itp = LinearInterpolation(
+                subgrid_level,
+                basin_level;
+                extrapolate = true,
+                cache_parameters = true,
+            )
+            # # These should only be pushed when the subgrid_id has changed
+            if subgrid_id_time[end] != subgrid_id
+                # Push the completed index_lookup of the previous subgrid_id
+                push_lookup!(current_interpolation_index, lookup_index, lookup_time)
+                # Push the new subgrid_id and basin_index
+                push!(subgrid_id_time, subgrid_id)
+                push!(basin_index_time, node_to_basin[node_id])
+                # Start new index_lookup contents
+                lookup_time = Float64[]
+                lookup_index = Int[]
+            end
+            push!(lookup_index, interpolation_index)
+            push!(lookup_time, time_group)
+            push!(interpolations_time, hh_itp)
+        else
+            has_error = true
+        end
+    end
+
+    # Push completed IndexLookup of the last group
+    if interpolation_index > 0
+        push_lookup!(current_interpolation_index, lookup_index, lookup_time)
+    end
+
+    has_error && error("Invalid Basin / subgrid_time table.")
+    level = fill(NaN, length(subgrid_id_static) + length(subgrid_id_time))
+
+    # Find the level indices
+    level_index_static = zeros(Int, length(subgrid_id_static))
+    level_index_time = zeros(Int, length(subgrid_id_time))
+    subgrid_ids = sort(vcat(subgrid_id_static, subgrid_id_time))
+    for (i, subgrid_id) in enumerate(subgrid_id_static)
+        level_index_static[i] = findsorted(subgrid_ids, subgrid_id)
+    end
+    for (i, subgrid_id) in enumerate(subgrid_id_time)
+        level_index_time[i] = findsorted(subgrid_ids, subgrid_id)
+    end
+
+    return Subgrid(;
+        level,
+        subgrid_id_static,
+        basin_index_static,
+        level_index_static,
+        interpolations_static,
+        subgrid_id_time,
+        basin_index_time,
+        level_index_time,
+        interpolations_time,
+        current_interpolation_index,
+    )
 end
 
 function Allocation(db::DB, config::Config, graph::MetaGraph)::Allocation

core/src/schema.jl

@@ -10,6 +10,7 @@
 @schema "ribasim.basin.profile" BasinProfile
 @schema "ribasim.basin.state" BasinState
 @schema "ribasim.basin.subgrid" BasinSubgrid
+@schema "ribasim.basin.subgridtime" BasinSubgridTime
 @schema "ribasim.basin.concentration" BasinConcentration
 @schema "ribasim.basin.concentrationexternal" BasinConcentrationExternal
 @schema "ribasim.basin.concentrationstate" BasinConcentrationState
@@ -58,9 +59,13 @@ function nodetype(
     type_string = string(T)
     elements = split(type_string, '.'; limit = 3)
     last_element = last(elements)
+    # Special case last elements that need an underscore
     if startswith(last_element, "concentration") && length(last_element) > 13
         elements[end] = "concentration_$(last_element[14:end])"
     end
+    if last_element == "subgridtime"
+        elements[end] = "subgrid_time"
+    end
     if isnode(sv)
         n = elements[2]
         k = Symbol(elements[3])
@@ -150,6 +155,14 @@ end
     subgrid_level::Float64
 end
 
+@version BasinSubgridTimeV1 begin
+    subgrid_id::Int32
+    node_id::Int32
+    time::DateTime
+    basin_level::Float64
+    subgrid_level::Float64
+end
+
 @version LevelBoundaryStaticV1 begin
     node_id::Int32
     active::Union{Missing, Bool}

core/src/validation.jl

@@ -321,8 +321,8 @@ Validate the entries for a single subgrid element.
 """
 function valid_subgrid(
     subgrid_id::Int32,
-    node_id::NodeID,
-    node_to_basin::Dict{NodeID, Int},
+    node_id::Int32,
+    node_to_basin::Dict{Int32, Int},
     basin_level::Vector{Float64},
     subgrid_level::Vector{Float64},
 )::Bool

core/src/write.jl

@@ -376,11 +376,14 @@
     (; t, saveval) = saved.subgrid_level
     subgrid = integrator.p.subgrid
 
-    nelem = length(subgrid.subgrid_id)
+    nelem = length(subgrid.level)
     ntsteps = length(t)
 
     time = repeat(datetime_since.(t, config.starttime); inner = nelem)
-    subgrid_id = repeat(subgrid.subgrid_id; outer = ntsteps)
+    subgrid_id = repeat(
+        sort(vcat(subgrid.subgrid_id_static, subgrid.subgrid_id_time));
+        outer = ntsteps,
+    )
     subgrid_level = FlatVector(saveval)
     return (; time, subgrid_id, subgrid_level)
 end
@@ -412,9 +415,9 @@
     mkpath(dirname(path))
     try
         Arrow.write(path, table; compress, metadata)
-    catch
+    catch e
         @error "Failed to write results, file may be locked." path
-        error("Failed to write results.")
+        rethrow(e)
     end
     return nothing
 end

core/test/run_models_test.jl

@@ -583,3 +583,29 @@ end
     @test all(isapprox.(Δinf[1:2:end], 25.0; atol = 1e-10))
     @test all(Δinf[2:2:end] .== 0.0)
 end
+
+@testitem "two_basin" begin
+    using DataFrames: DataFrame, nrow
+    using Dates: DateTime
+    import BasicModelInterface as BMI
+
+    toml_path = normpath(@__DIR__, "../../generated_testmodels/two_basin/ribasim.toml")
+    model = Ribasim.run(toml_path)
+    df = DataFrame(Ribasim.subgrid_level_table(model))
+
+    ntime = 367
+    @test nrow(df) == ntime * 2
+    @test df.subgrid_id == repeat(1:2; outer = ntime)
+    @test extrema(df.time) == (DateTime(2020), DateTime(2021))
+    @test all(df.subgrid_level[1:2] .== 0.01)
+
+    # After a month the h(h) of subgrid_id 2 increases by a meter
+    i_change = searchsortedfirst(df.time, DateTime(2020, 2))
+    @test df.subgrid_level[i_change + 1] - df.subgrid_level[i_change - 1] ≈ 1.0f0
+
+    # Besides the 1 meter shift the h(h) relations are 1:1
+    basin_level = copy(BMI.get_value_ptr(model, "basin.level"))
+    basin_level[2] += 1
+    @test basin_level ≈ df.subgrid_level[(end - 1):end]
+    @test basin_level ≈ model.integrator.p.subgrid.level
+end