patch_ib Memory Copy Reduction

src/simulation/m_data_output.fpp

@@ -889,6 +889,8 @@ contains
 
         integer, intent(in) :: time_step
 
+        $:GPU_UPDATE(host='[patch_ib(1:num_ibs)]')
+
         if (parallel_io) then
             call s_write_parallel_ib_data(time_step)
         else

src/simulation/m_ib_patches.fpp

@@ -456,9 +456,6 @@ contains
 
         center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg)
         center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg)
-        length(1) = patch_ib(patch_id)%length_x
-        length(2) = patch_ib(patch_id)%length_y
-        inverse_rotation(:,:) = patch_ib(patch_id)%rotation_matrix_inverse(:,:)
 
         ! encode the periodicity information into the patch_id
         call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, 0, encoded_patch_id)
@@ -468,21 +465,22 @@ contains
         jl = -gp_layers - 1
         ir = m + gp_layers + 1
         jr = n + gp_layers + 1
-        corner_distance = sqrt(dot_product(length, length))/2._wp  ! maximum distance any marker can be from the center
+        ! maximum distance any marker can be from the center
+        corner_distance = 0.5_wp*sqrt(patch_ib(patch_id)%length_x**2 + patch_ib(patch_id)%length_y**2)
         call get_bounding_indices(center(1) - corner_distance, center(1) + corner_distance, x_cc, il, ir)
         call get_bounding_indices(center(2) - corner_distance, center(2) + corner_distance, y_cc, jl, jr)
 
         ! Assign primitive variables if rectangle covers cell and patch has write permission
-        $:GPU_PARALLEL_LOOP(private='[i, j, xy_local]', copyin='[encoded_patch_id, center, length, inverse_rotation, x_cc, &
-                            & y_cc]', collapse=2)
+        $:GPU_PARALLEL_LOOP(private='[i, j, xy_local]', copyin='[encoded_patch_id, center]', collapse=2)
         do j = jl, jr
             do i = il, ir
                 ! get the x and y coordinates in the local IB frame
                 xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp]
-                xy_local = matmul(inverse_rotation, xy_local)
+                xy_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse, xy_local)
 
-                if (-0.5_wp*length(1) <= xy_local(1) .and. 0.5_wp*length(1) >= xy_local(1) .and. -0.5_wp*length(2) <= xy_local(2) &
-                    & .and. 0.5_wp*length(2) >= xy_local(2)) then
+                if (-0.5_wp*patch_ib(patch_id)%length_x <= xy_local(1) .and. 0.5_wp*patch_ib(patch_id)%length_x >= xy_local(1) &
+                    & .and. -0.5_wp*patch_ib(patch_id)%length_y <= xy_local(2) &
+                    & .and. 0.5_wp*patch_ib(patch_id)%length_y >= xy_local(2)) then
                     ! Updating the patch identities bookkeeping variable
                     ib_markers%sf(i, j, 0) = encoded_patch_id
                 end if
@@ -556,19 +554,14 @@ contains
         integer, intent(in)                :: xp, yp, zp  !< integers containing the periodicity projection information
         integer                            :: i, j, k, ir, il, jr, jl, kr, kl  !< Generic loop iterators
         integer                            :: encoded_patch_id
-        real(wp), dimension(1:3)           :: xyz_local, center, length  !< x and y coordinates in local IB frame
-        real(wp), dimension(1:3,1:3)       :: inverse_rotation
+        real(wp), dimension(1:3)           :: xyz_local, center  !< x and y coordinates in local IB frame
         real(wp)                           :: corner_distance
 
-        ! Transferring the cuboid's centroid and length information
+        ! Transferring the cuboid's centroid
 
         center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg)
         center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg)
         center(3) = patch_ib(patch_id)%z_centroid + real(zp, wp)*(z_domain%end - z_domain%beg)
-        length(1) = patch_ib(patch_id)%length_x
-        length(2) = patch_ib(patch_id)%length_y
-        length(3) = patch_ib(patch_id)%length_z
-        inverse_rotation(:,:) = patch_ib(patch_id)%rotation_matrix_inverse(:,:)
 
         ! encode the periodicity information into the patch_id
         call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, zp, encoded_patch_id)
