From 83939dd2014af338471c4678e16001fdf9a983ab Mon Sep 17 00:00:00 2001 From: Bart de Koning <74617371+SouthEndMusic@users.noreply.github.com> Date: Wed, 31 Jan 2024 11:02:37 +0100 Subject: [PATCH 01/40] Update allocation docs with main network (#1011) Fixes https://github.com/Deltares/Ribasim/issues/1005. --- docs/Project.toml | 2 + docs/core/allocation.qmd | 89 ++++++++++++++++++++++++++++++---------- pixi.toml | 4 +- 3 files changed, 72 insertions(+), 23 deletions(-) diff --git a/docs/Project.toml b/docs/Project.toml index a15dc604d..8174017e2 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -19,6 +19,7 @@ Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" MarkdownTables = "1862ce21-31c7-451e-824c-f20fa3f90fa2" OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" Ribasim = "aac5e3d9-0b8f-4d4f-8241-b1a7a9632635" +SQLite = "0aa819cd-b072-5ff4-a722-6bc24af294d9" [compat] Configurations = "0.17" @@ -33,4 +34,5 @@ Legolas = "0.5" Logging = "<0.0.1,1" MarkdownTables = "1" OrderedCollections = "1.6" +SQLite = "1.5.1" julia = "1.10" diff --git a/docs/core/allocation.qmd b/docs/core/allocation.qmd index f80c79a9c..14b16167e 100644 --- a/docs/core/allocation.qmd +++ b/docs/core/allocation.qmd @@ -2,15 +2,25 @@ title: "Allocation" --- -Allocation is the process of assigning an allocated abstraction flow rate to user nodes in the model based on information about sources, user demands over various priorities, constraints introduced by nodes, local water availability and graph topology. The allocation procedure implemented in Ribasim is heavily inspired by the [maximum flow problem](https://en.wikipedia.org/wiki/Maximum_flow_problem). +# Introduction +Allocation is the process of assigning an allocated abstraction flow rate to user nodes in the physical layer of the model based on information about sources, user demands over various priorities, constraints introduced by nodes, local water availability and graph topology. The allocation procedure implemented in Ribasim is heavily inspired by the [maximum flow problem](https://en.wikipedia.org/wiki/Maximum_flow_problem). -The allocation problem is solved per subnetwork of the Ribasim model. The subnetwork is used to formulate an optimization problem with the [JuMP](https://jump.dev/JuMP.jl/stable/) package, which is solved using the [HiGHS solver](https://highs.dev/). For more information see also the example of solving the maximum flow problem with `JuMP.jl` [here](https://jump.dev/JuMP.jl/stable/tutorials/linear/network_flows/#The-max-flow-problem). +The allocation problem is solved per subnetwork of the Ribasim model. The subnetwork is used to formulate an optimization problem with the [JuMP](https://jump.dev/JuMP.jl/stable/) package, which is solved using the [HiGHS solver](https://highs.dev/). For more in-depth information see also the example of solving the maximum flow problem with `JuMP.jl` [here](https://jump.dev/JuMP.jl/stable/tutorials/linear/network_flows/#The-max-flow-problem). -# The allocation problem +# The high level algorithm +The allocation algorithm consists of 3 types of optimization problems: + +1. **Subnetwork demand collection**: Collect the demands of a subnetwork from the main network by optimizing with unlimited capacity from the main network; +2. **Main network allocation**: Allocate to subnetworks with the above collected demand, and users in the main network; +3. **Subnetwork allocation**: Allocate to users in the subnetworks with the flows allocated to the subnetwork in the main network allocation. + +The total allocation algorithm consists of performing 1 for all subnetworks, then performing 2, then performing 3 for all subnetworks. Not having a main network is also supported, then 1 and 2 are skipped. + +# Elements of allocation The following data of the parameters and state of a Ribasim model are relevant for the allocation problem. -## Allocation problem input +## Schematisation input ### The subnetwork @@ -22,7 +32,7 @@ Sources are indicated by a set of edges in the subnetwork $$ E_S^\text{source} \subset \left(S \times S\right) \cap E. $$ -That is, if $(i,j) \in E_S^\text{source}$, then $Q_{ij}$ (see the [formal model description](equations.qmd#formal-model-description)) is treated as a source flow in the allocation problem. +That is, if $(i,j) \in E_S^\text{source}$, then $Q_{ij}$ (see the [formal model description](equations.qmd#formal-model-description)) is treated as a source flow in the allocation problem. These edges are either coming from a boundary/source node (e.g. a level or flow boundary) or connect the main network to a subnetwork. ### User demands @@ -35,6 +45,8 @@ $$ On this page we assume that the priorities are given by all integers from $1$ to some $p_{\max} \in \mathbb{N}$. However, in the Ribasim input this is not a requirement; some of these in between priority values can be missing, only the ordering of the given priorities is taken into account. ::: +## Simulation (physical layer) input + ### Vertical fluxes and local storage Apart from the source flows denoted by edges, there are other sources of water in the subnetwork, associated with the basins in the subnetwork $B_S = B \cap S$. Firstly there is the sum of the vertical fluxes (precipitation, evaporation, infiltration and drainage) for each basin @@ -48,20 +60,20 @@ $$ $$ Note that this value can be negative, which we interpret as a demand from the basin. -### Flow magnitude and direction constraints +### Constraining factors +#### Flow magnitude and direction constraints Nodes in the Ribasim model that have a `max_flow_rate`, i.e. pumps and outlets, put a constraint on the flow through that node. Some nodes only allow flow in one direction, like pumps, outlets and tabulated rating curves. -### Fractional flows and user return flows - +#### Fractional flows and user return flows Both fractional flow nodes and user nodes dictate proportional relationships between flows over edges in the subnetwork. Users have a return factor $0 \le r_i \le 1, i \in U_S$. -## The allocation optimization problem - -### The allocation subgraph +## The allocation graph A new graph is created from the subnetwork, which we call an allocation graph. The allocation graph is almost a subgraph of the main (flow) model, apart from the fact that an allocation graph can contain edges which are not in the main model. +### Nodes and edges + The allocation graph consists of: - Nodes $V'_S \subset V_S$, where each basin, source and user in the subnetwork get a node in the allocation graph. Also nodes that have fractional flow outneighbors get a node in the allocation graph. @@ -74,11 +86,11 @@ For notational convenience, we use the notation V^{\text{in}}_S(j) = \left\{i \in V'_S : (i,j) \in E_S\right\} \end{align} -for the set of in-neighbors and out-neighbors of a node in the allocation graph respectively +for the set of in-neighbors and out-neighbors of a node in the allocation graph respectively. -### The allocation graph capacities +### Capacities -The capacities of the edges of the allocation graph are collected in the sparse capacity matrix $C_S \in \overline{\mathbb{R}}_{\ge 0}^{n'\times n'}$ where $n' = \#V'_S$ is the number of nodes in the allocation graph. The capacities can be infinite. +Each edge in the allocation graph has an associated capacity. These capacities are collected in the sparse capacity matrix $C_S \in \overline{\mathbb{R}}_{\ge 0}^{n'\times n'}$ where $n' = \#V'_S$ is the number of nodes in the allocation graph. The capacities can be infinite, if there is nothing in the model constraining the capacity of the edge. The capacities are determined in 4 different ways: @@ -86,9 +98,13 @@ The capacities are determined in 4 different ways: - The capacity of the edge $e \in E_S$ is given by the smallest `max_flow_rate` of the nodes along the equivalent edges in the subnetwork. If there are no nodes with a `max_flow_rate`, the edge capacity is infinite; - If the edge is a source, the capacity of the edge is given by the flow rate of that source. -### The optimization variables +# The optimization problem -There are 2 types of variables whose value has to be determined to solve the allocation problem: +The optimization problem for a subnetwork is a linear optimization problem consisting of an objective function with associated constraints on a set of variables, all of which are introduced below. + +## The optimization variables + +There are 2 types of variable whose value has to be determined to solve the allocation problem: - The flows $F \in \mathbb{R}_{\ge 0}^{n'\times n'}$ over the edges in the allocation graph; - The allocations to the basins @@ -100,7 +116,7 @@ $$ Currently the basin allocations are not taken into account in the implementation. ::: -### The optimization objective +## The optimization objective The goal of allocation is to get the flow to the users as close as possible to their demand. To achieve this, the following objectives are supported: @@ -121,6 +137,11 @@ $$ \min \sum_{(i,j)\in E_S\;:\; i\in U_S} \left|1 - \frac{F_{ij}}{d_j^p(t)}\right| + c \sum_{e \in E_S} F_e $$ +:::{.callout-note} +When performing main network allocation, the connections to subnetworks are also interpreted as users with demands determined by subnetwork demand collection. +::: + + To avoid division by $0$ errors, if a `*_relative` objective is used and a demand is $0$, the coefficient of the flow $F_{ij}$ is set to $0$. For `*_absolute` objectives the optimizer cares about the actual amount of water allocated to a user, for `*_relative` objectives it cares about the fraction of the demand allocated to the user. For `quadratic_*` objectives the optimizer cares about avoiding large shortages, for `linear_*` objectives it treats all deviations equally. @@ -137,7 +158,7 @@ The absolute value applied here is not supported in a linear programming context In the future new optimization objectives will be introduced, for demands of basins and priorities over sources. These will be used in combination with the above, in the form of goal programming. ::: -### The optimization variable constraints +## The optimization constraints - Flow conservation: For the basins in the allocation graph we have that $$ \sum_{j=1}^{n'} F_{kj} \le \sum_{i=1}^{n'} F_{ik}, \quad \forall k \in B_S. @@ -148,6 +169,12 @@ $$ F_{ij} \le \left(C_S\right)_{ij}, \quad \forall(i,j) \in E_S. $$ {#eq-capacityconstraint} By the definition of $C_S$ this also includes the source flows. + +:::{.callout-note} +When performing subnetwork demand collection, these capacities are set to $\infty$ for edges which connect the main network to a subnetwork. +::: + + - User outflow: The outflow of the user is dictated by the inflow and the return factor: $$ F_{ik} = r_k \cdot F_{kj} \quad @@ -200,6 +227,26 @@ In the future there will be 2 more optimization solves: ## Example -:::{.callout-note} -An example with figures and data will be added here after addition of [allocation output files](https://github.com/Deltares/Ribasim/issues/659). -::: +The following is an example of an optimization problem for the example shown [here](../python/examples.ipynb#model-with-allocation): + + +```{julia} +# | code-fold: true +using Ribasim +using SQLite + +toml_path = normpath(@__DIR__, "../../generated_testmodels/allocation_example/ribasim.toml") +p = Ribasim.Model(toml_path).integrator.p + +allocation_model = p.allocation_models[1] +t = 0.0 +priority_idx = 1 + +Ribasim.set_flow!(p.graph, Ribasim.NodeID(1), Ribasim.NodeID(2), 1.0) + +Ribasim.adjust_source_flows!(allocation_model, p, priority_idx) +Ribasim.adjust_edge_capacities!(allocation_model, p, priority_idx) +Ribasim.set_objective_priority!(allocation_model, p, t, priority_idx) + +println(p.allocation_models[1].problem) +``` diff --git a/pixi.toml b/pixi.toml index 1f1d5183e..8732b2483 100644 --- a/pixi.toml +++ b/pixi.toml @@ -41,11 +41,11 @@ build-julia-docs = { cmd = "julia --project=docs docs/make.jl", depends_on = [ ] } quartodoc-build = { cmd = "quartodoc build && rm objects.json", cwd = "docs" } quarto-preview = { cmd = "quarto preview docs", depends_on = [ - "quartodoc-build", + "quartodoc-build", "generate-testmodels" ] } quarto-check = { cmd = "quarto check all", depends_on = ["quartodoc-build"] } quarto-render = { cmd = "julia --project=docs --eval 'using Pkg; Pkg.build(\"IJulia\")' && quarto render docs --to html --execute", depends_on = [ - "quartodoc-build", + "quartodoc-build", "generate-testmodels" ] } docs = { depends_on = ["build-julia-docs", "quarto-preview"] } # Lint From 58f46dd7f0910ee89461fdb2830843f9e64aecbd Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" <41898282+github-actions[bot]@users.noreply.github.com> Date: Wed, 31 Jan 2024 11:12:55 +0100 Subject: [PATCH 02/40] CompatHelper: bump compat for Dictionaries to 0.4 for package core, (keep existing compat) (#1013) This pull request changes the compat entry for the `Dictionaries` package from `0.3.25` to `0.3.25, 0.4` for package core. This keeps the compat entries for earlier versions. Note: I have not tested your package with this new compat entry. It is your responsibility to make sure that your package tests pass before you merge this pull request. Co-authored-by: CompatHelper Julia --- core/Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/core/Project.toml b/core/Project.toml index 2e22edab1..a8010212a 100644 --- a/core/Project.toml +++ b/core/Project.toml @@ -56,7 +56,7 @@ DataFrames = "1.4" DataInterpolations = "4.4" DataStructures = "0.18" Dates = "<0.0.1,1" -Dictionaries = "0.3.25" +Dictionaries = "0.3.25, 0.4" DiffEqCallbacks = "2.29.1" Documenter = "0.27,1" EnumX = "1.0" From 868f5c9ff7067b673ffd05821d667c770a0f6345 Mon Sep 17 00:00:00 2001 From: Martijn Visser Date: Wed, 31 Jan 2024 13:14:28 +0100 Subject: [PATCH 03/40] Fold docs/Project.toml into root/Project.toml (#1014) There isn't really added value to having it separately anymore, this is simpler. --- Manifest.toml | 10 ++-------- Project.toml | 28 +++++++++++++++++++++++++++- docs/Project.toml | 38 -------------------------------------- pixi.toml | 6 +++--- 4 files changed, 32 insertions(+), 50 deletions(-) delete mode 100644 docs/Project.toml diff --git a/Manifest.toml b/Manifest.toml index b7c0d6488..f161f91e0 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.10.0" manifest_format = "2.0" -project_hash = "a2e54351da64b7dac3109810459579d647b75bb7" +project_hash = "1a68754569dae741c9232c9ae2db3512f52c2b80" [[deps.ADTypes]] git-tree-sha1 = "41c37aa88889c171f1300ceac1313c06e891d245" @@ -1254,7 +1254,7 @@ version = "3.5.13" deps = ["Accessors", "Arrow", "BasicModelInterface", "CodecLz4", "CodecZstd", "ComponentArrays", "Configurations", "DBInterface", "DataInterpolations", "DataStructures", "Dates", "Dictionaries", "DiffEqCallbacks", "EnumX", "FiniteDiff", "ForwardDiff", "Graphs", "HiGHS", "IterTools", "JuMP", "Legolas", "Logging", "LoggingExtras", "MetaGraphsNext", "OrdinaryDiffEq", "PreallocationTools", "SQLite", "SciMLBase", "SparseArrays", "StructArrays", "Tables", "TerminalLoggers", "TimeZones", "TimerOutputs", "TranscodingStreams"] path = "core" uuid = "aac5e3d9-0b8f-4d4f-8241-b1a7a9632635" -version = "0.6.0" +version = "0.7.0" [[deps.RuntimeGeneratedFunctions]] deps = ["ExprTools", "SHA", "Serialization"] @@ -1662,12 +1662,6 @@ path = "build/create_binaries" uuid = "3cfb6a46-05f0-43df-bb16-bf763deb14b4" version = "0.1.0" -[[deps.docs]] -deps = ["Configurations", "DataFrames", "Dates", "Documenter", "DocumenterMarkdown", "IJulia", "InteractiveUtils", "JSON3", "Legolas", "Logging", "MarkdownTables", "OrderedCollections", "Ribasim"] -path = "docs" -uuid = "8daea9ca-fc6c-4731-aa85-717fa0b706bc" -version = "0.1.0" - [[deps.libblastrampoline_jll]] deps = ["Artifacts", "Libdl"] uuid = "8e850b90-86db-534c-a0d3-1478176c7d93" diff --git a/Project.toml b/Project.toml index 6df86f2b1..ec9276c94 100644 --- a/Project.toml +++ b/Project.toml @@ -5,9 +5,21 @@ description = "Meta-project used to share the Manifest of Ribasim and its depend [deps] Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45" BasicModelInterface = "59605e27-edc0-445a-b93d-c09a3a50b330" +Configurations = "5218b696-f38b-4ac9-8b61-a12ec717816d" +DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" +Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +DocumenterMarkdown = "997ab1e6-3595-5248-9280-8efb232c3433" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" +IJulia = "7073ff75-c697-5162-941a-fcdaad2a7d2a" Infiltrator = "5903a43b-9cc3-4c30-8d17-598619ec4e9b" +InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" +JSON3 = "0f8b85d8-7281-11e9-16c2-39a750bddbf1" +Legolas = "741b9549-f6ed-4911-9fbf-4a1c0c97f0cd" +Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" +MarkdownTables = "1862ce21-31c7-451e-824c-f20fa3f90fa2" MetaGraphsNext = "fa8bd995-216d-47f1-8a91-f3b68fbeb377" +OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" ReTestItems = "817f1d60-ba6b-4fd5-9520-3cf149f6a823" Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" @@ -18,6 +30,20 @@ Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" TestEnv = "1e6cf692-eddd-4d53-88a5-2d735e33781b" TimeZones = "f269a46b-ccf7-5d73-abea-4c690281aa53" create_binaries = "3cfb6a46-05f0-43df-bb16-bf763deb14b4" -docs = "8daea9ca-fc6c-4731-aa85-717fa0b706bc" libribasim = "f319f290-633d-4573-adfe-d6d5548b6388" ribasim_cli = "e45c999e-e944-4589-a8e6-7d0b7a60140a" + +[compat] +Configurations = "0.17" +DataFrames = "1" +Dates = "<0.0.1,1" +Documenter = "0.27,1" +DocumenterMarkdown = "0.2" +IJulia = "1" +InteractiveUtils = "<0.0.1,1" +JSON3 = "1.12" +Legolas = "0.5" +Logging = "<0.0.1,1" +MarkdownTables = "1" +OrderedCollections = "1.6" +julia = "1.10" diff --git a/docs/Project.toml b/docs/Project.toml deleted file mode 100644 index 8174017e2..000000000 --- a/docs/Project.toml +++ /dev/null @@ -1,38 +0,0 @@ -name = "docs" -uuid = "8daea9ca-fc6c-4731-aa85-717fa0b706bc" -authors = ["Deltares and contributors "] -description = "Build Ribasim docs and generate schemas" -manifest = "../Manifest.toml" -version = "0.1.0" - -[deps] -Configurations = "5218b696-f38b-4ac9-8b61-a12ec717816d" -DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" -Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" -Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -DocumenterMarkdown = "997ab1e6-3595-5248-9280-8efb232c3433" -IJulia = "7073ff75-c697-5162-941a-fcdaad2a7d2a" -InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" -JSON3 = "0f8b85d8-7281-11e9-16c2-39a750bddbf1" -Legolas = "741b9549-f6ed-4911-9fbf-4a1c0c97f0cd" -Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" -MarkdownTables = "1862ce21-31c7-451e-824c-f20fa3f90fa2" -OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" -Ribasim = "aac5e3d9-0b8f-4d4f-8241-b1a7a9632635" -SQLite = "0aa819cd-b072-5ff4-a722-6bc24af294d9" - -[compat] -Configurations = "0.17" -DataFrames = "1" -Dates = "<0.0.1,1" -Documenter = "0.27,1" -DocumenterMarkdown = "0.2" -IJulia = "1" -InteractiveUtils = "<0.0.1,1" -JSON3 = "1.12" -Legolas = "0.5" -Logging = "<0.0.1,1" -MarkdownTables = "1" -OrderedCollections = "1.6" -SQLite = "1.5.1" -julia = "1.10" diff --git a/pixi.toml b/pixi.toml index 8732b2483..fe84712ae 100644 --- a/pixi.toml +++ b/pixi.toml @@ -36,7 +36,7 @@ initialize-julia = { depends_on = [ "instantiate-julia", ] } # Docs -build-julia-docs = { cmd = "julia --project=docs docs/make.jl", depends_on = [ +build-julia-docs = { cmd = "julia --project docs/make.jl", depends_on = [ "initialize-julia", ] } quartodoc-build = { cmd = "quartodoc build && rm objects.json", cwd = "docs" } @@ -44,7 +44,7 @@ quarto-preview = { cmd = "quarto preview docs", depends_on = [ "quartodoc-build", "generate-testmodels" ] } quarto-check = { cmd = "quarto check all", depends_on = ["quartodoc-build"] } -quarto-render = { cmd = "julia --project=docs --eval 'using Pkg; Pkg.build(\"IJulia\")' && quarto render docs --to html --execute", depends_on = [ +quarto-render = { cmd = "julia --project --eval 'using Pkg; Pkg.build(\"IJulia\")' && quarto render docs --to html --execute", depends_on = [ "quartodoc-build", "generate-testmodels" ] } docs = { depends_on = ["build-julia-docs", "quarto-preview"] } @@ -88,7 +88,7 @@ test-ribasim-core-cov = { cmd = "julia --project=core --eval 'using Pkg; Pkg.tes generate-testmodels = "python utils/generate-testmodels.py" tests = { depends_on = ["lint", "test-ribasim-python", "test-ribasim-core"] } # Codegen -generate-schema = { cmd = "julia --project=docs docs/gen_schema.jl" } +generate-schema = { cmd = "julia --project docs/gen_schema.jl" } generate-python = "datamodel-codegen --output-model-type pydantic_v2.BaseModel --base-class ribasim.input_base.BaseModel --use-union-operator --use-title-as-name --use-double-quotes --disable-timestamp --use-default --strict-nullable --input-file-type=jsonschema --input docs/schema/root.schema.json --output python/ribasim/ribasim/models.py" codegen = { depends_on = ["generate-schema", "generate-python", "lint"] } # Publish From 5f5e09a4547fbea8b2c984316b96ce0e1616ac0a Mon Sep 17 00:00:00 2001 From: Bart de Koning <74617371+SouthEndMusic@users.noreply.github.com> Date: Wed, 31 Jan 2024 15:54:55 +0100 Subject: [PATCH 04/40] Remove using S for storage (use u) (#1003) Fixes https://github.com/Deltares/Ribasim/issues/968. --- docs/core/equations.qmd | 24 ++++++++++++------------ docs/core/numerics.qmd | 2 +- 2 files changed, 13 insertions(+), 13 deletions(-) diff --git a/docs/core/equations.qmd b/docs/core/equations.qmd index cc238d8a6..92d5d258f 100644 --- a/docs/core/equations.qmd +++ b/docs/core/equations.qmd @@ -42,7 +42,7 @@ $\mathbf{u}(t)$ is given by all the states of the model, which are (currently) t Given a single basin with node ID $i \in B$, the equation that dictates the change of its storage over time is given by $$ -\frac{\text{d}S_i}{\text{d}t} = +\frac{\text{d}u_i}{\text{d}t} = \sum_{(i',j') \in E | j' = i} Q_{i',j'} - \sum_{(i',j') \in E | i' = i} Q_{i',j'} + F_i(p,t). $$ @@ -74,15 +74,15 @@ $$ Most of these terms are always $0$, because a flow over an edge only depends on a small number of states. Therefore the matrix $J$ is very sparse. -For many contributions to the Jacobian the derivative of the level $l(S)$ of a basin with respect to its storage $S$ is required. To get an expression for this, we first look at the storage as a function of the level: +For many contributions to the Jacobian the derivative of the level $l(u)$ of a basin with respect to its storage $u$ is required. To get an expression for this, we first look at the storage as a function of the level: $$ -S(l) = \int_{l_0}^l A(\ell)d\ell. +u(l) = \int_{l_0}^l A(\ell)d\ell. $$ -From this we obtain $S'(l) = A(l)$ and thus +From this we obtain $u'(l) = A(l)$ and thus $$ -\frac{\text{d}l}{\text{d}S} = \frac{1}{A(S)}. +\frac{\text{d}l}{\text{d}u} = \frac{1}{A(u)}. $$ :::{.callout-note} @@ -149,11 +149,11 @@ plt.show() The precipitation term is given by $$ - Q_P = P \cdot A(S). + Q_P = P \cdot A(u). $$ {#eq-precip} Here $P = P(t)$ is the precipitation rate and $A$ is the wetted area. $A$ is a -function of the storage $S = S(t)$: as the volume of water changes, the area of the free water +function of the storage $u = u(t)$: as the volume of water changes, the area of the free water surface may change as well, depending on the slopes of the surface waters. ## Evaporation @@ -161,13 +161,13 @@ surface may change as well, depending on the slopes of the surface waters. The evaporation term is given by $$ - Q_E = E_\text{pot} \cdot A(S) \cdot \phi(d;0.1). + Q_E = E_\text{pot} \cdot A(u) \cdot \phi(d;0.1). $$ {#eq-evap} -Here $E_\text{pot} = E_\text{pot}(t)$ is the potential evaporation rate and $A$ is the wetted area. $\phi$ is the [reduction factor](equations.qmd#sec-reduction_factor) which depends on the depth $d$. It provides a smooth gradient as $S \rightarrow 0$. +Here $E_\text{pot} = E_\text{pot}(t)$ is the potential evaporation rate and $A$ is the wetted area. $\phi$ is the [reduction factor](equations.qmd#sec-reduction_factor) which depends on the depth $d$. It provides a smooth gradient as $u \rightarrow 0$. -A straightforward formulation $Q_E = \mathrm{max}(E_\text{pot} A(S), -0)$ is unsuitable, as $\frac{\mathrm{d}Q_E}{\mathrm{d}S}(S=0)$ is then not well-defined. +A straightforward formulation $Q_E = \mathrm{max}(E_\text{pot} A(u), +0)$ is unsuitable, as $\frac{\mathrm{d}Q_E}{\mathrm{d}u}(u=0)$ is then not well-defined.