Skip to content

Commit 3c43787

Browse files
MarcoArtianoranocha
MarcoArtiano
and
ranocha
authored
An optimized 2D kernel for non-conservative fluxes (#2653)
* format * fix names * remove some comments * change is to has_antisymmetric_flux * rename has_antisymmetric_flux to combine_conservative_and_nonconservative_fluxes * change comment for the trait * Apply suggestions from code review Co-authored-by: Hendrik Ranocha <[email protected]> * add tests * add test structured mesh 2d * fmt * fix test * fix fluxes in test * Apply suggestions from code review Co-authored-by: Hendrik Ranocha <[email protected]> * attempt to remove StructuredMesh specialization * add link to test for an example * change test with BC * add news entry * add 0.5 factor in the trait doc * update test * remove AMR from the test * fix test * fix test --------- Co-authored-by: Hendrik Ranocha <[email protected]>
1 parent 3fea822 commit 3c43787

File tree

5 files changed

+448
-1
lines changed

5 files changed

+448
-1
lines changed

‎NEWS.md‎

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -13,7 +13,11 @@ 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-
16+
- Optimized 2D kernel for nonconservative fluxes with `P4estMesh` was added ([#2653]).
17+
The optimized kernel can be enabled via the trait `Trixi.combine_conservative_and_nonconservative_fluxes(flux, equations)`.
18+
When the trait is set to `Trixi.True()`, a single method has to be defined, that computes and returns the tuple
19+
`flux_cons(u_ll, u_rr) + 0.5f0 * flux_noncons(u_ll, u_rr)` and
20+
`flux_cons(u_ll, u_rr) + 0.5f0 * flux_noncons(u_rr, u_ll)`.
1721

1822
## Changes when updating to v0.13 from v0.12.x
1923

‎src/equations/equations.jl‎

Lines changed: 36 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -326,6 +326,42 @@ combined with certain solvers (e.g., subcell limiting).
326326
"""
327327
function n_nonconservative_terms end
328328

329+
"""
330+
Trixi.combine_conservative_and_nonconservative_fluxes(flux, equations)
331+
332+
Trait function indicating whether the given `flux` and `equations` support
333+
fusing the computation of conservative fluxes with nonconservative fluxes.
334+
This is purely a performance optimization for equations with nonconservative
335+
terms (i.e., where `have_nonconservative_terms(equations)` is `Trixi.True()`).
336+
The default value is `Trixi.False()`, i.e., you have to pass a tuple of
337+
numerical fluxes for the conservative and the nonconservative terms, e.g.,
338+
to compute surface terms or the [`VolumeIntegralFluxDifferencing`](@ref).
339+
340+
For some systems and flux implementations, it is cheaper to compute
341+
342+
flux_noncons(u_ll, u_rr, orientation_or_normal_direction, equations)
343+
344+
and
345+
346+
flux_noncons(u_rr, u_ll, orientation_or_normal_direction, equations)
347+
348+
together, or to compute conservative and nonconservative flux contributions in
349+
a single fused kernel. In this case, you should set this trait to be `Trixi.True()`
350+
to take advantage of a more efficient implementation. In this case, you have to
351+
define a single method that computes
352+
353+
flux_cons(u_ll, u_rr, n, equations) + 0.5f0 * flux_noncons(u_ll, u_rr, n, equations)
354+
355+
and
356+
357+
flux_cons(u_ll, u_rr, n, equations) + 0.5f0 * flux_noncons(u_rr, u_ll, n, equations)
358+
359+
together and returns them as a tuple.
360+
See also the test section P4estMesh2D with combine_conservative_and_nonconservative_fluxes in
361+
[Test Performance](https://github.com/trixi-framework/Trixi.jl/blob/main/test/test_performance_specializations_2d.jl).
362+
"""
363+
combine_conservative_and_nonconservative_fluxes(flux, ::AbstractEquations) = False()
364+
329365
have_constant_speed(::AbstractEquations) = False()
330366

331367
"""

‎src/solvers/dgsem_p4est/dg_2d.jl‎

Lines changed: 108 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -224,6 +224,31 @@ end
224224
primary_element_index,
225225
secondary_node_index, secondary_direction_index,
226226
secondary_element_index)
227+
@unpack surface_flux = surface_integral
228+
calc_interface_flux!(surface_flux_values, mesh, have_nonconservative_terms,
229+
combine_conservative_and_nonconservative_fluxes(surface_flux,
230+
equations),
231+
equations,
232+
surface_integral, dg, cache,
233+
interface_index, normal_direction,
234+
primary_node_index, primary_direction_index,
235+
primary_element_index,
236+
secondary_node_index, secondary_direction_index,
237+
secondary_element_index)
238+
return nothing
239+
end
240+
241+
@inline function calc_interface_flux!(surface_flux_values,
242+
mesh::Union{P4estMesh{2}, T8codeMesh{2}},
243+
have_nonconservative_terms::True,
244+
combine_conservative_and_nonconservative_fluxes::False,
245+
equations,
246+
surface_integral, dg::DG, cache,
247+
interface_index, normal_direction,
248+
primary_node_index, primary_direction_index,
249+
primary_element_index,
250+
secondary_node_index, secondary_direction_index,
251+
secondary_element_index)
227252
@unpack u = cache.interfaces
228253
surface_flux, nonconservative_flux = surface_integral.surface_flux
229254

@@ -252,6 +277,33 @@ end
252277
return nothing
253278
end
254279

280+
@inline function calc_interface_flux!(surface_flux_values,
281+
mesh::Union{P4estMesh{2}, T8codeMesh{2}},
282+
have_nonconservative_terms::True,
283+
combine_conservative_and_nonconservative_fluxes::True,
284+
equations,
285+
surface_integral, dg::DG, cache,
286+
interface_index, normal_direction,
287+
primary_node_index, primary_direction_index,
288+
primary_element_index,
289+
secondary_node_index, secondary_direction_index,
290+
secondary_element_index)
291+
@unpack u = cache.interfaces
292+
@unpack surface_flux = surface_integral
293+
294+
u_ll, u_rr = get_surface_node_vars(u, equations, dg, primary_node_index,
295+
interface_index)
296+
297+
flux_left, flux_right = surface_flux(u_ll, u_rr, normal_direction, equations)
298+
299+
for v in eachvariable(equations)
300+
surface_flux_values[v, primary_node_index, primary_direction_index, primary_element_index] = flux_left[v]
301+
surface_flux_values[v, secondary_node_index, secondary_direction_index, secondary_element_index] = -flux_right[v]
302+
end
303+
304+
return nothing
305+
end
306+
255307
function prolong2boundaries!(cache, u,
256308
mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}},
257309
equations, surface_integral, dg::DG)
@@ -359,6 +411,26 @@ end
359411
i_index, j_index,
360412
node_index, direction_index, element_index,
361413
boundary_index)
414+
@unpack surface_flux = surface_integral
415+
calc_boundary_flux!(surface_flux_values, t, boundary_condition,
416+
mesh, have_nonconservative_terms,
417+
combine_conservative_and_nonconservative_fluxes(surface_flux,
418+
equations),
419+
equations, surface_integral, dg, cache,
420+
i_index, j_index,
421+
node_index, direction_index, element_index, boundary_index)
422+
return nothing
423+
end
424+
425+
@inline function calc_boundary_flux!(surface_flux_values, t, boundary_condition,
426+
mesh::Union{P4estMesh{2}, T8codeMesh{2}},
427+
have_nonconservative_terms::True,
428+
combine_conservative_and_nonconservative_fluxes::False,
429+
equations,
430+
surface_integral, dg::DG, cache,
431+
i_index, j_index,
432+
node_index, direction_index, element_index,
433+
boundary_index)
362434
@unpack boundaries = cache
363435
@unpack node_coordinates, contravariant_vectors = cache.elements
364436

