STLs appear to be working on the GPU with NVHPC!

Author: danieljvickers

src/common/m_derived_types.fpp

@@ -184,15 +184,15 @@ module m_derived_types
     end type t_model
 
     type :: t_model_array
-        ! --- Original CPU-side fields (unchanged) ---
+        ! Original CPU-side fields (unchanged)
         type(t_model), allocatable :: model
         real(wp), allocatable, dimension(:, :, :) :: boundary_v
         real(wp), allocatable, dimension(:, :) :: interpolated_boundary_v
         integer :: boundary_edge_count
         integer :: total_vertices
         integer :: interpolate
 
-        ! --- GPU-friendly flattened arrays ---
+        ! GPU-friendly flattened arrays
         integer :: ntrs  ! copy of model%ntrs
         real(wp), allocatable, dimension(:, :, :) :: trs_v  ! (3, 3, ntrs) - triangle vertices
         real(wp), allocatable, dimension(:, :) :: trs_n  ! (3, ntrs)    - triangle normals

src/common/m_model.fpp

@@ -17,22 +17,27 @@ module m_model
 
     private
 
-    public :: s_instantiate_STL_models, f_model_read, s_model_write, s_model_free, f_model_is_inside, models, gpu_ntrs, gpu_trs_v, gpu_trs_n
+    public :: s_instantiate_STL_models, f_model_read, s_model_write, s_model_free, f_model_is_inside, models, gpu_ntrs, &
+              gpu_trs_v, gpu_trs_n, gpu_boundary_v, gpu_interpolated_boundary_v, gpu_interpolate, gpu_boundary_edge_count, &
+              gpu_total_vertices
 
     ! Subroutines for STL immersed boundaries
     public :: f_check_boundary, f_register_edge, f_check_interpolation_2D, &
               f_check_interpolation_3D, f_interpolate_2D, f_interpolate_3D, &
               f_interpolated_distance, f_normals, f_distance, f_distance_normals_3D, f_tri_area, s_pack_model_for_gpu, &
-              f_model_is_inside_flat
+              f_model_is_inside_flat, f_distance_normals_3d_flat
 
     !! array of STL models that can be allocated and then used in IB marker and levelset compute
     type(t_model_array), allocatable, target :: models(:)
     !! GPU-friendly flat arrays for STL model data
     integer, allocatable :: gpu_ntrs(:)
     real(wp), allocatable, dimension(:, :, :, :) :: gpu_trs_v
     real(wp), allocatable, dimension(:, :, :) :: gpu_trs_n
-    real(wp), allocatable, dimension(:, :, :) :: gpu_boundary_v
+    real(wp), allocatable, dimension(:, :, :, :) :: gpu_boundary_v
     real(wp), allocatable, dimension(:, :, :) :: gpu_interpolated_boundary_v
+    integer, allocatable :: gpu_interpolate(:)
+    integer, allocatable :: gpu_boundary_edge_count(:)
+    integer, allocatable :: gpu_total_vertices(:)
 
 contains
 
@@ -174,34 +179,80 @@ contains
         ! Pack and upload flat arrays for GPU (AFTER the loop)
         block
             integer :: pid, max_ntrs
+            integer :: max_bv1, max_bv2, max_bv3, max_iv1, max_iv2
 
             max_ntrs = 0
+            max_bv1 = 0; max_bv2 = 0; max_bv3 = 0
+            max_iv1 = 0; max_iv2 = 0
+
             do pid = 1, num_ibs
                 if (allocated(models(pid)%model)) then
                     call s_pack_model_for_gpu(models(pid))
                     max_ntrs = max(max_ntrs, models(pid)%ntrs)
                 end if
+                if (allocated(models(pid)%boundary_v)) then
+                    max_bv1 = max(max_bv1, size(models(pid)%boundary_v, 1))
+                    max_bv2 = max(max_bv2, size(models(pid)%boundary_v, 2))
+                    max_bv3 = max(max_bv3, size(models(pid)%boundary_v, 3))
+                end if
+                if (allocated(models(pid)%interpolated_boundary_v)) then
+                    max_iv1 = max(max_iv1, size(models(pid)%interpolated_boundary_v, 1))
+                    max_iv2 = max(max_iv2, size(models(pid)%interpolated_boundary_v, 2))
+                end if
             end do
 
             if (max_ntrs > 0) then
                 allocate (gpu_ntrs(1:num_ibs))
                 allocate (gpu_trs_v(1:3, 1:3, 1:max_ntrs, 1:num_ibs))
                 allocate (gpu_trs_n(1:3, 1:max_ntrs, 1:num_ibs))
