Skip to content

Commit f0085c9

Browse files
ghost points are now computed on the GPU
1 parent 85889c0 commit f0085c9

4 files changed

Lines changed: 155 additions & 176 deletions

File tree

src/common/include/parallel_macros.fpp

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -174,6 +174,16 @@
174174
#endif
175175
#:enddef
176176

177+
#:def END_GPU_ATOMIC_CAPTURE()
178+
#:set acc_end_directive = '!$acc end atomic'
179+
#:set omp_end_directive = '!$omp end atomic'
180+
#if defined(MFC_OpenACC)
181+
$:acc_end_directive
182+
#elif defined(MFC_OpenMP)
183+
$:omp_end_directive
184+
#endif
185+
#:enddef
186+
177187
#:def GPU_UPDATE(host=None, device=None, extraAccArgs=None, extraOmpArgs=None)
178188
#:set acc_code = ACC_UPDATE(host=host, device=device, extraAccArgs=extraAccArgs)
179189
#:set omp_code = OMP_UPDATE(host=host, device=device, extraOmpArgs=extraOmpArgs)

src/simulation/m_compute_levelset.fpp

Lines changed: 0 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -33,8 +33,6 @@ contains
3333

3434
integer :: i, patch_id, patch_geometry
3535

36-
$:GPU_UPDATE(device='[gps(1:num_gps)]')
37-
3836
! 3D Patch Geometries
3937
if (p > 0) then
4038

@@ -92,8 +90,6 @@ contains
9290
end if
9391
end do
9492

95-
$:GPU_UPDATE(host='[gps(1:num_gps)]')
96-
9793
end subroutine s_apply_levelset
9894

9995
subroutine s_circle_levelset(gp)

src/simulation/m_ib_patches.fpp

Lines changed: 48 additions & 47 deletions
Original file line numberDiff line numberDiff line change
@@ -274,8 +274,8 @@ contains
274274
275275
$:GPU_PARALLEL_LOOP(private='[i,j]',&
276276
& copyin='[patch_id,center,radius]', collapse=2)
277-
do j = -gp_layers, n+gp_layers
278-
do i = -gp_layers, m+gp_layers
277+
do j = -gp_layers, n + gp_layers
278+
do i = -gp_layers, m + gp_layers
279279
if ((x_cc(i) - center(1))**2 &
280280
+ (y_cc(j) - center(2))**2 <= radius**2) &
281281
then
@@ -378,8 +378,8 @@ contains
378378
379379
$:GPU_PARALLEL_LOOP(private='[i,j,xy_local,k,f]', &
380380
& copyin='[patch_id,center,inverse_rotation,offset,ma,ca_in,airfoil_grid_u,airfoil_grid_l]', collapse=2)
381-
do j = -gp_layers, n+gp_layers
382-
do i = -gp_layers, m+gp_layers
381+
do j = -gp_layers, n + gp_layers
382+
do i = -gp_layers, m + gp_layers
383383
xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp] ! get coordinate frame centered on IB
384384
xy_local = matmul(inverse_rotation, xy_local) ! rotate the frame into the IB's coordinates
385385
xy_local = xy_local - offset ! airfoils are a patch that require a centroid offset
@@ -529,9 +529,9 @@ contains
529529

