Merge branch 'main' into main

Project.toml

@@ -1,7 +1,7 @@
 name = "Trixi"
 uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
 authors = ["Michael Schlottke-Lakemper <[email protected]>", "Gregor Gassner <[email protected]>", "Hendrik Ranocha <[email protected]>", "Andrew R. Winters <[email protected]>", "Jesse Chan <[email protected]>"]
-version = "0.7.15-pre"
+version = "0.7.16-pre"
 
 [deps]
 CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
@@ -111,4 +111,4 @@ julia = "1.8"
 [extras]
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 Convex = "f65535da-76fb-5f13-bab9-19810c17039a"
-ECOS = "e2685f51-7e38-5353-a97d-a921fd2c8199"
+ECOS = "e2685f51-7e38-5353-a97d-a921fd2c8199"

examples/tree_2d_dgsem/elixir_advection_restart_amr.jl

@@ -0,0 +1,58 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# 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, 3.0))
+
+###############################################################################
+# 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_000040.h5")
+mesh = load_mesh(restart_filename)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+tspan = (load_time(restart_filename), 5.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 = 3,
+                                      med_level = 4, med_threshold = 0.8,
+                                      max_level = 5, max_threshold = 1.2)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval = 5,
+                           adapt_initial_condition = false,
+                           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()...)
+
+# 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)
+
+###############################################################################
+# run the simulation
+
+sol = solve!(integrator)
+
+summary_callback() # print the timer summary

src/callbacks_step/alive.jl

@@ -68,7 +68,6 @@ function (alive_callback::AliveCallback)(u, t, integrator)
     # We need to check the number of accepted steps since callbacks are not
     # activated after a rejected step.
     return alive_interval > 0 && ((integrator.stats.naccept % alive_interval == 0 &&
-             !(integrator.stats.naccept == 0 && integrator.iter > 0) &&
              (analysis_interval == 0 ||
               integrator.stats.naccept % analysis_interval != 0)) ||
             isfinished(integrator))

src/callbacks_step/analysis.jl

@@ -26,8 +26,8 @@ or `extra_analysis_errors = (:conservation_error,)`.
 If you want to omit the computation (to safe compute-time) of the [`default_analysis_errors`](@ref), specify
 `analysis_errors = Symbol[]`.
 Note: `default_analysis_errors` are `:l2_error` and `:linf_error` for all equations.
-If you want to compute `extra_analysis_errors` such as `:conservation_error` solely, i.e., 
-without `:l2_error, :linf_error` you need to specify 
+If you want to compute `extra_analysis_errors` such as `:conservation_error` solely, i.e.,
+without `:l2_error, :linf_error` you need to specify
 `analysis_errors = [:conservation_error]` instead of `extra_analysis_errors = [:conservation_error]`.
 
 Further scalar functions `func` in `extra_analysis_integrals` are applied to the numerical
