diff --git a/src/solvers/dgsem_p4est/dg_2d.jl b/src/solvers/dgsem_p4est/dg_2d.jl index 3c86828918..1b38eaa2d4 100644 --- a/src/solvers/dgsem_p4est/dg_2d.jl +++ b/src/solvers/dgsem_p4est/dg_2d.jl @@ -13,11 +13,15 @@ function create_cache(mesh::Union{P4estMesh{2}, T8codeMesh{2}}, equations, MA2d = MArray{Tuple{nvariables(equations), nnodes(mortar_l2)}, uEltype, 2, nvariables(equations) * nnodes(mortar_l2)} - fstar_upper_threaded = MA2d[MA2d(undef) for _ in 1:Threads.nthreads()] - fstar_lower_threaded = MA2d[MA2d(undef) for _ in 1:Threads.nthreads()] + fstar_primary_upper_threaded = MA2d[MA2d(undef) for _ in 1:Threads.nthreads()] + fstar_primary_lower_threaded = MA2d[MA2d(undef) for _ in 1:Threads.nthreads()] + fstar_secondary_upper_threaded = MA2d[MA2d(undef) for _ in 1:Threads.nthreads()] + fstar_secondary_lower_threaded = MA2d[MA2d(undef) for _ in 1:Threads.nthreads()] u_threaded = MA2d[MA2d(undef) for _ in 1:Threads.nthreads()] - (; fstar_upper_threaded, fstar_lower_threaded, u_threaded) + (; fstar_primary_upper_threaded, fstar_primary_lower_threaded, + fstar_secondary_upper_threaded, fstar_secondary_lower_threaded, + u_threaded) end # index_to_start_step_2d(index::Symbol, index_range) @@ -451,13 +455,17 @@ function calc_mortar_flux!(surface_flux_values, surface_integral, dg::DG, cache) @unpack neighbor_ids, node_indices = cache.mortars @unpack contravariant_vectors = cache.elements - @unpack fstar_upper_threaded, fstar_lower_threaded = cache + @unpack (fstar_primary_upper_threaded, fstar_primary_lower_threaded, + fstar_secondary_upper_threaded, fstar_secondary_lower_threaded) = cache index_range = eachnode(dg) @threaded for mortar in eachmortar(dg, cache) # Choose thread-specific pre-allocated container - fstar = (fstar_lower_threaded[Threads.threadid()], - fstar_upper_threaded[Threads.threadid()]) + fstar_primary = (fstar_primary_lower_threaded[Threads.threadid()], + fstar_primary_upper_threaded[Threads.threadid()]) + + fstar_secondary = (fstar_secondary_lower_threaded[Threads.threadid()], + fstar_secondary_upper_threaded[Threads.threadid()]) # Get index information on the small elements small_indices = node_indices[1, mortar] @@ -480,7 +488,8 @@ function calc_mortar_flux!(surface_flux_values, contravariant_vectors, i_small, j_small, element) - calc_mortar_flux!(fstar, mesh, nonconservative_terms, equations, + calc_mortar_flux!(fstar_primary, fstar_secondary, + mesh, nonconservative_terms, equations, surface_integral, dg, cache, mortar, position, normal_direction, node) @@ -501,14 +510,15 @@ function calc_mortar_flux!(surface_flux_values, # "mortar_fluxes_to_elements!" instead. mortar_fluxes_to_elements!(surface_flux_values, mesh, equations, mortar_l2, dg, cache, - mortar, fstar, u_buffer) + mortar, fstar_primary, fstar_secondary, + u_buffer) end return nothing end # Inlined version of the mortar flux computation on small elements for conservation laws -@inline function calc_mortar_flux!(fstar, +@inline function calc_mortar_flux!(fstar_primary, fstar_secondary, mesh::Union{P4estMesh{2}, T8codeMesh{2}}, nonconservative_terms::False, equations, surface_integral, dg::DG, cache, @@ -523,12 +533,13 @@ end flux = surface_flux(u_ll, u_rr, normal_direction, equations) # Copy flux to buffer - set_node_vars!(fstar[position_index], flux, equations, dg, node_index) + set_node_vars!(fstar_primary[position_index], flux, equations, dg, node_index) + set_node_vars!(fstar_secondary[position_index], flux, equations, dg, node_index) end # Inlined version of the mortar flux computation on small elements for equations with conservative and # nonconservative terms -@inline function calc_mortar_flux!(fstar, +@inline function calc_mortar_flux!(fstar_primary, fstar_secondary, mesh::Union{P4estMesh{2}, T8codeMesh{2}}, nonconservative_terms::True, equations, surface_integral, dg::DG, cache, @@ -547,19 +558,25 @@ end # The nonconservative flux is scaled by a factor of 0.5 based on # the interpretation of global SBP operators coupled discontinuously via # central fluxes/SATs - noncons = nonconservative_flux(u_ll, u_rr, normal_direction, equations) + noncons_primary = nonconservative_flux(u_ll, u_rr, normal_direction, equations) + noncons_secondary = nonconservative_flux(u_rr, u_ll, normal_direction, equations) - flux_plus_noncons = flux + 0.5f0 * noncons + flux_plus_noncons_primary = flux + 0.5f0 * noncons_primary + flux_plus_noncons_secondary = flux + 0.5f0 * noncons_secondary # Copy to buffer - set_node_vars!(fstar[position_index], flux_plus_noncons, equations, dg, node_index) + set_node_vars!(fstar_primary[position_index], flux_plus_noncons_primary, equations, + dg, node_index) + set_node_vars!(fstar_secondary[position_index], flux_plus_noncons_secondary, + equations, dg, node_index) end @inline function mortar_fluxes_to_elements!(surface_flux_values, mesh::Union{P4estMesh{2}, T8codeMesh{2}}, equations, mortar_l2::LobattoLegendreMortarL2, - dg::DGSEM, cache, mortar, fstar, u_buffer) + dg::DGSEM, cache, mortar, fstar_primary, + fstar_secondary, u_buffer) @unpack neighbor_ids, node_indices = cache.mortars # Copy solution small to small @@ -570,16 +587,16 @@ end element = neighbor_ids[position, mortar] for i in eachnode(dg) for v in eachvariable(equations) - surface_flux_values[v, i, small_direction, element] = fstar[position][v, - i] + surface_flux_values[v, i, small_direction, element] = fstar_primary[position][v, + i] end end end # Project small fluxes to large element. multiply_dimensionwise!(u_buffer, - mortar_l2.reverse_upper, fstar[2], - mortar_l2.reverse_lower, fstar[1]) + mortar_l2.reverse_upper, fstar_secondary[2], + mortar_l2.reverse_lower, fstar_secondary[1]) # The flux is calculated in the outward direction of the small elements, # so the sign must be switched to get the flux in outward direction diff --git a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl index 23968a2bfa..baa9f7bf51 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parabolic.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parabolic.jl @@ -398,7 +398,8 @@ end mesh::Union{P4estMesh{2}, T8codeMesh{2}}, equations::AbstractEquationsParabolic, mortar_l2::LobattoLegendreMortarL2, - dg::DGSEM, cache, mortar, fstar, u_buffer) + dg::DGSEM, cache, mortar, fstar_primary, + fstar_secondary, u_buffer) @unpack neighbor_ids, node_indices = cache.mortars # Copy solution small to small small_indices = node_indices[1, mortar] @@ -408,16 +409,16 @@ end element = neighbor_ids[position, mortar] for i in eachnode(dg) for v in eachvariable(equations) - surface_flux_values[v, i, small_direction, element] = fstar[position][v, - i] + surface_flux_values[v, i, small_direction, element] = fstar_primary[position][v, + i] end end end # Project small fluxes to large element. multiply_dimensionwise!(u_buffer, - mortar_l2.reverse_upper, fstar[2], - mortar_l2.reverse_lower, fstar[1]) + mortar_l2.reverse_upper, fstar_secondary[2], + mortar_l2.reverse_lower, fstar_secondary[1]) # Copy interpolated flux values from buffer to large element face in the # correct orientation. @@ -772,13 +773,13 @@ function calc_mortar_flux_divergence!(surface_flux_values, surface_integral, dg::DG, cache) @unpack neighbor_ids, node_indices = cache.mortars @unpack contravariant_vectors = cache.elements - @unpack fstar_upper_threaded, fstar_lower_threaded = cache + @unpack fstar_primary_upper_threaded, fstar_primary_lower_threaded = cache index_range = eachnode(dg) @threaded for mortar in eachmortar(dg, cache) # Choose thread-specific pre-allocated container - fstar = (fstar_lower_threaded[Threads.threadid()], - fstar_upper_threaded[Threads.threadid()]) + fstar = (fstar_primary_lower_threaded[Threads.threadid()], + fstar_primary_upper_threaded[Threads.threadid()]) for position in 1:2 for node in eachnode(dg) @@ -802,7 +803,7 @@ function calc_mortar_flux_divergence!(surface_flux_values, # this reuses the hyperbolic version of `mortar_fluxes_to_elements!` mortar_fluxes_to_elements!(surface_flux_values, mesh, equations, mortar_l2, dg, cache, - mortar, fstar, u_buffer) + mortar, fstar, fstar, u_buffer) end return nothing @@ -813,7 +814,7 @@ end # The reasoning is that parabolic fluxes are treated like conservative # terms (e.g., we compute a viscous conservative "flux") and thus no # non-conservative terms are present. -@inline function calc_mortar_flux!(fstar, +@inline function calc_mortar_flux!(fstar_primary, fstar_secondary, mesh::Union{P4estMesh{2}, T8codeMesh{2}}, nonconservative_terms::False, equations::AbstractEquationsParabolic, @@ -830,7 +831,8 @@ end flux_ = 0.5f0 * (u_ll + u_rr) # Copy flux to buffer - set_node_vars!(fstar[position_index], flux_, equations, dg, node_index) + set_node_vars!(fstar_primary[position_index], flux_, equations, dg, node_index) + set_node_vars!(fstar_secondary[position_index], flux_, equations, dg, node_index) end # TODO: parabolic, finish implementing `calc_boundary_flux_gradients!` and `calc_boundary_flux_divergence!` diff --git a/src/solvers/dgsem_p4est/dg_2d_parallel.jl b/src/solvers/dgsem_p4est/dg_2d_parallel.jl index 3bf0cd0cab..633c6f1d06 100644 --- a/src/solvers/dgsem_p4est/dg_2d_parallel.jl +++ b/src/solvers/dgsem_p4est/dg_2d_parallel.jl @@ -209,13 +209,13 @@ function calc_mpi_mortar_flux!(surface_flux_values, surface_integral, dg::DG, cache) @unpack local_neighbor_ids, local_neighbor_positions, node_indices = cache.mpi_mortars @unpack contravariant_vectors = cache.elements - @unpack fstar_upper_threaded, fstar_lower_threaded = cache + @unpack fstar_primary_upper_threaded, fstar_primary_lower_threaded = cache index_range = eachnode(dg) @threaded for mortar in eachmpimortar(dg, cache) # Choose thread-specific pre-allocated container - fstar = (fstar_lower_threaded[Threads.threadid()], - fstar_upper_threaded[Threads.threadid()]) + fstar = (fstar_primary_lower_threaded[Threads.threadid()], + fstar_primary_upper_threaded[Threads.threadid()]) # Get index information on the small elements small_indices = node_indices[1, mortar] diff --git a/src/solvers/dgsem_tree/dg_2d.jl b/src/solvers/dgsem_tree/dg_2d.jl index 0585093b7a..069d4f09bc 100644 --- a/src/solvers/dgsem_tree/dg_2d.jl +++ b/src/solvers/dgsem_tree/dg_2d.jl @@ -95,14 +95,17 @@ function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMe # TODO: Taal performance using different types MA2d = MArray{Tuple{nvariables(equations), nnodes(mortar_l2)}, uEltype, 2, nvariables(equations) * nnodes(mortar_l2)} - fstar_upper_threaded = MA2d[MA2d(undef) for _ in 1:Threads.nthreads()] - fstar_lower_threaded = MA2d[MA2d(undef) for _ in 1:Threads.nthreads()] + fstar_primary_upper_threaded = MA2d[MA2d(undef) for _ in 1:Threads.nthreads()] + fstar_primary_lower_threaded = MA2d[MA2d(undef) for _ in 1:Threads.nthreads()] + fstar_secondary_upper_threaded = MA2d[MA2d(undef) for _ in 1:Threads.nthreads()] + fstar_secondary_lower_threaded = MA2d[MA2d(undef) for _ in 1:Threads.nthreads()] # A2d = Array{uEltype, 2} # fstar_upper_threaded = [A2d(undef, nvariables(equations), nnodes(mortar_l2)) for _ in 1:Threads.nthreads()] # fstar_lower_threaded = [A2d(undef, nvariables(equations), nnodes(mortar_l2)) for _ in 1:Threads.nthreads()] - (; fstar_upper_threaded, fstar_lower_threaded) + (; fstar_primary_upper_threaded, fstar_primary_lower_threaded, + fstar_secondary_upper_threaded, fstar_secondary_lower_threaded) end # TODO: Taal discuss/refactor timer, allowing users to pass a custom timer? @@ -900,23 +903,31 @@ function calc_mortar_flux!(surface_flux_values, surface_integral, dg::DG, cache) @unpack surface_flux = surface_integral @unpack u_lower, u_upper, orientations = cache.mortars - @unpack fstar_upper_threaded, fstar_lower_threaded = cache + @unpack (fstar_primary_upper_threaded, fstar_primary_lower_threaded, + fstar_secondary_upper_threaded, fstar_secondary_lower_threaded) = cache @threaded for mortar in eachmortar(dg, cache) # Choose thread-specific pre-allocated container - fstar_upper = fstar_upper_threaded[Threads.threadid()] - fstar_lower = fstar_lower_threaded[Threads.threadid()] + fstar_primary_upper = fstar_primary_upper_threaded[Threads.threadid()] + fstar_primary_lower = fstar_primary_lower_threaded[Threads.threadid()] + fstar_secondary_upper = fstar_secondary_upper_threaded[Threads.threadid()] + fstar_secondary_lower = fstar_secondary_lower_threaded[Threads.threadid()] # Calculate fluxes orientation = orientations[mortar] - calc_fstar!(fstar_upper, equations, surface_flux, dg, u_upper, mortar, + calc_fstar!(fstar_primary_upper, equations, surface_flux, dg, u_upper, mortar, orientation) - calc_fstar!(fstar_lower, equations, surface_flux, dg, u_lower, mortar, + calc_fstar!(fstar_primary_lower, equations, surface_flux, dg, u_lower, mortar, + orientation) + calc_fstar!(fstar_secondary_upper, equations, surface_flux, dg, u_upper, mortar, + orientation) + calc_fstar!(fstar_secondary_lower, equations, surface_flux, dg, u_lower, mortar, orientation) mortar_fluxes_to_elements!(surface_flux_values, mesh, equations, mortar_l2, dg, cache, - mortar, fstar_upper, fstar_lower) + mortar, fstar_primary_upper, fstar_primary_lower, + fstar_secondary_upper, fstar_secondary_lower) end return nothing @@ -929,18 +940,25 @@ function calc_mortar_flux!(surface_flux_values, surface_integral, dg::DG, cache) surface_flux, nonconservative_flux = surface_integral.surface_flux @unpack u_lower, u_upper, orientations, large_sides = cache.mortars - @unpack fstar_upper_threaded, fstar_lower_threaded = cache + @unpack (fstar_primary_upper_threaded, fstar_primary_lower_threaded, + fstar_secondary_upper_threaded, fstar_secondary_lower_threaded) = cache @threaded for mortar in eachmortar(dg, cache) # Choose thread-specific pre-allocated container - fstar_upper = fstar_upper_threaded[Threads.threadid()] - fstar_lower = fstar_lower_threaded[Threads.threadid()] + fstar_primary_upper = fstar_primary_upper_threaded[Threads.threadid()] + fstar_primary_lower = fstar_primary_lower_threaded[Threads.threadid()] + fstar_secondary_upper = fstar_secondary_upper_threaded[Threads.threadid()] + fstar_secondary_lower = fstar_secondary_lower_threaded[Threads.threadid()] # Calculate fluxes orientation = orientations[mortar] - calc_fstar!(fstar_upper, equations, surface_flux, dg, u_upper, mortar, + calc_fstar!(fstar_primary_upper, equations, surface_flux, dg, u_upper, mortar, + orientation) + calc_fstar!(fstar_primary_lower, equations, surface_flux, dg, u_lower, mortar, + orientation) + calc_fstar!(fstar_secondary_upper, equations, surface_flux, dg, u_upper, mortar, orientation) - calc_fstar!(fstar_lower, equations, surface_flux, dg, u_lower, mortar, + calc_fstar!(fstar_secondary_lower, equations, surface_flux, dg, u_lower, mortar, orientation) # Add nonconservative fluxes. @@ -958,14 +976,26 @@ function calc_mortar_flux!(surface_flux_values, u_lower_ll, u_lower_rr = get_surface_node_vars(u_lower, equations, dg, i, mortar) # Call pointwise nonconservative term - noncons_upper = nonconservative_flux(u_upper_ll, u_upper_rr, - orientation, equations) - noncons_lower = nonconservative_flux(u_lower_ll, u_lower_rr, - orientation, equations) + noncons_primary_upper = nonconservative_flux(u_upper_ll, u_upper_rr, + orientation, equations) + noncons_primary_lower = nonconservative_flux(u_lower_ll, u_lower_rr, + orientation, equations) + noncons_secondary_upper = nonconservative_flux(u_upper_rr, u_upper_ll, + orientation, equations) + noncons_secondary_lower = nonconservative_flux(u_lower_rr, u_lower_ll, + orientation, equations) # Add to primary and secondary temporary storage - multiply_add_to_node_vars!(fstar_upper, 0.5f0, noncons_upper, equations, + multiply_add_to_node_vars!(fstar_primary_upper, 0.5f0, + noncons_primary_upper, equations, dg, i) - multiply_add_to_node_vars!(fstar_lower, 0.5f0, noncons_lower, equations, + multiply_add_to_node_vars!(fstar_primary_lower, 0.5f0, + noncons_primary_lower, equations, + dg, i) + multiply_add_to_node_vars!(fstar_secondary_upper, 0.5f0, + noncons_secondary_upper, equations, + dg, i) + multiply_add_to_node_vars!(fstar_secondary_lower, 0.5f0, + noncons_secondary_lower, equations, dg, i) end else # large_sides[mortar] == 2 -> small elements on the left @@ -976,21 +1006,34 @@ function calc_mortar_flux!(surface_flux_values, u_lower_ll, u_lower_rr = get_surface_node_vars(u_lower, equations, dg, i, mortar) # Call pointwise nonconservative term - noncons_upper = nonconservative_flux(u_upper_rr, u_upper_ll, - orientation, equations) - noncons_lower = nonconservative_flux(u_lower_rr, u_lower_ll, - orientation, equations) + noncons_primary_upper = nonconservative_flux(u_upper_rr, u_upper_ll, + orientation, equations) + noncons_primary_lower = nonconservative_flux(u_lower_rr, u_lower_ll, + orientation, equations) + noncons_secondary_upper = nonconservative_flux(u_upper_ll, u_upper_rr, + orientation, equations) + noncons_secondary_lower = nonconservative_flux(u_lower_ll, u_lower_rr, + orientation, equations) # Add to primary and secondary temporary storage - multiply_add_to_node_vars!(fstar_upper, 0.5f0, noncons_upper, equations, + multiply_add_to_node_vars!(fstar_primary_upper, 0.5f0, + noncons_primary_upper, equations, + dg, i) + multiply_add_to_node_vars!(fstar_primary_lower, 0.5f0, + noncons_primary_lower, equations, + dg, i) + multiply_add_to_node_vars!(fstar_secondary_upper, 0.5f0, + noncons_secondary_upper, equations, dg, i) - multiply_add_to_node_vars!(fstar_lower, 0.5f0, noncons_lower, equations, + multiply_add_to_node_vars!(fstar_secondary_lower, 0.5f0, + noncons_secondary_lower, equations, dg, i) end end mortar_fluxes_to_elements!(surface_flux_values, mesh, equations, mortar_l2, dg, cache, - mortar, fstar_upper, fstar_lower) + mortar, fstar_primary_upper, fstar_primary_lower, + fstar_secondary_upper, fstar_secondary_lower) end return nothing @@ -1015,7 +1058,10 @@ end mesh::TreeMesh{2}, equations, mortar_l2::LobattoLegendreMortarL2, dg::DGSEM, cache, - mortar, fstar_upper, fstar_lower) + mortar, fstar_primary_upper, + fstar_primary_lower, + fstar_secondary_upper, + fstar_secondary_lower) large_element = cache.mortars.neighbor_ids[3, mortar] upper_element = cache.mortars.neighbor_ids[2, mortar] lower_element = cache.mortars.neighbor_ids[1, mortar] @@ -1038,8 +1084,8 @@ end direction = 4 end end - surface_flux_values[:, :, direction, upper_element] .= fstar_upper - surface_flux_values[:, :, direction, lower_element] .= fstar_lower + surface_flux_values[:, :, direction, upper_element] .= fstar_primary_upper + surface_flux_values[:, :, direction, lower_element] .= fstar_primary_lower # Project small fluxes to large element if cache.mortars.large_sides[mortar] == 1 # -> large element on left side @@ -1078,8 +1124,8 @@ end # depends on the types of fstar_upper/fstar_lower and dg.l2mortar_reverse_upper. # Using StaticArrays for both makes the code above faster for common test cases. multiply_dimensionwise!(view(surface_flux_values, :, :, direction, large_element), - mortar_l2.reverse_upper, fstar_upper, - mortar_l2.reverse_lower, fstar_lower) + mortar_l2.reverse_upper, fstar_secondary_upper, + mortar_l2.reverse_lower, fstar_secondary_lower) return nothing end diff --git a/src/solvers/dgsem_tree/dg_2d_parabolic.jl b/src/solvers/dgsem_tree/dg_2d_parabolic.jl index 9123a3bcd2..22c29eb472 100644 --- a/src/solvers/dgsem_tree/dg_2d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_2d_parabolic.jl @@ -623,12 +623,12 @@ function calc_mortar_flux!(surface_flux_values, surface_integral, dg::DG, cache) @unpack surface_flux = surface_integral @unpack u_lower, u_upper, orientations = cache.mortars - @unpack fstar_upper_threaded, fstar_lower_threaded = cache + @unpack fstar_primary_upper_threaded, fstar_primary_lower_threaded = cache @threaded for mortar in eachmortar(dg, cache) # Choose thread-specific pre-allocated container - fstar_upper = fstar_upper_threaded[Threads.threadid()] - fstar_lower = fstar_lower_threaded[Threads.threadid()] + fstar_upper = fstar_primary_upper_threaded[Threads.threadid()] + fstar_lower = fstar_primary_lower_threaded[Threads.threadid()] # Calculate fluxes orientation = orientations[mortar] diff --git a/src/solvers/dgsem_tree/dg_2d_parallel.jl b/src/solvers/dgsem_tree/dg_2d_parallel.jl index d0287e00b8..b21933f0a5 100644 --- a/src/solvers/dgsem_tree/dg_2d_parallel.jl +++ b/src/solvers/dgsem_tree/dg_2d_parallel.jl @@ -769,23 +769,31 @@ function calc_mpi_mortar_flux!(surface_flux_values, surface_integral, dg::DG, cache) @unpack surface_flux = surface_integral @unpack u_lower, u_upper, orientations = cache.mpi_mortars - @unpack fstar_upper_threaded, fstar_lower_threaded = cache + @unpack fstar_primary_upper_threaded, fstar_primary_lower_threaded, fstar_secondary_upper_threaded, fstar_secondary_lower_threaded = cache @threaded for mortar in eachmpimortar(dg, cache) # Choose thread-specific pre-allocated container - fstar_upper = fstar_upper_threaded[Threads.threadid()] - fstar_lower = fstar_lower_threaded[Threads.threadid()] + fstar_primary_upper = fstar_primary_upper_threaded[Threads.threadid()] + fstar_primary_lower = fstar_primary_lower_threaded[Threads.threadid()] + fstar_secondary_upper = fstar_secondary_upper_threaded[Threads.threadid()] + fstar_secondary_lower = fstar_secondary_lower_threaded[Threads.threadid()] - # Calculate fluxes + # Because `nonconservative_terms` is `False` the primary and secondary fluxes + # are identical. So, we could possibly save on computation and just pass two copies later. orientation = orientations[mortar] - calc_fstar!(fstar_upper, equations, surface_flux, dg, u_upper, mortar, + calc_fstar!(fstar_primary_upper, equations, surface_flux, dg, u_upper, mortar, orientation) - calc_fstar!(fstar_lower, equations, surface_flux, dg, u_lower, mortar, + calc_fstar!(fstar_primary_lower, equations, surface_flux, dg, u_lower, mortar, + orientation) + calc_fstar!(fstar_secondary_upper, equations, surface_flux, dg, u_upper, mortar, + orientation) + calc_fstar!(fstar_secondary_lower, equations, surface_flux, dg, u_lower, mortar, orientation) mpi_mortar_fluxes_to_elements!(surface_flux_values, mesh, equations, mortar_l2, dg, cache, - mortar, fstar_upper, fstar_lower) + mortar, fstar_primary_upper, fstar_primary_lower, + fstar_secondary_upper, fstar_secondary_lower) end return nothing @@ -795,7 +803,10 @@ end mesh::ParallelTreeMesh{2}, equations, mortar_l2::LobattoLegendreMortarL2, dg::DGSEM, cache, - mortar, fstar_upper, fstar_lower) + mortar, fstar_primary_upper, + fstar_primary_lower, + fstar_secondary_upper, + fstar_secondary_lower) local_neighbor_ids = cache.mpi_mortars.local_neighbor_ids[mortar] local_neighbor_positions = cache.mpi_mortars.local_neighbor_positions[mortar] @@ -821,9 +832,9 @@ end end if position == 1 - surface_flux_values[:, :, direction, element] .= fstar_lower + surface_flux_values[:, :, direction, element] .= fstar_primary_lower elseif position == 2 - surface_flux_values[:, :, direction, element] .= fstar_upper + surface_flux_values[:, :, direction, element] .= fstar_primary_upper end else # position == 3 -> current element is large # Project small fluxes to large element @@ -846,8 +857,8 @@ end end multiply_dimensionwise!(view(surface_flux_values, :, :, direction, element), - mortar_l2.reverse_upper, fstar_upper, - mortar_l2.reverse_lower, fstar_lower) + mortar_l2.reverse_upper, fstar_secondary_upper, + mortar_l2.reverse_lower, fstar_secondary_lower) end end diff --git a/src/solvers/dgsem_tree/dg_3d.jl b/src/solvers/dgsem_tree/dg_3d.jl index 4b3bc81ff8..381534d207 100644 --- a/src/solvers/dgsem_tree/dg_3d.jl +++ b/src/solvers/dgsem_tree/dg_3d.jl @@ -115,24 +115,44 @@ function create_cache(mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3}, equations, mortar_l2::LobattoLegendreMortarL2, uEltype) # TODO: Taal compare performance of different types A3d = Array{uEltype, 3} - fstar_upper_left_threaded = A3d[A3d(undef, nvariables(equations), nnodes(mortar_l2), - nnodes(mortar_l2)) - for _ in 1:Threads.nthreads()] - fstar_upper_right_threaded = A3d[A3d(undef, nvariables(equations), - nnodes(mortar_l2), nnodes(mortar_l2)) - for _ in 1:Threads.nthreads()] - fstar_lower_left_threaded = A3d[A3d(undef, nvariables(equations), nnodes(mortar_l2), - nnodes(mortar_l2)) - for _ in 1:Threads.nthreads()] - fstar_lower_right_threaded = A3d[A3d(undef, nvariables(equations), - nnodes(mortar_l2), nnodes(mortar_l2)) - for _ in 1:Threads.nthreads()] + fstar_primary_upper_left_threaded = A3d[A3d(undef, nvariables(equations), + nnodes(mortar_l2), + nnodes(mortar_l2)) + for _ in 1:Threads.nthreads()] + fstar_primary_upper_right_threaded = A3d[A3d(undef, nvariables(equations), + nnodes(mortar_l2), nnodes(mortar_l2)) + for _ in 1:Threads.nthreads()] + fstar_primary_lower_left_threaded = A3d[A3d(undef, nvariables(equations), + nnodes(mortar_l2), + nnodes(mortar_l2)) + for _ in 1:Threads.nthreads()] + fstar_primary_lower_right_threaded = A3d[A3d(undef, nvariables(equations), + nnodes(mortar_l2), nnodes(mortar_l2)) + for _ in 1:Threads.nthreads()] + fstar_secondary_upper_left_threaded = A3d[A3d(undef, nvariables(equations), + nnodes(mortar_l2), + nnodes(mortar_l2)) + for _ in 1:Threads.nthreads()] + fstar_secondary_upper_right_threaded = A3d[A3d(undef, nvariables(equations), + nnodes(mortar_l2), + nnodes(mortar_l2)) + for _ in 1:Threads.nthreads()] + fstar_secondary_lower_left_threaded = A3d[A3d(undef, nvariables(equations), + nnodes(mortar_l2), + nnodes(mortar_l2)) + for _ in 1:Threads.nthreads()] + fstar_secondary_lower_right_threaded = A3d[A3d(undef, nvariables(equations), + nnodes(mortar_l2), + nnodes(mortar_l2)) + for _ in 1:Threads.nthreads()] fstar_tmp1_threaded = A3d[A3d(undef, nvariables(equations), nnodes(mortar_l2), nnodes(mortar_l2)) for _ in 1:Threads.nthreads()] - (; fstar_upper_left_threaded, fstar_upper_right_threaded, - fstar_lower_left_threaded, fstar_lower_right_threaded, + (; fstar_primary_upper_left_threaded, fstar_primary_upper_right_threaded, + fstar_primary_lower_left_threaded, fstar_primary_lower_right_threaded, + fstar_secondary_upper_left_threaded, fstar_secondary_upper_right_threaded, + fstar_secondary_lower_left_threaded, fstar_secondary_lower_right_threaded, fstar_tmp1_threaded) end @@ -1038,33 +1058,51 @@ function calc_mortar_flux!(surface_flux_values, surface_integral, dg::DG, cache) @unpack surface_flux = surface_integral @unpack u_lower_left, u_lower_right, u_upper_left, u_upper_right, orientations = cache.mortars - @unpack (fstar_upper_left_threaded, fstar_upper_right_threaded, - fstar_lower_left_threaded, fstar_lower_right_threaded, + @unpack (fstar_primary_upper_left_threaded, fstar_primary_upper_right_threaded, + fstar_primary_lower_left_threaded, fstar_primary_lower_right_threaded, + fstar_secondary_upper_left_threaded, fstar_secondary_upper_right_threaded, + fstar_secondary_lower_left_threaded, fstar_secondary_lower_right_threaded, fstar_tmp1_threaded) = cache @threaded for mortar in eachmortar(dg, cache) # Choose thread-specific pre-allocated container - fstar_upper_left = fstar_upper_left_threaded[Threads.threadid()] - fstar_upper_right = fstar_upper_right_threaded[Threads.threadid()] - fstar_lower_left = fstar_lower_left_threaded[Threads.threadid()] - fstar_lower_right = fstar_lower_right_threaded[Threads.threadid()] + fstar_primary_upper_left = fstar_primary_upper_left_threaded[Threads.threadid()] + fstar_primary_upper_right = fstar_primary_upper_right_threaded[Threads.threadid()] + fstar_primary_lower_left = fstar_primary_lower_left_threaded[Threads.threadid()] + fstar_primary_lower_right = fstar_primary_lower_right_threaded[Threads.threadid()] + fstar_secondary_upper_left = fstar_secondary_upper_left_threaded[Threads.threadid()] + fstar_secondary_upper_right = fstar_secondary_upper_right_threaded[Threads.threadid()] + fstar_secondary_lower_left = fstar_secondary_lower_left_threaded[Threads.threadid()] + fstar_secondary_lower_right = fstar_secondary_lower_right_threaded[Threads.threadid()] fstar_tmp1 = fstar_tmp1_threaded[Threads.threadid()] # Calculate fluxes orientation = orientations[mortar] - calc_fstar!(fstar_upper_left, equations, surface_flux, dg, u_upper_left, mortar, - orientation) - calc_fstar!(fstar_upper_right, equations, surface_flux, dg, u_upper_right, - mortar, orientation) - calc_fstar!(fstar_lower_left, equations, surface_flux, dg, u_lower_left, mortar, - orientation) - calc_fstar!(fstar_lower_right, equations, surface_flux, dg, u_lower_right, - mortar, orientation) + calc_fstar!(fstar_primary_upper_left, equations, surface_flux, dg, + u_upper_left, mortar, orientation) + calc_fstar!(fstar_primary_upper_right, equations, surface_flux, dg, + u_upper_right, mortar, orientation) + calc_fstar!(fstar_primary_lower_left, equations, surface_flux, dg, + u_lower_left, mortar, orientation) + calc_fstar!(fstar_primary_lower_right, equations, surface_flux, dg, + u_lower_right, mortar, orientation) + calc_fstar!(fstar_secondary_upper_left, equations, surface_flux, dg, + u_upper_left, mortar, orientation) + calc_fstar!(fstar_secondary_upper_right, equations, surface_flux, dg, + u_upper_right, mortar, orientation) + calc_fstar!(fstar_secondary_lower_left, equations, surface_flux, dg, + u_lower_left, mortar, orientation) + calc_fstar!(fstar_secondary_lower_right, equations, surface_flux, dg, + u_lower_right, mortar, orientation) mortar_fluxes_to_elements!(surface_flux_values, mesh, equations, mortar_l2, dg, cache, mortar, - fstar_upper_left, fstar_upper_right, - fstar_lower_left, fstar_lower_right, + fstar_primary_upper_left, fstar_primary_upper_right, + fstar_primary_lower_left, fstar_primary_lower_right, + fstar_secondary_upper_left, + fstar_secondary_upper_right, + fstar_secondary_lower_left, + fstar_secondary_lower_right, fstar_tmp1) end @@ -1078,28 +1116,42 @@ function calc_mortar_flux!(surface_flux_values, surface_integral, dg::DG, cache) surface_flux, nonconservative_flux = surface_integral.surface_flux @unpack u_lower_left, u_lower_right, u_upper_left, u_upper_right, orientations, large_sides = cache.mortars - @unpack (fstar_upper_left_threaded, fstar_upper_right_threaded, - fstar_lower_left_threaded, fstar_lower_right_threaded, + @unpack (fstar_primary_upper_left_threaded, fstar_primary_upper_right_threaded, + fstar_primary_lower_left_threaded, fstar_primary_lower_right_threaded, + fstar_secondary_upper_left_threaded, fstar_secondary_upper_right_threaded, + fstar_secondary_lower_left_threaded, fstar_secondary_lower_right_threaded, fstar_tmp1_threaded) = cache @threaded for mortar in eachmortar(dg, cache) # Choose thread-specific pre-allocated container - fstar_upper_left = fstar_upper_left_threaded[Threads.threadid()] - fstar_upper_right = fstar_upper_right_threaded[Threads.threadid()] - fstar_lower_left = fstar_lower_left_threaded[Threads.threadid()] - fstar_lower_right = fstar_lower_right_threaded[Threads.threadid()] + fstar_primary_upper_left = fstar_primary_upper_left_threaded[Threads.threadid()] + fstar_primary_upper_right = fstar_primary_upper_right_threaded[Threads.threadid()] + fstar_primary_lower_left = fstar_primary_lower_left_threaded[Threads.threadid()] + fstar_primary_lower_right = fstar_primary_lower_right_threaded[Threads.threadid()] + fstar_secondary_upper_left = fstar_secondary_upper_left_threaded[Threads.threadid()] + fstar_secondary_upper_right = fstar_secondary_upper_right_threaded[Threads.threadid()] + fstar_secondary_lower_left = fstar_secondary_lower_left_threaded[Threads.threadid()] + fstar_secondary_lower_right = fstar_secondary_lower_right_threaded[Threads.threadid()] fstar_tmp1 = fstar_tmp1_threaded[Threads.threadid()] # Calculate fluxes orientation = orientations[mortar] - calc_fstar!(fstar_upper_left, equations, surface_flux, dg, u_upper_left, mortar, - orientation) - calc_fstar!(fstar_upper_right, equations, surface_flux, dg, u_upper_right, - mortar, orientation) - calc_fstar!(fstar_lower_left, equations, surface_flux, dg, u_lower_left, mortar, - orientation) - calc_fstar!(fstar_lower_right, equations, surface_flux, dg, u_lower_right, - mortar, orientation) + calc_fstar!(fstar_primary_upper_left, equations, surface_flux, dg, + u_upper_left, mortar, orientation) + calc_fstar!(fstar_primary_upper_right, equations, surface_flux, dg, + u_upper_right, mortar, orientation) + calc_fstar!(fstar_primary_lower_left, equations, surface_flux, dg, + u_lower_left, mortar, orientation) + calc_fstar!(fstar_primary_lower_right, equations, surface_flux, dg, + u_lower_right, mortar, orientation) + calc_fstar!(fstar_secondary_upper_left, equations, surface_flux, dg, + u_upper_left, mortar, orientation) + calc_fstar!(fstar_secondary_upper_right, equations, surface_flux, dg, + u_upper_right, mortar, orientation) + calc_fstar!(fstar_secondary_lower_left, equations, surface_flux, dg, + u_lower_left, mortar, orientation) + calc_fstar!(fstar_secondary_lower_right, equations, surface_flux, dg, + u_lower_right, mortar, orientation) # Add nonconservative fluxes. # These need to be adapted on the geometry (left/right) since the order of @@ -1126,28 +1178,62 @@ function calc_mortar_flux!(surface_flux_values, dg, i, j, mortar) # Call pointwise nonconservative term - noncons_upper_left = nonconservative_flux(u_upper_left_ll, - u_upper_left_rr, orientation, - equations) - noncons_upper_right = nonconservative_flux(u_upper_right_ll, - u_upper_right_rr, - orientation, equations) - noncons_lower_left = nonconservative_flux(u_lower_left_ll, - u_lower_left_rr, orientation, - equations) - noncons_lower_right = nonconservative_flux(u_lower_right_ll, - u_lower_right_rr, - orientation, equations) + noncons_primary_upper_left = nonconservative_flux(u_upper_left_ll, + u_upper_left_rr, + orientation, + equations) + noncons_primary_upper_right = nonconservative_flux(u_upper_right_ll, + u_upper_right_rr, + orientation, + equations) + noncons_primary_lower_left = nonconservative_flux(u_lower_left_ll, + u_lower_left_rr, + orientation, + equations) + noncons_primary_lower_right = nonconservative_flux(u_lower_right_ll, + u_lower_right_rr, + orientation, + equations) + noncons_secondary_upper_left = nonconservative_flux(u_upper_left_rr, + u_upper_left_ll, + orientation, + equations) + noncons_secondary_upper_right = nonconservative_flux(u_upper_right_rr, + u_upper_right_ll, + orientation, + equations) + noncons_secondary_lower_left = nonconservative_flux(u_lower_left_rr, + u_lower_left_ll, + orientation, + equations) + noncons_secondary_lower_right = nonconservative_flux(u_lower_right_rr, + u_lower_right_ll, + orientation, + equations) # Add to primary and secondary temporary storage - multiply_add_to_node_vars!(fstar_upper_left, 0.5f0, noncons_upper_left, + multiply_add_to_node_vars!(fstar_primary_upper_left, 0.5f0, + noncons_primary_upper_left, equations, dg, i, j) - multiply_add_to_node_vars!(fstar_upper_right, 0.5f0, - noncons_upper_right, + multiply_add_to_node_vars!(fstar_primary_upper_right, 0.5f0, + noncons_primary_upper_right, equations, dg, i, j) - multiply_add_to_node_vars!(fstar_lower_left, 0.5f0, noncons_lower_left, + multiply_add_to_node_vars!(fstar_primary_lower_left, 0.5f0, + noncons_primary_lower_left, equations, dg, i, j) - multiply_add_to_node_vars!(fstar_lower_right, 0.5f0, - noncons_lower_right, + multiply_add_to_node_vars!(fstar_primary_lower_right, 0.5f0, + noncons_primary_lower_right, + equations, dg, i, j) + multiply_add_to_node_vars!(fstar_secondary_upper_left, 0.5f0, + noncons_secondary_upper_left, + equations, dg, i, j) + multiply_add_to_node_vars!(fstar_secondary_upper_right, 0.5f0, + noncons_secondary_upper_right, + equations, dg, i, j) + multiply_add_to_node_vars!(fstar_secondary_lower_left, 0.5f0, + noncons_secondary_lower_left, + equations, dg, i, j) + multiply_add_to_node_vars!(fstar_secondary_lower_right, 0.5f0, + noncons_secondary_lower_right, equations, dg, i, j) end else # large_sides[mortar] == 2 -> small elements on the left @@ -1168,36 +1254,74 @@ function calc_mortar_flux!(surface_flux_values, dg, i, j, mortar) # Call pointwise nonconservative term - noncons_upper_left = nonconservative_flux(u_upper_left_rr, - u_upper_left_ll, orientation, - equations) - noncons_upper_right = nonconservative_flux(u_upper_right_rr, - u_upper_right_ll, - orientation, equations) - noncons_lower_left = nonconservative_flux(u_lower_left_rr, - u_lower_left_ll, orientation, - equations) - noncons_lower_right = nonconservative_flux(u_lower_right_rr, - u_lower_right_ll, - orientation, equations) + noncons_primary_upper_left = nonconservative_flux(u_upper_left_rr, + u_upper_left_ll, + orientation, + equations) + noncons_primary_upper_right = nonconservative_flux(u_upper_right_rr, + u_upper_right_ll, + orientation, + equations) + noncons_primary_lower_left = nonconservative_flux(u_lower_left_rr, + u_lower_left_ll, + orientation, + equations) + noncons_primary_lower_right = nonconservative_flux(u_lower_right_rr, + u_lower_right_ll, + orientation, + equations) + noncons_secondary_upper_left = nonconservative_flux(u_upper_left_ll, + u_upper_left_rr, + orientation, + equations) + noncons_secondary_upper_right = nonconservative_flux(u_upper_right_ll, + u_upper_right_rr, + orientation, + equations) + noncons_secondary_lower_left = nonconservative_flux(u_lower_left_ll, + u_lower_left_rr, + orientation, + equations) + noncons_secondary_lower_right = nonconservative_flux(u_lower_right_ll, + u_lower_right_rr, + orientation, + equations) # Add to primary and secondary temporary storage - multiply_add_to_node_vars!(fstar_upper_left, 0.5f0, noncons_upper_left, + multiply_add_to_node_vars!(fstar_primary_upper_left, 0.5f0, + noncons_primary_upper_left, + equations, dg, i, j) + multiply_add_to_node_vars!(fstar_primary_upper_right, 0.5f0, + noncons_primary_upper_right, + equations, dg, i, j) + multiply_add_to_node_vars!(fstar_primary_lower_left, 0.5f0, + noncons_primary_lower_left, + equations, dg, i, j) + multiply_add_to_node_vars!(fstar_primary_lower_right, 0.5f0, + noncons_primary_lower_right, + equations, dg, i, j) + multiply_add_to_node_vars!(fstar_secondary_upper_left, 0.5f0, + noncons_secondary_upper_left, equations, dg, i, j) - multiply_add_to_node_vars!(fstar_upper_right, 0.5f0, - noncons_upper_right, + multiply_add_to_node_vars!(fstar_secondary_upper_right, 0.5f0, + noncons_secondary_upper_right, equations, dg, i, j) - multiply_add_to_node_vars!(fstar_lower_left, 0.5f0, noncons_lower_left, + multiply_add_to_node_vars!(fstar_secondary_lower_left, 0.5f0, + noncons_secondary_lower_left, equations, dg, i, j) - multiply_add_to_node_vars!(fstar_lower_right, 0.5f0, - noncons_lower_right, + multiply_add_to_node_vars!(fstar_secondary_lower_right, 0.5f0, + noncons_secondary_lower_right, equations, dg, i, j) end end mortar_fluxes_to_elements!(surface_flux_values, mesh, equations, mortar_l2, dg, cache, mortar, - fstar_upper_left, fstar_upper_right, - fstar_lower_left, fstar_lower_right, + fstar_primary_upper_left, fstar_primary_upper_right, + fstar_primary_lower_left, fstar_primary_lower_right, + fstar_secondary_upper_left, + fstar_secondary_upper_right, + fstar_secondary_lower_left, + fstar_secondary_lower_right, fstar_tmp1) end @@ -1224,8 +1348,14 @@ end mortar_l2::LobattoLegendreMortarL2, dg::DGSEM, cache, mortar, - fstar_upper_left, fstar_upper_right, - fstar_lower_left, fstar_lower_right, + fstar_primary_upper_left, + fstar_primary_upper_right, + fstar_primary_lower_left, + fstar_primary_lower_right, + fstar_secondary_upper_left, + fstar_secondary_upper_right, + fstar_secondary_lower_left, + fstar_secondary_lower_right, fstar_tmp1) lower_left_element = cache.mortars.neighbor_ids[1, mortar] lower_right_element = cache.mortars.neighbor_ids[2, mortar] @@ -1257,10 +1387,10 @@ end direction = 6 end end - surface_flux_values[:, :, :, direction, upper_left_element] .= fstar_upper_left - surface_flux_values[:, :, :, direction, upper_right_element] .= fstar_upper_right - surface_flux_values[:, :, :, direction, lower_left_element] .= fstar_lower_left - surface_flux_values[:, :, :, direction, lower_right_element] .= fstar_lower_right + surface_flux_values[:, :, :, direction, upper_left_element] .= fstar_primary_upper_left + surface_flux_values[:, :, :, direction, upper_right_element] .= fstar_primary_upper_right + surface_flux_values[:, :, :, direction, lower_left_element] .= fstar_primary_lower_left + surface_flux_values[:, :, :, direction, lower_right_element] .= fstar_primary_lower_right # Project small fluxes to large element if cache.mortars.large_sides[mortar] == 1 # -> small elements on right side @@ -1290,19 +1420,19 @@ end multiply_dimensionwise!(view(surface_flux_values, :, :, :, direction, large_element), mortar_l2.reverse_lower, mortar_l2.reverse_upper, - fstar_upper_left, fstar_tmp1) + fstar_secondary_upper_left, fstar_tmp1) add_multiply_dimensionwise!(view(surface_flux_values, :, :, :, direction, large_element), mortar_l2.reverse_upper, mortar_l2.reverse_upper, - fstar_upper_right, fstar_tmp1) + fstar_secondary_upper_right, fstar_tmp1) add_multiply_dimensionwise!(view(surface_flux_values, :, :, :, direction, large_element), mortar_l2.reverse_lower, mortar_l2.reverse_lower, - fstar_lower_left, fstar_tmp1) + fstar_secondary_lower_left, fstar_tmp1) add_multiply_dimensionwise!(view(surface_flux_values, :, :, :, direction, large_element), mortar_l2.reverse_upper, mortar_l2.reverse_lower, - fstar_lower_right, fstar_tmp1) + fstar_secondary_lower_right, fstar_tmp1) return nothing end diff --git a/src/solvers/dgsem_tree/dg_3d_parabolic.jl b/src/solvers/dgsem_tree/dg_3d_parabolic.jl index 69956d5834..657f39bbff 100644 --- a/src/solvers/dgsem_tree/dg_3d_parabolic.jl +++ b/src/solvers/dgsem_tree/dg_3d_parabolic.jl @@ -726,37 +726,32 @@ function calc_mortar_flux!(surface_flux_values, surface_integral, dg::DG, cache) @unpack surface_flux = surface_integral @unpack u_lower_left, u_lower_right, u_upper_left, u_upper_right, orientations = cache.mortars - @unpack (fstar_upper_left_threaded, fstar_upper_right_threaded, - fstar_lower_left_threaded, fstar_lower_right_threaded, + @unpack (fstar_primary_upper_left_threaded, fstar_primary_upper_right_threaded, + fstar_primary_lower_left_threaded, fstar_primary_lower_right_threaded, fstar_tmp1_threaded) = cache @threaded for mortar in eachmortar(dg, cache) # Choose thread-specific pre-allocated container - fstar_upper_left = fstar_upper_left_threaded[Threads.threadid()] - fstar_upper_right = fstar_upper_right_threaded[Threads.threadid()] - fstar_lower_left = fstar_lower_left_threaded[Threads.threadid()] - fstar_lower_right = fstar_lower_right_threaded[Threads.threadid()] + fstar_upper_left = fstar_primary_upper_left_threaded[Threads.threadid()] + fstar_upper_right = fstar_primary_upper_right_threaded[Threads.threadid()] + fstar_lower_left = fstar_primary_lower_left_threaded[Threads.threadid()] + fstar_lower_right = fstar_primary_lower_right_threaded[Threads.threadid()] fstar_tmp1 = fstar_tmp1_threaded[Threads.threadid()] # Calculate fluxes orientation = orientations[mortar] calc_fstar!(fstar_upper_left, equations_parabolic, surface_flux, dg, - u_upper_left, mortar, - orientation) + u_upper_left, mortar, orientation) calc_fstar!(fstar_upper_right, equations_parabolic, surface_flux, dg, - u_upper_right, - mortar, orientation) + u_upper_right, mortar, orientation) calc_fstar!(fstar_lower_left, equations_parabolic, surface_flux, dg, - u_lower_left, mortar, - orientation) + u_lower_left, mortar, orientation) calc_fstar!(fstar_lower_right, equations_parabolic, surface_flux, dg, - u_lower_right, - mortar, orientation) + u_lower_right, mortar, orientation) mortar_fluxes_to_elements!(surface_flux_values, mesh, equations_parabolic, mortar_l2, dg, cache, - mortar, - fstar_upper_left, fstar_upper_right, + mortar, fstar_upper_left, fstar_upper_right, fstar_lower_left, fstar_lower_right, fstar_tmp1) end @@ -782,6 +777,95 @@ end return nothing end +@inline function mortar_fluxes_to_elements!(surface_flux_values, + mesh::TreeMesh{3}, + equations_parabolic::AbstractEquationsParabolic, + mortar_l2::LobattoLegendreMortarL2, + dg::DGSEM, cache, + mortar, + fstar_upper_left, fstar_upper_right, + fstar_lower_left, fstar_lower_right, + fstar_tmp1) + lower_left_element = cache.mortars.neighbor_ids[1, mortar] + lower_right_element = cache.mortars.neighbor_ids[2, mortar] + upper_left_element = cache.mortars.neighbor_ids[3, mortar] + upper_right_element = cache.mortars.neighbor_ids[4, mortar] + large_element = cache.mortars.neighbor_ids[5, mortar] + + # Copy flux small to small + if cache.mortars.large_sides[mortar] == 1 # -> small elements on right side + if cache.mortars.orientations[mortar] == 1 + # L2 mortars in x-direction + direction = 1 + elseif cache.mortars.orientations[mortar] == 2 + # L2 mortars in y-direction + direction = 3 + else # if cache.mortars.orientations[mortar] == 3 + # L2 mortars in z-direction + direction = 5 + end + else # large_sides[mortar] == 2 -> small elements on left side + if cache.mortars.orientations[mortar] == 1 + # L2 mortars in x-direction + direction = 2 + elseif cache.mortars.orientations[mortar] == 2 + # L2 mortars in y-direction + direction = 4 + else # if cache.mortars.orientations[mortar] == 3 + # L2 mortars in z-direction + direction = 6 + end + end + surface_flux_values[:, :, :, direction, upper_left_element] .= fstar_upper_left + surface_flux_values[:, :, :, direction, upper_right_element] .= fstar_upper_right + surface_flux_values[:, :, :, direction, lower_left_element] .= fstar_lower_left + surface_flux_values[:, :, :, direction, lower_right_element] .= fstar_lower_right + + # Project small fluxes to large element + if cache.mortars.large_sides[mortar] == 1 # -> small elements on right side + if cache.mortars.orientations[mortar] == 1 + # L2 mortars in x-direction + direction = 2 + elseif cache.mortars.orientations[mortar] == 2 + # L2 mortars in y-direction + direction = 4 + else # if cache.mortars.orientations[mortar] == 3 + # L2 mortars in z-direction + direction = 6 + end + else # large_sides[mortar] == 2 -> small elements on left side + if cache.mortars.orientations[mortar] == 1 + # L2 mortars in x-direction + direction = 1 + elseif cache.mortars.orientations[mortar] == 2 + # L2 mortars in y-direction + direction = 3 + else # if cache.mortars.orientations[mortar] == 3 + # L2 mortars in z-direction + direction = 5 + end + end + + multiply_dimensionwise!(view(surface_flux_values, :, :, :, direction, + large_element), + mortar_l2.reverse_lower, mortar_l2.reverse_upper, + fstar_upper_left, fstar_tmp1) + add_multiply_dimensionwise!(view(surface_flux_values, :, :, :, direction, + large_element), + mortar_l2.reverse_upper, mortar_l2.reverse_upper, + fstar_upper_right, fstar_tmp1) + add_multiply_dimensionwise!(view(surface_flux_values, :, :, :, direction, + large_element), + mortar_l2.reverse_lower, mortar_l2.reverse_lower, + fstar_lower_left, fstar_tmp1) + add_multiply_dimensionwise!(view(surface_flux_values, :, :, :, direction, + large_element), + mortar_l2.reverse_upper, mortar_l2.reverse_lower, + fstar_lower_right, fstar_tmp1) + + return nothing +end + # Calculate the gradient of the transformed variables function calc_gradient!(gradients, u_transformed, t, mesh::TreeMesh{3}, equations_parabolic, diff --git a/test/test_p4est_2d.jl b/test/test_p4est_2d.jl index 46dfa542f3..7a1a3f3ac0 100644 --- a/test/test_p4est_2d.jl +++ b/test/test_p4est_2d.jl @@ -296,8 +296,8 @@ end t = sol.t[end] u_ode = sol.u[end] du_ode = similar(u_ode) - # Larger values for allowed allocations due to usage of custom - # integrator which are not *recorded* for the methods from + # Larger values for allowed allocations due to usage of custom + # integrator which are not *recorded* for the methods from # OrdinaryDiffEq.jl # Corresponding issue: https://github.com/trixi-framework/Trixi.jl/issues/1877 @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 15000 @@ -514,8 +514,8 @@ end t = sol.t[end] u_ode = sol.u[end] du_ode = similar(u_ode) - # Larger values for allowed allocations due to usage of custom - # integrator which are not *recorded* for the methods from + # Larger values for allowed allocations due to usage of custom + # integrator which are not *recorded* for the methods from # OrdinaryDiffEq.jl # Corresponding issue: https://github.com/trixi-framework/Trixi.jl/issues/1877 @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 15000 @@ -617,18 +617,18 @@ end @trixi_testset "elixir_mhd_rotor.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_rotor.jl"), - l2=[0.45520942119628094, 0.8918052934746296, 0.8324715234698744, + l2=[0.4551839744017604, 0.8917986079085971, 0.832474072904728, 0.0, - 0.9801268321975156, 0.1047572273917542, - 0.15551326369065485, + 0.98015167453983, 0.10475978783943254, + 0.15551175906375883, 0.0, - 2.0604823850230503e-5], - linf=[10.194219681102396, 18.254409373691253, - 10.031954811617876, + 2.026208477271868e-5], + linf=[10.19496728149964, 18.23726813972206, + 10.04367783820621, 0.0, - 19.646870952133547, 1.39386796733875, 1.8725058402095736, + 19.63022306543678, 1.3952679820406384, 1.8716515525771589, 0.0, - 0.001619787048808356], + 0.0017266639582675424], tspan=(0.0, 0.02)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) diff --git a/test/test_t8code_2d.jl b/test/test_t8code_2d.jl index 1e50fba144..ddede0e737 100644 --- a/test/test_t8code_2d.jl +++ b/test/test_t8code_2d.jl @@ -336,17 +336,16 @@ end # This test is identical to the one in `test_p4est_2d.jl` besides minor # deviations in the expected error norms. @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_rotor.jl"), - l2=[0.4420732463420727, 0.8804644301158163, 0.8262542320734158, + l2=[0.4419337424073218, 0.8804938551016932, 0.8258185723720365, 0.0, - 0.9615023124248694, 0.10386709616933161, - 0.15403081916109138, + 0.961220188718187, 0.10397273631386837, 0.15408979488125943, 0.0, - 2.835066224683485e-5], - linf=[10.045486750338348, 17.998696851793447, 9.57580213608948, + 2.66769410449947e-5], + linf=[10.053140536236942, 18.17070117006211, 9.549208389448738, 0.0, - 19.431290734386764, 1.3821685025605288, 1.8186235976086789, + 19.676151923191583, 1.3896544044814965, 1.8153256887969416, 0.0, - 0.0023118793481168537], + 0.0022030404596184786], tspan=(0.0, 0.02)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) diff --git a/test/test_tree_2d_mhd.jl b/test/test_tree_2d_mhd.jl index cf70d0eef5..b7152cb8b7 100644 --- a/test/test_tree_2d_mhd.jl +++ b/test/test_tree_2d_mhd.jl @@ -82,26 +82,26 @@ end @trixi_testset "elixir_mhd_alfven_wave_mortar.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave_mortar.jl"), l2=[ - 3.7762324533854616e-6, - 1.5534623833573546e-6, - 1.4577234868196855e-6, - 1.7647724628707057e-6, - 1.4831911814574333e-6, - 1.456369119716533e-6, - 1.4115666913995062e-6, - 1.804758237422838e-6, - 8.320469738087189e-7 + 1.0896015330565795e-5, + 4.152763046029908e-6, + 3.851874655132384e-6, + 4.2295110232831874e-6, + 3.135859402264645e-6, + 3.29531401471973e-6, + 3.1347238307092746e-6, + 4.186230495566739e-6, + 1.670859989962532e-6 ], linf=[ - 3.670661330201774e-5, - 1.530289442645827e-5, - 1.3592183785327006e-5, - 1.5173897443654383e-5, - 9.43771379136038e-6, - 1.0906323046233624e-5, - 1.0603954940346938e-5, - 1.5900499596113726e-5, - 5.978772247650426e-6 + 5.3178641410078775e-5, + 3.09217107711951e-5, + 2.7722788709688695e-5, + 2.1631700804783383e-5, + 1.558520409727926e-5, + 1.73873627985488e-5, + 1.6635747942750356e-5, + 2.0751205947883156e-5, + 7.655540399230342e-6 ], tspan=(0.0, 1.0)) # Ensure that we do not have excessive memory allocations @@ -181,26 +181,26 @@ end @trixi_testset "elixir_mhd_orszag_tang.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_orszag_tang.jl"), l2=[ - 0.21967600768935716, - 0.2643126515795721, - 0.31488287201980875, + 0.21970081242543155, + 0.2643108041773596, + 0.31484243079966445, 0.0, - 0.5160141621186931, - 0.23028914748088603, - 0.34413527376463915, + 0.5159994161306146, + 0.23024218609799854, + 0.34413704351228147, 0.0, - 0.003178793090381426 + 0.003220120866497733 ], linf=[ - 1.2749969218080568, - 0.6737013368774057, - 0.8604154399895696, + 1.2753954566712156, + 0.6737923290533722, + 0.8574465081172007, 0.0, - 2.799342099887639, - 0.6473347557712643, - 0.9691773375490476, + 2.800507621357904, + 0.6472414758680339, + 0.9707631523292184, 0.0, - 0.05729832038724348 + 0.06528658804650658 ], tspan=(0.0, 0.09)) # Ensure that we do not have excessive memory allocations @@ -216,26 +216,26 @@ end @trixi_testset "elixir_mhd_orszag_tang.jl with flux_hlle" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_orszag_tang.jl"), l2=[ - 0.10806619664693064, - 0.20199136742199922, - 0.22984589847526207, + 0.10806640059794005, + 0.20199169830949384, + 0.22984626162122923, 0.0, - 0.29950152196422647, - 0.15688413207147794, - 0.24293641543490646, + 0.2995035634728381, + 0.1568851137962238, + 0.24293639539810255, 0.0, - 0.003246181006326598 + 0.003246131507524401 ], linf=[ - 0.560316034595759, - 0.5095520363866776, - 0.6536748458764621, + 0.5600698267839397, + 0.5095520220558266, + 0.6536747859174317, 0.0, - 0.9627447086204038, - 0.3981375420906146, - 0.673472146198816, + 0.9624343226044095, + 0.39814285051228965, + 0.6734722065677001, 0.0, - 0.04879208429337193 + 0.048789764358224214 ], tspan=(0.0, 0.06), surface_flux=(flux_hlle, @@ -290,26 +290,26 @@ end @trixi_testset "elixir_mhd_rotor.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_rotor.jl"), l2=[ - 1.2623319195262743, - 1.8273050553090515, - 1.7004151198284634, + 1.264189543599029, + 1.8320601832078407, + 1.700447583152504, 0.0, - 2.2978570581460818, - 0.2147235065899803, - 0.23558337696054493, + 2.3024199507805165, + 0.21477383173627232, + 0.23559923070707714, 0.0, - 0.0032515115395693483 + 0.0034025828879598176 ], linf=[ - 11.003677581472843, - 14.70614192714736, - 15.687648666952708, + 10.988505627764773, + 14.712395261659752, + 15.687199838635722, 0.0, - 17.098104835553823, - 1.3283750501377847, - 1.4365828094434892, + 17.095921959435447, + 1.335014119480973, + 1.4366904817630641, 0.0, - 0.07886241196068537 + 0.08464617851256993 ], tspan=(0.0, 0.05)) # Ensure that we do not have excessive memory allocations @@ -325,26 +325,26 @@ end @trixi_testset "elixir_mhd_blast_wave.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_blast_wave.jl"), l2=[ - 0.17646728395490927, - 3.866230215339417, - 2.4867304651291255, + 0.17638656371490055, + 3.8616031530927084, + 2.4810236453809127, 0.0, - 355.4562971958441, - 2.359493623565687, - 1.4030741420730297, + 354.6341111396657, + 2.353681534580767, + 1.3926633033090652, 0.0, - 0.029613599942667133 + 0.030696738560246576 ], linf=[ - 1.581630420824181, - 44.15725488910748, - 13.056964982196554, + 1.5823311254590813, + 44.156859286717044, + 13.036736942960012, 0.0, - 2244.875490238186, - 13.07679044647926, - 9.14612176426092, + 2187.5906984085345, + 12.552321899505023, + 9.147117303057248, 0.0, - 0.5154756722488522 + 0.5285917066723818 ], tspan=(0.0, 0.003), # Calling the AnalysisCallback before iteration 9 causes the interpolation @@ -391,8 +391,8 @@ end t = sol.t[end] u_ode = sol.u[end] du_ode = similar(u_ode) - # Larger values for allowed allocations due to usage of custom - # integrator which are not *recorded* for the methods from + # Larger values for allowed allocations due to usage of custom + # integrator which are not *recorded* for the methods from # OrdinaryDiffEq.jl # Corresponding issue: https://github.com/trixi-framework/Trixi.jl/issues/1877 @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 15000 diff --git a/test/test_tree_2d_mhdmulti.jl b/test/test_tree_2d_mhdmulti.jl index d36554a667..19eaa28cb6 100644 --- a/test/test_tree_2d_mhdmulti.jl +++ b/test/test_tree_2d_mhdmulti.jl @@ -108,16 +108,16 @@ end @trixi_testset "elixir_mhdmulti_rotor.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_rotor.jl"), - l2=[0.6574605535168556, 0.6623234319361953, 0.0, - 0.689806698245354, - 0.04883686128677976, 0.08382459729494686, 0.0, - 0.0021114516459281177, - 0.15909290019096098, 0.07954645009548049], - linf=[9.362339085941425, 9.169838118652539, 0.0, - 10.600957847359556, - 0.6628317732399827, 1.4185626901435056, 0.0, - 0.06914316292003836, - 3.328770801731456, 1.664385400865728], + l2=[0.6574459522153201, 0.6620356383023878, 0.0, + 0.6888912144519942, + 0.04882939911229928, 0.08366520368549821, 0.0, + 0.0021850987869278136, + 0.15909935226497424, 0.07954967613248712], + linf=[9.363623690550916, 9.178740037372911, 0.0, + 10.611054196904469, + 0.6628358023789442, 1.419291349928299, 0.0, + 0.0988733910381692, + 3.3287658922602334, 1.6643829461301167], tspan=(0.0, 0.01)) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) diff --git a/test/test_tree_3d_mhd.jl b/test/test_tree_3d_mhd.jl index 98016cb519..4822d2ff60 100644 --- a/test/test_tree_3d_mhd.jl +++ b/test/test_tree_3d_mhd.jl @@ -152,26 +152,26 @@ end @trixi_testset "elixir_mhd_alfven_wave_mortar.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave_mortar.jl"), l2=[ - 0.001879021634926363, - 0.007032724521848316, - 0.0032793932234187325, - 0.009056594733320348, - 0.007514150120617965, - 0.007328739509868727, - 0.00309794018112387, - 0.009026356949274878, - 0.0035732583778049776 + 0.002117092205724962, + 0.0082287162318041, + 0.0034356818644221947, + 0.009802676239657889, + 0.008065655848544878, + 0.00822223011240085, + 0.0033142782650662905, + 0.009782724705061424, + 0.003818346240751859 ], linf=[ - 0.013734346970999622, - 0.06173467158736011, - 0.02183946452704291, - 0.06258216169457917, - 0.03672304497348122, - 0.055120532123884625, - 0.018202716205672487, - 0.06133688282205586, - 0.019888161885935608 + 0.01498490966057986, + 0.1168609357063561, + 0.024026660552548984, + 0.09909160811731985, + 0.06507407924731945, + 0.09999894905805326, + 0.029292103110423517, + 0.09399116188535625, + 0.031263077562076205 ], tspan=(0.0, 0.25)) # Ensure that we do not have excessive memory allocations