Merge branch 'main' into mutable-2

.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,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.8.9-DEV"
+version = "0.8.10-DEV"
 
 [deps]
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
@@ -61,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/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_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_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/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

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/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

src/meshes/p4est_mesh.jl

@@ -6,12 +6,23 @@
 #! format: noindent
 
 """
-    P4estMesh{NDIMS} <: AbstractMesh{NDIMS}
+    P4estMesh{NDIMS, NDIMS_AMBIENT} <: AbstractMesh{NDIMS}
 
 An unstructured curved mesh based on trees that uses the C library `p4est`
 to manage trees and mesh refinement.
+
+The parameter `NDIMS` denotes the dimension of the spatial domain or manifold represented
+by the mesh itself, while `NDIMS_AMBIENT` denotes the dimension of the ambient space in
+which the mesh is embedded. For example, the type `P4estMesh{3, 3}` corresponds to a
+standard mesh for a three-dimensional volume, whereas `P4estMesh{2, 3}` corresponds to a
+mesh for a two-dimensional surface or shell in three-dimensional space.
+
+!!! warning "Experimental implementation"
+    The use of `NDIMS != NDIMS_AMBIENT` is an experimental feature and may change in future
+    releases.
 """
-mutable struct P4estMesh{NDIMS, RealT <: Real, IsParallel, P, Ghost, NDIMSP2, NNODES} <:
+mutable struct P4estMesh{NDIMS, NDIMS_AMBIENT, RealT <: Real, IsParallel, P, Ghost,
+                         NDIMSP2, NNODES} <:
                AbstractMesh{NDIMS}
     p4est       :: P # Either PointerWrapper{p4est_t} or PointerWrapper{p8est_t}
     is_parallel :: IsParallel
