1536 ib neighborhood radius automatic size checking

docs/documentation/case.md

@@ -414,7 +414,7 @@ Additional details on this specification can be found in [NACA airfoil](https://
 
 - `ib_coefficient_of_friction` is the coefficient of friction used in IB collisions.
 
-- `ib_neighborhood_radius` controls the size of the neighborhood size. This value defaults to 1, which indicates that any given rank is aware of IBs up to 1 ranks away. This parameter is required to strong-scale a case when IBs eventually grow to be larger than one full processor domain wide.
+- `ib_neighborhood_radius` controls the size of the neighborhood size. A value of $r$ indicates that any given rank is aware of IBs up to $r$ ranks away. This value defaults to 0, which leaves the radius unset so that it is selected automatically. This parameter is required to strong-scale a case when IBs eventually grow to be larger than one full processor domain wide.
 
 #### Particle Clouds
 

src/simulation/m_global_parameters.fpp

@@ -486,7 +486,7 @@ contains
         bub_pp%R_g = dflt_real; R_g = dflt_real
 
         ! Immersed Boundaries (sim-specific extras)
-        ib_neighborhood_radius = 1
+        ib_neighborhood_radius = 0
         collision_model = 0
         coefficient_of_restitution = dflt_real
         collision_time = dflt_real

src/simulation/m_ib_patches.fpp

@@ -23,7 +23,7 @@ module m_ib_patches
     implicit none
 
     private; public :: s_apply_ib_patches, s_update_ib_rotation_matrix, s_instantiate_STL_models, s_decode_patch_periodicity, &
-        & s_encode_patch_periodicity, s_initialize_ib_airfoils, s_get_periodicities
+        & s_encode_patch_periodicity, s_initialize_ib_airfoils, s_get_periodicities, s_get_ib_bound
 
 contains
 
@@ -65,7 +65,7 @@ contains
                             call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, zp, encoded_patch_id)
 
                             ! find the indices to the left and right of the IB in i, j, k
-                            call get_bounding_indices(patch_ib(patch_id), center, il, ir, jl, jr, kl, kr)
+                            call s_get_bounding_indices(patch_ib(patch_id), center, il, ir, jl, jr, kl, kr)
 
                             ! skip patches whose bounding box does not overlap this rank's domain
                             if (ir < il .or. jr < jl .or. kr < kl) cycle
@@ -137,7 +137,7 @@ contains
                         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
-                        call get_bounding_indices(patch_ib(patch_id), center, il, ir, jl, jr, kl, kr)
+                        call s_get_bounding_indices(patch_ib(patch_id), center, il, ir, jl, jr, kl, kr)
 
                         ! skip patches whose bounding box does not overlap this rank's domain
                         if (ir < il .or. jr < jl) cycle
@@ -218,7 +218,7 @@ contains
                             call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, zp, encoded_patch_id)
 
                             ! find the indices to the left and right of the IB in i, j, k
-                            call get_bounding_indices(patch_ib(patch_id), center, il, ir, jl, jr, kl, kr)
+                            call s_get_bounding_indices(patch_ib(patch_id), center, il, ir, jl, jr, kl, kr)
 
                             do k = kl, kr
                                 do j = jl, jr
@@ -287,7 +287,7 @@ contains
                         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
-                        call get_bounding_indices(patch_ib(patch_id), center, il, ir, jl, jr, kl, kr)
+                        call s_get_bounding_indices(patch_ib(patch_id), center, il, ir, jl, jr, kl, kr)
 
                         do j = jl, jr
                             do i = il, ir
@@ -458,31 +458,67 @@ contains
 
     end subroutine s_update_ib_rotation_matrix
 
-    subroutine get_bounding_indices(patch, center, il, ir, jl, jr, kl, kr)
+    subroutine s_get_ib_bound(patch, bound)
 
         $:GPU_ROUTINE(parallelism='[seq]')
 
         type(ib_patch_parameters), intent(in) :: patch
-        real(wp), dimension(3), intent(in)    :: center
-        integer, intent(out)                  :: il, ir, jl, jr, kl, kr
-        real(wp), dimension(3)                :: bbox_min, bbox_max, local_corner, world_corner
+        real(wp), intent(out)                 :: bound
         real(wp), dimension(2)                :: lx, ly, lz
-        integer                               :: cx, cy, cz
-        logical                               :: outside_domain
 
         if (patch%geometry == 2 .or. patch%geometry == 8) then
             ! circle and sphere geometries
-            bbox_min = center - patch%radius
-            bbox_max = center + patch%radius
+            bound = patch%radius
         else if (patch%geometry == 3) then
-            ! rectangular geometries
-            bbox_min = center - 0.5_wp*sqrt(patch%length_x**2 + patch%length_y**2)
-            bbox_max = center + 0.5_wp*sqrt(patch%length_x**2 + patch%length_y**2)
+            bound = 0.5_wp*sqrt(patch%length_x**2 + patch%length_y**2)
         else if (patch%geometry == 4 .or. patch%geometry == 11) then