@@ -580,25 +573,29 @@ contains
         ir = m + gp_layers + 1
         jr = n + gp_layers + 1
         kr = p + gp_layers + 1
-        corner_distance = sqrt(dot_product(length, length))/2._wp  ! maximum distance any marker can be from the center
+        corner_distance = 0.5_wp*sqrt(patch_ib(patch_id)%length_x**2 + patch_ib(patch_id)%length_y**2 &
+                                      & + patch_ib(patch_id)%length_z**2)  ! maximum distance any marker can be from the center
         call get_bounding_indices(center(1) - corner_distance, center(1) + corner_distance, x_cc, il, ir)
         call get_bounding_indices(center(2) - corner_distance, center(2) + corner_distance, y_cc, jl, jr)
         call get_bounding_indices(center(3) - corner_distance, center(3) + corner_distance, z_cc, kl, kr)
 
         ! Checking whether the cuboid covers a particular cell in the domain and verifying whether the current patch has permission
         ! to write to to that cell. If both queries check out, the primitive variables of the current patch are assigned to this
         ! cell.
-        $:GPU_PARALLEL_LOOP(private='[i, j, k, xyz_local]', copyin='[encoded_patch_id, center, length, inverse_rotation]', &
-                            & collapse=3)
+        $:GPU_PARALLEL_LOOP(private='[i, j, k, xyz_local]', copyin='[encoded_patch_id, center]', collapse=3)
         do k = kl, kr
             do j = jl, jr
                 do i = il, ir
                     xyz_local = [x_cc(i), y_cc(j), z_cc(k)] - center  ! get coordinate frame centered on IB
-                    xyz_local = matmul(inverse_rotation, xyz_local)  ! rotate the frame into the IB's coordinates
-
-                    if (-0.5*length(1) <= xyz_local(1) .and. 0.5*length(1) >= xyz_local(1) .and. -0.5*length(2) <= xyz_local(2) &
-                        & .and. 0.5*length(2) >= xyz_local(2) .and. -0.5*length(3) <= xyz_local(3) .and. 0.5*length(3) &
-                        & >= xyz_local(3)) then
+                    ! rotate the frame into the IB's coordinates
+                    xyz_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse, xyz_local)
+
+                    if (-0.5_wp*patch_ib(patch_id)%length_x <= xyz_local(1) &
+                        & .and. 0.5_wp*patch_ib(patch_id)%length_x >= xyz_local(1) .and. &
+                        & -0.5_wp*patch_ib(patch_id)%length_y <= xyz_local(2) &
+                        & .and. 0.5_wp*patch_ib(patch_id)%length_y >= xyz_local(2) .and. &
+                        & -0.5_wp*patch_ib(patch_id)%length_z <= xyz_local(3) &
+                        & .and. 0.5_wp*patch_ib(patch_id)%length_z >= xyz_local(3)) then
                         ! Updating the patch identities bookkeeping variable
                         ib_markers%sf(i, j, k) = encoded_patch_id
                     end if
@@ -617,21 +614,16 @@ contains
         integer, intent(in)                :: xp, yp, zp  !< integers containing the periodicity projection information
         integer                            :: i, j, k, il, ir, jl, jr, kl, kr  !< Generic loop iterators
         integer                            :: encoded_patch_id
-        real(wp)                           :: radius
         real(wp), dimension(1:3)           :: xyz_local, center, length  !< x and y coordinates in local IB frame
         real(wp), dimension(1:3,1:3)       :: inverse_rotation
         real(wp)                           :: corner_distance
 
-        ! Transferring the cylindrical patch's centroid, length, radius,
+        ! Transferring the cylindrical patch's centroid
 
         center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg)
         center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg)
         center(3) = patch_ib(patch_id)%z_centroid + real(zp, wp)*(z_domain%end - z_domain%beg)
-        length(1) = patch_ib(patch_id)%length_x
-        length(2) = patch_ib(patch_id)%length_y
-        length(3) = patch_ib(patch_id)%length_z
-        radius = patch_ib(patch_id)%radius
-        inverse_rotation(:,:) = patch_ib(patch_id)%rotation_matrix_inverse(:,:)
+        length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, patch_ib(patch_id)%length_z]
 
         ! encode the periodicity information into the patch_id
         call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, zp, encoded_patch_id)
