Handles restarting problem of adaptive time integration methods

examples/tree_2d_dgsem/elixir_advection_extended.jl

@@ -54,7 +54,7 @@ analysis_callback = AnalysisCallback(semi, interval=analysis_interval,
 alive_callback = AliveCallback(analysis_interval=analysis_interval)
 
 # The SaveRestartCallback allows to save a file from which a Trixi.jl simulation can be restarted
-save_restart = SaveRestartCallback(interval=100,
+save_restart = SaveRestartCallback(interval=40,
                                    save_final_restart=true)
 
 # The SaveSolutionCallback allows to save the solution to a file in regular intervals
@@ -77,9 +77,10 @@ callbacks = CallbackSet(summary_callback,
 # run the simulation
 
 # OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
-sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
+alg = CarpenterKennedy2N54(williamson_condition=false)
+sol = solve(ode, alg,
             dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
-            save_everystep=false, callback=callbacks);
+            save_everystep=false, callback=callbacks; ode_default_options()...);
 
 # Print the timer summary
 summary_callback()

examples/tree_2d_dgsem/elixir_advection_restart.jl

@@ -3,9 +3,10 @@ using OrdinaryDiffEq
 using Trixi
 
 ###############################################################################
-# create a restart file
-
-trixi_include(@__MODULE__, joinpath(@__DIR__, "elixir_advection_extended.jl"))
+# Define time integration algorithm
+alg = CarpenterKennedy2N54(williamson_condition=false)
+# Create a restart file
+trixi_include(@__MODULE__, joinpath(@__DIR__, "elixir_advection_extended.jl"), alg = alg, tspan = (0.0, 10.0))
 
 
 ###############################################################################
@@ -14,22 +15,26 @@ trixi_include(@__MODULE__, joinpath(@__DIR__, "elixir_advection_extended.jl"))
 # Note: If you get a restart file from somewhere else, you need to provide
 # appropriate setups in the elixir loading a restart file
 
-restart_filename = joinpath("out", "restart_000018.h5")
+restart_filename = joinpath("out", "restart_000040.h5")
 mesh = load_mesh(restart_filename)
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
 
-tspan = (load_time(restart_filename), 2.0)
+tspan = (load_time(restart_filename), 10.0)
 dt = load_dt(restart_filename)
 ode = semidiscretize(semi, tspan, restart_filename);
 
 # Do not overwrite the initial snapshot written by elixir_advection_extended.jl.
 save_solution.condition.save_initial_solution = false
 
-alg = CarpenterKennedy2N54(williamson_condition=false)
 integrator = init(ode, alg,
                   dt=dt, # solve needs some value here but it will be overwritten by the stepsize_callback
-                  save_everystep=false, callback=callbacks)
+                  save_everystep=false, callback=callbacks; ode_default_options()...)
+
+# Load saved context for adaptive time integrator
+if integrator.opts.adaptive
+    load_adaptive_time_integrator!(integrator, restart_filename)
+end
 
 # Get the last time index and work with that.
 load_timestep!(integrator, restart_filename)

src/Trixi.jl

@@ -255,7 +255,8 @@ export SummaryCallback, SteadyStateCallback, AnalysisCallback, AliveCallback,
        GlmSpeedCallback, LBMCollisionCallback, EulerAcousticsCouplingCallback,
        TrivialCallback, AnalysisCallbackCoupled
 
-export load_mesh, load_time, load_timestep, load_timestep!, load_dt
+export load_mesh, load_time, load_timestep, load_timestep!, load_dt,
+       load_adaptive_time_integrator!
 
 export ControllerThreeLevel, ControllerThreeLevelCombined,
        IndicatorLöhner, IndicatorLoehner, IndicatorMax,

src/callbacks_step/save_restart.jl

@@ -105,6 +105,11 @@ function (restart_callback::SaveRestartCallback)(integrator)
         end
 
         save_restart_file(u_ode, t, dt, iter, semi, restart_callback)
+        # If using an adaptive time stepping scheme, store controller values for restart
+        if integrator.opts.adaptive
+            save_adaptive_time_integrator(integrator, integrator.opts.controller,
+                                          restart_callback)
+        end
     end
 
     # avoid re-evaluating possible FSAL stages
@@ -168,5 +173,36 @@ function load_restart_file(semi::AbstractSemidiscretization, restart_file)
     load_restart_file(mesh_equations_solver_cache(semi)..., restart_file)
 end
 
+"""
+    load_adaptive_time_integrator!(integrator, restart_file::AbstractString)
+
+Load the context information for time integrators with error-based step size control
+saved in a `restart_file`.
+"""
+function load_adaptive_time_integrator!(integrator, restart_file::AbstractString)
+    controller = integrator.opts.controller
+    # Read context information for controller
+    h5open(restart_file, "r") do file
+        # Ensure that the necessary information was saved
+        if !("time_integrator_qold" in keys(attributes(file))) ||
+           !("time_integrator_dtpropose" in keys(attributes(file))) ||
+           (hasproperty(controller, :err) &&
+            !("time_integrator_controller_err" in keys(attributes(file))))
+            error("Missing data in restart file: check the consistency of adaptive time controller with initial setup!")
+        end
+        # Load data that is required both for PIController and PIDController
+        integrator.qold = read(attributes(file)["time_integrator_qold"])
+        integrator.dtpropose = read(attributes(file)["time_integrator_dtpropose"])
+        # Accept step to use dtpropose already in the first step
+        integrator.accept_step = true
+        # Reevaluate integrator.fsal_first on the first step
+        integrator.reeval_fsal = true
+        # Load additional parameters for PIDController
+        if hasproperty(controller, :err) # Distinguish PIDController from PIController 
+            controller.err[:] = read(attributes(file)["time_integrator_controller_err"])
+        end
+    end
+end
+
 include("save_restart_dg.jl")
 end # @muladd

src/callbacks_step/save_restart_dg.jl

@@ -327,4 +327,28 @@ function load_restart_file_on_root(mesh::Union{ParallelTreeMesh, ParallelP4estMe
 
     return u_ode
 end
+
+# Store controller values for an adaptive time stepping scheme
+function save_adaptive_time_integrator(integrator,
+                                       controller, restart_callback)
+    # Save only on root
+    if mpi_isroot()
+        @unpack output_directory = restart_callback
+        timestep = integrator.stats.naccept
+
+        # Filename based on current time step
+        filename = joinpath(output_directory, @sprintf("restart_%06d.h5", timestep))
+
+        # Open file (preserve existing content)
+        h5open(filename, "r+") do file
+            # Add context information as attributes both for PIController and PIDController
+            attributes(file)["time_integrator_qold"] = integrator.qold
+            attributes(file)["time_integrator_dtpropose"] = integrator.dtpropose
+            # For PIDController is necessary to save additional parameters
+            if hasproperty(controller, :err) # Distinguish PIDController from PIController
+                attributes(file)["time_integrator_controller_err"] = controller.err
+            end
+        end
+    end
+end
 end # @muladd

test/test_mpi_tree.jl

@@ -23,10 +23,22 @@ CI_ON_WINDOWS = (get(ENV, "GITHUB_ACTIONS", false) == "true") && Sys.iswindows()
   end
 
   @trixi_testset "elixir_advection_restart.jl" begin
-    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart.jl"),
-      # Expected errors are exactly the same as in the serial test!
-      l2   = [7.81674284320524e-6],
-      linf = [6.314906965243505e-5])
+    using OrdinaryDiffEq: RDPK3SpFSAL49
+    Trixi.mpi_isroot() && println("═"^100)
+    Trixi.mpi_isroot() && println(joinpath(EXAMPLES_DIR, "elixir_advection_extended.jl"))
+    trixi_include(@__MODULE__, joinpath(EXAMPLES_DIR, "elixir_advection_extended.jl"),
+      alg = RDPK3SpFSAL49(), tspan = (0.0, 10.0))
+    l2_expected, linf_expected = analysis_callback(sol)
+
+    Trixi.mpi_isroot() && println("═"^100)
+    Trixi.mpi_isroot() && println(joinpath(EXAMPLES_DIR, "elixir_advection_restart.jl"))
+    # Errors are exactly the same as in the elixir_advection_extended.jl
+    trixi_include(@__MODULE__, joinpath(EXAMPLES_DIR, "elixir_advection_restart.jl"),
+      alg = RDPK3SpFSAL49())
+    l2_actual, linf_actual = analysis_callback(sol)
+
+    Trixi.mpi_isroot() && @test l2_actual == l2_expected
+    Trixi.mpi_isroot() && @test linf_actual == linf_expected
   end
 
   @trixi_testset "elixir_advection_mortar.jl" begin

test/test_threaded.jl

@@ -12,27 +12,38 @@ Trixi.mpi_isroot() && isdir(outdir) && rm(outdir, recursive=true)
 @testset "Threaded tests" begin
   @testset "TreeMesh" begin
     @trixi_testset "elixir_advection_restart.jl" begin
-      @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_advection_restart.jl"),
-        # Expected errors are exactly the same as in the serial test!
-        l2   = [7.81674284320524e-6],
-        linf = [6.314906965243505e-5])
+      elixir = joinpath(examples_dir(), "tree_2d_dgsem", "elixir_advection_extended.jl")
+      Trixi.mpi_isroot() && println("═"^100)
+      Trixi.mpi_isroot() && println(elixir)
+      trixi_include(@__MODULE__, elixir, tspan = (0.0, 10.0))
+      l2_expected, linf_expected = analysis_callback(sol)
+
+      elixir = joinpath(examples_dir(), "tree_2d_dgsem", "elixir_advection_restart.jl")
+      Trixi.mpi_isroot() && println("═"^100)
+      Trixi.mpi_isroot() && println(elixir)
+      # Errors are exactly the same as in the elixir_advection_extended.jl
+      trixi_include(@__MODULE__, elixir)
+      l2_actual, linf_actual = analysis_callback(sol)
+
+      Trixi.mpi_isroot() && @test l2_actual == l2_expected
+      Trixi.mpi_isroot() && @test linf_actual == linf_expected
 
-        # 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)) < 5000
-        end
+      # 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)) < 5000
+      end
     end
 
     @trixi_testset "elixir_advection_restart.jl with threaded time integration" begin
       @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_advection_restart.jl"),
         alg = CarpenterKennedy2N54(williamson_condition = false, thread = OrdinaryDiffEq.True()),
         # Expected errors are exactly the same as in the serial test!
-        l2   = [7.81674284320524e-6],
-        linf = [6.314906965243505e-5])
+        l2   = [8.005068880114254e-6],
+        linf = [6.39093577996519e-5])
     end
 
     @trixi_testset "elixir_advection_amr_refine_twice.jl" begin

