Tries

core/src/Ribasim.jl

@@ -57,7 +57,11 @@ using MetaGraphsNext:
 using NonlinearSolve: solve
 using PreallocationTools: LazyBufferCache
 using OrdinaryDiffEq:
-    OrdinaryDiffEq, OrdinaryDiffEqRosenbrockAdaptiveAlgorithm, get_du, SteadyStateProblem
+    OrdinaryDiffEq,
+    OrdinaryDiffEqRosenbrockAdaptiveAlgorithm,
+    get_du,
+    SteadyStateProblem,
+    ImplicitEuler
 using SciMLBase:
     init,
     solve!,
@@ -72,7 +76,7 @@ using SciMLBase:
     ODESolution,
     VectorContinuousCallback,
     get_proposed_dt
-using SteadyStateDiffEq: DynamicSS
+using SteadyStateDiffEq: DynamicSS, SSRootfind
 using SQLite: SQLite, DB, Query, esc_id
 using StructArrays: StructVector
 using Tables: Tables, AbstractRow, columntable

core/src/model.jl

@@ -146,8 +146,14 @@ function Model(config::Config)::Model
             (; basin) = parameters
 
             # Compute steady state
+            println("Computing steady state...")
             parameters.fixed_t[] = 0.0
-            sol = solve(SteadyStateProblem(prob), DynamicSS())
+            sol = solve(
+                SteadyStateProblem(RHS, u0, parameters),
+                DynamicSS(ImplicitEuler());
+                abstol = 1e-10,
+                reltol = 1e-10,
+            )
             parameters.fixed_t[] = -1.0
 
             # Check steady state

core/src/solve.jl

@@ -54,6 +54,12 @@ function water_balance!(
     # Formulate du (controlled by PidControl)
     formulate_du_pid_controlled!(du, graph, pid_control)
 
+    if eltype(du) == Float64
+        println(p.basin.node_id[argmax(du.storage)])
+    end
+    #println(sum(u.storage))
+    #println("$(sqrt(sum(du .^ 2))),")
+
     return nothing
 end
 

core/src/util.jl

@@ -846,7 +846,7 @@ function get_jac_prototype(du0, u0, p, t0)
         TracerSparsityDetector(),
     )
     p.all_nodes_active[] = false
-    jac_prototype
+    Float64.(jac_prototype)
 end
 
 # Custom overloads