Restart HyperbolicParabolic with SplitODEProblem

examples/tree_1d_dgsem/elixir_advection_diffusion.jl

@@ -75,8 +75,12 @@ analysis_callback = AnalysisCallback(semi, interval = 100)
 # The AliveCallback prints short status information in regular intervals
 alive_callback = AliveCallback(analysis_interval = 100)
 
+# The SaveRestartCallback allows to save a file from which a Trixi.jl simulation can be restarted
+save_restart = SaveRestartCallback(interval = 100,
+                                   save_final_restart = true)
+
 # Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
-callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback)
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_restart)
 
 ###############################################################################
 # run the simulation

examples/tree_1d_dgsem/elixir_advection_diffusion_restart.jl

@@ -0,0 +1,30 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# create a restart file
+
+elixir_file = "elixir_advection_diffusion.jl"
+trixi_include(@__MODULE__, joinpath(@__DIR__, elixir_file))
+
+###############################################################################
+# initialize the ODE
+
+restart_file = "restart_000000018.h5"
+restart_filename = joinpath("out", restart_file)
+tspan = (load_time(restart_filename), 2.0)
+
+ode = semidiscretize(semi, tspan, restart_filename);
+
+# Do not save restart files here
+callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback)
+
+###############################################################################
+# run the simulation
+
+sol = solve(ode, KenCarp4(autodiff = false), abstol = time_abs_tol, reltol = time_int_tol,
+            save_everystep = false, callback = callbacks)
+
+# Print the timer summary
+summary_callback()

src/semidiscretization/semidiscretization.jl

@@ -102,7 +102,8 @@ function semidiscretize(semi::AbstractSemidiscretization, tspan;
 end
 
 """
-    semidiscretize(semi::AbstractSemidiscretization, tspan, restart_file::AbstractString)
+    semidiscretize(semi::AbstractSemidiscretization, tspan, 
+                   restart_file::AbstractString)
 
 Wrap the semidiscretization `semi` as an ODE problem in the time interval `tspan`
 that can be passed to `solve` from the [SciML ecosystem](https://diffeq.sciml.ai/latest/).

src/semidiscretization/semidiscretization_hyperbolic_parabolic.jl

@@ -295,6 +295,39 @@ function semidiscretize(semi::SemidiscretizationHyperbolicParabolic, tspan;
     return SplitODEProblem{iip}(rhs_parabolic!, rhs!, u0_ode, tspan, semi)
 end
 
+"""
+    semidiscretize(semi::SemidiscretizationHyperbolicParabolic, tspan,
+                   restart_file::AbstractString)
+
+Wrap the semidiscretization `semi` as a split ODE problem in the time interval `tspan`
+that can be passed to `solve` from the [SciML ecosystem](https://diffeq.sciml.ai/latest/).
+The parabolic right-hand side is the first function of the split ODE problem and
+will be used by default by the implicit part of IMEX methods from the
+SciML ecosystem.
+
+The initial condition etc. is taken from the `restart_file`.
+"""
+function semidiscretize(semi::SemidiscretizationHyperbolicParabolic, tspan,
+                        restart_file::AbstractString;
+                        reset_threads = true)
+    # Optionally reset Polyester.jl threads. See
+    # https://github.com/trixi-framework/Trixi.jl/issues/1583
+    # https://github.com/JuliaSIMD/Polyester.jl/issues/30
+    if reset_threads
+        Polyester.reset_threads!()
+    end
+
+    u0_ode = load_restart_file(semi, restart_file)
+    # TODO: MPI, do we want to synchronize loading and print debug statements, e.g. using
+    #       mpi_isparallel() && MPI.Barrier(mpi_comm())
+    #       See https://github.com/trixi-framework/Trixi.jl/issues/328
+    iip = true # is-inplace, i.e., we modify a vector when calling rhs_parabolic!, rhs!
+    # Note that the IMEX time integration methods of OrdinaryDiffEq.jl treat the
+    # first function implicitly and the second one explicitly. Thus, we pass the
+    # stiffer parabolic function first.
+    return SplitODEProblem{iip}(rhs_parabolic!, rhs!, u0_ode, tspan, semi)
+end
+
 function rhs!(du_ode, u_ode, semi::SemidiscretizationHyperbolicParabolic, t)
     @unpack mesh, equations, boundary_conditions, source_terms, solver, cache = semi
 

test/test_parabolic_1d.jl

@@ -28,6 +28,21 @@ isdir(outdir) && rm(outdir, recursive = true)
     end
 end
 
+@trixi_testset "TreeMesh1D: elixir_advection_diffusion_restart.jl" begin
+    @test_trixi_include(joinpath(examples_dir(), "tree_1d_dgsem",
+                                 "elixir_advection_diffusion_restart.jl"),
+                        l2=[1.0671615777620987e-5],
+                        linf=[3.861509422325993e-5])
+    # 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)) < 1000
+    end
+end
+
 @trixi_testset "TreeMesh1D: elixir_advection_diffusion.jl (AMR)" begin
     @test_trixi_include(joinpath(examples_dir(), "tree_1d_dgsem",
                                  "elixir_advection_diffusion.jl"),