diff --git a/docs/src/parallelization.md b/docs/src/parallelization.md index e55471bb256..610aa2cbd95 100644 --- a/docs/src/parallelization.md +++ b/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 diff --git a/src/meshes/mesh_io.jl b/src/meshes/mesh_io.jl index b9895e7d454..92e38ce1bf3 100644 --- a/src/meshes/mesh_io.jl +++ b/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, diff --git a/src/meshes/structured_mesh.jl b/src/meshes/structured_mesh.jl index df067db833d..553aabbbc20 100644 --- a/src/meshes/structured_mesh.jl +++ b/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) diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl index bc76aa1a9d2..182a486dce5 100644 --- a/src/solvers/dgmulti/dg.jl +++ b/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 diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index a665aa4b19d..36624f2ce8a 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/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) diff --git a/src/solvers/dgsem_p4est/dg_3d.jl b/src/solvers/dgsem_p4est/dg_3d.jl index dc69329474f..4c0845ba9af 100644 --- a/src/solvers/dgsem_p4est/dg_3d.jl +++ b/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)