Add unstable_check for custom SSP method

examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl

@@ -102,4 +102,5 @@ stage_callbacks = (SubcellLimiterIDPCorrection(),
 
 sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks = stage_callbacks);
                   dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
-                  ode_default_options()..., callback = callbacks);
+                  ode_default_options()..., unstable_check = unstable_check,
+                  callback = callbacks);

src/Trixi.jl

@@ -307,7 +307,7 @@ export PositivityPreservingLimiterZhangShu, EntropyBoundedLimiter
 export trixi_include, examples_dir, get_examples, default_example,
        default_example_unstructured, ode_default_options
 
-export ode_norm, ode_unstable_check
+export ode_norm, ode_unstable_check, unstable_check
 
 export convergence_test, jacobian_fd, jacobian_ad_forward, linear_structure
 

src/time_integration/methods_2N.jl

@@ -110,7 +110,8 @@ mutable struct SimpleIntegrator2N{RealT <: Real, uType, Params, Sol, F, Alg,
 end
 
 function init(ode::ODEProblem, alg::SimpleAlgorithm2N;
-              dt, callback::Union{CallbackSet, Nothing} = nothing, kwargs...)
+              dt, callback::Union{CallbackSet, Nothing} = nothing,
+              unstable_check = ode_unstable_check, kwargs...)
     u = copy(ode.u0)
     du = similar(u)
     u_tmp = similar(u)

src/time_integration/methods_3Sstar.jl

@@ -149,7 +149,8 @@ mutable struct SimpleIntegrator3Sstar{RealT <: Real, uType, Params, Sol, F, Alg,
 end
 
 function init(ode::ODEProblem, alg::SimpleAlgorithm3Sstar;
-              dt, callback::Union{CallbackSet, Nothing} = nothing, kwargs...)
+              dt, callback::Union{CallbackSet, Nothing} = nothing,
+              unstable_check = ode_unstable_check, kwargs...)
     u = copy(ode.u0)
     du = similar(u)
     u_tmp1 = similar(u)

src/time_integration/methods_SSP.jl

@@ -52,25 +52,27 @@ struct SimpleSSPRK33{StageCallbacks} <: SimpleAlgorithmSSP
 end
 
 # This struct is needed to fake https://github.com/SciML/OrdinaryDiffEq.jl/blob/0c2048a502101647ac35faabd80da8a5645beac7/src/integrators/type.jl#L1
-mutable struct SimpleIntegratorSSPOptions{Callback, TStops}
+mutable struct SimpleIntegratorSSPOptions{Callback, UnstableCheck, TStops}
     callback::Callback # callbacks; used in Trixi
+    unstable_check::UnstableCheck # unstable check function
     adaptive::Bool # whether the algorithm is adaptive; ignored
     dtmax::Float64 # ignored
     maxiters::Int # maximal number of time steps
     tstops::TStops # tstops from https://diffeq.sciml.ai/v6.8/basics/common_solver_opts/#Output-Control-1; ignored
 end
 
-function SimpleIntegratorSSPOptions(callback, tspan; maxiters = typemax(Int), kwargs...)
+function SimpleIntegratorSSPOptions(callback, tspan; maxiters = typemax(Int),
+                                    unstable_check = unstable_check, kwargs...)
     tstops_internal = BinaryHeap{eltype(tspan)}(FasterForward())
     # We add last(tspan) to make sure that the time integration stops at the end time
     push!(tstops_internal, last(tspan))
     # We add 2 * last(tspan) because add_tstop!(integrator, t) is only called by DiffEqCallbacks.jl if tstops contains a time that is larger than t
     # (https://github.com/SciML/DiffEqCallbacks.jl/blob/025dfe99029bd0f30a2e027582744528eb92cd24/src/iterative_and_periodic.jl#L92)
     push!(tstops_internal, 2 * last(tspan))
-    SimpleIntegratorSSPOptions{typeof(callback), typeof(tstops_internal)}(callback,
-                                                                          false, Inf,
-                                                                          maxiters,
-                                                                          tstops_internal)
+    SimpleIntegratorSSPOptions{typeof(callback), typeof(unstable_check),
+                               typeof(tstops_internal)}(callback, unstable_check,
+                                                        false, Inf, maxiters,
+                                                        tstops_internal)
 end
 
 # This struct is needed to fake https://github.com/SciML/OrdinaryDiffEq.jl/blob/0c2048a502101647ac35faabd80da8a5645beac7/src/integrators/type.jl#L77