@@ -48,7 +59,14 @@ mutable struct P4estMesh{NDIMS, RealT <: Real, IsParallel, P, Ghost, NDIMSP2, NN
         ghost = ghost_new_p4est(p4est)
         ghost_pw = PointerWrapper(ghost)
 
-        mesh = new{NDIMS, eltype(tree_node_coordinates), typeof(is_parallel),
+        # To enable the treatment of a manifold of dimension NDIMS embedded within an
+        # ambient space of dimension NDIMS_AMBIENT, we store both as type parameters and
+        # allow them to differ in the general case. This functionality is used for
+        # constructing discretizations on spherical shell domains for applications in
+        # global atmospheric modelling. The ambient dimension NDIMS_AMBIENT is therefore 
+        # set here in the inner constructor to size(tree_node_coordinates, 1).
+        mesh = new{NDIMS, size(tree_node_coordinates, 1),
+                   eltype(tree_node_coordinates), typeof(is_parallel),
                    typeof(p4est_pw), typeof(ghost_pw), NDIMS + 2, length(nodes)}(p4est_pw,
                                                                                  is_parallel,
                                                                                  ghost_pw,
@@ -66,8 +84,8 @@ mutable struct P4estMesh{NDIMS, RealT <: Real, IsParallel, P, Ghost, NDIMSP2, NN
     end
 end
 
-const SerialP4estMesh{NDIMS} = P4estMesh{NDIMS, <:Real, <:False}
-const ParallelP4estMesh{NDIMS} = P4estMesh{NDIMS, <:Real, <:True}
+const SerialP4estMesh{NDIMS} = P4estMesh{NDIMS, <:Any, <:Real, <:False}
+const ParallelP4estMesh{NDIMS} = P4estMesh{NDIMS, <:Any, <:Real, <:True}
 
 @inline mpi_parallel(mesh::SerialP4estMesh) = False()
 @inline mpi_parallel(mesh::ParallelP4estMesh) = True()
@@ -87,7 +105,8 @@ function destroy_mesh(mesh::P4estMesh{3})
 end
 
 @inline Base.ndims(::P4estMesh{NDIMS}) where {NDIMS} = NDIMS
-@inline Base.real(::P4estMesh{NDIMS, RealT}) where {NDIMS, RealT} = RealT
+@inline Base.real(::P4estMesh{NDIMS, NDIMS_AMBIENT, RealT}) where {NDIMS, NDIMS_AMBIENT, RealT} = RealT
+@inline ndims_ambient(::P4estMesh{NDIMS, NDIMS_AMBIENT}) where {NDIMS, NDIMS_AMBIENT} = NDIMS_AMBIENT
 
 @inline function ntrees(mesh::P4estMesh)
     return mesh.p4est.trees.elem_count[]
@@ -97,7 +116,8 @@ end
 @inline ncellsglobal(mesh::P4estMesh) = Int(mesh.p4est.global_num_quadrants[])
 
 function Base.show(io::IO, mesh::P4estMesh)
-    print(io, "P4estMesh{", ndims(mesh), ", ", real(mesh), "}")
+    print(io, "P4estMesh{", ndims(mesh), ", ", ndims_ambient(mesh), ", ", real(mesh),
+          "}")
 end
 
 function Base.show(io::IO, ::MIME"text/plain", mesh::P4estMesh)
@@ -107,11 +127,12 @@ function Base.show(io::IO, ::MIME"text/plain", mesh::P4estMesh)
         setup = [
             "#trees" => ntrees(mesh),
             "current #cells" => ncellsglobal(mesh),
-            "polydeg" => length(mesh.nodes) - 1,
+            "polydeg" => length(mesh.nodes) - 1
         ]
         summary_box(io,
-                    "P4estMesh{" * string(ndims(mesh)) * ", " * string(real(mesh)) *
-                    "}", setup)
+                    "P4estMesh{" * string(ndims(mesh)) * ", " *
+                    string(ndims_ambient(mesh)) *
+                    ", " * string(real(mesh)) * "}", setup)
     end
 end
 

src/meshes/t8code_mesh.jl

@@ -93,7 +93,7 @@ function Base.show(io::IO, ::MIME"text/plain", mesh::T8codeMesh)
         setup = [
             "#trees" => ntrees(mesh),
             "current #cells" => ncellsglobal(mesh),
-            "polydeg" => length(mesh.nodes) - 1,
+            "polydeg" => length(mesh.nodes) - 1
         ]
         summary_box(io,
                     "T8codeMesh{" * string(ndims(mesh)) * ", " * string(real(mesh)) * "}",
@@ -982,7 +982,7 @@ function fill_mesh_info!(mesh::T8codeMesh, interfaces, mortars, boundaries,
         [1, 2, 0, 0, 3, 4, 0, 0], # 2
         [0, 0, 1, 2, 0, 0, 3, 4], # 3
         [1, 2, 3, 4, 0, 0, 0, 0], # 4
-        [0, 0, 0, 0, 1, 2, 3, 4], # 5
+        [0, 0, 0, 0, 1, 2, 3, 4] # 5
     ]
 
     # Helper variables to compute unique global MPI interface/mortar ids.
@@ -1235,18 +1235,18 @@ function fill_mesh_info!(mesh::T8codeMesh, interfaces, mortars, boundaries,
                                 global_mortar_id_to_local[global_mortar_id] = local_mpi_mortar_id
 
                                 mpi_mesh_info.mpi_mortars.local_neighbor_ids[local_mpi_mortar_id] = [
-                                    current_index + 1,
+                                    current_index + 1
                                 ]
                                 mpi_mesh_info.mpi_mortars.local_neighbor_positions[local_mpi_mortar_id] = [
-                                    map_iface_to_ichild_to_position[iface + 1][t8_element_child_id(eclass_scheme, element) + 1],
+                                    map_iface_to_ichild_to_position[iface + 1][t8_element_child_id(eclass_scheme, element) + 1]
                                 ]
                                 init_mortar_node_indices!(mpi_mesh_info.mpi_mortars,
                                                           (iface, dual_faces[1]),
                                                           orientation, local_mpi_mortar_id)
 
                                 neighbor_ranks = [
                                     remotes[findlast(ghost_remote_first_elem .<=
-                                                     neighbor_ielements[1])],
+                                                     neighbor_ielements[1])]
                                 ]
                                 mpi_mesh_info.neighbor_ranks_mortar[local_mpi_mortar_id] = neighbor_ranks
 

src/meshes/tree_mesh.jl

@@ -200,7 +200,7 @@ function Base.show(io::IO, ::MIME"text/plain",
             "periodicity" => mesh.tree.periodicity,
             "current #cells" => mesh.tree.length,
             "#leaf-cells" => count_leaf_cells(mesh.tree),
-            "maximum #cells" => mesh.tree.capacity,
+            "maximum #cells" => mesh.tree.capacity
         ]
         summary_box(io, "TreeMesh{" * string(NDIMS) * ", " * string(TreeType) * "}",
                     setup)

src/semidiscretization/semidiscretization_euler_gravity.jl

@@ -59,7 +59,7 @@ function Base.show(io::IO, ::MIME"text/plain", parameters::ParametersEulerGravit
             "gravitational constant (G)" => parameters.gravitational_constant,
             "CFL (gravity)" => parameters.cfl,
             "max. #iterations" => parameters.n_iterations_max,
-            "time integrator" => parameters.timestep_gravity,
+            "time integrator" => parameters.timestep_gravity
         ]
         summary_box(io, "ParametersEulerGravity", setup)
     end