From 81efc91bfd6a9bdf80ef4977f58e3ab6e272ae27 Mon Sep 17 00:00:00 2001 From: Martijn Visser Date: Tue, 29 Aug 2023 16:43:21 +0200 Subject: [PATCH 01/13] add reduction_factor to TabulatedRatingCurve If the rating curve prescribes large Q even when the upstream Basin is nearly empty, this leads to an incorrect water balance. --- core/src/jac.jl | 4 ++-- core/src/solve.jl | 17 +++++++++++------ core/test/run_models.jl | 4 ++-- 3 files changed, 15 insertions(+), 10 deletions(-) diff --git a/core/src/jac.jl b/core/src/jac.jl index c474bd2cd..3e3e12143 100644 --- a/core/src/jac.jl +++ b/core/src/jac.jl @@ -404,7 +404,7 @@ function formulate_jac!( controlled_node_idx = findsorted(pump.node_id, controlled_node_id) listened_basin_storage = u.storage[listened_node_idx] - reduction_factor = min(listened_basin_storage, 10.0) / 10.0 + reduction_factor = clamp(listened_basin_storage, 0.0, 10.0) / 10.0 else controlled_node_idx = findsorted(outlet.node_id, controlled_node_id) @@ -414,7 +414,7 @@ function formulate_jac!( id_index(basin.node_id, upstream_node_id) if has_upstream_index upstream_basin_storage = u.storage[upstream_basin_idx] - reduction_factor = min(upstream_basin_storage, 10.0) / 10.0 + reduction_factor = clamp(upstream_basin_storage, 0.0, 10.0) / 10.0 else reduction_factor = 1.0 end diff --git a/core/src/solve.jl b/core/src/solve.jl index e61e68698..c1d040adf 100644 --- a/core/src/solve.jl +++ b/core/src/solve.jl @@ -572,7 +572,7 @@ function continuous_control!( controlled_node_idx = findsorted(pump.node_id, controlled_node_id) listened_basin_storage = u.storage[listened_node_idx] - reduction_factor = min(listened_basin_storage, 10.0) / 10.0 + reduction_factor = clamp(listened_basin_storage, 0.0, 10.0) / 10.0 else controlled_node_idx = findsorted(outlet.node_id, controlled_node_id) @@ -581,7 +581,7 @@ function continuous_control!( has_index, upstream_basin_idx = id_index(basin.node_id, upstream_node_id) if has_index upstream_basin_storage = u.storage[upstream_basin_idx] - reduction_factor = min(upstream_basin_storage, 10.0) / 10.0 + reduction_factor = clamp(upstream_basin_storage, 0.0, 10.0) / 10.0 else reduction_factor = 1.0 end @@ -704,9 +704,10 @@ Directed graph: outflow is positive! function formulate!( tabulated_rating_curve::TabulatedRatingCurve, p::Parameters, + storage::AbstractVector{Float64}, t::Float64, )::Nothing - (; connectivity) = p + (; basin, connectivity) = p (; graph_flow, flow) = connectivity (; node_id, active, tables) = tabulated_rating_curve for (i, id) in enumerate(node_id) @@ -714,7 +715,11 @@ function formulate!( downstream_ids = outneighbors(graph_flow, id) if active[i] - q = tables[i](get_level(p, upstream_basin_id, t)) + hasindex, basin_idx = id_index(basin.node_id, upstream_basin_id) + @assert hasindex "TabulatedRatingCurve must be downstream of a Basin" + s = storage[basin_idx] + reduction_factor = clamp(s, 0.0, 10.0) / 10.0 + q = reduction_factor * tables[i](get_level(p, upstream_basin_id, t)) else q = 0.0 end @@ -882,7 +887,7 @@ function formulate!( if hasindex # Pumping from basin s = storage[basin_idx] - reduction_factor = min(s, 10.0) / 10.0 + reduction_factor = clamp(s, 0.0, 10.0) / 10.0 q = reduction_factor * rate else # Pumping from level boundary @@ -932,7 +937,7 @@ function formulate_flows!( formulate!(linear_resistance, p, t) formulate!(manning_resistance, p, t) - formulate!(tabulated_rating_curve, p, t) + formulate!(tabulated_rating_curve, p, storage, t) formulate!(flow_boundary, p, t) formulate!(fractional_flow, p) formulate!(pump, p, storage) diff --git a/core/test/run_models.jl b/core/test/run_models.jl index c11b07abd..fa40bfda9 100644 --- a/core/test/run_models.jl +++ b/core/test/run_models.jl @@ -47,7 +47,7 @@ end @test model isa Ribasim.Model @test successful_retcode(model) @test length(model.integrator.p.basin.precipitation) == 4 - @test model.integrator.sol.u[end] ≈ Float32[469.8923, 469.89038, 410.71472, 1427.4194] skip = + @test model.integrator.sol.u[end] ≈ Float32[473.06985, 473.06796, 374.6535, 1428.3226] skip = Sys.isapple() end @@ -80,7 +80,7 @@ end model = Ribasim.run(toml_path) @test model isa Ribasim.Model @test successful_retcode(model) - @test model.integrator.sol.u[end] ≈ Float32[5.949285, 725.9446] skip = Sys.isapple() + @test model.integrator.sol.u[end] ≈ Float32[8.418005, 725.491] skip = Sys.isapple() # the highest level in the dynamic table is updated to 1.2 from the callback @test model.integrator.p.tabulated_rating_curve.tables[end].t[end] == 1.2 end From b5569f792640344013a684fdb2fdc55adaf476ad Mon Sep 17 00:00:00 2001 From: Hofer-Julian <30049909+Hofer-Julian@users.noreply.github.com> Date: Tue, 29 Aug 2023 17:06:22 +0200 Subject: [PATCH 02/13] Run pre-commit_auto_update even less often (#561) The commit waiting times are just too annoying --- .github/workflows/pre-commit_auto_update.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/pre-commit_auto_update.yml b/.github/workflows/pre-commit_auto_update.yml index c1b29fd76..bf4d418ae 100644 --- a/.github/workflows/pre-commit_auto_update.yml +++ b/.github/workflows/pre-commit_auto_update.yml @@ -2,8 +2,8 @@ name: Pre-commit on: schedule: - # every monday morning - - cron: "0 7 * * MON" + # At 03:00 on day 1 of the month + - cron: "0 3 1 * *" # on demand workflow_dispatch: From 91232b96b2ada1f151cedf4bf19874960c29eac5 Mon Sep 17 00:00:00 2001 From: Martijn Visser Date: Tue, 29 Aug 2023 17:37:06 +0200 Subject: [PATCH 03/13] bump ribasim-python version --- python/ribasim/ribasim/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/ribasim/ribasim/__init__.py b/python/ribasim/ribasim/__init__.py index 5695d01ec..0eb79db09 100644 --- a/python/ribasim/ribasim/__init__.py +++ b/python/ribasim/ribasim/__init__.py @@ -1,4 +1,4 @@ -__version__ = "0.2.0" +__version__ = "0.3.0" from ribasim import models, utils From 194c6424d035844de5d08f20699b758b55592cb9 Mon Sep 17 00:00:00 2001 From: Martijn Visser Date: Tue, 29 Aug 2023 17:41:11 +0200 Subject: [PATCH 04/13] set core version to 0.2.0 (#558) Won't be registered, but good to bump this for the release. Co-authored-by: Hofer-Julian <30049909+Hofer-Julian@users.noreply.github.com> --- core/Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/core/Project.toml b/core/Project.toml index 431f1b61a..0a3bfd55b 100644 --- a/core/Project.toml +++ b/core/Project.toml @@ -1,7 +1,7 @@ name = "Ribasim" uuid = "aac5e3d9-0b8f-4d4f-8241-b1a7a9632635" authors = ["Deltares and contributors"] -version = "0.1.0" +version = "0.2.0" [deps] Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45" From 74e0fba6ead688e1128b1a2688b7fe5e7b4b1317 Mon Sep 17 00:00:00 2001 From: Martijn Visser Date: Fri, 1 Sep 2023 13:02:46 +0200 Subject: [PATCH 05/13] allow overwriting open GeoPackage (#569) This works fine having the file open in QGIS on Windows. We write to another file, and at the end move this file over the destination file. Fixes #563 --- python/ribasim/ribasim/geometry/edge.py | 12 +++------ python/ribasim/ribasim/geometry/node.py | 16 +++-------- python/ribasim/ribasim/input_base.py | 30 +++++++++------------ python/ribasim/ribasim/model.py | 35 +++++++++++++++++++------ 4 files changed, 45 insertions(+), 48 deletions(-) diff --git a/python/ribasim/ribasim/geometry/edge.py b/python/ribasim/ribasim/geometry/edge.py index 14cd47dda..9f191a237 100644 --- a/python/ribasim/ribasim/geometry/edge.py +++ b/python/ribasim/ribasim/geometry/edge.py @@ -1,4 +1,3 @@ -from pathlib import Path from typing import Any, Dict, Union import geopandas as gpd @@ -43,20 +42,15 @@ class Config: def _layername(cls, field) -> str: return cls.get_input_type() - def write(self, directory: FilePath, modelname: str) -> None: + def write_layer(self, path: FilePath) -> None: """ Write the contents of the input to a GeoPackage. - The Geopackage will be written in ``directory`` and will be be named - ``{modelname}.gpkg``. - Parameters ---------- - directory : FilePath - modelname : str + path : FilePath """ self.sort() - directory = Path(directory) dataframe = self.static name = self._layername(dataframe) @@ -65,7 +59,7 @@ def write(self, directory: FilePath, modelname: str) -> None: gdf = gdf.set_geometry("geometry") else: gdf["geometry"] = None - gdf.to_file(directory / f"{modelname}.gpkg", layer=name, driver="GPKG") + gdf.to_file(path, layer=name, driver="GPKG") return diff --git a/python/ribasim/ribasim/geometry/node.py b/python/ribasim/ribasim/geometry/node.py index 9f3b23658..998fdaeda 100644 --- a/python/ribasim/ribasim/geometry/node.py +++ b/python/ribasim/ribasim/geometry/node.py @@ -1,4 +1,3 @@ -from pathlib import Path from typing import Any, Dict, Union import geopandas as gpd @@ -75,31 +74,22 @@ def get_node_ids_and_types(*nodes): return node_id, node_type - def write(self, directory: FilePath, modelname: str) -> None: + def write_layer(self, path: FilePath) -> None: """ Write the contents of the input to a GeoPackage. - The Geopackage will be written in ``directory`` and will be be named - ``{modelname}.gpkg``. - Parameters ---------- - directory : FilePath - modelname : str + path : FilePath """ self.sort() - directory = Path(directory) dataframe = self.static name = self._layername(dataframe) gdf = gpd.GeoDataFrame(data=dataframe) gdf = gdf.set_geometry("geometry") - gdf.to_file( - directory / f"{modelname}.gpkg", - layer=name, - driver="GPKG", - ) + gdf.to_file(path, layer=name, driver="GPKG") return diff --git a/python/ribasim/ribasim/input_base.py b/python/ribasim/ribasim/input_base.py index 05fe591f3..f4acc07d6 100644 --- a/python/ribasim/ribasim/input_base.py +++ b/python/ribasim/ribasim/input_base.py @@ -1,6 +1,6 @@ import re import textwrap -from pathlib import Path +from contextlib import closing from sqlite3 import Connection, connect from typing import Any, Dict, Set, Union @@ -22,11 +22,11 @@ def esc_id(identifier: str) -> str: def exists(connection: Connection, name: str) -> bool: """Check if a table exists in a SQLite database.""" - cursor = connection.cursor() - cursor.execute( - "SELECT name FROM sqlite_master WHERE type='table' AND name=?", (name,) - ) - result = cursor.fetchone() + with closing(connection.cursor()) as cursor: + cursor.execute( + "SELECT name FROM sqlite_master WHERE type='table' AND name=?", (name,) + ) + result = cursor.fetchone() return result is not None @@ -88,22 +88,16 @@ def get_node_IDs(self) -> Set[int]: def _layername(cls, field) -> str: return f"{cls.get_input_type()}{delimiter}{field}" - def write(self, directory: FilePath, modelname: str) -> None: + def write_table(self, connection: Connection) -> None: """ Write the contents of the input to a GeoPackage. - The Geopackage will be written in ``directory`` and will be be named - ``{modelname}.gpkg``. - Parameters ---------- - directory : FilePath - Path to the directory where to write the files. - modelname : str - Name of the model, used as a file name. + connection : Connection + SQLite connection to the GeoPackage. """ self.sort() - directory = Path(directory) sql = "INSERT INTO gpkg_contents (table_name, data_type, identifier) VALUES (?, ?, ?)" for field in self.fields(): dataframe = getattr(self, field) @@ -111,9 +105,9 @@ def write(self, directory: FilePath, modelname: str) -> None: continue name = self._layername(field) - with connect(directory / f"{modelname}.gpkg") as connection: - dataframe.to_sql(name, connection, index=False, if_exists="replace") - connection.execute(sql, (name, "attributes", name)) + dataframe.to_sql(name, connection, index=False, if_exists="replace") + with closing(connection.cursor()) as cursor: + cursor.execute(sql, (name, "attributes", name)) return diff --git a/python/ribasim/ribasim/model.py b/python/ribasim/ribasim/model.py index c08588dc7..ec4458219 100644 --- a/python/ribasim/ribasim/model.py +++ b/python/ribasim/ribasim/model.py @@ -1,7 +1,10 @@ import datetime import inspect +import shutil +from contextlib import closing from enum import Enum from pathlib import Path +from sqlite3 import connect from typing import Any, List, Optional, Type, Union, cast import matplotlib.pyplot as plt @@ -174,16 +177,32 @@ def _write_toml(self, directory: FilePath): return def _write_tables(self, directory: FilePath) -> None: - """Write the input to GeoPackage and Arrow tables.""" - # avoid adding tables to existing model + """Write the input to GeoPackage tables.""" + # We write all tables to a temporary GeoPackage with a dot prefix, + # and at the end move this over the target file. + # This does not throw a PermissionError if the file is open in QGIS. directory = Path(directory) gpkg_path = directory / f"{self.modelname}.gpkg" - gpkg_path.unlink(missing_ok=True) - - for name in self.fields(): - input_entry = getattr(self, name) - if isinstance(input_entry, TableModel): - input_entry.write(directory, self.modelname) + tempname = "." + self.modelname + temp_path = gpkg_path.with_stem(tempname) + # avoid adding tables to existing model + temp_path.unlink(missing_ok=True) + + # write to GeoPackage using geopandas + self.node.write_layer(temp_path) + self.edge.write_layer(temp_path) + + # write to GeoPackage using sqlite3 + with closing(connect(temp_path)) as connection: + for name in self.fields(): + input_entry = getattr(self, name) + is_geometry = isinstance(input_entry, Node) or isinstance( + input_entry, Edge + ) + if isinstance(input_entry, TableModel) and not is_geometry: + input_entry.write_table(connection) + + shutil.move(temp_path, gpkg_path) return @staticmethod From 6c70c5a4881ca267ba5e2961f3554b236a20fea1 Mon Sep 17 00:00:00 2001 From: Hofer-Julian <30049909+Hofer-Julian@users.noreply.github.com> Date: Fri, 1 Sep 2023 14:45:09 +0200 Subject: [PATCH 06/13] Add ribasim_cli tests (#560) Fixes #339 --- build/ribasim_cli/tests/test_models.py | 31 +++++++++ python/ribasim/tests/conftest.py | 67 +++---------------- .../ribasim_testmodels/basic.py | 3 +- 3 files changed, 41 insertions(+), 60 deletions(-) create mode 100644 build/ribasim_cli/tests/test_models.py diff --git a/build/ribasim_cli/tests/test_models.py b/build/ribasim_cli/tests/test_models.py new file mode 100644 index 000000000..b84600e33 --- /dev/null +++ b/build/ribasim_cli/tests/test_models.py @@ -0,0 +1,31 @@ +import subprocess +from pathlib import Path + +import pytest +import ribasim +import ribasim_testmodels + + +@pytest.mark.parametrize( + "model_constructor", + map(ribasim_testmodels.__dict__.get, ribasim_testmodels.__all__), +) +def test_ribasim_cli(model_constructor, tmp_path): + model = model_constructor() + assert isinstance(model, ribasim.Model) + model.write(tmp_path) + + executable = ( + Path(__file__).parents[2] + / "create_binaries" + / "ribasim_cli" + / "bin" + / "ribasim.exe" + ) + config_file = str(tmp_path / f"{model.modelname}.toml") + result = subprocess.run([executable, config_file]) + + if model.modelname.startswith("invalid"): + assert result.returncode != 0 + else: + assert result.returncode == 0 diff --git a/python/ribasim/tests/conftest.py b/python/ribasim/tests/conftest.py index fbdf80633..877cbfa4a 100644 --- a/python/ribasim/tests/conftest.py +++ b/python/ribasim/tests/conftest.py @@ -2,53 +2,28 @@ import pytest import ribasim -from ribasim_testmodels import ( - backwater_model, - basic_model, - basic_transient_model, - bucket_model, - discrete_control_of_pid_control_model, - dutch_waterways_model, - flow_boundary_time_model, - flow_condition_model, - invalid_discrete_control_model, - invalid_edge_types_model, - invalid_fractional_flow_model, - invalid_qh_model, - level_boundary_condition_model, - level_setpoint_with_minmax_model, - linear_resistance_model, - manning_resistance_model, - misc_nodes_model, - pid_control_equation_model, - pid_control_model, - pump_discrete_control_model, - rating_curve_model, - tabulated_rating_curve_control_model, - tabulated_rating_curve_model, - trivial_model, -) +import ribasim_testmodels # we can't call fixtures directly, so we keep separate versions @pytest.fixture() def basic() -> ribasim.Model: - return basic_model() + return ribasim_testmodels.basic_model() @pytest.fixture() -def basic_transient(basic) -> ribasim.Model: - return basic_transient_model(basic) +def basic_transient() -> ribasim.Model: + return ribasim_testmodels.basic_transient_model() @pytest.fixture() def tabulated_rating_curve() -> ribasim.Model: - return tabulated_rating_curve_model() + return ribasim_testmodels.tabulated_rating_curve_model() @pytest.fixture() def backwater() -> ribasim.Model: - return backwater_model() + return ribasim_testmodels.backwater_model() # write models to disk for Julia tests to use @@ -57,36 +32,10 @@ def backwater() -> ribasim.Model: models = [ model_generator() - for model_generator in ( - backwater_model, - basic_model, - bucket_model, - discrete_control_of_pid_control_model, - dutch_waterways_model, - flow_boundary_time_model, - flow_condition_model, - invalid_discrete_control_model, - invalid_edge_types_model, - invalid_fractional_flow_model, - invalid_qh_model, - level_boundary_condition_model, - level_setpoint_with_minmax_model, - linear_resistance_model, - manning_resistance_model, - misc_nodes_model, - pid_control_equation_model, - pid_control_model, - pump_discrete_control_model, - rating_curve_model, - tabulated_rating_curve_control_model, - tabulated_rating_curve_model, - trivial_model, + for model_generator in map( + ribasim_testmodels.__dict__.get, ribasim_testmodels.__all__ ) ] for model in models: model.write(datadir / model.modelname) - - if model.modelname == "basic": - model = basic_transient_model(model) - model.write(datadir / model.modelname) diff --git a/python/ribasim_testmodels/ribasim_testmodels/basic.py b/python/ribasim_testmodels/ribasim_testmodels/basic.py index 25f6b9797..8bbb72ccc 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/basic.py +++ b/python/ribasim_testmodels/ribasim_testmodels/basic.py @@ -210,9 +210,10 @@ def basic_model() -> ribasim.Model: return model -def basic_transient_model(model) -> ribasim.Model: +def basic_transient_model() -> ribasim.Model: """Update the basic model with transient forcing""" + model = basic_model() time = pd.date_range(model.starttime, model.endtime) day_of_year = time.day_of_year.to_numpy() seconds_per_day = 24 * 60 * 60 From b13162df46ab45a4412fad73b1069622834087b2 Mon Sep 17 00:00:00 2001 From: Bart de Koning <74617371+SouthEndMusic@users.noreply.github.com> Date: Fri, 1 Sep 2023 15:51:07 +0200 Subject: [PATCH 07/13] ForwardDiff.jl automatic Jacobian (#550) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit *Edit: this initial message is outdated, see the discussion below.* It turned out to be quite easy to get `ForwardDiff.jl` to work for the Jacobian of `water_balance!`, it was mainly a matter of replacing some instances of `Float64` with `Number` in function signatures or removing `Float64` altogether (I did this only where necessary). See the example below for POC. ```julia using ForwardDiff using SQLite using Ribasim toml_path = "C:/Users/konin_bt/Ribasim_development/Ribasim/data/basic_transient/basic_transient.toml" cfg = Config(toml_path) gpkg_path = input_path(cfg, cfg.geopackage) db = SQLite.DB(gpkg_path) p = Ribasim.Parameters(db, cfg) t = 0.0 state = load_structvector(db, cfg, BasinStateV1) storage, _ = get_storages_from_levels(p.basin, state.level) integral = zeros(length(p.pid_control.node_id)) u0 = ComponentVector{Float64}(; storage, integral) function water_balance(u::ComponentVector) du = zero(u) water_balance!(du,u,p,t) return du end jac_AD = ForwardDiff.jacobian(water_balance, u0) jac_AN = zero(jac_AD) water_balance_jac!(jac_AN,u0,p,t) println("AD Jacobian:") display(jac_AD) println("\nAN Jacobian:") display(transpose(jac_AN)) ``` output: ``` AD Jacobian: 4×4 Matrix{Float64}: -3.06216e-12 0.0 0.0 0.0 0.0 -2.25533e-7 0.0 0.0 0.0 2.48017e-8 -3.06216e-12 0.0 0.0 4.96033e-8 0.0 -8.26725e-7 AN Jacobian: 4×4 transpose(::Matrix{Float64}) with eltype Float64: 0.0 0.0 0.0 0.0 0.0 -2.2553e-7 0.0 0.0 0.0 2.48017e-8 0.0 0.0 0.0 4.96033e-8 0.0 -8.26722e-7 ``` ~~There is a remaining problem where the AD Jacobian contained NaNs unless I turned on NaN-safe mode as suggested [here](https://juliadiff.org/ForwardDiff.jl/latest/user/advanced.html#Fixing-NaN/Inf-Issues-1).~~ Regarding the current erroring test; why does `Vector{Float64}` not satisfy `Vector{Number}`? --- core/Project.toml | 4 + core/src/Ribasim.jl | 3 +- core/src/bmi.jl | 22 +- core/src/config.jl | 3 +- core/src/create.jl | 68 +- core/src/io.jl | 2 +- core/src/jac.jl | 601 ------------------ core/src/solve.jl | 240 ++++--- core/src/utils.jl | 28 +- core/src/validation.jl | 2 +- core/test/config.jl | 2 +- core/test/control.jl | 2 +- core/test/docs.toml | 3 +- core/test/equations.jl | 4 +- core/test/run_models.jl | 27 +- core/test/utils.jl | 8 +- core/test/validation.jl | 6 +- docs/core/usage.qmd | 6 +- libribasim/Manifest.toml | 7 + python/ribasim/ribasim/model.py | 1 - .../ribasim_testmodels/dutch_waterways.py | 3 - ribasim_cli/Manifest.toml | 7 + 22 files changed, 273 insertions(+), 776 deletions(-) delete mode 100644 core/src/jac.jl create mode 100644 libribasim/Manifest.toml create mode 100644 ribasim_cli/Manifest.toml diff --git a/core/Project.toml b/core/Project.toml index 0a3bfd55b..16d555d3f 100644 --- a/core/Project.toml +++ b/core/Project.toml @@ -17,12 +17,14 @@ Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" Dictionaries = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4" DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e" Legolas = "741b9549-f6ed-4911-9fbf-4a1c0c97f0cd" Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" LoggingExtras = "e6f89c97-d47a-5376-807f-9c37f3926c36" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46" SQLite = "0aa819cd-b072-5ff4-a722-6bc24af294d9" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" @@ -45,11 +47,13 @@ DataStructures = "0.18" Dictionaries = "0.3.25" DiffEqCallbacks = "2.29.1" FiniteDiff = "2.21" +ForwardDiff = "0.10" Graphs = "1.6" IterTools = "1.4" Legolas = "0.5" LoggingExtras = "1" OrdinaryDiffEq = "6.7" +PreallocationTools = "0.4" SQLite = "1.5.1" SciMLBase = "1.60" StructArrays = "0.6.13" diff --git a/core/src/Ribasim.jl b/core/src/Ribasim.jl index 37273b5b6..c1bc06678 100644 --- a/core/src/Ribasim.jl +++ b/core/src/Ribasim.jl @@ -13,12 +13,14 @@ using DataInterpolations: LinearInterpolation, derivative using Dates using DBInterface: execute, prepare using Dictionaries: Indices, Dictionary, gettoken, gettokenvalue, dictionary +using ForwardDiff: pickchunksize using DiffEqCallbacks using Graphs: DiGraph, add_edge!, adjacency_matrix, inneighbors, outneighbors using Legolas: Legolas, @schema, @version, validate, SchemaVersion, declared using Logging: current_logger, min_enabled_level, with_logger using LoggingExtras: EarlyFilteredLogger, LevelOverrideLogger using OrdinaryDiffEq +using PreallocationTools: DiffCache, get_tmp using SciMLBase using SparseArrays using SQLite: SQLite, DB, Query, esc_id @@ -32,7 +34,6 @@ TimerOutputs.complement!() include("validation.jl") include("solve.jl") -include("jac.jl") include("config.jl") using .config include("utils.jl") diff --git a/core/src/bmi.jl b/core/src/bmi.jl index 6837c3605..c8c7db092 100644 --- a/core/src/bmi.jl +++ b/core/src/bmi.jl @@ -78,15 +78,12 @@ function BMI.initialize(T::Type{Model}, config::Config)::Model # default to nearly empty basins, perhaps make required input fill(1.0, n) else - storages, errors = - get_storages_from_levels(parameters.basin, state.level) + storages, errors = get_storages_from_levels(parameters.basin, state.level) if errors - error( - "Encountered errors while parsing the initial levels of basins.", - ) + error("Encountered errors while parsing the initial levels of basins.") end storages - end::Vector{Float64} + end @assert length(storage) == n "Basin / state length differs from number of Basins" # Integrals for PID control integral = zeros(length(parameters.pid_control.node_id)) @@ -97,8 +94,7 @@ function BMI.initialize(T::Type{Model}, config::Config)::Model timespan = (zero(t_end), t_end) jac_prototype = config.solver.sparse ? get_jac_prototype(parameters) : nothing - jac = config.solver.jac ? water_balance_jac! : nothing - RHS = ODEFunction(water_balance!; jac_prototype, jac) + RHS = ODEFunction(water_balance!; jac_prototype) @timeit_debug to "Setup ODEProblem" begin prob = ODEProblem(RHS, u0, timespan, parameters) @@ -113,6 +109,7 @@ function BMI.initialize(T::Type{Model}, config::Config)::Model alg; progress = true, progress_name = "Simulating", + progress_steps = 100, callback, tstops, config.solver.saveat, @@ -249,7 +246,7 @@ function get_value( level_boundary_idx = findsorted(level_boundary.node_id, feature_id) if hasindex_basin - _, level, _ = get_area_and_level(basin, basin_idx, u[basin_idx]) + _, level = get_area_and_level(basin, basin_idx, u[basin_idx]) elseif !isnothing(level_boundary_idx) level = level_boundary.level[level_boundary_idx](t + Δt) else @@ -432,13 +429,16 @@ function set_control_params!(p::Parameters, node_id::Int, control_state::String) for (field, value) in zip(keys(new_state), new_state) if !ismissing(value) - getfield(node, field)[idx] = value + vec = get_tmp(getfield(node, field), 0) + vec[idx] = value end end end "Copy the current flow to the SavedValues" -save_flow(u, t, integrator) = copy(nonzeros(integrator.p.connectivity.flow)) +function save_flow(u, t, integrator) + copy(nonzeros(get_tmp(integrator.p.connectivity.flow, u))) +end "Load updates from 'Basin / time' into the parameters" function update_basin(integrator)::Nothing diff --git a/core/src/config.jl b/core/src/config.jl index 95edd1a01..fee738f45 100644 --- a/core/src/config.jl +++ b/core/src/config.jl @@ -74,8 +74,7 @@ const nodetypes = collect(keys(nodekinds)) reltol::Float64 = 1e-3 maxiters::Int = 1e9 sparse::Bool = true - jac::Bool = true - autodiff::Bool = false + autodiff::Bool = true end @enum Compression begin diff --git a/core/src/create.jl b/core/src/create.jl index 467129663..8221f4504 100644 --- a/core/src/create.jl +++ b/core/src/create.jl @@ -40,7 +40,7 @@ function parse_static_and_time( # If the type is a union, then the associated parameter is optional and # the type is of the form Union{Missing,ActualType} parameter_type = if parameter_name in time_interpolatables - LinearInterpolation + ScalarInterpolation elseif isa(parameter_type, Union) nonmissingtype(parameter_type) else @@ -190,7 +190,7 @@ function static_and_time_node_ids( return static_node_ids, time_node_ids, node_ids, !errors end -function Connectivity(db::DB)::Connectivity +function Connectivity(db::DB, config::Config, chunk_size::Int)::Connectivity if !valid_edge_types(db) error("Invalid edge types found.") end @@ -204,6 +204,10 @@ function Connectivity(db::DB)::Connectivity flow = adjacency_matrix(graph_flow, Float64) nonzeros(flow) .= 0.0 + if config.solver.autodiff + flow = DiffCache(flow, chunk_size) + end + return Connectivity( graph_flow, graph_control, @@ -228,7 +232,7 @@ function LinearResistance(db::DB, config::Config)::LinearResistance return LinearResistance( parsed_parameters.node_id, - parsed_parameters.active, + BitVector(parsed_parameters.active), parsed_parameters.resistance, parsed_parameters.control_mapping, ) @@ -309,7 +313,7 @@ function ManningResistance(db::DB, config::Config)::ManningResistance return ManningResistance( parsed_parameters.node_id, - parsed_parameters.active, + BitVector(parsed_parameters.active), parsed_parameters.length, parsed_parameters.manning_n, parsed_parameters.profile_width, @@ -398,7 +402,7 @@ function FlowBoundary(db::DB, config::Config)::FlowBoundary return FlowBoundary(node_ids, parsed_parameters.active, parsed_parameters.flow_rate) end -function Pump(db::DB, config::Config)::Pump +function Pump(db::DB, config::Config, chunk_size::Int)::Pump static = load_structvector(db, config, PumpStaticV1) defaults = (; min_flow_rate = 0.0, max_flow_rate = NaN, active = true) parsed_parameters, valid = parse_static_and_time(db, config, "Pump"; static, defaults) @@ -408,10 +412,17 @@ function Pump(db::DB, config::Config)::Pump error("Errors occurred when parsing Pump data.") end + # If flow rate is set by PID control, it is part of the AD Jacobian computations + flow_rate = if config.solver.autodiff + DiffCache(parsed_parameters.flow_rate, chunk_size) + else + parsed_parameters.flow_rate + end + return Pump( parsed_parameters.node_id, - parsed_parameters.active, - parsed_parameters.flow_rate, + BitVector(parsed_parameters.active), + flow_rate, parsed_parameters.min_flow_rate, parsed_parameters.max_flow_rate, parsed_parameters.control_mapping, @@ -419,7 +430,7 @@ function Pump(db::DB, config::Config)::Pump ) end -function Outlet(db::DB, config::Config)::Outlet +function Outlet(db::DB, config::Config, chunk_size::Int)::Outlet static = load_structvector(db, config, OutletStaticV1) defaults = (; min_flow_rate = 0.0, max_flow_rate = NaN, active = true) parsed_parameters, valid = parse_static_and_time(db, config, "Outlet"; static, defaults) @@ -429,10 +440,17 @@ function Outlet(db::DB, config::Config)::Outlet error("Errors occurred when parsing Outlet data.") end + # If flow rate is set by PID control, it is part of the AD Jacobian computations + flow_rate = if config.solver.autodiff + DiffCache(parsed_parameters.flow_rate, chunk_size) + else + parsed_parameters.flow_rate + end + return Outlet( parsed_parameters.node_id, - parsed_parameters.active, - parsed_parameters.flow_rate, + BitVector(parsed_parameters.active), + flow_rate, parsed_parameters.min_flow_rate, parsed_parameters.max_flow_rate, parsed_parameters.control_mapping, @@ -445,12 +463,16 @@ function Terminal(db::DB, config::Config)::Terminal return Terminal(static.node_id) end -function Basin(db::DB, config::Config)::Basin +function Basin(db::DB, config::Config, chunk_size::Int)::Basin node_id = get_ids(db, "Basin") n = length(node_id) current_level = zeros(n) current_area = zeros(n) - current_darea = zeros(n) + + if config.solver.autodiff + current_level = DiffCache(current_level, chunk_size) + current_area = DiffCache(current_area, chunk_size) + end precipitation = fill(NaN, length(node_id)) potential_evaporation = fill(NaN, length(node_id)) @@ -476,7 +498,6 @@ function Basin(db::DB, config::Config)::Basin infiltration, current_level, current_area, - current_darea, area, level, storage, @@ -529,7 +550,7 @@ function DiscreteControl(db::DB, config::Config)::DiscreteControl ) end -function PidControl(db::DB, config::Config)::PidControl +function PidControl(db::DB, config::Config, chunk_size::Int)::PidControl static = load_structvector(db, config, PidControlStaticV1) time = load_structvector(db, config, PidControlTimeV1) @@ -550,6 +571,10 @@ function PidControl(db::DB, config::Config)::PidControl pid_error = zeros(length(node_ids)) + if config.solver.autodiff + pid_error = DiffCache(pid_error, chunk_size) + end + # Combine PID parameters into one vector interpolation object pid_parameters = VectorInterpolation[] (; proportional, integral, derivative) = parsed_parameters @@ -578,7 +603,7 @@ function PidControl(db::DB, config::Config)::PidControl return PidControl( node_ids, - parsed_parameters.active, + BitVector(parsed_parameters.active), parsed_parameters.listen_node_id, parsed_parameters.target, pid_parameters, @@ -588,7 +613,10 @@ function PidControl(db::DB, config::Config)::PidControl end function Parameters(db::DB, config::Config)::Parameters - connectivity = Connectivity(db) + n_states = length(get_ids(db, "Basin")) + length(get_ids(db, "PidControl")) + chunk_size = pickchunksize(n_states) + + connectivity = Connectivity(db, config, chunk_size) linear_resistance = LinearResistance(db, config) manning_resistance = ManningResistance(db, config) @@ -596,13 +624,13 @@ function Parameters(db::DB, config::Config)::Parameters fractional_flow = FractionalFlow(db, config) level_boundary = LevelBoundary(db, config) flow_boundary = FlowBoundary(db, config) - pump = Pump(db, config) - outlet = Outlet(db, config) + pump = Pump(db, config, chunk_size) + outlet = Outlet(db, config, chunk_size) terminal = Terminal(db, config) discrete_control = DiscreteControl(db, config) - pid_control = PidControl(db, config) + pid_control = PidControl(db, config, chunk_size) - basin = Basin(db, config) + basin = Basin(db, config, chunk_size) p = Parameters( config.starttime, diff --git a/core/src/io.jl b/core/src/io.jl index 0cec18885..4be67f2db 100644 --- a/core/src/io.jl +++ b/core/src/io.jl @@ -179,7 +179,7 @@ function write_flow_output(model::Model, compress) (; t, saveval) = saved_flow (; connectivity) = integrator.p - I, J, _ = findnz(connectivity.flow) + I, J, _ = findnz(get_tmp(connectivity.flow, integrator.u)) unique_edge_ids = [connectivity.edge_ids_flow[ij] for ij in zip(I, J)] nflow = length(I) ntsteps = length(t) diff --git a/core/src/jac.jl b/core/src/jac.jl deleted file mode 100644 index 3e3e12143..000000000 --- a/core/src/jac.jl +++ /dev/null @@ -1,601 +0,0 @@ -""" -The Jacobian is a n x n matrix where n is the number of basins plus the number of -PidControl nodes. Each basin has a storage state, and each PidControl node has an error integral -state. If we write `water_balance!` as `f(u, p(t), t)` where u is the vector of all states, then -`J[i,j] = ∂f_j/∂u_i`. f_j dictates the time derivative of state j. - -For more on the sparsity see [`get_jac_prototype`](@ref). -""" -function water_balance_jac!( - J::AbstractMatrix, - u::ComponentVector{Float64}, - p::Parameters, - t, -)::Nothing - (; basin) = p - J .= 0.0 - - # Ensures current_* vectors are current - set_current_basin_properties!(basin, u.storage, t) - - for nodefield in nodefields(p) - if nodefield == :pid_control - continue - end - - formulate_jac!(getfield(p, nodefield), J, u, p, t) - end - - # PID control must be done last - formulate_jac!(p.pid_control, J, u, p, t) - - return nothing -end - -""" -The contributions of LinearResistance nodes to the Jacobian. -""" -function formulate_jac!( - linear_resistance::LinearResistance, - J::AbstractMatrix, - u::ComponentVector{Float64}, - p::Parameters, - t::Float64, -)::Nothing - (; basin, connectivity) = p - (; active, resistance, node_id) = linear_resistance - (; graph_flow) = connectivity - - for (id, isactive, R) in zip(node_id, active, resistance) - - # Inactive nodes do not contribute - if !isactive - continue - end - - id_in = only(inneighbors(graph_flow, id)) - id_out = only(outneighbors(graph_flow, id)) - - has_index_in, idx_in = id_index(basin.node_id, id_in) - has_index_out, idx_out = id_index(basin.node_id, id_out) - - if has_index_in - area_in = basin.current_area[idx_in] - term_in = 1 / (area_in * R) - J[idx_in, idx_in] -= term_in - end - - if has_index_out - area_out = basin.current_area[idx_out] - term_out = 1 / (area_out * R) - J[idx_out, idx_out] -= term_out - end - - if has_index_in && has_index_out - J[idx_in, idx_out] += term_out - J[idx_out, idx_in] += term_in - end - end - return nothing -end - -""" -The contributions of ManningResistance nodes to the Jacobian. -""" -function formulate_jac!( - manning_resistance::ManningResistance, - J::AbstractMatrix, - u::ComponentVector{Float64}, - p::Parameters, - t::Float64, -)::Nothing - (; basin, connectivity) = p - (; node_id, active, length, manning_n, profile_width, profile_slope) = - manning_resistance - (; graph_flow) = connectivity - - for (i, id) in enumerate(node_id) - - # Inactive nodes do not contribute - if !active[i] - continue - end - - #TODO: this was copied from formulate! for the manning_resistance, - # maybe put in separate function - basin_a_id = only(inneighbors(graph_flow, id)) - basin_b_id = only(outneighbors(graph_flow, id)) - - h_a = get_level(p, basin_a_id, t) - h_b = get_level(p, basin_b_id, t) - bottom_a, bottom_b = basin_bottoms(basin, basin_a_id, basin_b_id, id) - slope = profile_slope[i] - width = profile_width[i] - n = manning_n[i] - L = length[i] - - Δh = h_a - h_b - q_sign = sign(Δh) - - # Average d, A, R - d_a = h_a - bottom_a - d_b = h_b - bottom_b - d = 0.5 * (d_a + d_b) - - A_a = width * d + slope * d_a^2 - A_b = width * d + slope * d_b^2 - A = 0.5 * (A_a + A_b) - - slope_unit_length = sqrt(slope^2 + 1.0) - P_a = width + 2.0 * d_a * slope_unit_length - P_b = width + 2.0 * d_b * slope_unit_length - R_h_a = A_a / P_a - R_h_b = A_b / P_b - R_h = 0.5 * (R_h_a + R_h_b) - - k = 1000.0 - kΔh = k * Δh - atankΔh = atan(k * Δh) - ΔhatankΔh = Δh * atankΔh - R_hpow = R_h^(2 / 3) - root = sqrt(2 / π * ΔhatankΔh) - - id_in = only(inneighbors(graph_flow, id)) - id_out = only(outneighbors(graph_flow, id)) - - has_index_in, idx_in = id_index(basin.node_id, id_in) - has_index_out, idx_out = id_index(basin.node_id, id_out) - - if has_index_in - basin_in_area = basin.current_area[idx_in] - ∂A_a = (width + 2 * slope * d_a) / basin_in_area - ∂A = 0.5 * ∂A_a - ∂P_a = 2 * slope_unit_length / basin_in_area - ∂R_h_a = (P_a * ∂A_a - A_a * ∂P_a) / P_a^2 - ∂R_h_b = width / (2 * basin_in_area * P_b) - ∂R_h = 0.5 * (∂R_h_a + ∂R_h_b) - # This float exact comparison is deliberate since `sqrt_contribution` has a - # removable singularity, i.e. it doesn't exist at $\Delta h = 0$ because of - # division by zero but the limit Δh → 0 does exist and is equal to the given value. - if Δh == 0 - sqrt_contribution = 2 / (sqrt(2 * π) * basin_in_area) - else - sqrt_contribution = - (atankΔh + kΔh / (1 + kΔh^2)) / - (basin_in_area * sqrt(2 * π * ΔhatankΔh)) - end - # term_in = q * (∂A / A + ∂R_h / R_h + sqrt_contribution) - term_in = - q_sign * ( - ∂A * R_hpow * root + - A * R_hpow * ∂R_h / R_h * root + - A * R_hpow * sqrt_contribution - ) / (n * sqrt(L)) - J[idx_in, idx_in] -= term_in - end - - if has_index_out - basin_out_area = basin.current_area[idx_out] - ∂A_b = (width + 2 * slope * d_b) / basin_out_area - ∂A = 0.5 * ∂A_b - ∂P_b = 2 * slope_unit_length / basin_out_area - ∂R_h_b = (P_b * ∂A_b - A_b * ∂P_b) / P_b^2 - ∂R_h_b = width / (2 * basin_out_area * P_b) - ∂R_h = 0.5 * (∂R_h_b + ∂R_h_a) - # This float exact comparison is deliberate since `sqrt_contribution` has a - # removable singularity, i.e. it doesn't exist at $\Delta h = 0$ because of - # division by zero but the limit Δh → 0 does exist and is equal to the given value. - if Δh == 0 - sqrt_contribution = 2 / (sqrt(2 * π) * basin_out_area) - else - sqrt_contribution = - (atankΔh + kΔh / (1 + kΔh^2)) / - (basin_out_area * sqrt(2 * π * ΔhatankΔh)) - end - # term_out = q * (∂A / A + ∂R_h / R_h + sqrt_contribution) - term_out = - q_sign * ( - ∂A * R_hpow * root + - A * R_hpow * ∂R_h / R_h * root + - A * R_hpow * sqrt_contribution - ) / (n * sqrt(L)) - J[idx_out, idx_out] -= term_out - end - - if has_index_in && has_index_out - J[idx_in, idx_out] += term_out - J[idx_out, idx_in] += term_in - end - end - return nothing -end - -""" -The contributions of Pump and Outlet nodes to the Jacobian. -""" -function formulate_jac!( - node::Union{Pump, Outlet}, - J::AbstractMatrix, - u::ComponentVector{Float64}, - p::Parameters, - t::Float64, -)::Nothing - (; basin, fractional_flow, connectivity) = p - (; active, node_id, flow_rate, is_pid_controlled) = node - - (; graph_flow) = connectivity - - for (i, id) in enumerate(node_id) - - # Inactive nodes do not contribute - if !active[i] - continue - end - - if is_pid_controlled[i] - continue - end - - id_in = only(inneighbors(graph_flow, id)) - - # For inneighbors only directly connected basins give a contribution - has_index_in, idx_in = id_index(basin.node_id, id_in) - - # For outneighbors there can be directly connected basins - # or basins connected via a fractional flow - # (but not both at the same time!) - if has_index_in - s = u.storage[idx_in] - - if s < 10.0 - dq = flow_rate[i] / 10.0 - - J[idx_in, idx_in] -= dq - - has_index_out, idx_out = id_index(basin.node_id, id_in) - - idxs_fractional_flow, idxs_out = get_fractional_flow_connected_basins( - id, - basin, - fractional_flow, - graph_flow, - ) - - # Assumes either one outneighbor basin or one or more - # outneighbor fractional flows - if isempty(idxs_out) - id_out = only(outneighbors(graph_flow, id)) - has_index_out, idx_out = id_index(basin.node_id, id_out) - - if has_index_out - J[idx_in, idx_out] += dq - end - else - for (idx_fractional_flow, idx_out) in - zip(idxs_fractional_flow, idxs_out) - J[idx_in, idx_out] += - dq * fractional_flow.fraction[idx_fractional_flow] - end - end - end - end - end - return nothing -end - -""" -The contributions of TabulatedRatingCurve nodes to the Jacobian. -""" -function formulate_jac!( - tabulated_rating_curve::TabulatedRatingCurve, - J::AbstractMatrix, - u::ComponentVector{Float64}, - p::Parameters, - t::Float64, -)::Nothing - (; basin, fractional_flow, connectivity) = p - (; node_id, active, tables) = tabulated_rating_curve - (; graph_flow) = connectivity - - for (i, id) in enumerate(node_id) - if !active[i] - continue - end - - id_in = only(inneighbors(graph_flow, id)) - - # For inneighbors only directly connected basins give a contribution - has_index_in, idx_in = id_index(basin.node_id, id_in) - - # For outneighbors there can be directly connected basins - # or basins connected via a fractional flow - if has_index_in - # Computing this slope here is silly, - # should eventually be computed pre-simulation and cached! - table = tables[i] - level = basin.current_level[idx_in] - slope = scalar_interpolation_derivative( - table, - level; - extrapolate_up_constant = false, - ) - - dq = slope / basin.current_area[idx_in] - - J[idx_in, idx_in] -= dq - - idxs_fractional_flow, idxs_out = - get_fractional_flow_connected_basins(id, basin, fractional_flow, graph_flow) - - # Assumes either one outneighbor basin or one or more - # outneighbor fractional flows - if isempty(idxs_out) - id_out = only(outneighbors(graph_flow, id)) - has_index_out, idx_out = id_index(basin.node_id, id_out) - - if has_index_out - J[idx_in, idx_out] += dq - end - else - for (idx_fractional_flow, idx_out) in zip(idxs_fractional_flow, idxs_out) - J[idx_in, idx_out] += dq * fractional_flow.fraction[idx_fractional_flow] - end - end - end - end - return nothing -end - -""" -The contributions of PidControl nodes to the Jacobian. -""" -function formulate_jac!( - pid_control::PidControl, - J::AbstractMatrix, - u::ComponentVector{Float64}, - p::Parameters, - t::Float64, -)::Nothing - (; basin, connectivity, pump, outlet) = p - (; node_id, active, listen_node_id, pid_params, target, error) = pid_control - min_flow_rate_pump = pump.min_flow_rate - max_flow_rate_pump = pump.max_flow_rate - min_flow_rate_outlet = outlet.min_flow_rate - max_flow_rate_outlet = outlet.max_flow_rate - (; graph_flow, graph_control) = connectivity - - pid_params_interpolated = [params(t) for params in pid_params] - derivative = [params[3] for params in pid_params_interpolated] - - get_error!(pid_control, p, t) - - n_basins = length(basin.node_id) - integral_value = u.integral - - if any(.!isnan.(derivative)) - # Calling water_balance is expensive, but it is a sure way of getting - # the proper du for the pid controlled basins - # TODO: Do not allocate new memory here, make field of struct - du = zero(u) - water_balance!(du, u, p, t) - end - - for (i, id) in enumerate(node_id) - if !active[i] - continue - end - - controlled_node_id = only(outneighbors(graph_control, id)) - controls_pump = insorted(controlled_node_id, pump.node_id) - - if !controls_pump - if !insorted(controlled_node_id, outlet.node_id) - error( - "Node #$controlled_node_id controlled by PidControl #$id is neither a Pump nor an Outlet.", - ) - end - end - - listened_node_id = listen_node_id[i] - _, listened_node_idx = id_index(basin.node_id, listened_node_id) - listen_area = basin.current_area[listened_node_idx] - - if controls_pump - controlled_node_idx = findsorted(pump.node_id, controlled_node_id) - - listened_basin_storage = u.storage[listened_node_idx] - reduction_factor = clamp(listened_basin_storage, 0.0, 10.0) / 10.0 - else - controlled_node_idx = findsorted(outlet.node_id, controlled_node_id) - - # Upstream node of outlet does not have to be a basin - upstream_node_id = only(inneighbors(graph_flow, controlled_node_id)) - has_upstream_index, upstream_basin_idx = - id_index(basin.node_id, upstream_node_id) - if has_upstream_index - upstream_basin_storage = u.storage[upstream_basin_idx] - reduction_factor = clamp(upstream_basin_storage, 0.0, 10.0) / 10.0 - else - reduction_factor = 1.0 - end - end - - K_p, K_i, K_d = pid_params_interpolated[i] - - if !iszero(K_d) - if controls_pump - D = 1.0 - K_d * reduction_factor / listen_area - else - D = 1.0 + K_d * reduction_factor / listen_area - end - else - D = 1.0 - end - - E = 0.0 - - if !iszero(K_p) - E += K_p * error[i] - end - - if !iszero(K_i) - E += K_i * integral_value[i] - end - - if !iszero(K_d) - dtarget_level = scalar_interpolation_derivative(target[i], t) - du_listened_basin_old = du.storage[listened_node_idx] - E += K_d * (dtarget_level - du_listened_basin_old / listen_area) - end - - # Clip values outside pump flow rate bounds - if controls_pump - min_flow_rate = min_flow_rate_pump - max_flow_rate = max_flow_rate_pump - else - min_flow_rate = min_flow_rate_outlet - max_flow_rate = max_flow_rate_outlet - end - - flow_rate = reduction_factor * E / D - was_clipped = false - - if flow_rate < min_flow_rate[controlled_node_idx] - was_clipped = true - flow_rate = min_flow_rate[controlled_node_idx] - end - - if !isnan(max_flow_rate[controlled_node_idx]) - if flow_rate > max_flow_rate[controlled_node_idx] - was_clipped = true - flow_rate = max_flow_rate[controlled_node_idx] - end - end - - # PID control integral state - pid_state_idx = n_basins + i - J[pid_state_idx, listened_node_idx] -= 1 / listen_area - - # If the flow rate is clipped to one of the bounds it does - # not change with storages and thus doesn't contribute to the - # Jacobian - if was_clipped - continue - end - - # Only in this case the reduction factor has a non-zero derivative - reduction_factor_regime = if controls_pump - listened_basin_storage < 10.0 - else - if has_upstream_index - upstream_basin_storage < 10.0 - else - false - end - end - - # Computing D and E derivatives - if !iszero(K_d) - darea = basin.current_darea[listened_node_idx] - - dD_du_listen = reduction_factor * darea / (listen_area^2) - - if reduction_factor_regime - if controls_pump - dD_du_listen -= 0.1 / darea - else - dD_du_upstream = -0.1 * K_d / area - end - else - dD_du_upstream = 0.0 - end - - dD_du_listen *= K_d - - dE_du_listen = - -K_d * ( - listen_area * J[listened_node_idx, listened_node_idx] - - du.storage[listened_node_idx] * darea - ) / (listen_area^2) - else - dD_du_listen = 0.0 - dD_du_upstream = 0.0 - dE_du_listen = 0.0 - end - - if !iszero(K_p) - dE_du_listen -= K_p / listen_area - end - - dQ_du_listen = reduction_factor * (D * dE_du_listen - E * dD_du_listen) / (D^2) - - if controls_pump && reduction_factor_regime - dQ_du_listen += 0.1 * E / D - end - - if controls_pump - J[listened_node_idx, listened_node_idx] -= dQ_du_listen - - downstream_node_id = only(outneighbors(graph_flow, controlled_node_id)) - has_downstream_index, downstream_node_idx = - id_index(basin.node_id, downstream_node_id) - - if has_downstream_index - J[listened_node_idx, downstream_node_idx] += dQ_du_listen - end - else - J[listened_node_idx, listened_node_idx] += dQ_du_listen - - if has_upstream_index - J[listened_node_idx, upstream_basin_idx] -= dQ_du_listen - end - end - - if !controls_pump - if !isnan(K_d) && has_upstream_index - dE_du_upstream = -K_d * J[upstream_basin_idx, listened_node_idx] / area - - dQ_du_upstream = - reduction_factor * (D * dE_du_upstream - E * dD_du_upstream) / (D^2) - - if reduction_factor_regime - dQ_du_upstream += 0.1 * E / D - end - - J[upstream_basin_idx, listened_node_idx] += dQ_du_upstream - J[upstream_basin_idx, upstream_basin_idx] -= dQ_du_upstream - end - end - end - return nothing -end - -""" -Method for nodes that do not contribute to the Jacobian -""" -function formulate_jac!( - node::AbstractParameterNode, - J::AbstractMatrix, - u::ComponentVector{Float64}, - p::Parameters, - t::Float64, -)::Nothing - node_type = nameof(typeof(node)) - - if !isa( - node, - Union{ - Basin, - DiscreteControl, - FlowBoundary, - FractionalFlow, - LevelBoundary, - Terminal, - }, - ) - error( - "It is not specified how nodes of type $node_type contribute to the Jacobian.", - ) - end - return nothing -end diff --git a/core/src/solve.jl b/core/src/solve.jl index c1d040adf..cca2d16e1 100644 --- a/core/src/solve.jl +++ b/core/src/solve.jl @@ -11,11 +11,17 @@ graph_flow, graph_control: directed graph with vertices equal to ids flow: store the flow on every flow edge edge_ids_flow, edge_ids_control: get the external edge id from (src, dst) edge_connection_type_flow, edge_connection_types_control: get (src_node_type, dst_node_type) from edge id + +if autodiff + T = DiffCache{SparseArrays.SparseMatrixCSC{Float64, Int64}, Vector{Float64}} +else + T = SparseMatrixCSC{Float64, Int} +end """ -struct Connectivity +struct Connectivity{T} graph_flow::DiGraph{Int} graph_control::DiGraph{Int} - flow::SparseMatrixCSC{Float64, Int} + flow::T edge_ids_flow::Dictionary{Tuple{Int, Int}, Int} edge_ids_flow_inv::Dictionary{Int, Tuple{Int, Int}} edge_ids_control::Dictionary{Tuple{Int, Int}, Int} @@ -24,13 +30,13 @@ struct Connectivity function Connectivity( graph_flow, graph_control, - flow, + flow::T, edge_ids_flow, edge_ids_flow_inv, edge_ids_control, edge_connection_types_flow, edge_connection_types_control, - ) + ) where {T} invalid_networks = Vector{String}() if !valid_edges(edge_ids_flow, edge_connection_types_flow) @@ -42,7 +48,7 @@ struct Connectivity end if isempty(invalid_networks) - new( + new{T}( graph_flow, graph_control, flow, @@ -71,19 +77,22 @@ Requirements: Type parameter C indicates the content backing the StructVector, which can be a NamedTuple of vectors or Arrow Tables, and is added to avoid type instabilities. The node_id are Indices to support fast lookup of e.g. current_level using ID. + +if autodiff + T = DiffCache{Vector{Float64}} +else + T = Vector{Float64} +end """ -struct Basin{C} <: AbstractParameterNode +struct Basin{T, C} <: AbstractParameterNode node_id::Indices{Int} precipitation::Vector{Float64} potential_evaporation::Vector{Float64} drainage::Vector{Float64} infiltration::Vector{Float64} # cache this to avoid recomputation - current_level::Vector{Float64} - current_area::Vector{Float64} - # The derivative of the area with respect to the level - # used for the analytical Jacobian - current_darea::Vector{Float64} + current_level::T + current_area::T # Discrete values for interpolation area::Vector{Vector{Float64}} level::Vector{Vector{Float64}} @@ -97,17 +106,16 @@ struct Basin{C} <: AbstractParameterNode potential_evaporation, drainage, infiltration, - current_level, - current_area, - current_darea, + current_level::T, + current_area::T, area, level, storage, time::StructVector{BasinForcingV1, C, Int}, - ) where {C} + ) where {T, C} errors = valid_profiles(node_id, level, area) if isempty(errors) - return new{C}( + return new{T, C}( node_id, precipitation, potential_evaporation, @@ -115,7 +123,6 @@ struct Basin{C} <: AbstractParameterNode infiltration, current_level, current_area, - current_darea, area, level, storage, @@ -262,10 +269,10 @@ max_flow_rate: The maximum flow rate of the pump control_mapping: dictionary from (node_id, control_state) to target flow rate is_pid_controlled: whether the flow rate of this pump is governed by PID control """ -struct Pump <: AbstractParameterNode +struct Pump{T} <: AbstractParameterNode node_id::Vector{Int} active::BitVector - flow_rate::Vector{Float64} + flow_rate::T min_flow_rate::Vector{Float64} max_flow_rate::Vector{Float64} control_mapping::Dict{Tuple{Int, String}, NamedTuple} @@ -274,14 +281,14 @@ struct Pump <: AbstractParameterNode function Pump( node_id, active, - flow_rate, + flow_rate::T, min_flow_rate, max_flow_rate, control_mapping, is_pid_controlled, - ) - if valid_flow_rates(node_id, flow_rate, control_mapping, :Pump) - return new( + ) where {T} + if valid_flow_rates(node_id, get_tmp(flow_rate, 0), control_mapping, :Pump) + return new{T}( node_id, active, flow_rate, @@ -305,10 +312,10 @@ max_flow_rate: The maximum flow rate of the outlet control_mapping: dictionary from (node_id, control_state) to target flow rate is_pid_controlled: whether the flow rate of this outlet is governed by PID control """ -struct Outlet <: AbstractParameterNode +struct Outlet{T} <: AbstractParameterNode node_id::Vector{Int} active::BitVector - flow_rate::Vector{Float64} + flow_rate::T min_flow_rate::Vector{Float64} max_flow_rate::Vector{Float64} control_mapping::Dict{Tuple{Int, String}, NamedTuple} @@ -317,14 +324,14 @@ struct Outlet <: AbstractParameterNode function Outlet( node_id, active, - flow_rate, + flow_rate::T, min_flow_rate, max_flow_rate, control_mapping, is_pid_controlled, - ) - if valid_flow_rates(node_id, flow_rate, control_mapping, :Outlet) - return new( + ) where {T} + if valid_flow_rates(node_id, get_tmp(flow_rate, 0), control_mapping, :Outlet) + return new{T}( node_id, active, flow_rate, @@ -383,13 +390,13 @@ pid_params: a vector interpolation for parameters changing over time. where the last three are the coefficients for the PID equation. error: the current error; basin_target - current_level """ -struct PidControl <: AbstractParameterNode +struct PidControl{T} <: AbstractParameterNode node_id::Vector{Int} active::BitVector listen_node_id::Vector{Int} target::Vector{ScalarInterpolation} pid_params::Vector{VectorInterpolation} - error::Vector{Float64} + error::T control_mapping::Dict{Tuple{Int, String}, NamedTuple} end @@ -479,15 +486,17 @@ end function set_current_basin_properties!( basin::Basin, - storage::AbstractVector{Float64}, - t::Real, + current_area, + current_level, + storage::AbstractVector, + t::Float64, )::Nothing for i in eachindex(storage) s = storage[i] - area, level, darea = get_area_and_level(basin, i, s) - basin.current_level[i] = level - basin.current_area[i] = area - basin.current_darea[i] = darea + area, level = get_area_and_level(basin, i, s) + + current_level[i] = level + current_area[i] = area end end @@ -496,15 +505,17 @@ Linearize the evaporation flux when at small water depths Currently at less than 0.1 m. """ function formulate!( - du::AbstractVector{Float64}, + du::AbstractVector, basin::Basin, - storage::AbstractVector{Float64}, - t::Real, + storage::AbstractVector, + current_area, + current_level, + t::Float64, )::Nothing for i in eachindex(storage) # add all precipitation that falls within the profile - level = basin.current_level[i] - area = basin.current_area[i] + level = current_level[i] + area = current_area[i] bottom = basin.level[i][1] fixed_area = basin.area[i][end] @@ -521,26 +532,34 @@ function formulate!( return nothing end -function get_error!(pid_control::PidControl, p::Parameters, t::Float64) +function get_error!( + pid_control::PidControl, + p::Parameters, + current_level, + pid_error, + t::Float64, +) (; basin) = p (; listen_node_id, target) = pid_control - pid_error = pid_control.error - for i in eachindex(listen_node_id) listened_node_id = listen_node_id[i] has_index, listened_node_idx = id_index(basin.node_id, listened_node_id) @assert has_index "Listen node $listened_node_id is not a Basin." - pid_error[i] = target[i](t) - basin.current_level[listened_node_idx] + pid_error[i] = target[i](t) - current_level[listened_node_idx] end end function continuous_control!( - u::ComponentVector{Float64}, - du::ComponentVector{Float64}, + u::ComponentVector, + du::ComponentVector, + current_area, pid_control::PidControl, p::Parameters, - integral_value::SubArray{Float64}, + integral_value::SubArray, + current_level::AbstractVector, + flow::SparseMatrixCSC, + pid_error, t::Float64, )::Nothing (; connectivity, pump, outlet, basin, fractional_flow) = p @@ -548,10 +567,13 @@ function continuous_control!( max_flow_rate_pump = pump.max_flow_rate min_flow_rate_outlet = outlet.min_flow_rate max_flow_rate_outlet = outlet.max_flow_rate - (; graph_control, graph_flow, flow) = connectivity - (; node_id, active, target, pid_params, listen_node_id, error) = pid_control + (; graph_control, graph_flow) = connectivity + (; node_id, active, target, pid_params, listen_node_id) = pid_control - get_error!(pid_control, p, t) + outlet_flow_rate = get_tmp(outlet.flow_rate, u) + pump_flow_rate = get_tmp(pump.flow_rate, u) + + get_error!(pid_control, p, current_level, pid_error, t) for (i, id) in enumerate(node_id) if !active[i] @@ -560,7 +582,7 @@ function continuous_control!( return end - du.integral[i] = error[i] + du.integral[i] = pid_error[i] listened_node_id = listen_node_id[i] _, listened_node_idx = id_index(basin.node_id, listened_node_id) @@ -593,14 +615,14 @@ function continuous_control!( if !iszero(K_d) # dlevel/dstorage = 1/area - area = basin.current_area[listened_node_idx] + area = current_area[listened_node_idx] D = 1.0 - K_d * reduction_factor / area else D = 1.0 end if !iszero(K_p) - flow_rate += reduction_factor * K_p * error[i] / D + flow_rate += reduction_factor * K_p * pid_error[i] / D end if !iszero(K_i) @@ -635,10 +657,10 @@ function continuous_control!( # in formulate!(du, connectivity, basin), but in this function # flows are set so du has to be updated too. if controls_pump - pump.flow_rate[controlled_node_idx] = flow_rate + pump_flow_rate[controlled_node_idx] = flow_rate du.storage[listened_node_idx] -= flow_rate else - outlet.flow_rate[controlled_node_idx] = flow_rate + outlet_flow_rate[controlled_node_idx] = flow_rate du.storage[listened_node_idx] += flow_rate end @@ -678,16 +700,26 @@ end """ Directed graph: outflow is positive! """ -function formulate!(linear_resistance::LinearResistance, p::Parameters, t::Float64)::Nothing +function formulate!( + linear_resistance::LinearResistance, + p::Parameters, + current_level, + flow, + t::Float64, +)::Nothing (; connectivity) = p - (; graph_flow, flow) = connectivity + (; graph_flow) = connectivity (; node_id, active, resistance) = linear_resistance for (i, id) in enumerate(node_id) basin_a_id = only(inneighbors(graph_flow, id)) basin_b_id = only(outneighbors(graph_flow, id)) if active[i] - q = (get_level(p, basin_a_id, t) - get_level(p, basin_b_id, t)) / resistance[i] + q = + ( + get_level(p, basin_a_id, current_level, t) - + get_level(p, basin_b_id, current_level, t) + ) / resistance[i] flow[basin_a_id, id] = q flow[id, basin_b_id] = q else @@ -705,10 +737,13 @@ function formulate!( tabulated_rating_curve::TabulatedRatingCurve, p::Parameters, storage::AbstractVector{Float64}, + current_level, + flow::SparseMatrixCSC, t::Float64, )::Nothing (; basin, connectivity) = p (; graph_flow, flow) = connectivity + (; graph_flow, flow) = connectivity (; node_id, active, tables) = tabulated_rating_curve for (i, id) in enumerate(node_id) upstream_basin_id = only(inneighbors(graph_flow, id)) @@ -719,7 +754,9 @@ function formulate!( @assert hasindex "TabulatedRatingCurve must be downstream of a Basin" s = storage[basin_idx] reduction_factor = clamp(s, 0.0, 10.0) / 10.0 - q = reduction_factor * tables[i](get_level(p, upstream_basin_id, t)) + q = + reduction_factor * + tables[i](get_level(p, upstream_basin_id, current_level, t)) else q = 0.0 end @@ -774,10 +811,12 @@ dry. function formulate!( manning_resistance::ManningResistance, p::Parameters, + current_level, + flow, t::Float64, )::Nothing (; basin, connectivity) = p - (; graph_flow, flow) = connectivity + (; graph_flow) = connectivity (; node_id, active, length, manning_n, profile_width, profile_slope) = manning_resistance for (i, id) in enumerate(node_id) @@ -790,8 +829,8 @@ function formulate!( continue end - h_a = get_level(p, basin_a_id, t) - h_b = get_level(p, basin_b_id, t) + h_a = get_level(p, basin_a_id, current_level, t) + h_b = get_level(p, basin_b_id, current_level, t) bottom_a, bottom_b = basin_bottoms(basin, basin_a_id, basin_b_id, id) slope = profile_slope[i] width = profile_width[i] @@ -817,8 +856,10 @@ function formulate!( R_h_b = A_b / P_b R_h = 0.5 * (R_h_a + R_h_b) k = 1000.0 + # This epsilon makes sure the AD derivative at Δh = 0 does not give NaN + eps = 1e-200 - q = q_sign * A / n * R_h^(2 / 3) * sqrt(Δh / L * 2 / π * atan(k * Δh)) + q = q_sign * A / n * R_h^(2 / 3) * sqrt(Δh / L * 2 / π * atan(k * Δh) + eps) flow[basin_a_id, id] = q flow[id, basin_b_id] = q @@ -826,9 +867,9 @@ function formulate!( return nothing end -function formulate!(fractional_flow::FractionalFlow, p::Parameters)::Nothing +function formulate!(fractional_flow::FractionalFlow, flow, p::Parameters)::Nothing (; connectivity) = p - (; graph_flow, flow) = connectivity + (; graph_flow) = connectivity (; node_id, fraction) = fractional_flow for (i, id) in enumerate(node_id) downstream_id = only(outneighbors(graph_flow, id)) @@ -838,9 +879,9 @@ function formulate!(fractional_flow::FractionalFlow, p::Parameters)::Nothing return nothing end -function formulate!(flow_boundary::FlowBoundary, p::Parameters, t::Float64)::Nothing +function formulate!(flow_boundary::FlowBoundary, p::Parameters, flow, t::Float64)::Nothing (; connectivity) = p - (; graph_flow, flow) = connectivity + (; graph_flow) = connectivity (; node_id, active, flow_rate) = flow_boundary for (i, id) in enumerate(node_id) @@ -862,11 +903,13 @@ end function formulate!( node::Union{Pump, Outlet}, p::Parameters, - storage::AbstractVector{Float64}, + flow, + storage::AbstractVector, )::Nothing (; connectivity, basin) = p - (; graph_flow, flow) = connectivity + (; graph_flow) = connectivity (; node_id, active, flow_rate, is_pid_controlled) = node + flow_rate = get_tmp(flow_rate, storage) for (id, isactive, rate, pid_controlled) in zip(node_id, active, flow_rate, is_pid_controlled) src_id = only(inneighbors(graph_flow, id)) @@ -901,14 +944,15 @@ function formulate!( end function formulate!( - du::ComponentVector{Float64}, + du::ComponentVector, connectivity::Connectivity, + flow::SparseMatrixCSC, basin::Basin, )::Nothing # loop over basins # subtract all outgoing flows # add all ingoing flows - (; graph_flow, flow) = connectivity + (; graph_flow) = connectivity for (i, basin_id) in enumerate(basin.node_id) for in_id in inneighbors(graph_flow, basin_id) du[i] += flow[in_id, basin_id] @@ -922,7 +966,9 @@ end function formulate_flows!( p::Parameters, - storage::AbstractVector{Float64}, + storage::AbstractVector, + current_level, + flow::SparseMatrixCSC, t::Float64, )::Nothing (; @@ -935,13 +981,13 @@ function formulate_flows!( outlet, ) = p - formulate!(linear_resistance, p, t) - formulate!(manning_resistance, p, t) - formulate!(tabulated_rating_curve, p, storage, t) - formulate!(flow_boundary, p, t) - formulate!(fractional_flow, p) - formulate!(pump, p, storage) - formulate!(outlet, p, storage) + formulate!(linear_resistance, p, current_level, flow, t) + formulate!(manning_resistance, p, current_level, flow, t) + formulate!(tabulated_rating_curve, p, storage, current_level, flow, t) + formulate!(flow_boundary, p, flow, t) + formulate!(fractional_flow, flow, p) + formulate!(pump, p, flow, storage) + formulate!(outlet, p, flow, storage) return nothing end @@ -950,8 +996,8 @@ end The right hand side function of the system of ODEs set up by Ribasim. """ function water_balance!( - du::ComponentVector{Float64}, - u::ComponentVector{Float64}, + du::ComponentVector, + u::ComponentVector, p::Parameters, t::Float64, )::Nothing @@ -961,22 +1007,38 @@ function water_balance!( integral = u.integral du .= 0.0 - nonzeros(connectivity.flow) .= 0.0 + flow = get_tmp(connectivity.flow, u) + nonzeros(flow) .= 0.0 + + current_area = get_tmp(basin.current_area, u) + current_level = get_tmp(basin.current_level, u) + pid_error = get_tmp(pid_control.error, u) # Ensures current_* vectors are current - set_current_basin_properties!(basin, storage, t) + set_current_basin_properties!(basin, current_area, current_level, storage, t) # Basin forcings - formulate!(du, basin, storage, t) + formulate!(du, basin, storage, current_area, current_level, t) # First formulate intermediate flows - formulate_flows!(p, storage, t) + formulate_flows!(p, storage, current_level, flow, t) # Now formulate du - formulate!(du, connectivity, basin) + formulate!(du, connectivity, flow, basin) # PID control (changes the du of PID controlled basins) - continuous_control!(u, du, pid_control, p, integral, t) + continuous_control!( + u, + du, + current_area, + pid_control, + p, + integral, + current_level, + flow, + pid_error, + t, + ) # Negative storage musn't decrease, based on Shampine's et. al. advice # https://docs.sciml.ai/DiffEqCallbacks/stable/step_control/#DiffEqCallbacks.PositiveDomain diff --git a/core/src/utils.jl b/core/src/utils.jl index cdb7b4360..13da0e4b6 100644 --- a/core/src/utils.jl +++ b/core/src/utils.jl @@ -34,7 +34,7 @@ function create_graph( end "Calculate a profile storage by integrating the areas over the levels" -function profile_storage(levels::Vector{Float64}, areas::Vector{Float64})::Vector{Float64} +function profile_storage(levels::Vector, areas::Vector)::Vector{Float64} # profile starts at the bottom; first storage is 0 storages = zero(areas) n = length(storages) @@ -110,7 +110,7 @@ end """Compute the storages of the basins based on the water level of the basins.""" function get_storages_from_levels( basin::Basin, - levels::Vector{Float64}, + levels::Vector, )::Tuple{Vector{Float64}, Bool} storages = Float64[] @@ -124,11 +124,7 @@ end Compute the area and level of a basin given its storage. Also returns darea/dlevel as it is needed for the Jacobian. """ -function get_area_and_level( - basin::Basin, - state_idx::Int, - storage::Float64, -)::Tuple{Float64, Float64, Float64} +function get_area_and_level(basin::Basin, state_idx::Int, storage::Real)::Tuple{Real, Real} storage_discrete = basin.storage[state_idx] area_discrete = basin.area[state_idx] level_discrete = basin.level[state_idx] @@ -137,11 +133,11 @@ function get_area_and_level( end function get_area_and_level( - storage_discrete::Vector{Float64}, - area_discrete::Vector{Float64}, - level_discrete::Vector{Float64}, - storage::Float64, -)::Tuple{Float64, Float64, Float64} + storage_discrete::Vector, + area_discrete::Vector, + level_discrete::Vector, + storage::Real, +)::Tuple{Real, Real} # storage_idx: smallest index such that storage_discrete[storage_idx] >= storage storage_idx = searchsortedfirst(storage_discrete, storage) @@ -209,7 +205,7 @@ function get_area_and_level( end end - return area, level, darea + return area, level end """ @@ -410,10 +406,10 @@ end Get the current water level of a node ID. The ID can belong to either a Basin or a LevelBoundary. """ -function get_level(p::Parameters, node_id::Int, t::Float64)::Float64 +function get_level(p::Parameters, node_id::Int, current_level, t::Float64)::Real (; basin, level_boundary) = p # since the node_id fields are already Indices, Dictionary creation is instant - basin = Dictionary(basin.node_id, basin.current_level) + basin = Dictionary(basin.node_id, current_level) hasindex, token = gettoken(basin, node_id) return if hasindex gettokenvalue(basin, token) @@ -764,7 +760,7 @@ function update_jac_prototype!( has_index_out, idx_out = id_index(basin.node_id, id_out) if has_index_out - push!(idxs_out, idx_out) + jac_prototype[idx_in, idx_out] = 1.0 end else for idx_out in idxs_out diff --git a/core/src/validation.jl b/core/src/validation.jl index 32e99b100..2eb1919d4 100644 --- a/core/src/validation.jl +++ b/core/src/validation.jl @@ -411,7 +411,7 @@ Test whether static or discrete controlled flow rates are indeed non-negative. """ function valid_flow_rates( node_id::Vector{Int}, - flow_rate::Vector{Float64}, + flow_rate::Vector, control_mapping::Dict{Tuple{Int, String}, NamedTuple}, node_type::Symbol, )::Bool diff --git a/core/test/config.jl b/core/test/config.jl index 6333b64d1..5e61ef2c8 100644 --- a/core/test/config.jl +++ b/core/test/config.jl @@ -70,7 +70,7 @@ end Ribasim.algorithm(Ribasim.Solver(; algorithm = "QNDF", autodiff = false)), ) == AutoFiniteDiff() @test alg_autodiff(Ribasim.algorithm(Ribasim.Solver(; algorithm = "QNDF"))) == - AutoFiniteDiff() + AutoForwardDiff() # autodiff is not a kwargs for explicit algorithms, but we use try-catch to bypass Ribasim.algorithm(Ribasim.Solver(; algorithm = "Euler", autodiff = true)) end diff --git a/core/test/control.jl b/core/test/control.jl index 00da8160d..4480099fe 100644 --- a/core/test/control.jl +++ b/core/test/control.jl @@ -109,7 +109,7 @@ end a = abs(Δlevel / cos(phi)) # This bound is the exact envelope of the analytical solution bound = @. a * exp(alpha * timesteps[1:idx_target_change]) - eps = 3.5e-3 + eps = 5e-3 # Initial convergence to target level @test all(@. abs(level[1:idx_target_change] - target_level) < bound + eps) # Later closeness to target level diff --git a/core/test/docs.toml b/core/test/docs.toml index 576166f00..90e01e399 100644 --- a/core/test/docs.toml +++ b/core/test/docs.toml @@ -32,8 +32,7 @@ abstol = 1e-6 # optional, default 1e-6 reltol = 1e-3 # optional, default 1e-3 maxiters = 1e9 # optional, default 1e9 sparse = true # optional, default true -jac = true # optional, default true -autodiff = false # optional, default false +autodiff = true # optional, default true [logging] # defines the logging level of Ribasim diff --git a/core/test/equations.jl b/core/test/equations.jl index eacf49cce..50c61bff4 100644 --- a/core/test/equations.jl +++ b/core/test/equations.jl @@ -5,6 +5,7 @@ using Arrow import BasicModelInterface as BMI using SciMLBase: successful_retcode using TimerOutputs +using PreallocationTools: get_tmp include("../../utils/testdata.jl") @@ -203,7 +204,8 @@ end (; flow_boundary, fractional_flow, pump) = p q_boundary = flow_boundary.flow_rate[1].u[1] - q_pump = pump.flow_rate[1] + pump_flow_rate = get_tmp(pump.flow_rate, 0) + q_pump = pump_flow_rate[1] frac = fractional_flow.fraction[1] storage_both = Ribasim.get_storages_and_levels(model).storage diff --git a/core/test/run_models.jl b/core/test/run_models.jl index fa40bfda9..69721aea7 100644 --- a/core/test/run_models.jl +++ b/core/test/run_models.jl @@ -4,6 +4,7 @@ using Ribasim import BasicModelInterface as BMI using SciMLBase: successful_retcode import Tables +using PreallocationTools: get_tmp @testset "trivial model" begin toml_path = normpath(@__DIR__, "../../data/trivial/trivial.toml") @@ -51,26 +52,26 @@ end Sys.isapple() end -@testset "sparse and jac solver options" begin +@testset "sparse and AD/FDM jac solver options" begin toml_path = normpath(@__DIR__, "../../data/basic_transient/basic_transient.toml") - config = Ribasim.Config(toml_path; solver_sparse = true, solver_jac = true) - sparse_jac = Ribasim.run(config) - config = Ribasim.Config(toml_path; solver_sparse = false, solver_jac = true) - dense_jac = Ribasim.run(config) - config = Ribasim.Config(toml_path; solver_sparse = true, solver_jac = false) + config = Ribasim.Config(toml_path; solver_sparse = true, solver_autodiff = true) + sparse_ad = Ribasim.run(config) + config = Ribasim.Config(toml_path; solver_sparse = false, solver_autodiff = true) + dense_ad = Ribasim.run(config) + config = Ribasim.Config(toml_path; solver_sparse = true, solver_autodiff = false) sparse_fdm = Ribasim.run(config) - config = Ribasim.Config(toml_path; solver_sparse = false, solver_jac = false) + config = Ribasim.Config(toml_path; solver_sparse = false, solver_autodiff = false) dense_fdm = Ribasim.run(config) - @test successful_retcode(sparse_jac) - @test successful_retcode(dense_jac) + @test successful_retcode(sparse_ad) + @test successful_retcode(dense_ad) @test successful_retcode(sparse_fdm) @test successful_retcode(dense_fdm) - @test dense_jac.integrator.sol.u[end] ≈ sparse_jac.integrator.sol.u[end] - @test sparse_fdm.integrator.sol.u[end] ≈ sparse_jac.integrator.sol.u[end] atol = 1e-3 - @test dense_fdm.integrator.sol.u[end] ≈ sparse_jac.integrator.sol.u[end] atol = 1e-3 + @test dense_ad.integrator.sol.u[end] ≈ sparse_ad.integrator.sol.u[end] atol = 1e-3 + @test sparse_fdm.integrator.sol.u[end] ≈ sparse_ad.integrator.sol.u[end] + @test dense_fdm.integrator.sol.u[end] ≈ sparse_ad.integrator.sol.u[end] atol = 1e-3 end @testset "TabulatedRatingCurve model" begin @@ -194,7 +195,7 @@ end u = model.integrator.sol.u[end] p = model.integrator.p - h_actual = p.basin.current_level + h_actual = get_tmp(p.basin.current_level, u) x = collect(10.0:20.0:990.0) h_expected = standard_step_method(x, 5.0, 1.0, 0.04, h_actual[end], 1.0e-6) diff --git a/core/test/utils.jl b/core/test/utils.jl index 817df863b..a0ae96bd7 100644 --- a/core/test/utils.jl +++ b/core/test/utils.jl @@ -32,7 +32,6 @@ end [2.0, 3.0], [2.0, 3.0], [2.0, 3.0], - [2.0, 3.0], darea, area, level, @@ -89,7 +88,6 @@ end zeros(1), zeros(1), zeros(1), - zeros(1), [area], [level], [storage], @@ -173,9 +171,9 @@ end @test jac_prototype.m == 4 @test jac_prototype.n == 4 - @test jac_prototype.colptr == [1, 3, 5, 7, 9] - @test jac_prototype.rowval == [1, 2, 1, 2, 2, 3, 2, 4] - @test jac_prototype.nzval == ones(8) + @test jac_prototype.colptr == [1, 3, 5, 7, 10] + @test jac_prototype.rowval == [1, 2, 1, 2, 2, 3, 2, 3, 4] + @test jac_prototype.nzval == ones(9) toml_path = normpath(@__DIR__, "../../data/pid_control/pid_control.toml") diff --git a/core/test/validation.jl b/core/test/validation.jl index f3ff75a71..64720ea41 100644 --- a/core/test/validation.jl +++ b/core/test/validation.jl @@ -76,13 +76,13 @@ end @test length(logger.logs) == 3 @test logger.logs[1].level == Error @test logger.logs[1].message == - "Nodes of type Ribasim.Pump can have at most 1 flow inneighbor(s) (got 2 for node #1)." + "Nodes of type Ribasim.Pump{Vector{Float64}} can have at most 1 flow inneighbor(s) (got 2 for node #1)." @test logger.logs[2].level == Error @test logger.logs[2].message == - "Nodes of type Ribasim.Pump must have at least 1 flow outneighbor(s) (got 0 for node #1)." + "Nodes of type Ribasim.Pump{Vector{Float64}} must have at least 1 flow outneighbor(s) (got 0 for node #1)." @test logger.logs[3].level == Error @test logger.logs[3].message == - "Nodes of type Ribasim.Pump must have at least 1 flow inneighbor(s) (got 0 for node #6)." + "Nodes of type Ribasim.Pump{Vector{Float64}} must have at least 1 flow inneighbor(s) (got 0 for node #6)." add_edge!(graph_flow, 2, 5) add_edge!(graph_flow, 5, 3) diff --git a/docs/core/usage.qmd b/docs/core/usage.qmd index 009383f3d..83afbdc67 100644 --- a/docs/core/usage.qmd +++ b/docs/core/usage.qmd @@ -54,10 +54,8 @@ For instance, `saveat = 86400.0` will save output after every day that passed. The Jacobian matrix provides information about the local sensitivity of the model with respect to changes in the states. For implicit solvers it must be calculated often, which can be expensive to do. There are several methods to do this. -By default Ribasim uses an analytical Jacobian function written by the Ribasim developers. -This is enabled by default, with the solver setting `jac = true`. -If this is not used, the Jacobian is calculated with a finite difference method, which can be less accurate and more expensive. -In the future we want to additionally support `autodiff = true` to calculate the Jacobian with automatic differentiation. +By default Ribasim uses a Jacobian derived automatically using [ForwardDiff.jl](https://juliadiff.org/ForwardDiff.jl/stable/) with memory management provided by [PreallocationTools.jl](https://docs.sciml.ai/PreallocationTools/stable/). +If this is not used by setting `autodiff = false`, the Jacobian is calculated with a finite difference method, which can be less accurate and more expensive. By default the Jacobian matrix is a sparse matrix (`sparse = true`). Since each state typically only depends on a small number of other states, this is generally more efficient, especially for larger models. diff --git a/libribasim/Manifest.toml b/libribasim/Manifest.toml new file mode 100644 index 000000000..fd0914633 --- /dev/null +++ b/libribasim/Manifest.toml @@ -0,0 +1,7 @@ +# This file is machine-generated - editing it directly is not advised + +julia_version = "1.9.2" +manifest_format = "2.0" +project_hash = "da39a3ee5e6b4b0d3255bfef95601890afd80709" + +[deps] diff --git a/python/ribasim/ribasim/model.py b/python/ribasim/ribasim/model.py index ec4458219..66c2e3bdc 100644 --- a/python/ribasim/ribasim/model.py +++ b/python/ribasim/ribasim/model.py @@ -45,7 +45,6 @@ class Solver(BaseModel): reltol: Optional[float] maxiters: Optional[int] sparse: Optional[bool] - jac: Optional[bool] autodiff: Optional[bool] diff --git a/python/ribasim_testmodels/ribasim_testmodels/dutch_waterways.py b/python/ribasim_testmodels/ribasim_testmodels/dutch_waterways.py index 3f8a1c988..f0942de6e 100644 --- a/python/ribasim_testmodels/ribasim_testmodels/dutch_waterways.py +++ b/python/ribasim_testmodels/ribasim_testmodels/dutch_waterways.py @@ -275,8 +275,6 @@ def dutch_waterways_model(): ) ) - solver = ribasim.Solver(jac=False) - model = ribasim.Model( modelname="dutch_waterways", node=node, @@ -289,7 +287,6 @@ def dutch_waterways_model(): tabulated_rating_curve=rating_curve, pid_control=pid_control, discrete_control=discrete_control, - solver=solver, starttime="2020-01-01 00:00:00", endtime="2021-01-01 00:00:00", ) diff --git a/ribasim_cli/Manifest.toml b/ribasim_cli/Manifest.toml new file mode 100644 index 000000000..fd0914633 --- /dev/null +++ b/ribasim_cli/Manifest.toml @@ -0,0 +1,7 @@ +# This file is machine-generated - editing it directly is not advised + +julia_version = "1.9.2" +manifest_format = "2.0" +project_hash = "da39a3ee5e6b4b0d3255bfef95601890afd80709" + +[deps] From be0dad796ff8940642153a8700aa82a616bdfcaf Mon Sep 17 00:00:00 2001 From: Martijn Visser Date: Fri, 1 Sep 2023 16:54:47 +0200 Subject: [PATCH 08/13] resolve manifests (#571) ed18bb84d392e6f965d351d81d45d4b47baa3c66 accidentally added new manifests rather than resolving the existing ones. --- build/libribasim/Manifest.toml | 4 ++-- build/ribasim_cli/Manifest.toml | 4 ++-- libribasim/Manifest.toml | 7 ------- ribasim_cli/Manifest.toml | 7 ------- 4 files changed, 4 insertions(+), 18 deletions(-) delete mode 100644 libribasim/Manifest.toml delete mode 100644 ribasim_cli/Manifest.toml diff --git a/build/libribasim/Manifest.toml b/build/libribasim/Manifest.toml index 917829205..a55156dfb 100644 --- a/build/libribasim/Manifest.toml +++ b/build/libribasim/Manifest.toml @@ -855,10 +855,10 @@ uuid = "ae029012-a4dd-5104-9daa-d747884805df" version = "1.3.0" [[deps.Ribasim]] -deps = ["Arrow", "BasicModelInterface", "CodecLz4", "CodecZstd", "ComponentArrays", "Configurations", "DBInterface", "DataInterpolations", "DataStructures", "Dates", "Dictionaries", "DiffEqCallbacks", "FiniteDiff", "Graphs", "IterTools", "Legolas", "Logging", "LoggingExtras", "OrdinaryDiffEq", "SQLite", "SciMLBase", "SparseArrays", "StructArrays", "Tables", "TerminalLoggers", "TimerOutputs", "TranscodingStreams"] +deps = ["Arrow", "BasicModelInterface", "CodecLz4", "CodecZstd", "ComponentArrays", "Configurations", "DBInterface", "DataInterpolations", "DataStructures", "Dates", "Dictionaries", "DiffEqCallbacks", "FiniteDiff", "ForwardDiff", "Graphs", "IterTools", "Legolas", "Logging", "LoggingExtras", "OrdinaryDiffEq", "PreallocationTools", "SQLite", "SciMLBase", "SparseArrays", "StructArrays", "Tables", "TerminalLoggers", "TimerOutputs", "TranscodingStreams"] path = "../../core" uuid = "aac5e3d9-0b8f-4d4f-8241-b1a7a9632635" -version = "0.1.0" +version = "0.2.0" [[deps.RuntimeGeneratedFunctions]] deps = ["ExprTools", "SHA", "Serialization"] diff --git a/build/ribasim_cli/Manifest.toml b/build/ribasim_cli/Manifest.toml index 039ba14d1..ea21ec27a 100644 --- a/build/ribasim_cli/Manifest.toml +++ b/build/ribasim_cli/Manifest.toml @@ -855,10 +855,10 @@ uuid = "ae029012-a4dd-5104-9daa-d747884805df" version = "1.3.0" [[deps.Ribasim]] -deps = ["Arrow", "BasicModelInterface", "CodecLz4", "CodecZstd", "ComponentArrays", "Configurations", "DBInterface", "DataInterpolations", "DataStructures", "Dates", "Dictionaries", "DiffEqCallbacks", "FiniteDiff", "Graphs", "IterTools", "Legolas", "Logging", "LoggingExtras", "OrdinaryDiffEq", "SQLite", "SciMLBase", "SparseArrays", "StructArrays", "Tables", "TerminalLoggers", "TimerOutputs", "TranscodingStreams"] +deps = ["Arrow", "BasicModelInterface", "CodecLz4", "CodecZstd", "ComponentArrays", "Configurations", "DBInterface", "DataInterpolations", "DataStructures", "Dates", "Dictionaries", "DiffEqCallbacks", "FiniteDiff", "ForwardDiff", "Graphs", "IterTools", "Legolas", "Logging", "LoggingExtras", "OrdinaryDiffEq", "PreallocationTools", "SQLite", "SciMLBase", "SparseArrays", "StructArrays", "Tables", "TerminalLoggers", "TimerOutputs", "TranscodingStreams"] path = "../../core" uuid = "aac5e3d9-0b8f-4d4f-8241-b1a7a9632635" -version = "0.1.0" +version = "0.2.0" [[deps.RuntimeGeneratedFunctions]] deps = ["ExprTools", "SHA", "Serialization"] diff --git a/libribasim/Manifest.toml b/libribasim/Manifest.toml deleted file mode 100644 index fd0914633..000000000 --- a/libribasim/Manifest.toml +++ /dev/null @@ -1,7 +0,0 @@ -# This file is machine-generated - editing it directly is not advised - -julia_version = "1.9.2" -manifest_format = "2.0" -project_hash = "da39a3ee5e6b4b0d3255bfef95601890afd80709" - -[deps] diff --git a/ribasim_cli/Manifest.toml b/ribasim_cli/Manifest.toml deleted file mode 100644 index fd0914633..000000000 --- a/ribasim_cli/Manifest.toml +++ /dev/null @@ -1,7 +0,0 @@ -# This file is machine-generated - editing it directly is not advised - -julia_version = "1.9.2" -manifest_format = "2.0" -project_hash = "da39a3ee5e6b4b0d3255bfef95601890afd80709" - -[deps] From 12abc689a71c7f2d228a265ed6e56c21e01f4ad2 Mon Sep 17 00:00:00 2001 From: Martijn Visser Date: Tue, 5 Sep 2023 19:19:20 +0200 Subject: [PATCH 09/13] avoid returning DiffCache through BMI (#575) This should fix the TeamCity libribasim tests. `test_bmi.py` did not catch it since only getting volume is tested: `libribasim.get_value_ptr("volume")`. I'll add those tests later, but first let's unbreak CI. --- core/src/bmi.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/core/src/bmi.jl b/core/src/bmi.jl index c8c7db092..2890ead2c 100644 --- a/core/src/bmi.jl +++ b/core/src/bmi.jl @@ -508,7 +508,7 @@ function BMI.get_value_ptr(model::Model, name::AbstractString) if name == "volume" model.integrator.u.storage elseif name == "level" - model.integrator.p.basin.current_level + get_tmp(model.integrator.p.basin.current_level, 0) elseif name == "infiltration" model.integrator.p.basin.infiltration elseif name == "drainage" From 9a8b2a7f6b7758476a9af7b3dc387966698f5b81 Mon Sep 17 00:00:00 2001 From: visser_mn Date: Thu, 7 Sep 2023 15:35:43 +0200 Subject: [PATCH 10/13] TeamCity change in 'Ribasim / Ribasim' project: artifact dependencies of 'Test libribasim - Windows' build configuration were updated --- .../Ribasim_Ribasim_TestLibribasimWindows.xml | 19 ++--- .../Ribasim_Ribasim_TestRibasimCliWindows.xml | 82 +++++++++++++++++++ 2 files changed, 89 insertions(+), 12 deletions(-) create mode 100644 .teamcity/Ribasim_Ribasim/buildTypes/Ribasim_Ribasim_TestRibasimCliWindows.xml diff --git a/.teamcity/Ribasim_Ribasim/buildTypes/Ribasim_Ribasim_TestLibribasimWindows.xml b/.teamcity/Ribasim_Ribasim/buildTypes/Ribasim_Ribasim_TestLibribasimWindows.xml index 087d5109c..b628adbb6 100644 --- a/.teamcity/Ribasim_Ribasim/buildTypes/Ribasim_Ribasim_TestLibribasimWindows.xml +++ b/.teamcity/Ribasim_Ribasim/buildTypes/Ribasim_Ribasim_TestLibribasimWindows.xml @@ -6,9 +6,6 @@ - - - @@ -16,8 +13,7 @@ - +call conda env create --file environment.yml -p "%conda_env_path%"]]> @@ -25,9 +21,9 @@ call conda create --prefix "%conda_env_path%" python=3.10]]> +pip install "python/ribasim" +pip install "python/ribasim_testmodels" +pip install "python/ribasim_api"]]> @@ -35,7 +31,7 @@ call pip install --editable "python/ribasim_api[tests]"]]> +pytest tests --basetemp=tests/temp --junitxml="report.xml"]]> @@ -46,8 +42,7 @@ call pytest tests --basetemp=tests/temp --junitxml="report.xml"]]> - - + @@ -78,7 +73,7 @@ call pytest tests --basetemp=tests/temp --junitxml="report.xml"]]> - + diff --git a/.teamcity/Ribasim_Ribasim/buildTypes/Ribasim_Ribasim_TestRibasimCliWindows.xml b/.teamcity/Ribasim_Ribasim/buildTypes/Ribasim_Ribasim_TestRibasimCliWindows.xml new file mode 100644 index 000000000..ec3fdb75a --- /dev/null +++ b/.teamcity/Ribasim_Ribasim/buildTypes/Ribasim_Ribasim_TestRibasimCliWindows.xml @@ -0,0 +1,82 @@ + + + Test ribasim_cli - Windows + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + From 8bbab875cfdde822d3ee0f73be2020eced1ffecc Mon Sep 17 00:00:00 2001 From: Martijn Visser Date: Fri, 8 Sep 2023 20:26:10 +0200 Subject: [PATCH 11/13] Use FixedSizeDiffCache for flows (#581) When profiling runs with the default `autodiff=true`, this line was responsible for 35% of the time and almost all allocations: https://github.com/Deltares/Ribasim/blob/b3eb044a722d1655c5465bafe50951b75fe960d6/core/src/solve.jl#L1002 With this PR that drops down to 0%. `connectivity.flow` is a sparse matrix, but the DiffCache does not seem to like sparse matrixes. The `dual_du` field was a dense vector of length n x n x cache_size, and the `get_tmp` call led to further allocations trying to restructure the sparse matrix from the vector. Luckily there is the FixedSizeDiffCache that helps here: https://docs.sciml.ai/PreallocationTools/stable/#FixedSizeDiffCache This retains the sparsity in the dual, and returns a `ReinterpretArray` from `get_tmp` during autodiff. To avoid materializing this reinterpretarray I needed to additionally fill the parent array with zeros rather than the array itself. There is another unrelated performance fix here, and that is to concretely type the Parameter struct, by adding type parameters from its fields. Otherwise you have situations like ```julia struct A a::Vector end ``` where the compiler doesn't know the element type of the Vector, so it can perform less optimizations. The solution: ```julia struct A{T} a::Vector{T} end ``` Finally I consistently added AbstractVector/Matrix argument type annotations to ensure the ReinterpretArray could enter everywhere. And I renamed the functions to formulate flows to `formulate_flow`, to make it easier to separate them from the other `formulate!` methods. --- core/src/Ribasim.jl | 2 +- core/src/create.jl | 5 ++- core/src/solve.jl | 83 ++++++++++++++++++++----------------- core/test/run_models.jl | 5 +++ docs/contribute/addnode.qmd | 4 +- 5 files changed, 57 insertions(+), 42 deletions(-) diff --git a/core/src/Ribasim.jl b/core/src/Ribasim.jl index c1bc06678..78b615aaa 100644 --- a/core/src/Ribasim.jl +++ b/core/src/Ribasim.jl @@ -20,7 +20,7 @@ using Legolas: Legolas, @schema, @version, validate, SchemaVersion, declared using Logging: current_logger, min_enabled_level, with_logger using LoggingExtras: EarlyFilteredLogger, LevelOverrideLogger using OrdinaryDiffEq -using PreallocationTools: DiffCache, get_tmp +using PreallocationTools: DiffCache, FixedSizeDiffCache, get_tmp using SciMLBase using SparseArrays using SQLite: SQLite, DB, Query, esc_id diff --git a/core/src/create.jl b/core/src/create.jl index 8221f4504..c75a6c0e8 100644 --- a/core/src/create.jl +++ b/core/src/create.jl @@ -202,10 +202,11 @@ function Connectivity(db::DB, config::Config, chunk_size::Int)::Connectivity edge_ids_flow_inv = Dictionary(values(edge_ids_flow), keys(edge_ids_flow)) flow = adjacency_matrix(graph_flow, Float64) - nonzeros(flow) .= 0.0 + flow .= 0.0 if config.solver.autodiff - flow = DiffCache(flow, chunk_size) + # FixedSizeDiffCache performs better for sparse matrix + flow = FixedSizeDiffCache(flow, chunk_size) end return Connectivity( diff --git a/core/src/solve.jl b/core/src/solve.jl index cca2d16e1..a29c7d82f 100644 --- a/core/src/solve.jl +++ b/core/src/solve.jl @@ -401,21 +401,21 @@ struct PidControl{T} <: AbstractParameterNode end # TODO Automatically add all nodetypes here -struct Parameters +struct Parameters{T, TSparse, C1, C2} starttime::DateTime - connectivity::Connectivity - basin::Basin + connectivity::Connectivity{TSparse} + basin::Basin{T, C1} linear_resistance::LinearResistance manning_resistance::ManningResistance - tabulated_rating_curve::TabulatedRatingCurve + tabulated_rating_curve::TabulatedRatingCurve{C2} fractional_flow::FractionalFlow level_boundary::LevelBoundary flow_boundary::FlowBoundary - pump::Pump - outlet::Outlet + pump::Pump{T} + outlet::Outlet{T} terminal::Terminal discrete_control::DiscreteControl - pid_control::PidControl + pid_control::PidControl{T} lookup::Dict{Int, Symbol} end @@ -508,8 +508,8 @@ function formulate!( du::AbstractVector, basin::Basin, storage::AbstractVector, - current_area, - current_level, + current_area::AbstractVector, + current_level::AbstractVector, t::Float64, )::Nothing for i in eachindex(storage) @@ -553,13 +553,13 @@ end function continuous_control!( u::ComponentVector, du::ComponentVector, - current_area, + current_area::AbstractVector, pid_control::PidControl, p::Parameters, integral_value::SubArray, current_level::AbstractVector, - flow::SparseMatrixCSC, - pid_error, + flow::AbstractMatrix, + pid_error::AbstractVector, t::Float64, )::Nothing (; connectivity, pump, outlet, basin, fractional_flow) = p @@ -700,11 +700,11 @@ end """ Directed graph: outflow is positive! """ -function formulate!( +function formulate_flow!( linear_resistance::LinearResistance, p::Parameters, - current_level, - flow, + current_level::AbstractVector, + flow::AbstractMatrix, t::Float64, )::Nothing (; connectivity) = p @@ -733,17 +733,16 @@ end """ Directed graph: outflow is positive! """ -function formulate!( +function formulate_flow!( tabulated_rating_curve::TabulatedRatingCurve, p::Parameters, - storage::AbstractVector{Float64}, - current_level, + storage::AbstractVector, + current_level::AbstractVector, flow::SparseMatrixCSC, t::Float64, )::Nothing (; basin, connectivity) = p - (; graph_flow, flow) = connectivity - (; graph_flow, flow) = connectivity + (; graph_flow) = connectivity (; node_id, active, tables) = tabulated_rating_curve for (i, id) in enumerate(node_id) upstream_basin_id = only(inneighbors(graph_flow, id)) @@ -808,11 +807,11 @@ The average of the upstream and downstream water depth is used to compute cross- hydraulic radius. This ensures that a basin can receive water after it has gone dry. """ -function formulate!( +function formulate_flow!( manning_resistance::ManningResistance, p::Parameters, - current_level, - flow, + current_level::AbstractVector, + flow::AbstractMatrix, t::Float64, )::Nothing (; basin, connectivity) = p @@ -867,7 +866,11 @@ function formulate!( return nothing end -function formulate!(fractional_flow::FractionalFlow, flow, p::Parameters)::Nothing +function formulate_flow!( + fractional_flow::FractionalFlow, + flow::AbstractMatrix, + p::Parameters, +)::Nothing (; connectivity) = p (; graph_flow) = connectivity (; node_id, fraction) = fractional_flow @@ -879,7 +882,12 @@ function formulate!(fractional_flow::FractionalFlow, flow, p::Parameters)::Nothi return nothing end -function formulate!(flow_boundary::FlowBoundary, p::Parameters, flow, t::Float64)::Nothing +function formulate_flow!( + flow_boundary::FlowBoundary, + p::Parameters, + flow::AbstractMatrix, + t::Float64, +)::Nothing (; connectivity) = p (; graph_flow) = connectivity (; node_id, active, flow_rate) = flow_boundary @@ -900,10 +908,10 @@ function formulate!(flow_boundary::FlowBoundary, p::Parameters, flow, t::Float64 end end -function formulate!( +function formulate_flow!( node::Union{Pump, Outlet}, p::Parameters, - flow, + flow::AbstractMatrix, storage::AbstractVector, )::Nothing (; connectivity, basin) = p @@ -946,7 +954,7 @@ end function formulate!( du::ComponentVector, connectivity::Connectivity, - flow::SparseMatrixCSC, + flow::AbstractMatrix, basin::Basin, )::Nothing # loop over basins @@ -968,7 +976,7 @@ function formulate_flows!( p::Parameters, storage::AbstractVector, current_level, - flow::SparseMatrixCSC, + flow::AbstractMatrix, t::Float64, )::Nothing (; @@ -981,13 +989,13 @@ function formulate_flows!( outlet, ) = p - formulate!(linear_resistance, p, current_level, flow, t) - formulate!(manning_resistance, p, current_level, flow, t) - formulate!(tabulated_rating_curve, p, storage, current_level, flow, t) - formulate!(flow_boundary, p, flow, t) - formulate!(fractional_flow, flow, p) - formulate!(pump, p, flow, storage) - formulate!(outlet, p, flow, storage) + formulate_flow!(linear_resistance, p, current_level, flow, t) + formulate_flow!(manning_resistance, p, current_level, flow, t) + formulate_flow!(tabulated_rating_curve, p, storage, current_level, flow, t) + formulate_flow!(flow_boundary, p, flow, t) + formulate_flow!(fractional_flow, flow, p) + formulate_flow!(pump, p, flow, storage) + formulate_flow!(outlet, p, flow, storage) return nothing end @@ -1008,7 +1016,8 @@ function water_balance!( du .= 0.0 flow = get_tmp(connectivity.flow, u) - nonzeros(flow) .= 0.0 + # use parent to avoid materializing the ReinterpretArray from FixedSizeDiffCache + parent(flow) .= 0.0 current_area = get_tmp(basin.current_area, u) current_level = get_tmp(basin.current_level, u) diff --git a/core/test/run_models.jl b/core/test/run_models.jl index 69721aea7..581487201 100644 --- a/core/test/run_models.jl +++ b/core/test/run_models.jl @@ -32,6 +32,11 @@ end end @test model isa Ribasim.Model + p = model.integrator.p + @test p isa Ribasim.Parameters + @test isconcretetype(typeof(p)) + @test all(isconcretetype, fieldtypes(typeof(p))) + @test successful_retcode(model) @test model.integrator.sol.u[end] ≈ Float32[519.8817, 519.8798, 339.3959, 1418.4331] skip = Sys.isapple() atol = 1.5 diff --git a/docs/contribute/addnode.qmd b/docs/contribute/addnode.qmd index b925824ce..3c15e841f 100644 --- a/docs/contribute/addnode.qmd +++ b/docs/contribute/addnode.qmd @@ -61,10 +61,10 @@ end ## Node behavior -In general if the new node type dictates flow, the behaviour of the new node in the Ribasim core is defined in a method of the `formulate!` function, which is called within the `water_balance!` (both in `solve.jl`) function being the right hand side of the system of differential equations solved by Ribasim. Here the details depend highly on the specifics of the node type. An example structure of a `formulate!` method is given below. +In general if the new node type dictates flow, the behaviour of the new node in the Ribasim core is defined in a method of the `formulate_flow!` function, which is called within the `water_balance!` (both in `solve.jl`) function being the right hand side of the system of differential equations solved by Ribasim. Here the details depend highly on the specifics of the node type. An example structure of a `formulate_flow!` method is given below. ```julia -function formulate!(new_node_type::NewNodeType, p::Parameters)::Nothing +function formulate_flow!(new_node_type::NewNodeType, p::Parameters)::Nothing # Retrieve relevant parameters (; connectivity) = p (; flow) = connectivity From 196b50d61e5bd4ee114b6f668bf31f56e5dd5430 Mon Sep 17 00:00:00 2001 From: Martijn Visser Date: Fri, 8 Sep 2023 20:54:43 +0200 Subject: [PATCH 12/13] avoid allocating value vectors in get_level (#582) `Dictionary` uses `Indices{I}` as keys, and `Vector{T}` as values. The Parameters contain both, and therefore it was free to construct a `Dictionary` in a frequently called function like `get_level`. However with autodiff, the values could be a ReinterpretArray with Duals instead of just a Vector. This meant that on Dictionary creation it would convert the ReinterpretArray to a Vector, leading to many allocations. This is on top of https://github.com/Deltares/Ribasim/pull/581. After that, this was responsible for 94% of the time spent. With this PR that goes down to about 2%, leading to a nice little speedup. --------- Co-authored-by: Hofer-Julian <30049909+Hofer-Julian@users.noreply.github.com> --- core/src/Ribasim.jl | 2 +- core/src/solve.jl | 2 +- core/src/utils.jl | 32 +++++++++++++++++--------------- 3 files changed, 19 insertions(+), 17 deletions(-) diff --git a/core/src/Ribasim.jl b/core/src/Ribasim.jl index 78b615aaa..8821ef2c8 100644 --- a/core/src/Ribasim.jl +++ b/core/src/Ribasim.jl @@ -12,7 +12,7 @@ using ComponentArrays: ComponentVector using DataInterpolations: LinearInterpolation, derivative using Dates using DBInterface: execute, prepare -using Dictionaries: Indices, Dictionary, gettoken, gettokenvalue, dictionary +using Dictionaries: Indices, Dictionary, gettoken, dictionary using ForwardDiff: pickchunksize using DiffEqCallbacks using Graphs: DiGraph, add_edge!, adjacency_matrix, inneighbors, outneighbors diff --git a/core/src/solve.jl b/core/src/solve.jl index a29c7d82f..8992b9a9a 100644 --- a/core/src/solve.jl +++ b/core/src/solve.jl @@ -975,7 +975,7 @@ end function formulate_flows!( p::Parameters, storage::AbstractVector, - current_level, + current_level::AbstractVector, flow::AbstractMatrix, t::Float64, )::Nothing diff --git a/core/src/utils.jl b/core/src/utils.jl index 13da0e4b6..d4b80baaa 100644 --- a/core/src/utils.jl +++ b/core/src/utils.jl @@ -406,34 +406,36 @@ end Get the current water level of a node ID. The ID can belong to either a Basin or a LevelBoundary. """ -function get_level(p::Parameters, node_id::Int, current_level, t::Float64)::Real +function get_level( + p::Parameters, + node_id::Int, + current_level::AbstractVector, + t::Float64, +)::Real (; basin, level_boundary) = p - # since the node_id fields are already Indices, Dictionary creation is instant - basin = Dictionary(basin.node_id, current_level) - hasindex, token = gettoken(basin, node_id) + hasindex, i = id_index(basin.node_id, node_id) return if hasindex - gettokenvalue(basin, token) + current_level[i] else - boundary = Dictionary(level_boundary.node_id, level_boundary.level) - boundary[node_id](t) + i = findsorted(level_boundary.node_id, node_id)::Int + level_boundary.level[i](t) end end "Get the index of an ID in a set of indices." -function id_index(ids::Indices{Int}, id::Int) - # There might be a better approach for this, this feels too internal - # the second return is the token, a Tuple{Int, Int} - hasindex, (_, idx) = gettoken(ids, id) - return hasindex, idx +function id_index(ids::Indices{Int}, id::Int)::Tuple{Bool, Int} + # We avoid creating Dictionary here since it converts the values to a Vector, + # leading to allocations when used with PreallocationTools's ReinterpretArrays. + hasindex, (_, i) = gettoken(ids, id) + return hasindex, i end "Return the bottom elevation of the basin with index i, or nothing if it doesn't exist" function basin_bottom(basin::Basin, node_id::Int)::Union{Float64, Nothing} - basin = Dictionary(basin.node_id, basin.level) - hasindex, token = gettoken(basin, node_id) + hasindex, i = id_index(basin.node_id, node_id) return if hasindex # get level(storage) interpolation function - level_discrete = gettokenvalue(basin, token) + level_discrete = basin.level[i] # and return the first level in this vector, representing the bottom first(level_discrete) else From 2aecef634bf642d68a4986e6d2c815070d26c9c1 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Mon, 11 Sep 2023 08:57:49 +0200 Subject: [PATCH 13/13] Bump actions/checkout from 3 to 4 (#584) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Bumps [actions/checkout](https://github.com/actions/checkout) from 3 to 4.
Release notes

Sourced from actions/checkout's releases.

v4.0.0

What's Changed

New Contributors

Full Changelog: https://github.com/actions/checkout/compare/v3...v4.0.0

v3.6.0

What's Changed

New Contributors

Full Changelog: https://github.com/actions/checkout/compare/v3.5.3...v3.6.0

v3.5.3

What's Changed

New Contributors

Full Changelog: https://github.com/actions/checkout/compare/v3...v3.5.3

v3.5.2

What's Changed

Full Changelog: https://github.com/actions/checkout/compare/v3.5.1...v3.5.2

v3.5.1

What's Changed

New Contributors

... (truncated)

Changelog

Sourced from actions/checkout's changelog.

Changelog

v4.0.0

v3.6.0

v3.5.3

v3.5.2

v3.5.1

v3.5.0

v3.4.0

v3.3.0

v3.2.0

v3.1.0

v3.0.2

v3.0.1

... (truncated)

Commits

[![Dependabot compatibility score](https://dependabot-badges.githubapp.com/badges/compatibility_score?dependency-name=actions/checkout&package-manager=github_actions&previous-version=3&new-version=4)](https://docs.github.com/en/github/managing-security-vulnerabilities/about-dependabot-security-updates#about-compatibility-scores) Dependabot will resolve any conflicts with this PR as long as you don't alter it yourself. You can also trigger a rebase manually by commenting `@dependabot rebase`. [//]: # (dependabot-automerge-start) [//]: # (dependabot-automerge-end) ---
Dependabot commands and options
You can trigger Dependabot actions by commenting on this PR: - `@dependabot rebase` will rebase this PR - `@dependabot recreate` will recreate this PR, overwriting any edits that have been made to it - `@dependabot merge` will merge this PR after your CI passes on it - `@dependabot squash and merge` will squash and merge this PR after your CI passes on it - `@dependabot cancel merge` will cancel a previously requested merge and block automerging - `@dependabot reopen` will reopen this PR if it is closed - `@dependabot close` will close this PR and stop Dependabot recreating it. You can achieve the same result by closing it manually - `@dependabot show ignore conditions` will show all of the ignore conditions of the specified dependency - `@dependabot ignore this major version` will close this PR and stop Dependabot creating any more for this major version (unless you reopen the PR or upgrade to it yourself) - `@dependabot ignore this minor version` will close this PR and stop Dependabot creating any more for this minor version (unless you reopen the PR or upgrade to it yourself) - `@dependabot ignore this dependency` will close this PR and stop Dependabot creating any more for this dependency (unless you reopen the PR or upgrade to it yourself)
Signed-off-by: dependabot[bot] Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com> --- .github/workflows/core_tests.yml | 2 +- .github/workflows/docs.yml | 2 +- .github/workflows/pre-commit_auto_update.yml | 2 +- .github/workflows/pre-commit_check.yml | 2 +- .github/workflows/python_lint.yml | 2 +- .github/workflows/python_tests.yml | 2 +- .github/workflows/qgis.yml | 2 +- 7 files changed, 7 insertions(+), 7 deletions(-) diff --git a/.github/workflows/core_tests.yml b/.github/workflows/core_tests.yml index 8b3367920..04822cfae 100644 --- a/.github/workflows/core_tests.yml +++ b/.github/workflows/core_tests.yml @@ -28,7 +28,7 @@ jobs: arch: - x64 steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - name: Setup Micromamba uses: mamba-org/setup-micromamba@v1 with: diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml index e4dc493ed..26a58be4e 100644 --- a/.github/workflows/docs.yml +++ b/.github/workflows/docs.yml @@ -24,7 +24,7 @@ jobs: arch: - x64 steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - uses: julia-actions/setup-julia@v1 with: diff --git a/.github/workflows/pre-commit_auto_update.yml b/.github/workflows/pre-commit_auto_update.yml index bf4d418ae..bfe8f0621 100644 --- a/.github/workflows/pre-commit_auto_update.yml +++ b/.github/workflows/pre-commit_auto_update.yml @@ -11,7 +11,7 @@ jobs: auto-update: runs-on: ubuntu-latest steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 with: ssh-key: ${{ secrets.SSH_PRIVATE_KEY }} - uses: actions/setup-python@v4 diff --git a/.github/workflows/pre-commit_check.yml b/.github/workflows/pre-commit_check.yml index fbad75fd6..2876b910a 100644 --- a/.github/workflows/pre-commit_check.yml +++ b/.github/workflows/pre-commit_check.yml @@ -10,7 +10,7 @@ jobs: check: runs-on: ubuntu-latest steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - uses: actions/setup-python@v4 with: python-version: "3.11" diff --git a/.github/workflows/python_lint.yml b/.github/workflows/python_lint.yml index de6c69f28..859391309 100644 --- a/.github/workflows/python_lint.yml +++ b/.github/workflows/python_lint.yml @@ -18,7 +18,7 @@ jobs: run: shell: bash -l {0} steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - name: Setup Micromamba uses: mamba-org/setup-micromamba@v1 with: diff --git a/.github/workflows/python_tests.yml b/.github/workflows/python_tests.yml index fcf173ba8..1400ea31b 100644 --- a/.github/workflows/python_tests.yml +++ b/.github/workflows/python_tests.yml @@ -30,7 +30,7 @@ jobs: arch: - x86 steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - name: Setup Micromamba uses: mamba-org/setup-micromamba@v1 diff --git a/.github/workflows/qgis.yml b/.github/workflows/qgis.yml index 6863ecd4e..788a7e059 100644 --- a/.github/workflows/qgis.yml +++ b/.github/workflows/qgis.yml @@ -17,7 +17,7 @@ jobs: steps: - name: Check out repository - uses: actions/checkout@v3 + uses: actions/checkout@v4 - name: Launching docker compose run: ./start.sh - name: Running tests