+                allocate (gpu_interpolate(1:num_ibs))
+                allocate (gpu_boundary_edge_count(1:num_ibs))
+                allocate (gpu_total_vertices(1:num_ibs))
 
                 gpu_ntrs = 0
                 gpu_trs_v = 0._wp
                 gpu_trs_n = 0._wp
+                gpu_interpolate = 0
+                gpu_boundary_edge_count = 0
+                gpu_total_vertices = 0
+
+                if (max_bv1 > 0) then
+                    allocate (gpu_boundary_v(1:max_bv1, 1:max_bv2, 1:max_bv3, 1:num_ibs))
+                    gpu_boundary_v = 0._wp
+                end if
+
+                if (max_iv1 > 0) then
+                    allocate (gpu_interpolated_boundary_v(1:max_iv1, 1:max_iv2, 1:num_ibs))
+                    gpu_interpolated_boundary_v = 0._wp
+                end if
 
                 do pid = 1, num_ibs
                     if (allocated(models(pid)%model)) then
                         gpu_ntrs(pid) = models(pid)%ntrs
                         gpu_trs_v(:, :, 1:models(pid)%ntrs, pid) = models(pid)%trs_v
                         gpu_trs_n(:, 1:models(pid)%ntrs, pid) = models(pid)%trs_n
+                        gpu_interpolate(pid) = models(pid)%interpolate
+                        gpu_boundary_edge_count(pid) = models(pid)%boundary_edge_count
+                        gpu_total_vertices(pid) = models(pid)%total_vertices
+                    end if
+                    if (allocated(models(pid)%boundary_v)) then
+                        gpu_boundary_v(1:size(models(pid)%boundary_v, 1), &
+                                       1:size(models(pid)%boundary_v, 2), &
+                                       1:size(models(pid)%boundary_v, 3), pid) = models(pid)%boundary_v
+                    end if
+                    if (allocated(models(pid)%interpolated_boundary_v)) then
+                        gpu_interpolated_boundary_v(1:size(models(pid)%interpolated_boundary_v, 1), &
+                                                    1:size(models(pid)%interpolated_boundary_v, 2), pid) = models(pid)%interpolated_boundary_v
                     end if
                 end do
 
-                $:GPU_ENTER_DATA(copyin='[gpu_ntrs, gpu_trs_v, gpu_trs_n]')
-
+                $:GPU_ENTER_DATA(copyin='[gpu_ntrs, gpu_trs_v, gpu_trs_n, gpu_interpolate, gpu_boundary_edge_count, gpu_total_vertices]')
+                if (allocated(gpu_boundary_v)) then
+                    $:GPU_ENTER_DATA(copyin='[gpu_boundary_v]')
+                end if
+                if (allocated(gpu_interpolated_boundary_v)) then
+                    $:GPU_ENTER_DATA(copyin='[gpu_interpolated_boundary_v]')
+                end if
             end if
         end block
 
@@ -731,7 +782,7 @@ contains
                 end if
             end do
 
-            ! if thje ray hits an odd number of triangles on its way out, then
+            ! if the ray hits an odd number of triangles on its way out, then
             ! it must be on the inside of the model
             nInOrOut = nInOrOut + mod(nHits, 2)
         end do
@@ -1371,6 +1422,66 @@ contains
 
     end subroutine f_distance_normals_3D
 