-            ! airfoil geometries TODO :: This can be better optimized, since airfoils are typically very long in one dimension
-            bbox_min = center - ib_airfoil(patch%airfoil_id)%c
-            bbox_max = center + ib_airfoil(patch%airfoil_id)%c
+            ! rectangular geometries
+            bound = ib_airfoil(patch%airfoil_id)%c
         else if (patch%geometry == 5) then
+            ! STL model geometry
+            lx(1) = stl_bounding_boxes(patch%model_id, 1, 1)
+            lx(2) = stl_bounding_boxes(patch%model_id, 1, 3)
+            ly(1) = stl_bounding_boxes(patch%model_id, 2, 1)
+            ly(2) = stl_bounding_boxes(patch%model_id, 2, 3)
+
+            bound = 0.5_wp*sqrt((lx(2) - lx(1))**2 + (ly(2) - ly(1))**2)
+        else if (patch%geometry == 6) then
+            ! ellipse geometry
+            bound = 0.5_wp*max(patch%length_x, patch%length_y)
+        else if (patch%geometry == 9) then
+            ! cuboid geometries
+            bound = 0.5_wp*sqrt(patch%length_x**2 + patch%length_y**2 + patch%length_z**2)
+        else if (patch%geometry == 10) then
+            ! cylinder geometry
+            bound = sqrt(patch%radius**2 + patch%length_x**2)
+        else if (patch%geometry == 12) then
+            ! Local-space bounding box extents (min=1, max=2 in the third index)
+            lx(1) = stl_bounding_boxes(patch%model_id, 1, 1) + patch%centroid_offset(1)
+            lx(2) = stl_bounding_boxes(patch%model_id, 1, 3) + patch%centroid_offset(1)
+            ly(1) = stl_bounding_boxes(patch%model_id, 2, 1) + patch%centroid_offset(2)
+            ly(2) = stl_bounding_boxes(patch%model_id, 2, 3) + patch%centroid_offset(2)
+            lz(1) = stl_bounding_boxes(patch%model_id, 3, 1) + patch%centroid_offset(3)
+            lz(2) = stl_bounding_boxes(patch%model_id, 3, 3) + patch%centroid_offset(3)
+
+            bound = 0.5_wp*sqrt((lx(2) - lx(1))**2 + (ly(2) - ly(1))**2 + (lz(2) - lz(1))**2)
+        end if
+
+    end subroutine s_get_ib_bound
+
+    subroutine s_get_bounding_indices(patch, center, il, ir, jl, jr, kl, kr)
+
+        $:GPU_ROUTINE(parallelism='[seq]')
+
+        type(ib_patch_parameters), intent(in) :: patch
+        real(wp), dimension(3), intent(in)    :: center
+        integer, intent(out)                  :: il, ir, jl, jr, kl, kr
+        real(wp), dimension(3)                :: bbox_min, bbox_max, local_corner, world_corner
+        real(wp), dimension(2)                :: lx, ly, lz
+        real(wp)                              :: bound
+        integer                               :: cx, cy, cz
+        logical                               :: outside_domain
+
+        if (patch%geometry == 5) then
             ! STL model geometry
             lx(1) = stl_bounding_boxes(patch%model_id, 1, 1) + patch%centroid_offset(1)
             lx(2) = stl_bounding_boxes(patch%model_id, 1, 3) + patch%centroid_offset(1)
@@ -502,18 +538,6 @@ contains
                     bbox_max(2) = max(bbox_max(2), world_corner(2))
                 end do
             end do
-        else if (patch%geometry == 6) then
-            ! ellipse geometry
-            bbox_min = center - 0.5_wp*max(patch%length_x, patch%length_y)
-            bbox_max = center + 0.5_wp*max(patch%length_x, patch%length_y)
-        else if (patch%geometry == 9) then
-            ! cuboid geometries
-            bbox_min = center - 0.5_wp*sqrt(patch%length_x**2 + patch%length_y**2 + patch%length_z**2)
-            bbox_max = center + 0.5_wp*sqrt(patch%length_x**2 + patch%length_y**2 + patch%length_z**2)
-        else if (patch%geometry == 10) then
-            ! cylinder geometry
-            bbox_min = center - sqrt(patch%radius**2 + patch%length_x**2)
-            bbox_max = center + sqrt(patch%radius**2 + patch%length_x**2)
         else if (patch%geometry == 12) then
             ! Local-space bounding box extents (min=1, max=2 in the third index)
             lx(1) = stl_bounding_boxes(patch%model_id, 1, 1) + patch%centroid_offset(1)
@@ -540,6 +564,11 @@ contains
                     end do
                 end do
             end do
+        else
+            ! All other IBs
+            call s_get_ib_bound(patch, bound)
+            bbox_min = center - bound
+            bbox_max = center + bound
         end if
 
         ! completely skip patches whose bounding box does not overlap this rank's domain
