Some fixes

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"

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"

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"

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
 

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

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,20 +27,20 @@ $$
 	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),
 $$
 
 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).

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