Skip to content

Commit

Permalink
1D, 3D
Browse files Browse the repository at this point in the history
  • Loading branch information
@DanielDoehring
DanielDoehring committed Oct 30, 2024
1 parent b5bc90c commit c4394ea
Show file tree
Hide file tree
Showing 3 changed files with 42 additions and 85 deletions.
27 changes: 0 additions & 27 deletions src/solvers/dgsem_tree/dg.jl
View file
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
View file
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))
@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
50 changes: 21 additions & 29 deletions src/solvers/dgsem_tree/dg_3d.jl
View file
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))
@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

0 comments on commit c4394ea

Please sign in to comment.