Skip to content

Commit

Permalink
Complete 1D
Browse files Browse the repository at this point in the history
  • Loading branch information
@huiyuxie
huiyuxie committed Sep 7, 2024
1 parent b208541 commit b61d915
Show file tree
Hide file tree
Showing 5 changed files with 158 additions and 46 deletions.
87 changes: 85 additions & 2 deletions src/solvers/dg_1d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -283,6 +283,7 @@ function boundary_flux_kernel!(surface_flux_values, boundaries_u, node_coordinat

@inbounds begin
for ii in axes(surface_flux_values, 1)
# `boundary_flux_node` can be nothing if periodic boundary condition is applied
surface_flux_values[ii, direction, neighbor] = boundary_flux_node === nothing ? # bad
surface_flux_values[ii, direction,
neighbor] :
Expand All @@ -294,6 +295,46 @@ function boundary_flux_kernel!(surface_flux_values, boundaries_u, node_coordinat
return nothing
end

function boundary_flux_kernel!(surface_flux_values, boundaries_u, node_coordinates, t, boundary_arr,
indices_arr, neighbor_ids, neighbor_sides, orientations,
boundary_conditions::NamedTuple, equations::AbstractEquations{1},
surface_flux::Any, nonconservative_flux::Any)
k = (blockIdx().x - 1) * blockDim().x + threadIdx().x

if (k <= length(boundary_arr))
boundary = boundary_arr[k]
direction = (indices_arr[1] <= boundary) + (indices_arr[2] <= boundary)

neighbor = neighbor_ids[boundary]
side = neighbor_sides[boundary]
orientation = orientations[boundary]

u_ll, u_rr = get_surface_node_vars(boundaries_u, equations, boundary)
u_inner = (2 - side) * u_ll + (side - 1) * u_rr
x = get_node_coords(node_coordinates, equations, boundary)

# TODO: Improve this part
if direction == 1
u_boundary = boundary_conditions[1].boundary_value_function(x, t, equations)
flux_node = surface_flux(u_boundary, u_inner, orientation, equations)
noncons_flux_node = nonconservative_flux(u_boundary, u_inner, orientation, equations)
else
u_boundary = boundary_conditions[2].boundary_value_function(x, t, equations)
flux_node = surface_flux(u_inner, u_boundary, orientation, equations)
noncons_flux_node = nonconservative_flux(u_inner, u_boundary, orientation, equations)
end

@inbounds begin
for ii in axes(surface_flux_values, 1)
surface_flux_values[ii, direction, neighbor] = flux_node[ii] +
0.5 * noncons_flux_node[ii]
end
end
end

return nothing
end

# Kernel for calculating surface integrals
function surface_integral_kernel!(du, factor_arr, surface_flux_values,
equations::AbstractEquations{1})
Expand Down Expand Up @@ -594,6 +635,48 @@ end

function cuda_boundary_flux!(t, mesh::TreeMesh{1}, boundary_conditions::NamedTuple,
nonconservative_terms::True, equations, dg::DGSEM, cache)
surface_flux, nonconservative_flux = dg.surface_integral.surface_flux

n_boundaries_per_direction = CuArray{Int64}(cache.boundaries.n_boundaries_per_direction)
neighbor_ids = CuArray{Int64}(cache.boundaries.neighbor_ids)
neighbor_sides = CuArray{Int64}(cache.boundaries.neighbor_sides)
orientations = CuArray{Int64}(cache.boundaries.orientations)
boundaries_u = CuArray{Float64}(cache.boundaries.u)
node_coordinates = CuArray{Float64}(cache.boundaries.node_coordinates)
surface_flux_values = CuArray{Float64}(cache.elements.surface_flux_values)

lasts = zero(n_boundaries_per_direction)
firsts = zero(n_boundaries_per_direction)

last_first_indices_kernel = @cuda launch=false last_first_indices_kernel!(lasts, firsts,
n_boundaries_per_direction)
last_first_indices_kernel(lasts, firsts, n_boundaries_per_direction;
configurator_1d(last_first_indices_kernel, lasts)...)

lasts, firsts = Array(lasts), Array(firsts)
boundary_arr = CuArray{Int64}(firsts[1]:lasts[2])
indices_arr = CuArray{Int64}([firsts[1], firsts[2]])

# Replace with callable functions (not necessary here)
# boundary_conditions_callable = replace_boundary_conditions(boundary_conditions)

boundary_flux_kernel = @cuda launch=false boundary_flux_kernel!(surface_flux_values,
boundaries_u, node_coordinates,
t, boundary_arr, indices_arr,
neighbor_ids, neighbor_sides,
orientations,
boundary_conditions,
equations,
surface_flux,
nonconservative_flux)
boundary_flux_kernel(surface_flux_values, boundaries_u, node_coordinates, t, boundary_arr,
indices_arr, neighbor_ids, neighbor_sides, orientations,
boundary_conditions, equations, surface_flux, nonconservative_flux;
configurator_1d(boundary_flux_kernel, boundary_arr)...)

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

return nothing
end

# Pack kernels for calculating surface integrals
Expand Down Expand Up @@ -642,9 +725,9 @@ function cuda_sources!(du, u, t, source_terms, equations::AbstractEquations{1},
return nothing
end

# Put everything together into a single function
# Ref: `rhs!` function in Trixi.jl
# Put everything together into a single function.

# Ref: `rhs!` function in Trixi.jl
function rhs_gpu!(du_cpu, u_cpu, t, mesh::TreeMesh{1}, equations, boundary_conditions,
source_terms::Source, dg::DGSEM, cache) where {Source}
du, u = copy_to_device!(du_cpu, u_cpu)
Expand Down
5 changes: 3 additions & 2 deletions src/solvers/dg_2d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -376,6 +376,7 @@ function boundary_flux_kernel!(surface_flux_values, boundaries_u, node_coordinat

@inbounds begin
for ii in axes(surface_flux_values, 1)
# `boundary_flux_node` can be nothing if periodic boundary condition is applied
surface_flux_values[ii, j, direction, neighbor] = boundary_flux_node === nothing ? # bad
surface_flux_values[ii, j,
direction,
Expand Down Expand Up @@ -1067,9 +1068,9 @@ function cuda_sources!(du, u, t, source_terms, equations::AbstractEquations{2},
return nothing
end

# Put everything together into a single function
# Ref: `rhs!` function in Trixi.jl
# Put everything together into a single function.

# Ref: `rhs!` function in Trixi.jl
function rhs_gpu!(du_cpu, u_cpu, t, mesh::TreeMesh{2}, equations, boundary_conditions,
source_terms::Source, dg::DGSEM, cache) where {Source}
du, u = copy_to_device!(du_cpu, u_cpu)
Expand Down
5 changes: 3 additions & 2 deletions src/solvers/dg_3d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -297,6 +297,7 @@ function boundary_flux_kernel!(surface_flux_values, boundaries_u, node_coordinat

@inbounds begin
for ii in axes(surface_flux_values, 1)
# `boundary_flux_node` can be nothing if periodic boundary condition is applied
surface_flux_values[ii, j1, j2, direction, neighbor] = boundary_flux_node ===
nothing ? # bad
surface_flux_values[ii, j1,
Expand Down Expand Up @@ -1141,9 +1142,9 @@ function cuda_sources!(du, u, t, source_terms, equations::AbstractEquations{3},
return nothing
end

# Put everything together into a single function
# Ref: `rhs!` function in Trixi.jl
# Put everything together into a single function.

# Ref: `rhs!` function in Trixi.jl
function rhs_gpu!(du_cpu, u_cpu, t, mesh::TreeMesh{3}, equations, boundary_conditions,
source_terms::Source, dg::DGSEM, cache) where {Source}
du, u = copy_to_device!(du_cpu, u_cpu)
Expand Down
27 changes: 27 additions & 0 deletions test.jl
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,27 @@
using Trixi

equations = IdealGlmMhdMulticomponentEquations1D(gammas = (2.0, 2.0, 2.0),
gas_constants = (2.0, 2.0, 2.0))

initial_condition = initial_condition_weak_blast_wave

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

coordinates_min = 0.0
coordinates_max = 1.0
mesh = TreeMesh(coordinates_min, coordinates_max,
initial_refinement_level = 4,
n_cells_max = 10_000)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

(; mesh, equations, initial_condition, boundary_conditions, source_terms, solver, cache) = semi

surface_flux_values = cache.elements.surface_flux_values
n_boundaries_per_direction = cache.boundaries.n_boundaries_per_direction

# Calculate indices
lasts = accumulate(+, n_boundaries_per_direction)
firsts = lasts - n_boundaries_per_direction .+ 1
80 changes: 40 additions & 40 deletions test/test_script.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -115,43 +115,43 @@ surface_flux_values_gpu = replace(cache_gpu.elements.surface_flux_values, NaN =>
surface_flux_values = replace(cache.elements.surface_flux_values, NaN => 0.0)
@test surface_flux_values_gpu ≈ 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)
# # u_upper_gpu = replace(cache_gpu.mortars.u_upper, NaN => 0.0)
# # u_lower_gpu = replace(cache_gpu.mortars.u_lower, NaN => 0.0)
# # u_upper = replace(cache.mortars.u_upper, NaN => 0.0)
# # u_lower = replace(cache.mortars.u_lower, NaN => 0.0)
# # @test u_upper_gpu ≈ u_upper
# # @test u_lower_gpu ≈ u_lower

# # # 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)
# # surface_flux_values_gpu = replace(cache_gpu.elements.surface_flux_values, NaN => 0.0)
# # surface_flux_values = replace(cache.elements.surface_flux_values, NaN => 0.0)
# # @test surface_flux_values_gpu ≈ 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 CUDA.@allowscalar 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 CUDA.@allowscalar 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 CUDA.@allowscalar du_gpu ≈ du

# # Copy data back to host
# # du_cpu, u_cpu = TrixiGPU.copy_to_host!(du_gpu, u_gpu)
# # 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)
# u_upper_gpu = replace(cache_gpu.mortars.u_upper, NaN => 0.0)
# u_lower_gpu = replace(cache_gpu.mortars.u_lower, NaN => 0.0)
# u_upper = replace(cache.mortars.u_upper, NaN => 0.0)
# u_lower = replace(cache.mortars.u_lower, NaN => 0.0)
# @test u_upper_gpu ≈ u_upper
# @test u_lower_gpu ≈ u_lower

# # 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)
# surface_flux_values_gpu = replace(cache_gpu.elements.surface_flux_values, NaN => 0.0)
# surface_flux_values = replace(cache.elements.surface_flux_values, NaN => 0.0)
# @test surface_flux_values_gpu ≈ 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 CUDA.@allowscalar 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 CUDA.@allowscalar 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 CUDA.@allowscalar du_gpu ≈ du

# Copy data back to host
# du_cpu, u_cpu = TrixiGPU.copy_to_host!(du_gpu, u_gpu)

0 comments on commit b61d915

Please sign in to comment.