@@ -566,7 +595,7 @@ contains
         call get_indices_from_bounds(bbox_min(2), bbox_max(2), y_cc, jl, jr)
         if (num_dims == 3) call get_indices_from_bounds(bbox_min(3), bbox_max(3), z_cc, kl, kr)
 
-    end subroutine get_bounding_indices
+    end subroutine s_get_bounding_indices
 
     subroutine get_indices_from_bounds(left_bound, right_bound, cell_centers, left_index, right_index)
 

src/simulation/m_start_up.fpp

@@ -898,7 +898,9 @@ contains
                 type(ib_patch_parameters), allocatable :: particle_cloud_ibs(:)
                 integer                                :: num_particle_cloud_ibs
 
-                call get_neighbor_bounds()
+                call s_instantiate_STL_models()
+                call s_initialize_ib_airfoils()
+                call s_get_neighbor_bounds()
 
                 if (cfl_dt .and. n_start > 0) then
                     call s_read_ib_restart_data(n_start)
@@ -911,8 +913,6 @@ contains
                 else
                     call s_generate_particle_clouds(particle_cloud_ibs, num_particle_cloud_ibs)
                 end if
-                call s_instantiate_STL_models()
-                call s_initialize_ib_airfoils()
                 call s_reduce_ib_patch_array(particle_cloud_ibs, num_particle_cloud_ibs)
                 deallocate (particle_cloud_ibs)
             end block
@@ -1492,10 +1492,10 @@ contains
 
     end subroutine s_compute_ib_neighbor_ranks
 
-    subroutine get_neighbor_bounds()
+    subroutine s_get_neighbor_bounds()
 
-        real(wp) :: beg_val, end_val, recv_val
-        integer  :: k, send_neighbor, recv_neighbor, ierr
+        real(wp) :: beg_val, end_val, recv_val, bound, max_ib_bound, local_rank_width, min_rank_width
+        integer  :: k, send_neighbor, recv_neighbor, ierr, temporary_radius
 
         ! Default: unbounded in all directions (covers single-rank and no-MPI cases)
 
@@ -1507,6 +1507,32 @@ contains
         neighbor_domain_z%end = huge(0._wp)
 
 #ifdef MFC_MPI
+        ! perform setup if we are doing automatic radius checking
+        if (ib_neighborhood_radius < 1) then
+            ib_neighborhood_radius = 0  ! ensure we are starting with 0 neighborhood radius
+
+            ! determine the maximum length of space that needs to be contained by the neighborhood
+            max_ib_bound = -1._wp
+            do k = 1, num_ibs
+                call s_get_ib_bound(patch_ib(k), bound)
+                max_ib_bound = max(max_ib_bound, bound)
+            end do
+            do k = 1, num_particle_clouds
+                max_ib_bound = max(max_ib_bound, particle_cloud(k)%radius)
+            end do
+
+            ! determine the upper bound on the size
+            local_rank_width = -1._wp
+            #:for X, ID, DIM in [('x', 1, 'm'), ('y', 2, 'n'), ('z', 3, 'p')]
+                if (num_dims >= ${ID}$) local_rank_width = max(local_rank_width, abs(${X}$_cb(${DIM}$) - ${X}$_cb(-1)))
+            #:endfor
+            call s_mpi_allreduce_min(local_rank_width, min_rank_width)
+
+            ! approximate the size of the neighborhood with a local 1.1x fudge factor for safety, lower bound of 1
+            ib_neighborhood_radius = max(1, ceiling(1.1_wp*max_ib_bound/(min_rank_width)))
+            if (proc_rank == 0) print *, "Automatic choice of ib_neighborhood_radius selected: ", ib_neighborhood_radius
+        end if
+
         ! For each direction, propagate the left/right boundary edges outward ib_neighborhood_radius hops. After k rounds: beg_val =
         ! left edge of the rank k hops to the left; end_val = right edge of the rank k hops to the right.
         #:for X, ID, TAG, DIM in [('x', 1, 100, 'm'), ('y', 2, 102, 'n'), ('z', 3, 104, 'p')]
@@ -1541,6 +1567,6 @@ contains
         #:endfor
 #endif
 
-    end subroutine get_neighbor_bounds
+    end subroutine s_get_neighbor_bounds
 
 end module m_start_up

toolchain/mfc/params/definitions.py

@@ -376,7 +376,7 @@ def get_value_label(param_name: str, value: int) -> str:
     "num_fluids": {"min": 1, "max": NF},
     "num_patches": {"min": 0, "max": NUM_PATCHES_MAX},
     "num_ibs": {"min": 0},
-    "ib_neighborhood_radius": {"min": 1},
+    "ib_neighborhood_radius": {"min": 0},
     "num_source": {"min": 1},
     "num_probes": {"min": 1},
     "nb": {"min": 1},