Skip to content

Commit

Permalink
Add interface flux 3D
Browse files Browse the repository at this point in the history
  • Loading branch information
@huiyuxie
huiyuxie committed Sep 15, 2024
1 parent 862b658 commit 1e9cd92
Show file tree
Hide file tree
Showing 3 changed files with 159 additions and 55 deletions.
8 changes: 5 additions & 3 deletions src/solvers/dg_2d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -419,7 +419,8 @@ end

# Kernel for setting interface fluxes
function interface_flux_kernel!(surface_flux_values, surface_flux_arr, noncons_left_arr,
noncons_right_arr, neighbor_ids, orientations)
noncons_right_arr, neighbor_ids, orientations,
equations::AbstractEquations{2})
i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
j = (blockIdx().y - 1) * blockDim().y + threadIdx().y
k = (blockIdx().z - 1) * blockDim().z + threadIdx().z
Expand Down Expand Up @@ -1098,9 +1099,10 @@ function cuda_interface_flux!(mesh::TreeMesh{2}, nonconservative_terms::True, eq
surface_flux_arr,
noncons_left_arr,
noncons_right_arr,
neighbor_ids, orientations)
neighbor_ids, orientations,
equations)
interface_flux_kernel(surface_flux_values, surface_flux_arr, noncons_left_arr,
noncons_right_arr, neighbor_ids, orientations;
noncons_right_arr, neighbor_ids, orientations, equations;
configurator_3d(interface_flux_kernel, size_arr)...)

