Skip to content

Commit

Permalink
Replace used variables
Browse files Browse the repository at this point in the history
  • Loading branch information
@huiyuxie
huiyuxie committed Sep 4, 2024
1 parent ba4c64b commit eab1329
Show file tree
Hide file tree
Showing 2 changed files with 89 additions and 59 deletions.
40 changes: 23 additions & 17 deletions src/solvers/dg_2d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,9 +57,11 @@ function volume_flux_kernel!(volume_flux_arr1, volume_flux_arr2, u, equations::A
k = (blockIdx().y - 1) * blockDim().y + threadIdx().y

if (j <= size(u, 2)^3 && k <= size(u, 4))
j1 = div(j - 1, size(u, 2)^2) + 1
j2 = div(rem(j - 1, size(u, 2)^2), size(u, 2)) + 1
j3 = rem(rem(j - 1, size(u, 2)^2), size(u, 2)) + 1
u2 = size(u, 2)

j1 = div(j - 1, u2^2) + 1
j2 = div(rem(j - 1, u2^2), u2) + 1
j3 = rem(rem(j - 1, u2^2), u2) + 1

u_node = get_node_vars(u, equations, j1, j2, k)
u_node1 = get_node_vars(u, equations, j3, j2, k)
Expand Down Expand Up @@ -88,9 +90,11 @@ function symmetric_noncons_flux_kernel!(symmetric_flux_arr1, symmetric_flux_arr2
k = (blockIdx().y - 1) * blockDim().y + threadIdx().y

if (j <= size(u, 2)^3 && k <= size(u, 4))
j1 = div(j - 1, size(u, 2)^2) + 1
j2 = div(rem(j - 1, size(u, 2)^2), size(u, 2)) + 1
j3 = rem(rem(j - 1, size(u, 2)^2), size(u, 2)) + 1
u2 = size(u, 2)

j1 = div(j - 1, u2^2) + 1
j2 = div(rem(j - 1, u2^2), u2) + 1
j3 = rem(rem(j - 1, u2^2), u2) + 1

u_node = get_node_vars(u, equations, j1, j2, k)
u_node1 = get_node_vars(u, equations, j3, j2, k)
Expand Down Expand Up @@ -176,15 +180,16 @@ function prolong_interfaces_kernel!(interfaces_u, u, neighbor_ids, orientations,
k = (blockIdx().y - 1) * blockDim().y + threadIdx().y

if (j <= size(interfaces_u, 2) * size(interfaces_u, 3) && k <= size(interfaces_u, 4))
j1 = div(j - 1, size(interfaces_u, 3)) + 1
j2 = rem(j - 1, size(interfaces_u, 3)) + 1
# size(interfaces_u, 3) == size(u, 2)
u2 = size(u, 2)

j1 = div(j - 1, u2) + 1
j2 = rem(j - 1, u2) + 1

orientation = orientations[k]
left_element = neighbor_ids[1, k]
right_element = neighbor_ids[2, k]

u2 = size(u, 2)

@inbounds begin
interfaces_u[1, j1, j2, k] = u[j1,
(2 - orientation) * u2 + (orientation - 1) * j2,
Expand Down Expand Up @@ -310,15 +315,16 @@ function prolong_boundaries_kernel!(boundaries_u, u, neighbor_ids, neighbor_side
k = (blockIdx().y - 1) * blockDim().y + threadIdx().y

if (j <= size(boundaries_u, 2) * size(boundaries_u, 3) && k <= size(boundaries_u, 4))
j1 = div(j - 1, size(boundaries_u, 3)) + 1
j2 = rem(j - 1, size(boundaries_u, 3)) + 1
# size(boundaries_u, 3) == size(u, 2)
u2 = size(u, 2)

j1 = div(j - 1, u2) + 1
j2 = rem(j - 1, u2) + 1

element = neighbor_ids[k]
side = neighbor_sides[k]
orientation = orientations[k]

u2 = size(u, 2)

@inbounds begin
boundaries_u[1, j1, j2, k] = u[j1,
(2 - orientation) * u2 + (orientation - 1) * j2,
Expand Down Expand Up @@ -557,11 +563,11 @@ function surface_integral_kernel!(du, factor_arr, surface_flux_values,
k = (blockIdx().z - 1) * blockDim().z + threadIdx().z

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

u2 = size(du, 2)

j1 = div(j - 1, u2) + 1
j2 = rem(j - 1, u2) + 1

@inbounds begin
du[i, j1, j2, k] -= (surface_flux_values[i, j2, 1, k] * isequal(j1, 1) +
surface_flux_values[i, j1, 3, k] * isequal(j2, 1)) * factor_arr[1]
Expand Down
108 changes: 66 additions & 42 deletions src/solvers/dg_3d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,9 +11,11 @@ function flux_kernel!(flux_arr1, flux_arr2, flux_arr3, u, equations::AbstractEqu
k = (blockIdx().y - 1) * blockDim().y + threadIdx().y

if (j <= size(u, 2)^3 && k <= size(u, 5))
j1 = div(j - 1, size(u, 2)^2) + 1
j2 = div(rem(j - 1, size(u, 2)^2), size(u, 2)) + 1
j3 = rem(rem(j - 1, size(u, 2)^2), size(u, 2)) + 1
u2 = size(u, 2)

j1 = div(j - 1, u2^2) + 1
j2 = div(rem(j - 1, u2^2), u2) + 1
j3 = rem(rem(j - 1, u2^2), u2) + 1

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

Expand All @@ -40,9 +42,12 @@ function weak_form_kernel!(du, derivative_dhat, flux_arr1, flux_arr2, flux_arr3)
k = (blockIdx().z - 1) * blockDim().z + threadIdx().z

if (i <= size(du, 1) && j <= size(du, 2)^3 && k <= size(du, 5))
j1 = div(j - 1, size(du, 2)^2) + 1
j2 = div(rem(j - 1, size(du, 2)^2), size(du, 2)) + 1
j3 = rem(rem(j - 1, size(du, 2)^2), size(du, 2)) + 1
# size(du, 2) == size(u, 2)
u2 = size(du, 2)

j1 = div(j - 1, u2^2) + 1
j2 = div(rem(j - 1, u2^2), u2) + 1
j3 = rem(rem(j - 1, u2^2), u2) + 1

@inbounds begin
for ii in axes(du, 2)
Expand All @@ -63,10 +68,12 @@ function volume_flux_kernel!(volume_flux_arr1, volume_flux_arr2, volume_flux_arr
k = (blockIdx().y - 1) * blockDim().y + threadIdx().y

if (j <= size(u, 2)^4 && k <= size(u, 5))
j1 = div(j - 1, size(u, 2)^3) + 1
j2 = div(rem(j - 1, size(u, 2)^3), size(u, 2)^2) + 1
j3 = div(rem(j - 1, size(u, 2)^2), size(u, 2)) + 1
j4 = rem(j - 1, size(u, 2)) + 1
u2 = size(u, 2)

j1 = div(j - 1, u2^3) + 1
j2 = div(rem(j - 1, u2^3), u2^2) + 1
j3 = div(rem(j - 1, u2^2), u2) + 1
j4 = rem(j - 1, u2) + 1

u_node = get_node_vars(u, equations, j1, j2, j3, k)
u_node1 = get_node_vars(u, equations, j4, j2, j3, k)
Expand Down Expand Up @@ -97,9 +104,12 @@ function volume_integral_kernel!(du, derivative_split, volume_flux_arr1, volume_
k = (blockIdx().z - 1) * blockDim().z + threadIdx().z

if (i <= size(du, 1) && j <= size(du, 2)^3 && k <= size(du, 5))
j1 = div(j - 1, size(du, 2)^2) + 1
j2 = div(rem(j - 1, size(du, 2)^2), size(du, 2)) + 1
j3 = rem(rem(j - 1, size(du, 2)^2), size(du, 2)) + 1
# size(du, 2) == size(u, 2)
u2 = size(du, 2)

j1 = div(j - 1, u2^2) + 1
j2 = div(rem(j - 1, u2^2), u2) + 1
j3 = rem(rem(j - 1, u2^2), u2) + 1

@inbounds begin
for ii in axes(du, 2)
Expand All @@ -123,16 +133,17 @@ function prolong_interfaces_kernel!(interfaces_u, u, neighbor_ids, orientations,
k = (blockIdx().y - 1) * blockDim().y + threadIdx().y

if (j <= size(interfaces_u, 2) * size(interfaces_u, 3)^2 && k <= size(interfaces_u, 5))
j1 = div(j - 1, size(interfaces_u, 3)^2) + 1
j2 = div(rem(j - 1, size(interfaces_u, 3)^2), size(interfaces_u, 3)) + 1
j3 = rem(rem(j - 1, size(interfaces_u, 3)^2), size(interfaces_u, 3)) + 1
# size(interfaces_u, 3) == size(u, 2)
u2 = size(u, 2)

j1 = div(j - 1, u2^2) + 1
j2 = div(rem(j - 1, u2^2), u2) + 1
j3 = rem(rem(j - 1, u2^2), u2) + 1

orientation = orientations[k]
left_element = neighbor_ids[1, k]
right_element = neighbor_ids[2, k]

u2 = size(u, 2)

@inbounds begin
interfaces_u[1, j1, j2, j3, k] = u[j1,
isequal(orientation, 1) * u2 + isequal(orientation, 2) * j2 + isequal(orientation, 3) * j2,
Expand Down Expand Up @@ -212,16 +223,17 @@ function prolong_boundaries_kernel!(boundaries_u, u, neighbor_ids, neighbor_side
k = (blockIdx().y - 1) * blockDim().y + threadIdx().y

if (j <= size(boundaries_u, 2) * size(boundaries_u, 3)^2 && k <= size(boundaries_u, 5))
j1 = div(j - 1, size(boundaries_u, 3)^2) + 1
j2 = div(rem(j - 1, size(boundaries_u, 3)^2), size(boundaries_u, 3)) + 1
j3 = rem(rem(j - 1, size(boundaries_u, 3)^2), size(boundaries_u, 3)) + 1
# size(boundaries_u, 3) == size(u, 2)
u2 = size(u, 2)

j1 = div(j - 1, u2^2) + 1
j2 = div(rem(j - 1, u2^2), u2) + 1
j3 = rem(rem(j - 1, u2^2), u2) + 1

element = neighbor_ids[k]
side = neighbor_sides[k]
orientation = orientations[k]

u2 = size(u, 2)

@inbounds begin
boundaries_u[1, j1, j2, j3, k] = u[j1,
isequal(orientation, 1) * u2 + isequal(orientation, 2) * j2 + isequal(orientation, 3) * j2,
Expand Down Expand Up @@ -312,8 +324,11 @@ function prolong_mortars_small2small_kernel!(u_upper_left, u_upper_right, u_lowe
k = (blockIdx().z - 1) * blockDim().z + threadIdx().z

if (i <= size(u_upper_left, 2) && j <= size(u_upper_left, 3)^2 && k <= size(u_upper_left, 5))
j1 = div(j - 1, size(u_upper_left, 3)) + 1
j2 = rem(j - 1, size(u_upper_left, 3)) + 1
# size(u_upper_left, 3) == size(u, 2)
u2 = size(u, 2)

j1 = div(j - 1, u2) + 1
j2 = rem(j - 1, u2) + 1

large_side = large_sides[k]
orientation = orientations[k]
Expand All @@ -323,8 +338,6 @@ function prolong_mortars_small2small_kernel!(u_upper_left, u_upper_right, u_lowe
upper_left_element = neighbor_ids[3, k]
upper_right_element = neighbor_ids[4, k]

u2 = size(u, 2)

@inbounds begin
u_upper_left[2, i, j1, j2, k] = u[i,
isequal(orientation, 1) + isequal(orientation, 2) * j1 + isequal(orientation, 3) * j1,
Expand Down Expand Up @@ -397,15 +410,17 @@ function prolong_mortars_large2small_kernel!(u_upper_left, u_upper_right, u_lowe
k = (blockIdx().z - 1) * blockDim().z + threadIdx().z

if (i <= size(u_upper_left, 2) && j <= size(u_upper_left, 3)^2 && k <= size(u_upper_left, 5))
j1 = div(j - 1, size(u_upper_left, 3)) + 1
j2 = rem(j - 1, size(u_upper_left, 3)) + 1
# size(u_upper_left, 3) == size(u, 2)
u2 = size(u, 2)

j1 = div(j - 1, u2) + 1
j2 = rem(j - 1, u2) + 1

large_side = large_sides[k]
orientation = orientations[k]
large_element = neighbor_ids[5, k]

leftright = large_side
u2 = size(u, 2)

@inbounds begin
for j1j1 in axes(forward_lower, 2)
Expand Down Expand Up @@ -655,20 +670,23 @@ function surface_integral_kernel!(du, factor_arr, surface_flux_values,
k = (blockIdx().z - 1) * blockDim().z + threadIdx().z

if (i <= size(du, 1) && j <= size(du, 2)^3 && k <= size(du, 5))
j1 = div(j - 1, size(du, 2)^2) + 1
j2 = div(rem(j - 1, size(du, 2)^2), size(du, 2)) + 1
j3 = rem(rem(j - 1, size(du, 2)^2), size(du, 2)) + 1
# size(du, 2) == size(u, 2)
u2 = size(du, 2)

j1 = div(j - 1, u2^2) + 1
j2 = div(rem(j - 1, u2^2), u2) + 1
j3 = rem(rem(j - 1, u2^2), u2) + 1

@inbounds begin
du[i, j1, j2, j3, k] -= (surface_flux_values[i, j2, j3, 1, k] * isequal(j1, 1) +
surface_flux_values[i, j1, j3, 3, k] * isequal(j2, 1) +
surface_flux_values[i, j1, j2, 5, k] * isequal(j3, 1)) *
factor_arr[1]
du[i, j1, j2, j3, k] += (surface_flux_values[i, j2, j3, 2, k] *
isequal(j1, size(du, 2)) +
isequal(j1, u2) +
surface_flux_values[i, j1, j3, 4, k] *
isequal(j2, size(du, 2)) +
surface_flux_values[i, j1, j2, 6, k] * isequal(j3, size(du, 2))) *
isequal(j2, u2) +
surface_flux_values[i, j1, j2, 6, k] * isequal(j3, u2)) *
factor_arr[2]
end
end
Expand All @@ -683,9 +701,12 @@ function jacobian_kernel!(du, inverse_jacobian, equations::AbstractEquations{3})
k = (blockIdx().z - 1) * blockDim().z + threadIdx().z

if (i <= size(du, 1) && j <= size(du, 2)^3 && k <= size(du, 5))
j1 = div(j - 1, size(du, 2)^2) + 1
j2 = div(rem(j - 1, size(du, 2)^2), size(du, 2)) + 1
j3 = rem(rem(j - 1, size(du, 2)^2), size(du, 2)) + 1
# size(du, 2) == size(u, 2)
u2 = size(du, 2)

j1 = div(j - 1, u2^2) + 1
j2 = div(rem(j - 1, u2^2), u2) + 1
j3 = rem(rem(j - 1, u2^2), u2) + 1

@inbounds du[i, j1, j2, j3, k] *= -inverse_jacobian[k]
end
Expand All @@ -700,9 +721,12 @@ function source_terms_kernel!(du, u, node_coordinates, t, equations::AbstractEqu
k = (blockIdx().y - 1) * blockDim().y + threadIdx().y

if (j <= size(du, 2)^3 && k <= size(du, 5))
j1 = div(j - 1, size(du, 2)^2) + 1
j2 = div(rem(j - 1, size(du, 2)^2), size(du, 2)) + 1
j3 = rem(rem(j - 1, size(du, 2)^2), size(du, 2)) + 1
# size(du, 2) == size(u, 2)
u2 = size(u, 2)

j1 = div(j - 1, u2^2) + 1
j2 = div(rem(j - 1, u2^2), u2) + 1
j3 = rem(rem(j - 1, u2^2), u2) + 1

u_local = get_node_vars(u, equations, j1, j2, j3, k)
x_local = get_node_coords(node_coordinates, equations, j1, j2, j3, k)
Expand Down

0 comments on commit eab1329

Please sign in to comment.