diff --git a/examples/p4est_2d_dgsem/elixir_advection_flag_gauss_legendre.jl b/examples/p4est_2d_dgsem/elixir_advection_flag_gauss_legendre.jl new file mode 100644 index 00000000000..815317ae26c --- /dev/null +++ b/examples/p4est_2d_dgsem/elixir_advection_flag_gauss_legendre.jl @@ -0,0 +1,49 @@ +using OrdinaryDiffEqLowStorageRK +using Trixi + +############################################################################### +# semidiscretization of the linear advection equation + +advection_velocity = (0.2, -0.7) +equations = LinearScalarAdvectionEquation2D(advection_velocity) + +solver = DGSEM(polydeg = 4, surface_flux = flux_lax_friedrichs, + basis_type = GaussLegendreBasis) + +# Deformed rectangle that looks like a waving flag, +# lower and upper faces are sinus curves, left and right are vertical lines. +f1(s) = SVector(-1.0, s - 1.0) +f2(s) = SVector(1.0, s + 1.0) +f3(s) = SVector(s, -1.0 + sinpi(0.5 * s)) +f4(s) = SVector(s, 1.0 + sinpi(0.5 * s)) + +# Create P4estMesh with 3 x 2 trees and 6 x 4 elements, +# approximate the geometry with a smaller polydeg for testing. +trees_per_dimension = (3, 2) +mesh = P4estMesh(trees_per_dimension, polydeg = 3, + faces = (f1, f2, f3, f4), + periodicity = true, + initial_refinement_level = 1) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, + solver; + boundary_conditions = boundary_condition_periodic) + +############################################################################### +# ODE solvers, callbacks etc. + +ode = semidiscretize(semi, (0.0, 0.2)) + +summary_callback = SummaryCallback() +analysis_callback = AnalysisCallback(semi, interval = 100) +stepsize_callback = StepsizeCallback(cfl = 0.5) + +callbacks = CallbackSet(summary_callback, analysis_callback, + stepsize_callback) + +############################################################################### +# 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 + ode_default_options()..., callback = callbacks); diff --git a/src/meshes/p4est_mesh_view.jl b/src/meshes/p4est_mesh_view.jl index 4f0bce4c69b..3d23a49f2bf 100644 --- a/src/meshes/p4est_mesh_view.jl +++ b/src/meshes/p4est_mesh_view.jl @@ -97,6 +97,10 @@ function extract_interfaces(mesh::P4estMeshView, interfaces_parent) # Copy relevant entries from parent mesh @views interfaces.u .= interfaces_parent.u[.., mask] + if interfaces.normal_directions !== nothing # Needed for Gauss-Legendre basis + @views interfaces.normal_directions .= interfaces_parent.normal_directions[.., + mask] + end @views interfaces.node_indices .= interfaces_parent.node_indices[.., mask] @views neighbor_ids = interfaces_parent.neighbor_ids[.., mask] diff --git a/src/solvers/dgsem_p4est/containers.jl b/src/solvers/dgsem_p4est/containers.jl index 2dd360b50cc..c86c2ea3d24 100644 --- a/src/solvers/dgsem_p4est/containers.jl +++ b/src/solvers/dgsem_p4est/containers.jl @@ -190,11 +190,8 @@ function Adapt.adapt_structure(to, size(elements.surface_flux_values)) new_type_params = (NDIMS, - RealT, - uEltype, - NDIMS + 1, - NDIMS + 2, - NDIMS + 3, + RealT, uEltype, + NDIMS + 1, NDIMS + 2, NDIMS + 3, typeof(inverse_jacobian), # ArrayRealTNDIMSP1 typeof(node_coordinates), # ArrayRealTNDIMSP2 typeof(jacobian_matrix), # ArrayRealTNDIMSP3 @@ -213,38 +210,48 @@ function Adapt.adapt_structure(to, _surface_flux_values) end -mutable struct P4estInterfaceContainer{NDIMS, uEltype <: Real, NDIMSP2, +mutable struct P4estInterfaceContainer{NDIMS, RealT <: Real, uEltype <: Real, + NDIMSP1, NDIMSP2, uArray <: DenseArray{uEltype, NDIMSP2}, + NormalArray <: + Union{DenseArray{RealT, NDIMSP1}, Nothing}, IdsMatrix <: DenseMatrix{Int}, IndicesMatrix <: DenseMatrix{NTuple{NDIMS, Symbol}}, uVector <: DenseVector{uEltype}, + NormalVector <: + Union{DenseVector{RealT}, Nothing}, IdsVector <: DenseVector{Int}, IndicesVector <: DenseVector{NTuple{NDIMS, Symbol}}} <: AbstractInterfaceContainer - u::uArray # [primary/secondary, variable, i, j, interface] - neighbor_ids::IdsMatrix # [primary/secondary, interface] - node_indices::IndicesMatrix # [primary/secondary, interface] + u::uArray # [primary/secondary, variable, i, j, interface] + normal_directions::NormalArray # [dimension, i, j, interface] + neighbor_ids::IdsMatrix # [primary/secondary, interface] + node_indices::IndicesMatrix # [primary/secondary, interface] # internal `resize!`able storage _u::uVector + _normal_directions::NormalVector _neighbor_ids::IdsVector _node_indices::IndicesVector end +# `trivial` means that the interface normals can be taken from +# the outer/surface element nodes +@inline trivial_interface_normals(::LobattoLegendreBasis) = true +# For Gauss-Legendre basis, the interface normals need to be interpolated, +# analogous to the surface values. +@inline trivial_interface_normals(::GaussLegendreBasis) = false + @inline function ninterfaces(interfaces::P4estInterfaceContainer) return size(interfaces.neighbor_ids, 2) end @inline Base.ndims(::P4estInterfaceContainer{NDIMS}) where {NDIMS} = NDIMS -@inline function Base.eltype(::P4estInterfaceContainer{NDIMS, uEltype}) where {NDIMS, - uEltype} - return uEltype -end # See explanation of Base.resize! for the element container function Base.resize!(interfaces::P4estInterfaceContainer, capacity) - @unpack _u, _neighbor_ids, _node_indices = interfaces + @unpack _u, _normal_directions, _neighbor_ids, _node_indices = interfaces n_dims = ndims(interfaces) n_nodes = size(interfaces.u, 3) @@ -256,6 +263,18 @@ function Base.resize!(interfaces::P4estInterfaceContainer, capacity) (2, n_variables, ntuple(_ -> n_nodes, n_dims - 1)..., capacity)) + if _normal_directions === nothing # trivial interface normals, e.g. for Lobatto-Legendre basis + interfaces.normal_directions = nothing + else # non-trivial interface normals, e.g. for Gauss-Legendre basis + resize!(_normal_directions, n_dims * n_nodes^(n_dims - 1) * capacity) + interfaces.normal_directions = unsafe_wrap(ArrayType, + pointer(_normal_directions), + (n_dims, + ntuple(_ -> n_nodes, + n_dims - 1)..., + capacity)) + end + resize!(_neighbor_ids, 2 * capacity) interfaces.neighbor_ids = unsafe_wrap(ArrayType, pointer(_neighbor_ids), (2, capacity)) @@ -271,6 +290,7 @@ end function init_interfaces(mesh::Union{P4estMesh, P4estMeshView, T8codeMesh}, equations, basis, elements) NDIMS = ndims(elements) + RealT = real(mesh) uEltype = eltype(elements) # Initialize container @@ -283,59 +303,94 @@ function init_interfaces(mesh::Union{P4estMesh, P4estMeshView, T8codeMesh}, equa (2, nvariables(equations), ntuple(_ -> nnodes(basis), NDIMS - 1)..., n_interfaces)) + if !trivial_interface_normals(basis) + _normal_directions = Vector{RealT}(undef, + NDIMS * nnodes(basis)^(NDIMS - 1) * + n_interfaces) + normal_directions = unsafe_wrap(Array, pointer(_normal_directions), + (NDIMS, + ntuple(_ -> nnodes(basis), NDIMS - 1)..., + n_interfaces)) + else + _normal_directions = nothing + normal_directions = nothing + end + _neighbor_ids = Vector{Int}(undef, 2 * n_interfaces) neighbor_ids = unsafe_wrap(Array, pointer(_neighbor_ids), (2, n_interfaces)) _node_indices = Vector{NTuple{NDIMS, Symbol}}(undef, 2 * n_interfaces) node_indices = unsafe_wrap(Array, pointer(_node_indices), (2, n_interfaces)) - interfaces = P4estInterfaceContainer{NDIMS, uEltype, NDIMS + 2, - typeof(u), typeof(neighbor_ids), - typeof(node_indices), typeof(_u), + interfaces = P4estInterfaceContainer{NDIMS, RealT, uEltype, + NDIMS + 1, NDIMS + 2, + typeof(u), typeof(normal_directions), + typeof(neighbor_ids), typeof(node_indices), + typeof(_u), typeof(_normal_directions), typeof(_neighbor_ids), typeof(_node_indices)}(u, + normal_directions, neighbor_ids, node_indices, _u, + _normal_directions, _neighbor_ids, _node_indices) - init_interfaces!(interfaces, mesh) + init_interfaces!(interfaces, elements, mesh, basis) return interfaces end -function init_interfaces!(interfaces, mesh::Union{P4estMesh, P4estMeshView}) +function init_interfaces!(interfaces, elements, + mesh::Union{P4estMesh, P4estMeshView}, basis) init_surfaces!(interfaces, nothing, nothing, mesh) + init_normal_directions!(interfaces, basis, elements) return interfaces end function Adapt.parent_type(::Type{<:P4estInterfaceContainer{<:Any, <:Any, <:Any, + <:Any, <:Any, ArrayT}}) where {ArrayT} return ArrayT end # Manual adapt_structure since we have aliasing memory -function Adapt.adapt_structure(to, interfaces::P4estInterfaceContainer) +function Adapt.adapt_structure(to, + interfaces::P4estInterfaceContainer{NDIMS, + RealT, + uEltype}) where { + NDIMS, + uEltype, + RealT + } # Adapt underlying storage _u = adapt(to, interfaces._u) + _normal_directions = interfaces._normal_directions === nothing ? nothing : + adapt(to, interfaces._normal_directions) _neighbor_ids = adapt(to, interfaces._neighbor_ids) _node_indices = adapt(to, interfaces._node_indices) # Wrap arrays again u = unsafe_wrap_or_alloc(to, _u, size(interfaces.u)) + normal_directions = _normal_directions === nothing ? nothing : # Check if interface normals are trivial + unsafe_wrap_or_alloc(to, _normal_directions, + size(interfaces.normal_directions)) neighbor_ids = unsafe_wrap_or_alloc(to, _neighbor_ids, size(interfaces.neighbor_ids)) node_indices = unsafe_wrap_or_alloc(to, _node_indices, size(interfaces.node_indices)) - NDIMS = ndims(interfaces) new_type_params = (NDIMS, - eltype(_u), - NDIMS + 2, - typeof(u), typeof(neighbor_ids), typeof(node_indices), - typeof(_u), typeof(_neighbor_ids), typeof(_node_indices)) - return P4estInterfaceContainer{new_type_params...}(u, neighbor_ids, node_indices, - _u, _neighbor_ids, _node_indices) + RealT, eltype(_u), + NDIMS + 1, NDIMS + 2, + typeof(u), typeof(normal_directions), + typeof(neighbor_ids), typeof(node_indices), + typeof(_u), typeof(_normal_directions), + typeof(_neighbor_ids), typeof(_node_indices)) + return P4estInterfaceContainer{new_type_params...}(u, normal_directions, + neighbor_ids, node_indices, + _u, _normal_directions, + _neighbor_ids, _node_indices) end mutable struct P4estBoundaryContainer{NDIMS, uEltype <: Real, NDIMSP1, @@ -618,8 +673,7 @@ function Adapt.adapt_structure(to, mortars::P4estMortarContainer) NDIMS = ndims(mortars) new_type_params = (NDIMS, eltype(_u), - NDIMS + 1, - NDIMS + 3, + NDIMS + 1, NDIMS + 3, typeof(u), typeof(neighbor_ids), typeof(node_indices), typeof(_u), typeof(_neighbor_ids), typeof(_node_indices)) return P4estMortarContainer{new_type_params...}(u, neighbor_ids, node_indices, @@ -655,7 +709,12 @@ function reinitialize_containers!(mesh::P4estMesh, equations, dg::DGSEM, cache) # re-initialize containers together to reduce # the number of iterations over the mesh in `p4est` - return init_surfaces!(interfaces, mortars, boundaries, mesh) + init_surfaces!(interfaces, mortars, boundaries, mesh) + + # init_normal_directions! requires that `node_indices` have been initialized + init_normal_directions!(interfaces, dg.basis, elements) + + return nothing end # A helper struct used in initialization methods below diff --git a/src/solvers/dgsem_p4est/containers_2d.jl b/src/solvers/dgsem_p4est/containers_2d.jl index 0cd04609d4a..a0f7accf6c3 100644 --- a/src/solvers/dgsem_p4est/containers_2d.jl +++ b/src/solvers/dgsem_p4est/containers_2d.jl @@ -8,7 +8,7 @@ # Initialize data structures in element container function init_elements!(elements, mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, - basis::LobattoLegendreBasis) + basis::AbstractBasisSBP) @unpack node_coordinates, jacobian_matrix, contravariant_vectors, inverse_jacobian = elements @@ -29,7 +29,7 @@ end function calc_node_coordinates!(node_coordinates, mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, - basis::LobattoLegendreBasis) + basis::AbstractBasisSBP) # Hanging nodes will cause holes in the mesh if its polydeg is higher # than the polydeg of the solver. @assert length(basis.nodes)>=length(mesh.nodes) "The solver can't have a lower polydeg than the mesh" @@ -87,6 +87,89 @@ function calc_node_coordinates!(node_coordinates, return node_coordinates end +# For Gauss-Lobatto-Legendre (GLL) nodes, interface normals are computed on-the-fly +# during flux evaluation and do not need dedicated storage in the interface container. +function init_normal_directions!(interfaces::P4estInterfaceContainer{2}, + basis::LobattoLegendreBasis, elements) + return nothing +end + +# For Gauss-Legendre (GL) nodes, the interface normals are +# computed from interpolation of the volume node normals to the surface nodes. +function init_normal_directions!(interfaces::P4estInterfaceContainer{2}, + basis::GaussLegendreBasis, elements) + @unpack neighbor_ids, node_indices, normal_directions = interfaces + @unpack contravariant_vectors = elements + @unpack boundary_interpolation = basis + index_range = eachnode(basis) + + for interface in axes(neighbor_ids, 2) + primary_element = neighbor_ids[1, interface] + primary_indices = node_indices[1, interface] + primary_direction = indices2direction(primary_indices) + + i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1], + index_range) + j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2], + index_range) + + # The index direction is identified based on `{i,j}_{primary, secondary}_step`. + # For step = 0, the direction identified by this index is normal to the face. + # For step != 0 (1 or -1), the direction identified by this index is tangential to the face. + + i_primary = i_primary_start + j_primary = j_primary_start + + if i_primary_step == 0 + # i is the normal direction (constant), j varies along the surface + # => Interpolate in first/normal direction + # Interpolation side is governed by element orientation + side = interpolation_side(primary_indices[1]) + for i in index_range + normal_1 = zero(eltype(normal_directions)) + normal_2 = zero(eltype(normal_directions)) + for ii in index_range + factor = boundary_interpolation[ii, side] + # Retrieve normal directions at element/volume nodes + normal_direction = get_normal_direction(primary_direction, + contravariant_vectors, + ii, j_primary, + primary_element) + normal_1 = normal_1 + normal_direction[1] * factor + normal_2 = normal_2 + normal_direction[2] * factor + end + normal_directions[1, i, interface] = normal_1 + normal_directions[2, i, interface] = normal_2 + j_primary += j_primary_step # incrementing j_primary suffices (i_primary_step = 0) + end + else # j_primary_step == 0 + # j is the normal direction (constant), i varies along the surface + # => Interpolate in second/normal direction + # Interpolation side is governed by element orientation + side = interpolation_side(primary_indices[2]) + for i in index_range + normal_1 = zero(eltype(normal_directions)) + normal_2 = zero(eltype(normal_directions)) + for jj in index_range + factor = boundary_interpolation[jj, side] + # Retrieve normal directions at element/volume nodes + normal_direction = get_normal_direction(primary_direction, + contravariant_vectors, + i_primary, jj, + primary_element) + normal_1 = normal_1 + normal_direction[1] * factor + normal_2 = normal_2 + normal_direction[2] * factor + end + normal_directions[1, i, interface] = normal_1 + normal_directions[2, i, interface] = normal_2 + i_primary += i_primary_step # incrementing i_primary suffices (j_primary_step = 0) + end + end + end + + return nothing +end + # Initialize node_indices of interface container @inline function init_interface_node_indices!(interfaces::P4estInterfaceContainer{2}, faces, orientation, interface_id) diff --git a/src/solvers/dgsem_p4est/containers_3d.jl b/src/solvers/dgsem_p4est/containers_3d.jl index 9d343c76f02..c0e295153c7 100644 --- a/src/solvers/dgsem_p4est/containers_3d.jl +++ b/src/solvers/dgsem_p4est/containers_3d.jl @@ -76,6 +76,12 @@ function calc_node_coordinates!(node_coordinates, return node_coordinates end +# Not yet implemented and needed for 3D +function init_normal_directions!(interfaces::P4estInterfaceContainer{3}, + basis::LobattoLegendreBasis, elements) + return nothing +end + # Initialize node_indices of interface container @inline function init_interface_node_indices!(interfaces::P4estInterfaceContainer{3}, faces, orientation, interface_id) diff --git a/src/solvers/dgsem_p4est/containers_parallel.jl b/src/solvers/dgsem_p4est/containers_parallel.jl index 82cf2bcbf56..f4cb9a31ede 100644 --- a/src/solvers/dgsem_p4est/containers_parallel.jl +++ b/src/solvers/dgsem_p4est/containers_parallel.jl @@ -252,8 +252,10 @@ end # Overload init! function for regular interfaces, regular mortars and boundaries since they must # call the appropriate init_surfaces! function for parallel p4est meshes -function init_interfaces!(interfaces, mesh::P4estMeshParallel) +function init_interfaces!(interfaces, elements, + mesh::P4estMeshParallel, basis) init_surfaces!(interfaces, nothing, nothing, nothing, nothing, mesh) + init_normal_directions!(interfaces, basis, elements) return interfaces end @@ -312,6 +314,9 @@ function reinitialize_containers!(mesh::P4estMeshParallel, equations, dg::DGSEM, # the number of iterations over the mesh in p4est init_surfaces!(interfaces, mortars, boundaries, mpi_interfaces, mpi_mortars, mesh) + # init_normal_directions! requires that `node_indices` have been initialized + init_normal_directions!(interfaces, dg.basis, elements) + # re-initialize MPI cache @unpack mpi_cache = cache init_mpi_cache!(mpi_cache, mesh, mpi_interfaces, mpi_mortars, @@ -320,7 +325,9 @@ function reinitialize_containers!(mesh::P4estMeshParallel, equations, dg::DGSEM, # re-initialize and distribute normal directions of MPI mortars; requires MPI communication, so # the MPI cache must be re-initialized before init_normal_directions!(mpi_mortars, dg.basis, elements) - return exchange_normal_directions!(mpi_mortars, mpi_cache, mesh, nnodes(dg)) + exchange_normal_directions!(mpi_mortars, mpi_cache, mesh, nnodes(dg)) + + return nothing end # A helper struct used in initialization methods below diff --git a/src/solvers/dgsem_p4est/containers_parallel_2d.jl b/src/solvers/dgsem_p4est/containers_parallel_2d.jl index d531d33821b..003f97aab5f 100644 --- a/src/solvers/dgsem_p4est/containers_parallel_2d.jl +++ b/src/solvers/dgsem_p4est/containers_parallel_2d.jl @@ -39,8 +39,8 @@ end # Normal directions of small element surfaces are needed to calculate the mortar fluxes. Initialize # them for locally available small elements. -function init_normal_directions!(mpi_mortars::P4estMPIMortarContainer{2}, basis, - elements) +function init_normal_directions!(mpi_mortars::P4estMPIMortarContainer{2}, + basis, elements) @unpack local_neighbor_ids, local_neighbor_positions, node_indices = mpi_mortars @unpack contravariant_vectors = elements index_range = eachnode(basis) @@ -77,5 +77,7 @@ function init_normal_directions!(mpi_mortars::P4estMPIMortarContainer{2}, basis, end end end + + return nothing end end # muladd diff --git a/src/solvers/dgsem_p4est/containers_parallel_3d.jl b/src/solvers/dgsem_p4est/containers_parallel_3d.jl index 56f0a543b97..8a2c2939679 100644 --- a/src/solvers/dgsem_p4est/containers_parallel_3d.jl +++ b/src/solvers/dgsem_p4est/containers_parallel_3d.jl @@ -97,8 +97,8 @@ end # Normal directions of small element surfaces are needed to calculate the mortar fluxes. Initialize # them for locally available small elements. -function init_normal_directions!(mpi_mortars::P4estMPIMortarContainer{3}, basis, - elements) +function init_normal_directions!(mpi_mortars::P4estMPIMortarContainer{3}, + basis, elements) @unpack local_neighbor_ids, local_neighbor_positions, node_indices = mpi_mortars @unpack contravariant_vectors = elements index_range = eachnode(basis) @@ -145,5 +145,7 @@ function init_normal_directions!(mpi_mortars::P4estMPIMortarContainer{3}, basis, end end end + + return nothing end end # muladd diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index c00bf87358e..f60ece112c1 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -62,18 +62,24 @@ end end end +# Infer interpolation side, i.e., left (1) or right (2) for an element. +# Required for boundary interpolation with Gauss-Legendre nodes. +@inline function interpolation_side(index) + return (index === :begin) ? 1 : 2 +end + function prolong2interfaces!(backend::Nothing, cache, u, mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, - equations, dg::DG) + equations, dg::DGSEM{<:LobattoLegendreBasis}) @unpack interfaces = cache @unpack neighbor_ids, node_indices = cache.interfaces index_range = eachnode(dg) @threaded for interface in eachinterface(dg, cache) - prolong2interfaces_per_interface!(interfaces.u, u, interface, typeof(mesh), - equations, neighbor_ids, node_indices, - index_range) + prolong2interfaces_per_interface!(interfaces.u, u, interface, + typeof(mesh), equations, + neighbor_ids, node_indices, index_range) end return nothing end @@ -81,7 +87,7 @@ end function prolong2interfaces!(backend::Backend, cache, u, mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, - equations, dg::DG) + equations, dg::DGSEM{<:LobattoLegendreBasis}) @unpack interfaces = cache ninterfaces(interfaces) == 0 && return nothing @unpack neighbor_ids, node_indices = cache.interfaces @@ -104,13 +110,13 @@ end neighbor_ids, node_indices, index_range) end +# Version for Gauss-Lobatto-Legendre @inline function prolong2interfaces_per_interface!(interfaces_u, u, interface, ::Type{<:Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}}, equations, neighbor_ids, - node_indices, - index_range) + node_indices, index_range) primary_element = neighbor_ids[1, interface] primary_indices = node_indices[1, interface] @@ -154,12 +160,139 @@ end return nothing end +function prolong2interfaces!(backend::Nothing, cache, u, + mesh::Union{P4estMesh{2}, P4estMeshView{2}}, + equations, dg::DGSEM{<:GaussLegendreBasis}) + @unpack interfaces = cache + @unpack neighbor_ids, node_indices = cache.interfaces + @unpack boundary_interpolation = dg.basis + index_range = eachnode(dg) + + @threaded for interface in eachinterface(dg, cache) + prolong2interfaces_per_interface!(interfaces.u, u, interface, + typeof(mesh), equations, + neighbor_ids, node_indices, index_range, + boundary_interpolation) + end + + return nothing +end + +# Version for Gauss-Legendre, which requires passing in the boundary interpolation matrix +@inline function prolong2interfaces_per_interface!(interfaces_u, u, interface, + ::Type{<:Union{P4estMesh{2}, + P4estMeshView{2}}}, + equations, neighbor_ids, + node_indices, index_range, + boundary_interpolation) + # Interpolate solution data from the primary element to the interface. + primary_element = neighbor_ids[1, interface] + primary_indices = node_indices[1, interface] + + i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1], + index_range) + j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2], + index_range) + # The index direction is identified based on `{i,j}_{primary, secondary}_step`. + # For step = 0, the direction identified by this index is normal to the face. + # For step != 0 (1 or -1), the direction identified by this index is tangential to the face. + + # Note that in the current implementation, the interface will be + # "aligned at the primary element", i.e., the index of the primary side + # will always run forwards. + + i_primary = i_primary_start + j_primary = j_primary_start + + if i_primary_step == 0 + # i is the normal direction (constant), j varies along the surface + # => Interpolate in first/normal direction + # Interpolation side is governed by element orientation + side = interpolation_side(primary_indices[1]) + for i in index_range + for v in eachvariable(equations) + u_primary = zero(eltype(interfaces_u)) + for ii in index_range + u_primary = (u_primary + + u[v, ii, j_primary, primary_element] * + boundary_interpolation[ii, side]) + end + interfaces_u[1, v, i, interface] = u_primary + end + j_primary += j_primary_step # incrementing j_primary suffices + end + else # j_primary_step == 0 + # j is the normal direction (constant), i varies along the surface + # => Interpolate in second/normal direction + # Interpolation side is governed by element orientation + side = interpolation_side(primary_indices[2]) + for i in index_range + for v in eachvariable(equations) + u_primary = zero(eltype(interfaces_u)) + for jj in index_range + u_primary = (u_primary + + u[v, i_primary, jj, primary_element] * + boundary_interpolation[jj, side]) + end + interfaces_u[1, v, i, interface] = u_primary + end + i_primary += i_primary_step # incrementing i_primary suffices + end + end + + # Interpolate solution data from the secondary element to the interface. + secondary_element = neighbor_ids[2, interface] + secondary_indices = node_indices[2, interface] + + i_secondary_start, i_secondary_step = index_to_start_step_2d(secondary_indices[1], + index_range) + j_secondary_start, j_secondary_step = index_to_start_step_2d(secondary_indices[2], + index_range) + + i_secondary = i_secondary_start + j_secondary = j_secondary_start + if i_secondary_step == 0 + side = interpolation_side(secondary_indices[1]) + for i in index_range + for v in eachvariable(equations) + u_secondary = zero(eltype(interfaces_u)) + for ii in index_range + u_secondary = (u_secondary + + u[v, ii, j_secondary, secondary_element] * + boundary_interpolation[ii, side]) + end + interfaces_u[2, v, i, interface] = u_secondary + end + j_secondary += j_secondary_step # incrementing j_secondary suffices + end + else # j_secondary_step == 0 + side = interpolation_side(secondary_indices[2]) + for i in index_range + for v in eachvariable(equations) + u_secondary = zero(eltype(interfaces_u)) + for jj in index_range + u_secondary = (u_secondary + + u[v, i_secondary, jj, secondary_element] * + boundary_interpolation[jj, side]) + end + interfaces_u[2, v, i, interface] = u_secondary + end + i_secondary += i_secondary_step # incrementing i_secondary suffices + end + end + + return nothing +end + function calc_interface_flux!(backend::Nothing, surface_flux_values, mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, have_nonconservative_terms, - equations, surface_integral, dg::DG, cache) + equations, surface_integral, + dg::DGSEM{<:LobattoLegendreBasis}, cache) @unpack neighbor_ids, node_indices = cache.interfaces + # Take for Gauss-Lobatto-Legendre (GLL) the interface normals from the outer volume nodes, i.e., + # element data. @unpack contravariant_vectors = cache.elements index_range = eachnode(dg) @@ -179,7 +312,8 @@ function calc_interface_flux!(backend::Backend, surface_flux_values, mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, have_nonconservative_terms, - equations, surface_integral, dg::DG, cache) + equations, surface_integral, + dg::DGSEM{<:LobattoLegendreBasis}, cache) ninterfaces(cache.interfaces) == 0 && return nothing @unpack neighbor_ids, node_indices = cache.interfaces @unpack contravariant_vectors = cache.elements @@ -217,7 +351,7 @@ end T8codeMesh{2}}}, have_nonconservative_terms, equations, surface_integral, - SolverT::Type{<:DG}, + SolverT::Type{<:DGSEM{<:LobattoLegendreBasis}}, u_interface, interface, neighbor_ids, node_indices, contravariant_vectors, @@ -265,8 +399,9 @@ end calc_interface_flux!(surface_flux_values, MeshT, have_nonconservative_terms, equations, surface_integral, SolverT, u_interface, - interface, normal_direction, node, primary_direction, - primary_element, node_secondary, + interface, normal_direction, node, + primary_direction, primary_element, + node_secondary, secondary_direction, secondary_element) # Increment primary element indices to pull the normal direction @@ -279,6 +414,81 @@ end return nothing end +function calc_interface_flux!(backend::Nothing, surface_flux_values, + mesh::Union{P4estMesh{2}, P4estMeshView{2}}, + have_nonconservative_terms, + equations, surface_integral, + dg::DGSEM{<:GaussLegendreBasis}, cache) + @unpack neighbor_ids, node_indices = cache.interfaces + # Take for Gauss-Legendre (GL) the interface normals from the interfaces, i.e., + # interface data. + @unpack normal_directions = cache.interfaces + index_range = eachnode(dg) + + @threaded for interface in eachinterface(dg, cache) + calc_interface_flux_per_interface!(surface_flux_values, typeof(mesh), + have_nonconservative_terms, + equations, surface_integral, typeof(dg), + cache.interfaces.u, interface, + neighbor_ids, node_indices, + normal_directions, index_range) + end + + return nothing +end + +@inline function calc_interface_flux_per_interface!(surface_flux_values, + MeshT::Type{<:Union{P4estMesh{2}, + P4estMeshView{2}}}, + have_nonconservative_terms, + equations, surface_integral, + SolverT::Type{<:DGSEM{<:GaussLegendreBasis}}, + u_interface, interface, + neighbor_ids, + node_indices, normal_directions, + index_range) + index_end = last(index_range) + + # Get element and side index information on the primary element + primary_element = neighbor_ids[1, interface] + primary_indices = node_indices[1, interface] + primary_direction = indices2direction(primary_indices) + + # Get element and side index information on the secondary element + secondary_element = neighbor_ids[2, interface] + secondary_indices = node_indices[2, interface] + secondary_direction = indices2direction(secondary_indices) + + # Initiate the secondary index to be used in the surface for loop. + # This index on the primary side will always run forward but + # the secondary index might need to run backwards for flipped sides. + if :i_backward in secondary_indices + node_secondary = index_end + node_secondary_step = -1 + else + node_secondary = 1 + node_secondary_step = 1 + end + + for node in index_range + # This seems to be faster than `@views normal_direction = normal_directions[:, node, interface]` + normal_direction = SVector(normal_directions[1, node, interface], + normal_directions[2, node, interface]) + + calc_interface_flux!(surface_flux_values, MeshT, have_nonconservative_terms, + equations, surface_integral, SolverT, u_interface, + interface, normal_direction, node, + primary_direction, primary_element, + node_secondary, + secondary_direction, secondary_element) + + # Increment the surface node index along the secondary element + node_secondary += node_secondary_step + end + + return nothing +end + # Inlined version of the interface flux computation for conservation laws @inline function calc_interface_flux!(surface_flux_values, ::Type{<:Union{P4estMesh{2}, @@ -906,7 +1116,7 @@ function calc_surface_integral!(backend::Nothing, du, u, mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, equations, surface_integral::SurfaceIntegralWeakForm, - dg::DGSEM, cache) + dg::DGSEM{<:LobattoLegendreBasis}, cache) @unpack inverse_weights = dg.basis @unpack surface_flux_values = cache.elements @@ -917,12 +1127,27 @@ function calc_surface_integral!(backend::Nothing, du, u, end end +function calc_surface_integral!(backend::Nothing, du, u, + mesh::Union{P4estMesh{2}, P4estMeshView{2}}, + equations, surface_integral::SurfaceIntegralWeakForm, + dg::DGSEM{<:GaussLegendreBasis}, cache) + @unpack boundary_interpolation_inverse_weights = dg.basis + @unpack surface_flux_values = cache.elements + + @threaded for element in eachelement(dg, cache) + calc_surface_integral_per_element!(du, typeof(mesh), equations, + surface_integral, dg, + boundary_interpolation_inverse_weights, + surface_flux_values, element) + end +end + function calc_surface_integral!(backend::Backend, du, u, mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}}, equations, surface_integral::SurfaceIntegralWeakForm, - dg::DGSEM, cache) + dg::DGSEM{<:LobattoLegendreBasis}, cache) nelements(dg, cache) == 0 && return nothing @unpack inverse_weights = dg.basis @unpack surface_flux_values = cache.elements @@ -939,7 +1164,8 @@ end T8codeMesh{2}}}, equations, surface_integral::SurfaceIntegralWeakForm, - dg::DGSEM, factor, + dg::DGSEM{<:LobattoLegendreBasis}, + factor, surface_flux_values) element = @index(Global) calc_surface_integral_per_element!(du, MeshT, equations, surface_integral, @@ -952,9 +1178,9 @@ end T8codeMesh{2}}}, equations, surface_integral::SurfaceIntegralWeakForm, - dg::DGSEM, factor, - surface_flux_values, - element) + dg::DGSEM{<:LobattoLegendreBasis}, + factor, + surface_flux_values, element) # Note that all fluxes have been computed with outward-pointing normal vectors. # This computes the **negative** surface integral contribution, # i.e., M^{-1} * boundary_interpolation^T (which is for Gauss-Lobatto DGSEM just M^{-1} * B) @@ -991,6 +1217,59 @@ end return nothing end +function calc_surface_integral_per_element!(du, + ::Type{<:Union{P4estMesh{2}, + P4estMeshView{2}}}, + equations, + surface_integral::SurfaceIntegralWeakForm, + dg::DGSEM{<:GaussLegendreBasis}, + boundary_interpolation_inverse_weights, + surface_flux_values, element) + # Note that all fluxes have been computed with outward-pointing normal vectors. + # This computes the **negative** surface integral contribution, + # i.e., M^{-1} * boundary_interpolation^T + # and the missing "-" is taken care of by `apply_jacobian!`. + # + # We also use explicit assignments instead of `+=` to let `@muladd` turn these + # into FMAs (see comment at the top of the file). + for l in eachnode(dg) + for v in eachvariable(equations) + # Aliases for repeatedly accessed variables + surface_flux_minus_x = surface_flux_values[v, l, 1, element] + surface_flux_plus_x = surface_flux_values[v, l, 2, element] + for ii in eachnode(dg) + # surface at -x + du[v, ii, l, element] = (du[v, ii, l, element] + + surface_flux_minus_x * + boundary_interpolation_inverse_weights[ii, + 1]) + # surface at +x + du[v, ii, l, element] = (du[v, ii, l, element] + + surface_flux_plus_x * + boundary_interpolation_inverse_weights[ii, + 2]) + end + + surface_flux_minus_y = surface_flux_values[v, l, 3, element] + surface_flux_plus_y = surface_flux_values[v, l, 4, element] + for jj in eachnode(dg) + # surface at -y + du[v, l, jj, element] = (du[v, l, jj, element] + + surface_flux_minus_y * + boundary_interpolation_inverse_weights[jj, + 1]) + # surface at +y + du[v, l, jj, element] = (du[v, l, jj, element] + + surface_flux_plus_y * + boundary_interpolation_inverse_weights[jj, + 2]) + end + end + end + + return nothing +end + # Call this for coupled P4estMeshView simulations. # The coupling calculations (especially boundary conditions) require data from the parent mesh, which is why # the additional variable u_parent is needed, compared to non-coupled systems. diff --git a/src/solvers/dgsem_structured/containers_2d.jl b/src/solvers/dgsem_structured/containers_2d.jl index 71151a785d2..03241705266 100644 --- a/src/solvers/dgsem_structured/containers_2d.jl +++ b/src/solvers/dgsem_structured/containers_2d.jl @@ -60,7 +60,7 @@ end # Calculate Jacobian matrix of the mapping from the reference element to the element in the physical domain function calc_jacobian_matrix!(jacobian_matrix, element, node_coordinates::AbstractArray{<:Any, 4}, - basis::LobattoLegendreBasis) + basis::AbstractBasisSBP) @unpack derivative_matrix = basis # The code below is equivalent to the following matrix multiplications, which diff --git a/src/solvers/dgsem_t8code/containers.jl b/src/solvers/dgsem_t8code/containers.jl index dea0bf1a9f1..56db011f3b1 100644 --- a/src/solvers/dgsem_t8code/containers.jl +++ b/src/solvers/dgsem_t8code/containers.jl @@ -54,7 +54,8 @@ function count_required_surfaces(mesh::T8codeMesh) end # Compatibility to `dgsem_p4est/containers.jl`. -function init_interfaces!(interfaces, mesh::T8codeMesh) +function init_interfaces!(interfaces, elements, + mesh::T8codeMesh, basis) # Already computed. Do nothing. return nothing end diff --git a/src/solvers/dgsem_t8code/dg.jl b/src/solvers/dgsem_t8code/dg.jl index ce35d524f8b..47169c5d8e7 100644 --- a/src/solvers/dgsem_t8code/dg.jl +++ b/src/solvers/dgsem_t8code/dg.jl @@ -33,5 +33,6 @@ include("containers_2d.jl") include("containers_3d.jl") include("containers_parallel.jl") + include("dg_parallel.jl") end # @muladd diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 739706d209b..f9be047fba6 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -1278,8 +1278,8 @@ function calc_surface_integral!(backend::Nothing, du, u, end function calc_surface_integral!(backend::Nothing, du, u, - mesh::Union{TreeMesh{2}, StructuredMesh{2}, - StructuredMeshView{2}}, + mesh::Union{TreeMesh{2}, + StructuredMesh{2}, StructuredMeshView{2}}, equations, surface_integral::SurfaceIntegralWeakForm, dg::DGSEM{<:GaussLegendreBasis}, cache) @unpack boundary_interpolation_inverse_weights = dg.basis diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index 4ebd2599314..5147ce29aa3 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -30,6 +30,17 @@ isdir(outdir) && rm(outdir, recursive = true) @test real(semi32.mesh) == Float64 end +@trixi_testset "elixir_advection_basic.jl (Gauss-Legendre)" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"), + solver=DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs, + basis_type = GaussLegendreBasis), + cfl=0.8, + l2=[3.721398353159235e-6], linf=[1.8621131703255855e-5]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + @test_allocations(Trixi.rhs!, semi, sol, 1000) +end + @trixi_testset "elixir_advection_basic.jl (Float32)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic_gpu.jl"), # Expected errors are exactly the same as with TreeMesh! @@ -62,6 +73,15 @@ end @test_allocations(Trixi.rhs!, semi, sol, 1000) end +@trixi_testset "elixir_advection_flag_gauss_legendre.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_advection_flag_gauss_legendre.jl"), + l2=[0.0004734270965231037], linf=[0.002206239481024719]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + @test_allocations(Trixi.rhs!, semi, sol, 1000) +end + @trixi_testset "elixir_advection_unstructured_flag.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_unstructured_flag.jl"), l2=[0.0005379687442422346],