Merge branch 'main' into HLL_2_Wave_Improvements_NonBreaking

.github/workflows/SpellCheck.yml

@@ -10,4 +10,4 @@ jobs:
       - name: Checkout Actions Repository
         uses: actions/checkout@v3
       - name: Check spelling
-        uses: crate-ci/typos@v1.15.10
+        uses: crate-ci/typos@v1.16.0

AUTHORS.md

@@ -24,6 +24,7 @@ are listed in alphabetical order:
 
 * Maximilian D. Bertrand
 * Benjamin Bolm
+* Simon Candelaresi
 * Jesse Chan
 * Lars Christmann
 * Christof Czernik

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.5.31-pre"
+version = "0.5.32-pre"
 
 [deps]
 CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"

docs/make.jl

@@ -92,6 +92,7 @@ makedocs(
         "Getting started" => [
             "Overview" => "overview.md",
             "Visualization" => "visualization.md",
+            "Restart simulation" => "restart.md",
         ],
         "Tutorials" => tutorials,
         "Basic building blocks" => [

docs/src/restart.md

@@ -0,0 +1,89 @@
+# [Restart simulation](@id restart)
+
+You can continue running an already finished simulation by first
+preparing the simulation for the restart and then performing the restart.
+Here we suppose that in the first run your simulation stops at time 1.0
+and then you want it to run further to time 2.0.
+
+## [Prepare the simulation for a restart](@id restart_preparation)
+In you original elixir you need to specify to write out restart files.
+Those will later be read for the restart of your simulation.
+This is done almost the same way as writing the snapshots using the
+[`SaveSolutionCallback`](@ref) callback.
+For the restart files it is called [`SaveRestartCallback`](@ref):
+```julia
+save_restart = SaveRestartCallback(interval=100,
+                                   save_final_restart=true)
+```
+Make this part of your `CallbackSet`.
+
+An example is
+[```examples/examples/structured_2d_dgsem/elixir_advection_extended.jl```](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/structured_2d_dgsem/elixir_advection_extended.jl).
+
+
+## [Perform the simulation restart](@id restart_perform)
+Since all of the information about the simulation can be obtained from the
+last snapshot, the restart can be done with relatively few lines
+in an extra elixir file.
+However, some might prefer to keep everything in one elixir and
+conditionals like ```if restart``` with a boolean variable ```restart``` that is user defined.
+
+First we need to define from which file we want to restart, e.g.
+```julia
+restart_file = "restart_000021.h5"
+restart_filename = joinpath("out", restart_file)
+```
+
+Then we load the mesh file:
+```julia
+mesh = load_mesh(restart_filename)
+```
+
+This is then needed for the semidiscretization:
+```julia
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+```
+
+We then define a new time span for the simulation that takes as starting
+time the one form the snapshot:
+```julia
+tspan = (load_time(restart_filename), 2.0)
+```
+
+We now also take the last ```dt```, so that our solver does not need to first find
+one to fulfill the CFL condition:
+```julia
+dt = load_dt(restart_filename)
+```
+
+The ODE that we will pass to the solver is now:
+```julia
+ode = semidiscretize(semi, tspan, restart_filename)
+```
+
+You should now define a [`SaveSolutionCallback`](@ref) similar to the
+[original simulation](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/structured_2d_dgsem/elixir_advection_extended.jl),
+but with ```save_initial_solution=false```, otherwise our initial snapshot will be overwritten.
+If you are using one file for the original simulation and the restart
+you can reuse your [`SaveSolutionCallback`](@ref), but need to set
+```julia
+save_solution.condition.save_initial_solution = false
+```
+
+Before we compute the solution using 
+[OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl)
+we need to set the integrator
+and its time step number, e.g.:
+```julia
+integrator = init(ode, CarpenterKennedy2N54(williamson_condition=false),
+                  dt=dt, save_everystep=false, callback=callbacks);
+integrator.iter = load_timestep(restart_filename)
+integrator.stats.naccept = integrator.iter
+```
+
+Now we can compute the solution:
+```julia
+sol = solve!(integrator)
+```
+
+An example is in `[``examples/structured_2d_dgsem/elixir_advection_restart.jl```](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/structured_2d_dgsem/elixir_advection_restart.jl).

examples/p4est_2d_dgsem/elixir_advection_restart.jl

@@ -24,13 +24,23 @@ 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
+
+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);
+
+# Get the last time index and work with that.
+integrator.iter = load_timestep(restart_filename)
+integrator.stats.naccept = integrator.iter
+
 
 ###############################################################################
 # run the simulation
 
-sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
-            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
-            save_everystep=false, callback=callbacks);
+sol = solve!(integrator)
 summary_callback() # print the timer summary

examples/p4est_2d_dgsem/elixir_linearizedeuler_gaussian_source.jl

