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

Thread-Parallelize blended DG-FV #2138

Merged
merged 7 commits into from
Oct 31, 2024
Merged

Thread-Parallelize blended DG-FV #2138

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
0919b88
Prototype SCHG 2D thread-para
DanielDoehring Oct 30, 2024
5fdc4cc
remove not-needed datastructures
DanielDoehring Oct 30, 2024
b5bc90c
optimize
DanielDoehring Oct 30, 2024
c4394ea
1D, 3D
DanielDoehring Oct 30, 2024
f535fbf
dgmulti
DanielDoehring Oct 31, 2024
f6ec17e
Apply suggestions from code review
DanielDoehring Oct 31, 2024
8881a9d
Merge branch 'main' into ThreadParaBlendedDG_FV
ranocha Oct 31, 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
77 changes: 22 additions & 55 deletions src/solvers/dgmulti/shock_capturing.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,9 +2,6 @@
function create_cache(mesh::DGMultiMesh{NDIMS}, equations,
volume_integral::VolumeIntegralShockCapturingHG,
dg::DGMultiFluxDiff{<:GaussSBP}, RealT, uEltype) where {NDIMS}
element_ids_dg = Int[]
element_ids_dgfv = Int[]

# build element to element (element_to_element_connectivity) connectivity for smoothing of
# shock capturing parameters.
face_to_face_connectivity = mesh.md.FToF # num_faces x num_elements matrix
Expand All @@ -31,8 +28,7 @@ function create_cache(mesh::DGMultiMesh{NDIMS}, equations,
sparsity_pattern = sum(map(A -> abs.(A)', sparse_hybridized_SBP_operators)) .>
100 * eps()

return (; element_ids_dg, element_ids_dgfv,
sparse_hybridized_SBP_operators, sparsity_pattern,
return (; sparse_hybridized_SBP_operators, sparsity_pattern,
element_to_element_connectivity)
end

Expand Down Expand Up @@ -151,70 +147,41 @@ function apply_smoothing!(mesh::DGMultiMesh, alpha, alpha_tmp, dg::DGMulti, cach
end
end

# pure_and_blended_element_ids!(element_ids_dg, element_ids_dgfv, alpha, dg, cache)
#
# Given blending factors `alpha` and the solver `dg`, fill
# `element_ids_dg` with the IDs of elements using a pure DG scheme and
# `element_ids_dgfv` with the IDs of elements using a blended DG-FV scheme.
function pure_and_blended_element_ids!(element_ids_dg, element_ids_dgfv, alpha,
mesh::DGMultiMesh, dg::DGMulti)
empty!(element_ids_dg)
empty!(element_ids_dgfv)
# For `Float64`, this gives 1.8189894035458565e-12
# For `Float32`, this gives 1.1920929f-5
RealT = eltype(alpha)
atol = max(100 * eps(RealT), eps(RealT)^convert(RealT, 0.75f0))

for element in eachelement(mesh, dg)
# Clip blending factor for values close to zero (-> pure DG)
dg_only = isapprox(alpha[element], 0, atol = atol)
if dg_only
push!(element_ids_dg, element)
else
push!(element_ids_dgfv, element)
end
end

return nothing
end

function calc_volume_integral!(du, u,
mesh::DGMultiMesh,
have_nonconservative_terms, equations,
volume_integral::VolumeIntegralShockCapturingHG,
dg::DGMultiFluxDiff, cache)
(; element_ids_dg, element_ids_dgfv) = cache
(; volume_flux_dg, volume_flux_fv, indicator) = volume_integral

# Calculate blending factors α: u = u_DG * (1 - α) + u_FV * α
alpha = @trixi_timeit timer() "blending factors" indicator(u, mesh, equations, dg,
cache)

# Determine element ids for DG-only and blended DG-FV volume integral
pure_and_blended_element_ids!(element_ids_dg, element_ids_dgfv, alpha, mesh, dg)

# Loop over pure DG elements
@trixi_timeit timer() "pure DG" @threaded for idx_element in eachindex(element_ids_dg)
element = element_ids_dg[idx_element]
flux_differencing_kernel!(du, u, element, mesh, have_nonconservative_terms,
equations, volume_flux_dg, dg, cache)
end
# For `Float64`, this gives 1.8189894035458565e-12
# For `Float32`, this gives 1.1920929f-5
RealT = eltype(alpha)
atol = max(100 * eps(RealT), eps(RealT)^convert(RealT, 0.75f0))

# Loop over blended DG-FV elements, blend the high and low order RHS contributions
# via `rhs_high * (1 - alpha) + rhs_low * (alpha)`.
@trixi_timeit timer() "blended DG-FV" @threaded for idx_element in eachindex(element_ids_dgfv)
element = element_ids_dgfv[idx_element]
@threaded for element in eachelement(mesh, dg)
alpha_element = alpha[element]
# Clip blending factor for values close to zero (-> pure DG)
dg_only = isapprox(alpha_element, 0, atol = atol)

# Calculate DG volume integral contribution
flux_differencing_kernel!(du, u, element, mesh,
have_nonconservative_terms, equations,
volume_flux_dg, dg, cache, 1 - alpha_element)

# Calculate "FV" low order volume integral contribution
low_order_flux_differencing_kernel!(du, u, element, mesh,
have_nonconservative_terms, equations,
volume_flux_fv, dg, cache, alpha_element)
if dg_only
flux_differencing_kernel!(du, u, element, mesh, have_nonconservative_terms,
equations, volume_flux_dg, dg, cache)
else
# Calculate DG volume integral contribution
flux_differencing_kernel!(du, u, element, mesh,
have_nonconservative_terms, equations,
volume_flux_dg, dg, cache, 1 - alpha_element)

# Calculate "FV" low order volume integral contribution
low_order_flux_differencing_kernel!(du, u, element, mesh,
have_nonconservative_terms, equations,
volume_flux_fv, dg, cache, alpha_element)
end
end

return nothing
Expand Down
27 changes: 0 additions & 27 deletions src/solvers/dgsem_tree/dg.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,33 +15,6 @@ function reset_du!(du, dg, cache)
return du
end

# pure_and_blended_element_ids!(element_ids_dg, element_ids_dgfv, alpha, dg, cache)
#
# Given blending factors `alpha` and the solver `dg`, fill
# `element_ids_dg` with the IDs of elements using a pure DG scheme and
# `element_ids_dgfv` with the IDs of elements using a blended DG-FV scheme.
function pure_and_blended_element_ids!(element_ids_dg, element_ids_dgfv, alpha, dg::DG,
cache)
empty!(element_ids_dg)
empty!(element_ids_dgfv)
# For `Float64`, this gives 1.8189894035458565e-12
# For `Float32`, this gives 1.1920929f-5
RealT = eltype(alpha)
atol = max(100 * eps(RealT), eps(RealT)^convert(RealT, 0.75f0))

for element in eachelement(dg, cache)
# Clip blending factor for values close to zero (-> pure DG)
dg_only = isapprox(alpha[element], 0, atol = atol)
if dg_only
push!(element_ids_dg, element)
else
push!(element_ids_dgfv, element)
end
end

return nothing
end

function volume_jacobian(element, mesh::TreeMesh, cache)
return inv(cache.elements.inverse_jacobian[element])^ndims(mesh)
end
Expand Down
52 changes: 23 additions & 29 deletions src/solvers/dgsem_tree/dg_1d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,9 +42,6 @@ end
function create_cache(mesh::Union{TreeMesh{1}, StructuredMesh{1}, P4estMesh{1}},
equations,
volume_integral::VolumeIntegralShockCapturingHG, dg::DG, uEltype)
element_ids_dg = Int[]
element_ids_dgfv = Int[]

cache = create_cache(mesh, equations,
VolumeIntegralFluxDifferencing(volume_integral.volume_flux_dg),
dg, uEltype)
Expand All @@ -55,8 +52,7 @@ function create_cache(mesh::Union{TreeMesh{1}, StructuredMesh{1}, P4estMesh{1}},
fstar1_R_threaded = A2dp1_x[A2dp1_x(undef, nvariables(equations), nnodes(dg) + 1)
for _ in 1:Threads.nthreads()]

return (; cache..., element_ids_dg, element_ids_dgfv, fstar1_L_threaded,
fstar1_R_threaded)
return (; cache..., fstar1_L_threaded, fstar1_R_threaded)
end

function create_cache(mesh::Union{TreeMesh{1}, StructuredMesh{1}, P4estMesh{1}},
Expand Down Expand Up @@ -256,37 +252,35 @@ function calc_volume_integral!(du, u,
nonconservative_terms, equations,
volume_integral::VolumeIntegralShockCapturingHG,
dg::DGSEM, cache)
@unpack element_ids_dg, element_ids_dgfv = cache
@unpack volume_flux_dg, volume_flux_fv, indicator = volume_integral

# Calculate blending factors α: u = u_DG * (1 - α) + u_FV * α
alpha = @trixi_timeit timer() "blending factors" indicator(u, mesh, equations, dg,
cache)

# Determine element ids for DG-only and blended DG-FV volume integral
pure_and_blended_element_ids!(element_ids_dg, element_ids_dgfv, alpha, dg, cache)

# Loop over pure DG elements
@trixi_timeit timer() "pure DG" @threaded for idx_element in eachindex(element_ids_dg)
element = element_ids_dg[idx_element]
flux_differencing_kernel!(du, u, element, mesh, nonconservative_terms,
equations,
volume_flux_dg, dg, cache)
end

# Loop over blended DG-FV elements
@trixi_timeit timer() "blended DG-FV" @threaded for idx_element in eachindex(element_ids_dgfv)
element = element_ids_dgfv[idx_element]
# For `Float64`, this gives 1.8189894035458565e-12
# For `Float32`, this gives 1.1920929f-5
RealT = eltype(alpha)
atol = max(100 * eps(RealT), eps(RealT)^convert(RealT, 0.75f0))
DanielDoehring marked this conversation as resolved.
Show resolved Hide resolved
@threaded for element in eachelement(dg, cache)
alpha_element = alpha[element]

# Calculate DG volume integral contribution
flux_differencing_kernel!(du, u, element, mesh, nonconservative_terms,
equations,
volume_flux_dg, dg, cache, 1 - alpha_element)

# Calculate FV volume integral contribution
fv_kernel!(du, u, mesh, nonconservative_terms, equations, volume_flux_fv,
dg, cache, element, alpha_element)
# Clip blending factor for values close to zero (-> pure DG)
dg_only = isapprox(alpha_element, 0, atol = atol)

if dg_only
flux_differencing_kernel!(du, u, element, mesh,
nonconservative_terms, equations,
volume_flux_dg, dg, cache)
else
# Calculate DG volume integral contribution
flux_differencing_kernel!(du, u, element, mesh,
nonconservative_terms, equations,
volume_flux_dg, dg, cache, 1 - alpha_element)

# Calculate FV volume integral contribution
fv_kernel!(du, u, mesh, nonconservative_terms, equations, volume_flux_fv,
dg, cache, element, alpha_element)
end
end

return nothing
Expand Down
49 changes: 22 additions & 27 deletions src/solvers/dgsem_tree/dg_2d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -47,9 +47,6 @@ end
function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMesh2D,
P4estMesh{2}, T8codeMesh{2}}, equations,
volume_integral::VolumeIntegralShockCapturingHG, dg::DG, uEltype)
element_ids_dg = Int[]
element_ids_dgfv = Int[]

cache = create_cache(mesh, equations,
VolumeIntegralFluxDifferencing(volume_integral.volume_flux_dg),
dg, uEltype)
Expand All @@ -66,7 +63,7 @@ function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}, UnstructuredMe
fstar2_R_threaded = A3dp1_y[A3dp1_y(undef, nvariables(equations), nnodes(dg),
nnodes(dg) + 1) for _ in 1:Threads.nthreads()]

return (; cache..., element_ids_dg, element_ids_dgfv,
return (; cache...,
fstar1_L_threaded, fstar1_R_threaded, fstar2_L_threaded, fstar2_R_threaded)
end

Expand Down Expand Up @@ -338,37 +335,35 @@ function calc_volume_integral!(du, u,
nonconservative_terms, equations,
volume_integral::VolumeIntegralShockCapturingHG,
dg::DGSEM, cache)
@unpack element_ids_dg, element_ids_dgfv = cache
@unpack volume_flux_dg, volume_flux_fv, indicator = volume_integral

# Calculate blending factors α: u = u_DG * (1 - α) + u_FV * α
alpha = @trixi_timeit timer() "blending factors" indicator(u, mesh, equations, dg,
cache)

# Determine element ids for DG-only and blended DG-FV volume integral
pure_and_blended_element_ids!(element_ids_dg, element_ids_dgfv, alpha, dg, cache)

# Loop over pure DG elements
@trixi_timeit timer() "pure DG" @threaded for idx_element in eachindex(element_ids_dg)
element = element_ids_dg[idx_element]
flux_differencing_kernel!(du, u, element, mesh,
nonconservative_terms, equations,
volume_flux_dg, dg, cache)
end

# Loop over blended DG-FV elements
@trixi_timeit timer() "blended DG-FV" @threaded for idx_element in eachindex(element_ids_dgfv)
element = element_ids_dgfv[idx_element]
# For `Float64`, this gives 1.8189894035458565e-12
# For `Float32`, this gives 1.1920929f-5
RealT = eltype(alpha)
atol = max(100 * eps(RealT), eps(RealT)^convert(RealT, 0.75f0))
DanielDoehring marked this conversation as resolved.
Show resolved Hide resolved
@threaded for element in eachelement(dg, cache)
alpha_element = alpha[element]
# Clip blending factor for values close to zero (-> pure DG)
dg_only = isapprox(alpha_element, 0, atol = atol)

# Calculate DG volume integral contribution
flux_differencing_kernel!(du, u, element, mesh,
nonconservative_terms, equations,
volume_flux_dg, dg, cache, 1 - alpha_element)

# Calculate FV volume integral contribution
fv_kernel!(du, u, mesh, nonconservative_terms, equations, volume_flux_fv,
dg, cache, element, alpha_element)
if dg_only
flux_differencing_kernel!(du, u, element, mesh,
nonconservative_terms, equations,
volume_flux_dg, dg, cache)
else
# Calculate DG volume integral contribution
flux_differencing_kernel!(du, u, element, mesh,
nonconservative_terms, equations,
volume_flux_dg, dg, cache, 1 - alpha_element)

# Calculate FV volume integral contribution
fv_kernel!(du, u, mesh, nonconservative_terms, equations, volume_flux_fv,
dg, cache, element, alpha_element)
end
end

return nothing
Expand Down
52 changes: 23 additions & 29 deletions src/solvers/dgsem_tree/dg_3d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -47,9 +47,6 @@ function create_cache(mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3},
T8codeMesh{3}},
equations,
volume_integral::VolumeIntegralShockCapturingHG, dg::DG, uEltype)
element_ids_dg = Int[]
element_ids_dgfv = Int[]

cache = create_cache(mesh, equations,
VolumeIntegralFluxDifferencing(volume_integral.volume_flux_dg),
dg, uEltype)
Expand All @@ -76,8 +73,7 @@ function create_cache(mesh::Union{TreeMesh{3}, StructuredMesh{3}, P4estMesh{3},
nnodes(dg), nnodes(dg) + 1)
for _ in 1:Threads.nthreads()]

return (; cache..., element_ids_dg, element_ids_dgfv, fstar1_L_threaded,
fstar1_R_threaded,
return (; cache..., fstar1_L_threaded, fstar1_R_threaded,
fstar2_L_threaded, fstar2_R_threaded, fstar3_L_threaded, fstar3_R_threaded)
end

Expand Down Expand Up @@ -386,37 +382,35 @@ function calc_volume_integral!(du, u,
nonconservative_terms, equations,
volume_integral::VolumeIntegralShockCapturingHG,
dg::DGSEM, cache)
@unpack element_ids_dg, element_ids_dgfv = cache
@unpack volume_flux_dg, volume_flux_fv, indicator = volume_integral

# Calculate blending factors α: u = u_DG * (1 - α) + u_FV * α
alpha = @trixi_timeit timer() "blending factors" indicator(u, mesh, equations, dg,
cache)

# Determine element ids for DG-only and blended DG-FV volume integral
pure_and_blended_element_ids!(element_ids_dg, element_ids_dgfv, alpha, dg, cache)

# Loop over pure DG elements
@trixi_timeit timer() "pure DG" @threaded for idx_element in eachindex(element_ids_dg)
element = element_ids_dg[idx_element]
flux_differencing_kernel!(du, u, element, mesh,
nonconservative_terms, equations,
volume_flux_dg, dg, cache)
end

# Loop over blended DG-FV elements
@trixi_timeit timer() "blended DG-FV" @threaded for idx_element in eachindex(element_ids_dgfv)
element = element_ids_dgfv[idx_element]
# For `Float64`, this gives 1.8189894035458565e-12
# For `Float32`, this gives 1.1920929f-5
RealT = eltype(alpha)
atol = max(100 * eps(RealT), eps(RealT)^convert(RealT, 0.75f0))
DanielDoehring marked this conversation as resolved.
Show resolved Hide resolved
@threaded for element in eachelement(dg, cache)
alpha_element = alpha[element]

# Calculate DG volume integral contribution
flux_differencing_kernel!(du, u, element, mesh,
nonconservative_terms, equations,
volume_flux_dg, dg, cache, 1 - alpha_element)

# Calculate FV volume integral contribution
fv_kernel!(du, u, mesh, nonconservative_terms, equations, volume_flux_fv,
dg, cache, element, alpha_element)
# Clip blending factor for values close to zero (-> pure DG)
dg_only = isapprox(alpha_element, 0, atol = atol)

if dg_only
flux_differencing_kernel!(du, u, element, mesh,
nonconservative_terms, equations,
volume_flux_dg, dg, cache)
else
# Calculate DG volume integral contribution
flux_differencing_kernel!(du, u, element, mesh,
nonconservative_terms, equations,
volume_flux_dg, dg, cache, 1 - alpha_element)

# Calculate FV volume integral contribution
fv_kernel!(du, u, mesh, nonconservative_terms, equations, volume_flux_fv,
dg, cache, element, alpha_element)
end
end

return nothing
Expand Down
Loading