Steady state callback

src/TrixiParticles.jl

@@ -20,7 +20,7 @@ using Polyester: Polyester, @batch
 using Printf: @printf, @sprintf
 using RecipesBase: RecipesBase, @series
 using SciMLBase: CallbackSet, DiscreteCallback, DynamicalODEProblem, u_modified!,
- get_tmp_cache, set_proposed_dt!, ODESolution, ODEProblem
+ get_tmp_cache, set_proposed_dt!, ODESolution, ODEProblem, terminate!
 @reexport using StaticArrays: SVector
 using StaticArrays: @SMatrix, SMatrix, setindex
 using StrideArrays: PtrArray, StaticInt
@@ -57,7 +57,7 @@ export WeaklyCompressibleSPHSystem, EntropicallyDampedSPHSystem, TotalLagrangian
  BoundarySPHSystem, DEMSystem, BoundaryDEMSystem, OpenBoundarySPHSystem, InFlow,
  OutFlow
 export InfoCallback, SolutionSavingCallback, DensityReinitializationCallback,
- PostprocessCallback, StepsizeCallback, UpdateCallback
+ PostprocessCallback, StepsizeCallback, UpdateCallback, SteadyStateCallback
 export ContinuityDensity, SummationDensity
 export PenaltyForceGanzenmueller
 export SchoenbergCubicSplineKernel, SchoenbergQuarticSplineKernel,

src/callbacks/callbacks.jl

@@ -31,3 +31,4 @@ include("density_reinit.jl")
 include("post_process.jl")
 include("stepsize.jl")
 include("update.jl")
+include("steady_state.jl")

src/callbacks/info.jl

@@ -57,18 +57,7 @@ end
 # affect!
 function (info_callback::InfoCallback)(integrator)
  if isfinished(integrator)
- println("─"^100)
- println("Trixi simulation finished. Final time: ", integrator.t,
- " Time steps: ", integrator.stats.naccept, " (accepted), ",
- integrator.iter, " (total)")
- println("─"^100)
- println()
-
- # Print timer
- TimerOutputs.complement!(timer())
- print_timer(timer(), title="TrixiParticles.jl",
- allocations=true, linechars=:unicode, compact=false)
- println()
+ print_summary(integrator)
  else
  t = integrator.t
  t_initial = first(integrator.sol.prob.tspan)
@@ -266,3 +255,18 @@ function summary_footer(io; total_width=100, indentation_level=0)
 
  print(io, s)
 end
+
+function print_summary(integrator)
+ println("─"^100)
+ println("Trixi simulation finished. Final time: ", integrator.t,
+ " Time steps: ", integrator.stats.naccept, " (accepted), ",
+ integrator.iter, " (total)")
+ println("─"^100)
+ println()
+
+ # Print timer
+ TimerOutputs.complement!(timer())
+ print_timer(timer(), title="TrixiParticles.jl",
+ allocations=true, linechars=:unicode, compact=false)
+ println()
+end

src/callbacks/steady_state.jl