@@ -642,28 +634,31 @@ contains
         ir = m + gp_layers + 1
         jr = n + gp_layers + 1
         kr = p + gp_layers + 1
-        corner_distance = sqrt(radius**2 + maxval(length)**2)  ! distance to rim of cylinder
+        corner_distance = sqrt(patch_ib(patch_id)%radius**2 + maxval(length)**2)  ! distance to rim of cylinder
         call get_bounding_indices(center(1) - corner_distance, center(1) + corner_distance, x_cc, il, ir)
         call get_bounding_indices(center(2) - corner_distance, center(2) + corner_distance, y_cc, jl, jr)
         call get_bounding_indices(center(3) - corner_distance, center(3) + corner_distance, z_cc, kl, kr)
 
         ! Checking whether the cylinder covers a particular cell in the domain and verifying whether the current patch has the
         ! permission to write to that cell. If both queries check out, the primitive variables of the current patch are assigned to
         ! this cell.
-        $:GPU_PARALLEL_LOOP(private='[i, j, k, xyz_local]', copyin='[encoded_patch_id, center, length, radius, &
-                            & inverse_rotation]', collapse=3)
+        $:GPU_PARALLEL_LOOP(private='[i, j, k, xyz_local]', copyin='[encoded_patch_id, center]', collapse=3)
         do k = kl, kr
             do j = jl, jr
                 do i = il, ir
                     xyz_local = [x_cc(i), y_cc(j), z_cc(k)] - center  ! get coordinate frame centered on IB
-                    xyz_local = matmul(inverse_rotation, xyz_local)  ! rotate the frame into the IB's coordinates
-
-                    if (((.not. f_is_default(length(1)) .and. xyz_local(2)**2 + xyz_local(3)**2 <= radius**2 .and. &
-                        & -0.5_wp*length(1) <= xyz_local(1) .and. 0.5_wp*length(1) >= xyz_local(1)) &
-                        & .or. (.not. f_is_default(length(2)) .and. xyz_local(1)**2 + xyz_local(3)**2 <= radius**2 .and. &
-                        & -0.5_wp*length(2) <= xyz_local(2) .and. 0.5_wp*length(2) >= xyz_local(2)) &
-                        & .or. (.not. f_is_default(length(3)) .and. xyz_local(1)**2 + xyz_local(2)**2 <= radius**2 .and. &
-                        & -0.5_wp*length(3) <= xyz_local(3) .and. 0.5_wp*length(3) >= xyz_local(3)))) then
+                    ! rotate the frame into the IB's coordinates
+                    xyz_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse, xyz_local)
+
+                    if (((.not. f_is_default(patch_ib(patch_id)%length_x) .and. xyz_local(2)**2 + xyz_local(3) &
+                        & **2 <= patch_ib(patch_id)%radius**2 .and. -0.5_wp*patch_ib(patch_id)%length_x <= xyz_local(1) &
+                        & .and. 0.5_wp*patch_ib(patch_id)%length_x >= xyz_local(1)) &
+                        & .or. (.not. f_is_default(patch_ib(patch_id)%length_y) .and. xyz_local(1)**2 + xyz_local(3) &
+                        & **2 <= patch_ib(patch_id)%radius**2 .and. -0.5_wp*patch_ib(patch_id)%length_y <= xyz_local(2) &
+                        & .and. 0.5_wp*patch_ib(patch_id)%length_y >= xyz_local(2)) &
+                        & .or. (.not. f_is_default(patch_ib(patch_id)%length_z) .and. xyz_local(1)**2 + xyz_local(2) &
+                        & **2 <= patch_ib(patch_id)%radius**2 .and. -0.5_wp*patch_ib(patch_id)%length_z <= xyz_local(3) &
+                        & .and. 0.5_wp*patch_ib(patch_id)%length_z >= xyz_local(3)))) then
                         ! Updating the patch identities bookkeeping variable
                         ib_markers%sf(i, j, k) = encoded_patch_id
                     end if
@@ -683,40 +678,37 @@ contains
         integer                            :: i, j, il, ir, jl, jr  !< Generic loop iterators
         integer                            :: encoded_patch_id
         real(wp), dimension(1:3)           :: xy_local              !< x and y coordinates in local IB frame
