Skip to content

Commit df9013a

Browse files
amruedaamrueda
authored
Fix non-conservative mortars for P4estMesh 3D (trixi-framework#2127)
* 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
1 parent 4a81e00 commit df9013a

File tree

4 files changed

+123
-74
lines changed

4 files changed

+123
-74
lines changed

‎src/solvers/dgsem_p4est/dg_3d.jl

+42-22
Original file line numberDiff line numberDiff line change
@@ -10,9 +10,14 @@
1010
function create_cache(mesh::Union{P4estMesh{3}, T8codeMesh{3}}, equations,
1111
mortar_l2::LobattoLegendreMortarL2, uEltype)
1212
# TODO: Taal compare performance of different types
13-
fstar_threaded = [Array{uEltype, 4}(undef, nvariables(equations), nnodes(mortar_l2),
14-
nnodes(mortar_l2), 4)
15-
for _ in 1:Threads.nthreads()]
13+
fstar_primary_threaded = [Array{uEltype, 4}(undef, nvariables(equations),
14+
nnodes(mortar_l2),
15+
nnodes(mortar_l2), 4)
16+
for _ in 1:Threads.nthreads()]
17+
fstar_secondary_threaded = [Array{uEltype, 4}(undef, nvariables(equations),
18+
nnodes(mortar_l2),
19+
nnodes(mortar_l2), 4)
20+
for _ in 1:Threads.nthreads()]
1621

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

24-
(; fstar_threaded, fstar_tmp_threaded, u_threaded)
29+
(; fstar_primary_threaded, fstar_secondary_threaded, fstar_tmp_threaded, u_threaded)
2530
end
2631

2732
# index_to_start_step_3d(index::Symbol, index_range)
@@ -521,12 +526,13 @@ function calc_mortar_flux!(surface_flux_values,
521526
surface_integral, dg::DG, cache)
522527
@unpack neighbor_ids, node_indices = cache.mortars
523528
@unpack contravariant_vectors = cache.elements
524-
@unpack fstar_threaded, fstar_tmp_threaded = cache
529+
@unpack fstar_primary_threaded, fstar_secondary_threaded, fstar_tmp_threaded = cache
525530
index_range = eachnode(dg)
526531

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

532538
# Get index information on the small elements
@@ -555,7 +561,8 @@ function calc_mortar_flux!(surface_flux_values,
555561
i_small, j_small, k_small,
556562
element)
557563

558-
calc_mortar_flux!(fstar, mesh, nonconservative_terms, equations,
564+
calc_mortar_flux!(fstar_primary, fstar_secondary, mesh,
565+
nonconservative_terms, equations,
559566
surface_integral, dg, cache,
560567
mortar, position, normal_direction,
561568
i, j)
@@ -581,14 +588,15 @@ function calc_mortar_flux!(surface_flux_values,
581588
# "mortar_fluxes_to_elements!" instead.
582589
mortar_fluxes_to_elements!(surface_flux_values,
583590
mesh, equations, mortar_l2, dg, cache,
584-
mortar, fstar, u_buffer, fstar_tmp)
591+
mortar, fstar_primary, fstar_secondary, u_buffer,
592+
fstar_tmp)
585593
end
586594

587595
return nothing
588596
end
589597

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

605613
# Copy flux to buffer
606-
set_node_vars!(fstar, flux, equations, dg, i_node_index, j_node_index,
614+
set_node_vars!(fstar_primary, flux, equations, dg, i_node_index, j_node_index,
615+
position_index)
616+
set_node_vars!(fstar_secondary, flux, equations, dg, i_node_index, j_node_index,
607617
position_index)
608618
end
609619

610620
# Inlined version of the mortar flux computation on small elements for conservation fluxes
611621
# with nonconservative terms
612-
@inline function calc_mortar_flux!(fstar,
622+
@inline function calc_mortar_flux!(fstar_primary, fstar_secondary,
613623
mesh::Union{P4estMesh{3}, T8codeMesh{3}},
614624
nonconservative_terms::True, equations,
615625
surface_integral, dg::DG, cache,
@@ -627,20 +637,28 @@ end
627637
# Compute nonconservative flux and add it to the flux scaled by a factor of 0.5 based on
628638
# the interpretation of global SBP operators coupled discontinuously via
629639
# central fluxes/SATs
630-
noncons = nonconservative_flux(u_ll, u_rr, normal_direction, equations)
631-
flux_plus_noncons = flux + 0.5f0 * noncons
640+
noncons_primary = nonconservative_flux(u_ll, u_rr, normal_direction, equations)
641+
noncons_secondary = nonconservative_flux(u_rr, u_ll, normal_direction, equations)
642+
flux_plus_noncons_primary = flux + 0.5f0 * noncons_primary
643+
flux_plus_noncons_secondary = flux + 0.5f0 * noncons_secondary
632644

633645
# Copy to buffer
634-
set_node_vars!(fstar, flux_plus_noncons, equations, dg, i_node_index, j_node_index,
646+
set_node_vars!(fstar_primary, flux_plus_noncons_primary, equations, dg,
647+
i_node_index,
648+
j_node_index,
649+
position_index)
650+
set_node_vars!(fstar_secondary, flux_plus_noncons_secondary, equations, dg,
651+
i_node_index,
652+
j_node_index,
635653
position_index)
636654
end
637655

638656
@inline function mortar_fluxes_to_elements!(surface_flux_values,
639657
mesh::Union{P4estMesh{3}, T8codeMesh{3}},
640658
equations,
641659
mortar_l2::LobattoLegendreMortarL2,
642-
dg::DGSEM, cache, mortar, fstar, u_buffer,
643-
fstar_tmp)
660+
dg::DGSEM, cache, mortar, fstar_primary,
661+
fstar_secondary, u_buffer, fstar_tmp)
644662
@unpack neighbor_ids, node_indices = cache.mortars
645663
index_range = eachnode(dg)
646664

@@ -652,28 +670,30 @@ end
652670
element = neighbor_ids[position, mortar]
653671
for j in eachnode(dg), i in eachnode(dg)
654672
for v in eachvariable(equations)
655-
surface_flux_values[v, i, j, small_direction, element] = fstar[v, i, j,
656-
position]
673+
surface_flux_values[v, i, j, small_direction, element] = fstar_primary[v,
674+
i,
675+
j,
676+
position]
657677
end
658678
end
659679
end
660680

661681
# Project small fluxes to large element.
662682
multiply_dimensionwise!(u_buffer,
663683
mortar_l2.reverse_lower, mortar_l2.reverse_lower,
664-
view(fstar, .., 1),
684+
view(fstar_secondary, .., 1),
665685
fstar_tmp)
666686
add_multiply_dimensionwise!(u_buffer,
667687
mortar_l2.reverse_upper, mortar_l2.reverse_lower,
668-
view(fstar, .., 2),
688+
view(fstar_secondary, .., 2),
669689
fstar_tmp)
670690
add_multiply_dimensionwise!(u_buffer,
671691
mortar_l2.reverse_lower, mortar_l2.reverse_upper,
672-
view(fstar, .., 3),
692+
view(fstar_secondary, .., 3),
673693
fstar_tmp)
674694
add_multiply_dimensionwise!(u_buffer,
675695
mortar_l2.reverse_upper, mortar_l2.reverse_upper,
676-
view(fstar, .., 4),
696+
view(fstar_secondary, .., 4),
677697
fstar_tmp)
678698

679699
# The flux is calculated in the outward direction of the small elements,

‎src/solvers/dgsem_p4est/dg_3d_parabolic.jl

+17-12
Original file line numberDiff line numberDiff line change
@@ -271,7 +271,8 @@ end
271271
mesh::P4estMesh{3},
272272
equations::AbstractEquationsParabolic,
273273
mortar_l2::LobattoLegendreMortarL2,
274-
dg::DGSEM, cache, mortar, fstar, u_buffer,
274+
dg::DGSEM, cache, mortar, fstar_primary,
275+
fstar_secondary, u_buffer,
275276
fstar_tmp)
276277
@unpack neighbor_ids, node_indices = cache.mortars
277278
index_range = eachnode(dg)
@@ -283,28 +284,30 @@ end
283284
element = neighbor_ids[position, mortar]
284285
for j in eachnode(dg), i in eachnode(dg)
285286
for v in eachvariable(equations)
286-
surface_flux_values[v, i, j, small_direction, element] = fstar[v, i, j,
287-
position]
287+
surface_flux_values[v, i, j, small_direction, element] = fstar_primary[v,
288+
i,
289+
j,
290+
position]
288291
end
289292
end
290293
end
291294

292295
# Project small fluxes to large element.
293296
multiply_dimensionwise!(u_buffer,
294297
mortar_l2.reverse_lower, mortar_l2.reverse_lower,
295-
view(fstar, .., 1),
298+
view(fstar_secondary, .., 1),
296299
fstar_tmp)
297300
add_multiply_dimensionwise!(u_buffer,
298301
mortar_l2.reverse_upper, mortar_l2.reverse_lower,
299-
view(fstar, .., 2),
302+
view(fstar_secondary, .., 2),
300303
fstar_tmp)
301304
add_multiply_dimensionwise!(u_buffer,
302305
mortar_l2.reverse_lower, mortar_l2.reverse_upper,
303-
view(fstar, .., 3),
306+
view(fstar_secondary, .., 3),
304307
fstar_tmp)
305308
add_multiply_dimensionwise!(u_buffer,
306309
mortar_l2.reverse_upper, mortar_l2.reverse_upper,
307-
view(fstar, .., 4),
310+
view(fstar_secondary, .., 4),
308311
fstar_tmp)
309312

310313
# The flux is calculated in the outward direction of the small elements,
@@ -788,12 +791,12 @@ function calc_mortar_flux_divergence!(surface_flux_values,
788791
surface_integral, dg::DG, cache)
789792
@unpack neighbor_ids, node_indices = cache.mortars
790793
@unpack contravariant_vectors = cache.elements
791-
@unpack fstar_threaded, fstar_tmp_threaded = cache
794+
@unpack fstar_primary_threaded, fstar_tmp_threaded = cache
792795
index_range = eachnode(dg)
793796

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

799802
# Get index information on the small elements
@@ -842,7 +845,7 @@ function calc_mortar_flux_divergence!(surface_flux_values,
842845
# this reuses the hyperbolic version of `mortar_fluxes_to_elements!`
843846
mortar_fluxes_to_elements!(surface_flux_values,
844847
mesh, equations, mortar_l2, dg, cache,
845-
mortar, fstar, u_buffer, fstar_tmp)
848+
mortar, fstar, fstar, u_buffer, fstar_tmp)
846849
end
847850

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

‎src/solvers/dgsem_p4est/dg_3d_parallel.jl

+2-2
Original file line numberDiff line numberDiff line change
@@ -384,12 +384,12 @@ function calc_mpi_mortar_flux!(surface_flux_values,
384384
surface_integral, dg::DG, cache)
385385
@unpack local_neighbor_ids, local_neighbor_positions, node_indices = cache.mpi_mortars
386386
@unpack contravariant_vectors = cache.elements
387-
@unpack fstar_threaded, fstar_tmp_threaded = cache
387+
@unpack fstar_primary_threaded, fstar_tmp_threaded = cache
388388
index_range = eachnode(dg)
389389

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

395395
# Get index information on the small elements

0 commit comments

Comments
 (0)