Initial move of code from simulation to pre_process

Author: danieljvickers

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
 
     type(scalar_field), allocatable, dimension(:) :: q_cons_temp
 
@@ -733,6 +733,50 @@ contains
 
     end subroutine s_initialize_data_output_module
 
+    !> @brief Writes restart_data/ib_state_0.dat: the initial IB layout (namelist patch_ib entries, then any generated
+    !! particle-cloud beds). Rank-0-only. Read back by simulation at startup via s_read_ib_restart_data(0, ...)
+    !! (src/simulation/m_start_up.fpp). Uses the same 20-field record layout s_write_serial_ib_state
+    !! (src/simulation/m_data_output.fpp) writes, so simulation's own writers can freely overwrite this file later at t_step_start
+    !! == 0 without changing its layout - only position (fields 17:19) and radius (field 20) are populated here; everything else
+    !! (time, force, torque, vel, angular_vel, angles) is zero for a freshly generated IB.
+    impure subroutine s_write_ib_state_0(particle_cloud_ibs, num_particle_cloud_ibs)
+
+        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
+        integer, parameter                                  :: NFIELDS_PER_IB = 20
+        real(wp)                                            :: ib_buf(NFIELDS_PER_IB)
+
+        call s_create_directory(trim(case_dir) // '/restart_data')
+
+        file_loc = trim(case_dir) // '/restart_data/ib_state_0.dat'
+
+        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))
+
+        ib_buf = 0._wp
+
+        do i = 1, num_ibs
+            ib_buf(17) = patch_ib(i)%x_centroid
+            ib_buf(18) = patch_ib(i)%y_centroid
+            ib_buf(19) = patch_ib(i)%z_centroid
+            ib_buf(20) = patch_ib(i)%radius
+            write (file_unit) ib_buf
+        end do
+
+        do i = 1, num_particle_cloud_ibs
+            ib_buf(17) = particle_cloud_ibs(i)%x_centroid
+            ib_buf(18) = particle_cloud_ibs(i)%y_centroid
+            ib_buf(19) = particle_cloud_ibs(i)%z_centroid
+            ib_buf(20) = particle_cloud_ibs(i)%radius
+            write (file_unit) ib_buf
+        end do
+
+        close (file_unit)
+
+    end subroutine s_write_ib_state_0
+
     !> Resets s_write_data_files pointer
     impure subroutine s_finalize_data_output_module
 

src/pre_process/m_global_parameters.fpp

@@ -361,6 +361,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/simulation/m_particle_cloud.fpp src/pre_process/m_particle_cloud.fppsrc/simulation/m_particle_cloud.fpp renamed to src/pre_process/m_particle_cloud.fpp

@@ -1,18 +1,19 @@
 !>
 !! @file m_particle_cloud.fpp
 !! @brief Generates particle beds: converts particle_cloud specifications into
-!!        individual sphere/circle particle_cloud_ibs entries before reduction.
+!!        individual sphere/circle particle_cloud_ibs entries, written to
+!!        restart_data/ib_state_0.dat for simulation to read at startup.
 
 #:include 'macros.fpp'
 
-!> @brief Generates particle beds by converting particle_cloud patch specifications into individual immersed boundary patches before
-!! domain reduction. Each rank runs the same deterministic placement so no MPI broadcast of particle positions is needed.
+!> @brief Generates particle beds by converting particle_cloud patch specifications into individual immersed boundary patches. Runs
+!! once, on rank 0 only, producing the full global (unfiltered) particle list - simulation applies its own neighborhood filtering
+!! when it reads the result back from restart_data/ib_state_0.dat.
 module m_particle_cloud
 
     use m_global_parameters
     use m_constants
     use m_mpi_common
-    use m_collisions
 
     implicit none
 
@@ -22,11 +23,10 @@ module m_particle_cloud
 
 contains
 
-    !> Generate all particle beds and fill particle_cloud_ibs. Called on all ranks before s_reduce_ib_patch_array. Each packing
-    !! method owns and allocates its own per-cloud working array (see s_particle_cloud_lattice / s_particle_cloud_rejection_pack)
-    !! and hands back only the entries that fall within this rank's IB neighborhood. Only the first num_particle_cloud_ibs of them
-    !! are actually written - callers must use that count, not size(particle_cloud_ibs), since the remainder of the array is left
-    !! uninitialized.
+    !> Generate all particle beds and fill particle_cloud_ibs with the full, global (unfiltered) particle list. Called once, on rank
+    !! 0 only. Each packing method owns and allocates its own per-cloud working array (see s_particle_cloud_lattice /
+    !! s_particle_cloud_rejection_pack). Only the first num_particle_cloud_ibs of them are actually written - callers must use that
+    !! count, not size(particle_cloud_ibs), since the remainder of the array is left uninitialized.
     impure subroutine s_generate_particle_clouds(particle_cloud_ibs, num_particle_cloud_ibs)
 
         type(ib_patch_parameters), allocatable, intent(out), dimension(:) :: particle_cloud_ibs