530530
$:GPU_PARALLEL_LOOP(private='[i,j,l,xyz_local,k,f]',&
531531
& copyin='[patch_id,center,inverse_rotation,offset,ma,ca_in,airfoil_grid_u,airfoil_grid_l,z_min,z_max]', collapse=3)
532-
do l = -gp_layers, p+gp_layers
533-
do j = -gp_layers, n+gp_layers
534-
do i = -gp_layers, m+gp_layers
532+
do l = -gp_layers, p + gp_layers
533+
do j = -gp_layers, n + gp_layers
534+
do i = -gp_layers, m + gp_layers
535535
xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), z_cc(l) - center(3)] ! get coordinate frame centered on IB
536536
xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinates
537537
xyz_local = xyz_local - offset ! airfoils are a patch that require a centroid offset
@@ -623,8 +623,8 @@ contains
623623
! variables of the current patch are assigned to this cell.
624624
$:GPU_PARALLEL_LOOP(private='[i,j, xy_local]',&
625625
& copyin='[patch_id,center,length,inverse_rotation,x_cc,y_cc]', collapse=2)
626-
do j = -gp_layers, n+gp_layers
627-
do i = -gp_layers, m+gp_layers
626+
do j = -gp_layers, n + gp_layers
627+
do i = -gp_layers, m + gp_layers
628628
! get the x and y coordinates in the local IB frame
629629
xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp]
630630
xy_local = matmul(inverse_rotation, xy_local)
@@ -673,6 +673,7 @@ contains
673673
center(3) = patch_ib(patch_id)%z_centroid
674674
radius = patch_ib(patch_id)%radius
675675
676+
! find the indices to the left and right of the IB in i, j, k
676677
il = -gp_layers
677678
jl = -gp_layers
678679
kl = -gp_layers
@@ -689,10 +690,10 @@ contains
689690
! the current patch are assigned to this cell.
690691
$:GPU_PARALLEL_LOOP(private='[i,j,k,cart_y,cart_z]',&
691692
& copyin='[patch_id,center,radius]', collapse=3)
692-
do k = kl-1, kr+1
693-
do j = jl-1, jr+1
694-
do i = il-1, ir+1
695-
! do i = -gp_layers, m+gp_layers
693+
do k = kl - 1, kr + 1
694+
do j = jl - 1, jr + 1
695+
do i = il - 1, ir + 1
696+
! do i = -gp_layers, m+gp_layers
696697
if (grid_geometry == 3) then
697698
call s_convert_cylindrical_to_cartesian_coord(y_cc(j), z_cc(k))
698699
else
@@ -746,9 +747,9 @@ contains
746747
! of the current patch are assigned to this cell.
747748
$:GPU_PARALLEL_LOOP(private='[i,j,k,xyz_local,cart_y,cart_z]',&
748749
& copyin='[patch_id,center,length,inverse_rotation]', collapse=3)
749-
do k = -gp_layers, p+gp_layers
750-
do j = -gp_layers, n+gp_layers
751-
do i = -gp_layers, m+gp_layers
750+
do k = -gp_layers, p + gp_layers
751+
do j = -gp_layers, n + gp_layers
752+
do i = -gp_layers, m + gp_layers
752753

