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 12, 2024
1 parent 9838536 commit fb45900
Show file tree
Hide file tree
Showing 5 changed files with 258 additions and 203 deletions.
242 changes: 127 additions & 115 deletions src/solvers/dg_1d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -181,9 +181,63 @@ function volume_flux_dgfv_kernel!(volume_flux_arr, fstar1_L, fstar1_R, u,
return nothing
end

# Kernel for calculating pure DG and DG-FV volume fluxes
function volume_flux_dgfv_kernel!(volume_flux_arr, noncons_flux_arr, fstar1_L, fstar1_R, u,
element_ids_dgfv, derivative_split,
equations::AbstractEquations{1},
volume_flux_dg::Any, nonconservative_flux_dg::Any,
volume_flux_fv::Any, nonconservative_flux_fv::Any)
j = (blockIdx().x - 1) * blockDim().x + threadIdx().x
k = (blockIdx().y - 1) * blockDim().y + threadIdx().y

if (j <= size(u, 2)^2 && k <= size(u, 3))
# length(element_ids_dgfv) == size(u, 3)

j1 = div(j - 1, size(u, 2)) + 1
j2 = rem(j - 1, size(u, 2)) + 1

element_dgfv = element_ids_dgfv[k] # check if `element_dgfv` is zero

# The sets of `get_node_vars` operations may be combined
# into a single set of operation for better performance (to be explored)

u_node = get_node_vars(u, equations, j1, k)
u_node1 = get_node_vars(u, equations, j2, k)

volume_flux_node = volume_flux_dg(u_node, u_node1, 1, equations)
noncons_flux_node = nonconservative_flux_dg(u_node, u_node1, 1, equations)

@inbounds begin
for ii in axes(u, 1)
volume_flux_arr[ii, j1, j2, k] = derivative_split[j1, j2] * volume_flux_node[ii]
noncons_flux_arr[ii, j1, j2, k] = noncons_flux_node[ii]
end
end

if j1 !== 1 && j2 === 1 && element_dgfv !== 0 # bad
u_ll = get_node_vars(u, equations, j1 - 1, element_dgfv)
u_rr = get_node_vars(u, equations, j1, element_dgfv)

f1_node = volume_flux_fv(u_ll, u_rr, 1, equations)

f1_L_node = nonconservative_flux_fv(u_ll, u_rr, 1, equations)
f1_R_node = nonconservative_flux_fv(u_rr, u_ll, 1, equations)

@inbounds begin
for ii in axes(u, 1)
fstar1_L[ii, j1, element_dgfv] = f1_node[ii] + 0.5 * f1_L_node[ii]
fstar1_R[ii, j1, element_dgfv] = f1_node[ii] + 0.5 * f1_R_node[ii]
end
end
end
end

return nothing
end

# Kernel for calculating DG volume integral contribution
function volume_integral_dg_kernel!(du, element_ids_dg, element_ids_dgfv, alpha, derivative_split,
volume_flux_arr)
volume_flux_arr, equations::AbstractEquations{1})
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 @@ -216,6 +270,53 @@ function volume_integral_dg_kernel!(du, element_ids_dg, element_ids_dgfv, alpha,
return nothing
end

# Kernel for calculating DG volume integral contribution
function volume_integral_dg_kernel!(du, element_ids_dg, element_ids_dgfv, alpha, derivative_split,
volume_flux_arr, noncons_flux_arr,
equations::AbstractEquations{1})
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(du, 1) && j <= size(du, 2) && k <= size(du, 3))
# length(element_ids_dg) == size(du, 3)
# length(element_ids_dgfv) == size(du, 3)

element_dg = element_ids_dg[k] # check if `element_dg` is zero
element_dgfv = element_ids_dgfv[k] # check if `element_dgfv` is zero
alpha_element = alpha[k]

@inbounds begin
if element_dg !== 0 # bad
integral_contribution = 0.0

for ii in axes(du, 2)
du[i, j, element_dg] += volume_flux_arr[i, j, ii, element_dg]
integral_contribution += derivative_split[j, ii] *
noncons_flux_arr[i, j, ii, element_dg]
end

du[i, j, element_dg] += 0.5 * integral_contribution
end

if element_dgfv !== 0 # bad
integral_contribution = 0.0

for ii in axes(du, 2)
du[i, j, element_dgfv] += (1 - alpha_element) *
volume_flux_arr[i, j, ii, element_dgfv]
integral_contribution += derivative_split[j, ii] *
noncons_flux_arr[i, j, ii, element_dgfv]
end

du[i, j, element_dgfv] += 0.5 * (1 - alpha_element) * integral_contribution
end
end
end

return nothing
end

