Merge branch 'main' into AMR_Parabolic_2D3D_Tree

docs/src/parallelization.md

@@ -186,6 +186,11 @@ julia> set_preferences!(
            "libhdf5" => "/path/to/your/libhdf5.so",
            "libhdf5_hl" => "/path/to/your/libhdf5_hl.so", force = true)
 ```
+Alternatively, with HDF5.jl v0.17.1 or higher you can use
+```julia
+julia> using HDF5
+julia> HDF5.API.set_libraries!("/path/to/your/libhdf5.so", "/path/to/your/libhdf5_hl.so")
+```
 For more information see also the
 [documentation of HDF5.jl](https://juliaio.github.io/HDF5.jl/stable/mpi/). In total, you should
 have a file called LocalPreferences.toml in the project directory that contains a section

src/meshes/mesh_io.jl

@@ -263,32 +263,27 @@ function load_mesh_serial(mesh_file::AbstractString; n_cells_max, RealT)
         size = Tuple(size_)
 
         # TODO: `@eval` is evil
-        # A temporary workaround to evaluate the code that defines the domain mapping in a local scope.
-        # This prevents errors when multiple restart elixirs are executed in one session, where one
-        # defines `mapping` as a variable, while the other defines it as a function.
         #
         # This should be replaced with something more robust and secure,
         # see https://github.com/trixi-framework/Trixi.jl/issues/541).
-        expr = Meta.parse(mapping_as_string)
-        if expr.head == :toplevel
-            expr.head = :block
-        end
-
         if ndims == 1
-            mapping = @eval function (xi)
-                $expr
+            mapping = eval(Meta.parse("""function (xi)
+                $mapping_as_string
                 mapping(xi)
             end
+            """))
         elseif ndims == 2
-            mapping = @eval function (xi, eta)
-                $expr
+            mapping = eval(Meta.parse("""function (xi, eta)
+                $mapping_as_string
                 mapping(xi, eta)
             end
+            """))
         else # ndims == 3
-            mapping = @eval function (xi, eta, zeta)
-                $expr
+            mapping = eval(Meta.parse("""function (xi, eta, zeta)
+                $mapping_as_string
                 mapping(xi, eta, zeta)
             end
+            """))
         end
 
         mesh = StructuredMesh(size, mapping; RealT = RealT, unsaved_changes = false,

src/meshes/structured_mesh.jl

@@ -96,13 +96,17 @@ function StructuredMesh(cells_per_dimension, faces::Tuple; RealT = Float64,
 
     # Collect definitions of face functions in one string (separated by semicolons)
     face2substring(face) = code_string(face, ntuple(_ -> Float64, NDIMS - 1))
-    join_semicolon(strings) = join(strings, "; ")
+    join_newline(strings) = join(strings, "\n")
 
-    faces_definition = faces .|> face2substring .|> string |> join_semicolon
+    faces_definition = faces .|> face2substring .|> string |> join_newline
 
     # Include faces definition in `mapping_as_string` to allow for evaluation
     # without knowing the face functions
-    mapping_as_string = "$faces_definition; faces = $(string(faces)); mapping = transfinite_mapping(faces)"
+    mapping_as_string = """
+        $faces_definition
+        faces = $(string(faces))
+        mapping = transfinite_mapping(faces)
+        """
 
     return StructuredMesh(cells_per_dimension, mapping; RealT = RealT,
                           periodicity = periodicity,
@@ -123,13 +127,14 @@ Create a StructuredMesh that represents a uncurved structured mesh with a rectan
 """
 function StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max;
                         periodicity = true)
-    NDIMS = length(cells_per_dimension)
     RealT = promote_type(eltype(coordinates_min), eltype(coordinates_max))
 
     mapping = coordinates2mapping(coordinates_min, coordinates_max)
-    mapping_as_string = "coordinates_min = $coordinates_min; " *
-                        "coordinates_max = $coordinates_max; " *
-                        "mapping = coordinates2mapping(coordinates_min, coordinates_max)"
+    mapping_as_string = """
+        coordinates_min = $coordinates_min
+        coordinates_max = $coordinates_max
+        mapping = coordinates2mapping(coordinates_min, coordinates_max)
+        """
     return StructuredMesh(cells_per_dimension, mapping; RealT = RealT,
                           periodicity = periodicity,
                           mapping_as_string = mapping_as_string)

src/solvers/dgmulti/dg.jl

@@ -302,7 +302,12 @@ function calc_volume_integral!(du, u, mesh::DGMultiMesh,
     @threaded for e in eachelement(mesh, dg, cache)
         flux_values = local_values_threaded[Threads.threadid()]
         for i in eachdim(mesh)
-            flux_values .= flux.(view(u_values, :, e), i, equations)
+            # Here, the broadcasting operation does allocate
+            #flux_values .= flux.(view(u_values, :, e), i, equations)
+            # Use loop instead
+            for j in eachindex(flux_values)
+                flux_values[j] = flux(u_values[j, e], i, equations)
+            end
             for j in eachdim(mesh)
                 apply_to_each_field(mul_by_accum!(weak_differentiation_matrices[j],
                                                   dxidxhatj[i, j][1, e]),
@@ -327,6 +332,7 @@ function calc_volume_integral!(du, u, mesh::DGMultiMesh{NDIMS, <:NonAffine},
     @threaded for e in eachelement(mesh, dg, cache)
         flux_values = cache.flux_threaded[Threads.threadid()]
         for i in eachdim(mesh)
+            # Here, the broadcasting operation does not allocate
             flux_values[i] .= flux.(view(u_values, :, e), i, equations)
         end
 

src/solvers/dgsem_p4est/dg_2d.jl

@@ -501,6 +501,11 @@ function calc_mortar_flux!(surface_flux_values,
         # copying in the correct orientation
         u_buffer = cache.u_threaded[Threads.threadid()]
 
+        # in calc_interface_flux!, the interface flux is computed once over each 
+        # interface using the normal from the "primary" element. The result is then 
+        # passed back to the "secondary" element, flipping the sign to account for the 
+        # change in the normal direction. For mortars, this sign flip occurs in 
+        # "mortar_fluxes_to_elements!" instead.
         mortar_fluxes_to_elements!(surface_flux_values,
                                    mesh, equations, mortar_l2, dg, cache,
                                    mortar, fstar, u_buffer)

src/solvers/dgsem_p4est/dg_3d.jl

@@ -580,6 +580,11 @@ function calc_mortar_flux!(surface_flux_values,
         # copying in the correct orientation
         u_buffer = cache.u_threaded[Threads.threadid()]
 
+        # in calc_interface_flux!, the interface flux is computed once over each 
+        # interface using the normal from the "primary" element. The result is then 
+        # passed back to the "secondary" element, flipping the sign to account for the 
+        # change in the normal direction. For mortars, this sign flip occurs in 
+        # "mortar_fluxes_to_elements!" instead.
         mortar_fluxes_to_elements!(surface_flux_values,
                                    mesh, equations, mortar_l2, dg, cache,
                                    mortar, fstar, u_buffer, fstar_tmp)