fix: lazy-allocate patch_ib to restore FP-stability CI

num_ib_patches_max = 2050000 caused a ~2.25 GB unconditional heap
allocation + full init loop in s_assign_default_values_to_user_inputs
on every startup, crashing the CI debug build even for cases with no IBM.
Also, s_reduce_ib_patch_array had a 2.25 GB stack-allocated local array
that caused SIGSEGV for IBM cases.

Introduce num_ib_patches_max_namelist = 50000 (restoring the pre-particle-bed
budget) for the initial allocation. s_generate_particle_beds grows patch_ib
to num_ib_patches_max via MOVE_ALLOC only when particle beds are actually
being generated. s_reduce_ib_patch_array now uses a heap-allocated local
array sized to num_ibs instead of the full num_ib_patches_max.

Author: sbryngelson

src/common/m_constants.fpp

@@ -16,15 +16,17 @@ module m_constants
     real(wp), parameter  :: verysmall = 1.e-12_wp       !< Very small number
     !> Radius cutoff to avoid division by zero for 3D spherical harmonic patch (geometry 14)
     real(wp), parameter :: small_radius = 1.e-32_wp
-    integer, parameter  :: num_stcls_min = 5                   !< Minimum # of stencils
-    integer, parameter  :: path_len = 400                      !< Maximum path length
-    integer, parameter  :: name_len = 50                       !< Maximum name length
-    integer, parameter  :: dflt_int = -100                     !< Default integer value
-    integer, parameter  :: fourier_rings = 5                   !< Fourier filter ring limit
-    integer, parameter  :: num_fluids_max = 10                 !< Maximum number of fluids in the simulation
-    integer, parameter  :: num_probes_max = 10                 !< Maximum number of flow probes in the simulation
-    integer, parameter  :: num_patches_max = 10                !< Maximum number of IC patches
-    integer, parameter  :: num_ib_patches_max = 2050000        !< Maximum number of immersed boundary patches (patch_ib)
+    integer, parameter  :: num_stcls_min = 5             !< Minimum # of stencils
+    integer, parameter  :: path_len = 400                !< Maximum path length
+    integer, parameter  :: name_len = 50                 !< Maximum name length
+    integer, parameter  :: dflt_int = -100               !< Default integer value
+    integer, parameter  :: fourier_rings = 5             !< Fourier filter ring limit
+    integer, parameter  :: num_fluids_max = 10           !< Maximum number of fluids in the simulation
+    integer, parameter  :: num_probes_max = 10           !< Maximum number of flow probes in the simulation
+    integer, parameter  :: num_patches_max = 10          !< Maximum number of IC patches
+    integer, parameter  :: num_ib_patches_max = 2050000  !< Maximum number of immersed boundary patches (patch_ib)
+    !> Max patches readable from the namelist; patch_ib grows to num_ib_patches_max only when particle beds are used
+    integer, parameter  :: num_ib_patches_max_namelist = 50000
     integer, parameter  :: num_local_ibs_max = 2000            !< Maximum number of immersed boundary patches (patch_ib)
     integer, parameter  :: num_particle_beds_max = 10          !< Maximum number of particle bed patch specifications
     integer, parameter  :: num_bc_patches_max = 10             !< Maximum number of boundary condition patches

src/simulation/m_global_parameters.fpp

@@ -829,8 +829,8 @@ contains
             particle_bed(i)%seed = 0
         end do
 
-        allocate (patch_ib(num_ib_patches_max))
-        do i = 1, num_ib_patches_max
+        allocate (patch_ib(num_ib_patches_max_namelist))
+        do i = 1, num_ib_patches_max_namelist
             patch_ib(i)%gbl_patch_id = i
             patch_ib(i)%geometry = dflt_int
             patch_ib(i)%x_centroid = 0._wp

src/simulation/m_particle_bed.fpp

@@ -40,6 +40,18 @@ contains
 
         if (num_particle_beds == 0) return
 
+        ! Grow patch_ib from the namelist-sized allocation to the full capacity needed for particle beds
+        if (size(patch_ib) < num_ib_patches_max) then
+            block
+                type(ib_patch_parameters), allocatable :: tmp(:)
+                integer                                :: n
+                n = size(patch_ib)
+                call move_alloc(patch_ib, tmp)
+                allocate (patch_ib(num_ib_patches_max))
+                patch_ib(1:n) = tmp
+            end block
+        end if
+
         do b = 1, num_particle_beds
             xmin = particle_bed(b)%x_centroid - 0.5_wp*particle_bed(b)%length_x
             xmax = particle_bed(b)%x_centroid + 0.5_wp*particle_bed(b)%length_x

src/simulation/m_start_up.fpp

@@ -1022,7 +1022,7 @@ contains
             "on CPUs"
 #endif
         else
-            allocate (patch_ib(num_ib_patches_max))
+            allocate (patch_ib(num_ib_patches_max_namelist))
         end if
 
         call s_mpi_bcast_user_inputs()
@@ -1210,11 +1210,11 @@ contains
     !! the local computational domain.
     subroutine s_reduce_ib_patch_array()
 
-        type(ib_patch_parameters), dimension(num_ib_patches_max) :: patch_ib_gbl
-        real(wp), dimension(3)                                   :: centroid
-        integer                                                  :: i, j
-        integer                                                  :: num_aware_ibs
-        logical                                                  :: is_in_neighborhood, is_local
+        type(ib_patch_parameters), allocatable :: patch_ib_gbl(:)
+        real(wp), dimension(3)                 :: centroid
+        integer                                :: i, j
+        integer                                :: num_aware_ibs
+        logical                                :: is_in_neighborhood, is_local
 
         ! do all set up for moving immersed boundaries
 
@@ -1226,6 +1226,7 @@ contains
             end if
         end do
 
+        allocate (patch_ib_gbl(num_ibs))
         patch_ib_gbl(1:num_ibs) = patch_ib(1:num_ibs)
         call get_neighbor_bounds()  ! make sure the bounds of the neighbors are correctly set up
         call s_compute_ib_neighbor_ranks()  ! build lookup of all neighbor MPI ranks
@@ -1244,7 +1245,7 @@ contains
 #ifdef MFC_MPI
         ! fallback for 1-rank case
         if (num_procs == 1) then
-            patch_ib(:) = patch_ib_gbl(1:num_aware_ibs)
+            patch_ib(1:num_gbl_ibs) = patch_ib_gbl(1:num_gbl_ibs)
         else
             ! determine the set of patches owned by local rank
             num_local_ibs = 0
@@ -1269,9 +1270,11 @@ contains
         end if
 #else
         ! reduce the size of the array for local simulation in no-MPI case
-        patch_ib(:) = patch_ib_gbl(1:num_aware_ibs)
+        patch_ib(1:num_gbl_ibs) = patch_ib_gbl(1:num_gbl_ibs)
 #endif
 
+        deallocate (patch_ib_gbl)
+
         @:ALLOCATE(ib_gbl_idx_lookup(1:num_gbl_ibs))
 
     end subroutine s_reduce_ib_patch_array