Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix noncons mortars for remaining mesh types #2134

Merged
merged 12 commits into from
Nov 14, 2024
55 changes: 36 additions & 19 deletions src/solvers/dgsem_p4est/dg_2d.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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]
Expand All @@ -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)
Expand All @@ -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,
Expand All @@ -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,
Expand All @@ -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
amrueda marked this conversation as resolved.
Show resolved Hide resolved

# 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
Expand All @@ -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
Expand Down
24 changes: 13 additions & 11 deletions src/solvers/dgsem_p4est/dg_2d_parabolic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand All @@ -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.
Expand Down Expand Up @@ -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)
Expand All @@ -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
Expand All @@ -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,
Expand All @@ -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!`
Expand Down
6 changes: 3 additions & 3 deletions src/solvers/dgsem_p4est/dg_2d_parallel.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down
Loading
Loading