753754
if (grid_geometry == 3) then
754755
! TODO :: This does not work and is not covered by any tests. This should be fixed
@@ -816,9 +817,9 @@ contains
816817
! variables of the current patch are assigned to this cell.
817818
$:GPU_PARALLEL_LOOP(private='[i,j,k,xyz_local,cart_y,cart_z]',&
818819
& copyin='[patch_id,center,length,radius,inverse_rotation]', collapse=3)
819-
do k = -gp_layers, p+gp_layers
820-
do j = -gp_layers, n+gp_layers
821-
do i = -gp_layers, m+gp_layers
820+
do k = -gp_layers, p + gp_layers
821+
do j = -gp_layers, n + gp_layers
822+
do i = -gp_layers, m + gp_layers
822823
823824
if (grid_geometry == 3) then
824825
call s_convert_cylindrical_to_cartesian_coord(y_cc(j), z_cc(k))
@@ -879,8 +880,8 @@ contains
879880
! domain
880881
$:GPU_PARALLEL_LOOP(private='[i,j, xy_local]',&
881882
& copyin='[patch_id,center,ellipse_coeffs,inverse_rotation,x_cc,y_cc]', collapse=2)
882-
do j = -gp_layers, n+gp_layers
883-
do i = -gp_layers, m+gp_layers
883+
do j = -gp_layers, n + gp_layers
884+
do i = -gp_layers, m + gp_layers
884885
! get the x and y coordinates in the local IB frame
885886
xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp]
886887
xy_local = matmul(inverse_rotation, xy_local)
@@ -920,23 +921,23 @@ contains
920921
inverse_rotation(:, :) = patch_ib(patch_id)%rotation_matrix_inverse(:, :)
921922
offset(:) = patch_ib(patch_id)%centroid_offset(:)
922923
923-
do i = 0-gp_layers, m+gp_layers
924-
do j = -gp_layers, n+gp_layers
924+
do i = 0 - gp_layers, m + gp_layers
925+
do j = -gp_layers, n + gp_layers
925926
926-
xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp]
927-
xy_local = matmul(inverse_rotation, xy_local)
928-
xy_local = xy_local - offset
927+
xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp]
928+
xy_local = matmul(inverse_rotation, xy_local)
929+
xy_local = xy_local - offset
929930
930-
if (grid_geometry == 3) then
931-
xy_local = f_convert_cyl_to_cart(xy_local)
932-
end if
931+
if (grid_geometry == 3) then
932+
xy_local = f_convert_cyl_to_cart(xy_local)
933+
end if
933934
934-
eta = f_model_is_inside(model, xy_local, (/dx(i), dy(j), 0._wp/), patch_ib(patch_id)%model_spc)
935+
eta = f_model_is_inside(model, xy_local, (/dx(i), dy(j), 0._wp/), patch_ib(patch_id)%model_spc)
935936
936-
! Reading STL boundary vertices and compute the levelset and levelset_norm
937-
if (eta > patch_ib(patch_id)%model_threshold) then
938-
ib_markers%sf(i, j, 0) = patch_id
939-
end if
937+
! Reading STL boundary vertices and compute the levelset and levelset_norm
938+
if (eta > patch_ib(patch_id)%model_threshold) then
939+
ib_markers%sf(i, j, 0) = patch_id
940+
end if
940941
941942
end do
942943
end do
@@ -968,9 +969,9 @@ contains
968969
inverse_rotation(:, :) = patch_ib(patch_id)%rotation_matrix_inverse(:, :)
969970
offset(:) = patch_ib(patch_id)%centroid_offset(:)
970971
971-
do i = 0-gp_layers, m+gp_layers
972-
do j = -gp_layers, n+gp_layers
973-
do k = -gp_layers, p+gp_layers
972+
do i = 0 - gp_layers, m + gp_layers
973+
do j = -gp_layers, n + gp_layers
974+
do k = -gp_layers, p + gp_layers
974975
975976
xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), z_cc(k) - center(3)]
976977
xyz_local = matmul(inverse_rotation, xyz_local)
@@ -1072,21 +1073,21 @@ contains
10721073
10731074
subroutine get_bounding_indices(left_bound, right_bound, cell_centers, left_index, right_index)
10741075
1075-
real(wp), intent(in) :: left_bound, right_bound
1076-
integer, intent(inout) :: left_index, right_index
1077-
real(wp), dimension(-buff_size:), intent(in) :: cell_centers
1076+
real(wp), intent(in) :: left_bound, right_bound
1077+
integer, intent(inout) :: left_index, right_index
1078+
real(wp), dimension(-buff_size:), intent(in) :: cell_centers
10781079
1079-
integer :: itr_left, itr_middle, itr_right
1080+
integer :: itr_left, itr_middle, itr_right
10801081
10811082
itr_left = left_index
10821083
itr_right = right_index
10831084
10841085
do while (itr_left + 1 < itr_right)
1085-
itr_middle = (itr_left + itr_right) / 2
1086+
itr_middle = (itr_left + itr_right)/2
10861087
if (cell_centers(itr_middle) < left_bound) then
1087-
itr_left = itr_middle
1088+
itr_left = itr_middle
10881089
else if (cell_centers(itr_middle) > left_bound) then
1089-
itr_right = itr_middle
1090+
itr_right = itr_middle
10901091
else
10911092
itr_left = itr_middle
10921093
exit
@@ -1096,11 +1097,11 @@ contains
10961097
10971098
itr_right = right_index
10981099
do while (itr_left + 1 < itr_right)
1099-
itr_middle = (itr_left + itr_right) / 2
1100+
itr_middle = (itr_left + itr_right)/2
11001101
if (cell_centers(itr_middle) < right_bound) then
1101-
itr_left = itr_middle
1102+
itr_left = itr_middle
11021103
else if (cell_centers(itr_middle) > right_bound) then
1103-
itr_right = itr_middle
1104+
itr_right = itr_middle
11041105
else
11051106
itr_right = itr_middle
11061107
exit

0 commit comments

Comments
 (0)