Move Particle Cloud Generation to Pre-Process

docs/documentation/contributing.md

@@ -218,7 +218,7 @@ Both human reviewers and AI code reviewers reference this section.
 ### Parameter Plumbing
 
 - **Derived-type parameters are not auto-broadcast.** `generated_bcast.fpp` covers namelist *scalars* only. Each derived type (`chem_params`, `lag_params`, `rburn`) needs a hand-written `_emit_<name>` in `toolchain/mfc/params/generators/fortran_gen.py` plus its call site in that generator's simulation branch, and, if it is read on device, an explicit ``$:GPU_UPDATE(device='[name]')`` in both the target's `m_global_parameters.fpp` and `src/simulation/m_start_up.fpp` — `GPU_DECLARE` alone does not make it device-resident. Regrouping existing scalars into a derived type silently drops their broadcast, leaving every non-root rank holding the `dflt_real` sentinel. Single-rank golden files cannot catch this, so pair such a change with a `ppn=2` test and confirm it fails without the emitter.
-- **A `patch_ib` member that immersed-boundary ghost-point code reads must also be set in `s_add_cloud_particle`** (`src/simulation/m_particle_cloud.fpp`). `particle_cloud_ibs` is allocated without default initialization, and `s_reduce_ib_patch_array` copies the whole struct into `patch_ib`, overwriting the defaults assigned in `s_assign_default_values_to_user_inputs`. Anything left unset reaches the solver as uninitialized memory, and only where the allocation is not already zero-filled. A platform-only NaN is the signature of this class: a garbage `v_blow` once failed an AMD lane with `ICFL is NaN` while every NVIDIA lane and all local runs passed.
+- **A `patch_ib` member that immersed-boundary ghost-point code reads must also be set for particle-cloud IBs in `s_assign_particle_cloud_ib_defaults`** (`src/simulation/m_start_up.fpp`). Pre-process writes only position, kinematics and radius to the IB state file; simulation builds every other property there, writing into a reused `patch_ib` slot. Anything it leaves unset keeps whatever that slot held, which may be a namelist patch's value or uninitialized memory, and shows up only where that memory is not already zero-filled. A platform-only NaN is the signature of this class: a garbage `v_blow` once failed an AMD lane with `ICFL is NaN` while every NVIDIA lane and all local runs passed.
 - **Runtime checks go where they run.** Shared constraints belong in `src/common/m_checker_common.fpp`, simulation-only ones in `src/simulation/m_checker.fpp`, and pre- and post-process ones in their own `m_checker.fpp`. Those two `s_check_inputs` are currently empty; that is still the correct home for their checks, not `m_checker_common`.
 - **Analytic initial conditions are compiled into the binary** and their expressions are AST-validated at case load, so syntax errors and unknown variables surface immediately and by name. Each IC variable maps to an `eqn_idx` expression in `QPVF_IDX_VARS` (`toolchain/mfc/case.py`); adding a patch-settable conserved variable means updating that map and the Fortran `eqn_idx` builder together, because a mismatch is a silent wrong index.
 - **Under `--case-optimization` the baked-in constants are dropped from the namelist**, so changing one requires a rebuild rather than a case-file edit.

src/common/m_helper.fpp

@@ -10,6 +10,7 @@ module m_helper
 
     use m_derived_types
     use m_global_parameters
+    use m_constants, only: BC_PERIODIC
     use ieee_arithmetic  !< For checking NaN
 
     implicit none
@@ -18,7 +19,8 @@ module m_helper
     public :: s_comp_n_from_prim, s_comp_n_from_cons, s_initialize_bubbles_model, s_initialize_nonpoly, s_simpson, s_transcoeff, &
         & s_int_to_str, s_transform_vec, s_transform_triangle, s_transform_model, s_swap, f_cross, f_create_transform_matrix, &
         & f_create_bbox, s_print_2D_array, f_xor, f_logical_to_int, associated_legendre, real_ylm, double_factorial, factorial, &
