Deltares · SouthEndMusic · Apr 16, 2026 · Apr 16, 2026 · Apr 16, 2026 · Apr 17, 2026
diff --git a/Manifest.toml b/Manifest.toml
diff --git a/Project.toml b/Project.toml
@@ -9,6 +9,7 @@ projects = ["core"]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
 Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
 BasicModelInterface = "59605e27-edc0-445a-b93d-c09a3a50b330"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"

diff --git a/core/Project.toml b/core/Project.toml
@@ -9,6 +9,7 @@ projects = ["test"]
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
+ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
 BasicModelInterface = "59605e27-edc0-445a-b93d-c09a3a50b330"
 Configurations = "5218b696-f38b-4ac9-8b61-a12ec717816d"
 DBInterface = "a10d1c49-ce27-4219-8d33-6db1a4562965"
@@ -60,6 +61,7 @@ TranscodingStreams = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
 [compat]
 ADTypes = "1.13"
 Accessors = "0.1"
+ArrayInterface = "7.27.0"
 BasicModelInterface = "0.1"
 Configurations = "0.17"
 DBInterface = "2.4"

diff --git a/core/src/Ribasim.jl b/core/src/Ribasim.jl
@@ -22,23 +22,28 @@ using PrecompileTools: @setup_workload, @compile_workload
 using DifferentiationInterface:
     AutoSparse,
     Constant,
-    Cache,
     prepare_jacobian,
+    value_and_jacobian!,
     jacobian!,
     prepare_derivative,
     derivative!,
     second_derivative
 
-using ForwardDiff: derivative as forward_diff
+using ForwardDiff: ForwardDiff, Dual, Partials, seed!, partials, valtype, derivative as forward_diff
+
+using ArrayInterface: ArrayInterface
 
 # Algorithms for solving ODEs.
-using OrdinaryDiffEqCore: OrdinaryDiffEqCore, get_du
-using OrdinaryDiffEqDifferentiation:
-    OrdinaryDiffEqDifferentiation, dolinsolve, jacobian2W!
-using SciMLOperators: WOperator, MatrixOperator
+using OrdinaryDiffEqCore:
+    OrdinaryDiffEqCore,
+    OrdinaryDiffEqAdaptiveImplicitAlgorithm,
+    OrdinaryDiffEqImplicitAlgorithm,
+    get_du
+using DiffEqBase: DiffEqBase
 import ADTypes
 using ADTypes: AutoForwardDiff
 import ForwardDiff
+import NaNMath
 import OrdinaryDiffEqBDF
 
 # Interface for defining and solving the ODE problem of the physical layer.
@@ -54,21 +59,22 @@ using SciMLBase:
     DEIntegrator,
     FullSpecialize,
     NoSpecialize,
-    SciMLOperators,
-    AbstractSciMLOperator,
     LinearProblem,
     LinearSolution
 
+using OrdinaryDiffEqDifferentiation:
+    OrdinaryDiffEqDifferentiation, do_newJW, jacobian2W!, dolinsolve
+
 # Automatically detecting the sparsity pattern of the Jacobian of water_balance!
 # through operator overloading
-using SparseConnectivityTracer: GradientTracer, TracerSparsityDetector
+using SparseConnectivityTracer: TracerSparsityDetector
 using SparseMatrixColorings: GreedyColoringAlgorithm, sparsity_pattern
 
 # For efficient sparse computations
-using SparseArrays: SparseMatrixCSC, sparse, nzrange
+using SparseArrays: SparseMatrixCSC, spzeros, sparse, findnz
 
 # Linear algebra
-using LinearAlgebra: LinearAlgebra, mul!, UniformScaling
+using LinearAlgebra: LinearAlgebra, mul!, dot, I
 
 # Interpolation functionality, used for e.g.
 # basin profiles and TabulatedRatingCurve. See also the node
@@ -108,7 +114,7 @@ import BasicModelInterface as BMI
 using DelimitedFiles: writedlm
 
 # Reading GeoPackage files, which are SQLite databases with spatial data
-using SQLite: SQLite, DB, Query, esc_id
+using SQLite: SQLite, DB, esc_id
 using DBInterface: execute, prepare
 
 # Logging to both the console and a file
@@ -169,13 +175,17 @@ using Printf: @sprintf
 
 using Base.Threads: nthreads
 
