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.

Sign up for GitHub

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

Jump to bottom

Fix non-conservative mortars for P4estMesh 3D #2127

Open
wants to merge 6 commits into
base: main
Choose a base branch
from
+123 −74
Open

Fix non-conservative mortars for P4estMesh 3D #2127

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
82c7ff1
Added non-unique mortar fluxes to be able to handle non-conservative …
amrueda Oct 22, 2024
0553a80
format
amrueda Oct 22, 2024
9fd429a
Updated p4est 3d tests
amrueda Oct 22, 2024
2ab09df
Updated parabolic and parallel p4est 3d implementations for new namin…
amrueda Oct 22, 2024
35aaf4e
format
amrueda Oct 22, 2024
f19d4d8
Fixed problem with parabolic non-conforming discretization
amrueda Oct 23, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
64 changes: 42 additions & 22 deletions src/solvers/dgsem_p4est/dg_3d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,9 +10,14 @@
function create_cache(mesh::Union{P4estMesh{3}, T8codeMesh{3}}, equations,
mortar_l2::LobattoLegendreMortarL2, uEltype)
# TODO: Taal compare performance of different types
fstar_threaded = [Array{uEltype, 4}(undef, nvariables(equations), nnodes(mortar_l2),
nnodes(mortar_l2), 4)
for _ in 1:Threads.nthreads()]
fstar_primary_threaded = [Array{uEltype, 4}(undef, nvariables(equations),
nnodes(mortar_l2),
nnodes(mortar_l2), 4)
for _ in 1:Threads.nthreads()]
fstar_secondary_threaded = [Array{uEltype, 4}(undef, nvariables(equations),
nnodes(mortar_l2),
nnodes(mortar_l2), 4)
for _ in 1:Threads.nthreads()]

fstar_tmp_threaded = [Array{uEltype, 3}(undef, nvariables(equations),
nnodes(mortar_l2), nnodes(mortar_l2))
Expand All @@ -21,7 +26,7 @@ function create_cache(mesh::Union{P4estMesh{3}, T8codeMesh{3}}, equations,
nnodes(mortar_l2))
for _ in 1:Threads.nthreads()]

(; fstar_threaded, fstar_tmp_threaded, u_threaded)
(; fstar_primary_threaded, fstar_secondary_threaded, fstar_tmp_threaded, u_threaded)
end

# index_to_start_step_3d(index::Symbol, index_range)
Expand Down Expand Up @@ -521,12 +526,13 @@ 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_threaded, fstar_tmp_threaded = cache
@unpack fstar_primary_threaded, fstar_secondary_threaded, fstar_tmp_threaded = cache
index_range = eachnode(dg)

@threaded for mortar in eachmortar(dg, cache)
# Choose thread-specific pre-allocated container
fstar = fstar_threaded[Threads.threadid()]
fstar_primary = fstar_primary_threaded[Threads.threadid()]
fstar_secondary = fstar_secondary_threaded[Threads.threadid()]
fstar_tmp = fstar_tmp_threaded[Threads.threadid()]

# Get index information on the small elements
Expand Down Expand Up @@ -555,7 +561,8 @@ function calc_mortar_flux!(surface_flux_values,
i_small, j_small, k_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,
i, j)
Expand All @@ -581,14 +588,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, fstar_tmp)
mortar, fstar_primary, fstar_secondary, u_buffer,
fstar_tmp)
end

return nothing
end

