Skip to content

Commit

Permalink
Some multi-threading improvements (#1630)
Browse files Browse the repository at this point in the history 
* fix multi-threaded parabolic terms on ARM

On ARM, the previous versions resulted in
  cfunction: closures are not supported on this platform
With this change, everything seems to work fine locally.
At least test/test_threaded.jl runs fine with two threads.

* reduce alloactions of multi-threaded parabolic terms a bit

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.

* format
  • Loading branch information
@ranocha
ranocha authored Sep 12, 2023
1 parent 953f88a commit 7791faa
Show file tree
Hide file tree
Showing 4 changed files with 86 additions and 69 deletions.
5 changes: 3 additions & 2 deletions src/solvers/dgmulti/dg_parabolic.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -62,9 +62,10 @@ end
function transform_variables!(u_transformed, u, mesh,
equations_parabolic::AbstractEquationsParabolic,
dg::DGMulti, parabolic_scheme, cache, cache_parabolic)
transformation = gradient_variable_transformation(equations_parabolic)

@threaded for i in eachindex(u)
u_transformed[i] = gradient_variable_transformation(equations_parabolic)(u[i],
equations_parabolic)
u_transformed[i] = transformation(u[i], equations_parabolic)
end
end

Expand Down
34 changes: 20 additions & 14 deletions src/solvers/dgsem_tree/dg_1d_parabolic.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -105,12 +105,13 @@ end
function transform_variables!(u_transformed, u, mesh::TreeMesh{1},
equations_parabolic::AbstractEquationsParabolic,
dg::DG, parabolic_scheme, cache, cache_parabolic)
transformation = gradient_variable_transformation(equations_parabolic)

@threaded for element in eachelement(dg, cache)
# Calculate volume terms in one element
for i in eachnode(dg)
u_node = get_node_vars(u, equations_parabolic, dg, i, element)
u_transformed_node = gradient_variable_transformation(equations_parabolic)(u_node,
equations_parabolic)
u_transformed_node = transformation(u_node, equations_parabolic)
set_node_vars!(u_transformed, u_transformed_node, equations_parabolic, dg,
i, element)
end
Expand Down Expand Up @@ -147,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 @@ -204,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 @@ -552,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
51 changes: 28 additions & 23 deletions src/solvers/dgsem_tree/dg_2d_parabolic.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -118,12 +118,13 @@ end
function transform_variables!(u_transformed, u, mesh::Union{TreeMesh{2}, P4estMesh{2}},
equations_parabolic::AbstractEquationsParabolic,
dg::DG, parabolic_scheme, cache, cache_parabolic)
transformation = gradient_variable_transformation(equations_parabolic)

@threaded for element in eachelement(dg, cache)
# Calculate volume terms in one element
for j in eachnode(dg), i in eachnode(dg)
u_node = get_node_vars(u, equations_parabolic, dg, i, j, element)
u_transformed_node = gradient_variable_transformation(equations_parabolic)(u_node,
equations_parabolic)
u_transformed_node = transformation(u_node, equations_parabolic)
set_node_vars!(u_transformed, u_transformed_node, equations_parabolic, dg,
i, j, element)
end
Expand Down Expand Up @@ -168,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 @@ -244,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 @@ -608,7 +611,7 @@ function prolong2mortars!(cache, flux_viscous::Tuple{AbstractArray, AbstractArra
end

# NOTE: Use analogy to "calc_mortar_flux!" for hyperbolic eqs with no nonconservative terms.
# Reasoning: "calc_interface_flux!" for parabolic part is implemented as the version for
# Reasoning: "calc_interface_flux!" for parabolic part is implemented as the version for
# hyperbolic terms with conserved terms only, i.e., no nonconservative terms.
function calc_mortar_flux!(surface_flux_values,
mesh::TreeMesh{2},
Expand Down Expand Up @@ -934,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
Loading

0 comments on commit 7791faa

Please sign in to comment.