@@ -78,12 +78,10 @@ contains
 
     end subroutine s_generate_particle_clouds
 
-    !> Rejection-samples particle centres into a box or hemisphere-shell region with a minimum centre-to-centre spacing. Rejection
-    !! sampling needs every placed particle tracked (regardless of which rank's neighborhood it falls in) to detect overlaps
-    !! deterministically, so cloud_ibs is allocated here to the cloud's full requested particle count and only pared down to this
-    !! rank's neighborhood afterwards, via s_reduce_particle_cloud_ibs. Only the per-candidate geometry sampling differs between box
+    !> Rejection-samples particle centres into a box or hemisphere-shell region with a minimum centre-to-centre spacing. cloud_ibs
+    !! is allocated here to the cloud's full requested particle count. Only the per-candidate geometry sampling differs between box
     !! and hemisphere shell; it is delegated to s_sample_cloud_candidate, and every other step (overlap rejection via the spatial
-    !! hash, acceptance, reduction) is geometry-independent.
+    !! hash, acceptance) is geometry-independent.
     subroutine s_particle_cloud_rejection_pack(cloud_idx, glbl_idx, cloud_ibs, num_cloud_ibs)
 
         integer, intent(in)                                               :: cloud_idx
@@ -174,7 +172,6 @@ contains
 
         deallocate (placed, hash_head, chain_next)
 
-        call s_reduce_particle_cloud_ibs(cloud_ibs, ib_idx)
         num_cloud_ibs = ib_idx
 
     end subroutine s_particle_cloud_rejection_pack
@@ -243,9 +240,7 @@ contains
     !> Places particles on the optimally dense lattice for the cloud region: a triangular lattice in 2D, a face-centered cubic
     !! lattice in 3D. The lattice spacing is set by the particle density (num_particles over the region area/volume); if that
     !! spacing falls below the required centre-to-centre distance (2*radius + min_spacing), the region is too dense and the run is
-    !! aborted. No two lattice sites can overlap, so unlike rejection packing each site's IB neighborhood membership
-    !! (get_neighbor_bounds() must already have run) is checked as it is generated and only in-neighborhood sites are stored;
-    !! cloud_ibs is therefore allocated to the neighborhood-sized cap rather than the cloud's full particle count.
+    !! aborted. No two lattice sites can overlap, so cloud_ibs is allocated to the cloud's full particle count.
     subroutine s_particle_cloud_lattice(cloud_idx, glbl_idx, cloud_ibs, num_cloud_ibs)
 
         integer, intent(in)                                               :: cloud_idx
@@ -259,7 +254,7 @@ contains
         real(wp), dimension(4)                                            :: bx_off, by_off, bz_off
         real(wp), dimension(3)                                            :: centroid
 
-        allocate (cloud_ibs(min(num_ib_patches_max_namelist, particle_cloud(cloud_idx)%num_particles)))
+        allocate (cloud_ibs(particle_cloud(cloud_idx)%num_particles))
         ib_idx = 0
 
         xmin = particle_cloud(cloud_idx)%x_centroid - 0.5_wp*particle_cloud(cloud_idx)%length_x
@@ -301,10 +296,7 @@ contains
                 do while (px <= xmax .and. n_placed < n_target)
                     glbl_idx = glbl_idx + 1
                     centroid = [px, py, particle_cloud(cloud_idx)%z_centroid]
-                    if (f_neighborhood_ranks_own_location(centroid)) then
-                        call s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, centroid(1), centroid(2), centroid(3), &
-                                                  & cloud_ibs)
-                    end if
+                    call s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, centroid(1), centroid(2), centroid(3), cloud_ibs)
                     n_placed = n_placed + 1
                     col = col + 1
                     px = x0 + real(col, wp)*spacing
@@ -328,10 +320,8 @@ contains
                             centroid = [xmin + real(ix, wp)*cell + bx_off(b), ymin + real(jy, wp)*cell + by_off(b), &
                                                     & zmin + real(kz, wp)*cell + bz_off(b)]
                             glbl_idx = glbl_idx + 1
-                            if (f_neighborhood_ranks_own_location(centroid)) then
-                                call s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, centroid(1), centroid(2), &
-                                                          & centroid(3), cloud_ibs)
-                            end if
+                            call s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, centroid(1), centroid(2), centroid(3), &
+                                                      & cloud_ibs)
                             n_placed = n_placed + 1
                         end do
                     end do
@@ -344,11 +334,11 @@ contains
 
     end subroutine s_particle_cloud_lattice
 