test/test_tree_2d_advection.jl

@@ -25,10 +25,22 @@ EXAMPLES_DIR = pkgdir(Trixi, "examples", "tree_2d_dgsem")
   end
 
   @trixi_testset "elixir_advection_restart.jl" begin
-    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart.jl"),
-      # Expected errors are exactly the same as in the parallel test!
-      l2   = [7.81674284320524e-6],
-      linf = [6.314906965243505e-5])
+    using OrdinaryDiffEq: SSPRK43
+    println("═"^100)
+    println(joinpath(EXAMPLES_DIR, "elixir_advection_extended.jl"))
+    trixi_include(@__MODULE__, joinpath(EXAMPLES_DIR, "elixir_advection_extended.jl"),
+      alg = SSPRK43(), tspan = (0.0, 10.0))
+    l2_expected, linf_expected = analysis_callback(sol)
+
+    println("═"^100)
+    println(joinpath(EXAMPLES_DIR, "elixir_advection_restart.jl"))
+    # Errors are exactly the same as in the elixir_advection_extended.jl
+    trixi_include(@__MODULE__, joinpath(EXAMPLES_DIR, "elixir_advection_restart.jl"),
+      alg = SSPRK43())
+    l2_actual, linf_actual = analysis_callback(sol)
+
+    @test l2_actual == l2_expected
+    @test linf_actual == linf_expected
   end
 
   @trixi_testset "elixir_advection_mortar.jl" begin