@@ -391,6 +463,42 @@ end
391463
return nothing
392464
end
393465

466+
@inline function calc_boundary_flux!(surface_flux_values, t, boundary_condition,
467+
mesh::Union{P4estMesh{2}, T8codeMesh{2}},
468+
have_nonconservative_terms::True,
469+
combine_conservative_and_nonconservative_fluxes::True,
470+
equations,
471+
surface_integral, dg::DG, cache,
472+
i_index, j_index,
473+
node_index, direction_index, element_index,
474+
boundary_index)
475+
@unpack boundaries = cache
476+
@unpack node_coordinates, contravariant_vectors = cache.elements
477+
478+
# Extract solution data from boundary container
479+
u_inner = get_node_vars(boundaries.u, equations, dg, node_index, boundary_index)
480+
481+
# Outward-pointing normal direction (not normalized)
482+
normal_direction = get_normal_direction(direction_index, contravariant_vectors,
483+
i_index, j_index, element_index)
484+
485+
# Coordinates at boundary node
486+
x = get_node_coords(node_coordinates, equations, dg,
487+
i_index, j_index, element_index)
488+
489+
# Call pointwise numerical flux functions for the conservative and nonconservative part
490+
# in the normal direction on the boundary
491+
flux = boundary_condition(u_inner, normal_direction, x, t,
492+
surface_integral.surface_flux, equations)
493+
494+
# Copy flux to element storage in the correct orientation
495+
for v in eachvariable(equations)
496+
surface_flux_values[v, node_index, direction_index, element_index] = flux[v]
497+
end
498+
499+
return nothing
500+
end
501+
394502
function prolong2mortars!(cache, u,
395503
mesh::Union{P4estMesh{2}, P4estMeshView{2}, T8codeMesh{2}},
396504
equations,

‎src/solvers/dgsem_structured/dg_2d.jl‎

Lines changed: 109 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -126,6 +126,25 @@ end
126126
T8codeMesh{2}},
127127
have_nonconservative_terms::True, equations,
128128
volume_flux, dg::DGSEM, cache, alpha = true)
129+
flux_differencing_kernel!(du, u, element, mesh, have_nonconservative_terms,
130+
combine_conservative_and_nonconservative_fluxes(volume_flux,
131+
equations),
132+
equations,
133+
volume_flux,
134+
dg, cache, alpha)
135+
return nothing
136+
end
137+
138+
@inline function flux_differencing_kernel!(du, u,
139+
element,
140+
mesh::Union{StructuredMesh{2},
141+
StructuredMeshView{2},
142+
UnstructuredMesh2D, P4estMesh{2},
143+
T8codeMesh{2}},
144+
have_nonconservative_terms::True,
145+
combine_conservative_and_nonconservative_fluxes::False,
146+
equations,
147+
volume_flux, dg::DGSEM, cache, alpha = true)
129148
@unpack derivative_split = dg.basis
130149
@unpack contravariant_vectors = cache.elements
131150
symmetric_flux, nonconservative_flux = volume_flux
@@ -190,6 +209,78 @@ end
190209
return nothing
191210
end
192211

