diff --git a/core/src/Ribasim.jl b/core/src/Ribasim.jl index 03edb77fc..99c0124fd 100644 --- a/core/src/Ribasim.jl +++ b/core/src/Ribasim.jl @@ -63,8 +63,9 @@ using DataInterpolations: invert_integral, derivative, integral, - AbstractInterpolation -using DataInterpolations.ExtrapolationType: Constant, Extension, Linear + AbstractInterpolation, + ExtrapolationType +using DataInterpolations.ExtrapolationType: Constant, Periodic, Extension, Linear # Modeling language for Mathematical Optimization. # Used for allocation, see the docs: https://ribasim.org/dev/allocation.html diff --git a/core/src/model.jl b/core/src/model.jl index 58e3a354c..18ab13b69 100644 --- a/core/src/model.jl +++ b/core/src/model.jl @@ -55,6 +55,12 @@ function Model(config::Config)::Model error("Invalid link types found.") end + # for Float32 this method allows max ~1000 year simulations without accuracy issues + t_end = seconds_since(config.endtime, config.starttime) + @assert eps(t_end) < 3600 "Simulation time too long" + t0 = zero(t_end) + timespan = (t0, t_end) + local parameters, tstops try parameters = Parameters(db, config) @@ -63,7 +69,16 @@ function Model(config::Config)::Model error("Invalid discrete control state definition(s).") end - (; pid_control, graph, outlet, pump, basin, tabulated_rating_curve) = parameters + (; + pid_control, + graph, + outlet, + pump, + basin, + tabulated_rating_curve, + level_boundary, + flow_boundary, + ) = parameters if !valid_pid_connectivity(pid_control.node_id, pid_control.listen_node_id, graph) error("Invalid PidControl connectivity.") end @@ -83,11 +98,8 @@ function Model(config::Config)::Model # tell the solver to stop when new data comes in tstops = Vector{Float64}[] for schema_version in [ - BasinTimeV1, DiscreteControlConditionV1, - FlowBoundaryTimeV1, FlowDemandTimeV1, - LevelBoundaryTimeV1, LevelDemandTimeV1, PidControlTimeV1, TabulatedRatingCurveTimeV1, @@ -97,6 +109,21 @@ function Model(config::Config)::Model push!(tstops, get_tstops(time_schema.time, config.starttime)) end + # Tell the solver to stop at all data points from periodically extrapolated + # if applicable, otherwise just the original timeseries + for interpolations in [ + basin.forcing.precipitation, + basin.forcing.potential_evaporation, + basin.forcing.drainage, + basin.forcing.infiltration, + level_boundary.level, + flow_boundary.flow_rate, + ] + for itp in interpolations + push!(tstops, get_cyclic_tstops(itp, t_end)) + end + end + finally # always close the database, also in case of an error close(db) @@ -113,12 +140,6 @@ function Model(config::Config)::Model # The Solver algorithm alg = algorithm(config.solver; u0) - # for Float32 this method allows max ~1000 year simulations without accuracy issues - t_end = seconds_since(config.endtime, config.starttime) - @assert eps(t_end) < 3600 "Simulation time too long" - t0 = zero(t_end) - timespan = (t0, t_end) - # Synchronize level with storage set_current_basin_properties!(du0, u0, parameters, t0) diff --git a/core/src/read.jl b/core/src/read.jl index c4dc02b1b..29f83c8a3 100644 --- a/core/src/read.jl +++ b/core/src/read.jl @@ -1,16 +1,23 @@ """ Process the data in the static and time tables for a given node type. -The 'defaults' named tuple dictates how missing data is filled in. -'time_interpolatables' is a vector of Symbols of parameter names -for which a time interpolation (linear) object must be constructed. -The control mapping for DiscreteControl is also constructed in this function. -This function currently does not support node states that are defined by more -than one row in a table, as is the case for TabulatedRatingCurve. + +Inputs: +- `db`: The Ribasim GeoPackage database +- `config`: The configuration object +- `node_type_string`: A string representation of a node type, e.g. "LinearResistance" + +Optional inputs: +- `static`: The StructVector with the static node data +- `time`: The StructVector with the transient node data +- `time_interpolatables`: The names of the variables which can be interpolated over time +- `interpolation_type`: The type of interpolation used for the time interpolatable variables +- `is_complete`: Toggles whether it should be checked that an id that exists in the Node table should be + in either the static or time table. Defaults to `true`. """ function parse_static_and_time( db::DB, config::Config, - node_type::Type; + node_type_string::String; static::Union{StructVector, Nothing} = nothing, time::Union{StructVector, Nothing} = nothing, defaults::NamedTuple = (; active = true), @@ -30,12 +37,12 @@ function parse_static_and_time( # The types of the variables that can define a control state parameter_types = collect(fieldtypes(static_type))[mask] - # A vector of vectors, for each parameter the (initial) values for all nodes + # A vector of vectors and dicts, for each parameter the (initial) values for all nodes # of the current type vals_out = [] - node_type_string = String(split(string(node_type), '.')[end]) node_ids = get_node_ids(db, node_type_string) + cyclic_times = get_cyclic_time(db, node_type_string) ids = Int32.(node_ids) n_nodes = length(node_ids) @@ -109,7 +116,7 @@ function parse_static_and_time( errors = false trivial_timespan = [0.0, prevfloat(Inf)] - for node_id in node_ids + for (cyclic_time, node_id) in zip(cyclic_times, node_ids) if node_id in static_node_ids # The interval of rows of the static table that have the current node_id rows = searchsorted(static_node_id_vec, node_id) @@ -167,6 +174,7 @@ function parse_static_and_time( default_value = hasproperty(defaults, parameter_name) ? defaults[parameter_name] : NaN, interpolation_type, + extrapolation = cyclic_time ? Periodic : Constant, ) else # Activity of transient nodes is assumed to be true @@ -188,6 +196,7 @@ function parse_static_and_time( [val, val], trivial_timespan; cache_parameters = true, + extrapolation = Constant, ) end getfield(out, parameter_name)[node_id.idx] = val @@ -407,7 +416,7 @@ function LinearResistance(db::DB, config::Config, graph::MetaGraph)::LinearResis static = load_structvector(db, config, LinearResistanceStaticV1) defaults = (; max_flow_rate = Inf, active = true) parsed_parameters, valid = - parse_static_and_time(db, config, LinearResistance; static, defaults) + parse_static_and_time(db, config, "LinearResistance"; static, defaults) if !valid error( @@ -565,7 +574,8 @@ function ManningResistance( basin::Basin, )::ManningResistance static = load_structvector(db, config, ManningResistanceStaticV1) - parsed_parameters, valid = parse_static_and_time(db, config, ManningResistance; static) + parsed_parameters, valid = + parse_static_and_time(db, config, "ManningResistance"; static) if !valid error("Errors occurred when parsing ManningResistance data.") @@ -603,8 +613,14 @@ function LevelBoundary(db::DB, config::Config)::LevelBoundary end time_interpolatables = [:level] - parsed_parameters, valid = - parse_static_and_time(db, config, LevelBoundary; static, time, time_interpolatables) + parsed_parameters, valid = parse_static_and_time( + db, + config, + "LevelBoundary"; + static, + time, + time_interpolatables, + ) substances = get_substances(db, config) concentration = zeros(length(node_ids), length(substances)) @@ -638,8 +654,14 @@ function FlowBoundary(db::DB, config::Config, graph::MetaGraph)::FlowBoundary end time_interpolatables = [:flow_rate] - parsed_parameters, valid = - parse_static_and_time(db, config, FlowBoundary; static, time, time_interpolatables) + parsed_parameters, valid = parse_static_and_time( + db, + config, + "FlowBoundary"; + static, + time, + time_interpolatables, + ) for itp in parsed_parameters.flow_rate if any(itp.u .< 0.0) @@ -679,7 +701,7 @@ function Pump(db::DB, config::Config, graph::MetaGraph)::Pump max_downstream_level = Inf, active = true, ) - parsed_parameters, valid = parse_static_and_time(db, config, Pump; static, defaults) + parsed_parameters, valid = parse_static_and_time(db, config, "Pump"; static, defaults) if !valid error("Errors occurred when parsing Pump data.") @@ -714,7 +736,7 @@ function Outlet(db::DB, config::Config, graph::MetaGraph)::Outlet max_downstream_level = Inf, active = true, ) - parsed_parameters, valid = parse_static_and_time(db, config, Outlet; static, defaults) + parsed_parameters, valid = parse_static_and_time(db, config, "Outlet"; static, defaults) if !valid error("Errors occurred when parsing Outlet data.") @@ -852,7 +874,7 @@ function Basin(db::DB, config::Config, graph::MetaGraph)::Basin parsed_parameters, valid = parse_static_and_time( db, config, - Basin; + "Basin"; static, time, time_interpolatables, @@ -1182,7 +1204,7 @@ function PidControl(db::DB, config::Config, graph::MetaGraph)::PidControl time_interpolatables = [:target, :proportional, :integral, :derivative] parsed_parameters, valid = - parse_static_and_time(db, config, PidControl; static, time, time_interpolatables) + parse_static_and_time(db, config, "PidControl"; static, time, time_interpolatables) if !valid error("Errors occurred when parsing PidControl data.") @@ -1401,7 +1423,7 @@ function LevelDemand(db::DB, config::Config)::LevelDemand parsed_parameters, valid = parse_static_and_time( db, config, - LevelDemand; + "LevelDemand"; static, time, time_interpolatables = [:min_level, :max_level], @@ -1429,7 +1451,7 @@ function FlowDemand(db::DB, config::Config)::FlowDemand parsed_parameters, valid = parse_static_and_time( db, config, - FlowDemand; + "FlowDemand"; static, time, time_interpolatables = [:demand], @@ -1796,6 +1818,11 @@ function get_node_ids(db::DB, node_type)::Vector{NodeID} return node_ids end +function get_cyclic_time(db::DB, node_type)::Vector{Bool} + sql = "SELECT cyclic_time FROM Node WHERE node_type = $(esc_id(node_type)) ORDER BY node_id" + return only(execute(columntable, db, sql)) +end + function exists(db::DB, tablename::String) query = execute( db, diff --git a/core/src/util.jl b/core/src/util.jl index de2de7dd0..70ceaf457 100644 --- a/core/src/util.jl +++ b/core/src/util.jl @@ -65,7 +65,6 @@ function get_level_from_storage(basin::Basin, state_idx::Int, storage) end end -"Linear interpolation of a scalar with constant extrapolation." function get_scalar_interpolation( starttime::DateTime, time::AbstractVector, @@ -73,22 +72,34 @@ function get_scalar_interpolation( param::Symbol; default_value::Float64 = 0.0, interpolation_type::Type{<:AbstractInterpolation}, + extrapolation::ExtrapolationType.T = Constant, )::interpolation_type rows = searchsorted(time.node_id, node_id) parameter = getproperty(time, param)[rows] parameter = coalesce.(parameter, default_value) times = seconds_since.(time.time[rows], starttime) + errors = false + if !allunique(times) - @error "The time series for $node_id has repeated times, this can not be interpolated." - error("Invalid time series.") + errors = true + @error "(One of) the time series for $node_id has repeated times, this can not be interpolated." end - return interpolation_type( - parameter, - times; - extrapolation = Constant, - cache_parameters = true, - ) + + if extrapolation == Periodic + if !(all(isnan, parameter) || (first(parameter) == last(parameter))) + errors = true + @error "$node_id is denoted as cyclic but in (one of) its time series the first and last value are not the same." + end + + if length(times) < 2 + errors = true + @error "$node_id is denoted as cyclic but (one of) its time series has fewer than 2 data points." + end + end + + errors && error("Invalid time series.") + return interpolation_type(parameter, times; extrapolation, cache_parameters = true) end """ @@ -1134,3 +1145,32 @@ function find_index(x::Symbol, s::OrderedSet{Symbol}) end error(lazy"$x not found in $s.") end + +function get_cyclic_tstops(itp::AbstractInterpolation, endtime::Float64)::Vector{Float64} + # The length of the period + T = last(itp.t) - first(itp.t) + + # How many periods back from first(itp.t) are needed + nT_back = itp.extrapolation_left == Periodic ? Int(ceil((first(itp.t)) / T)) : 0 + + # How many periods forward from first(itp.t) are needed + nT_forward = + itp.extrapolation_right == Periodic ? Int(ceil((endtime - first(itp.t)) / T)) : 0 + + tstops = Float64[] + + for i in (-nT_back):nT_forward + # Append the timepoints of the interpolation shifted by an integer amount of + # periods to the tstops, filtering out values outside the simulation period + if i == nT_forward + append!(tstops, filter(t -> 0 ≤ t ≤ endtime, itp.t .+ i * T)) + else + # Because of floating point errors last(itp.t) = first(itp.t) + T + # does not always hold exactly, so to prevent that these become separate + # very close tstops we only use the last time point of the period in the last period + append!(tstops, filter(t -> 0 ≤ t ≤ endtime, itp.t[1:(end - 1)] .+ i * T)) + end + end + + return tstops +end diff --git a/core/test/time_test.jl b/core/test/time_test.jl index bcc2eac16..61445025a 100644 --- a/core/test/time_test.jl +++ b/core/test/time_test.jl @@ -72,3 +72,37 @@ end # Given that precipitation stops after 15 of the 20 days @test mean_precipitation ≈ 3 / 4 * starting_precipitation end + +@testitem "get_cyclic_tstops" begin + using Ribasim: get_cyclic_tstops + using DataInterpolations: LinearInterpolation, ConstantInterpolation + using DataInterpolations.ExtrapolationType: Periodic + + itp = LinearInterpolation(zeros(3), [0.5, 1.0, 1.5]) + @test get_cyclic_tstops(itp, 5.0) == itp.t + + itp = LinearInterpolation(zeros(3), [0.5, 1.0, 1.5]; extrapolation = Periodic) + @test get_cyclic_tstops(itp, 5.0) == 0:0.5:5 + + itp = ConstantInterpolation(zeros(2), [0.3, 0.5]; extrapolation = Periodic) + @test get_cyclic_tstops(itp, 1.0) ≈ 0.1:0.2:0.9 +end + +@testitem "cyclic time" begin + using DataInterpolations.ExtrapolationType: Periodic + + toml_path = normpath(@__DIR__, "../../generated_testmodels/cyclic_time/ribasim.toml") + @test ispath(toml_path) + + model = Ribasim.Model(toml_path) + (; level_boundary, flow_boundary, basin) = model.integrator.p + + function test_extrapolation(itp) + @test itp.extrapolation_left == Periodic + @test itp.extrapolation_right == Periodic + end + + test_extrapolation(basin.forcing.precipitation[1]) + test_extrapolation(level_boundary.level[1]) + test_extrapolation(flow_boundary.flow_rate[1]) +end diff --git a/docs/reference/usage.qmd b/docs/reference/usage.qmd index f95753dd7..4059a1bed 100644 --- a/docs/reference/usage.qmd +++ b/docs/reference/usage.qmd @@ -172,10 +172,17 @@ geom | Point | (optional) name | String | (optional, does not have to be unique) subnetwork_id | Int32 | (optional) source_priority | Int32 | (optional, does not have to be unique) +cyclic_time | Bool | (optional, defaults to false) Adding a point geometry to the node table can be helpful to examine models in [QGIS](https://qgis.org/en/site/), as it will show the location of the nodes on the map. The geometry is not used by Ribasim. +## cyclic time series + +When cyclic forcing is set to true for a node in the Node table, every time series associated with that node in the corresponding table(s) will be interpreted as cyclic. That is: the time series is exactly repeated left and right of the original time interval to cover the whole simulation period. For this it is validated that the first and last data values in the timeseries are the same. For instance, quarterly precipitation requires giving values for every quarter at the start of the quarter, and then the value for the first quarter again at the start of the next year. + +Note that periods like months or years are not of constant length in the calendar, so over long simulation periods the timeseries can get out of sync with these periods on the calendar. + # Link {#sec-link} Links define connections between nodes. The only thing that defines an link is the nodes it connects, and in what direction. diff --git a/python/ribasim/ribasim/config.py b/python/ribasim/ribasim/config.py index 67af8991e..5d7bc05c4 100644 --- a/python/ribasim/ribasim/config.py +++ b/python/ribasim/ribasim/config.py @@ -190,6 +190,8 @@ class Node(pydantic.BaseModel): Optionally adds this node to a subnetwork, which is input for the allocation algorithm. source_priority : int Optionally overrides the source priority for this node, which is used in the allocation algorithm. + cyclic_time : bool + Optionally extrapolate forcing timeseries periodically. Defaults to False. """ node_id: NonNegativeInt | None = None @@ -197,6 +199,7 @@ class Node(pydantic.BaseModel): name: str = "" subnetwork_id: int | None = None source_priority: int | None = None + cyclic_time: bool = False model_config = ConfigDict(arbitrary_types_allowed=True, extra="allow") @@ -224,6 +227,7 @@ def into_geodataframe(self, node_type: str, node_id: int) -> GeoDataFrame: "source_priority": pd.Series( [self.source_priority], dtype=pd.Int32Dtype() ), + "cyclic_time": pd.Series([self.cyclic_time], dtype=bool), **extra, }, geometry=[self.geometry], diff --git a/python/ribasim/ribasim/geometry/node.py b/python/ribasim/ribasim/geometry/node.py index 3b2eac19e..753c4cce9 100644 --- a/python/ribasim/ribasim/geometry/node.py +++ b/python/ribasim/ribasim/geometry/node.py @@ -28,6 +28,7 @@ class NodeSchema(_GeoBaseSchema): source_priority: Series[pd.Int32Dtype] = pa.Field( default=pd.NA, nullable=True, coerce=True ) + cyclic_time: Series[bool] = pa.Field(default=False) geometry: GeoSeries[Point] = pa.Field(default=None, nullable=True) @classmethod diff --git a/python/ribasim_testmodels/ribasim_testmodels/__init__.py b/python/ribasim_testmodels/ribasim_testmodels/__init__.py index ed781ef9f..18a933aed 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/__init__.py +++ b/python/ribasim_testmodels/ribasim_testmodels/__init__.py @@ -25,6 +25,7 @@ basic_arrow_model, basic_model, basic_transient_model, + cyclic_time_model, outlet_model, tabulated_rating_curve_model, ) @@ -78,6 +79,7 @@ "concentration_condition_model", "continuous_concentration_condition_model", "connector_node_flow_condition_model", + "cyclic_time_model", "discrete_control_of_pid_control_model", "fair_distribution_model", "flow_boundary_time_model", diff --git a/python/ribasim_testmodels/ribasim_testmodels/basic.py b/python/ribasim_testmodels/ribasim_testmodels/basic.py index 5fd796ba3..0249e5c25 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/basic.py +++ b/python/ribasim_testmodels/ribasim_testmodels/basic.py @@ -5,6 +5,7 @@ import numpy as np import pandas as pd import ribasim +from ribasim import Model from ribasim.config import Experimental, Node from ribasim.input_base import TableModel from ribasim.nodes import ( @@ -20,9 +21,9 @@ from shapely.geometry import MultiPolygon, Point -def basic_model() -> ribasim.Model: +def basic_model() -> Model: # Setup model - model = ribasim.Model( + model = Model( starttime="2020-01-01", endtime="2021-01-01", crs="EPSG:28992", @@ -205,14 +206,14 @@ def basic_model() -> ribasim.Model: return model -def basic_arrow_model() -> ribasim.Model: +def basic_arrow_model() -> Model: model = basic_model() model.basin.profile.set_filepath(Path("profile.arrow")) model.input_dir = Path("input") return model -def basic_transient_model() -> ribasim.Model: +def basic_transient_model() -> Model: """Update the basic model with transient forcing.""" model = basic_model() time = pd.date_range(model.starttime, model.endtime) @@ -262,7 +263,7 @@ def basic_transient_model() -> ribasim.Model: return model -def tabulated_rating_curve_model() -> ribasim.Model: +def tabulated_rating_curve_model() -> Model: """ Set up a model where the upstream Basin has two TabulatedRatingCurve attached. @@ -271,7 +272,7 @@ def tabulated_rating_curve_model() -> ribasim.Model: Only the upstream Basin receives a (constant) precipitation. """ # Setup a model: - model = ribasim.Model( + model = Model( starttime="2020-01-01", endtime="2021-01-01", crs="EPSG:28992", @@ -350,7 +351,7 @@ def tabulated_rating_curve_model() -> ribasim.Model: def outlet_model(): """Set up a basic model with an outlet that encounters various physical constraints.""" - model = ribasim.Model( + model = Model( starttime="2020-01-01", endtime="2021-01-01", crs="EPSG:28992", @@ -392,3 +393,59 @@ def outlet_model(): model.link.add(model.outlet[2], model.basin[3]) return model + + +def cyclic_time_model() -> Model: + model = Model( + starttime="2020-01-01", + endtime="2021-01-01", + crs="EPSG:28992", + ) + + bsn = model.basin.add( + Node(1, Point(0, 0), cyclic_time=True), + [ + basin.Profile(level=[0.0, 1.0], area=100.0), + basin.Time( + time=[ + "2020-01-01", + "2020-04-01", + "2020-07-01", + "2020-10-01", + "2021-01-01", + ], + precipitation=[1.0, 2.0, 1.0, 2.0, 1.0], + ), + basin.State(level=[5.0]), + ], + ) + + lr = model.linear_resistance.add( + Node(2, Point(1, 0)), [linear_resistance.Static(resistance=[1.0])] + ) + + lb = model.level_boundary.add( + Node(3, Point(2, 0), cyclic_time=True), + [ + level_boundary.Time( + time=["2020-01-01", "2020-05-01", "2020-10-01"], + level=[2.0, 3.0, 2.0], + ) + ], + ) + + fb = model.flow_boundary.add( + Node(4, Point(0, 1), cyclic_time=True), + [ + flow_boundary.Time( + time=["2020-01-01", "2020-07-01", "2020-08-01"], + flow_rate=[1.0, 2.0, 1.0], + ) + ], + ) + + model.edge.add(bsn, lr) + model.edge.add(lr, lb) + model.edge.add(fb, bsn) + + return model