+    subroutine f_distance_normals_3D_flat(ntrs, trs_v, trs_n, pid, point, normals, distance)
+
+        $:GPU_ROUTINE(parallelism='[seq]')
+
+        integer, intent(in) :: ntrs
+        real(wp), dimension(:, :, :, :), intent(in) :: trs_v
+        real(wp), dimension(:, :, :), intent(in) :: trs_n
+        integer, intent(in) :: pid
+        real(wp), dimension(1:3), intent(in) :: point
+        real(wp), dimension(1:3), intent(out) :: normals
+        real(wp), intent(out) :: distance
+
+        real(wp), dimension(1:3, 1:3) :: tri
+        real(wp) :: dist_min, dist_t_min
+        real(wp) :: dist_min_normal, dist_buffer_normal
+        real(wp), dimension(1:3) :: midp
+        real(wp), dimension(1:3) :: dist_buffer
+        integer :: i, j, tri_idx
+
+        dist_min = 1.e12_wp
+        dist_min_normal = 1.e12_wp
+        distance = 0._wp
+
+        tri_idx = 0
+        do i = 1, ntrs
+            do j = 1, 3
+                tri(j, 1) = trs_v(j, 1, i, pid)
+                tri(j, 2) = trs_v(j, 2, i, pid)
+                tri(j, 3) = trs_v(j, 3, i, pid)
+                dist_buffer(j) = sqrt((point(1) - tri(j, 1))**2 + &
+                                      (point(2) - tri(j, 2))**2 + &
+                                      (point(3) - tri(j, 3))**2)
+            end do
+
+            do j = 1, 3
+                midp(j) = (tri(1, j) + tri(2, j) + tri(3, j))/3
+            end do
+
+            dist_t_min = minval(dist_buffer(1:3))
+            dist_buffer_normal = sqrt((point(1) - midp(1))**2 + &
+                                      (point(2) - midp(2))**2 + &
+                                      (point(3) - midp(3))**2)
+
+            if (dist_t_min < dist_min) then
+                dist_min = dist_t_min
+            end if
+
+            if (dist_buffer_normal < dist_min_normal) then
+                dist_min_normal = dist_buffer_normal
+                tri_idx = i
+            end if
+        end do
+
+        normals(1) = trs_n(1, tri_idx, pid)
+        normals(2) = trs_n(2, tri_idx, pid)
+        normals(3) = trs_n(3, tri_idx, pid)
+        distance = dist_min
+
+    end subroutine f_distance_normals_3D_flat
+
     !> This procedure determines the levelset distance of 2D models without interpolation.
     !! @param boundary_v                   Group of all the boundary vertices of the 2D model without interpolation
     !! @param boundary_vertex_count        Output the total number of boundary vertices

src/simulation/m_compute_levelset.fpp

@@ -50,6 +50,8 @@ contains
                     call s_cylinder_levelset(gps(i))
                 elseif (patch_geometry == 11) then
                     call s_3d_airfoil_levelset(gps(i))
+                elseif (patch_geometry == 12) then
+                    call s_model_levelset(gps(i))
                 end if
             end do
             $:END_GPU_PARALLEL_LOOP()
@@ -70,6 +72,8 @@ contains
                     call s_rectangle_levelset(gps(i))
                 elseif (patch_geometry == 4) then
                     call s_airfoil_levelset(gps(i))
+                elseif (patch_geometry == 5) then
+                    call s_model_levelset(gps(i))
                 elseif (patch_geometry == 6) then
                     call s_ellipse_levelset(gps(i))
                 end if
@@ -79,17 +83,6 @@ contains
 
         end if
 
-        ! STL models computed on the CPU for now
-        do i = 1, num_gps
-            patch_id = gps(i)%ib_patch_id
-            patch_geometry = patch_ib(patch_id)%geometry
-
-            if (patch_geometry == 5 .or. patch_geometry == 12) then
-                call s_model_levelset(gps(i))
-                $:GPU_UPDATE(device='[gps(i)]')
-            end if
-        end do
-
     end subroutine s_apply_levelset
 
     subroutine s_circle_levelset(gp)
@@ -651,7 +644,7 @@ contains
         type(ghost_point), intent(inout) :: gp
 
         integer :: i, j, k, patch_id, boundary_edge_count, total_vertices
-        logical :: interpolate
+        integer :: interpolate
         real(wp), dimension(1:3) :: center, xyz_local
         real(wp) :: normals(1:3) !< Boundary normal buffer
         real(wp) :: distance