-    !> Writes a single placed particle into particle_cloud_ibs at the next free slot, advancing ib_idx. The caller decides whether
-    !! this particle belongs in the array (neighborhood membership, for lattice packing, or unconditionally for rejection packing -
-    !! see s_particle_cloud_lattice / s_particle_cloud_rejection_pack) and supplies its already-assigned, absolute global patch id
-    !! via glbl_idx - s_reduce_ib_patch_array copies gbl_patch_id as-is. Shared by all packing methods so the per-particle
-    !! ib_patch_parameters setup stays in one place.
+    !> Writes a single placed particle into particle_cloud_ibs at the next free slot, advancing ib_idx, tagged with its
+    !! already-assigned, absolute global patch id via glbl_idx. Shared by all packing methods so the per-particle
+    !! ib_patch_parameters setup stays in one place. Only x/y/z_centroid and radius are ever read back out of this struct (by
+    !! s_write_ib_state_0); the rest is filled in for completeness and because leaving fields uninitialized here would otherwise
+    !! surface as garbage if this struct is ever read further downstream.
     subroutine s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, px, py, pz, particle_cloud_ibs)
 
         integer, intent(in)                                    :: cloud_idx, glbl_idx, geom
@@ -358,7 +348,7 @@ contains
 
         ib_idx = ib_idx + 1
         @:PROHIBIT(ib_idx > size(particle_cloud_ibs), &
-                   & "Too many particle-cloud IBs in one rank's neighborhood. Modify case file or increase num_ib_patches_max_namelist.")
+                   & "Too many particle-cloud IBs. Modify case file or increase num_ib_patches_max_namelist.")
 
         particle_cloud_ibs(ib_idx)%gbl_patch_id = glbl_idx
         particle_cloud_ibs(ib_idx)%geometry = geom
@@ -388,12 +378,9 @@ contains
         particle_cloud_ibs(ib_idx)%moving_ibm = particle_cloud(cloud_idx)%moving_ibm
         particle_cloud_ibs(ib_idx)%slip = .false.
 
-        ! Particles are inert surfaces. These must be set explicitly: particle_cloud_ibs is
-        ! allocated (not default-initialized) and s_reduce_ib_patch_array copies the whole
-        ! struct into patch_ib, overwriting the defaults from
-        ! s_assign_default_values_to_user_inputs -- so anything left unset here reaches the
-        ! solver as uninitialized memory (a nonzero v_blow injects a garbage wall-normal
-        ! velocity and NaNs the field).
+        ! Particles are inert surfaces. particle_cloud_ibs is allocated (not default-initialized),
+        ! so these must be set explicitly even though only x/y/z_centroid and radius are read back
+        ! out of this struct downstream.
         particle_cloud_ibs(ib_idx)%v_blow = 0._wp
         particle_cloud_ibs(ib_idx)%inj_species = 0
         particle_cloud_ibs(ib_idx)%burn_rate_exp = 0._wp
@@ -500,29 +487,6 @@ contains
 
     end subroutine s_check_cloud_particle_overlap
 
-    !> Compacts cloud_ibs(1:num_ibs) in place, discarding entries outside this rank's IB neighborhood (get_neighbor_bounds() must
-    !! already have run) and updating num_ibs to the retained count. Used by rejection packing, which cannot filter as it places
-    !! particles (see s_particle_cloud_rejection_pack), to pare its full, unfiltered placement down to this rank's neighborhood.
-    subroutine s_reduce_particle_cloud_ibs(cloud_ibs, num_cloud_ibs)
-
-        type(ib_patch_parameters), intent(inout), dimension(:) :: cloud_ibs
-        integer, intent(inout)                                 :: num_cloud_ibs
-        integer                                                :: i, write_idx
-        real(wp), dimension(3)                                 :: centroid
-
-        write_idx = 0
-        do i = 1, num_cloud_ibs
-            centroid = [cloud_ibs(i)%x_centroid, cloud_ibs(i)%y_centroid, 0._wp]
-            if (num_dims == 3) centroid(3) = cloud_ibs(i)%z_centroid
-            if (f_neighborhood_ranks_own_location(centroid)) then
-                write_idx = write_idx + 1
-                if (write_idx /= i) cloud_ibs(write_idx) = cloud_ibs(i)
-            end if
-        end do
-        num_cloud_ibs = write_idx
-
-    end subroutine s_reduce_particle_cloud_ibs
-
     !> Xorshift PRNG. Advances seed in-place and returns a value in [0, 1).
     function f_xorshift(seed) result(rval)
 

src/pre_process/m_start_up.fpp

@@ -28,6 +28,7 @@ module m_start_up
 
     use m_check_patches
     use m_check_ib_patches
+    use m_particle_cloud
     use m_helper
     use m_checker_common
     use m_checker
@@ -133,7 +134,21 @@ contains
 
         call s_check_patches()
 
-        if (ib) call s_check_ib_patches()
+        if (ib) then
+            call s_check_ib_patches()
+
+            if (proc_rank == 0) then
+                block
+                    type(ib_patch_parameters), allocatable :: particle_cloud_ibs(:)
+                    integer                                :: num_particle_cloud_ibs
+
+                    call s_generate_particle_clouds(particle_cloud_ibs, num_particle_cloud_ibs)
+                    call s_write_ib_state_0(particle_cloud_ibs, num_particle_cloud_ibs)
+                    deallocate (particle_cloud_ibs)
+                end block
+            end if
+            call s_mpi_barrier()
+        end if
 
     end subroutine s_check_input_file