-include("carrays.jl")
-using .CArrays: CVector, getaxes, getdata
+using SciMLOperators: SciMLOperators, WOperator, AbstractSciMLOperator
+
+include("cvectors.jl")
+using .CVectors: CVector, getaxes, getdata
 include("schema.jl")
 include("config.jl")
 using .config
+using .config: with_mass_matrix
 include("parameter.jl")
 include("validation.jl")
+include("formulate_flows.jl")
 include("solve.jl")
 include("logo.jl")
 include("logging.jl")
@@ -184,7 +194,6 @@ include("allocation_init.jl")
 include("allocation_optim.jl")
 include("util.jl")
 include("graph.jl")
-include("differentiation.jl")
 include("model.jl")
 include("read.jl")
 include("write.jl")

diff --git a/core/src/allocation_init.jl b/core/src/allocation_init.jl
@@ -167,7 +167,7 @@ function add_user_demand!(
         allocation_model::AllocationModel,
         p_independent::ParametersIndependent,
     )::Nothing
-    (; problem, cumulative_supplied_volume, node_ids_in_subnetwork) = allocation_model
+    (; problem, node_ids_in_subnetwork) = allocation_model
     (; user_demand_ids_subnetwork) = node_ids_in_subnetwork
     (; user_demand) = p_independent
     (; inflow_links, outflow_link) = user_demand
@@ -235,14 +235,6 @@ function add_user_demand!(
             return_factor * sum(flow[lm.link] for lm in inflow_links[node_id.idx]),
         base_name = "user_demand_return_flow"
     )
-
-    # Add the links for which the supplied volume is required for output
-    for node_id in user_demand_ids_subnetwork
-        for link_metadata in inflow_links[node_id.idx]
-            cumulative_supplied_volume[link_metadata.link] = 0.0
-        end
-    end
-
     return nothing
 end
 
@@ -253,11 +245,11 @@ function add_flow_demand!(
         allocation_model::AllocationModel,
         p_independent::ParametersIndependent,
     )::Nothing
-    (; problem, cumulative_supplied_volume, scaling, node_ids_in_subnetwork) =
+    (; problem, scaling, node_ids_in_subnetwork) =
         allocation_model
-    (; node_ids_subnetwork_with_flow_demand, flow_demand_ids_subnetwork) =
+    (; node_ids_subnetwork_with_flow_demand) =
         node_ids_in_subnetwork
-    (; graph, flow_demand) = p_independent
+    (; graph) = p_independent
     flow = problem[:flow]
 
     # Define decision variables: flow allocated to FlowDemand node per demand priority
@@ -327,11 +319,6 @@ function add_flow_demand!(
             d - flow_demand_allocated[node_id, demand_priority],
         base_name = "flow_demand_relative_error_constraint"
     )
-
-    # Add the links for which the supplied volume is required for output
-    for node_id in flow_demand_ids_subnetwork
-        cumulative_supplied_volume[flow_demand.inflow_link[node_id.idx].link] = 0.0
-    end
     return nothing
 end
 
@@ -433,15 +420,15 @@ function add_linearized_connector_node!(
 
             # Only linearize if the level comes from a Basin
             upstream_node = inflow_link[node_id.idx].link[1]
-            if upstream_node.type == NodeType.Basin
+            if upstream_node.is_basin
                 JuMP.add_to_expression!(
                     linearization,
                     ∂q∂h_upstream * storage_change[upstream_node] / A,
                 )
             end
 
             downstream_node = outflow_link[node_id.idx].link[2]
-            if downstream_node.type == NodeType.Basin
+            if downstream_node.is_basin
                 JuMP.add_to_expression!(
                     linearization,
                     ∂q∂h_downstream * storage_change[downstream_node] / A,
@@ -926,13 +913,13 @@ function NodeIDsInSubnetwork(
         # basin_ids_subnetwork_with_level_demand
         get_nodes(
             node_id ->
-            node_id.type == NodeType.Basin &&
+            node_id.is_basin &&
                 !isnothing(get_external_demand_id(p_independent, node_id)),
         ),
         # node_ids_subnetwork_with_flow_demand
         get_nodes(
             node_id ->
-            node_id.type != NodeType.Basin &&
+            !node_id.is_basin &&
                 !isnothing(get_external_demand_id(p_independent, node_id)),
         ),
     )