Save and load user_data for p4est (#1915)

* save and load user_data for p4est

this fixes restating when using AMR

* use p4est_reset_data instead

* add test for p4est, restart, and AMR

* adapt errors in test

* 3D = p8!

* return something

* remove coverage_override

* remove coverage_override

* use ode_default_options()

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,8 @@ 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()...);
+            callback = callbacks;
+            ode_default_options()...); # default options because an adaptive time stepping method is used in test_mpi_tree.jl
 
 # Print the timer summary
 summary_callback()

examples/tree_2d_dgsem/elixir_advection_restart.jl

@@ -29,7 +29,8 @@ 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()...)
+                  callback = callbacks;
+                  ode_default_options()...); # default options because an adaptive time stepping method is used in test_mpi_tree.jl
 
 # 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

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

test/test_p4est_2d.jl

@@ -97,7 +97,24 @@ end
                         linf=[6.21489667023134e-5],
                         # With the default `maxiters = 1` in coverage tests,
                         # there would be no time steps after the restart.
-                        coverage_override=(maxiters = 100_000,))
+                        coverage_override=(maxiters = 25,))
+    # 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_advection_restart_amr.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart_amr.jl"),
+                        l2=[2.869137983727866e-6],
+                        linf=[3.8353423270964804e-5],
+                        # With the default `maxiters = 1` in coverage tests,
+                        # there would be no time steps after the restart.
+                        coverage_override=(maxiters = 25,))
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
     let