@@ -117,7 +119,8 @@ has_tstop(integrator::SimpleIntegratorSSP) = !isempty(integrator.opts.tstops)
 first_tstop(integrator::SimpleIntegratorSSP) = first(integrator.opts.tstops)
 
 function init(ode::ODEProblem, alg::SimpleAlgorithmSSP;
-              dt, callback::Union{CallbackSet, Nothing} = nothing, kwargs...)
+              dt, callback::Union{CallbackSet, Nothing} = nothing,
+              unstable_check = ode_unstable_check, kwargs...)
     u = copy(ode.u0)
     du = similar(u)
     u_tmp = similar(u)
@@ -127,6 +130,7 @@ function init(ode::ODEProblem, alg::SimpleAlgorithmSSP;
     integrator = SimpleIntegratorSSP(u, du, u_tmp, t, tdir, dt, dt, iter, ode.p,
                                      (prob = ode,), ode.f, alg,
                                      SimpleIntegratorSSPOptions(callback, ode.tspan;
+                                                                unstable_check = unstable_check,
                                                                 kwargs...),
                                      false, true, false)
 
@@ -156,6 +160,7 @@ function solve!(integrator::SimpleIntegratorSSP)
     @unpack alg = integrator
     t_end = last(prob.tspan)
     callbacks = integrator.opts.callback
+    (; unstable_check) = integrator.opts
 
     integrator.finalstep = false
     @trixi_timeit timer() "main loop" while !integrator.finalstep
@@ -206,6 +211,11 @@ function solve!(integrator::SimpleIntegratorSSP)
             end
         end
 
+        if unstable_check(integrator.dt, integrator.u, integrator, integrator.t)
+            @warn "Instability detected. Aborting"
+            terminate!(integrator)
+        end
+
         # respect maximum number of iterations
         if integrator.iter >= integrator.opts.maxiters && !integrator.finalstep
             @warn "Interrupted. Larger maxiters is needed."

src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl

@@ -115,10 +115,10 @@ The original paper is
 
 # Arguments
 - `num_stages` (`Int`): Number of stages in the PERK method.
-- `base_path_monomial_coeffs` (`AbstractString`): Path to a file containing 
+- `base_path_monomial_coeffs` (`AbstractString`): Path to a file containing
     monomial coefficients of the stability polynomial of PERK method.
     The coefficients should be stored in a text file at `joinpath(base_path_monomial_coeffs, "gamma_$(num_stages).txt")` and separated by line breaks.
-- `dt_opt` (`Float64`, optional): Optimal time step size for the simulation setup. Can be `nothing` if it is unknown. 
+- `dt_opt` (`Float64`, optional): Optimal time step size for the simulation setup. Can be `nothing` if it is unknown.
     In this case the optimal CFL number cannot be computed and the [`StepsizeCallback`](@ref) cannot be used.
 - `tspan`: Time span of the simulation.
 - `semi` (`AbstractSemidiscretization`): Semidiscretization setup.
@@ -132,7 +132,7 @@ The original paper is
 
 !!! note
     To use this integrator, the user must import the
-    [Convex.jl](https://github.com/jump-dev/Convex.jl) and 
+    [Convex.jl](https://github.com/jump-dev/Convex.jl) and
     [ECOS.jl](https://github.com/jump-dev/ECOS.jl) packages
     unless the coefficients are provided in a `gamma_<num_stages>.txt` file.
 """
@@ -214,7 +214,8 @@ mutable struct PairedExplicitRK2Integrator{RealT <: Real, uType, Params, Sol, F,
 end
 
 function init(ode::ODEProblem, alg::PairedExplicitRK2;
-              dt, callback::Union{CallbackSet, Nothing} = nothing, kwargs...)
+              dt, callback::Union{CallbackSet, Nothing} = nothing,
+              unstable_check = ode_unstable_check, kwargs...)
     u0 = copy(ode.u0)
     du = zero(u0)
     u_tmp = zero(u0)

src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl

@@ -111,17 +111,17 @@ The original paper is
   Third-order Paired Explicit Runge-Kutta schemes for stiff systems of equations
   [DOI: 10.1016/j.jcp.2022.111470](https://doi.org/10.1016/j.jcp.2022.111470)
 
-  While the changes to SSPRK33 base-scheme are described in 
+  While the changes to SSPRK33 base-scheme are described in
 - Doehring, Schlottke-Lakemper, Gassner, Torrilhon (2024)
   Multirate Time-Integration based on Dynamic ODE Partitioning through Adaptively Refined Meshes for Compressible Fluid Dynamics
-  [DOI: 10.1016/j.jcp.2024.113223](https://doi.org/10.1016/j.jcp.2024.113223)                      
+  [DOI: 10.1016/j.jcp.2024.113223](https://doi.org/10.1016/j.jcp.2024.113223)
 
 # Arguments
 - `num_stages` (`Int`): Number of stages in the paired explicit Runge-Kutta (PERK) method.
-- `base_path_a_coeffs` (`AbstractString`): Path to a file containing some coefficients in the A-matrix in 
+- `base_path_a_coeffs` (`AbstractString`): Path to a file containing some coefficients in the A-matrix in
     the Butcher tableau of the Runge Kutta method.
     The matrix should be stored in a text file at `joinpath(base_path_a_coeffs, "a_$(num_stages).txt")` and separated by line breaks.
-- `dt_opt` (`Float64`, optional): Optimal time step size for the simulation setup. Can be `nothing` if it is unknown. 
+- `dt_opt` (`Float64`, optional): Optimal time step size for the simulation setup. Can be `nothing` if it is unknown.
     In this case the optimal CFL number cannot be computed and the [`StepsizeCallback`](@ref) cannot be used.
 - `tspan`: Time span of the simulation.
 - `semi` (`AbstractSemidiscretization`): Semidiscretization setup.
@@ -132,8 +132,8 @@ The original paper is
     $S$ is the number of stages. Default is `1.0f0`.
 
 !!! note
-    To use this integrator, the user must import the 
-    [Convex.jl](https://github.com/jump-dev/Convex.jl), 
+    To use this integrator, the user must import the
+    [Convex.jl](https://github.com/jump-dev/Convex.jl),
     [ECOS.jl](https://github.com/jump-dev/ECOS.jl), and
     [NLSolve.jl](https://github.com/JuliaNLSolvers/NLsolve.jl) packages
     unless the A-matrix coefficients are provided in a `a_<num_stages>.txt` file.
@@ -208,7 +208,8 @@ mutable struct PairedExplicitRK3Integrator{RealT <: Real, uType, Params, Sol, F,
 end
 
 function init(ode::ODEProblem, alg::PairedExplicitRK3;
-              dt, callback::Union{CallbackSet, Nothing} = nothing, kwargs...)
+              dt, callback::Union{CallbackSet, Nothing} = nothing,
+              unstable_check = ode_unstable_check, kwargs...)
     u0 = copy(ode.u0)
     du = zero(u0)
     u_tmp = zero(u0)

src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl

@@ -108,18 +108,18 @@ end
 The following structures and methods provide an implementation of
 the fourth-order paired explicit Runge-Kutta (P-ERK) method
 optimized for a certain simulation setup (PDE, IC & BCs, Riemann Solver, DG Solver).
-The method has been proposed in 
+The method has been proposed in
 
 - D. Doehring, L. Christmann, M. Schlottke-Lakemper, G. J. Gassner and M. Torrilhon (2024).
   Fourth-Order Paired-Explicit Runge-Kutta Methods
   [DOI:10.48550/arXiv.2408.05470](https://doi.org/10.48550/arXiv.2408.05470)
 
 # Arguments
 - `num_stages` (`Int`): Number of stages in the paired explicit Runge-Kutta (P-ERK) method.
-- `base_path_a_coeffs` (`AbstractString`): Path to a file containing some coefficients in the A-matrix in 
+- `base_path_a_coeffs` (`AbstractString`): Path to a file containing some coefficients in the A-matrix in
     the Butcher tableau of the Runge Kutta method.
     The matrix should be stored in a text file at `joinpath(base_path_a_coeffs, "a_$(num_stages).txt")` and separated by line breaks.
-- `dt_opt` (`Float64`, optional): Optimal time step size for the simulation setup. Can be `nothing` if it is unknown. 
+- `dt_opt` (`Float64`, optional): Optimal time step size for the simulation setup. Can be `nothing` if it is unknown.
     In this case the optimal CFL number cannot be computed and the [`StepsizeCallback`](@ref) cannot be used.
 - `tspan`: Time span of the simulation.
 - `semi` (`AbstractSemidiscretization`): Semidiscretization setup.
@@ -130,16 +130,16 @@ The method has been proposed in
 
 !!! note
     To use this integrator, the user must import the
-    [Convex.jl](https://github.com/jump-dev/Convex.jl) and 
+    [Convex.jl](https://github.com/jump-dev/Convex.jl) and
     [ECOS.jl](https://github.com/jump-dev/ECOS.jl) packages
-    unless the coefficients are provided in a `a_<num_stages>.txt` file.   
+    unless the coefficients are provided in a `a_<num_stages>.txt` file.
 """
 struct PairedExplicitRK4 <: AbstractPairedExplicitRKSingle
     num_stages::Int # S
 
     # Optimized coefficients, i.e., flexible part of the Butcher array matrix A.
     a_matrix::Union{Matrix{Float64}, Nothing}
-    # This part of the Butcher array matrix A is constant for all PERK methods, i.e., 
+    # This part of the Butcher array matrix A is constant for all PERK methods, i.e.,
     # regardless of the optimized coefficients.
     a_matrix_constant::SMatrix{2, 3, Float64}
     c::Vector{Float64}
@@ -208,7 +208,8 @@ mutable struct PairedExplicitRK4Integrator{RealT <: Real, uType, Params, Sol, F,
 end
 
 function init(ode::ODEProblem, alg::PairedExplicitRK4;
-              dt, callback::Union{CallbackSet, Nothing} = nothing, kwargs...)
+              dt, callback::Union{CallbackSet, Nothing} = nothing,
+              unstable_check = ode_unstable_check, kwargs...)
     u0 = copy(ode.u0)
     du = zero(u0)
     u_tmp = zero(u0)

src/time_integration/time_integration.jl

@@ -36,20 +36,22 @@ function set_proposed_dt!(integrator::AbstractTimeIntegrator, dt)
     (integrator.dt = dt; integrator.dtcache = dt)
 end
 
-# Required e.g. for `glm_speed_callback` 
+# Required e.g. for `glm_speed_callback`
 function get_proposed_dt(integrator::AbstractTimeIntegrator)
     return integrator.dt
 end
 
 """
-    Trixi.solve(ode::ODEProblem, alg::AbstractTimeIntegrationAlgorithm; 
+    Trixi.solve(ode::ODEProblem, alg::AbstractTimeIntegrationAlgorithm;
                 dt, callbacks, kwargs...)
 
 Fakes `solve` from https://diffeq.sciml.ai/v6.8/basics/overview/#Solving-the-Problems-1
 """
 function solve(ode::ODEProblem, alg::AbstractTimeIntegrationAlgorithm;
-               dt, callback = nothing, kwargs...)
-    integrator = init(ode, alg, dt = dt, callback = callback; kwargs...)
+               dt, callback = nothing,
+               unstable_check = ode_unstable_check, kwargs...)
+    integrator = init(ode, alg, dt = dt, callback = callback,
+                      unstable_check = unstable_check; kwargs...)
 
     # Start actual solve
     solve!(integrator)
@@ -82,6 +84,15 @@ function DiffEqBase.get_tstops_max(integrator::AbstractTimeIntegrator)
     return maximum(get_tstops_array(integrator))
 end
 
+function unstable_check(dt, u_ode, integrator::AbstractTimeIntegrator, t)
+    if mpi_isparallel()
+        u_isfinite = MPI.Allreduce!(Ref(all(isfinite, u_ode)), Base.min, mpi_comm())[]
+    else
+        u_isfinite = all(isfinite, u_ode)
+    end
+    return !isfinite(dt) || !u_isfinite
+end
+
 function finalize_callbacks(integrator::AbstractTimeIntegrator)
     callbacks = integrator.opts.callback