From 4fe2824c94e6364bb05bd2e86071833ddbee493f Mon Sep 17 00:00:00 2001
From: Benjamin Bolm <74359358+bennibolm@users.noreply.github.com>
Date: Mon, 6 May 2024 09:35:08 +0200
Subject: [PATCH] Fix intern time integration method with fixed time step
 (#1919)

* Take minimum with current time step

* Fix bug using `dtcache`

* Reset not needed change

* Adjust comment on `dtcache` [skip ci]

* Modify time step before checking for end time

* Add test for SSP method with fixed time step and tstops

* fmt

* Add alive callback

* Move solution callback into test

* Adapt fixed step size to exactly hit 0.1
---
 src/time_integration/methods_SSP.jl | 16 +++++++--------
 test/test_tree_2d_euler.jl          | 32 +++++++++++++++++++++++++++++
 2 files changed, 40 insertions(+), 8 deletions(-)

diff --git a/src/time_integration/methods_SSP.jl b/src/time_integration/methods_SSP.jl
index 9d1e06488b4..18f9c613249 100644
--- a/src/time_integration/methods_SSP.jl
+++ b/src/time_integration/methods_SSP.jl
@@ -88,7 +88,7 @@ mutable struct SimpleIntegratorSSP{RealT <: Real, uType, Params, Sol, F, Alg,
     t::RealT
     tdir::RealT
     dt::RealT # current time step
-    dtcache::RealT # ignored
+    dtcache::RealT # manually set time step
     iter::Int # current number of time steps (iteration)
     p::Params # will be the semidiscretization from Trixi
     sol::Sol # faked
@@ -102,7 +102,7 @@ end
 
 """
     add_tstop!(integrator::SimpleIntegratorSSP, t)
-Add a time stop during the time integration process. 
+Add a time stop during the time integration process.
 This function is called after the periodic SaveSolutionCallback to specify the next stop to save the solution.
 """
 function add_tstop!(integrator::SimpleIntegratorSSP, t)
@@ -145,7 +145,7 @@ function solve(ode::ODEProblem, alg = SimpleSSPRK33()::SimpleAlgorithmSSP;
     t = first(ode.tspan)
     tdir = sign(ode.tspan[end] - ode.tspan[1])
     iter = 0
-    integrator = SimpleIntegratorSSP(u, du, r0, t, tdir, dt, zero(dt), iter, ode.p,
+    integrator = SimpleIntegratorSSP(u, du, r0, t, tdir, dt, dt, iter, ode.p,
                                      (prob = ode,), ode.f, alg,
                                      SimpleIntegratorSSPOptions(callback, ode.tspan;
                                                                 kwargs...),
@@ -185,6 +185,8 @@ function solve!(integrator::SimpleIntegratorSSP)
             error("time step size `dt` is NaN")
         end
 
+        modify_dt_for_tstops!(integrator)
+
         # if the next iteration would push the simulation beyond the end time, set dt accordingly
         if integrator.t + integrator.dt > t_end ||
            isapprox(integrator.t + integrator.dt, t_end)
@@ -192,8 +194,6 @@ 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]
@@ -232,7 +232,7 @@ function solve!(integrator::SimpleIntegratorSSP)
         end
     end
 
-    # Empty the tstops array. 
+    # Empty the tstops array.
     # This cannot be done in terminate!(integrator::SimpleIntegratorSSP) because DiffEqCallbacks.PeriodicCallbackAffect would return at error.
     extract_all!(integrator.opts.tstops)
 
@@ -253,12 +253,12 @@ u_modified!(integrator::SimpleIntegratorSSP, ::Bool) = false
 
 # used by adaptive timestepping algorithms in DiffEq
 function set_proposed_dt!(integrator::SimpleIntegratorSSP, dt)
-    integrator.dt = dt
+    (integrator.dt = dt; integrator.dtcache = dt)
 end
 
 # used by adaptive timestepping algorithms in DiffEq
 function get_proposed_dt(integrator::SimpleIntegratorSSP)
-    return integrator.dt
+    return ifelse(integrator.opts.adaptive, integrator.dt, integrator.dtcache)
 end
 
 # stop the time integration
diff --git a/test/test_tree_2d_euler.jl b/test/test_tree_2d_euler.jl
index 97a549ca4b8..d593dc540f7 100644
--- a/test/test_tree_2d_euler.jl
+++ b/test/test_tree_2d_euler.jl
@@ -235,6 +235,38 @@ end
     end
 end
 
+@trixi_testset "elixir_euler_shockcapturing_subcell.jl (fixed time step)" begin
+    # Testing local SSP method without stepsize callback
+    # Additionally, tests combination with SaveSolutionCallback using time interval
+    @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                 "elixir_euler_shockcapturing_subcell.jl"),
+                        dt=2.0e-3,
+                        tspan=(0.0, 0.25),
+                        save_solution=SaveSolutionCallback(dt = 0.1 + 1.0e-8),
+                        callbacks=CallbackSet(summary_callback, save_solution,
+                                              analysis_callback, alive_callback),
+                        l2=[
+                            0.05624855363458103,
+                            0.06931288786158463,
+                            0.06931283188960778,
+                            0.6200535829842072,
+                        ],
+                        linf=[
+                            0.29029967648805566,
+                            0.6494728865862608,
+                            0.6494729363533714,
+                            3.0949621505674787,
+                        ])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+        t = sol.t[end]
+        u_ode = sol.u[end]
+        du_ode = similar(u_ode)
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 15000
+    end
+end
+
 @trixi_testset "elixir_euler_blast_wave.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blast_wave.jl"),
                         l2=[