diff --git a/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl b/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl index d63e3cd06ea..875c3eae434 100644 --- a/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl +++ b/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); diff --git a/src/Trixi.jl b/src/Trixi.jl index 6f31873a721..40ad73c4578 100644 --- a/src/Trixi.jl +++ b/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 diff --git a/src/time_integration/methods_2N.jl b/src/time_integration/methods_2N.jl index d626163d1a0..c38d754c0d1 100644 --- a/src/time_integration/methods_2N.jl +++ b/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) diff --git a/src/time_integration/methods_3Sstar.jl b/src/time_integration/methods_3Sstar.jl index 7862a273544..4f87b0e6f26 100644 --- a/src/time_integration/methods_3Sstar.jl +++ b/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) diff --git a/src/time_integration/methods_SSP.jl b/src/time_integration/methods_SSP.jl index 3286178d1c4..0dadb6f1522 100644 --- a/src/time_integration/methods_SSP.jl +++ b/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." diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl index 3f9460dd433..4e087b3fc60 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK2.jl +++ b/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_.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) diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl index 4f7344678d9..dde29a3ed37 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK3.jl +++ b/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_.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) diff --git a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl index d0571b9b0ae..abbd4b006c8 100644 --- a/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl +++ b/src/time_integration/paired_explicit_runge_kutta/methods_PERK4.jl @@ -108,7 +108,7 @@ 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 @@ -116,10 +116,10 @@ The method has been proposed in # 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_.txt` file. + unless the coefficients are provided in a `a_.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) diff --git a/src/time_integration/time_integration.jl b/src/time_integration/time_integration.jl index 9e3d802935b..f2c005cffb1 100644 --- a/src/time_integration/time_integration.jl +++ b/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