Merge branch 'main' into twolayer_swe

trixi-framework · Aug 9, 2023 · 30d73b2 · 30d73b2
2 parents fae20c6 + ddf0892
commit 30d73b2
Show file tree

Hide file tree

Showing 140 changed files with 8,566 additions and 366 deletions.
diff --git a/.github/workflows/SpellCheck.yml b/.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.6
+        uses: crate-ci/typos@v1.16.2
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
@@ -69,6 +69,7 @@ jobs:
           - structured
           - p4est_part1
           - p4est_part2
+          - t8code_part1
           - unstructured_dgmulti
           - parabolic
           - paper_self_gravitating_gas_dynamics

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

diff --git a/NEWS.md b/NEWS.md
@@ -10,6 +10,7 @@ for human readability.
 
 - Experimental support for 3D parabolic diffusion terms has been added.
 - Capability to set truly discontinuous initial conditions in 1D.
+- Wetting and drying feature and examples for 1D and 2D shallow water equations
 
 #### Changed
 

diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "Trixi"
 uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
 authors = ["Michael Schlottke-Lakemper <m.schlottke-lakemper@hlrs.de>", "Gregor Gassner <ggassner@uni-koeln.de>", "Hendrik Ranocha <mail@ranocha.de>", "Andrew R. Winters <andrew.ross.winters@liu.se>", "Jesse Chan <jesse.chan@rice.edu>"]
-version = "0.5.31-pre"
+version = "0.5.38-pre"
 
 [deps]
 CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
@@ -37,6 +37,7 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 StrideArrays = "d1fa6d79-ef01-42a6-86c9-f7c551f8593b"
 StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
 SummationByPartsOperators = "9f78cca6-572e-554e-b819-917d2f1cf240"
+T8code = "d0cc0030-9a40-4274-8435-baadcfd54fa1"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 Triangulate = "f7e6ffb2-c36d-4f8f-a77e-16e897189344"
 TriplotBase = "981d1d27-644d-49a2-9326-4793e63143c3"
@@ -59,13 +60,13 @@ HDF5 = "0.14, 0.15, 0.16"
 IfElse = "0.1"
 LinearMaps = "2.7, 3.0"
 LoopVectorization = "0.12.118"
-Makie = "0.19"
 MPI = "0.20"
+Makie = "0.19"
 MuladdMacro = "0.2.2"
 Octavian = "0.3.5"
 OffsetArrays = "1.3"
 P4est = "0.4"
-Polyester = "0.3.4, 0.5, 0.6, 0.7"
+Polyester = "0.7.5"
 PrecompileTools = "1.1"
 RecipesBase = "1.1"
 Reexport = "1.0"
@@ -79,7 +80,8 @@ StaticArrayInterface = "1.4"
 StaticArrays = "1"
 StrideArrays = "0.1.18"
 StructArrays = "0.6"
-SummationByPartsOperators = "0.5.25"
+SummationByPartsOperators = "0.5.41"
+T8code = "0.4.1"
 TimerOutputs = "0.5"
 Triangulate = "2.0"
 TriplotBase = "0.1"

diff --git a/README.md b/README.md
@@ -17,6 +17,16 @@
   <img width="300px" src="https://trixi-framework.github.io/assets/logo.png">
 </p>
 