@@ -119,9 +119,7 @@ function AnalysisCallback(mesh, equations::AbstractEquations, solver, cache;
     # We need to check the number of accepted steps since callbacks are not
     # activated after a rejected step.
     condition = (u, t, integrator) -> interval > 0 &&
-        ((integrator.stats.naccept % interval == 0 &&
-          !(integrator.stats.naccept == 0 && integrator.iter > 0)) ||
-         isfinished(integrator))
+        (integrator.stats.naccept % interval == 0 || isfinished(integrator))
 
     analyzer = SolutionAnalyzer(solver; kwargs...)
     cache_analysis = create_cache_analysis(analyzer, mesh, equations, solver, cache,
@@ -696,7 +694,7 @@ include("analysis_dg3d_parallel.jl")
 
 # Special analyze for `SemidiscretizationHyperbolicParabolic` such that
 # precomputed gradients are available. Required for `enstrophy` (see above) and viscous forces.
-# Note that this needs to be included after `analysis_surface_integral_2d.jl` to 
+# Note that this needs to be included after `analysis_surface_integral_2d.jl` to
 # have `VariableViscous` available.
 function analyze(quantity::AnalysisSurfaceIntegral{Variable},
                  du, u, t,

src/callbacks_step/save_restart.jl

@@ -83,8 +83,7 @@ function (restart_callback::SaveRestartCallback)(u, t, integrator)
     #    (total #steps)       (#accepted steps)
     # We need to check the number of accepted steps since callbacks are not
     # activated after a rejected step.
-    return interval > 0 && (((integrator.stats.naccept % interval == 0) &&
-             !(integrator.stats.naccept == 0 && integrator.iter > 0)) ||
+    return interval > 0 && (integrator.stats.naccept % interval == 0 ||
             (save_final_restart && isfinished(integrator)))
 end
 
@@ -198,7 +197,7 @@ function load_adaptive_time_integrator!(integrator, restart_file::AbstractString
         # 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 
+        if hasproperty(controller, :err) # Distinguish PIDController from PIController
             controller.err[:] = read(attributes(file)["time_integrator_controller_err"])
         end
     end

src/callbacks_step/save_solution.jl

@@ -174,8 +174,7 @@ function (solution_callback::SaveSolutionCallback)(u, t, integrator)
     #    (total #steps)       (#accepted steps)
     # We need to check the number of accepted steps since callbacks are not
     # activated after a rejected step.
-    return interval_or_dt > 0 && (((integrator.stats.naccept % interval_or_dt == 0) &&
-             !(integrator.stats.naccept == 0 && integrator.iter > 0)) ||
+    return interval_or_dt > 0 && (integrator.stats.naccept % interval_or_dt == 0 ||
             (save_final_solution && isfinished(integrator)))
 end
 

src/callbacks_step/visualization.jl

@@ -137,8 +137,7 @@ function (visualization_callback::VisualizationCallback)(u, t, integrator)
     #    (total #steps)       (#accepted steps)
     # We need to check the number of accepted steps since callbacks are not
     # activated after a rejected step.
-    return interval > 0 && ((integrator.stats.naccept % interval == 0 &&
-             !(integrator.stats.naccept == 0 && integrator.iter > 0)) ||
+    return interval > 0 && (integrator.stats.naccept % interval == 0 ||
             isfinished(integrator))
 end
 

src/meshes/mesh_io.jl

@@ -370,17 +370,18 @@ function load_mesh_parallel(mesh_file::AbstractString; n_cells_max, RealT)
 
     if mesh_type == "TreeMesh"
         if mpi_isroot()
-            n_cells = h5open(mesh_file, "r") do file
-                read(attributes(file)["n_cells"])
+            n_cells, capacity = h5open(mesh_file, "r") do file
+                return read(attributes(file)["n_cells"]),
+                       read(attributes(file)["capacity"])
             end
-            MPI.Bcast!(Ref(ndims_), mpi_root(), mpi_comm())
             MPI.Bcast!(Ref(n_cells), mpi_root(), mpi_comm())
+            MPI.Bcast!(Ref(capacity), mpi_root(), mpi_comm())
         else
-            ndims_ = MPI.Bcast!(Ref(0), mpi_root(), mpi_comm())[]
             n_cells = MPI.Bcast!(Ref(0), mpi_root(), mpi_comm())[]
+            capacity = MPI.Bcast!(Ref(0), mpi_root(), mpi_comm())[]
         end
 
-        mesh = TreeMesh(ParallelTree{ndims_}, max(n_cells, n_cells_max))
+        mesh = TreeMesh(ParallelTree{ndims_}, max(n_cells, n_cells_max, capacity))
         load_mesh!(mesh, mesh_file)
     elseif mesh_type == "P4estMesh"
         if mpi_isroot()

test/test_mpi_tree.jl

@@ -70,6 +70,14 @@ CI_ON_WINDOWS = (get(ENV, "GITHUB_ACTIONS", false) == "true") && Sys.iswindows()
                             coverage_override=(maxiters = 6,))
     end
 
+    @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,))
+    end
+
     # Linear scalar advection with AMR
     # These example files are only for testing purposes and have no practical use
     @trixi_testset "elixir_advection_amr_refine_twice.jl" begin