AA/adaptive-dt

CHANGELOG.md

@@ -6,7 +6,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## Version [v0.5.3] - 2026-02-19
 
 ### Added
-*  No new functionality added
+*  Added optional adaptive time stepping based on Courant number control (`AdaptiveTimeStepping`) [#98](@ref)
 
 
 ### Fixed

examples/2D_flatplate_fixedT_KOmega.jl

@@ -99,7 +99,8 @@ solvers = (
         solver      = Gmres(), # Bicgstab(), Gmres()
         preconditioner = Jacobi(),
         convergence = 1e-7,
-        relax       = 0.3,
+        # relax       = 0.3,
+        relax       = 0.2,
         atol = 1e-2
     ),
     h = SolverSetup(

src/Calculate/Calculate_0_gradient.jl

@@ -255,4 +255,9 @@ function grad!(grad::Grad{Midpoint,F,R,I,M}, phif, phi, BCs, time, config) where
         correct_interpolation!(grad, phif, phi, config)
         green_gauss!(grad, phif, config)
     end
+end
+
+# This is a gradient method that does not take BCs as arguments (implicitly assuming the the gradient at boundaries is 0)
+function grad!(grad::Grad{Gauss,F,R,I,M}, phif, phi, time, config) where {F,R<:VectorField,I,M}
+    green_gauss!(grad, phif, config)
 end

src/Discretise/Discretise_1_schemes.jl

@@ -30,7 +30,7 @@ end
 @inline scheme_source!(
     term::Operator{F,P,I,Time{Euler}}, cell, cID, cIndex, prev, runtime)  where {F,P,I} = begin
         volume = cell.volume
-        vol_rdt = volume/runtime.dt
+        vol_rdt = volume/runtime.dt[1]
 
         # Increment sparse and b arrays 
         ac = vol_rdt
@@ -48,7 +48,7 @@ end
 @inline scheme_source!(
     term::Operator{F,P,I,Time{CrankNicolson}}, cell, cID, cIndex, prev, runtime)  where {F,P,I} = begin
         volume = cell.volume
-        vol_rdt = volume/runtime.dt
+        vol_rdt = volume/runtime.dt[1]
 
         # Increment sparse and b arrays 
         ac = vol_rdt

src/ModelPhysics/0_type_definition.jl

@@ -4,6 +4,7 @@ export Momentum
 export AbstractTimeModel
 export Transient, Steady
 
+
 """
     struct Physics{T,F,SO,M,Tu,E,D,BI}
         time::T
@@ -159,4 +160,4 @@ end
 mutable struct ModelState{T1}
     residuals::T1
     converged::Bool
-end
+end

src/ModelPhysics/Energy/Energy.jl

@@ -19,7 +19,6 @@ include("energy_types.jl")
 
 # Energy models
 include("Sensible_Enthalpy.jl")
-# include("Laplace_Energy.jl")
 include("Conduction.jl")
 
 # Property Models

src/ModelPhysics/Energy/Sensible_Enthalpy.jl

@@ -159,7 +159,7 @@ function energy!(
     # Kbounded = ScalarField(mesh)
     Pr = model.fluid.Pr
 
-    dt = runtime.dt
+    dt = runtime.dt[1]
 
     # Pre-allocate auxiliary variables
     TF = _get_float(mesh)

src/ModelPhysics/ModelPhysics.jl

@@ -34,4 +34,5 @@ include("Turbulence/Turbulence.jl")
 include("FluidProperties/FluidProperties.jl")
 
 
+
 end # end module

src/ModelPhysics/Turbulence/RANS_laminar.jl

@@ -102,6 +102,19 @@ function save_output(model::Physics{T,F,SO,M,Tu,E,D,BI}, outputWriter, iteration
     write_results(iteration, time, model.domain, outputWriter, config.boundaries, args...)
 end
 
+# function save_output(model::Physics{T,F,SO,M,Tu,E,D,BI}, outputWriter, iteration, time, config
+#     ) where {T,F<:Multiphase,SO,M,Tu<:Laminar,E<:Nothing,D,BI}
+
+#     args = (
+#         ("U", model.momentum.U), 
+#         ("p", model.momentum.p),
+#         ("alpha", model.fluid.alpha),
+#         ("rho", model.fluid.rho),
+#         # ("p_rgh", model.fluid.p_rgh)
+#     )
+#     write_results(iteration, time, model.domain, outputWriter, config.boundaries, args...)
+# end
+
 function save_output(model::Physics{T,F,SO,M,Tu,E,D,BI}, outputWriter, iteration, time, config
     ) where {T,F,SO,M,Tu<:Laminar,E<:Nothing,D,BI}
     args = (

src/Simulate/Simulate_0_types.jl

@@ -2,12 +2,12 @@ export Configuration
 export Hardware
 
 """
-    @kwdef struct Configuration{SC,SL,RT,HW,PP}
+    struct Configuration{SC,SL,RT,HW,PP}
         schemes::SC
         solvers::SL
         runtime::RT
         hardware::HW
-        postprocess::PP = nothing
+        postprocess::PP
     end
 
 The `Configuration` type is passed to all flow solvers and provides all the relevant information to run a simulation. 
@@ -29,16 +29,29 @@ config = Configuration(
     solvers=solvers, schemes=schemes, runtime=runtime, hardware=hardware, boundaries=BCs)
 ```
 """
-@kwdef struct Configuration{SC,SL,RT,HW,BC,PP}
+struct Configuration{SC,SL,RT,HW,BC,PP}
     schemes::SC
     solvers::SL
     runtime::RT
     hardware::HW
     boundaries::BC
-    postprocess::PP = nothing
+    postprocess::PP
 end
 Adapt.@adapt_structure Configuration
 
+Configuration(; schemes, solvers, runtime, hardware, boundaries, postprocess=nothing) = begin
+Configuration(
+    schemes,
+    solvers,
+    adapt(hardware.backend, runtime),
+    hardware,
+    boundaries,
+    postprocess
+)
+end
+
+
+
 
 """
     hardware = Hardware(backend, workgroup)

src/Solve/Solve_1_api.jl

@@ -3,6 +3,7 @@ export explicit_relaxation!, implicit_relaxation!, implicit_relaxation_diagdom!,
 export solve_system!
 export solve_equation!
 export residual!
+export AdaptiveTimeStepping
 
 struct SolverSetup{
     F<:AbstractFloat,
@@ -81,27 +82,68 @@ SolverSetup(;
             float_type(atol),
             float_type(rtol))
 
-struct Runtime{I<:Integer,F<:AbstractFloat}
+struct AdaptiveTimeStepping{F<:AbstractFloat}
+    maxCo::F
+    minShrink::F
+    maxGrow::F
+end
+Adapt.@adapt_structure AdaptiveTimeStepping
+
+"""
+    AdaptiveTimeStepping(; 
+        # keyword arguments
+
+        maxCo=0.75,
+        minShrink=0.1,
+        maxGrow=1.2
+    )
+
+Constructs an `AdaptiveTimeStepping` object used to control automatic time-step adjustment
+based on the Courant number.
+
+This struct is passed optionally to `Runtime` and enables adaptive time stepping in transient
+simulations. If not provided, a fixed time step is used.
+
+# Input arguments
+
+- `maxCo::AbstractFloat`: target maximum Courant number. The time step will be adjusted
+  such that the computed Courant number approaches this value.
+- `minShrink::AbstractFloat`: lower bound on the multiplicative factor applied to the
+  current time step. Prevents excessively large reductions in a single update.
+- `maxGrow::AbstractFloat`: upper bound on the multiplicative factor applied to the
+  current time step. Prevents excessive time-step growth.
+"""
+AdaptiveTimeStepping(;
+    maxCo=0.75,
+    minShrink=0.1,
+    maxGrow=1.2
+) = AdaptiveTimeStepping(float(maxCo), float(minShrink), float(maxGrow))
+
+struct Runtime{I<:Integer,F<:AbstractFloat, V<:AbstractVector{F}, A<:Union{Nothing, AdaptiveTimeStepping}}
     iterations::I
-    dt::F
+    dt::V
     write_interval::I
+    adaptive::A
 end
+Adapt.@adapt_structure Runtime
 
 """
     Runtime(; 
             # keyword arguments
 
             iterations::I, 
             write_interval::I, 
-            time_step::N
+            time_step::N,
+            adaptive::A
         ) where {I<:Integer,N<:Number} = begin
 
         # returned Runtime struct
         Runtime{I<:Integer,F<:AbstractFloat}
             (
                 iterations=iterations, 
                 dt=time_step, 
-                write_interval=write_interval
+                write_interval=write_interval,
+                adaptive=adaptive
             )
     end
 
@@ -112,15 +154,21 @@ This is a convenience function to set the top-level runtime information. The inp
 - `iterations::Integer`: specifies the number of iterations in a simulation run.
 - `write_interval::Integer`: defines how often simulation results are written to file (on the current working directory). The interval is currently based on number of iterations. Set to `-1` to run without writing results to file.
 - `time_step::AbstractFloat`: the time step to use in the simulation. Notice that for steady solvers this is simply a counter and it is recommended to simply use `1`.
+- `adaptive::Union{Nothing, AdaptiveTimeStepping}`: optionally enables adaptive time stepping. Pass an `AdaptiveTimeStepping` object to automatically adjust `dt` based on the Courant number during transient simulations. Defaults to `nothing`, meaning a fixed time step is used.
 
 # Example
 
 ```julia
 runtime = Runtime(iterations=2000, time_step=1, write_interval=2000)
 ```
 """
-Runtime(; iterations::I, write_interval::I, time_step::N) where {I<:Integer,N<:Number} = begin
-    Runtime(iterations, float(time_step), write_interval)
+Runtime(; iterations::I,
+          write_interval::I,
+          time_step::N,
+          adaptive=nothing) where {I<:Integer,N<:Number} = begin
+
+    val = float(time_step)
+    Runtime(iterations, [val], write_interval, adaptive)
 end
 
 # Set schemes function definition with default set variables

src/Solvers/Solvers_0_functions.jl

@@ -267,7 +267,7 @@ end
 @kernel function _max_courant_number!(cellsCourant, U, runtime, mesh::Mesh3)
     i = @index(Global)
     @uniform cells = mesh.cells
-    dt = runtime.dt
+    dt = runtime.dt[1]
     umag = norm(U[i])
     volume = cells[i].volume
     dx = volume^0.333333
@@ -277,9 +277,20 @@ end
 @kernel function _max_courant_number!(cellsCourant, U, runtime, mesh::Mesh2)
     i = @index(Global)
     @uniform cells = mesh.cells
-    dt = runtime.dt
+    dt = runtime.dt[1]
     umag = norm(U[i])
     volume = cells[i].volume
     dx = volume^0.5
     cellsCourant[i] = umag * dt / dx
+end
+
+
+update_dt!(runtime::Runtime{<:Any,<:Any,<:Any,Nothing}, courant) = nothing
+
+function update_dt!(runtime::Runtime{<:Any,<:Any,<:Any,<:AdaptiveTimeStepping}, courant)
+    (; maxCo, maxGrow, minShrink) = runtime.adaptive
+
+    courant_factor = maxCo / (courant + eps())
+    new_dt_factor = clamp(courant_factor, minShrink, maxGrow)
+    runtime.dt .= runtime.dt .* new_dt_factor
 end

src/Solvers/Solvers_1_CSIMPLE.jl

@@ -131,8 +131,11 @@ function CSIMPLE(
     (; iterations, write_interval,dt) = runtime
     (; backend) = hardware
 
+    dt_cpu = zeros(_get_float(mesh), 1)
+    copyto!(dt_cpu, config.runtime.dt)
+
     # rho = get_flux(U_eqn, 1)
-    postprocess = convert_time_to_iterations(postprocess,model,dt,iterations)
+    postprocess = convert_time_to_iterations(postprocess,model,dt_cpu[1],iterations)
     mdotf = get_flux(U_eqn, 2)
     mueff = get_flux(U_eqn, 3)
     mueffgradUt = get_source(U_eqn, 2)

src/Solvers/Solvers_1_LAPLACE.jl

@@ -125,8 +125,10 @@ function LAPLACE(
     (; iterations, write_interval, dt) = runtime
     (; backend) = hardware
 
+    dt_cpu = zeros(_get_float(mesh), 1)
+    copyto!(dt_cpu, config.runtime.dt)
 
-    postprocess = convert_time_to_iterations(postprocess,model,dt,iterations)
+    postprocess = convert_time_to_iterations(postprocess,model,dt_cpu[1],iterations)
     @info "Starting LAPLACE loops..."
     progress = Progress(iterations; dt=1.0, showspeed=true)
 
@@ -154,7 +156,7 @@ function LAPLACE(
 
         ProgressMeter.next!(
             progress, showvalues = [
-                (:time, iteration*runtime.dt),
+                (:time, iteration*dt_cpu[1]),
                 (:T_residual, R_T[iteration])
                 ]
             )

src/Solvers/Solvers_1_SIMPLE.jl

@@ -112,8 +112,11 @@ function SIMPLE(
     (; solvers, schemes, runtime, hardware, boundaries, postprocess) = config
     (; iterations, write_interval,dt) = runtime
     (; backend) = hardware
+
+    dt_cpu = zeros(_get_float(mesh), 1)
+    copyto!(dt_cpu, config.runtime.dt)
 
-    postprocess = convert_time_to_iterations(postprocess,model,dt,iterations)
+    postprocess = convert_time_to_iterations(postprocess,model,dt_cpu[1],iterations)
     mdotf = get_flux(U_eqn, 2)
     nueff = get_flux(U_eqn, 3)
     rDf = get_flux(p_eqn, 1)
@@ -150,7 +153,6 @@ function SIMPLE(
     grad!(∇p, pf, p, boundaries.p, time, config)
     limit_gradient!(schemes.p.limiter, ∇p, p, config)
 
-
     update_nueff!(nueff, nu, model.turbulence, config)
 
     @info "Starting SIMPLE loops..."