-        real(wp), dimension(1:2)           :: ellipse_coeffs        !< a and b in the ellipse coefficients
         real(wp), dimension(1:2)           :: center                !< x and y coordinates in local IB frame
-        real(wp), dimension(1:3,1:3)       :: inverse_rotation
+        real(wp)                           :: bounding_radius
 
-        ! Transferring the ellipse's centroid and length information
+        ! Transferring the ellipse's centroid
 
         center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg)
         center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg)
-        ellipse_coeffs(1) = 0.5_wp*patch_ib(patch_id)%length_x
-        ellipse_coeffs(2) = 0.5_wp*patch_ib(patch_id)%length_y
-        inverse_rotation(:,:) = patch_ib(patch_id)%rotation_matrix_inverse(:,:)
 
         ! encode the periodicity information into the patch_id
         call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, 0, encoded_patch_id)
 
         ! find the indices to the left and right of the IB in i, j, k
+        bounding_radius = 0.5_wp*max(patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y)
         il = -gp_layers - 1
         jl = -gp_layers - 1
         ir = m + gp_layers + 1
         jr = n + gp_layers + 1
-        call get_bounding_indices(center(1) - maxval(ellipse_coeffs)*2._wp, center(1) + maxval(ellipse_coeffs)*2._wp, x_cc, il, ir)
-        call get_bounding_indices(center(2) - maxval(ellipse_coeffs)*2._wp, center(2) + maxval(ellipse_coeffs)*2._wp, y_cc, jl, jr)
+        call get_bounding_indices(center(1) - bounding_radius*2._wp, center(1) + bounding_radius*2._wp, x_cc, il, ir)
+        call get_bounding_indices(center(2) - bounding_radius*2._wp, center(2) + bounding_radius*2._wp, y_cc, jl, jr)
 
         ! Checking whether the ellipse covers a particular cell in the domain
-        $:GPU_PARALLEL_LOOP(private='[i, j, xy_local]', copyin='[encoded_patch_id, center, ellipse_coeffs, inverse_rotation, &
-                            & x_cc, y_cc]', collapse=2)
+        $:GPU_PARALLEL_LOOP(private='[i, j, xy_local]', copyin='[encoded_patch_id, center]', collapse=2)
         do j = jl, jr
             do i = il, ir
                 ! get the x and y coordinates in the local IB frame
                 xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp]
-                xy_local = matmul(inverse_rotation, xy_local)
+                xy_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse(:,:), xy_local)
 
                 ! Ellipse condition (x/a)^2 + (y/b)^2 <= 1
-                if ((xy_local(1)/ellipse_coeffs(1))**2 + (xy_local(2)/ellipse_coeffs(2))**2 <= 1._wp) then
+                if ((xy_local(1)/(0.5_wp*patch_ib(patch_id)%length_x))**2 + (xy_local(2)/(0.5_wp*patch_ib(patch_id)%length_y)) &
+                    & **2 <= 1._wp) then
                     ! Updating the patch identities bookkeeping variable
                     ib_markers%sf(i, j, 0) = encoded_patch_id
                 end if
@@ -895,6 +887,8 @@ contains
     !> Compute a rotation matrix for converting to the rotating frame of the boundary
     subroutine s_update_ib_rotation_matrix(patch_id)
 
+        $:GPU_ROUTINE(parallelism='[seq]')
+
         integer, intent(in)          :: patch_id
         real(wp), dimension(3, 3, 3) :: rotation
         real(wp)                     :: angle

src/simulation/m_ibm.fpp

@@ -76,14 +76,17 @@ contains
 
         call nvtxStartRange("SETUP-IBM-MODULE")
 
+        $:GPU_UPDATE(device='[patch_ib(1:num_ibs)]')
+
         ! do all set up for moving immersed boundaries
+        $:GPU_PARALLEL_LOOP(private='[i]')
         do i = 1, num_ibs
             if (patch_ib(i)%moving_ibm /= 0) then
                 call s_compute_moment_of_inertia(i, patch_ib(i)%angular_vel)
             end if
             call s_update_ib_rotation_matrix(i)
         end do
