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

Remove interfaces and boundaries from cache_parabolic #1754

Draft
wants to merge 10 commits into
base: main
Choose a base branch
from
8 changes: 4 additions & 4 deletions src/callbacks_step/amr_dg1d.jl
Original file line number Diff line number Diff line change
Expand Up @@ -89,9 +89,9 @@ function refine!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1},
reinitialize_containers!(mesh, equations, dg, cache_parabolic)

# Sanity check
@unpack interfaces = cache_parabolic
(; interfaces) = cache
if isperiodic(mesh.tree)
@assert ninterfaces(interfaces)==1 * nelements(dg, cache_parabolic) ("For 1D and periodic domains, the number of interfaces must be the same as the number of elements")
@assert ninterfaces(interfaces)==1 * nelements(dg, cache) ("For 1D and periodic domains, the number of interfaces must be the same as the number of elements")
end

return nothing
Expand Down Expand Up @@ -235,9 +235,9 @@ function coarsen!(u_ode::AbstractVector, adaptor, mesh::TreeMesh{1},
reinitialize_containers!(mesh, equations, dg, cache_parabolic)

# Sanity check
@unpack interfaces = cache_parabolic
@unpack interfaces = cache
if isperiodic(mesh.tree)
@assert ninterfaces(interfaces)==1 * nelements(dg, cache_parabolic) ("For 1D and periodic domains, the number of interfaces must be the same as the number of elements")
@assert ninterfaces(interfaces)==1 * nelements(dg, cache) ("For 1D and periodic domains, the number of interfaces must be the same as the number of elements")
end

return nothing
Expand Down
51 changes: 25 additions & 26 deletions src/solvers/dgsem_tree/dg_1d_parabolic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -62,39 +62,40 @@ function rhs_parabolic!(du, u, t, mesh::TreeMesh{1},

# Prolong solution to interfaces
@trixi_timeit timer() "prolong2interfaces" begin
prolong2interfaces!(cache_parabolic, flux_viscous, mesh, equations_parabolic,
prolong2interfaces!(cache.interfaces,
flux_viscous, mesh, equations_parabolic,
dg.surface_integral, dg, cache)
end

# Calculate interface fluxes
@trixi_timeit timer() "interface flux" begin
calc_interface_flux!(cache_parabolic.elements.surface_flux_values, mesh,
equations_parabolic, dg, cache_parabolic)
calc_interface_flux!(cache.elements.surface_flux_values, mesh,
equations_parabolic, dg, cache)
end

# Prolong solution to boundaries
@trixi_timeit timer() "prolong2boundaries" begin
prolong2boundaries!(cache_parabolic, flux_viscous, mesh, equations_parabolic,
prolong2boundaries!(cache.boundaries, flux_viscous, mesh, equations_parabolic,
dg.surface_integral, dg, cache)
end

# Calculate boundary fluxes
@trixi_timeit timer() "boundary flux" begin
calc_boundary_flux_divergence!(cache_parabolic, t,
calc_boundary_flux_divergence!(cache, t,
boundary_conditions_parabolic, mesh,
equations_parabolic,
dg.surface_integral, dg)
end

# Calculate surface integrals
@trixi_timeit timer() "surface integral" begin
calc_surface_integral!(du, u, mesh, equations_parabolic,
dg.surface_integral, dg, cache_parabolic)
calc_surface_integral!(du, u, mesh, equations_parabolic,
dg.surface_integral, dg, cache)
end

# Apply Jacobian from mapping to reference element
@trixi_timeit timer() "Jacobian" begin
apply_jacobian_parabolic!(du, mesh, equations_parabolic, dg, cache_parabolic)
apply_jacobian_parabolic!(du, mesh, equations_parabolic, dg, cache)
end

return nothing
Expand Down Expand Up @@ -144,11 +145,10 @@ end

# This is the version used when calculating the divergence of the viscous fluxes
# We pass the `surface_integral` argument solely for dispatch
function prolong2interfaces!(cache_parabolic, flux_viscous,
function prolong2interfaces!(interfaces, flux_viscous,
mesh::TreeMesh{1},
equations_parabolic::AbstractEquationsParabolic,
surface_integral, dg::DG, cache)
@unpack interfaces = cache_parabolic
@unpack neighbor_ids = interfaces
interfaces_u = interfaces.u

Expand All @@ -170,10 +170,10 @@ end
# This is the version used when calculating the divergence of the viscous fluxes
function calc_interface_flux!(surface_flux_values,
mesh::TreeMesh{1}, equations_parabolic,
dg::DG, cache_parabolic)
@unpack neighbor_ids, orientations = cache_parabolic.interfaces
dg::DG, cache)
@unpack neighbor_ids, orientations = cache.interfaces

@threaded for interface in eachinterface(dg, cache_parabolic)
@threaded for interface in eachinterface(dg, cache)
# Get neighboring elements
left_id = neighbor_ids[1, interface]
right_id = neighbor_ids[2, interface]
Expand All @@ -184,7 +184,7 @@ function calc_interface_flux!(surface_flux_values,
right_direction = 2 * orientations[interface] - 1

# Get precomputed fluxes at interfaces
flux_ll, flux_rr = get_surface_node_vars(cache_parabolic.interfaces.u,
flux_ll, flux_rr = get_surface_node_vars(cache.interfaces.u,
equations_parabolic,
dg, interface)

Expand All @@ -203,15 +203,14 @@ function calc_interface_flux!(surface_flux_values,
end

# This is the version used when calculating the divergence of the viscous fluxes
function prolong2boundaries!(cache_parabolic, flux_viscous,
function prolong2boundaries!(boundaries, flux_viscous,
mesh::TreeMesh{1},
equations_parabolic::AbstractEquationsParabolic,
surface_integral, dg::DG, cache)
@unpack boundaries = cache_parabolic
@unpack neighbor_sides, neighbor_ids = boundaries
boundaries_u = boundaries.u

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

if neighbor_sides[boundary] == 1
Expand Down Expand Up @@ -432,18 +431,18 @@ function calc_gradient!(gradients, u_transformed, t,
end

# Prolong solution to interfaces
@trixi_timeit timer() "prolong2interfaces" prolong2interfaces!(cache_parabolic,
@trixi_timeit timer() "prolong2interfaces" prolong2interfaces!(cache,
u_transformed, mesh,
equations_parabolic,
dg.surface_integral,
dg)

# Calculate interface fluxes
@trixi_timeit timer() "interface flux" begin
@unpack surface_flux_values = cache_parabolic.elements
@unpack neighbor_ids, orientations = cache_parabolic.interfaces
@unpack surface_flux_values = cache.elements
@unpack neighbor_ids, orientations = cache.interfaces

@threaded for interface in eachinterface(dg, cache_parabolic)
@threaded for interface in eachinterface(dg, cache)
# Get neighboring elements
left_id = neighbor_ids[1, interface]
right_id = neighbor_ids[2, interface]
Expand All @@ -454,7 +453,7 @@ function calc_gradient!(gradients, u_transformed, t,
right_direction = 2 * orientations[interface] - 1

# Call pointwise Riemann solver
u_ll, u_rr = get_surface_node_vars(cache_parabolic.interfaces.u,
u_ll, u_rr = get_surface_node_vars(cache.interfaces.u,
equations_parabolic, dg, interface)
flux = 0.5 * (u_ll + u_rr)

Expand All @@ -467,14 +466,14 @@ function calc_gradient!(gradients, u_transformed, t,
end

# Prolong solution to boundaries
@trixi_timeit timer() "prolong2boundaries" prolong2boundaries!(cache_parabolic,
@trixi_timeit timer() "prolong2boundaries" prolong2boundaries!(cache,
u_transformed, mesh,
equations_parabolic,
dg.surface_integral,
dg)

# Calculate boundary fluxes
@trixi_timeit timer() "boundary flux" calc_boundary_flux_gradients!(cache_parabolic,
@trixi_timeit timer() "boundary flux" calc_boundary_flux_gradients!(cache,
t,
boundary_conditions_parabolic,
mesh,
Expand All @@ -485,7 +484,7 @@ function calc_gradient!(gradients, u_transformed, t,
# Calculate surface integrals
@trixi_timeit timer() "surface integral" begin
@unpack boundary_interpolation = dg.basis
@unpack surface_flux_values = cache_parabolic.elements
@unpack surface_flux_values = cache.elements

# Note that all fluxes have been computed with outward-pointing normal vectors.
# Access the factors only once before beginning the loop to increase performance.
Expand All @@ -512,7 +511,7 @@ function calc_gradient!(gradients, u_transformed, t,
# Apply Jacobian from mapping to reference element
@trixi_timeit timer() "Jacobian" begin
apply_jacobian_parabolic!(gradients, mesh, equations_parabolic, dg,
cache_parabolic)
cache)
end

return nothing
Expand Down
Loading
Loading