cache.elements.surface_flux_values = surface_flux_values # copy back to host automatically
Expand Down
104 changes: 100 additions & 4 deletions src/solvers/dg_3d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -463,12 +463,13 @@ function interface_flux_kernel!(surface_flux_values, surface_flux_arr, neighbor_
left_id = neighbor_ids[1, k]
right_id = neighbor_ids[2, k]

left_dir = 2 * orientations[k]
right_dir = 2 * orientations[k] - 1
left_direction = 2 * orientations[k]
right_direction = 2 * orientations[k] - 1

@inbounds begin
surface_flux_values[i, j1, j2, left_dir, left_id] = surface_flux_arr[i, j1, j2, k]
surface_flux_values[i, j1, j2, right_dir, right_id] = surface_flux_arr[i, j1, j2, k]
surface_flux_values[i, j1, j2, left_direction, left_id] = surface_flux_arr[i, j1, j2, k]
surface_flux_values[i, j1, j2, right_direction, right_id] = surface_flux_arr[i, j1, j2,
k]
end
end

Expand Down Expand Up @@ -1283,6 +1284,70 @@ function cuda_interface_flux!(mesh::TreeMesh{3}, nonconservative_terms::False, e
return nothing
end

# Kernel for calculating surface and both nonconservative fluxes
function surface_noncons_flux_kernel!(surface_flux_arr, noncons_left_arr, noncons_right_arr,
interfaces_u, orientations, equations::AbstractEquations{3},
surface_flux::Any, nonconservative_flux::Any)
j1 = (blockIdx().x - 1) * blockDim().x + threadIdx().x
j2 = (blockIdx().y - 1) * blockDim().y + threadIdx().y
k = (blockIdx().z - 1) * blockDim().z + threadIdx().z

if (j1 <= size(surface_flux_arr, 2) && j2 <= size(surface_flux_arr, 3) &&
k <= size(surface_flux_arr, 4))
u_ll, u_rr = get_surface_node_vars(interfaces_u, equations, j1, j2, k)
orientation = orientations[k]

surface_flux_node = surface_flux(u_ll, u_rr, orientation, equations)
noncons_left_node = nonconservative_flux(u_ll, u_rr, orientation, equations)
noncons_right_node = nonconservative_flux(u_rr, u_ll, orientation, equations)

@inbounds begin
for ii in axes(surface_flux_arr, 1)
surface_flux_arr[ii, j1, j2, k] = surface_flux_node[ii]
noncons_left_arr[ii, j1, j2, k] = noncons_left_node[ii]
noncons_right_arr[ii, j1, j2, k] = noncons_right_node[ii]
end
end
end

return nothing
end

# Kernel for setting interface fluxes
function interface_flux_kernel!(surface_flux_values, surface_flux_arr, noncons_left_arr,
noncons_right_arr, neighbor_ids, orientations,
equations::AbstractEquations{3})
i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
j = (blockIdx().y - 1) * blockDim().y + threadIdx().y
k = (blockIdx().z - 1) * blockDim().z + threadIdx().z

if (i <= size(surface_flux_values, 1) && j <= size(surface_flux_arr, 2)^2 &&
k <= size(surface_flux_arr, 4))
j1 = div(j - 1, size(surface_flux_arr, 2)) + 1
j2 = rem(j - 1, size(surface_flux_arr, 2)) + 1

left_id = neighbor_ids[1, k]
right_id = neighbor_ids[2, k]

left_direction = 2 * orientations[k]
right_direction = 2 * orientations[k] - 1

@inbounds begin
surface_flux_values[i, j1, j2, left_direction, left_id] = surface_flux_arr[i, j1, j2,
k] +
0.5 *
noncons_left_arr[i, j1, j2, k]
surface_flux_values[i, j1, j2, right_direction, right_id] = surface_flux_arr[i, j1, j2,
k] +
0.5 *
noncons_right_arr[i, j1, j2,
k]
end
end

return nothing
end

# Pack kernels for calculating interface fluxes
function cuda_interface_flux!(mesh::TreeMesh{3}, nonconservative_terms::True, equations, dg::DGSEM,
cache)
Expand All @@ -1296,6 +1361,37 @@ function cuda_interface_flux!(mesh::TreeMesh{3}, nonconservative_terms::True, eq
noncons_right_arr = CuArray{Float64}(undef, size(interfaces_u)[2:end]...)
surface_flux_values = CuArray{Float64}(cache.elements.surface_flux_values)

size_arr = CuArray{Float64}(undef, size(interfaces_u, 3), size(interfaces_u, 4),
size(interfaces_u, 5))

surface_noncons_flux_kernel = @cuda launch=false surface_noncons_flux_kernel!(surface_flux_arr,
noncons_left_arr,
noncons_right_arr,
interfaces_u,
orientations,
equations,
surface_flux,
nonconservative_flux)
surface_noncons_flux_kernel(surface_flux_arr, noncons_left_arr, noncons_right_arr, interfaces_u,
orientations, equations, surface_flux, nonconservative_flux;
configurator_3d(surface_noncons_flux_kernel, size_arr)...)

size_arr = CuArray{Float64}(undef, size(surface_flux_values, 1), size(interfaces_u, 3)^2,
size(interfaces_u, 5))

interface_flux_kernel = @cuda launch=false interface_flux_kernel!(surface_flux_values,
surface_flux_arr,
noncons_left_arr,
noncons_right_arr,
neighbor_ids, orientations,
equations)
interface_flux_kernel(surface_flux_values, surface_flux_arr, noncons_left_arr,
noncons_right_arr,
neighbor_ids, orientations, equations;
configurator_3d(interface_flux_kernel, size_arr)...)

cache.elements.surface_flux_values = surface_flux_values # copy back to host automatically

return nothing
end

Expand Down
102 changes: 54 additions & 48 deletions test/test_script.jl
View file
Original file line number Diff line number Diff line change
@@ -1,18 +1,24 @@
include("test_trixigpu.jl")

advection_velocity = (0.2, -0.7, 0.5)
equations = LinearScalarAdvectionEquation3D(advection_velocity)
equations = IdealGlmMhdEquations3D(5 / 3)

solver = DGSEM(polydeg = 3, surface_flux = flux_lax_friedrichs)
initial_condition = initial_condition_convergence_test

volume_flux = (flux_hindenlang_gassner, flux_nonconservative_powell)
solver = DGSEM(polydeg = 3,
surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell),
volume_integral = VolumeIntegralFluxDifferencing(volume_flux))

coordinates_min = (-1.0, -1.0, -1.0)
coordinates_max = (1.0, 1.0, 1.0)
mesh = TreeMesh(coordinates_min, coordinates_max,
initial_refinement_level = 2,
n_cells_max = 10_000)

mesh = TreeMesh(coordinates_min, coordinates_max, initial_refinement_level = 3,
n_cells_max = 30_000)
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

semi = SemidiscretizationHyperbolic(mesh, equations,
initial_condition_convergence_test, solver)
###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 1.0)

Expand Down Expand Up @@ -68,44 +74,44 @@ TrixiGPU.cuda_prolong2boundaries!(u_gpu, mesh_gpu, boundary_conditions_gpu, equa
Trixi.prolong2boundaries!(cache, u, mesh, equations, solver.surface_integral, solver)
@test_approx cache_gpu.boundaries.u ≈ cache.boundaries.u

# Test `cuda_boundary_flux!`
TrixiGPU.cuda_boundary_flux!(t_gpu, mesh_gpu, boundary_conditions_gpu,
Trixi.have_nonconservative_terms(equations_gpu), equations_gpu,
solver_gpu, cache_gpu)
Trixi.calc_boundary_flux!(cache, t, boundary_conditions, mesh, equations,
solver.surface_integral, solver)
@test_approx cache_gpu.elements.surface_flux_values ≈ cache.elements.surface_flux_values

# Test `cuda_prolong2mortars!`
TrixiGPU.cuda_prolong2mortars!(u_gpu, mesh_gpu, TrixiGPU.check_cache_mortars(cache_gpu),
solver_gpu, cache_gpu)
Trixi.prolong2mortars!(cache, u, mesh, equations,
solver.mortar, solver.surface_integral, solver)
@test_approx cache_gpu.mortars.u_upper_left ≈ cache.mortars.u_upper_left
@test_approx cache_gpu.mortars.u_upper_right ≈ cache.mortars.u_upper_right
@test_approx cache_gpu.mortars.u_lower_left ≈ cache.mortars.u_lower_left
@test_approx cache_gpu.mortars.u_lower_right ≈ cache.mortars.u_lower_right

# Test `cuda_mortar_flux!`
TrixiGPU.cuda_mortar_flux!(mesh_gpu, TrixiGPU.check_cache_mortars(cache_gpu),
Trixi.have_nonconservative_terms(equations_gpu), equations_gpu,
solver_gpu, cache_gpu)
Trixi.calc_mortar_flux!(cache.elements.surface_flux_values, mesh,
Trixi.have_nonconservative_terms(equations), equations,
solver.mortar, solver.surface_integral, solver, cache)
@test_approx cache_gpu.elements.surface_flux_values ≈ cache.elements.surface_flux_values

# Test `cuda_surface_integral!`
TrixiGPU.cuda_surface_integral!(du_gpu, mesh_gpu, equations_gpu, solver_gpu, cache_gpu)
Trixi.calc_surface_integral!(du, u, mesh, equations, solver.surface_integral, solver, cache)
@test_approx du_gpu ≈ du

# Test `cuda_jacobian!`
TrixiGPU.cuda_jacobian!(du_gpu, mesh_gpu, equations_gpu, cache_gpu)
Trixi.apply_jacobian!(du, mesh, equations, solver, cache)
@test_approx du_gpu ≈ du

# Test `cuda_sources!`
TrixiGPU.cuda_sources!(du_gpu, u_gpu, t_gpu, source_terms_gpu, equations_gpu, cache_gpu)
Trixi.calc_sources!(du, u, t, source_terms, equations, solver, cache)
@test_approx du_gpu ≈ du
# # Test `cuda_boundary_flux!`
# TrixiGPU.cuda_boundary_flux!(t_gpu, mesh_gpu, boundary_conditions_gpu,
# Trixi.have_nonconservative_terms(equations_gpu), equations_gpu,
# solver_gpu, cache_gpu)
# Trixi.calc_boundary_flux!(cache, t, boundary_conditions, mesh, equations,
# solver.surface_integral, solver)
# @test_approx cache_gpu.elements.surface_flux_values ≈ cache.elements.surface_flux_values

# # Test `cuda_prolong2mortars!`
# TrixiGPU.cuda_prolong2mortars!(u_gpu, mesh_gpu, TrixiGPU.check_cache_mortars(cache_gpu),
# solver_gpu, cache_gpu)
# Trixi.prolong2mortars!(cache, u, mesh, equations,
# solver.mortar, solver.surface_integral, solver)
# @test_approx cache_gpu.mortars.u_upper_left ≈ cache.mortars.u_upper_left
# @test_approx cache_gpu.mortars.u_upper_right ≈ cache.mortars.u_upper_right
# @test_approx cache_gpu.mortars.u_lower_left ≈ cache.mortars.u_lower_left
# @test_approx cache_gpu.mortars.u_lower_right ≈ cache.mortars.u_lower_right

# # Test `cuda_mortar_flux!`
# TrixiGPU.cuda_mortar_flux!(mesh_gpu, TrixiGPU.check_cache_mortars(cache_gpu),
# Trixi.have_nonconservative_terms(equations_gpu), equations_gpu,
# solver_gpu, cache_gpu)
# Trixi.calc_mortar_flux!(cache.elements.surface_flux_values, mesh,
# Trixi.have_nonconservative_terms(equations), equations,
# solver.mortar, solver.surface_integral, solver, cache)
# @test_approx cache_gpu.elements.surface_flux_values ≈ cache.elements.surface_flux_values

# # Test `cuda_surface_integral!`
# TrixiGPU.cuda_surface_integral!(du_gpu, mesh_gpu, equations_gpu, solver_gpu, cache_gpu)
# Trixi.calc_surface_integral!(du, u, mesh, equations, solver.surface_integral, solver, cache)
# @test_approx du_gpu ≈ du

# # Test `cuda_jacobian!`
# TrixiGPU.cuda_jacobian!(du_gpu, mesh_gpu, equations_gpu, cache_gpu)
# Trixi.apply_jacobian!(du, mesh, equations, solver, cache)
# @test_approx du_gpu ≈ du

# # Test `cuda_sources!`
# TrixiGPU.cuda_sources!(du_gpu, u_gpu, t_gpu, source_terms_gpu, equations_gpu, cache_gpu)
# Trixi.calc_sources!(du, u, t, source_terms, equations, solver, cache)
# @test_approx du_gpu ≈ du

0 comments on commit 1e9cd92

Please sign in to comment.