Skip to content

Commit

Permalink
Fix non-conservative mortars for P4estMesh 3D (trixi-framework#2127)
Browse files Browse the repository at this point in the history 
* Added non-unique mortar fluxes to be able to handle non-conservative equations

* format

* Updated p4est 3d tests

* Updated parabolic and parallel p4est 3d implementations for new naming convention

* format

* Fixed problem with parabolic non-conforming discretization
  • Loading branch information
@amrueda
amrueda authored Oct 25, 2024
1 parent 4a81e00 commit df9013a
Show file tree
Hide file tree
Showing 4 changed files with 123 additions and 74 deletions.
64 changes: 42 additions & 22 deletions src/solvers/dgsem_p4est/dg_3d.jl
View file
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
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
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

0 comments on commit df9013a

Please sign in to comment.