From 380ad01ee451423c6eb074c01f2172a013cd0afb Mon Sep 17 00:00:00 2001 From: Bart de Koning Date: Tue, 17 Sep 2024 19:18:11 +0200 Subject: [PATCH] Some fixes --- Manifest.toml | 8 ++++---- Project.toml | 2 +- core/Project.toml | 4 ++-- core/src/Ribasim.jl | 1 + core/src/parameter.jl | 5 +++++ docs/concept/equations.qmd | 25 +++++++++++++------------ docs/concept/modelconcept.qmd | 2 +- 7 files changed, 27 insertions(+), 20 deletions(-) diff --git a/Manifest.toml b/Manifest.toml index d03ba3a69..425ca14c5 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.10.5" manifest_format = "2.0" -project_hash = "d694251e41cfe037d14e0d129b2e2e49035c0990" +project_hash = "f8b2dea16bf0ccef5eef2111258ff045598662c9" [[deps.ADTypes]] git-tree-sha1 = "99a6f5d0ce1c7c6afdb759daa30226f71c54f6b0" @@ -375,9 +375,9 @@ uuid = "ade2ca70-3891-5945-98fb-dc099432e06a" [[deps.DiffEqBase]] deps = ["ArrayInterface", "ConcreteStructs", "DataStructures", "DocStringExtensions", "EnumX", "EnzymeCore", "FastBroadcast", "FastClosures", "ForwardDiff", "FunctionWrappers", "FunctionWrappersWrappers", "LinearAlgebra", "Logging", "Markdown", "MuladdMacro", "Parameters", "PreallocationTools", "PrecompileTools", "Printf", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "SciMLStructures", "Setfield", "Static", "StaticArraysCore", "Statistics", "Tricks", "TruncatedStacktraces"] -git-tree-sha1 = "d7d43a1cc11dc4e4e5378816ae720fccd75a77c8" +git-tree-sha1 = "fa7d580038451a8df4434ef5b079ac9b2d486194" uuid = "2b5f629d-d688-5b77-993f-72d75c75574e" -version = "6.155.0" +version = "6.155.1" [deps.DiffEqBase.extensions] DiffEqBaseCUDAExt = "CUDA" @@ -1344,7 +1344,7 @@ uuid = "295af30f-e4ad-537b-8983-00126c2a3abe" version = "3.5.18" [[deps.Ribasim]] -deps = ["Accessors", "Arrow", "BasicModelInterface", "CodecZstd", "ComponentArrays", "Configurations", "DBInterface", "DataInterpolations", "DataStructures", "Dates", "DiffEqCallbacks", "EnumX", "FiniteDiff", "ForwardDiff", "Graphs", "HiGHS", "IterTools", "JuMP", "Legolas", "LineSearches", "LinearSolve", "Logging", "LoggingExtras", "MetaGraphsNext", "OrdinaryDiffEqBDF", "OrdinaryDiffEqCore", "OrdinaryDiffEqLowOrderRK", "OrdinaryDiffEqNonlinearSolve", "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqSDIRK", "OrdinaryDiffEqTsit5", "PreallocationTools", "SQLite", "SciMLBase", "SparseArrays", "SparseConnectivityTracer", "StructArrays", "Tables", "TerminalLoggers", "TranscodingStreams"] +deps = ["Accessors", "Arrow", "BasicModelInterface", "CodecZstd", "ComponentArrays", "Configurations", "DBInterface", "DataInterpolations", "DataStructures", "Dates", "DiffEqBase", "DiffEqCallbacks", "EnumX", "FiniteDiff", "Graphs", "HiGHS", "IterTools", "JuMP", "Legolas", "LineSearches", "LinearSolve", "Logging", "LoggingExtras", "MetaGraphsNext", "OrdinaryDiffEqBDF", "OrdinaryDiffEqCore", "OrdinaryDiffEqLowOrderRK", "OrdinaryDiffEqNonlinearSolve", "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqSDIRK", "OrdinaryDiffEqTsit5", "PreallocationTools", "SQLite", "SciMLBase", "SparseArrays", "SparseConnectivityTracer", "StructArrays", "Tables", "TerminalLoggers", "TranscodingStreams"] path = "core" uuid = "aac5e3d9-0b8f-4d4f-8241-b1a7a9632635" version = "2024.10.0" diff --git a/Project.toml b/Project.toml index f8fd65409..ef7bcf793 100644 --- a/Project.toml +++ b/Project.toml @@ -15,11 +15,11 @@ DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" +DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" EnumX = "4e289a0a-7415-4d19-859d-a7e5c4648b56" FindFirstFunctions = "64ca27bc-2ba2-4a57-88aa-44e436879224" FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" -ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" Glob = "c27321d9-0574-5035-807b-f59d2c89b15c" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" diff --git a/core/Project.toml b/core/Project.toml index 521d86e91..76be2b5cd 100644 --- a/core/Project.toml +++ b/core/Project.toml @@ -15,10 +15,10 @@ DBInterface = "a10d1c49-ce27-4219-8d33-6db1a4562965" DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" +DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def" EnumX = "4e289a0a-7415-4d19-859d-a7e5c4648b56" FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" -ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6" HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e" @@ -67,10 +67,10 @@ DataFrames = "1.4" DataInterpolations = "6" DataStructures = "0.18" Dates = "<0.0.1, 1" +DiffEqBase = "6.155" DiffEqCallbacks = "3.6" EnumX = "1.0" FiniteDiff = "2.21" -ForwardDiff = "0.10" Graphs = "1.9" HiGHS = "1.7" IOCapture = "0.2" diff --git a/core/src/Ribasim.jl b/core/src/Ribasim.jl index 06c1bf7c3..0c0ab18f5 100644 --- a/core/src/Ribasim.jl +++ b/core/src/Ribasim.jl @@ -21,6 +21,7 @@ using OrdinaryDiffEqCore: get_du, AbstractNLSolver, calculate_residuals! +using DiffEqBase: DiffEqBase using OrdinaryDiffEqNonlinearSolve: OrdinaryDiffEqNonlinearSolve, relax!, _compute_rhs! using LineSearches: BackTracking diff --git a/core/src/parameter.jl b/core/src/parameter.jl index f867e46c6..bfaeb1fbe 100644 --- a/core/src/parameter.jl +++ b/core/src/parameter.jl @@ -842,3 +842,8 @@ const ModelGraph = MetaGraph{ water_balance_abstol::Float64 water_balance_reltol::Float64 end + +# To opt-out of type checking for ForwardDiff +function DiffEqBase.anyeltypedual(::Parameters, ::Type{Val{counter}}) where {counter} + Any +end diff --git a/docs/concept/equations.qmd b/docs/concept/equations.qmd index 2f5c45776..ed03e89c6 100644 --- a/docs/concept/equations.qmd +++ b/docs/concept/equations.qmd @@ -11,13 +11,13 @@ $$ which is a system of coupled first order differential equations. The model is given by a directed graph, consisting of a set of node IDs (vertices) $V$ and edges $E$, consisting of ordered pairs of node IDs. -We denote the subset of the nodes given by the basins $B \subset V$, and the subset of nodes that prescribe flow $N \subset V$. +We denote the subset of the nodes given by the Basins $B \subset V$, and the subset of nodes that prescribe flow $N \subset V$. The states $\mathbf{u}$ of the model are given by cumulative flows since the start of the simulation as prescribed by the nodes $N$: $$ u_n(t) = \int_{t_0}^t q_n\text{d}t' \quad \forall n \in N, $$ -as well as by the basin forcings: +as well as by the Basin forcings: $$ u_b^\text{forcing}(t) = \int_{t_0}^t q_b^\text{forcing}\text{d}t' \quad \forall b \in B. $$ @@ -27,7 +27,7 @@ $$ u_i(t_0) = 0 \quad \forall i. $$ -From these cumulative flows, the storage in each basin can be determined at each point in time: +From these cumulative flows, the storage in each Basin can be determined at each point in time: $$ S_b(t) = S_i(0) + S^\text{exact}(t) - u_b^\text{forcing}(t) + \sum_{n\;|\;(n,b)\in E} u(t) - \sum_{n\;|\;(b,n)\in E} u(t), $$ @@ -35,12 +35,12 @@ $$ i.e. the storage is given by: - the initial storage; -- Plus the exactly integrated flows (more on that below); +- plus the exactly integrated flows (more on that below); - minus the cumulative outgoing forcings; - plus the cumulative horizontal inflows; - minus the cumulative horizontal outflows. -From these storages in combination with the basin profiles the basin levels $h$ are computed. +From these storages in combination with the Basin profiles the Basin levels $h$ are computed. The relationship between the profile and the storage is given by $$ S_b = \int_{h_0}^h A_b(\ell)\text{d}\ell, @@ -52,7 +52,7 @@ These levels are then inputs for determining the flows prescribed by the nodes $ $$ \frac{\text{d}h}{\text{d}t} = \frac{1}{A_b}, $$ -and so areas of zero are not allowed in the basin profiles. +and so areas of zero are not allowed in the Basin profiles. ## The PID control integral state @@ -64,7 +64,7 @@ This is the error integral for PID control, further explained in [PID equations] The more states the problem has, the more time it takes to solve it. Therefore we want to minimize the number of states. Flows do not have to be states when they can be integrated over time exactly because they do not depend on the other states. -This is true for flow boundaries, precipitation (which uses a fixed basin area) and drainage. +This is true for FlowBoundary nodes, and Basin precipitation (which uses a fixed basin area) and drainage. ## The Jacobian @@ -81,16 +81,17 @@ i.e. $J_{i,j}$ quantifies how $f_j$. the time derivative of state $j$, changes w ## Why this formulation -You might wonder why in the above explanation the states are given by the cumulative flows and not by the basin storages, which is arguably conceptually simpler. -The reason is that we do not just want to model the storages in the basins over time, but we also want accurate output of each individual flow, e.g. to model the spread of pollutants. +You might wonder why in the above explanation the states are given by the cumulative flows and not by the Basin storages, which is arguably conceptually simpler. +The reason is that we do not just want to model the storages in the Basins over time, but we also want accurate output of each individual flow, e.g. to model the spread of pollutants. -When the states are given by the storages, generally the individual flows can not accurately be computed from that as a post processing step. +When the states are given by the storages, generally the individual flows can not accurately be computed from that as a post processing step, because there are more flows than storages. +Also, we can only compute flows at individual points in time explicitly, not over a whole interval. When the states are given by the cumulative flows however, the output of the problem solve gives these flows directly, and from those the storage over time can be computed accurately. -Hence in short, the formulation above gives more information than a formulation with basin storages as states. +Hence in short, the formulation above gives more information than a formulation with Basin storages as states. ## Numerical solution -Ribasim uses `OrdinaryDiffEq.jl` to provide a numerical solution to the water balance equations. +Ribasim uses [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl/) to provide a numerical solution to the water balance equations. Changes to forcings or parameters such as precipitation, but also the allocated water abstraction is managed through the use of callback functions [@callbacks]. In a coupled run, the exchanges with MODFLOW 6 are also managed via the use of a callback function. For more a more in-depth discussion of numerical computations see [Numerical considerations](numerics.qmd). diff --git a/docs/concept/modelconcept.qmd b/docs/concept/modelconcept.qmd index bb0c83af2..c9e46bea9 100644 --- a/docs/concept/modelconcept.qmd +++ b/docs/concept/modelconcept.qmd @@ -27,7 +27,7 @@ $$ $$ We don't use these equations directly. -Rather, we use an equivalent formulation where solve for the cumulative flows instead of the basin storages. +Rather, we use an equivalent formulation where solve for the cumulative flows instead of the Basin storages. For more details on this see [Equations](equations.qmd). ## Time