Skip to content

Commit

Permalink
reduce alloactions of multi-threaded parabolic terms a bit
Browse files Browse the repository at this point in the history 
Polyester.jl passes arrays as pointer arrays to the closures without requiring allocations.
More complicated structs may still require allocations, so unpacking some arrays before entering a threaded loop can reduce allocations.
  • Loading branch information
@ranocha
ranocha committed Sep 12, 2023
1 parent dffec18 commit 935b4c2
Show file tree
Hide file tree
Showing 3 changed files with 72 additions and 59 deletions.
29 changes: 17 additions & 12 deletions src/solvers/dgsem_tree/dg_1d_parabolic.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -148,16 +148,18 @@ function prolong2interfaces!(cache_parabolic, flux_viscous,
equations_parabolic::AbstractEquationsParabolic,
surface_integral, dg::DG, cache)
@unpack interfaces = cache_parabolic
@unpack neighbor_ids = interfaces
interfaces_u = interfaces.u

@threaded for interface in eachinterface(dg, cache)
left_element = interfaces.neighbor_ids[1, interface]
right_element = interfaces.neighbor_ids[2, interface]
left_element = neighbor_ids[1, interface]
right_element = neighbor_ids[2, interface]

# interface in x-direction
for v in eachvariable(equations_parabolic)
# OBS! `interfaces.u` stores the interpolated *fluxes* and *not the solution*!
interfaces.u[1, v, interface] = flux_viscous[v, nnodes(dg), left_element]
interfaces.u[2, v, interface] = flux_viscous[v, 1, right_element]
# OBS! `interfaces_u` stores the interpolated *fluxes* and *not the solution*!
interfaces_u[1, v, interface] = flux_viscous[v, nnodes(dg), left_element]
interfaces_u[2, v, interface] = flux_viscous[v, 1, right_element]
end
end

Expand Down Expand Up @@ -205,21 +207,22 @@ function prolong2boundaries!(cache_parabolic, flux_viscous,
equations_parabolic::AbstractEquationsParabolic,
surface_integral, dg::DG, cache)
@unpack boundaries = cache_parabolic
@unpack neighbor_sides = boundaries
@unpack neighbor_sides, neighbor_ids = boundaries
boundaries_u = boundaries.u

@threaded for boundary in eachboundary(dg, cache_parabolic)
element = boundaries.neighbor_ids[boundary]
element = neighbor_ids[boundary]

if neighbor_sides[boundary] == 1
# element in -x direction of boundary
for v in eachvariable(equations_parabolic)
# OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*!
boundaries.u[1, v, boundary] = flux_viscous[v, nnodes(dg), element]
# OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*!
boundaries_u[1, v, boundary] = flux_viscous[v, nnodes(dg), element]
end
else # Element in +x direction of boundary
for v in eachvariable(equations_parabolic)
# OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*!
boundaries.u[2, v, boundary] = flux_viscous[v, 1, element]
# OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*!
boundaries_u[2, v, boundary] = flux_viscous[v, 1, element]
end
end
end
Expand Down Expand Up @@ -553,8 +556,10 @@ end
# where f(u) is the inviscid flux and g(u) is the viscous flux.
function apply_jacobian_parabolic!(du, mesh::TreeMesh{1},
equations::AbstractEquationsParabolic, dg::DG, cache)
@unpack inverse_jacobian = cache.elements

@threaded for element in eachelement(dg, cache)
factor = cache.elements.inverse_jacobian[element]
factor = inverse_jacobian[element]

for i in eachnode(dg)
for v in eachvariable(equations)
Expand Down
44 changes: 24 additions & 20 deletions src/solvers/dgsem_tree/dg_2d_parabolic.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -169,30 +169,31 @@ function prolong2interfaces!(cache_parabolic, flux_viscous,
equations_parabolic::AbstractEquationsParabolic,
surface_integral, dg::DG, cache)
@unpack interfaces = cache_parabolic
@unpack orientations = interfaces
@unpack orientations, neighbor_ids = interfaces
interfaces_u = interfaces.u