@@ -663,16 +656,17 @@ contains
         k = gp%loc(3)
 
         ! load in model values
-        interpolate = models(patch_id)%interpolate
-        boundary_edge_count = models(patch_id)%boundary_edge_count
-        total_vertices = models(patch_id)%total_vertices
+        interpolate = gpu_interpolate(patch_id)
+        boundary_edge_count = gpu_boundary_edge_count(patch_id)
+        total_vertices = gpu_total_vertices(patch_id)
 
         center = 0._wp
         if (.not. f_is_default(patch_ib(patch_id)%x_centroid)) center(1) = patch_ib(patch_id)%x_centroid
         if (.not. f_is_default(patch_ib(patch_id)%y_centroid)) center(2) = patch_ib(patch_id)%y_centroid
         if (p > 0) then
             if (.not. f_is_default(patch_ib(patch_id)%z_centroid)) center(3) = patch_ib(patch_id)%z_centroid
         end if
+
         inverse_rotation(:, :) = patch_ib(patch_id)%rotation_matrix_inverse(:, :)
         rotation(:, :) = patch_ib(patch_id)%rotation_matrix(:, :)
 
@@ -691,11 +685,11 @@ contains
         if (p > 0) then
 
             ! Get the boundary normals and shortest distance between the cell center and the model boundary
-            call f_distance_normals_3D(models(patch_id)%model, xyz_local, normals, distance)
+            call f_distance_normals_3D_flat(gpu_ntrs(patch_id), gpu_trs_v, gpu_trs_n, patch_id, xyz_local, normals, distance)
 
             ! Get the shortest distance between the cell center and the interpolated model boundary
             if (interpolate) then
-                gp%levelset = f_interpolated_distance(models(patch_id)%interpolated_boundary_v, total_vertices, xyz_local)
+                gp%levelset = f_interpolated_distance(gpu_interpolated_boundary_v(:, :, patch_id), total_vertices, xyz_local)
             else
                 gp%levelset = distance
             end if
@@ -707,19 +701,18 @@ contains
             gp%levelset_norm = matmul(rotation, normals(1:3))
         else
             ! 2D models
-            if (interpolate) then
+            if (interpolate == 1) then
                 ! Get the shortest distance between the cell center and the model boundary
-                gp%levelset = f_interpolated_distance(models(patch_id)%interpolated_boundary_v, total_vertices, xyz_local)
+                gp%levelset = f_interpolated_distance(gpu_interpolated_boundary_v(:, :, patch_id), total_vertices, xyz_local)  ! Change 7
             else
-                ! Get the shortest distance between the cell center and the interpolated model boundary
-                gp%levelset = f_distance(models(patch_id)%boundary_v, boundary_edge_count, xyz_local)
+                gp%levelset = f_distance(gpu_boundary_v(:, :, :, patch_id), boundary_edge_count, xyz_local)  ! Change 8
             end if
 
             ! Correct the sign of the levelset
             gp%levelset = -abs(gp%levelset)
 
             ! Get the boundary normals
-            call f_normals(models(patch_id)%boundary_v, &
+            call f_normals(gpu_boundary_v(:, :, :, patch_id), &
                            boundary_edge_count, &
                            xyz_local, &
                            normals)
@@ -729,8 +722,6 @@ contains
 
         end if
 
-        ! print gp%levelset, gp%levelset_norm(1), gp%levelset_norm(2), gp%levelset_norm(3)
-
     end subroutine s_model_levelset
 
 end module m_compute_levelset

src/simulation/m_ib_patches.fpp

@@ -877,8 +877,6 @@ contains
         end do
         $:END_GPU_PARALLEL_LOOP()
 
-        $:GPU_UPDATE(host='[ib_markers%sf]')
-
     end subroutine s_ib_model
 
     !> The STL patch is a 3D geometry that is imported from an STL file.
@@ -935,8 +933,6 @@ contains
             end do
         end do
 
-        $:GPU_UPDATE(host='[ib_markers%sf]')
-
     end subroutine s_ib_3d_model
 
     !> Subroutine that computes a rotation matrix for converting to the rotating frame of the boundary