# Kernel for calculating FV volume integral contribution
function volume_integral_fv_kernel!(du, fstar1_L, fstar1_R, inverse_weights, element_ids_dgfv,
alpha)
Expand Down Expand Up @@ -477,7 +578,7 @@ end

# Kernel for calculating source terms
function source_terms_kernel!(du, u, node_coordinates, t, equations::AbstractEquations{1},
source_terms::Function)
source_terms::Any)
j = (blockIdx().x - 1) * blockDim().x + threadIdx().x
k = (blockIdx().y - 1) * blockDim().y + threadIdx().y

Expand Down Expand Up @@ -631,9 +732,10 @@ function cuda_volume_integral!(du, u, mesh::TreeMesh{1}, nonconservative_terms::
element_ids_dgfv,
alpha,
derivative_split,
volume_flux_arr)
volume_flux_arr,
equations)
volume_integral_dg_kernel(du, element_ids_dg, element_ids_dgfv, alpha, derivative_split,
volume_flux_arr;
volume_flux_arr, equations;
configurator_3d(volume_integral_dg_kernel, du)...)

size_arr = CuArray{Float64}(undef, size(u, 2), size(u, 3))
Expand All @@ -649,117 +751,6 @@ function cuda_volume_integral!(du, u, mesh::TreeMesh{1}, nonconservative_terms::
return nothing
end

# Kernel for calculating pure DG and DG-FV volume fluxes
function volume_flux_dgfv_kernel!(volume_flux_arr, noncons_flux_arr, fstar1_L, fstar1_R, u,
element_ids_dgfv, derivative_split,
equations::AbstractEquations{1},
volume_flux_dg::Any, nonconservative_flux_dg::Any,
volume_flux_fv::Any, nonconservative_flux_fv::Any)
j = (blockIdx().x - 1) * blockDim().x + threadIdx().x
k = (blockIdx().y - 1) * blockDim().y + threadIdx().y

if (j <= size(u, 2)^2 && k <= size(u, 3))
# length(element_ids_dgfv) == size(u, 3)

j1 = div(j - 1, size(u, 2)) + 1
j2 = rem(j - 1, size(u, 2)) + 1

element_dgfv = element_ids_dgfv[k] # check if `element_dgfv` is zero

# The sets of `get_node_vars` operations may be combined
# into a single set of operation for better performance (to be explored)

u_node = get_node_vars(u, equations, j1, k)
u_node1 = get_node_vars(u, equations, j2, k)

volume_flux_node = volume_flux_dg(u_node, u_node1, 1, equations)
noncons_flux_node = nonconservative_flux_dg(u_node, u_node1, 1, equations)

@inbounds begin
for ii in axes(u, 1)
volume_flux_arr[ii, j1, j2, k] = volume_flux_node[ii]
noncons_flux_arr[ii, j1, j2, k] = noncons_flux_node[ii] * derivative_split[j1, j2]
end
end

if j1 !== 1 && j2 === 1 && element_dgfv !== 0 # bad
u_ll = get_node_vars(u, equations, j1 - 1, element_dgfv)
u_rr = get_node_vars(u, equations, j1, element_dgfv)

f1_node = volume_flux_fv(u_ll, u_rr, 1, equations)

f1_L_node = nonconservative_flux_fv(u_ll, u_rr, 1, equations)
f1_R_node = nonconservative_flux_fv(u_rr, u_ll, 1, equations)

@inbounds begin
for ii in axes(u, 1)
fstar1_L[ii, j1, element_dgfv] = f1_node[ii] + 0.5 * f1_L_node[ii]
fstar1_R[ii, j1, element_dgfv] = f1_node[ii] + 0.5 * f1_R_node[ii]
end
end
end
end

return nothing
end

# # Kernel for calculating symmetric and nonconservative volume integrals
# function volume_integral_kernel!(du, derivative_split, symmetric_flux_arr, noncons_flux_arr)
# 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(du, 1) && j <= size(du, 2) && k <= size(du, 3))
# @inbounds begin
# integral_contribution = 0.0 # change back to `Float32`

# for ii in axes(du, 2)
# du[i, j, k] += derivative_split[j, ii] * symmetric_flux_arr[i, j, ii, k]
# integral_contribution += noncons_flux_arr[i, j, ii, k]
# end

# du[i, j, k] += 0.5 * integral_contribution # change back to `Float32`
# end
# end

# return nothing
# end

# # Kernel for calculating DG volume integral contribution
# function volume_integral_dg_kernel!(du, element_ids_dg, element_ids_dgfv, alpha, derivative_split,
# volume_flux_arr, noncons_flux_arr)
# 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(du, 1) && j <= size(du, 2) && k <= size(du, 3))
# # length(element_ids_dg) == size(du, 3)
# # length(element_ids_dgfv) == size(du, 3)

# element_dg = element_ids_dg[k] # check if `element_dg` is zero
# element_dgfv = element_ids_dgfv[k] # check if `element_dgfv` is zero
# alpha_element = alpha[k]

# @inbounds begin
# if element_dg !== 0 # bad
# for ii in axes(du, 2)
# du[i, j, element_dg] += derivative_split[j, ii] *
# volume_flux_arr[i, j, ii, element_dg]
# end
# end

# if element_dgfv !== 0 # bad
# for ii in axes(du, 2)
# du[i, j, element_dgfv] += (1 - alpha_element) * derivative_split[j, ii] *
# volume_flux_arr[i, j, ii, element_dgfv]
# end
# end
# end
# end

# return nothing
# end

# Pack kernels for calculating volume integrals
function cuda_volume_integral!(du, u, mesh::TreeMesh{1}, nonconservative_terms::True, equations,
volume_integral::VolumeIntegralShockCapturingHG, dg::DGSEM, cache)
Expand Down Expand Up @@ -812,6 +803,27 @@ function cuda_volume_integral!(du, u, mesh::TreeMesh{1}, nonconservative_terms::
nonconservative_flux_dg, volume_flux_fv, nonconservative_flux_fv;
configurator_2d(volume_flux_dgfv_kernel, size_arr)...)

volume_integral_dg_kernel = @cuda launch=false volume_integral_dg_kernel!(du, element_ids_dg,
element_ids_dgfv,
alpha,
derivative_split,
volume_flux_arr,
noncons_flux_arr,
equations)
volume_integral_dg_kernel(du, element_ids_dg, element_ids_dgfv, alpha, derivative_split,
volume_flux_arr, noncons_flux_arr, equations;
configurator_3d(volume_integral_dg_kernel, du)...)

size_arr = CuArray{Float64}(undef, size(u, 2), size(u, 3))

volume_integral_fv_kernel = @cuda launch=false volume_integral_fv_kernel!(du, fstar1_L,
fstar1_R,
inverse_weights,
element_ids_dgfv,
alpha)
volume_integral_fv_kernel(du, fstar1_L, fstar1_R, inverse_weights, element_ids_dgfv, alpha;
configurator_2d(volume_integral_fv_kernel, size_arr)...)

return nothing
end

Expand Down
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 @@ -240,7 +240,8 @@ end

# Kernel for calculating DG volume integral contribution
function volume_integral_dg_kernel!(du, element_ids_dg, element_ids_dgfv, alpha, derivative_split,
volume_flux_arr1, volume_flux_arr2)
volume_flux_arr1, volume_flux_arr2,
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 @@ -979,9 +980,10 @@ function cuda_volume_integral!(du, u, mesh::TreeMesh{2}, nonconservative_terms::
alpha,
derivative_split,
volume_flux_arr1,
volume_flux_arr2)
volume_flux_arr2,
equations)
volume_integral_dg_kernel(du, element_ids_dg, element_ids_dgfv, alpha, derivative_split,
volume_flux_arr1, volume_flux_arr2;
volume_flux_arr1, volume_flux_arr2, equations;
configurator_3d(volume_integral_dg_kernel, size_arr)...)

