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.

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
Show file tree
Hide file tree
Changes from 4 commits
Commits
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
27 changes: 0 additions & 27 deletions src/solvers/dgsem_tree/dg.jl
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
50 changes: 21 additions & 29 deletions src/solvers/dgsem_tree/dg_1d.jl
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,33 @@ 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]
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
47 changes: 20 additions & 27 deletions src/solvers/dgsem_tree/dg_2d.jl
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,33 @@ 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]
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
50 changes: 21 additions & 29 deletions src/solvers/dgsem_tree/dg_3d.jl
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,33 @@ 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]
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