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

p4est surface performance 3D #783

Merged
merged 8 commits into from
Aug 13, 2021
Merged

p4est surface performance 3D #783

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
403b3b0
P4estmesh: 3D prolong2interfaces
ranocha Aug 13, 2021
62b2a54
P4estmesh: 3D calc_interface_flux
ranocha Aug 13, 2021
585e7d8
P4estmesh: 3D prolong2boundaries
ranocha Aug 13, 2021
be3c9f8
P4estmesh: 3D calc_boundary_flux
ranocha Aug 13, 2021
bdbe42d
P4estmesh: 3D prolong2mortars
ranocha Aug 13, 2021
5b427d7
P4estmesh: 3D calc_mortar_flux
ranocha Aug 13, 2021
77a2c39
document index_to_start_step_2d/index_to_start_step_3d
ranocha Aug 13, 2021
1b6b974
improve documentation of new auxiliary functions
ranocha Aug 13, 2021
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
3 changes: 2 additions & 1 deletion examples/p4est_3d_dgsem/elixir_advection_cubed_sphere.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,8 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
# ODE solvers, callbacks etc.

# Create ODE problem with time span from 0.0 to 1.0
ode = semidiscretize(semi, (0.0, 1.0));
tspan = (0.0, 1.0)
ode = semidiscretize(semi, tspan)

# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
# and resets the timers
Expand Down
3 changes: 2 additions & 1 deletion examples/p4est_3d_dgsem/elixir_advection_nonconforming.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,8 @@ semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergen
# ODE solvers, callbacks etc.

# Create ODE problem with time span from 0.0 to 1.0
ode = semidiscretize(semi, (0.0, 1.0));
tspan = (0.0, 1.0)
ode = semidiscretize(semi, tspan)

# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
# and resets the timers
Expand Down
4 changes: 2 additions & 2 deletions examples/tree_3d_dgsem/elixir_advection_extended.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,8 +21,8 @@ boundary_conditions = boundary_condition_periodic
# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)

coordinates_min = (-1, -1, -1) # minimum coordinates (min(x), min(y), min(z))
coordinates_max = ( 1, 1, 1) # maximum coordinates (max(x), max(y), max(z))
coordinates_min = (-1.0, -1.0, -1.0) # minimum coordinates (min(x), min(y), min(z))
coordinates_max = ( 1.0, 1.0, 1.0) # maximum coordinates (max(x), max(y), max(z))

# Create a uniformly refined mesh with periodic boundaries
mesh = TreeMesh(coordinates_min, coordinates_max,
Expand Down
17 changes: 0 additions & 17 deletions src/solvers/dgsem_p4est/dg.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,23 +27,6 @@ function create_cache(mesh::P4estMesh, equations::AbstractEquations, dg::DG, ::A
end


# TODO: p4est interface performance, remove
# Extract outward-pointing normal vector (contravariant vector ±Ja^i, i = index) as SVector
# Note that this vector is not normalized
@inline function get_normal_vector(direction, cache, indices...)
@unpack contravariant_vectors = cache.elements

orientation = div(direction + 1, 2)
normal = get_contravariant_vector(orientation, contravariant_vectors, indices...)

# Contravariant vectors at interfaces in negative coordinate direction are pointing inwards
if direction in (1, 3, 5)
normal *= -1
end

return normal
end

# Extract outward-pointing normal direction
# (contravariant vector ±Ja^i, i = index)
# Note that this vector is not normalized
Expand Down
97 changes: 57 additions & 40 deletions src/solvers/dgsem_p4est/dg_2d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,10 +20,31 @@ function create_cache(mesh::P4estMesh{2}, equations, mortar_l2::LobattoLegendreM
end


# TODO: p4est interface performance, move and generalize this function for 3D
@inline function index_to_start_step(index::Symbol, index_range)
# index_to_start_step_2d(index::Symbol, index_range)
#
# Given a symbolic `index` and an `indexrange` (usually `eachnode(dg)`),
# return `index_start, index_step`, i.e., a tuple containing
# - `index_start`, an index value to begin a loop
# - `index_step`, an index step to update during a loop
# The resulting indices translate surface indices to volume indices.
#
# !!! warning
# This assumes that loops using the return values are written as
#
# i_volume_start, i_volume_step = index_to_start_step_2d(symbolic_index_i, index_range)
# j_volume_start, j_volume_step = index_to_start_step_2d(symbolic_index_j, index_range)
#
# i_volume, j_volume = i_volume_start, j_volume_start
# for i_surface in index_range
# # do stuff with `i_surface` and `(i_volume, j_volume)`
#
# i_volume += i_volume_step
# j_volume += j_volume_step
# end
@inline function index_to_start_step_2d(index::Symbol, index_range)
index_begin = first(index_range)
index_end = last(index_range)

if index === :one
return index_begin, 0
elseif index === :end
Expand All @@ -42,16 +63,16 @@ function prolong2interfaces!(cache, u,
index_range = eachnode(dg)

@threaded for interface in eachinterface(dg, cache)
# Copy solution data from the primary element on a case-by-case basis to get
# the correct face and orientation.
# Copy solution data from the primary element on a case-by-case basis
# to get the correct face and orientation.
# Note that in the current implementation, the interface will be
# "aligned at the primary element", i.e., the index of the primary side
# will always run forwards.
primary_element = interfaces.element_ids[1, interface]
primary_indices = interfaces.node_indices[1, interface]

i_primary_start, i_primary_step = index_to_start_step(primary_indices[1], index_range)
j_primary_start, j_primary_step = index_to_start_step(primary_indices[2], index_range)
i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1], index_range)
j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2], index_range)