@@ -0,0 +1,89 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+# Based on the TreeMesh example `elixir_acoustics_gaussian_source.jl`.
+# The acoustic perturbation equations have been replaced with the linearized Euler
+# equations and instead of the Cartesian `TreeMesh` a rotated `P4estMesh` is used
+
+# Oscillating Gaussian-shaped source terms
+function source_terms_gauss(u, x, t, equations::LinearizedEulerEquations2D)
+  r = 0.1
+  A = 1.0
+  f = 2.0
+
+  # Velocity sources
+  s2 = 0.0
+  s3 = 0.0
+  # Density and pressure source
+  s1 = s4 = exp(-(x[1]^2 + x[2]^2) / (2 * r^2)) * A * sin(2 * pi * f * t)
+
+  return SVector(s1, s2, s3, s4)
+end
+
+initial_condition_zero(x, t, equations::LinearizedEulerEquations2D) = SVector(0.0, 0.0, 0.0, 0.0)
+
+###############################################################################
+# semidiscretization of the linearized Euler equations
+
+# Create a domain that is a 30° rotated version of [-3, 3]^2
+c = cospi(2 * 30.0 / 360.0)
+s = sinpi(2 * 30.0 / 360.0)
+rot_mat = Trixi.SMatrix{2, 2}([c -s; s c])
+mapping(xi, eta) = rot_mat * SVector(3.0*xi, 3.0*eta)
+
+# Mean density and speed of sound are slightly off from 1.0 to allow proper verification of
+# curved LEE implementation using this elixir (some things in the LEE cancel if both are 1.0)
+equations = LinearizedEulerEquations2D(v_mean_global=Tuple(rot_mat * SVector(-0.5, 0.25)),
+                                       c_mean_global=1.02, rho_mean_global=1.01)
+
+initial_condition = initial_condition_zero
+
+# Create DG solver with polynomial degree = 3 and upwind flux as surface flux
+solver = DGSEM(polydeg=3, surface_flux=flux_godunov)
+
+# Create a uniformly refined mesh with periodic boundaries
+trees_per_dimension = (4, 4)
+mesh = P4estMesh(trees_per_dimension, polydeg=1,
+                 mapping=mapping,
+                 periodicity=true, initial_refinement_level=2)
+
+# A semidiscretization collects data structures and functions for the spatial discretization
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms=source_terms_gauss)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+# Create ODE problem with time span from 0.0 to 2.0
+tspan = (0.0, 2.0)
+ode = semidiscretize(semi, tspan)
+
+# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
+# and resets the timers
+summary_callback = SummaryCallback()
+
+# The AnalysisCallback allows to analyse the solution in regular intervals and prints the results
+analysis_callback = AnalysisCallback(semi, interval=100)
+
+# The SaveSolutionCallback allows to save the solution to a file in regular intervals
+save_solution = SaveSolutionCallback(interval=100)
+
+# The StepsizeCallback handles the re-calculation of the maximum Δt after each time step
+stepsize_callback = StepsizeCallback(cfl=0.5)
+
+# Create a CallbackSet to collect all callbacks such that they can be passed to the ODE solver
+callbacks = CallbackSet(summary_callback, analysis_callback, save_solution, stepsize_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),
+            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep=false, callback=callbacks);
+
+# Print the timer summary
+summary_callback()

examples/p4est_3d_dgsem/elixir_advection_restart.jl

@@ -21,13 +21,23 @@ mesh = load_mesh(restart_filename)
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, solver)
 
 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
+
+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);
+
+# Get the last time index and work with that.
+integrator.iter = load_timestep(restart_filename)
+integrator.stats.naccept = integrator.iter
+
 
 ###############################################################################
 # run the simulation
 
-sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
-            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
-            save_everystep=false, callback=callbacks);
+sol = solve!(integrator)
 summary_callback() # print the timer summary

examples/structured_2d_dgsem/elixir_advection_restart.jl

@@ -23,13 +23,22 @@ mesh = load_mesh(restart_filename)
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
 
 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
+
+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);
+
+# Get the last time index and work with that.
+integrator.iter = load_timestep(restart_filename)
+integrator.stats.naccept = integrator.iter
 
 ###############################################################################
 # run the simulation
 
-sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
-            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
-            save_everystep=false, callback=callbacks);
+sol = solve!(integrator)
 summary_callback() # print the timer summary

examples/structured_3d_dgsem/elixir_advection_restart.jl

@@ -21,13 +21,23 @@ mesh = load_mesh(restart_filename)
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, solver)
 
 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
+
+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);
+
+# Get the last time index and work with that.
+integrator.iter = load_timestep(restart_filename)
+integrator.stats.naccept = integrator.iter
+
 
 ###############################################################################
 # run the simulation
 
-sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
-            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
-            save_everystep=false, callback=callbacks);
+sol = solve!(integrator)
 summary_callback() # print the timer summary

examples/tree_2d_dgsem/elixir_advection_restart.jl

@@ -20,13 +20,23 @@ mesh = load_mesh(restart_filename)
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
 
 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
+
+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)
+
+# Get the last time index and work with that.
+integrator.iter = load_timestep(restart_filename)
+integrator.stats.naccept = integrator.iter
 
 ###############################################################################
 # run the simulation
 
-sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
-            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
-            save_everystep=false, callback=callbacks);
+sol = solve!(integrator)
+
 summary_callback() # print the timer summary