@@ -0,0 +1,175 @@
+"""
+ SteadyStateCallback(; interval::Integer=0, dt=0.0, interval_size::Integer=10,
+ abstol=1.0e-8, reltol=1.0e-6)
+
+Terminates the integration when the residual of the change in kinetic energy
+falls below the threshold specified by `abstol, reltol`.
+
+# Keywords
+- `interval=0`: Check steady state condition every `interval` time steps.
+- `dt=0.0`: Check steady state condition in regular intervals of `dt` in terms
+ of integration time by adding additional `tstops`.
+- `interval_size`: The interval in which the change of the kinetic energy is considered.
+ `interval_size` is a (integer) multiple of `interval` or `dt`.
+- `abstol`: Absolute tolerance.
+- `reltol`: Relative tolerance.
+"""
+mutable struct SteadyStateCallback{I, RealT <: Real}
+ interval :: I
+ abstol :: RealT
+ reltol :: RealT
+ previous_ekin :: Vector{Float64}
+ interval_size :: Int
+end
+
+function SteadyStateCallback(; interval::Integer=0, dt=0.0, interval_size::Integer=10,
+ abstol=1.0e-8, reltol=1.0e-6)
+ abstol, reltol = promote(abstol, reltol)
+
+ if dt > 0 && interval > 0
+ throw(ArgumentError("setting both `interval` and `dt` is not supported"))
+ end
+
+ if dt > 0
+ interval = Float64(dt)
+ end
+
+ steady_state_callback = SteadyStateCallback(interval, abstol, reltol, [Inf64],
+ interval_size)
+
+ if dt > 0
+ return PeriodicCallback(steady_state_callback, dt, save_positions=(false, false),
+ final_affect=true)
+ else
+ return DiscreteCallback(steady_state_callback, steady_state_callback,
+ save_positions=(false, false))
+ end
+end
+
+function Base.show(io::IO, cb::DiscreteCallback{<:Any, <:SteadyStateCallback})
+ @nospecialize cb # reduce precompilation time
+
+ cb_ = cb.affect!
+
+ print(io, "SteadyStateCallback(abstol=", cb_.abstol, ", ", "reltol=", cb_.reltol, ")")
+end
+
+function Base.show(io::IO,
+ cb::DiscreteCallback{<:Any,
+ <:PeriodicCallbackAffect{<:SteadyStateCallback}})
+ @nospecialize cb # reduce precompilation time
+
+ cb_ = cb.affect!.affect!
+
+ print(io, "SteadyStateCallback(abstol=", cb_.abstol, ", reltol=", cb_.reltol, ")")
+end
+
+function Base.show(io::IO, ::MIME"text/plain",
+ cb::DiscreteCallback{<:Any, <:SteadyStateCallback})
+ @nospecialize cb # reduce precompilation time
+
+ if get(io, :compact, false)
+ show(io, cb)
+ else
+ cb_ = cb.affect!
+
+ setup = ["absolute tolerance" => cb_.abstol,
+ "relative tolerance" => cb_.reltol,
+ "interval" => cb_.interval,
+ "interval size" => cb_.interval_size]
+ summary_box(io, "SteadyStateCallback", setup)
+ end
+end
+
+function Base.show(io::IO, ::MIME"text/plain",
+ cb::DiscreteCallback{<:Any,
+ <:PeriodicCallbackAffect{<:SteadyStateCallback}})
+ @nospecialize cb # reduce precompilation time
+
+ if get(io, :compact, false)
+ show(io, cb)
+ else
+ cb_ = cb.affect!.affect!
+
+ setup = ["absolute tolerance" => cb_.abstol,
+ "relative tolerance" => cb_.reltol,
+ "interval" => cb_.interval,
+ "interval_size" => cb_.interval_size]
+ summary_box(io, "SteadyStateCallback", setup)
+ end
+end
+
+# `affect!` (`PeriodicCallback`)
+function (cb::SteadyStateCallback)(integrator)
+ steady_state_condition(cb, integrator) || return nothing
+
+ print_summary(integrator)
+
+ terminate!(integrator)
+end
+
+# `affect!` (`DiscreteCallback`)
+function (cb::SteadyStateCallback{Int})(integrator)
+ print_summary(integrator)
+
+ terminate!(integrator)
+end
+
+# `condition` (`DiscreteCallback`)
+function (steady_state_callback::SteadyStateCallback)(vu_ode, t, integrator)
+ return steady_state_condition(steady_state_callback, integrator)
+end
+
+@inline function steady_state_condition(cb, integrator)
+ (; abstol, reltol, previous_ekin, interval_size) = cb
+
+ vu_ode = integrator.u
+ v_ode, u_ode = vu_ode.x
+ semi = integrator.p
+
+ # Calculate kinetic energy
+ ekin = 0.0
+ foreach_system(semi) do system
+ v = wrap_v(v_ode, system, semi)
+
+ # Calculate kintetic energy
+ for particle in each_moving_particle(system)
+ velocity = current_velocity(v, system, particle)
+ ekin += 0.5 * system.mass[particle] * dot(velocity, velocity)
+ end
+ end
+
+ push!(previous_ekin, ekin)
+
+ terminate = false
+
+ if interval_size_condition(cb, integrator)
+ popfirst!(previous_ekin)
+
+ # Calculate MSE only over the `interval_size`
+ mse = 0.0
+ for index in 1:interval_size
+ mse += (previous_ekin[index] - ekin)^2
+ end
+
+ mse /= interval_size
+
+ threshold = abstol + reltol * ekin
+
+ terminate = mse <= threshold
+ end
+
+ return terminate
+end
+
+# `DiscreteCallback`
+@inline function interval_size_condition(cb::SteadyStateCallback{Int}, integrator)
+ return integrator.stats.naccept > 0 &&
+ round(integrator.stats.naccept / cb.interval) > cb.interval_size
+end
+
+# `PeriodicCallback`
+@inline function interval_size_condition(cb::SteadyStateCallback, integrator)
+ return integrator.stats.naccept > 0 &&
+ round(Int, integrator.t / cb.interval) > cb.interval_size
+end