-        $:GPU_UPDATE(device='[patch_ib(1:num_ibs)]')
+        $:END_GPU_PARALLEL_LOOP()
 
         ! allocate some arrays for MPI communication, if required by this simulation
 #ifdef MFC_MPI
@@ -868,13 +871,13 @@ contains
         $:END_GPU_PARALLEL_LOOP()
 
         ! recalulcate the rotation matrix based upon the new angles
+        $:GPU_PARALLEL_LOOP(private='[i]')
         do i = 1, num_ibs
             if (patch_ib(i)%moving_ibm /= 0) then
                 call s_update_ib_rotation_matrix(i)
             end if
         end do
-
-        $:GPU_UPDATE(device='[patch_ib]')
+        $:END_GPU_PARALLEL_LOOP()
 
         ! recompute the new ib_patch locations and broadcast them.
         call nvtxStartRange("APPLY-IB-PATCHES")
@@ -933,7 +936,7 @@ contains
         $:GPU_PARALLEL_LOOP(private='[ib_idx, ib_idx_temp, encoded_ib_idx, fluid_idx, radial_vector, local_force_contribution, &
                             & cell_volume, local_torque_contribution, dynamic_viscosity, viscous_stress_div, &
                             & viscous_stress_div_1, viscous_stress_div_2, dx, dy, dz]', copy='[forces, torques]', &
-                            & copyin='[patch_ib, dynamic_viscosities]', collapse=3)
+                            & copyin='[dynamic_viscosities]', collapse=3)
         do i = 0, m
             do j = 0, n
                 do k = 0, p
@@ -1122,6 +1125,8 @@ contains
     !> Computes the moment of inertia for an immersed boundary
     subroutine s_compute_moment_of_inertia(ib_idx, axis)
 
+        $:GPU_ROUTINE(parallelism='[seq]')
+
         real(wp), dimension(3), intent(in) :: axis  !< the axis about which we compute the moment. Only required in 3D.
         integer, intent(in)                :: ib_idx
         real(wp)                           :: moment, distance_to_axis, cell_volume
@@ -1158,13 +1163,10 @@ contains
 
             ib_marker = patch_ib(ib_idx)%gbl_patch_id
 
-            $:GPU_PARALLEL_LOOP(private='[position, closest_point_along_axis, vector_to_axis, distance_to_axis]', copy='[moment, &
-                                & count]', copyin='[ib_marker, cell_volume, normal_axis]', collapse=3)
             do i = 0, m
                 do j = 0, n
                     do k = 0, p
                         if (ib_markers%sf(i, j, k) == ib_marker) then
-                            $:GPU_ATOMIC(atomic='update')
                             count = count + 1  ! increment the count of total cells in the boundary
 
                             ! get the position in local coordinates so that the axis passes through 0, 0, 0
@@ -1182,13 +1184,11 @@ contains
                             distance_to_axis = dot_product(vector_to_axis, vector_to_axis)  ! saves the distance to the axis squared
 
                             ! compute the position component of the moment
-                            $:GPU_ATOMIC(atomic='update')
                             moment = moment + distance_to_axis
                         end if
                     end do
                 end do
             end do
-            $:END_GPU_PARALLEL_LOOP()
 
             ! write the final moment assuming the points are all uniform density
             patch_ib(ib_idx)%moment = moment*patch_ib(ib_idx)%mass/(count*cell_volume)

src/simulation/m_time_steppers.fpp

@@ -714,6 +714,7 @@ contains
 
         if (moving_immersed_boundary_flag) call s_compute_ib_forces(q_prim_vf, fluid_pp)
 
+        $:GPU_PARALLEL_LOOP(private='[i, gbl_id]', copyin='[s]')
         do i = 1, num_ibs
             if (s == 1) then
                 patch_ib(i)%step_vel = patch_ib(i)%vel
@@ -759,8 +760,8 @@ contains
                          & 2)*patch_ib(i)%z_centroid + rk_coef(s, 3)*patch_ib(i)%vel(3)*dt)/rk_coef(s, 4)
             end if
         end do
+        $:END_GPU_PARALLEL_LOOP()
 
-        $:GPU_UPDATE(device='[patch_ib]')
         call s_update_mib(num_ibs)
 
         call nvtxEndRange