Enable save solution with time intervals for SimpleIntegratorSSP

Project.toml

@@ -42,6 +42,8 @@ TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 Triangulate = "f7e6ffb2-c36d-4f8f-a77e-16e897189344"
 TriplotBase = "981d1d27-644d-49a2-9326-4793e63143c3"
 TriplotRecipes = "808ab39a-a642-4abf-81ff-4cb34ebbffa3"
+DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 
 [weakdeps]
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"

examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl

@@ -68,7 +68,7 @@ analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
 
 alive_callback = AliveCallback(analysis_interval=analysis_interval)
 
-save_solution = SaveSolutionCallback(interval=100,
+save_solution = SaveSolutionCallback(dt=0.1,
                                      save_initial_solution=true,
                                      save_final_solution=true,
                                      solution_variables=cons2prim)

src/Trixi.jl

@@ -37,7 +37,8 @@ using SciMLBase: CallbackSet, DiscreteCallback,
 import SciMLBase: get_du, get_tmp_cache, u_modified!,
                   AbstractODEIntegrator, init, step!, check_error,
                   get_proposed_dt, set_proposed_dt!,
-                  terminate!, remake
+                  terminate!, remake, add_tstop!, has_tstop, first_tstop
+using OrdinaryDiffEq: modify_dt_for_tstops!
 using CodeTracking: CodeTracking
 using ConstructionBase: ConstructionBase
 using DiffEqCallbacks: PeriodicCallback, PeriodicCallbackAffect
@@ -70,6 +71,7 @@ using TriplotBase: TriplotBase
 using TriplotRecipes: DGTriPseudocolor
 @reexport using SimpleUnPack: @unpack
 using SimpleUnPack: @pack!
+using DataStructures: BinaryHeap, FasterForward
 
 # finite difference SBP operators
 using SummationByPartsOperators: AbstractDerivativeOperator,

src/time_integration/methods_SSP.jl

@@ -50,17 +50,22 @@ 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}
+mutable struct SimpleIntegratorSSPOptions{Callback, TStops}
     callback::Callback # callbacks; used in Trixi
     adaptive::Bool # whether the algorithm is adaptive; ignored
     dtmax::Float64 # ignored
     maxiters::Int # maximal number of time steps
-    tstops::Vector{Float64} # tstops from https://diffeq.sciml.ai/v6.8/basics/common_solver_opts/#Output-Control-1; ignored
+    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...)
-    SimpleIntegratorSSPOptions{typeof(callback)}(callback, false, Inf, maxiters,
-                                                 [last(tspan)])
+    tstops_internal = BinaryHeap{eltype(tspan)}(FasterForward())
+    push!(tstops_internal, last(tspan))
+    push!(tstops_internal, 2 * last(tspan))
+    SimpleIntegratorSSPOptions{typeof(callback), typeof(tstops_internal)}(callback,
+                                                                          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
@@ -73,6 +78,7 @@ mutable struct SimpleIntegratorSSP{RealT <: Real, uType, Params, Sol, F, Alg,
     du::uType
     r0::uType
     t::RealT
+    tdir::RealT
     dt::RealT # current time step
     dtcache::RealT # ignored
     iter::Int # current number of time steps (iteration)
@@ -82,8 +88,26 @@ mutable struct SimpleIntegratorSSP{RealT <: Real, uType, Params, Sol, F, Alg,
     alg::Alg
     opts::SimpleIntegratorSSPOptions
     finalstep::Bool # added for convenience
+    dtchangeable::Bool
+    force_stepfail::Bool
 end
 
+"""
+    add_tstop!(integrator::SimpleIntegratorSSP, t)
+Add
+"""
+function add_tstop!(integrator::SimpleIntegratorSSP, t)
+    integrator.tdir * (t - integrator.t) < zero(integrator.t) &&
+        error("Tried to add a tstop that is behind the current time. This is strictly forbidden")
+    if length(integrator.opts.tstops) > 1
+        pop!(integrator.opts.tstops)
+    end
+    push!(integrator.opts.tstops, integrator.tdir * t)
+end
+
+has_tstop(integrator::SimpleIntegratorSSP) = !isempty(integrator.opts.tstops)
+first_tstop(integrator::SimpleIntegratorSSP) = first(integrator.opts.tstops)
+
 # Forward integrator.stats.naccept to integrator.iter (see GitHub PR#771)
 function Base.getproperty(integrator::SimpleIntegratorSSP, field::Symbol)
     if field === :stats
@@ -108,11 +132,13 @@ function solve(ode::ODEProblem, alg = SimpleSSPRK33()::SimpleAlgorithmSSP;
     du = similar(u)
     r0 = similar(u)
     t = first(ode.tspan)
+    tdir = sign(ode.tspan[end] - ode.tspan[1])
     iter = 0
-    integrator = SimpleIntegratorSSP(u, du, r0, t, dt, zero(dt), iter, ode.p,
+    integrator = SimpleIntegratorSSP(u, du, r0, t, tdir, dt, zero(dt), iter, ode.p,
                                      (prob = ode,), ode.f, alg,
                                      SimpleIntegratorSSPOptions(callback, ode.tspan;
-                                                                kwargs...), false)
+                                                                kwargs...),
+                                     false, true, false)
 
     # resize container
     resize!(integrator.p, nelements(integrator.p.solver, integrator.p.cache))
@@ -155,6 +181,8 @@ function solve!(integrator::SimpleIntegratorSSP)
             terminate!(integrator)
         end
 
+        modify_dt_for_tstops!(integrator)
+
         @. integrator.r0 = integrator.u
         for stage in eachindex(alg.c)
             t_stage = integrator.t + integrator.dt * alg.c[stage]
@@ -216,7 +244,7 @@ end
 # stop the time integration
 function terminate!(integrator::SimpleIntegratorSSP)
     integrator.finalstep = true
-    empty!(integrator.opts.tstops)
+    #empty!(integrator.opts.tstops)
 end
 
 # used for AMR

test/test_tree_2d_euler.jl

@@ -65,8 +65,8 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem")
 
   @trixi_testset "elixir_euler_shockcapturing_subcell.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_shockcapturing_subcell.jl"),
-      l2   = [0.08508147906199143, 0.04510299017724501, 0.045103019801950375, 0.6930704343869766],
-      linf = [0.31123546471463326, 0.5616274869594462, 0.5619692712224448, 2.88670199345138])
+      l2   = [0.08508152653623638, 0.04510301725066843, 0.04510304668512745, 0.6930705064715306],
+      linf = [0.31136518019691406, 0.5617651935473419, 0.5621200790240503, 2.8866869108596056])
   end
 
   @trixi_testset "elixir_euler_blast_wave.jl" begin