Skip to content

Commit 3580e95

Browse files
MarcoArtianoranocha
MarcoArtiano
and
ranocha
authored
An optimized 3D kernel for non-conservative fluxes (#2663)
* first commit * fmt * add tests * fix typo * add news entry * fix typo * fix test * add muladd and fix kernel * fix import ln_mean * add entry in AUTHORS.md * Apply suggestions from code review Co-authored-by: Hendrik Ranocha <[email protected]> --------- Co-authored-by: Hendrik Ranocha <[email protected]>
1 parent 2614a5f commit 3580e95

File tree

6 files changed

+423
-1
lines changed

6 files changed

+423
-1
lines changed

‎AUTHORS.md‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -24,6 +24,7 @@ provided substantial additions or modifications. Together, these two groups form
2424
The following people contributed major additions or modifications to Trixi.jl and
2525
are listed in alphabetical order:
2626

27+
* Marco Artiano
2728
* Maximilian D. Bertrand
2829
* Benjamin Bolm
2930
* Simon Candelaresi

‎NEWS.md‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -13,7 +13,7 @@ for human readability.
1313
In the new version, IDP positivity limiting for conservative variables (using
1414
the keyword `positivity_variables_cons` in `SubcellLimiterIDP()`) is supported.
1515
`BoundsCheckCallback` is not supported in 3D yet.
16-
- Optimized 2D kernel for nonconservative fluxes with `P4estMesh` was added ([#2653]).
16+
- Optimized 2D and 3D kernels for nonconservative fluxes with `P4estMesh` were added ([#2653], [#2663]).
1717
The optimized kernel can be enabled via the trait `Trixi.combine_conservative_and_nonconservative_fluxes(flux, equations)`.
1818
When the trait is set to `Trixi.True()`, a single method has to be defined, that computes and returns the tuple
1919
`flux_cons(u_ll, u_rr) + 0.5f0 * flux_noncons(u_ll, u_rr)` and

‎src/solvers/dgsem_p4est/dg_3d.jl‎

Lines changed: 112 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -291,6 +291,29 @@ end
291291
secondary_i_node_index, secondary_j_node_index,
292292
secondary_direction_index,
293293
secondary_element_index)
294+
calc_interface_flux!(surface_flux_values, mesh, have_nonconservative_terms,
295+
combine_conservative_and_nonconservative_fluxes(surface_integral.surface_flux,
296+
equations),
297+
equations, surface_integral, dg, cache, interface_index,
298+
normal_direction, primary_i_node_index, primary_j_node_index,
299+
primary_direction_index, primary_element_index,
300+
secondary_i_node_index, secondary_j_node_index,
301+
secondary_direction_index, secondary_element_index)
302+
return nothing
303+
end
304+
305+
@inline function calc_interface_flux!(surface_flux_values,
306+
mesh::Union{P4estMesh{3}, T8codeMesh{3}},
307+
have_nonconservative_terms::True,
308+
combine_conservative_and_nonconservative_fluxes::False,
309+
equations,
310+
surface_integral, dg::DG, cache,
311+
interface_index, normal_direction,
312+
primary_i_node_index, primary_j_node_index,
313+
primary_direction_index, primary_element_index,
314+
secondary_i_node_index, secondary_j_node_index,
315+
secondary_direction_index,
316+
secondary_element_index)
294317
@unpack u = cache.interfaces
295318
surface_flux, nonconservative_flux = surface_integral.surface_flux
296319

@@ -320,6 +343,36 @@ end
320343
return nothing
321344
end
322345

346+
@inline function calc_interface_flux!(surface_flux_values,
347+
mesh::Union{P4estMesh{3}, T8codeMesh{3}},
348+
have_nonconservative_terms::True,
349+
combine_conservative_and_nonconservative_fluxes::True,
350+
equations,
351+
surface_integral, dg::DG, cache,
352+
interface_index, normal_direction,
353+
primary_i_node_index, primary_j_node_index,
354+
primary_direction_index, primary_element_index,
355+
secondary_i_node_index, secondary_j_node_index,
356+
secondary_direction_index,
357+
secondary_element_index)
358+
@unpack u = cache.interfaces
359+
@unpack surface_flux = surface_integral
360+
u_ll, u_rr = get_surface_node_vars(u, equations, dg, primary_i_node_index,
361+
primary_j_node_index, interface_index)
362+
363+
flux_left, flux_right = surface_flux(u_ll, u_rr, normal_direction, equations)
364+
365+
# Store the flux with nonconservative terms on the primary and secondary elements
366+
for v in eachvariable(equations)
367+
surface_flux_values[v, primary_i_node_index, primary_j_node_index,
368+
primary_direction_index, primary_element_index] = flux_left[v]
369+
surface_flux_values[v, secondary_i_node_index, secondary_j_node_index,
370+
secondary_direction_index, secondary_element_index] = -flux_right[v]
371+
end
372+
373+
return nothing
374+
end
375+
323376
function prolong2boundaries!(cache, u,
324377
mesh::Union{P4estMesh{3}, T8codeMesh{3}},
325378
equations, surface_integral, dg::DG)
@@ -452,6 +505,26 @@ end
452505
k_index, i_node_index, j_node_index,
453506
direction_index,
454507
element_index, boundary_index)
508+
calc_boundary_flux!(surface_flux_values, t, boundary_condition, mesh,
509+
nonconservative_terms,
510+
combine_conservative_and_nonconservative_fluxes(surface_integral.surface_flux,
511+
equations),
512+
equations,
513+
surface_integral, dg, cache,
514+
i_index, j_index, k_index, i_node_index, j_node_index,
515+
direction_index, element_index, boundary_index)
516+
return nothing
517+
end
518+
519+
@inline function calc_boundary_flux!(surface_flux_values, t, boundary_condition,
520+
mesh::Union{P4estMesh{3}, T8codeMesh{3}},
521+
nonconservative_terms::True,
522+
combine_conservative_and_nonconservative_fluxes::False,
523+
equations,
524+
surface_integral, dg::DG, cache, i_index, j_index,
525+
k_index, i_node_index, j_node_index,
526+
direction_index,
527+
element_index, boundary_index)
455528
@unpack boundaries = cache
456529
@unpack node_coordinates, contravariant_vectors = cache.elements
457530
@unpack surface_flux = surface_integral
@@ -486,6 +559,45 @@ end
486559
return nothing
487560
end
488561

562+
@inline function calc_boundary_flux!(surface_flux_values, t, boundary_condition,
563+
mesh::Union{P4estMesh{3}, T8codeMesh{3}},
564+
nonconservative_terms::True,
565+
combine_conservative_and_nonconservative_fluxes::True,
566+
equations,
567+
surface_integral, dg::DG, cache, i_index, j_index,
568+
k_index, i_node_index, j_node_index,
569+
direction_index,
570+
element_index, boundary_index)
571+
@unpack boundaries = cache
572+
@unpack node_coordinates, contravariant_vectors = cache.elements
573+
@unpack surface_flux = surface_integral
574+
575+
# Extract solution data from boundary container
576+
u_inner = get_node_vars(boundaries.u, equations, dg, i_node_index, j_node_index,
577+
boundary_index)
578+
579+
# Outward-pointing normal direction (not normalized)
580+
normal_direction = get_normal_direction(direction_index, contravariant_vectors,
581+
i_index, j_index, k_index, element_index)
582+
583+
# Coordinates at boundary node
584+
x = get_node_coords(node_coordinates, equations, dg,
585+
i_index, j_index, k_index, element_index)
586+
587+
# Call pointwise numerical flux functions for the conservative and nonconservative part
588+
# in the normal direction on the boundary
589+
flux = boundary_condition(u_inner, normal_direction, x, t,
590+
surface_flux, equations)
591+
592+
# Copy flux to element storage in the correct orientation
593+
for v in eachvariable(equations)
594+
surface_flux_values[v, i_node_index, j_node_index,
595+
direction_index, element_index] = flux[v]
596+
end
597+
598+
return nothing
599+
end
600+
489601
function prolong2mortars!(cache, u,
490602
mesh::Union{P4estMesh{3}, T8codeMesh{3}}, equations,
491603
mortar_l2::LobattoLegendreMortarL2,

‎src/solvers/dgsem_structured/dg_3d.jl‎

Lines changed: 107 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -159,6 +159,22 @@ end
159159
T8codeMesh{3}},
160160
have_nonconservative_terms::True, equations,
161161
volume_flux, dg::DGSEM, cache, alpha = true)
162+
flux_differencing_kernel!(du, u, element, mesh, have_nonconservative_terms,
163+
combine_conservative_and_nonconservative_fluxes(volume_flux,
164+
equations),
165+
equations, volume_flux, dg, cache, alpha)
166+
167+
return nothing
168+
end
169+
170+
@inline function flux_differencing_kernel!(du, u,
171+
element,
172+
mesh::Union{StructuredMesh{3}, P4estMesh{3},
173+
T8codeMesh{3}},
174+
have_nonconservative_terms::True,
175+
combine_conservative_and_nonconservative_fluxes::False,
176+
equations,
177+
volume_flux, dg::DGSEM, cache, alpha = true)
162178
@unpack derivative_split = dg.basis
163179
@unpack contravariant_vectors = cache.elements
164180
symmetric_flux, nonconservative_flux = volume_flux
@@ -242,6 +258,97 @@ end
242258
return nothing
243259
end
244260

261+
@inline function flux_differencing_kernel!(du, u,
262+
element,
263+
mesh::Union{StructuredMesh{3}, P4estMesh{3},
264+
T8codeMesh{3}},
265+
have_nonconservative_terms::True,
266+
combine_conservative_and_nonconservative_fluxes::True,
267+
equations,
268+
volume_flux, dg::DGSEM, cache, alpha = true)
269+
@unpack derivative_split = dg.basis
270+
@unpack contravariant_vectors = cache.elements
271+
272+
# Calculate volume integral in one element
273+
for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
274+
u_node = get_node_vars(u, equations, dg, i, j, k, element)
275+
276+
# pull the contravariant vectors in each coordinate direction
277+
Ja1_node = get_contravariant_vector(1, contravariant_vectors, i, j, k, element)
278+
Ja2_node = get_contravariant_vector(2, contravariant_vectors, i, j, k, element)
279+
Ja3_node = get_contravariant_vector(3, contravariant_vectors, i, j, k, element)
280+
281+
# All diagonal entries of `derivative_split` are zero. Thus, we can skip
282+
# the computation of the diagonal terms. In addition, we use the symmetry
283+
# of the `volume_flux` to save half of the possible two-point flux
284+
# computations.
285+
286+
# x direction
287+
for ii in (i + 1):nnodes(dg)
288+
u_node_ii = get_node_vars(u, equations, dg, ii, j, k, element)
289+
# pull the contravariant vectors and compute the average
290+
Ja1_node_ii = get_contravariant_vector(1, contravariant_vectors,
291+
ii, j, k, element)
292+
# average mapping terms in first coordinate direction,
293+
# used as normal vector in the flux computation
294+
Ja1_avg = 0.5f0 * (Ja1_node + Ja1_node_ii)
295+
# compute the contravariant sharp flux in the direction of the
296+
# averaged contravariant vector
297+
fluxtilde1_left, fluxtilde1_right = volume_flux(u_node, u_node_ii, Ja1_avg,
298+
equations)
299+
multiply_add_to_node_vars!(du, alpha * derivative_split[i, ii],
300+
fluxtilde1_left,
301+
equations, dg, i, j, k, element)
302+
multiply_add_to_node_vars!(du, alpha * derivative_split[ii, i],
303+
fluxtilde1_right,
304+
equations, dg, ii, j, k, element)
305+
end
306+
307+
# y direction
308+
for jj in (j + 1):nnodes(dg)
309+
u_node_jj = get_node_vars(u, equations, dg, i, jj, k, element)
310+
# pull the contravariant vectors and compute the average
311+
Ja2_node_jj = get_contravariant_vector(2, contravariant_vectors,
312+
i, jj, k, element)
313+
# average mapping terms in second coordinate direction,
314+
# used as normal vector in the flux computation
315+
Ja2_avg = 0.5f0 * (Ja2_node + Ja2_node_jj)
316+
# compute the contravariant sharp flux in the direction of the
317+
# averaged contravariant vector
318+
fluxtilde2_left, fluxtilde2_right = volume_flux(u_node, u_node_jj, Ja2_avg,
319+
equations)
320+
multiply_add_to_node_vars!(du, alpha * derivative_split[j, jj],
321+
fluxtilde2_left,
322+
equations, dg, i, j, k, element)
323+
multiply_add_to_node_vars!(du, alpha * derivative_split[jj, j],
324+
fluxtilde2_right,
325+
equations, dg, i, jj, k, element)
326+
end
327+
328+
# z direction
329+
for kk in (k + 1):nnodes(dg)
330+
u_node_kk = get_node_vars(u, equations, dg, i, j, kk, element)
331+
# pull the contravariant vectors and compute the average
332+
Ja3_node_kk = get_contravariant_vector(3, contravariant_vectors,
333+
i, j, kk, element)
334+
# average mapping terms in third coordinate direction,
335+
# used as normal vector in the flux computation
336+
Ja3_avg = 0.5f0 * (Ja3_node + Ja3_node_kk)
337+
# compute the contravariant sharp flux in the direction of the
338+
# averaged contravariant vector
339+
fluxtilde3_left, fluxtilde3_right = volume_flux(u_node, u_node_kk, Ja3_avg,
340+
equations)
341+
multiply_add_to_node_vars!(du, alpha * derivative_split[k, kk],
342+
fluxtilde3_left,
343+
equations, dg, i, j, k, element)
344+
multiply_add_to_node_vars!(du, alpha * derivative_split[kk, k],
345+
fluxtilde3_right,
346+
equations, dg, i, j, kk, element)
347+
end
348+
end
349+
return nothing
350+
end
351+
245352
# Computing the normal vector for the FV method on curvilinear subcells.
246353
# To fulfill free-stream preservation we use the explicit formula B.53 in Appendix B.4
247354
# by Hennemann, Rueda-Ramirez, Hindenlang, Gassner (2020)

‎test/test_performance_specializations_2d.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -172,6 +172,7 @@ end
172172
@test du_specialized ≈ du_baseline
173173
end
174174
end
175+
175176
@timed_testset "P4estMesh2D, combine_conservative_and_nonconservative_fluxes" begin
176177
trixi_include(@__MODULE__,
177178
joinpath(EXAMPLES_DIR, "p4est_2d_dgsem", "elixir_mhd_rotor.jl"),

0 commit comments

Comments
 (0)