+***
+**Trixi.jl at JuliaCon 2023**<br/>
+At this year's JuliaCon, we will be present with an online contribution that involves Trixi.jl:
+
+* [Scaling Trixi.jl to more than 10,000 cores using MPI](https://pretalx.com/juliacon2023/talk/PC8PZ8/),
+  27th July 2023, 10:30–11:30 (US/Eastern), 32-G449 (Kiva)
+
+We are looking forward to seeing you there ♥️
+***
+
 **Trixi.jl** is a numerical simulation framework for hyperbolic conservation
 laws written in [Julia](https://julialang.org). A key objective for the
 framework is to be useful to both scientists and students. Therefore, next to

diff --git a/docs/literate/src/files/DGSEM_FluxDiff.jl b/docs/literate/src/files/DGSEM_FluxDiff.jl
@@ -96,13 +96,13 @@
 # \begin{align*}
 # J \underline{\dot{u}}(t) &= - M^{-1} B (\underline{f}^* - \underline{f}) - 2D \underline{f}_{vol}(u^-, u^+)\\[5pt]
 # &= - M^{-1} B (\underline{f}^* - \underline{f}_{vol}(\underline{u}, \underline{u})) - 2D \underline{f}_{vol}(u^-, u^+)\\[5pt]
-# &= - M^{-1} B \underline{f}_{sur}^* - (2D - M^{-1} B) \underline{f}_{vol}\\[5pt]
-# &= - M^{-1} B \underline{f}_{sur}^* - D_{split} \underline{f}_{vol}
+# &= - M^{-1} B \underline{f}_{surface}^* - (2D - M^{-1} B) \underline{f}_{vol}\\[5pt]
+# &= - M^{-1} B \underline{f}_{surface}^* - D_{split} \underline{f}_{vol}
 # \end{align*}
 # ```
 # This formulation is in a weak form type formulation and can be implemented by using the derivative
 # split matrix $D_{split}=(2D-M^{-1}B)$ and two different fluxes. We divide between the surface
-# flux $f=f_{sur}$ used for the numerical flux $f_{sur}^*$ and the already mentioned volume
+# flux $f=f_{surface}$ used for the numerical flux $f_{surface}^*$ and the already mentioned volume
 # flux $f_{vol}$ especially for this formulation.
 
 

diff --git a/docs/literate/src/files/shock_capturing.jl b/docs/literate/src/files/shock_capturing.jl
@@ -48,7 +48,7 @@
 # with the total energy $\mathbb{E}=\max\big(\frac{m_N^2}{\sum_{j=0}^N m_j^2}, \frac{m_{N-1}^2}{\sum_{j=0}^{N-1} m_j^2}\big)$,
 # threshold $\mathbb{T}= 0.5 * 10^{-1.8*(N+1)^{1/4}}$ and parameter $s=ln\big(\frac{1-0.0001}{0.0001}\big)\approx 9.21024$.
 
-# For computational efficiency, $\alpha_{min}$ is introduced und used for
+# For computational efficiency, $\alpha_{min}$ is introduced and used for
 # ```math
 # \tilde{\alpha} = \begin{cases}
 # 0, & \text{if } \alpha<\alpha_{min}\\

diff --git a/docs/make.jl b/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" => [

diff --git a/docs/src/restart.md b/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).
diff --git a/examples/dgmulti_2d/elixir_euler_bilinear.jl b/examples/dgmulti_2d/elixir_euler_bilinear.jl
@@ -2,7 +2,7 @@
 using Trixi, OrdinaryDiffEq
 
 dg = DGMulti(polydeg = 3, element_type = Quad(), approximation_type = SBP(),
-             surface_integral = SurfaceIntegralWeakForm(FluxHLL()),
+             surface_integral = SurfaceIntegralWeakForm(flux_hll),
              volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha))
 
 equations = CompressibleEulerEquations2D(1.4)

diff --git a/examples/dgmulti_2d/elixir_euler_curved.jl b/examples/dgmulti_2d/elixir_euler_curved.jl
@@ -2,7 +2,7 @@
 using Trixi, OrdinaryDiffEq
 
 dg = DGMulti(polydeg = 3, element_type = Quad(), approximation_type = SBP(),
-             surface_integral = SurfaceIntegralWeakForm(FluxHLL()),
+             surface_integral = SurfaceIntegralWeakForm(flux_hll),
              volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha))
 
 equations = CompressibleEulerEquations2D(1.4)
@@ -42,7 +42,8 @@ callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback)
 ###############################################################################
 # run the simulation
 
