@@ -1107,35 +1107,33 @@ contains
11071107 end subroutine s_compute_centroid_offset
11081108
11091109 !> Computes the moment of inertia for an immersed boundary
1110- subroutine s_compute_moment_of_inertia (ib_marker , axis )
1110+ subroutine s_compute_moment_of_inertia (ib_idx , axis )
11111111
11121112 real (wp), dimension (3 ), intent (in ) :: axis !< the axis about which we compute the moment. Only required in 3D .
1113- integer , intent (in ) :: ib_marker
1113+ integer , intent (in ) :: ib_idx
11141114 real (wp) :: moment, distance_to_axis, cell_volume
11151115 real (wp), dimension (3 ) :: position, closest_point_along_axis, vector_to_axis, normal_axis
1116- integer :: i, j, k, count
1116+ integer :: i, j, k, count, ib_marker
11171117
11181118 if (p == 0 ) then
11191119 normal_axis = [0 , 0 , 1 ]
11201120 else if (sqrt (sum (axis** 2 )) < sgm_eps) then
11211121 ! if the object is not actually rotating at this time, return a dummy value and exit
1122- patch_ib(ib_marker )%moment = 1._wp
1122+ patch_ib(ib_idx )%moment = 1._wp
11231123 return
11241124 else
11251125 normal_axis = axis/ sqrt (sum (axis))
11261126 end if
11271127
11281128 ! if the IB is in 2D or a 3D sphere, we can compute this exactly
1129- if (patch_ib(ib_marker)%geometry == 2 ) then ! circle
1130- patch_ib(ib_marker)%moment = 0.5_wp * patch_ib(ib_marker)%mass* (patch_ib(ib_marker)%radius)** 2
1131- else if (patch_ib(ib_marker)%geometry == 3 ) then ! rectangle
1132- patch_ib(ib_marker)%moment = patch_ib(ib_marker)%mass* (patch_ib(ib_marker)%length_x** 2 + patch_ib(ib_marker) &
1133- & %length_y** 2 )/ 6._wp
1134- else if (patch_ib(ib_marker)%geometry == 6 ) then ! ellipse
1135- patch_ib(ib_marker)%moment = 0.0625_wp * patch_ib(ib_marker)%mass* (patch_ib(ib_marker)%length_x** 2 + patch_ib(ib_marker) &
1136- & %length_y** 2 )
1137- else if (patch_ib(ib_marker)%geometry == 8 ) then ! sphere
1138- patch_ib(ib_marker)%moment = 0.4 * patch_ib(ib_marker)%mass* (patch_ib(ib_marker)%radius)** 2
1129+ if (patch_ib(ib_idx)%geometry == 2 ) then ! circle
1130+ patch_ib(ib_idx)%moment = 0.5_wp * patch_ib(ib_idx)%mass* (patch_ib(ib_idx)%radius)** 2
1131+ else if (patch_ib(ib_idx)%geometry == 3 ) then ! rectangle
1132+ patch_ib(ib_idx)%moment = patch_ib(ib_idx)%mass* (patch_ib(ib_idx)%length_x** 2 + patch_ib(ib_idx) %length_y** 2 )/ 6._wp
1133+ else if (patch_ib(ib_idx)%geometry == 6 ) then ! ellipse
1134+ patch_ib(ib_idx)%moment = 0.0625_wp * patch_ib(ib_idx)%mass* (patch_ib(ib_idx)%length_x** 2 + patch_ib(ib_idx) %length_y** 2 )
1135+ else if (patch_ib(ib_idx)%geometry == 8 ) then ! sphere
1136+ patch_ib(ib_idx)%moment = 0.4 * patch_ib(ib_idx)%mass* (patch_ib(ib_idx)%radius)** 2
11391137 else ! we do not have an analytic moment of inertia calculation and need to approximate it directly via a sum
11401138 count = 0
11411139 moment = 0._wp
@@ -1145,6 +1143,8 @@ contains
11451143 cell_volume = cell_volume* (z_cc(1 ) - z_cc(0 ))
11461144 end if
11471145
1146+ ib_marker = patch_ib(ib_idx)%gbl_patch_id
1147+
11481148 $:GPU_PARALLEL_LOOP(private= ' [position, closest_point_along_axis, vector_to_axis, distance_to_axis]' , copy= ' [moment, &
11491149 & count]' , copyin= ' [ib_marker, cell_volume, normal_axis]' , collapse= 3 )
11501150 do i = 0 , m
@@ -1155,12 +1155,12 @@ contains
11551155 count = count + 1 ! increment the count of total cells in the boundary
11561156
11571157 ! get the position in local coordinates so that the axis passes through 0 , 0 , 0
1158- if (p == 0 ) then
1159- position = [x_cc(i), y_cc(j), 0._wp ] - [patch_ib(ib_marker )%x_centroid, &
1160- & patch_ib(ib_marker)%y_centroid, 0._wp ]
1158+ if (num_dims < 3 ) then
1159+ position = [x_cc(i), y_cc(j), 0._wp ] - [patch_ib(ib_idx )%x_centroid, patch_ib(ib_idx)%y_centroid , &
1160+ & 0._wp ]
11611161 else
1162- position = [x_cc(i), y_cc(j), z_cc(k)] - [patch_ib(ib_marker )%x_centroid, &
1163- & patch_ib(ib_marker )%y_centroid, patch_ib(ib_marker )%z_centroid]
1162+ position = [x_cc(i), y_cc(j), z_cc(k)] - [patch_ib(ib_idx )%x_centroid, &
1163+ & patch_ib(ib_idx )%y_centroid, patch_ib(ib_idx )%z_centroid]
11641164 end if
11651165
11661166 ! project the position along the axis to find the closest distance to the rotation axis
@@ -1178,8 +1178,7 @@ contains
11781178 $:END_GPU_PARALLEL_LOOP()
11791179
11801180 ! write the final moment assuming the points are all uniform density
1181- patch_ib(ib_marker)%moment = moment* patch_ib(ib_marker)%mass/ (count* cell_volume)
1182- $:GPU_UPDATE(device= ' [patch_ib(ib_marker)%moment]' )
1181+ patch_ib(ib_idx)%moment = moment* patch_ib(ib_idx)%mass/ (count* cell_volume)
11831182 end if
11841183
11851184 end subroutine s_compute_moment_of_inertia
@@ -1506,4 +1505,21 @@ contains
15061505
15071506 end subroutine s_handoff_ib_ownership
15081507
1508+ subroutine get_neighborhood_idx(gbl_idx, neighborhood_idx)
1509+
1510+ integer, intent(in) :: gbl_idx
1511+ integer, intent(out) :: neighborhood_idx
1512+ integer :: i
1513+
1514+ neighborhood_idx = -1
1515+
1516+ do i = 1, num_ibs
1517+ if (patch_ib(i)%gbl_patch_id == gbl_idx) then
1518+ neighborhood_idx = i
1519+ exit
1520+ end if
1521+ end do
1522+
1523+ end subroutine get_neighborhood_idx
1524+
15091525end module m_ibm
0 commit comments