212+
@inline function flux_differencing_kernel!(du, u,
213+
element,
214+
mesh::Union{StructuredMesh{2},
215+
StructuredMeshView{2},
216+
UnstructuredMesh2D, P4estMesh{2},
217+
T8codeMesh{2}},
218+
have_nonconservative_terms::True,
219+
combine_conservative_and_nonconservative_fluxes::True,
220+
equations,
221+
volume_flux, dg::DGSEM, cache, alpha = true)
222+
@unpack derivative_split = dg.basis
223+
@unpack contravariant_vectors = cache.elements
224+
225+
# Calculate volume integral in one element
226+
for j in eachnode(dg), i in eachnode(dg)
227+
u_node = get_node_vars(u, equations, dg, i, j, element)
228+
229+
# pull the contravariant vectors in each coordinate direction
230+
Ja1_node = get_contravariant_vector(1, contravariant_vectors, i, j, element)
231+
Ja2_node = get_contravariant_vector(2, contravariant_vectors, i, j, element)
232+
233+
# All diagonal entries of `derivative_split` are zero. Thus, we can skip
234+
# the computation of the diagonal terms. In addition, we use the symmetry
235+
# of the `volume_flux` to save half of the possible two-point flux
236+
# computations.
237+
238+
# x direction
239+
for ii in (i + 1):nnodes(dg)
240+
u_node_ii = get_node_vars(u, equations, dg, ii, j, element)
241+
# pull the contravariant vectors and compute the average
242+
Ja1_node_ii = get_contravariant_vector(1, contravariant_vectors,
243+
ii, j, element)
244+
# average mapping terms in first coordinate direction,
245+
# used as normal vector in the flux computation
246+
Ja1_avg = 0.5f0 * (Ja1_node + Ja1_node_ii)
247+
# compute the contravariant sharp flux in the direction of the
248+
# averaged contravariant vector
249+
fluxtilde1_left, fluxtilde1_right = volume_flux(u_node, u_node_ii, Ja1_avg,
250+
equations)
251+
multiply_add_to_node_vars!(du, alpha * derivative_split[i, ii],
252+
fluxtilde1_left,
253+
equations, dg, i, j, element)
254+
multiply_add_to_node_vars!(du, alpha * derivative_split[ii, i],
255+
fluxtilde1_right,
256+
equations, dg, ii, j, element)
257+
end
258+
259+
# y direction
260+
for jj in (j + 1):nnodes(dg)
261+
u_node_jj = get_node_vars(u, equations, dg, i, jj, element)
262+
# pull the contravariant vectors and compute the average
263+
Ja2_node_jj = get_contravariant_vector(2, contravariant_vectors,
264+
i, jj, element)
265+
# average mapping terms in second coordinate direction,
266+
# used as normal vector in the flux computation
267+
Ja2_avg = 0.5f0 * (Ja2_node + Ja2_node_jj)
268+
# compute the contravariant sharp flux in the direction of the
269+
# averaged contravariant vector
270+
fluxtilde2_left, fluxtilde2_right = volume_flux(u_node, u_node_jj, Ja2_avg,
271+
equations)
272+
multiply_add_to_node_vars!(du, alpha * derivative_split[j, jj],
273+
fluxtilde2_left,
274+
equations, dg, i, j, element)
275+
multiply_add_to_node_vars!(du, alpha * derivative_split[jj, j],
276+
fluxtilde2_right,
277+
equations, dg, i, jj, element)
278+
end
279+
end
280+
281+
return nothing
282+
end
283+
193284
# Computing the normal vector for the FV method on curvilinear subcells.
194285
# To fulfill free-stream preservation we use the explicit formula B.53 in Appendix B.4
195286
# by Hennemann, Rueda-Ramirez, Hindenlang, Gassner (2020)
@@ -451,6 +542,24 @@ end
451542
StructuredMeshView{2}},
452543
have_nonconservative_terms::True, equations,
453544
surface_integral, dg::DG, cache)
545+
calc_interface_flux!(surface_flux_values, left_element, right_element, orientation,
546+
u, mesh, have_nonconservative_terms,
547+
combine_conservative_and_nonconservative_fluxes(surface_integral.surface_flux,
548+
equations),
549+
equations, surface_integral,
550+
dg, cache)
551+
return nothing
552+
end
553+
554+
@inline function calc_interface_flux!(surface_flux_values, left_element, right_element,
555+
orientation, u,
556+
mesh::Union{StructuredMesh{2},
557+
StructuredMeshView{2}},
558+
have_nonconservative_terms::True,
559+
combine_conservative_and_nonconservative_fluxes::False,
560+
equations,
561+
surface_integral, dg::DG, cache)
562+
454563
# See comment on `calc_interface_flux!` with `have_nonconservative_terms::False`
455564
if left_element <= 0 # left_element = 0 at boundaries
456565
return nothing

0 commit comments

Comments
 (0)