From b1002b95b9e455871c10c89aa5e4af9ac5659b8b Mon Sep 17 00:00:00 2001
From: Brian Groenke <brian.groenke@awi.de>
Date: Thu, 26 Dec 2024 21:02:01 +0100
Subject: [PATCH] Fix CGEuler not respecting dtmax and dtmin

---
 Project.toml                             |  2 +-
 examples/heat_sfcc_samoylov_implicitT.jl |  2 +-
 src/Solvers/LiteImplicit/cglite_types.jl | 13 +++++++++++--
 src/Solvers/basic_solvers.jl             |  2 +-
 src/Solvers/integrator.jl                |  2 +-
 src/problem.jl                           |  7 +++++--
 6 files changed, 20 insertions(+), 8 deletions(-)

diff --git a/Project.toml b/Project.toml
index 8a45647e..e6a046dd 100755
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "CryoGrid"
 uuid = "a535b82e-5f3d-4d97-8b0b-d6483f5bebd5"
 authors = ["Brian Groenke <brian.groenke@awi.de>", "Jan Nitzbon <jan.nitzbon@awi.de>", "Moritz Langer <moritz.langer@awi.de>"]
-version = "0.23.1"
+version = "0.23.2"
 
 [deps]
 Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
diff --git a/examples/heat_sfcc_samoylov_implicitT.jl b/examples/heat_sfcc_samoylov_implicitT.jl
index 469d71c3..6c355aa3 100644
--- a/examples/heat_sfcc_samoylov_implicitT.jl
+++ b/examples/heat_sfcc_samoylov_implicitT.jl
@@ -38,7 +38,7 @@ tile = Tile(strat, modelgrid, forcings, ssinit);
 tspan = (DateTime(2000,10,1), DateTime(2002,10,1))
 u0, du0 = initialcondition!(tile, tspan);
 prob = CryoGridProblem(tile, u0, tspan, saveat=24*3600, savevars=(:T,:θw,:∂H∂T), step_limiter=nothing)
-integrator = init(prob, ImplicitEuler(autodiff=false))
+integrator = init(prob, ImplicitEuler())
 # test one step
 @time step!(integrator)
 last_t = integrator.t
diff --git a/src/Solvers/LiteImplicit/cglite_types.jl b/src/Solvers/LiteImplicit/cglite_types.jl
index 47c993f9..11d3253f 100644
--- a/src/Solvers/LiteImplicit/cglite_types.jl
+++ b/src/Solvers/LiteImplicit/cglite_types.jl
@@ -22,7 +22,16 @@ struct LiteImplicitEulerCache{Tu,TA} <: SciMLBase.DECache
     D::TA
 end
 
-function DiffEqBase.__init(prob::CryoGridProblem, alg::LiteImplicitEuler, args...; dt=24*3600.0, saveat=dt, kwargs...)
+function DiffEqBase.__init(
+    prob::CryoGridProblem,
+    alg::LiteImplicitEuler,
+    args...;
+    dt=24*3600.0,
+    dtmax=dt,
+    dtmin=1.0,
+    saveat=dt,
+    kwargs...
+)
     tile = Tile(prob.f)
     grid = tile.grid
     u0 = copy(prob.u0)
@@ -54,6 +63,6 @@ function DiffEqBase.__init(prob::CryoGridProblem, alg::LiteImplicitEuler, args..
         similar(u0, eltype(u0), length(prob.u0.H)),
     )
     p = isnothing(prob.p) ? prob.p : collect(prob.p)
-    opts = CryoGridIntegratorOptions(; dtmax=dt, saveat=CryoGrid.expandtstep(saveat, prob.tspan), kwargs...)
+    opts = CryoGridIntegratorOptions(; saveat=CryoGrid.expandtstep(saveat, prob.tspan), dtmax, dtmin, kwargs...)
     return CryoGridIntegrator(alg, cache, opts, sol, copy(u0), p, t0, convert(eltype(prob.tspan), dt), 1, 1)
 end
diff --git a/src/Solvers/basic_solvers.jl b/src/Solvers/basic_solvers.jl
index 5e9afe2d..c03a4b99 100644
--- a/src/Solvers/basic_solvers.jl
+++ b/src/Solvers/basic_solvers.jl
@@ -74,6 +74,6 @@ function perform_step!(integrator::CryoGridIntegrator{CGEuler})
     # set next dt
     # compute maximum timestep
     dtmax = CryoGrid.timestep(tile, du, u, p, t₀)
-    integrator.dt = dtmax
+    integrator.dt = max(min(dtmax, integrator.opts.dtmax), integrator.opts.dtmin)
     return nothing
 end
diff --git a/src/Solvers/integrator.jl b/src/Solvers/integrator.jl
index 53ed2732..195bab2b 100755
--- a/src/Solvers/integrator.jl
+++ b/src/Solvers/integrator.jl
@@ -52,7 +52,7 @@ function DiffEqBase.sensitivity_solution(sol::CryoGridSolution, u, t)
 end
 
 Base.@kwdef mutable struct CryoGridIntegratorOptions{Tt,Tsaveat,Ttstops}
-    dtmax::Tt = 24*3600.0
+    dtmax::Tt = 3600.0
     dtmin::Tt = one(typeof(dtmax))
     saveat::Tsaveat = typeof(dtmax)[]
     tstops::Ttstops = SortedSet{typeof(dtmin)}()
diff --git a/src/problem.jl b/src/problem.jl
index 6c05881f..db3ebf19 100644
--- a/src/problem.jl
+++ b/src/problem.jl
@@ -47,7 +47,8 @@ Constructor for `CryoGridProblem` that automatically generates all necessary cal
 function CryoGridProblem(
     tile::Tile,
     u0::ComponentVector,
-    tspan::NTuple{2,Float64};
+    tspan::NTuple{2,Float64},
+    p::Union{Nothing,AbstractVector}=nothing;
     diagnostic_stepsize=3600.0,
     saveat=3600.0,
     savevars=(),
@@ -68,7 +69,7 @@ function CryoGridProblem(
     # strip all "fixed" parameters
     tile = stripparams(FixedParam, tile)
     # retrieve variable parameters
-    p = length(ModelParameters.params(tile)) > 0 ? parameters(tile) : nothing
+    tilepara = length(ModelParameters.params(tile)) > 0 ? parameters(tile) : nothing
     du0 = zero(u0)
     # remove units
     tile = stripunits(tile)
@@ -95,6 +96,8 @@ function CryoGridProblem(
     tile.data.outputs = savevals
     # build mass matrix
     mass_matrix = Numerics.build_mass_matrix(tile.state)
+    # get params
+    p = isnothing(p) && !isnothing(tilepara) ? ustrip.(vec(tilepara)) : p
 	func = odefunction(tile, u0, p, tspan; mass_matrix, specialization, function_kwargs...)
 	return CryoGridProblem{true}(func, u0, tspan, p, callbacks, saveat, getsavestate, isoutofdomain, prob_kwargs)
 end