# Inlined version of the mortar flux computation on small elements for conservation fluxes
@inline function calc_mortar_flux!(fstar,
@inline function calc_mortar_flux!(fstar_primary, fstar_secondary,
mesh::Union{P4estMesh{3}, T8codeMesh{3}},
nonconservative_terms::False, equations,
surface_integral, dg::DG, cache,
Expand All @@ -603,13 +611,15 @@ end
flux = surface_flux(u_ll, u_rr, normal_direction, equations)

# Copy flux to buffer
set_node_vars!(fstar, flux, equations, dg, i_node_index, j_node_index,
set_node_vars!(fstar_primary, flux, equations, dg, i_node_index, j_node_index,
position_index)
set_node_vars!(fstar_secondary, flux, equations, dg, i_node_index, j_node_index,
position_index)
end

# Inlined version of the mortar flux computation on small elements for conservation fluxes
# with nonconservative terms
@inline function calc_mortar_flux!(fstar,
@inline function calc_mortar_flux!(fstar_primary, fstar_secondary,
mesh::Union{P4estMesh{3}, T8codeMesh{3}},
nonconservative_terms::True, equations,
surface_integral, dg::DG, cache,
Expand All @@ -627,20 +637,28 @@ end
# Compute nonconservative flux and add it to the flux 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)
flux_plus_noncons = flux + 0.5f0 * noncons
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_primary = flux + 0.5f0 * noncons_primary
flux_plus_noncons_secondary = flux + 0.5f0 * noncons_secondary

# Copy to buffer
set_node_vars!(fstar, flux_plus_noncons, equations, dg, i_node_index, j_node_index,
set_node_vars!(fstar_primary, flux_plus_noncons_primary, equations, dg,
i_node_index,
j_node_index,
position_index)
set_node_vars!(fstar_secondary, flux_plus_noncons_secondary, equations, dg,
i_node_index,
j_node_index,
position_index)
end

@inline function mortar_fluxes_to_elements!(surface_flux_values,
mesh::Union{P4estMesh{3}, T8codeMesh{3}},
equations,
mortar_l2::LobattoLegendreMortarL2,
dg::DGSEM, cache, mortar, fstar, u_buffer,
fstar_tmp)
dg::DGSEM, cache, mortar, fstar_primary,
fstar_secondary, u_buffer, fstar_tmp)
@unpack neighbor_ids, node_indices = cache.mortars
index_range = eachnode(dg)

Expand All @@ -652,28 +670,30 @@ end
element = neighbor_ids[position, mortar]
for j in eachnode(dg), i in eachnode(dg)
for v in eachvariable(equations)
surface_flux_values[v, i, j, small_direction, element] = fstar[v, i, j,
position]
surface_flux_values[v, i, j, small_direction, element] = fstar_primary[v,
i,
j,
position]
end
end
end

# Project small fluxes to large element.
multiply_dimensionwise!(u_buffer,
mortar_l2.reverse_lower, mortar_l2.reverse_lower,
view(fstar, .., 1),
view(fstar_secondary, .., 1),
fstar_tmp)
add_multiply_dimensionwise!(u_buffer,
mortar_l2.reverse_upper, mortar_l2.reverse_lower,
view(fstar, .., 2),
view(fstar_secondary, .., 2),
fstar_tmp)
add_multiply_dimensionwise!(u_buffer,
mortar_l2.reverse_lower, mortar_l2.reverse_upper,
view(fstar, .., 3),
view(fstar_secondary, .., 3),
fstar_tmp)
add_multiply_dimensionwise!(u_buffer,
mortar_l2.reverse_upper, mortar_l2.reverse_upper,
view(fstar, .., 4),
view(fstar_secondary, .., 4),
fstar_tmp)

# The flux is calculated in the outward direction of the small elements,
Expand Down
29 changes: 17 additions & 12 deletions src/solvers/dgsem_p4est/dg_3d_parabolic.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -271,7 +271,8 @@ end
mesh::P4estMesh{3},
equations::AbstractEquationsParabolic,
mortar_l2::LobattoLegendreMortarL2,
dg::DGSEM, cache, mortar, fstar, u_buffer,
dg::DGSEM, cache, mortar, fstar_primary,
fstar_secondary, u_buffer,
fstar_tmp)
@unpack neighbor_ids, node_indices = cache.mortars
index_range = eachnode(dg)
Expand All @@ -283,28 +284,30 @@ end
element = neighbor_ids[position, mortar]
for j in eachnode(dg), i in eachnode(dg)
for v in eachvariable(equations)
surface_flux_values[v, i, j, small_direction, element] = fstar[v, i, j,
position]
surface_flux_values[v, i, j, small_direction, element] = fstar_primary[v,
i,
j,
position]
end
end
end

# Project small fluxes to large element.
multiply_dimensionwise!(u_buffer,
mortar_l2.reverse_lower, mortar_l2.reverse_lower,
view(fstar, .., 1),
view(fstar_secondary, .., 1),
fstar_tmp)
add_multiply_dimensionwise!(u_buffer,
mortar_l2.reverse_upper, mortar_l2.reverse_lower,
view(fstar, .., 2),
view(fstar_secondary, .., 2),
fstar_tmp)
add_multiply_dimensionwise!(u_buffer,
mortar_l2.reverse_lower, mortar_l2.reverse_upper,
view(fstar, .., 3),
view(fstar_secondary, .., 3),
fstar_tmp)
add_multiply_dimensionwise!(u_buffer,
mortar_l2.reverse_upper, mortar_l2.reverse_upper,
view(fstar, .., 4),
view(fstar_secondary, .., 4),
fstar_tmp)

# The flux is calculated in the outward direction of the small elements,
Expand Down Expand Up @@ -788,12 +791,12 @@ 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_threaded, fstar_tmp_threaded = cache
@unpack fstar_primary_threaded, fstar_tmp_threaded = cache
index_range = eachnode(dg)

@threaded for mortar in eachmortar(dg, cache)
# Choose thread-specific pre-allocated container
fstar = fstar_threaded[Threads.threadid()]
fstar = fstar_primary_threaded[Threads.threadid()]
fstar_tmp = fstar_tmp_threaded[Threads.threadid()]

# Get index information on the small elements
Expand Down Expand Up @@ -842,7 +845,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, fstar_tmp)
mortar, fstar, fstar, u_buffer, fstar_tmp)
end

return nothing
Expand All @@ -851,7 +854,7 @@ end
# NOTE: Use analogy to "calc_mortar_flux!" for hyperbolic eqs with no nonconservative terms.
# Reasoning: "calc_interface_flux!" for parabolic part is implemented as the version for
# hyperbolic terms with conserved terms only, i.e., no nonconservative terms.
@inline function calc_mortar_flux!(fstar,
@inline function calc_mortar_flux!(fstar_primary, fstar_secondary,
mesh::P4estMesh{3},
nonconservative_terms::False,
equations::AbstractEquationsParabolic,
Expand All @@ -867,7 +870,9 @@ end
# TODO: parabolic; only BR1 at the moment
flux_ = 0.5f0 * (u_ll + u_rr)
# Copy flux to buffer
set_node_vars!(fstar, flux_, equations, dg, i_node_index, j_node_index,
set_node_vars!(fstar_primary, flux_, equations, dg, i_node_index, j_node_index,
position_index)
set_node_vars!(fstar_secondary, flux_, equations, dg, i_node_index, j_node_index,
position_index)
end

Expand Down
4 changes: 2 additions & 2 deletions src/solvers/dgsem_p4est/dg_3d_parallel.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -384,12 +384,12 @@ 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_threaded, fstar_tmp_threaded = cache
@unpack fstar_primary_threaded, fstar_tmp_threaded = cache
index_range = eachnode(dg)

@threaded for mortar in eachmpimortar(dg, cache)
# Choose thread-specific pre-allocated container
fstar = fstar_threaded[Threads.threadid()]
fstar = fstar_primary_threaded[Threads.threadid()]
fstar_tmp = fstar_tmp_threaded[Threads.threadid()]

# Get index information on the small elements
Expand Down
Loading
Loading