-        & f_cut_on, f_cut_off, s_downsample_data, s_upsample_data, s_cross_product, f_unit_vector, s_prng, modmul
+        & f_cut_on, f_cut_off, s_downsample_data, s_upsample_data, s_cross_product, f_unit_vector, s_prng, modmul, &
+        & f_local_rank_owns_location
 
 contains
 
@@ -699,4 +701,45 @@ contains
 
     end subroutine s_upsample_data
 
+    !> @brief True if `location` falls within this rank's own subdomain (a strict partition - each location is owned by exactly one
+    !! rank, unlike an overlapping multi-rank ghost-stencil neighborhood). Used by both pre_process (to decide which
+    !! generated/namelist IBs a rank writes to its own restart_data/ib_state_0.dat chunk) and simulation (to decide which IBs a rank
+    !! owns for its own periodic IB-state writes) so the two stay consistent. glb_bounds_in is the global domain extent, used only
+    !! to project a location that falls just outside the domain (floating-point edge case) onto the domain so some rank still claims
+    !! it.
+    function f_local_rank_owns_location(location, glb_bounds_in) result(owns_location)
+
+        $:GPU_ROUTINE(parallelism='[seq]')
+
+        real(wp), dimension(3), intent(in)          :: location
+        type(bounds_info), dimension(3), intent(in) :: glb_bounds_in
+        logical                                     :: owns_location
+        real(wp), dimension(3)                      :: projected_location
+
+        owns_location = .true.
+
+#ifdef MFC_MPI
+        if (num_procs > 1) then
+            projected_location(:) = location(:)
+
+            ! catch the edge case where the location lies just outside the computational domain
+            #:for X, ID, DIM in [('x', 1, 'm'), ('y', 2, 'n'), ('z', 3, 'p')]
+                if (num_dims >= ${ID}$) then
+                    if (ib_bc_${X}$%beg /= BC_PERIODIC) then
+                        ! if it is outside the domain in one direction, project it somewhere inside so at least one rank owns it
+                        if (location(${ID}$) < glb_bounds_in(${ID}$)%beg) then
+                            projected_location(${ID}$) = glb_bounds_in(${ID}$)%beg
+                        else if (glb_bounds_in(${ID}$)%end < location(${ID}$)) then
+                            projected_location(${ID}$) = glb_bounds_in(${ID}$)%end - 1.0e-10_wp
+                        end if
+                    end if
+                    owns_location = owns_location .and. ${X}$_cb(-1) <= projected_location(${ID}$) &
+                        & .and. projected_location(${ID}$) < ${X}$_cb(${DIM}$)
+                end if
+            #:endfor
+        end if
+#endif
+
+    end function f_local_rank_owns_location
+
 end module m_helper

src/post_process/m_global_parameters.fpp

@@ -93,6 +93,7 @@ module m_global_parameters
     !> @name Boundary conditions in the x-, y- and z-coordinate directions
     !> @{
     type(int_bounds_info) :: bc_x, bc_y, bc_z
+    type(int_bounds_info) :: ib_bc_x, ib_bc_y, ib_bc_z  !< bc_x/y/z before decomposition overwrites them with MPI neighbor ranks
     type(bc_xyz_info)     :: bc
     !> @}
 
@@ -364,7 +365,7 @@ contains
 
         ! Particle clouds expand into individual IB patches at simulation startup, so num_ibs as read
         ! from the case file counts only the namelist patches. Match the global count the simulation
-        ! arrives at (s_reduce_ib_patch_array) so the IB state records can be read back.
+        ! arrives at (s_read_ib_restart_data) so the IB state records can be read back.
         do i = 1, num_particle_clouds
             num_ibs = num_ibs + particle_cloud(i)%num_particles
         end do

src/post_process/m_start_up.fpp

@@ -952,6 +952,12 @@ contains
         end if
 
         call s_mpi_bcast_user_inputs()
