Save and load user_data for p4est

examples/p4est_2d_dgsem/elixir_advection_restart_amr.jl

@@ -0,0 +1,54 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# create a restart file
+
+elixir_file = "elixir_advection_extended.jl"
+restart_file = "restart_000021.h5"
+
+trixi_include(@__MODULE__, joinpath(@__DIR__, elixir_file))
+
+###############################################################################
+# adapt the parameters that have changed compared to "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_file)
+mesh = load_mesh(restart_filename)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_conditions)
+
+tspan = (load_time(restart_filename), 2.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
+
+# Add AMR callback
+amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable = first),
+                                      base_level = 0,
+                                      med_level = 0, med_threshold = 0.8,
+                                      max_level = 1, max_threshold = 1.2)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval = 5,
+                           adapt_initial_condition = true,
+                           adapt_initial_condition_only_refine = true)
+callbacks_ext = CallbackSet(amr_callback, callbacks.discrete_callbacks...)
+
+integrator = init(ode, CarpenterKennedy2N54(williamson_condition = false),
+                  dt = dt, # solve needs some value here but it will be overwritten by the stepsize_callback
+                  save_everystep = false, callback = callbacks_ext, maxiters = 100_000);
+
+# Get the last time index and work with that.
+load_timestep!(integrator, restart_filename)
+
+###############################################################################
+# run the simulation
+
+sol = solve!(integrator)
+summary_callback() # print the timer summary

examples/tree_2d_dgsem/elixir_advection_extended.jl

@@ -78,7 +78,7 @@ callbacks = CallbackSet(summary_callback,
 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; ode_default_options()...);
+            save_everystep = false, callback = callbacks);
 
 # Print the timer summary
 summary_callback()

examples/tree_2d_dgsem/elixir_advection_restart.jl

@@ -29,7 +29,7 @@ save_solution.condition.save_initial_solution = false
 
 integrator = init(ode, alg,
                   dt = dt, # solve needs some value here but it will be overwritten by the stepsize_callback
-                  callback = callbacks, maxiters = 100_000; ode_default_options()...)
+                  save_everystep = false, callback = callbacks, maxiters = 100_000)
 
 # Load saved context for adaptive time integrator
 if integrator.opts.adaptive

examples/tree_2d_dgsem/elixir_advection_restart_amr.jl

@@ -34,13 +34,13 @@ amr_controller = ControllerThreeLevel(semi, IndicatorMax(semi, variable = first)
                                       max_level = 5, max_threshold = 1.2)
 amr_callback = AMRCallback(semi, amr_controller,
                            interval = 5,
-                           adapt_initial_condition = false,
+                           adapt_initial_condition = true,
                            adapt_initial_condition_only_refine = true)
 callbacks_ext = CallbackSet(amr_callback, callbacks.discrete_callbacks...)
 
 integrator = init(ode, alg,
                   dt = dt, # solve needs some value here but it will be overwritten by the stepsize_callback
-                  callback = callbacks_ext, maxiters = 100_000; ode_default_options()...)
+                  save_everystep = false, callback = callbacks_ext, maxiters = 100_000)
 
 # Load saved context for adaptive time integrator
 if integrator.opts.adaptive

src/auxiliary/p4est.jl

@@ -99,14 +99,30 @@ end
 function load_p4est(file, ::Val{2})
     conn_vec = Vector{Ptr{p4est_connectivity_t}}(undef, 1)
     comm = P4est.uses_mpi() ? mpi_comm() : C_NULL # Use Trixi.jl's MPI communicator if p4est supports MPI
-    p4est_load_ext(file, comm, 0, 0, 1, 0, C_NULL, pointer(conn_vec))
+    p4est = p4est_load_ext(file,
+                           comm,
+                           0,         # Size of user data
+                           0,         # Flag to load user data
+                           1,         # Autopartition: ignore saved partition
+                           0,         # Have only rank 0 read headers and bcast them
+                           C_NULL,    # No pointer to user data
+                           pointer(conn_vec))
+    # p4est_load_ext only allocates memory when also data is read
+    # use p4est_reset_data to allocate uninitialized memory
+    p4est_reset_data(p4est,
+                     2 * sizeof(Int), # Use Int-Vector of size 2 as quadrant user data
+                     C_NULL,          # No init function
+                     C_NULL)          # No pointer to user data
+    return p4est
 end
 
 # 3D
 function load_p4est(file, ::Val{3})
     conn_vec = Vector{Ptr{p8est_connectivity_t}}(undef, 1)
     comm = P4est.uses_mpi() ? mpi_comm() : C_NULL # Use Trixi.jl's MPI communicator if p4est supports MPI
-    p8est_load_ext(file, comm, 0, 0, 1, 0, C_NULL, pointer(conn_vec))
+    p4est = p8est_load_ext(file, comm, 0, 0, 1, 0, C_NULL, pointer(conn_vec))
+    p8est_reset_data(p4est, 2 * sizeof(Int), C_NULL, C_NULL)
+    return p4est
 end
 
 # Read `p4est` connectivity from Abaqus mesh file (.inp)

test/test_mpi_tree.jl

@@ -74,8 +74,7 @@ CI_ON_WINDOWS = (get(ENV, "GITHUB_ACTIONS", false) == "true") && Sys.iswindows()
         @test_trixi_include(joinpath(EXAMPLES_DIR,
                                      "elixir_advection_restart_amr.jl"),
                             l2=[7.870371848717432e-5],
-                            linf=[0.0007374081713964475],
-                            coverage_override=(maxiters = 50,))
+                            linf=[0.0007374081713964475])
     end
 
     # Linear scalar advection with AMR

test/test_p4est_2d.jl

@@ -108,6 +108,20 @@ end
     end
 end
 
+@trixi_testset "elixir_advection_restart_amr.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart_amr.jl"),
+                        l2=[3.2039786506780588e-6],
+                        linf=[3.848344339896226e-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 "elixir_euler_source_terms_nonconforming_unstructured_flag.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR,
                                  "elixir_euler_source_terms_nonconforming_unstructured_flag.jl"),