-sol = solve(ode, RDPK3SpFSAL49(); abstol=1.0e-6, reltol=1.0e-6,
+alg = RDPK3SpFSAL49()
+sol = solve(ode, alg; abstol=1.0e-6, reltol=1.0e-6,
             ode_default_options()..., callback=callbacks);
 
 summary_callback() # print the timer summary
diff --git a/examples/dgmulti_2d/elixir_euler_triangulate_pkg_mesh.jl b/examples/dgmulti_2d/elixir_euler_triangulate_pkg_mesh.jl
@@ -1,7 +1,7 @@
 using Trixi, OrdinaryDiffEq
 
 dg = DGMulti(polydeg = 3, element_type = Tri(),
-             surface_integral = SurfaceIntegralWeakForm(FluxHLL()),
+             surface_integral = SurfaceIntegralWeakForm(flux_hll),
              volume_integral = VolumeIntegralWeakForm())
 
 equations = CompressibleEulerEquations2D(1.4)

diff --git a/examples/dgmulti_2d/elixir_euler_weakform.jl b/examples/dgmulti_2d/elixir_euler_weakform.jl
@@ -2,7 +2,7 @@
 using Trixi, OrdinaryDiffEq
 
 dg = DGMulti(polydeg = 3, element_type = Tri(), approximation_type = Polynomial(),
-             surface_integral = SurfaceIntegralWeakForm(FluxHLL()),
+             surface_integral = SurfaceIntegralWeakForm(flux_hll),
              volume_integral = VolumeIntegralWeakForm())
 
 equations = CompressibleEulerEquations2D(1.4)

diff --git a/examples/dgmulti_2d/elixir_euler_weakform_periodic.jl b/examples/dgmulti_2d/elixir_euler_weakform_periodic.jl
@@ -2,7 +2,7 @@
 using Trixi, OrdinaryDiffEq
 
 dg = DGMulti(polydeg = 3, element_type = Tri(), approximation_type = Polynomial(),
-             surface_integral = SurfaceIntegralWeakForm(FluxHLL()),
+             surface_integral = SurfaceIntegralWeakForm(flux_hll),
              volume_integral = VolumeIntegralWeakForm())
 
 equations = CompressibleEulerEquations2D(1.4)

diff --git a/examples/dgmulti_3d/elixir_euler_curved.jl b/examples/dgmulti_3d/elixir_euler_curved.jl
@@ -2,7 +2,7 @@
 using Trixi, OrdinaryDiffEq
 
 dg = DGMulti(polydeg = 3, element_type = Hex(), approximation_type=SBP(),
-             surface_integral = SurfaceIntegralWeakForm(FluxHLL()),
+             surface_integral = SurfaceIntegralWeakForm(flux_hll),
              volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha))
 
 equations = CompressibleEulerEquations3D(1.4)

diff --git a/examples/dgmulti_3d/elixir_euler_weakform.jl b/examples/dgmulti_3d/elixir_euler_weakform.jl
@@ -2,7 +2,7 @@
 using Trixi, OrdinaryDiffEq
 
 dg = DGMulti(polydeg = 3, element_type = Tet(),
-             surface_integral = SurfaceIntegralWeakForm(FluxHLL()),
+             surface_integral = SurfaceIntegralWeakForm(flux_hll),
              volume_integral = VolumeIntegralWeakForm())
 
 equations = CompressibleEulerEquations3D(1.4)

diff --git a/examples/dgmulti_3d/elixir_euler_weakform_periodic.jl b/examples/dgmulti_3d/elixir_euler_weakform_periodic.jl
@@ -2,7 +2,7 @@
 using Trixi, OrdinaryDiffEq
 
 dg = DGMulti(polydeg = 3, element_type = Tet(), approximation_type = Polynomial(),
-             surface_integral = SurfaceIntegralWeakForm(FluxHLL()),
+             surface_integral = SurfaceIntegralWeakForm(flux_hll),
              volume_integral = VolumeIntegralWeakForm())
 
 equations = CompressibleEulerEquations3D(1.4)

diff --git a/examples/p4est_2d_dgsem/elixir_advection_restart.jl b/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