+
+        ! Save original BCs before decomposition overwrites them with MPI neighbor ranks
+        ib_bc_x = bc_x
+        ib_bc_y = bc_y
+        ib_bc_z = bc_z
+
         call s_initialize_parallel_io()
         output_offsets = (/offset_x, offset_y, offset_z/)
         call s_mpi_decompose_computational_domain(write_silo_ghost_offsets=format == format_silo, adjust_local_domains=.false., &

src/pre_process/m_data_output.fpp

@@ -28,7 +28,7 @@ module m_data_output
 
     private
     public :: s_write_serial_data_files, s_write_parallel_data_files, s_write_data_files, s_initialize_data_output_module, &
-        & s_finalize_data_output_module
+        & s_finalize_data_output_module, s_write_ib_state_0_file
 
     type(scalar_field), allocatable, dimension(:) :: q_cons_temp
 
@@ -733,6 +733,75 @@ contains
 
     end subroutine s_initialize_data_output_module
 
+    !> @brief Writes the initial IB layout (namelist patch_ib entries, then generated particle-cloud beds) that simulation reads
+    !! back at startup (s_read_ib_restart_data, src/simulation/m_start_up.fpp), in the layouts simulation's own IB state writers
+    !! use. Under file_per_process each rank writes only the IBs it owns to restart_data/lustre_0/ib_state_0_<rank>.dat; otherwise
+    !! only rank 0 calls this and writes every IB, in global-id order, to restart_data/ib_state_0.dat.
+    impure subroutine s_write_ib_state_0_file(glb_bounds, particle_cloud_ibs, num_particle_cloud_ibs)
+
+        type(bounds_info), dimension(3), intent(in)         :: glb_bounds
+        type(ib_patch_parameters), dimension(:), intent(in) :: particle_cloud_ibs
+        integer, intent(in)                                 :: num_particle_cloud_ibs
+        character(LEN=len_trim(case_dir) + 2*name_len)      :: file_loc
+        integer                                             :: i, ios, file_unit
+        logical, dimension(num_ibs)                         :: owned
+        real(wp), dimension(3)                              :: centroid
+
+        if (file_per_process) then
+            do i = 1, num_ibs
+                centroid = [patch_ib(i)%x_centroid, patch_ib(i)%y_centroid, patch_ib(i)%z_centroid]
+                owned(i) = f_local_rank_owns_location(centroid, glb_bounds)
+            end do
+
+            if (proc_rank == 0) call s_create_directory(trim(case_dir) // '/restart_data/lustre_0')
+            call s_mpi_barrier()
+            call s_delay_file_access(proc_rank)
+            write (file_loc, '(A,i7.7,A)') '/restart_data/lustre_0/ib_state_0_', proc_rank, '.dat'
+        else
+            owned = .true.
+            call s_create_directory(trim(case_dir) // '/restart_data')
+            file_loc = '/restart_data/ib_state_0.dat'
+        end if
+        file_loc = trim(case_dir) // trim(file_loc)
+
+        open (newunit=file_unit, file=trim(file_loc), form='unformatted', access='stream', status='replace', iostat=ios)
+        if (ios /= 0) call s_mpi_abort('Cannot open IB state output file: ' // trim(file_loc))
+
+        if (file_per_process) write (file_unit) count(owned) + num_particle_cloud_ibs
+        do i = 1, num_ibs
+            if (owned(i)) call s_write_ib_state_record(patch_ib(i), i)
+        end do
+        do i = 1, num_particle_cloud_ibs
+            call s_write_ib_state_record(particle_cloud_ibs(i), particle_cloud_ibs(i)%gbl_patch_id)
+        end do
+
+        close (file_unit)
+
+    contains
+
+        !> Writes one 20-field IB state record, prefixed by its global id under file_per_process.
+        subroutine s_write_ib_state_record(ib_patch, gbl_id)
+
+            type(ib_patch_parameters), intent(in) :: ib_patch
+            integer, intent(in)                   :: gbl_id
+            real(wp), dimension(20)               :: ib_buf
+
+            ib_buf = 0._wp
+            ib_buf(8:10) = ib_patch%vel
+            ib_buf(11:13) = ib_patch%angular_vel
+            ib_buf(14:16) = ib_patch%angles
+            ib_buf(17) = ib_patch%x_centroid
+            ib_buf(18) = ib_patch%y_centroid
+            ib_buf(19) = ib_patch%z_centroid
+            ib_buf(20) = ib_patch%radius
+
+            if (file_per_process) write (file_unit) gbl_id
+            write (file_unit) ib_buf
+
+        end subroutine s_write_ib_state_record
+
+    end subroutine s_write_ib_state_0_file
+
     !> Resets s_write_data_files pointer
     impure subroutine s_finalize_data_output_module
 

src/pre_process/m_global_parameters.fpp

@@ -52,7 +52,8 @@ module m_global_parameters
     ! Cell indices (InDices With BUFFer): includes buffer except in pre_process
     type(int_bounds_info) :: idwbuff(1:3)
     type(int_bounds_info) :: bc_x, bc_y, bc_z  !< Boundary conditions in the x-, y- and z-coordinate directions
-    type(bc_xyz_info)     :: bc                !< Combined BC storage (used by the shared beta-buffer routines; pre-process-local)
+    type(int_bounds_info) :: ib_bc_x, ib_bc_y, ib_bc_z  !< bc_x/y/z before decomposition overwrites them with MPI neighbor ranks
+    type(bc_xyz_info)     :: bc  !< Combined BC storage (used by the shared beta-buffer routines; pre-process-local)
     ! simplex_params: auto-generated in generated_decls.fpp
     ! shear_num/shear_indices/shear_BC_flip_*, bc: in m_global_parameters_common
     integer                           :: fd_order    !< Finite-difference order for CoM/probe derivative approximations
@@ -361,6 +362,27 @@ contains
             patch_ib(i)%rotation_matrix_inverse = patch_ib(i)%rotation_matrix
         end do
 
+        num_particle_clouds = 0
+        do i = 1, num_particle_clouds_max
+            particle_cloud(i)%x_centroid = 0._wp
+            particle_cloud(i)%y_centroid = 0._wp
+            particle_cloud(i)%z_centroid = 0._wp
+            particle_cloud(i)%length_x = dflt_real
+            particle_cloud(i)%length_y = dflt_real
+            particle_cloud(i)%length_z = dflt_real
+            particle_cloud(i)%num_particles = 0
+            particle_cloud(i)%radius = dflt_real
+            particle_cloud(i)%mass = dflt_real
+            particle_cloud(i)%min_spacing = 0._wp
+            particle_cloud(i)%shell_inner_radius = dflt_real
+            particle_cloud(i)%shell_outer_radius = dflt_real
+            particle_cloud(i)%moving_ibm = 0
+            particle_cloud(i)%seed = 0
+            particle_cloud(i)%cloud_geometry = 1
+            particle_cloud(i)%packing_method = dflt_int
+            particle_cloud(i)%periodic = 0
+        end do
+
         do i = 1, num_ib_airfoils_max
             ib_airfoil(i)%c = dflt_real
             ib_airfoil(i)%p = dflt_real

src/pre_process/m_mpi_proxy.fpp

@@ -135,6 +135,20 @@ contains
             call MPI_BCAST(patch_ib(i)%slip, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
         end do
 
+        ! manual: particle_cloud (runtime loop to num_particle_clouds; irregular member subset)
+        do i = 1, num_particle_clouds
+            #:for VAR in ['x_centroid', 'y_centroid', 'z_centroid', 'length_x', 'length_y', 'length_z', &
+                & 'radius', 'mass', 'min_spacing', 'shell_inner_radius', 'shell_outer_radius']
+                call MPI_BCAST(particle_cloud(i)%${VAR}$, 1, mpi_p, 0, MPI_COMM_WORLD, ierr)
+            #:endfor
+            call MPI_BCAST(particle_cloud(i)%num_particles, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
+            call MPI_BCAST(particle_cloud(i)%moving_ibm, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
+            call MPI_BCAST(particle_cloud(i)%seed, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
+            call MPI_BCAST(particle_cloud(i)%cloud_geometry, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
+            call MPI_BCAST(particle_cloud(i)%packing_method, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
+            call MPI_BCAST(particle_cloud(i)%periodic, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
+        end do
+
         ! manual: ib_airfoil (kept manual alongside patch_ib)
         do i = 1, num_ib_airfoils_max
             #:for VAR in ['c', 'p', 't', 'm']