Merge branch 'main' into CommentsAllocs

.github/workflows/FormatCheck.yml

@@ -29,7 +29,7 @@ jobs:
         # TODO: Change the call below to
         #       format(".")
         run: |
-          julia  -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter", version="1.0.45"))'
+          julia  -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter", version="1.0.60"))'
           julia  -e 'using JuliaFormatter; format(["benchmark", "examples", "ext", "src", "test", "utils"])'
       - name: Format check
         run: |

Project.toml

@@ -1,9 +1,10 @@
 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.8.9-DEV"
+version = "0.8.11-DEV"
 
 [deps]
+Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
 CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
 ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
@@ -60,6 +61,7 @@ TrixiConvexECOSExt = ["Convex", "ECOS"]
 TrixiMakieExt = "Makie"
 
 [compat]
+Accessors = "0.1.12"
 CodeTracking = "1.0.5"
 ConstructionBase = "1.3"
 Convex = "0.16"

docs/src/styleguide.md

@@ -51,7 +51,7 @@ PRs that verify that running JuliaFormatter.jl again will not change the source
 To format your contributions before created a PR (or, at least, before requesting a review
 of your PR), you need to install JuliaFormatter.jl first by running
 ```shell
-julia -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter", version="1.0.45"))'
+julia -e 'using Pkg; Pkg.add(PackageSpec(name = "JuliaFormatter", version="1.0.60"))'
 ```
 You can then recursively format the core Julia files in the Trixi.jl repo by executing
 ```shell

src/Trixi.jl

@@ -18,6 +18,7 @@ module Trixi
 # Include other packages that are used in Trixi.jl
 # (standard library packages first, other packages next, all of them sorted alphabetically)
 
+using Accessors: @reset
 using LinearAlgebra: LinearAlgebra, Diagonal, diag, dot, mul!, norm, cross, normalize, I,
                      UniformScaling, det
 using Printf: @printf, @sprintf, println

src/auxiliary/math.jl

@@ -284,6 +284,11 @@ end
 # when using `@fastmath`, which we also get from
 # [Fortran](https://godbolt.org/z/Yrsa1js7P)
 # or [C++](https://godbolt.org/z/674G7Pccv).
+#
+# Note however that such a custom reimplementation can cause incompatibilities with other
+# packages. Currently we are affected by an issue with MPI.jl on ARM, see
+# https://github.com/trixi-framework/Trixi.jl/issues/1922
+# The workaround is to resort to Base.min / Base.max when using MPI reductions.
 """
     Trixi.max(x, y, ...)
 

src/callbacks_step/alive.jl

@@ -45,7 +45,7 @@ function Base.show(io::IO, ::MIME"text/plain",
         alive_callback = cb.affect!
 
         setup = [
-            "interval" => alive_callback.alive_interval,
+            "interval" => alive_callback.alive_interval
         ]
         summary_box(io, "AliveCallback", setup)
     end

src/callbacks_step/analysis.jl

@@ -434,7 +434,8 @@ function (analysis_callback::AnalysisCallback)(io, du, u, u_ode, t, semi)
             res = maximum(abs, view(du, v, ..))
             if mpi_isparallel()
                 # TODO: Debugging, here is a type instability
-                global_res = MPI.Reduce!(Ref(res), max, mpi_root(), mpi_comm())
+                # Base.max instead of max needed, see comment in src/auxiliary/math.jl
+                global_res = MPI.Reduce!(Ref(res), Base.max, mpi_root(), mpi_comm())
                 if mpi_isroot()
                     res::eltype(du) = global_res[]
                 end

src/callbacks_step/analysis_dg2d.jl

@@ -29,10 +29,39 @@ function create_cache_analysis(analyzer, mesh::TreeMesh{2},
     return (; u_local, u_tmp1, x_local, x_tmp1)
 end
 
+# Specialized cache for P4estMesh to allow for different ambient dimension from mesh dimension
+function create_cache_analysis(analyzer, mesh::P4estMesh{2, NDIMS_AMBIENT},
+                               equations, dg::DG, cache,
+                               RealT, uEltype) where {NDIMS_AMBIENT}
+
+    # pre-allocate buffers
+    # We use `StrideArray`s here since these buffers are used in performance-critical
+    # places and the additional information passed to the compiler makes them faster
+    # than native `Array`s.
+    u_local = StrideArray(undef, uEltype,
+                          StaticInt(nvariables(equations)), StaticInt(nnodes(analyzer)),
+                          StaticInt(nnodes(analyzer)))
+    u_tmp1 = StrideArray(undef, uEltype,
+                         StaticInt(nvariables(equations)), StaticInt(nnodes(analyzer)),
+                         StaticInt(nnodes(dg)))
+    x_local = StrideArray(undef, RealT,
+                          StaticInt(NDIMS_AMBIENT), StaticInt(nnodes(analyzer)),
+                          StaticInt(nnodes(analyzer)))
+    x_tmp1 = StrideArray(undef, RealT,
+                         StaticInt(NDIMS_AMBIENT), StaticInt(nnodes(analyzer)),
+                         StaticInt(nnodes(dg)))
+    jacobian_local = StrideArray(undef, RealT,
+                                 StaticInt(nnodes(analyzer)),
+                                 StaticInt(nnodes(analyzer)))
+    jacobian_tmp1 = StrideArray(undef, RealT,
+                                StaticInt(nnodes(analyzer)), StaticInt(nnodes(dg)))
+
+    return (; u_local, u_tmp1, x_local, x_tmp1, jacobian_local, jacobian_tmp1)
+end
+
 function create_cache_analysis(analyzer,
                                mesh::Union{StructuredMesh{2}, StructuredMeshView{2},
-                                           UnstructuredMesh2D,
-                                           P4estMesh{2}, T8codeMesh{2}},
+                                           UnstructuredMesh2D, T8codeMesh{2}},
                                equations, dg::DG, cache,
                                RealT, uEltype)
 

src/callbacks_step/analysis_dg2d_parallel.jl

@@ -131,7 +131,8 @@ function calc_error_norms(func, u, t, analyzer,
     global_l2_error = Vector(l2_error)
     global_linf_error = Vector(linf_error)
     MPI.Reduce!(global_l2_error, +, mpi_root(), mpi_comm())
-    MPI.Reduce!(global_linf_error, max, mpi_root(), mpi_comm())
+    # Base.max instead of max needed, see comment in src/auxiliary/math.jl
+    MPI.Reduce!(global_linf_error, Base.max, mpi_root(), mpi_comm())
     total_volume = MPI.Reduce(volume, +, mpi_root(), mpi_comm())
     if mpi_isroot()
         l2_error = convert(typeof(l2_error), global_l2_error)

src/callbacks_step/analysis_dg3d.jl

@@ -35,9 +35,51 @@ function create_cache_analysis(analyzer, mesh::TreeMesh{3},
     return (; u_local, u_tmp1, u_tmp2, x_local, x_tmp1, x_tmp2)
 end
 
+# Specialized cache for P4estMesh to allow for different ambient dimension from mesh dimension
 function create_cache_analysis(analyzer,
-                               mesh::Union{StructuredMesh{3}, P4estMesh{3},
-                                           T8codeMesh{3}},
+                               mesh::P4estMesh{3, NDIMS_AMBIENT},
+                               equations, dg::DG, cache,
+                               RealT, uEltype) where {NDIMS_AMBIENT}
+
+    # pre-allocate buffers
+    # We use `StrideArray`s here since these buffers are used in performance-critical
+    # places and the additional information passed to the compiler makes them faster
+    # than native `Array`s.
+    u_local = StrideArray(undef, uEltype,
+                          StaticInt(nvariables(equations)), StaticInt(nnodes(analyzer)),
+                          StaticInt(nnodes(analyzer)), StaticInt(nnodes(analyzer)))
+    u_tmp1 = StrideArray(undef, uEltype,
+                         StaticInt(nvariables(equations)), StaticInt(nnodes(analyzer)),
+                         StaticInt(nnodes(dg)), StaticInt(nnodes(dg)))
+    u_tmp2 = StrideArray(undef, uEltype,
+                         StaticInt(nvariables(equations)), StaticInt(nnodes(analyzer)),
+                         StaticInt(nnodes(analyzer)), StaticInt(nnodes(dg)))
+    x_local = StrideArray(undef, RealT,
+                          StaticInt(NDIMS_AMBIENT), StaticInt(nnodes(analyzer)),
+                          StaticInt(nnodes(analyzer)), StaticInt(nnodes(analyzer)))
+    x_tmp1 = StrideArray(undef, RealT,
+                         StaticInt(NDIMS_AMBIENT), StaticInt(nnodes(analyzer)),
+                         StaticInt(nnodes(dg)), StaticInt(nnodes(dg)))
+    x_tmp2 = StrideArray(undef, RealT,
+                         StaticInt(NDIMS_AMBIENT), StaticInt(nnodes(analyzer)),
+                         StaticInt(nnodes(analyzer)), StaticInt(nnodes(dg)))
+    jacobian_local = StrideArray(undef, RealT,
+                                 StaticInt(nnodes(analyzer)),
+                                 StaticInt(nnodes(analyzer)),
+                                 StaticInt(nnodes(analyzer)))
+    jacobian_tmp1 = StrideArray(undef, RealT,
+                                StaticInt(nnodes(analyzer)), StaticInt(nnodes(dg)),
+                                StaticInt(nnodes(dg)))
+    jacobian_tmp2 = StrideArray(undef, RealT,
+                                StaticInt(nnodes(analyzer)),
+                                StaticInt(nnodes(analyzer)), StaticInt(nnodes(dg)))
+
+    return (; u_local, u_tmp1, u_tmp2, x_local, x_tmp1, x_tmp2, jacobian_local,
+            jacobian_tmp1, jacobian_tmp2)
+end
+
+function create_cache_analysis(analyzer,
+                               mesh::Union{StructuredMesh{3}, T8codeMesh{3}},
                                equations, dg::DG, cache,
                                RealT, uEltype)
 

src/callbacks_step/analysis_dg3d_parallel.jl

@@ -49,7 +49,8 @@ function calc_error_norms(func, u, t, analyzer,
     global_l2_error = Vector(l2_error)
     global_linf_error = Vector(linf_error)
     MPI.Reduce!(global_l2_error, +, mpi_root(), mpi_comm())
-    MPI.Reduce!(global_linf_error, max, mpi_root(), mpi_comm())
+    # Base.max instead of max needed, see comment in src/auxiliary/math.jl
+    MPI.Reduce!(global_linf_error, Base.max, mpi_root(), mpi_comm())
     total_volume = MPI.Reduce(volume, +, mpi_root(), mpi_comm())
     if mpi_isroot()
         l2_error = convert(typeof(l2_error), global_l2_error)

src/callbacks_step/averaging.jl

@@ -45,7 +45,7 @@ function Base.show(io::IO, ::MIME"text/plain",
 
         setup = [
             "Start time" => first(averaging_callback.tspan),
-            "Final time" => last(averaging_callback.tspan),
+            "Final time" => last(averaging_callback.tspan)
         ]
         summary_box(io, "AveragingCallback", setup)
     end

src/callbacks_step/glm_speed.jl

@@ -48,7 +48,7 @@ function Base.show(io::IO, ::MIME"text/plain",
         setup = [
             "GLM wave speed scaling" => glm_speed_callback.glm_scale,
             "Expected CFL number" => glm_speed_callback.cfl,
-            "Selected semidiscretizations" => glm_speed_callback.semi_indices,
+            "Selected semidiscretizations" => glm_speed_callback.semi_indices
         ]
         summary_box(io, "GlmSpeedCallback", setup)
     end
@@ -83,7 +83,10 @@ function update_cleaning_speed!(semi, glm_speed_callback, dt)
     c_h_deltat = calc_dt_for_cleaning_speed(cfl, mesh, equations, solver, cache)
 
     # c_h is proportional to its own time step divided by the complete MHD time step
-    equations.c_h = glm_scale * c_h_deltat / dt
+    # We use @reset here since the equations are immutable (to work on GPUs etc.).
+    # Thus, we need to modify the equations field of the semidiscretization.
+    @reset equations.c_h = glm_scale * c_h_deltat / dt
+    semi.equations = equations
 
     return semi
 end

src/callbacks_step/save_restart.jl

@@ -38,7 +38,7 @@ function Base.show(io::IO, ::MIME"text/plain",
             "interval" => save_restart_callback.interval,
             "save final solution" => save_restart_callback.save_final_restart ? "yes" :
                                      "no",
-            "output directory" => abspath(normpath(save_restart_callback.output_directory)),
+            "output directory" => abspath(normpath(save_restart_callback.output_directory))
         ]
         summary_box(io, "SaveRestartCallback", setup)
     end

src/callbacks_step/save_solution.jl

@@ -62,7 +62,7 @@ function Base.show(io::IO, ::MIME"text/plain",
                                        "yes" : "no",
             "save final solution" => save_solution_callback.save_final_solution ?
                                      "yes" : "no",
-            "output directory" => abspath(normpath(save_solution_callback.output_directory)),
+            "output directory" => abspath(normpath(save_solution_callback.output_directory))
         ]
         summary_box(io, "SaveSolutionCallback", setup)
     end
@@ -85,7 +85,7 @@ function Base.show(io::IO, ::MIME"text/plain",
                                        "yes" : "no",
             "save final solution" => save_solution_callback.save_final_solution ?
                                      "yes" : "no",
-            "output directory" => abspath(normpath(save_solution_callback.output_directory)),
+            "output directory" => abspath(normpath(save_solution_callback.output_directory))
         ]
         summary_box(io, "SaveSolutionCallback", setup)
     end

src/callbacks_step/steady_state.jl

@@ -43,7 +43,7 @@ function Base.show(io::IO, ::MIME"text/plain",
 
         setup = [
             "absolute tolerance" => steady_state_callback.abstol,
-            "relative tolerance" => steady_state_callback.reltol,
+            "relative tolerance" => steady_state_callback.reltol
         ]
         summary_box(io, "SteadyStateCallback", setup)
     end

src/callbacks_step/stepsize.jl

@@ -33,7 +33,7 @@ function Base.show(io::IO, ::MIME"text/plain",
         stepsize_callback = cb.affect!
 
         setup = [
-            "CFL number" => stepsize_callback.cfl_number,
+            "CFL number" => stepsize_callback.cfl_number
         ]
         summary_box(io, "StepsizeCallback", setup)
     end

src/callbacks_step/stepsize_dg2d.jl

@@ -54,7 +54,8 @@ function max_dt(u, t, mesh::ParallelTreeMesh{2},
                       typeof(constant_speed), typeof(equations), typeof(dg),
                       typeof(cache)},
                 u, t, mesh, constant_speed, equations, dg, cache)
-    dt = MPI.Allreduce!(Ref(dt), min, mpi_comm())[]
+    # Base.min instead of min needed, see comment in src/auxiliary/math.jl
+    dt = MPI.Allreduce!(Ref(dt), Base.min, mpi_comm())[]
 
     return dt
 end
@@ -70,7 +71,8 @@ function max_dt(u, t, mesh::ParallelTreeMesh{2},
                       typeof(constant_speed), typeof(equations), typeof(dg),
                       typeof(cache)},
                 u, t, mesh, constant_speed, equations, dg, cache)
-    dt = MPI.Allreduce!(Ref(dt), min, mpi_comm())[]
+    # Base.min instead of min needed, see comment in src/auxiliary/math.jl
+    dt = MPI.Allreduce!(Ref(dt), Base.min, mpi_comm())[]
 
     return dt
 end
@@ -154,7 +156,8 @@ function max_dt(u, t, mesh::ParallelP4estMesh{2},
                       typeof(constant_speed), typeof(equations), typeof(dg),
                       typeof(cache)},
                 u, t, mesh, constant_speed, equations, dg, cache)
-    dt = MPI.Allreduce!(Ref(dt), min, mpi_comm())[]
+    # Base.min instead of min needed, see comment in src/auxiliary/math.jl
+    dt = MPI.Allreduce!(Ref(dt), Base.min, mpi_comm())[]
 
     return dt
 end
@@ -170,7 +173,8 @@ function max_dt(u, t, mesh::ParallelP4estMesh{2},
                       typeof(constant_speed), typeof(equations), typeof(dg),
                       typeof(cache)},
                 u, t, mesh, constant_speed, equations, dg, cache)
-    dt = MPI.Allreduce!(Ref(dt), min, mpi_comm())[]
+    # Base.min instead of min needed, see comment in src/auxiliary/math.jl
+    dt = MPI.Allreduce!(Ref(dt), Base.min, mpi_comm())[]
 
     return dt
 end
@@ -186,7 +190,8 @@ function max_dt(u, t, mesh::ParallelT8codeMesh{2},
                       typeof(constant_speed), typeof(equations), typeof(dg),
                       typeof(cache)},
                 u, t, mesh, constant_speed, equations, dg, cache)
-    dt = MPI.Allreduce!(Ref(dt), min, mpi_comm())[]
+    # Base.min instead of min needed, see comment in src/auxiliary/math.jl
+    dt = MPI.Allreduce!(Ref(dt), Base.min, mpi_comm())[]
 
     return dt
 end
@@ -202,7 +207,8 @@ function max_dt(u, t, mesh::ParallelT8codeMesh{2},
                       typeof(constant_speed), typeof(equations), typeof(dg),
                       typeof(cache)},
                 u, t, mesh, constant_speed, equations, dg, cache)
-    dt = MPI.Allreduce!(Ref(dt), min, mpi_comm())[]
+    # Base.min instead of min needed, see comment in src/auxiliary/math.jl
+    dt = MPI.Allreduce!(Ref(dt), Base.min, mpi_comm())[]
 
     return dt
 end

src/callbacks_step/stepsize_dg3d.jl

@@ -130,7 +130,8 @@ function max_dt(u, t, mesh::ParallelP4estMesh{3},
                       typeof(constant_speed), typeof(equations), typeof(dg),
                       typeof(cache)},
                 u, t, mesh, constant_speed, equations, dg, cache)
-    dt = MPI.Allreduce!(Ref(dt), min, mpi_comm())[]
+    # Base.min instead of min needed, see comment in src/auxiliary/math.jl
+    dt = MPI.Allreduce!(Ref(dt), Base.min, mpi_comm())[]
 
     return dt
 end
@@ -146,7 +147,8 @@ function max_dt(u, t, mesh::ParallelP4estMesh{3},
                       typeof(constant_speed), typeof(equations), typeof(dg),
                       typeof(cache)},
                 u, t, mesh, constant_speed, equations, dg, cache)
-    dt = MPI.Allreduce!(Ref(dt), min, mpi_comm())[]
+    # Base.min instead of min needed, see comment in src/auxiliary/math.jl
+    dt = MPI.Allreduce!(Ref(dt), Base.min, mpi_comm())[]
 
     return dt
 end
@@ -162,7 +164,8 @@ function max_dt(u, t, mesh::ParallelT8codeMesh{3},
                       typeof(constant_speed), typeof(equations), typeof(dg),
                       typeof(cache)},
                 u, t, mesh, constant_speed, equations, dg, cache)
-    dt = MPI.Allreduce!(Ref(dt), min, mpi_comm())[]
+    # Base.min instead of min needed, see comment in src/auxiliary/math.jl
+    dt = MPI.Allreduce!(Ref(dt), Base.min, mpi_comm())[]
 
     return dt
 end
@@ -178,7 +181,8 @@ function max_dt(u, t, mesh::ParallelT8codeMesh{3},
                       typeof(constant_speed), typeof(equations), typeof(dg),
                       typeof(cache)},
                 u, t, mesh, constant_speed, equations, dg, cache)
-    dt = MPI.Allreduce!(Ref(dt), min, mpi_comm())[]
+    # Base.min instead of min needed, see comment in src/auxiliary/math.jl
+    dt = MPI.Allreduce!(Ref(dt), Base.min, mpi_comm())[]
 
     return dt
 end

src/callbacks_step/time_series.jl

@@ -72,7 +72,7 @@ function Base.show(io::IO, ::MIME"text/plain",
             "interval" => time_series_callback.interval,
             "solution_variables" => time_series_callback.solution_variables,
             "output_directory" => time_series_callback.output_directory,
-            "filename" => time_series_callback.filename,
+            "filename" => time_series_callback.filename
         ]
         summary_box(io, "TimeSeriesCallback", setup)
     end

src/callbacks_step/visualization.jl

@@ -50,7 +50,7 @@ function Base.show(io::IO, ::MIME"text/plain",
             "variable names" => visualization_callback.variable_names,
             "show mesh" => visualization_callback.show_mesh,
             "plot creator" => visualization_callback.plot_creator,
-            "plot data creator" => visualization_callback.plot_data_creator,
+            "plot data creator" => visualization_callback.plot_data_creator
         ]
         summary_box(io, "VisualizationCallback", setup)
     end