Fix CGEuler not respecting dtmax and dtmin

Project.toml

 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"

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

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

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

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)}()

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