size_arr = CuArray{Float64}(undef, size(u, 2)^2, size(u, 4))
Expand Down
8 changes: 5 additions & 3 deletions src/solvers/dg_3d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -208,7 +208,8 @@ end

# Kernel for calculating DG volume integral contribution
function volume_integral_dg_kernel!(du, element_ids_dg, element_ids_dgfv, alpha, derivative_split,
volume_flux_arr1, volume_flux_arr2, volume_flux_arr3)
volume_flux_arr1, volume_flux_arr2, volume_flux_arr3,
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
Expand Down Expand Up @@ -1061,9 +1062,10 @@ function cuda_volume_integral!(du, u, mesh::TreeMesh{3}, nonconservative_terms::
derivative_split,
volume_flux_arr1,
volume_flux_arr2,
volume_flux_arr3)
volume_flux_arr3,
equations)
volume_integral_dg_kernel(du, element_ids_dg, element_ids_dgfv, alpha, derivative_split,
volume_flux_arr1, volume_flux_arr2, volume_flux_arr3;
volume_flux_arr1, volume_flux_arr2, volume_flux_arr3, equations;
configurator_3d(volume_integral_dg_kernel, size_arr)...)

size_arr = CuArray{Float64}(undef, size(u, 2)^3, size(u, 5))
Expand Down
Loading

0 comments on commit fb45900

Please sign in to comment.