flux_viscous_x, flux_viscous_y = flux_viscous

@threaded for interface in eachinterface(dg, cache)
left_element = interfaces.neighbor_ids[1, interface]
right_element = interfaces.neighbor_ids[2, interface]
left_element = neighbor_ids[1, interface]
right_element = neighbor_ids[2, interface]

if orientations[interface] == 1
# interface in x-direction
for j in eachnode(dg), v in eachvariable(equations_parabolic)
# OBS! `interfaces.u` stores the interpolated *fluxes* and *not the solution*!
interfaces.u[1, v, j, interface] = flux_viscous_x[v, nnodes(dg), j,
# OBS! `interfaces_u` stores the interpolated *fluxes* and *not the solution*!
interfaces_u[1, v, j, interface] = flux_viscous_x[v, nnodes(dg), j,
left_element]
interfaces.u[2, v, j, interface] = flux_viscous_x[v, 1, j,
interfaces_u[2, v, j, interface] = flux_viscous_x[v, 1, j,
right_element]
end
else # if orientations[interface] == 2
# interface in y-direction
for i in eachnode(dg), v in eachvariable(equations_parabolic)
# OBS! `interfaces.u` stores the interpolated *fluxes* and *not the solution*!
interfaces.u[1, v, i, interface] = flux_viscous_y[v, i, nnodes(dg),
# OBS! `interfaces_u` stores the interpolated *fluxes* and *not the solution*!
interfaces_u[1, v, i, interface] = flux_viscous_y[v, i, nnodes(dg),
left_element]
interfaces.u[2, v, i, interface] = flux_viscous_y[v, i, 1,
interfaces_u[2, v, i, interface] = flux_viscous_y[v, i, 1,
right_element]
end
end
Expand Down Expand Up @@ -245,41 +246,42 @@ function prolong2boundaries!(cache_parabolic, flux_viscous,
equations_parabolic::AbstractEquationsParabolic,
surface_integral, dg::DG, cache)
@unpack boundaries = cache_parabolic
@unpack orientations, neighbor_sides = boundaries
@unpack orientations, neighbor_sides, neighbor_ids = boundaries
boundaries_u = boundaries.u
flux_viscous_x, flux_viscous_y = flux_viscous

@threaded for boundary in eachboundary(dg, cache_parabolic)
element = boundaries.neighbor_ids[boundary]
element = neighbor_ids[boundary]

if orientations[boundary] == 1
# boundary in x-direction
if neighbor_sides[boundary] == 1
# element in -x direction of boundary
for l in eachnode(dg), v in eachvariable(equations_parabolic)
# OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*!
boundaries.u[1, v, l, boundary] = flux_viscous_x[v, nnodes(dg), l,
# OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*!
boundaries_u[1, v, l, boundary] = flux_viscous_x[v, nnodes(dg), l,
element]
end
else # Element in +x direction of boundary
for l in eachnode(dg), v in eachvariable(equations_parabolic)
# OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*!
boundaries.u[2, v, l, boundary] = flux_viscous_x[v, 1, l, element]
# OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*!
boundaries_u[2, v, l, boundary] = flux_viscous_x[v, 1, l, element]
end
end
else # if orientations[boundary] == 2
# boundary in y-direction
if neighbor_sides[boundary] == 1
# element in -y direction of boundary
for l in eachnode(dg), v in eachvariable(equations_parabolic)
# OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*!
boundaries.u[1, v, l, boundary] = flux_viscous_y[v, l, nnodes(dg),
# OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*!
boundaries_u[1, v, l, boundary] = flux_viscous_y[v, l, nnodes(dg),
element]
end
else
# element in +y direction of boundary
for l in eachnode(dg), v in eachvariable(equations_parabolic)
# OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*!
boundaries.u[2, v, l, boundary] = flux_viscous_y[v, l, 1, element]
# OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*!
boundaries_u[2, v, l, boundary] = flux_viscous_y[v, l, 1, element]
end
end
end
Expand Down Expand Up @@ -935,8 +937,10 @@ end
# where f(u) is the inviscid flux and g(u) is the viscous flux.
function apply_jacobian_parabolic!(du, mesh::Union{TreeMesh{2}, P4estMesh{2}},
equations::AbstractEquationsParabolic, dg::DG, cache)
@unpack inverse_jacobian = cache.elements

@threaded for element in eachelement(dg, cache)
factor = cache.elements.inverse_jacobian[element]
factor = inverse_jacobian[element]

for j in eachnode(dg), i in eachnode(dg)
for v in eachvariable(equations)
Expand Down
58 changes: 31 additions & 27 deletions src/solvers/dgsem_tree/dg_3d_parabolic.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -176,43 +176,44 @@ function prolong2interfaces!(cache_parabolic, flux_viscous,
equations_parabolic::AbstractEquationsParabolic,
surface_integral, dg::DG, cache)
@unpack interfaces = cache_parabolic
@unpack orientations = interfaces
@unpack orientations, neighbor_ids = interfaces
interfaces_u = interfaces.u

flux_viscous_x, flux_viscous_y, flux_viscous_z = flux_viscous

@threaded for interface in eachinterface(dg, cache)
left_element = interfaces.neighbor_ids[1, interface]
right_element = interfaces.neighbor_ids[2, interface]
left_element = neighbor_ids[1, interface]
right_element = neighbor_ids[2, interface]

if orientations[interface] == 1
# interface in x-direction
for k in eachnode(dg), j in eachnode(dg),
v in eachvariable(equations_parabolic)
# OBS! `interfaces.u` stores the interpolated *fluxes* and *not the solution*!
interfaces.u[1, v, j, k, interface] = flux_viscous_x[v, nnodes(dg), j,
# OBS! `interfaces_u` stores the interpolated *fluxes* and *not the solution*!
interfaces_u[1, v, j, k, interface] = flux_viscous_x[v, nnodes(dg), j,
k, left_element]
interfaces.u[2, v, j, k, interface] = flux_viscous_x[v, 1, j, k,
interfaces_u[2, v, j, k, interface] = flux_viscous_x[v, 1, j, k,
right_element]
end
elseif orientations[interface] == 2
# interface in y-direction
for k in eachnode(dg), i in eachnode(dg),
v in eachvariable(equations_parabolic)
# OBS! `interfaces.u` stores the interpolated *fluxes* and *not the solution*!
interfaces.u[1, v, i, k, interface] = flux_viscous_y[v, i, nnodes(dg),
# OBS! `interfaces_u` stores the interpolated *fluxes* and *not the solution*!
interfaces_u[1, v, i, k, interface] = flux_viscous_y[v, i, nnodes(dg),
k, left_element]
interfaces.u[2, v, i, k, interface] = flux_viscous_y[v, i, 1, k,
interfaces_u[2, v, i, k, interface] = flux_viscous_y[v, i, 1, k,
right_element]
end
else # if orientations[interface] == 3
# interface in z-direction
for j in eachnode(dg), i in eachnode(dg),
v in eachvariable(equations_parabolic)
# OBS! `interfaces.u` stores the interpolated *fluxes* and *not the solution*!
interfaces.u[1, v, i, j, interface] = flux_viscous_z[v, i, j,
# OBS! `interfaces_u` stores the interpolated *fluxes* and *not the solution*!
interfaces_u[1, v, i, j, interface] = flux_viscous_z[v, i, j,
nnodes(dg),
left_element]
interfaces.u[2, v, i, j, interface] = flux_viscous_z[v, i, j, 1,
interfaces_u[2, v, i, j, interface] = flux_viscous_z[v, i, j, 1,
right_element]
end
end
Expand Down Expand Up @@ -266,27 +267,28 @@ function prolong2boundaries!(cache_parabolic, flux_viscous,
equations_parabolic::AbstractEquationsParabolic,
surface_integral, dg::DG, cache)
@unpack boundaries = cache_parabolic
@unpack orientations, neighbor_sides = boundaries
@unpack orientations, neighbor_sides, neighbor_ids = boundaries
boundaries_u = boundaries.u
flux_viscous_x, flux_viscous_y, flux_viscous_z = flux_viscous

@threaded for boundary in eachboundary(dg, cache_parabolic)
element = boundaries.neighbor_ids[boundary]
element = neighbor_ids[boundary]

if orientations[boundary] == 1
# boundary in x-direction
if neighbor_sides[boundary] == 1
# element in -x direction of boundary
for k in eachnode(dg), j in eachnode(dg),
v in eachvariable(equations_parabolic)
# OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*!
boundaries.u[1, v, j, k, boundary] = flux_viscous_x[v, nnodes(dg),
# OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*!
boundaries_u[1, v, j, k, boundary] = flux_viscous_x[v, nnodes(dg),
j, k, element]
end
else # Element in +x direction of boundary
for k in eachnode(dg), j in eachnode(dg),
v in eachvariable(equations_parabolic)
# OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*!
boundaries.u[2, v, j, k, boundary] = flux_viscous_x[v, 1, j, k,
# OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*!
boundaries_u[2, v, j, k, boundary] = flux_viscous_x[v, 1, j, k,
element]
end
end
Expand All @@ -296,17 +298,17 @@ function prolong2boundaries!(cache_parabolic, flux_viscous,
# element in -y direction of boundary
for k in eachnode(dg), i in eachnode(dg),
v in eachvariable(equations_parabolic)
# OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*!
boundaries.u[1, v, i, k, boundary] = flux_viscous_y[v, i,
# OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*!
boundaries_u[1, v, i, k, boundary] = flux_viscous_y[v, i,
nnodes(dg), k,
element]
end
else
# element in +y direction of boundary
for k in eachnode(dg), i in eachnode(dg),
v in eachvariable(equations_parabolic)
# OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*!
boundaries.u[2, v, i, k, boundary] = flux_viscous_y[v, i, 1, k,
# OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*!
boundaries_u[2, v, i, k, boundary] = flux_viscous_y[v, i, 1, k,
element]
end
end
Expand All @@ -316,17 +318,17 @@ function prolong2boundaries!(cache_parabolic, flux_viscous,
# element in -z direction of boundary
for j in eachnode(dg), i in eachnode(dg),
v in eachvariable(equations_parabolic)
# OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*!
boundaries.u[1, v, i, j, boundary] = flux_viscous_z[v, i, j,
# OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*!
boundaries_u[1, v, i, j, boundary] = flux_viscous_z[v, i, j,
nnodes(dg),
element]
end
else
# element in +z direction of boundary
for j in eachnode(dg), i in eachnode(dg),
v in eachvariable(equations_parabolic)
# OBS! `boundaries.u` stores the interpolated *fluxes* and *not the solution*!
boundaries.u[2, v, i, j, boundary] = flux_viscous_z[v, i, j, 1,
# OBS! `boundaries_u` stores the interpolated *fluxes* and *not the solution*!
boundaries_u[2, v, i, j, boundary] = flux_viscous_z[v, i, j, 1,
element]
end
end
Expand Down Expand Up @@ -1125,8 +1127,10 @@ end
# where f(u) is the inviscid flux and g(u) is the viscous flux.
function apply_jacobian_parabolic!(du, mesh::Union{TreeMesh{3}, P4estMesh{3}},
equations::AbstractEquationsParabolic, dg::DG, cache)
@unpack inverse_jacobian = cache.elements

@threaded for element in eachelement(dg, cache)
factor = cache.elements.inverse_jacobian[element]
factor = inverse_jacobian[element]

for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
for v in eachvariable(equations)
Expand Down

0 comments on commit 935b4c2

Please sign in to comment.