i_primary = i_primary_start
j_primary = j_primary_start
Expand All @@ -63,13 +84,13 @@ function prolong2interfaces!(cache, u,
j_primary += j_primary_step
end

# Copy solution data from the secondary element on a case-by-case basis to get
# the correct face and orientation.
# Copy solution data from the secondary element on a case-by-case basis
# to get the correct face and orientation.
secondary_element = interfaces.element_ids[2, interface]
secondary_indices = interfaces.node_indices[2, interface]

i_secondary_start, i_secondary_step = index_to_start_step(secondary_indices[1], index_range)
j_secondary_start, j_secondary_step = index_to_start_step(secondary_indices[2], index_range)
i_secondary_start, i_secondary_step = index_to_start_step_2d(secondary_indices[1], index_range)
j_secondary_start, j_secondary_step = index_to_start_step_2d(secondary_indices[2], index_range)

i_secondary = i_secondary_start
j_secondary = j_secondary_start
Expand Down Expand Up @@ -102,8 +123,8 @@ function calc_interface_flux!(surface_flux_values,
primary_indices = node_indices[1, interface]
primary_direction = indices2direction(primary_indices)

i_primary_start, i_primary_step = index_to_start_step(primary_indices[1], index_range)
j_primary_start, j_primary_step = index_to_start_step(primary_indices[2], index_range)
i_primary_start, i_primary_step = index_to_start_step_2d(primary_indices[1], index_range)
j_primary_start, j_primary_step = index_to_start_step_2d(primary_indices[2], index_range)

i_primary = i_primary_start
j_primary = j_primary_start
Expand Down Expand Up @@ -132,8 +153,6 @@ function calc_interface_flux!(surface_flux_values,

# Note that the index of the primary side will always run forward but
# the secondary index might need to run backwards for flipped sides.
# TODO: p4est interface performance; see whether this can be made simpler and
# more general when working on the 3D version
if :i_backwards in secondary_indices
for i in eachnode(dg)
for v in eachvariable(equations)
Expand Down Expand Up @@ -167,8 +186,8 @@ function prolong2boundaries!(cache, u,
element = boundaries.element_ids[boundary]
node_indices = boundaries.node_indices[boundary]

i_node_start, i_node_step = index_to_start_step(node_indices[1], index_range)
j_node_start, j_node_step = index_to_start_step(node_indices[2], index_range)
i_node_start, i_node_step = index_to_start_step_2d(node_indices[1], index_range)
j_node_start, j_node_step = index_to_start_step_2d(node_indices[2], index_range)

i_node = i_node_start
j_node = j_node_start
Expand Down Expand Up @@ -203,8 +222,8 @@ function calc_boundary_flux!(cache, t, boundary_condition, boundary_indexing,
node_indices = boundaries.node_indices[boundary]
direction = indices2direction(node_indices)

i_node_start, i_node_step = index_to_start_step(node_indices[1], index_range)
j_node_start, j_node_step = index_to_start_step(node_indices[2], index_range)
i_node_start, i_node_step = index_to_start_step_2d(node_indices[1], index_range)
j_node_start, j_node_step = index_to_start_step_2d(node_indices[2], index_range)

i_node = i_node_start
j_node = j_node_start
Expand Down Expand Up @@ -240,12 +259,12 @@ function prolong2mortars!(cache, u,
index_range = eachnode(dg)

@threaded for mortar in eachmortar(dg, cache)
# Copy solution data from the small elements on a case-by-case basis to get
# the correct face and orientation.
# Copy solution data from the small elements on a case-by-case basis
# to get the correct face and orientation.
small_indices = node_indices[1, mortar]

i_small_start, i_small_step = index_to_start_step(small_indices[1], index_range)
j_small_start, j_small_step = index_to_start_step(small_indices[2], index_range)
i_small_start, i_small_step = index_to_start_step_2d(small_indices[1], index_range)
j_small_start, j_small_step = index_to_start_step_2d(small_indices[2], index_range)

for position in 1:2
i_small = i_small_start
Expand All @@ -265,11 +284,12 @@ function prolong2mortars!(cache, u,
# before interpolating
u_buffer = cache.u_threaded[Threads.threadid()]

# Copy solution of large element face to buffer in the correct orientation
# Copy solution of large element face to buffer in the
# correct orientation
large_indices = node_indices[2, mortar]

i_large_start, i_large_step = index_to_start_step(large_indices[1], index_range)
j_large_start, j_large_step = index_to_start_step(large_indices[2], index_range)
i_large_start, i_large_step = index_to_start_step_2d(large_indices[1], index_range)
j_large_start, j_large_step = index_to_start_step_2d(large_indices[2], index_range)

i_large = i_large_start
j_large = j_large_start
Expand Down Expand Up @@ -316,8 +336,8 @@ function calc_mortar_flux!(surface_flux_values,
small_indices = node_indices[1, mortar]
small_direction = indices2direction(small_indices)

i_small_start, i_small_step = index_to_start_step(small_indices[1], index_range)
j_small_start, j_small_step = index_to_start_step(small_indices[2], index_range)
i_small_start, i_small_step = index_to_start_step_2d(small_indices[1], index_range)
j_small_start, j_small_step = index_to_start_step_2d(small_indices[2], index_range)

# Contravariant vectors at interfaces in negative coordinate direction
# are pointing inwards. This is handled by `get_normal_direction`.
Expand All @@ -341,8 +361,8 @@ function calc_mortar_flux!(surface_flux_values,
end
end

# Buffer to interpolate flux values of the large element to before copying
# in the correct orientation
# Buffer to interpolate flux values of the large element to before
# copying in the correct orientation
u_buffer = cache.u_threaded[Threads.threadid()]

mortar_fluxes_to_elements!(surface_flux_values,
Expand All @@ -360,9 +380,6 @@ end
dg::DGSEM, cache, mortar, fstar, u_buffer)
@unpack element_ids, node_indices = cache.mortars

large_indices = node_indices[2, mortar]
large_direction = indices2direction(large_indices)

# Copy solution small to small
small_indices = node_indices[1, mortar]
small_direction = indices2direction(small_indices)
Expand All @@ -376,8 +393,6 @@ end
end
end

large_element = element_ids[3, mortar]

# Project small fluxes to large element.
multiply_dimensionwise!(u_buffer,
mortar_l2.reverse_upper, fstar[2],
Expand All @@ -386,18 +401,20 @@ end
# The flux is calculated in the outward direction of the small elements,
# so the sign must be switched to get the flux in outward direction
# of the large element.
# The contravariant vectors of the large element (and therefore the normal vectors
# of the large element as well) are twice as large as the contravariant vectors
# of the small elements. Therefore, the flux need to be scaled by a factor of 2
# to obtain the flux of the large element.
# The contravariant vectors of the large element (and therefore the normal
# vectors of the large element as well) are twice as large as the
# contravariant vectors of the small elements. Therefore, the flux needs
# to be scaled by a factor of 2 to obtain the flux of the large element.
u_buffer .*= -2

# Copy interpolated flux values from buffer to large element face in the
# correct orientation.
# Note that the index of the small sides will always run forward but
# the index of the large side might need to run backwards for flipped sides.
# TODO: p4est interface performance; see whether this can be made simpler and
# more general when working on the 3D version
large_element = element_ids[3, mortar]
large_indices = node_indices[2, mortar]
large_direction = indices2direction(large_indices)

if :i_backwards in large_indices
for i in eachnode(dg)
for v in eachvariable(equations)
Expand Down
Loading