From 9c1904e8edb5c595c1c2399718312e131dcb744e Mon Sep 17 00:00:00 2001 From: danieljvickers Date: Thu, 4 Jun 2026 21:52:00 -0400 Subject: [PATCH 1/4] Renamed particle beds to particle clouds --- docs/documentation/case.md | 2 +- docs/module_categories.json | 2 +- .../case.py | 27 ++-- src/common/m_constants.fpp | 2 +- src/common/m_derived_types.fpp | 5 +- src/simulation/m_checker.fpp | 21 +++ src/simulation/m_global_parameters.fpp | 31 +++-- src/simulation/m_mpi_proxy.fpp | 13 +- ..._particle_bed.fpp => m_particle_cloud.fpp} | 128 +++++++++--------- src/simulation/m_start_up.fpp | 40 +++--- .../golden-metadata.txt | 0 tests/{BA186DDF => A7E7A245}/golden.txt | 0 toolchain/mfc/case_validator.py | 8 +- toolchain/mfc/params/definitions.py | 13 +- toolchain/mfc/params/descriptions.py | 2 +- 15 files changed, 160 insertions(+), 134 deletions(-) rename examples/{2D_mibm_particle_bed => 2D_mibm_particle_cloud}/case.py (85%) rename src/simulation/{m_particle_bed.fpp => m_particle_cloud.fpp} (56%) rename tests/{BA186DDF => A7E7A245}/golden-metadata.txt (100%) rename tests/{BA186DDF => A7E7A245}/golden.txt (100%) diff --git a/docs/documentation/case.md b/docs/documentation/case.md index 3fd9c143ef..dd5f2692de 100644 --- a/docs/documentation/case.md +++ b/docs/documentation/case.md @@ -313,7 +313,7 @@ This is enabled by adding ``'elliptic_smoothing': "T",`` and ``'elliptic_smoothi | ---: | :----: | :--- | | `num_ibs` | Integer | Number of immersed boundary patches | | `num_stl_models` | Integer | Number of STL/OBJ model entries in the `stl_models` array | -| `num_particle_beds` | Integer | Number of particle bed specifications to generate immersed boundary patches from | +| `num_particle_clouds` | Integer | Number of particle bed specifications to generate immersed boundary patches from | | `ib_neighborhood_radius` | Integer | Parameter that controls the neighborhood size for IB detection. | | `geometry` | Integer | Geometry configuration of the patch.| | `x[y,z]_centroid` | Real | Centroid of the applied geometry in the [x,y,z]-direction. | diff --git a/docs/module_categories.json b/docs/module_categories.json index 21de015c71..7cf0ccc995 100644 --- a/docs/module_categories.json +++ b/docs/module_categories.json @@ -38,7 +38,7 @@ "m_compute_cbc", "m_boundary_common", "m_ibm", - "m_particle_bed", + "m_particle_cloud", "m_igr", "m_ib_patches", "m_compute_levelset" diff --git a/examples/2D_mibm_particle_bed/case.py b/examples/2D_mibm_particle_cloud/case.py similarity index 85% rename from examples/2D_mibm_particle_bed/case.py rename to examples/2D_mibm_particle_cloud/case.py index 5db1e707c7..68abf5badc 100644 --- a/examples/2D_mibm_particle_bed/case.py +++ b/examples/2D_mibm_particle_cloud/case.py @@ -85,19 +85,20 @@ "collision_time": collision_time, "ib_coefficient_of_friction": 0.1, # Particle bed: 20 free-floating circles placed randomly in region - "num_particle_beds": 1, - "particle_bed(1)%x_centroid": bed_x, - "particle_bed(1)%y_centroid": bed_y, - "particle_bed(1)%z_centroid": 0.0, - "particle_bed(1)%length_x": bed_lx, - "particle_bed(1)%length_y": bed_ly, - "particle_bed(1)%length_z": 0.0, - "particle_bed(1)%num_particles": 20, - "particle_bed(1)%radius": particle_radius, - "particle_bed(1)%mass": particle_mass, - "particle_bed(1)%min_spacing": particle_min_spacing, - "particle_bed(1)%moving_ibm": 2, - "particle_bed(1)%seed": 42, + "num_particle_clouds": 1, + "particle_cloud(1)%x_centroid": bed_x, + "particle_cloud(1)%y_centroid": bed_y, + "particle_cloud(1)%z_centroid": 0.0, + "particle_cloud(1)%length_x": bed_lx, + "particle_cloud(1)%length_y": bed_ly, + "particle_cloud(1)%length_z": 0.0, + "particle_cloud(1)%num_particles": 20, + "particle_cloud(1)%radius": particle_radius, + "particle_cloud(1)%mass": particle_mass, + "particle_cloud(1)%min_spacing": particle_min_spacing, + "particle_cloud(1)%moving_ibm": 2, + "particle_cloud(1)%seed": 42, + "particle_cloud(1)%packing_method": 1, # Output "format": 1, "precision": 2, diff --git a/src/common/m_constants.fpp b/src/common/m_constants.fpp index 8ed324f363..b2259ae44d 100644 --- a/src/common/m_constants.fpp +++ b/src/common/m_constants.fpp @@ -31,7 +31,7 @@ module m_constants !> Fixed capacity of patch_ib (namelist patches + local particle bed subset after reduction) integer, parameter :: num_ib_patches_max_namelist = 54000 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_particle_clouds_max = 10 !< Maximum number of particle bed patch specifications integer, parameter :: num_bc_patches_max = 10 !< Maximum number of boundary condition patches integer, parameter :: max_2d_fourier_modes = 10 !< Max Fourier mode index for 2D modal patch (geometry 13) integer, parameter :: max_sph_harm_degree = 5 !< Max degree L for 3D spherical harmonic patch (geometry 14) diff --git a/src/common/m_derived_types.fpp b/src/common/m_derived_types.fpp index d2ce0189f1..5dd132a6bc 100644 --- a/src/common/m_derived_types.fpp +++ b/src/common/m_derived_types.fpp @@ -320,7 +320,7 @@ module m_derived_types real(wp), dimension(1:3) :: step_angular_vel !< velocity array used to store intermediate steps in the time_stepper module end type ib_patch_parameters - type particle_bed_parameters + type particle_cloud_parameters real(wp) :: x_centroid, y_centroid, z_centroid !< Center of the particle bed region real(wp) :: length_x, length_y, length_z !< Dimensions of the particle bed region integer :: num_particles !< Number of particles to generate @@ -329,7 +329,8 @@ module m_derived_types real(wp) :: min_spacing !< Minimum surface-to-surface gap (particle centers are 2*radius + min_spacing apart) integer :: moving_ibm !< Motion flag: 0=static, 1=moving (forces), 2=forced path integer :: seed !< Random seed for reproducible placement - end type particle_bed_parameters + integer :: packing_method !< Packing algorithm: 1=rejection sampling + end type particle_cloud_parameters !> Derived type annexing the physical parameters (PP) of the fluids. These include the specific heat ratio function and liquid !! stiffness function. diff --git a/src/simulation/m_checker.fpp b/src/simulation/m_checker.fpp index 4cce28299b..2ed8aca3c3 100644 --- a/src/simulation/m_checker.fpp +++ b/src/simulation/m_checker.fpp @@ -38,6 +38,10 @@ contains @:PROHIBIT(ib_state_wrt .and. .not. ib, "ib_state_wrt requires ib to be enabled") + if (num_particle_clouds > 0) then + call s_check_inputs_particle_clouds + end if + end subroutine s_check_inputs !> Checks constraints on compiler options @@ -106,4 +110,21 @@ contains end subroutine s_check_inputs_nvidia_uvm + !> Checks that each active particle cloud has a valid packing_method specified + impure subroutine s_check_inputs_particle_clouds + + integer :: i + character(len=5) :: idxStr + + do i = 1, num_particle_clouds + call s_int_to_str(i, idxStr) + @:PROHIBIT(particle_cloud(i)%packing_method == dflt_int, & + & "particle_cloud("//trim(idxStr)//")%packing_method must be specified (1 = rejection sampling)") + @:PROHIBIT(particle_cloud(i)%packing_method /= 1, & + & "particle_cloud("//trim(idxStr) & + & //")%packing_method must be 1 (rejection sampling is the only supported method)") + end do + + end subroutine s_check_inputs_particle_clouds + end module m_checker diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index f4a886fbfb..1bd6ea26f8 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/src/simulation/m_global_parameters.fpp @@ -263,7 +263,7 @@ module m_global_parameters !> @{ type(ib_patch_parameters), dimension(num_ib_patches_max_namelist) :: patch_ib !< Immersed boundary patch parameters integer, dimension(num_local_ibs_max) :: local_ib_patch_ids !< lookup table of IBs in the local compute domain - type(particle_bed_parameters), dimension(num_particle_beds_max) :: particle_bed !< Particle bed specifications + type(particle_cloud_parameters), dimension(num_particle_clouds_max) :: particle_cloud !< Particle bed specifications integer, allocatable, dimension(:,:,:) :: ib_neighbor_ranks !< MPI ranks of neighborhood domains, indexed (-N:N,-N:N,-N:N) type(ib_airfoil_parameters), dimension(num_ib_airfoils_max) :: ib_airfoil !< Per-airfoil NACA user inputs (namelist) type(ib_airfoil_grid), dimension(num_ib_airfoils_max) :: ib_airfoil_grids !< Per-airfoil computed surface grids @@ -689,20 +689,21 @@ contains ib_airfoil_grids(i)%Np = 0 end do - num_particle_beds = 0 - do i = 1, num_particle_beds_max - particle_bed(i)%x_centroid = 0._wp - particle_bed(i)%y_centroid = 0._wp - particle_bed(i)%z_centroid = 0._wp - particle_bed(i)%length_x = dflt_real - particle_bed(i)%length_y = dflt_real - particle_bed(i)%length_z = dflt_real - particle_bed(i)%num_particles = 0 - particle_bed(i)%radius = dflt_real - particle_bed(i)%mass = dflt_real - particle_bed(i)%min_spacing = 0._wp - particle_bed(i)%moving_ibm = 0 - particle_bed(i)%seed = 0 + 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)%moving_ibm = 0 + particle_cloud(i)%seed = 0 + particle_cloud(i)%packing_method = dflt_int end do do i = 1, num_ib_patches_max_namelist diff --git a/src/simulation/m_mpi_proxy.fpp b/src/simulation/m_mpi_proxy.fpp index f754cd8d26..ae955cef55 100644 --- a/src/simulation/m_mpi_proxy.fpp +++ b/src/simulation/m_mpi_proxy.fpp @@ -74,7 +74,7 @@ contains & 'wave_speeds', 'avg_state', 'precision', 'bc_x%beg', 'bc_x%end', & & 'bc_y%beg', 'bc_y%end', 'bc_z%beg', 'bc_z%end', 'fd_order', & & 'num_probes', 'num_integrals', 'bubble_model', 'thermal', & - & 'num_source', 'relax_model', 'num_ibs', 'num_particle_beds', 'n_start', & + & 'num_source', 'relax_model', 'num_ibs', 'num_particle_clouds', 'n_start', & & 'num_bc_patches', 'num_igr_iters', 'num_igr_warm_start_iters', & & 'adap_dt_max_iters', 'collision_model', 'ib_neighborhood_radius', & & 'int_comp' ] @@ -222,14 +222,15 @@ contains #:endfor end do - do i = 1, num_particle_beds + 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'] - call MPI_BCAST(particle_bed(i)%${VAR}$, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(particle_cloud(i)%${VAR}$, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) #:endfor - call MPI_BCAST(particle_bed(i)%num_particles, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) - call MPI_BCAST(particle_bed(i)%moving_ibm, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) - call MPI_BCAST(particle_bed(i)%seed, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + 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)%packing_method, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) end do do j = 1, num_probes_max diff --git a/src/simulation/m_particle_bed.fpp b/src/simulation/m_particle_cloud.fpp similarity index 56% rename from src/simulation/m_particle_bed.fpp rename to src/simulation/m_particle_cloud.fpp index 7c30e74fa5..b50e188510 100644 --- a/src/simulation/m_particle_bed.fpp +++ b/src/simulation/m_particle_cloud.fpp @@ -1,11 +1,11 @@ !> -!! @file m_particle_bed.fpp -!! @brief Generates particle beds: converts particle_bed specifications into -!! individual sphere/circle particle_bed_ibs entries before reduction. +!! @file m_particle_cloud.fpp +!! @brief Generates particle beds: converts particle_cloud specifications into +!! individual sphere/circle particle_cloud_ibs entries before reduction. -!> @brief Generates particle beds by converting particle_bed patch specifications into individual immersed boundary patches before +!> @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. -module m_particle_bed +module m_particle_cloud use m_global_parameters use m_constants @@ -15,16 +15,16 @@ module m_particle_bed private - public :: s_generate_particle_beds + public :: s_generate_particle_clouds contains - !> Generate all particle beds and fill particle_bed_ibs. Called on all ranks before s_reduce_ib_patch_array. Uses a spatial hash - !! grid (cell size = min_dist) so each candidate requires only 3^dim distance checks on average instead of O(n). The placement - !! is fully deterministic given the per-bed seed, so all ranks produce an identical result without MPI. - impure subroutine s_generate_particle_beds(particle_bed_ibs) + !> Generate all particle beds and fill particle_cloud_ibs. Called on all ranks before s_reduce_ib_patch_array. Uses a spatial + !! hash grid (cell size = min_dist) so each candidate requires only 3^dim distance checks on average instead of O(n). The + !! placement is fully deterministic given the per-bed seed, so all ranks produce an identical result without MPI. + impure subroutine s_generate_particle_clouds(particle_cloud_ibs) - type(ib_patch_parameters), allocatable, intent(out), dimension(:) :: particle_bed_ibs + type(ib_patch_parameters), allocatable, intent(out), dimension(:) :: particle_cloud_ibs integer :: b, ib_idx, geom integer :: n_placed, n_total_placed, n_total_particles integer(8) :: n_attempts, max_attempts @@ -41,32 +41,32 @@ contains integer :: dx_b, dy_b, dz_b, dz_lo, dz_hi, j integer, allocatable :: hash_head(:), chain_next(:) - if (num_particle_beds == 0) then - allocate (particle_bed_ibs(0)) + if (num_particle_clouds == 0) then + allocate (particle_cloud_ibs(0)) return end if call cpu_time(t_start) n_total_placed = 0 - ! Pre-count total particles across all beds so particle_bed_ibs can be allocated exactly once. + ! Pre-count total particles across all beds so particle_cloud_ibs can be allocated exactly once. n_total_particles = 0 - do b = 1, num_particle_beds - n_total_particles = n_total_particles + particle_bed(b)%num_particles + do b = 1, num_particle_clouds + n_total_particles = n_total_particles + particle_cloud(b)%num_particles end do - allocate (particle_bed_ibs(n_total_particles)) + allocate (particle_cloud_ibs(n_total_particles)) - ib_idx = 0 ! index into particle_bed_ibs + ib_idx = 0 ! index into particle_cloud_ibs - 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 - ymin = particle_bed(b)%y_centroid - 0.5_wp*particle_bed(b)%length_y - ymax = particle_bed(b)%y_centroid + 0.5_wp*particle_bed(b)%length_y - zmin = particle_bed(b)%z_centroid - 0.5_wp*particle_bed(b)%length_z - zmax = particle_bed(b)%z_centroid + 0.5_wp*particle_bed(b)%length_z + do b = 1, num_particle_clouds + xmin = particle_cloud(b)%x_centroid - 0.5_wp*particle_cloud(b)%length_x + xmax = particle_cloud(b)%x_centroid + 0.5_wp*particle_cloud(b)%length_x + ymin = particle_cloud(b)%y_centroid - 0.5_wp*particle_cloud(b)%length_y + ymax = particle_cloud(b)%y_centroid + 0.5_wp*particle_cloud(b)%length_y + zmin = particle_cloud(b)%z_centroid - 0.5_wp*particle_cloud(b)%length_z + zmax = particle_cloud(b)%z_centroid + 0.5_wp*particle_cloud(b)%length_z - min_dist = 2._wp*particle_bed(b)%radius + particle_bed(b)%min_spacing + min_dist = 2._wp*particle_cloud(b)%radius + particle_cloud(b)%min_spacing if (p == 0) then geom = 2 ! circle for 2D @@ -78,29 +78,29 @@ contains dz_hi = 1 end if - max_attempts = int(particle_bed(b)%num_particles, 8)*1000_8 + max_attempts = int(particle_cloud(b)%num_particles, 8)*1000_8 n_placed = 0 n_attempts = 0 - seed = particle_bed(b)%seed + seed = particle_cloud(b)%seed if (seed == 0) seed = 1 + b*1013904223 - allocate (placed(3, particle_bed(b)%num_particles)) + allocate (placed(3, particle_cloud(b)%num_particles)) ! Hash table: 4x overprovisioned for ~25% load factor, minimum 16 buckets. chain_next(i) links placed particle i to the ! previous occupant of its bucket. - hash_size = max(16, 4*particle_bed(b)%num_particles) + hash_size = max(16, 4*particle_cloud(b)%num_particles) allocate (hash_head(hash_size)) - allocate (chain_next(particle_bed(b)%num_particles)) + allocate (chain_next(particle_cloud(b)%num_particles)) hash_head = -1 chain_next = -1 - do while (n_placed < particle_bed(b)%num_particles .and. n_attempts < max_attempts) + do while (n_placed < particle_cloud(b)%num_particles .and. n_attempts < max_attempts) n_attempts = n_attempts + 1 rx = xmin + f_xorshift(seed)*(xmax - xmin) ry = ymin + f_xorshift(seed)*(ymax - ymin) if (p == 0) then - rz = particle_bed(b)%z_centroid + rz = particle_cloud(b)%z_centroid else rz = zmin + f_xorshift(seed)*(zmax - zmin) end if @@ -149,38 +149,38 @@ contains ib_idx = ib_idx + 1 - ! gbl_patch_id is relative within particle_bed_ibs here; s_reduce_ib_patch_array adjusts to global indexing. - particle_bed_ibs(ib_idx)%gbl_patch_id = ib_idx - particle_bed_ibs(ib_idx)%geometry = geom - particle_bed_ibs(ib_idx)%x_centroid = rx - particle_bed_ibs(ib_idx)%y_centroid = ry - particle_bed_ibs(ib_idx)%z_centroid = rz - particle_bed_ibs(ib_idx)%step_x_centroid = rx - particle_bed_ibs(ib_idx)%step_y_centroid = ry - particle_bed_ibs(ib_idx)%step_z_centroid = rz - particle_bed_ibs(ib_idx)%angles(:) = 0._wp - particle_bed_ibs(ib_idx)%step_angles(:) = 0._wp - particle_bed_ibs(ib_idx)%vel(:) = 0._wp - particle_bed_ibs(ib_idx)%step_vel(:) = 0._wp - particle_bed_ibs(ib_idx)%angular_vel(:) = 0._wp - particle_bed_ibs(ib_idx)%step_angular_vel(:) = 0._wp - particle_bed_ibs(ib_idx)%force(:) = 0._wp - particle_bed_ibs(ib_idx)%torque(:) = 0._wp - particle_bed_ibs(ib_idx)%centroid_offset(:) = 0._wp - particle_bed_ibs(ib_idx)%rotation_matrix = 0._wp - particle_bed_ibs(ib_idx)%rotation_matrix(1, 1) = 1._wp - particle_bed_ibs(ib_idx)%rotation_matrix(2, 2) = 1._wp - particle_bed_ibs(ib_idx)%rotation_matrix(3, 3) = 1._wp - particle_bed_ibs(ib_idx)%rotation_matrix_inverse = particle_bed_ibs(ib_idx)%rotation_matrix - particle_bed_ibs(ib_idx)%radius = particle_bed(b)%radius - particle_bed_ibs(ib_idx)%mass = particle_bed(b)%mass - particle_bed_ibs(ib_idx)%moment = dflt_real - particle_bed_ibs(ib_idx)%moving_ibm = particle_bed(b)%moving_ibm - particle_bed_ibs(ib_idx)%slip = .false. + ! gbl_patch_id is relative within particle_cloud_ibs here; s_reduce_ib_patch_array adjusts to global indexing. + particle_cloud_ibs(ib_idx)%gbl_patch_id = ib_idx + particle_cloud_ibs(ib_idx)%geometry = geom + particle_cloud_ibs(ib_idx)%x_centroid = rx + particle_cloud_ibs(ib_idx)%y_centroid = ry + particle_cloud_ibs(ib_idx)%z_centroid = rz + particle_cloud_ibs(ib_idx)%step_x_centroid = rx + particle_cloud_ibs(ib_idx)%step_y_centroid = ry + particle_cloud_ibs(ib_idx)%step_z_centroid = rz + particle_cloud_ibs(ib_idx)%angles(:) = 0._wp + particle_cloud_ibs(ib_idx)%step_angles(:) = 0._wp + particle_cloud_ibs(ib_idx)%vel(:) = 0._wp + particle_cloud_ibs(ib_idx)%step_vel(:) = 0._wp + particle_cloud_ibs(ib_idx)%angular_vel(:) = 0._wp + particle_cloud_ibs(ib_idx)%step_angular_vel(:) = 0._wp + particle_cloud_ibs(ib_idx)%force(:) = 0._wp + particle_cloud_ibs(ib_idx)%torque(:) = 0._wp + particle_cloud_ibs(ib_idx)%centroid_offset(:) = 0._wp + particle_cloud_ibs(ib_idx)%rotation_matrix = 0._wp + particle_cloud_ibs(ib_idx)%rotation_matrix(1, 1) = 1._wp + particle_cloud_ibs(ib_idx)%rotation_matrix(2, 2) = 1._wp + particle_cloud_ibs(ib_idx)%rotation_matrix(3, 3) = 1._wp + particle_cloud_ibs(ib_idx)%rotation_matrix_inverse = particle_cloud_ibs(ib_idx)%rotation_matrix + particle_cloud_ibs(ib_idx)%radius = particle_cloud(b)%radius + particle_cloud_ibs(ib_idx)%mass = particle_cloud(b)%mass + particle_cloud_ibs(ib_idx)%moment = dflt_real + particle_cloud_ibs(ib_idx)%moving_ibm = particle_cloud(b)%moving_ibm + particle_cloud_ibs(ib_idx)%slip = .false. end if end do - if (n_placed < particle_bed(b)%num_particles) then + if (n_placed < particle_cloud(b)%num_particles) then call s_mpi_abort("Error :: Failed to place all particles in particle bed") end if @@ -192,7 +192,7 @@ contains if (proc_rank == 0) print '(a,i0,a,f0.3,a)', 'Particle beds placed ', n_total_placed, ' particles in ', t_end - t_start, & & ' seconds.' - end subroutine s_generate_particle_beds + end subroutine s_generate_particle_clouds !> Xorshift PRNG. Advances seed in-place and returns a value in [0, 1). function f_xorshift(seed) result(rval) @@ -221,4 +221,4 @@ contains end function f_bin_hash -end module m_particle_bed +end module m_particle_cloud diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 7de68ea9b7..b695694070 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -42,7 +42,7 @@ module m_start_up use m_ibm use m_ib_patches use m_model - use m_particle_bed + use m_particle_cloud use m_collisions use m_compile_specific use m_checker_common @@ -872,21 +872,21 @@ contains if (model_eqns == 3) call s_initialize_internal_energy_equations(q_cons_ts(1)%vf) if (ib) then block - type(ib_patch_parameters), allocatable :: particle_bed_ibs(:) + type(ib_patch_parameters), allocatable :: particle_cloud_ibs(:) if (cfl_dt .and. n_start > 0) then call s_read_ib_restart_data(n_start) - allocate (particle_bed_ibs(0)) + allocate (particle_cloud_ibs(0)) else if (t_step_start > 0) then call s_read_ib_restart_data(t_step_start) - allocate (particle_bed_ibs(0)) + allocate (particle_cloud_ibs(0)) else - call s_generate_particle_beds(particle_bed_ibs) + call s_generate_particle_clouds(particle_cloud_ibs) end if call s_instantiate_STL_models() call s_initialize_ib_airfoils() - call s_reduce_ib_patch_array(particle_bed_ibs) - deallocate (particle_bed_ibs) + call s_reduce_ib_patch_array(particle_cloud_ibs) + deallocate (particle_cloud_ibs) end block call s_ibm_setup() if (t_step_start == 0 .or. (cfl_dt .and. n_start == 0)) then @@ -1158,20 +1158,20 @@ contains end subroutine s_read_ib_restart_data - !> @brief Merges patch_ib (namelist patches, fixed at num_ib_patches_max_namelist) with particle_bed_ibs (CPU-only, exact size) - !! and reduces to only the patches in or near the local computational domain. patch_ib is never reallocated; the local subset is - !! written in-place from the front. particle_bed_ibs is owned by the caller and freed there after this returns. - subroutine s_reduce_ib_patch_array(particle_bed_ibs) + !> @brief Merges patch_ib (namelist patches, fixed at num_ib_patches_max_namelist) with particle_cloud_ibs (CPU-only, exact + !! size) and reduces to only the patches in or near the local computational domain. patch_ib is never reallocated; the local + !! subset is written in-place from the front. particle_cloud_ibs is owned by the caller and freed there after this returns. + subroutine s_reduce_ib_patch_array(particle_cloud_ibs) - type(ib_patch_parameters), intent(in), dimension(:) :: particle_bed_ibs + type(ib_patch_parameters), intent(in), dimension(:) :: particle_cloud_ibs real(wp), dimension(3) :: centroid integer :: i integer :: num_namelist_ibs, num_bed_ibs num_namelist_ibs = num_ibs num_bed_ibs = 0 - do i = 1, num_particle_beds - num_bed_ibs = num_bed_ibs + particle_bed(i)%num_particles + do i = 1, num_particle_clouds + num_bed_ibs = num_bed_ibs + particle_cloud(i)%num_particles end do ! Check for moving IBs across both namelist and particle bed patches. @@ -1184,7 +1184,7 @@ contains end do if (.not. moving_immersed_boundary_flag) then do i = 1, num_bed_ibs - if (particle_bed_ibs(i)%moving_ibm /= 0) then + if (particle_cloud_ibs(i)%moving_ibm /= 0) then moving_immersed_boundary_flag = .true. exit end if @@ -1202,7 +1202,7 @@ contains @:PROHIBIT(num_gbl_ibs > num_ib_patches_max_namelist, & & "Total IB count exceeds patch_ib capacity. Increase num_ib_patches_max_namelist.") do i = 1, num_bed_ibs - patch_ib(num_namelist_ibs + i) = particle_bed_ibs(i) + patch_ib(num_namelist_ibs + i) = particle_cloud_ibs(i) patch_ib(num_namelist_ibs + i)%gbl_patch_id = num_namelist_ibs + i end do num_ibs = num_gbl_ibs @@ -1228,13 +1228,13 @@ contains end if end do do i = 1, num_bed_ibs - centroid = [particle_bed_ibs(i)%x_centroid, particle_bed_ibs(i)%y_centroid, 0._wp] - if (num_dims == 3) centroid(3) = particle_bed_ibs(i)%z_centroid + centroid = [particle_cloud_ibs(i)%x_centroid, particle_cloud_ibs(i)%y_centroid, 0._wp] + if (num_dims == 3) centroid(3) = particle_cloud_ibs(i)%z_centroid if (f_neighborhood_ranks_own_location(centroid)) then num_ibs = num_ibs + 1 @:PROHIBIT(num_ibs > num_ib_patches_max_namelist, & & "Local IB count exceeds patch_ib capacity. Increase num_ib_patches_max_namelist.") - patch_ib(num_ibs) = particle_bed_ibs(i) + patch_ib(num_ibs) = particle_cloud_ibs(i) patch_ib(num_ibs)%gbl_patch_id = num_namelist_ibs + i if (f_local_rank_owns_location(centroid)) then num_local_ibs = num_local_ibs + 1 @@ -1250,7 +1250,7 @@ contains @:PROHIBIT(num_gbl_ibs > num_ib_patches_max_namelist, & & "Total IB count exceeds patch_ib capacity. Increase num_ib_patches_max_namelist.") do i = 1, num_bed_ibs - patch_ib(num_namelist_ibs + i) = particle_bed_ibs(i) + patch_ib(num_namelist_ibs + i) = particle_cloud_ibs(i) patch_ib(num_namelist_ibs + i)%gbl_patch_id = num_namelist_ibs + i end do num_ibs = num_gbl_ibs diff --git a/tests/BA186DDF/golden-metadata.txt b/tests/A7E7A245/golden-metadata.txt similarity index 100% rename from tests/BA186DDF/golden-metadata.txt rename to tests/A7E7A245/golden-metadata.txt diff --git a/tests/BA186DDF/golden.txt b/tests/A7E7A245/golden.txt similarity index 100% rename from tests/BA186DDF/golden.txt rename to tests/A7E7A245/golden.txt diff --git a/toolchain/mfc/case_validator.py b/toolchain/mfc/case_validator.py index 98c10f7834..4848f77892 100644 --- a/toolchain/mfc/case_validator.py +++ b/toolchain/mfc/case_validator.py @@ -577,15 +577,15 @@ def check_ibm(self): ib = self.get("ib", "F") == "T" n = self.get("n", 0) num_ibs = self.get("num_ibs", 0) - num_particle_beds = self.get("num_particle_beds", 0) or 0 + num_particle_clouds = self.get("num_particle_clouds", 0) or 0 ib_state_wrt = self.get("ib_state_wrt", "F") == "T" self.prohibit(ib and n <= 0, "Immersed Boundaries do not work in 1D (requires n > 0)") - has_particle_beds = num_particle_beds > 0 and any((self.get(f"particle_bed({i})%num_particles", 0) or 0) > 0 for i in range(1, num_particle_beds + 1)) + has_particle_clouds = num_particle_clouds > 0 and any((self.get(f"particle_cloud({i})%num_particles", 0) or 0) > 0 for i in range(1, num_particle_clouds + 1)) self.prohibit( - ib and num_ibs <= 0 and not has_particle_beds, - "num_ibs must be >= 1 when ib is enabled (or specify at least one particle_bed with num_particles > 0)", + ib and num_ibs <= 0 and not has_particle_clouds, + "num_ibs must be >= 1 when ib is enabled (or specify at least one particle_cloud with num_particles > 0)", ) num_ib_patches_max = get_fortran_constants().get("num_ib_patches_max_namelist", 50000) self.prohibit( diff --git a/toolchain/mfc/params/definitions.py b/toolchain/mfc/params/definitions.py index 86fd51e688..e22787089b 100644 --- a/toolchain/mfc/params/definitions.py +++ b/toolchain/mfc/params/definitions.py @@ -34,7 +34,7 @@ def _fc(name: str, default: int) -> int: NIB = _fc("num_ib_patches_max", 54000) # patch_ib (Fortran array bound) NAF = _fc("num_ib_airfoils_max", 5) # ib_airfoil (Fortran array bound) NSM = _fc("num_stl_models_max", 10) # stl_models (Fortran array bound) -NPB = _fc("num_particle_beds_max", 10) # particle_bed (Fortran array bound) +NPB = _fc("num_particle_clouds_max", 10) # particle_cloud (Fortran array bound) # Enumeration limits for families not yet converted to IndexedFamily. # These are smaller than the Fortran array bounds to keep the registry compact. # The CONSTRAINTS dict below uses the Fortran constants for validation. @@ -610,7 +610,7 @@ def _load(): # Immersed boundary _r("num_ibs", INT, {"ib"}) _r("num_stl_models", INT, {"ib"}) - _r("num_particle_beds", INT, {"ib"}) + _r("num_particle_clouds", INT, {"ib"}) _r("ib_neighborhood_radius", INT, {"ib"}) _r("ib", LOG, {"ib"}) _r("collision_model", INT, {"ib"}) @@ -913,7 +913,7 @@ def _load(): ) ) - # particle_bed — compact bed specification that expands into individual patch_ib spheres/circles at startup + # particle_cloud — compact bed specification that expands into individual patch_ib spheres/circles at startup _pb_tags = {"ib"} _pb_attrs: Dict[str, tuple] = {} for _d in ["x", "y", "z"]: @@ -925,9 +925,10 @@ def _load(): _pb_attrs["min_spacing"] = (REAL, _pb_tags) _pb_attrs["moving_ibm"] = (INT, _pb_tags) _pb_attrs["seed"] = (INT, _pb_tags) + _pb_attrs["packing_method"] = (INT, _pb_tags) REGISTRY.register_family( IndexedFamily( - base_name="particle_bed", + base_name="particle_cloud", attrs=_pb_attrs, tags=_pb_tags, max_index=NPB, @@ -1251,9 +1252,9 @@ def _nv(targets: set, *names: str) -> None: "coefficient_of_restitution", "collision_time", "ib_coefficient_of_friction", - "num_particle_beds", + "num_particle_clouds", "ib_neighborhood_radius", - "particle_bed", + "particle_cloud", "tau_star", "cont_damage_s", "alpha_bar", diff --git a/toolchain/mfc/params/descriptions.py b/toolchain/mfc/params/descriptions.py index 9627751e70..c165577854 100644 --- a/toolchain/mfc/params/descriptions.py +++ b/toolchain/mfc/params/descriptions.py @@ -134,7 +134,7 @@ "ib": "Enable immersed boundary method", "num_ibs": "Number of immersed boundary patches", "num_stl_models": "Number of STL/OBJ model entries in the stl_models array", - "num_particle_beds": "Number of particle bed specifications to generate immersed boundary patches from", + "num_particle_clouds": "Number of particle bed specifications to generate immersed boundary patches from", "ib_neighborhood_radius": "Neighborhood radius in ranks for IB awareness", # Acoustic sources "acoustic_source": "Enable acoustic source terms", From d8673233c5373901177056ca55ad5289da66a685 Mon Sep 17 00:00:00 2001 From: danieljvickers Date: Thu, 4 Jun 2026 22:16:33 -0400 Subject: [PATCH 2/4] Separated the particle cloud generation to have a master routine that calls downstream routines --- src/simulation/m_particle_cloud.fpp | 287 ++++++++++++++-------------- 1 file changed, 147 insertions(+), 140 deletions(-) diff --git a/src/simulation/m_particle_cloud.fpp b/src/simulation/m_particle_cloud.fpp index b50e188510..7c8f24bd36 100644 --- a/src/simulation/m_particle_cloud.fpp +++ b/src/simulation/m_particle_cloud.fpp @@ -25,21 +25,8 @@ contains impure subroutine s_generate_particle_clouds(particle_cloud_ibs) type(ib_patch_parameters), allocatable, intent(out), dimension(:) :: particle_cloud_ibs - integer :: b, ib_idx, geom - integer :: n_placed, n_total_placed, n_total_particles - integer(8) :: n_attempts, max_attempts - real(wp) :: xmin, xmax, ymin, ymax, zmin, zmax, min_dist - real(wp) :: rx, ry, rz, dist + integer :: cloud_idx, ib_idx, n_total_particles real(wp) :: t_start, t_end - integer :: seed - logical :: overlaps - real(wp), allocatable :: placed(:,:) - - ! Spatial hash grid - integer :: hash_size, slot - integer :: bx, by, bz, nbx, nby, nbz - integer :: dx_b, dy_b, dz_b, dz_lo, dz_hi, j - integer, allocatable :: hash_head(:), chain_next(:) if (num_particle_clouds == 0) then allocate (particle_cloud_ibs(0)) @@ -47,152 +34,172 @@ contains end if call cpu_time(t_start) - n_total_placed = 0 ! Pre-count total particles across all beds so particle_cloud_ibs can be allocated exactly once. n_total_particles = 0 - do b = 1, num_particle_clouds - n_total_particles = n_total_particles + particle_cloud(b)%num_particles + do cloud_idx = 1, num_particle_clouds + n_total_particles = n_total_particles + particle_cloud(cloud_idx)%num_particles end do allocate (particle_cloud_ibs(n_total_particles)) ib_idx = 0 ! index into particle_cloud_ibs - do b = 1, num_particle_clouds - xmin = particle_cloud(b)%x_centroid - 0.5_wp*particle_cloud(b)%length_x - xmax = particle_cloud(b)%x_centroid + 0.5_wp*particle_cloud(b)%length_x - ymin = particle_cloud(b)%y_centroid - 0.5_wp*particle_cloud(b)%length_y - ymax = particle_cloud(b)%y_centroid + 0.5_wp*particle_cloud(b)%length_y - zmin = particle_cloud(b)%z_centroid - 0.5_wp*particle_cloud(b)%length_z - zmax = particle_cloud(b)%z_centroid + 0.5_wp*particle_cloud(b)%length_z + do cloud_idx = 1, num_particle_clouds + select case (particle_cloud(cloud_idx)%packing_method) + case (1) ! random box packing method + call s_particle_cloud_random_box(cloud_idx, ib_idx, particle_cloud_ibs) + end select + end do + + call cpu_time(t_end) + if (proc_rank == 0) print '(a,i0,a,f0.3,a)', 'Particle beds placed ', ib_idx, ' particles in ', t_end - t_start, ' seconds.' + + end subroutine s_generate_particle_clouds + + subroutine s_particle_cloud_random_box(cloud_idx, ib_idx, particle_cloud_ibs) + + integer, intent(in) :: cloud_idx + integer, intent(inout) :: ib_idx + type(ib_patch_parameters), intent(inout), dimension(:) :: particle_cloud_ibs + integer :: n_placed, geom, seed + integer(8) :: n_attempts, max_attempts + real(wp) :: xmin, xmax, ymin, ymax, zmin, zmax, min_dist + real(wp) :: rx, ry, rz, dist + logical :: overlaps + real(wp), allocatable :: placed(:,:) + integer :: hash_size, slot + integer :: bx, by, bz, nbx, nby, nbz + integer :: dx_b, dy_b, dz_b, dz_lo, dz_hi, j + integer, allocatable :: hash_head(:), chain_next(:) + + xmin = particle_cloud(cloud_idx)%x_centroid - 0.5_wp*particle_cloud(cloud_idx)%length_x + xmax = particle_cloud(cloud_idx)%x_centroid + 0.5_wp*particle_cloud(cloud_idx)%length_x + ymin = particle_cloud(cloud_idx)%y_centroid - 0.5_wp*particle_cloud(cloud_idx)%length_y + ymax = particle_cloud(cloud_idx)%y_centroid + 0.5_wp*particle_cloud(cloud_idx)%length_y + zmin = particle_cloud(cloud_idx)%z_centroid - 0.5_wp*particle_cloud(cloud_idx)%length_z + zmax = particle_cloud(cloud_idx)%z_centroid + 0.5_wp*particle_cloud(cloud_idx)%length_z + + min_dist = 2._wp*particle_cloud(cloud_idx)%radius + particle_cloud(cloud_idx)%min_spacing + + if (p == 0) then + geom = 2 ! circle for 2D + dz_lo = 0 + dz_hi = 0 + else + geom = 8 ! sphere for 3D + dz_lo = -1 + dz_hi = 1 + end if + + max_attempts = int(particle_cloud(cloud_idx)%num_particles, 8)*1000_8 + n_placed = 0 + n_attempts = 0 + seed = particle_cloud(cloud_idx)%seed + if (seed == 0) seed = 1 + cloud_idx*1013904223 + + allocate (placed(3, particle_cloud(cloud_idx)%num_particles)) - min_dist = 2._wp*particle_cloud(b)%radius + particle_cloud(b)%min_spacing + ! Hash table: 4x overprovisioned for ~25% load factor, minimum 16 buckets. chain_next(i) links placed particle i to the + ! previous occupant of its bucket. + hash_size = max(16, 4*particle_cloud(cloud_idx)%num_particles) + allocate (hash_head(hash_size)) + allocate (chain_next(particle_cloud(cloud_idx)%num_particles)) + hash_head = -1 + chain_next = -1 + do while (n_placed < particle_cloud(cloud_idx)%num_particles .and. n_attempts < max_attempts) + n_attempts = n_attempts + 1 + + rx = xmin + f_xorshift(seed)*(xmax - xmin) + ry = ymin + f_xorshift(seed)*(ymax - ymin) if (p == 0) then - geom = 2 ! circle for 2D - dz_lo = 0 - dz_hi = 0 + rz = particle_cloud(cloud_idx)%z_centroid else - geom = 8 ! sphere for 3D - dz_lo = -1 - dz_hi = 1 + rz = zmin + f_xorshift(seed)*(zmax - zmin) end if - max_attempts = int(particle_cloud(b)%num_particles, 8)*1000_8 - n_placed = 0 - n_attempts = 0 - seed = particle_cloud(b)%seed - if (seed == 0) seed = 1 + b*1013904223 - - allocate (placed(3, particle_cloud(b)%num_particles)) - - ! Hash table: 4x overprovisioned for ~25% load factor, minimum 16 buckets. chain_next(i) links placed particle i to the - ! previous occupant of its bucket. - hash_size = max(16, 4*particle_cloud(b)%num_particles) - allocate (hash_head(hash_size)) - allocate (chain_next(particle_cloud(b)%num_particles)) - hash_head = -1 - chain_next = -1 - - do while (n_placed < particle_cloud(b)%num_particles .and. n_attempts < max_attempts) - n_attempts = n_attempts + 1 - - rx = xmin + f_xorshift(seed)*(xmax - xmin) - ry = ymin + f_xorshift(seed)*(ymax - ymin) - if (p == 0) then - rz = particle_cloud(b)%z_centroid - else - rz = zmin + f_xorshift(seed)*(zmax - zmin) - end if - - bx = int(floor(rx/min_dist)) - by = int(floor(ry/min_dist)) - bz = 0 - if (p /= 0) bz = int(floor(rz/min_dist)) - - ! Check 3x3(x3) neighboring bins - O(1) average via hash lookup - overlaps = .false. - outer: do dx_b = -1, 1 - do dy_b = -1, 1 - do dz_b = dz_lo, dz_hi - nbx = bx + dx_b - nby = by + dy_b - nbz = bz + dz_b - slot = f_bin_hash(nbx, nby, nbz, hash_size) - j = hash_head(slot) - do while (j > 0) - if (p == 0) then - dist = sqrt((rx - placed(1, j))**2 + (ry - placed(2, j))**2) - else - dist = sqrt((rx - placed(1, j))**2 + (ry - placed(2, j))**2 + (rz - placed(3, j))**2) - end if - if (dist < min_dist) then - overlaps = .true. - exit outer - end if - j = chain_next(j) - end do + bx = int(floor(rx/min_dist)) + by = int(floor(ry/min_dist)) + bz = 0 + if (p /= 0) bz = int(floor(rz/min_dist)) + + ! Check 3x3(x3) neighboring bins - O(1) average via hash lookup + overlaps = .false. + outer: do dx_b = -1, 1 + do dy_b = -1, 1 + do dz_b = dz_lo, dz_hi + nbx = bx + dx_b + nby = by + dy_b + nbz = bz + dz_b + slot = f_bin_hash(nbx, nby, nbz, hash_size) + j = hash_head(slot) + do while (j > 0) + if (p == 0) then + dist = sqrt((rx - placed(1, j))**2 + (ry - placed(2, j))**2) + else + dist = sqrt((rx - placed(1, j))**2 + (ry - placed(2, j))**2 + (rz - placed(3, j))**2) + end if + if (dist < min_dist) then + overlaps = .true. + exit outer + end if + j = chain_next(j) end do end do - end do outer - - if (.not. overlaps) then - n_placed = n_placed + 1 - placed(1, n_placed) = rx - placed(2, n_placed) = ry - placed(3, n_placed) = rz - - ! Insert into hash grid as head of bucket chain - slot = f_bin_hash(bx, by, bz, hash_size) - chain_next(n_placed) = hash_head(slot) - hash_head(slot) = n_placed - - ib_idx = ib_idx + 1 - - ! gbl_patch_id is relative within particle_cloud_ibs here; s_reduce_ib_patch_array adjusts to global indexing. - particle_cloud_ibs(ib_idx)%gbl_patch_id = ib_idx - particle_cloud_ibs(ib_idx)%geometry = geom - particle_cloud_ibs(ib_idx)%x_centroid = rx - particle_cloud_ibs(ib_idx)%y_centroid = ry - particle_cloud_ibs(ib_idx)%z_centroid = rz - particle_cloud_ibs(ib_idx)%step_x_centroid = rx - particle_cloud_ibs(ib_idx)%step_y_centroid = ry - particle_cloud_ibs(ib_idx)%step_z_centroid = rz - particle_cloud_ibs(ib_idx)%angles(:) = 0._wp - particle_cloud_ibs(ib_idx)%step_angles(:) = 0._wp - particle_cloud_ibs(ib_idx)%vel(:) = 0._wp - particle_cloud_ibs(ib_idx)%step_vel(:) = 0._wp - particle_cloud_ibs(ib_idx)%angular_vel(:) = 0._wp - particle_cloud_ibs(ib_idx)%step_angular_vel(:) = 0._wp - particle_cloud_ibs(ib_idx)%force(:) = 0._wp - particle_cloud_ibs(ib_idx)%torque(:) = 0._wp - particle_cloud_ibs(ib_idx)%centroid_offset(:) = 0._wp - particle_cloud_ibs(ib_idx)%rotation_matrix = 0._wp - particle_cloud_ibs(ib_idx)%rotation_matrix(1, 1) = 1._wp - particle_cloud_ibs(ib_idx)%rotation_matrix(2, 2) = 1._wp - particle_cloud_ibs(ib_idx)%rotation_matrix(3, 3) = 1._wp - particle_cloud_ibs(ib_idx)%rotation_matrix_inverse = particle_cloud_ibs(ib_idx)%rotation_matrix - particle_cloud_ibs(ib_idx)%radius = particle_cloud(b)%radius - particle_cloud_ibs(ib_idx)%mass = particle_cloud(b)%mass - particle_cloud_ibs(ib_idx)%moment = dflt_real - particle_cloud_ibs(ib_idx)%moving_ibm = particle_cloud(b)%moving_ibm - particle_cloud_ibs(ib_idx)%slip = .false. - end if - end do - - if (n_placed < particle_cloud(b)%num_particles) then - call s_mpi_abort("Error :: Failed to place all particles in particle bed") + end do + end do outer + + if (.not. overlaps) then + n_placed = n_placed + 1 + placed(1, n_placed) = rx + placed(2, n_placed) = ry + placed(3, n_placed) = rz + + ! Insert into hash grid as head of bucket chain + slot = f_bin_hash(bx, by, bz, hash_size) + chain_next(n_placed) = hash_head(slot) + hash_head(slot) = n_placed + + ib_idx = ib_idx + 1 + + ! gbl_patch_id is relative within particle_cloud_ibs here; s_reduce_ib_patch_array adjusts to global indexing. + particle_cloud_ibs(ib_idx)%gbl_patch_id = ib_idx + particle_cloud_ibs(ib_idx)%geometry = geom + particle_cloud_ibs(ib_idx)%x_centroid = rx + particle_cloud_ibs(ib_idx)%y_centroid = ry + particle_cloud_ibs(ib_idx)%z_centroid = rz + particle_cloud_ibs(ib_idx)%step_x_centroid = rx + particle_cloud_ibs(ib_idx)%step_y_centroid = ry + particle_cloud_ibs(ib_idx)%step_z_centroid = rz + particle_cloud_ibs(ib_idx)%angles(:) = 0._wp + particle_cloud_ibs(ib_idx)%step_angles(:) = 0._wp + particle_cloud_ibs(ib_idx)%vel(:) = 0._wp + particle_cloud_ibs(ib_idx)%step_vel(:) = 0._wp + particle_cloud_ibs(ib_idx)%angular_vel(:) = 0._wp + particle_cloud_ibs(ib_idx)%step_angular_vel(:) = 0._wp + particle_cloud_ibs(ib_idx)%force(:) = 0._wp + particle_cloud_ibs(ib_idx)%torque(:) = 0._wp + particle_cloud_ibs(ib_idx)%centroid_offset(:) = 0._wp + particle_cloud_ibs(ib_idx)%rotation_matrix = 0._wp + particle_cloud_ibs(ib_idx)%rotation_matrix(1, 1) = 1._wp + particle_cloud_ibs(ib_idx)%rotation_matrix(2, 2) = 1._wp + particle_cloud_ibs(ib_idx)%rotation_matrix(3, 3) = 1._wp + particle_cloud_ibs(ib_idx)%rotation_matrix_inverse = particle_cloud_ibs(ib_idx)%rotation_matrix + particle_cloud_ibs(ib_idx)%radius = particle_cloud(cloud_idx)%radius + particle_cloud_ibs(ib_idx)%mass = particle_cloud(cloud_idx)%mass + particle_cloud_ibs(ib_idx)%moment = dflt_real + particle_cloud_ibs(ib_idx)%moving_ibm = particle_cloud(cloud_idx)%moving_ibm + particle_cloud_ibs(ib_idx)%slip = .false. end if - - n_total_placed = n_total_placed + n_placed - deallocate (placed, hash_head, chain_next) end do - call cpu_time(t_end) - if (proc_rank == 0) print '(a,i0,a,f0.3,a)', 'Particle beds placed ', n_total_placed, ' particles in ', t_end - t_start, & - & ' seconds.' + if (n_placed < particle_cloud(cloud_idx)%num_particles) then + call s_mpi_abort("Error :: Failed to place all particles in particle bed") + end if - end subroutine s_generate_particle_clouds + deallocate (placed, hash_head, chain_next) + + end subroutine s_particle_cloud_random_box !> Xorshift PRNG. Advances seed in-place and returns a value in [0, 1). function f_xorshift(seed) result(rval) From 4fbeacd2376b169b6747f2699a9fd2bb3254caee Mon Sep 17 00:00:00 2001 From: danieljvickers Date: Thu, 4 Jun 2026 22:18:35 -0400 Subject: [PATCH 3/4] Adjusted comments --- src/simulation/m_particle_cloud.fpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/simulation/m_particle_cloud.fpp b/src/simulation/m_particle_cloud.fpp index 7c8f24bd36..1b023e6e41 100644 --- a/src/simulation/m_particle_cloud.fpp +++ b/src/simulation/m_particle_cloud.fpp @@ -19,9 +19,7 @@ module m_particle_cloud contains - !> Generate all particle beds and fill particle_cloud_ibs. Called on all ranks before s_reduce_ib_patch_array. Uses a spatial - !! hash grid (cell size = min_dist) so each candidate requires only 3^dim distance checks on average instead of O(n). The - !! placement is fully deterministic given the per-bed seed, so all ranks produce an identical result without MPI. + !> Generate all particle beds and fill particle_cloud_ibs. Called on all ranks before s_reduce_ib_patch_array. impure subroutine s_generate_particle_clouds(particle_cloud_ibs) type(ib_patch_parameters), allocatable, intent(out), dimension(:) :: particle_cloud_ibs @@ -56,6 +54,7 @@ contains end subroutine s_generate_particle_clouds + !> Generates a random distributions of particles in a box with a minimum spacing subroutine s_particle_cloud_random_box(cloud_idx, ib_idx, particle_cloud_ibs) integer, intent(in) :: cloud_idx From 8d0583707f0b7d239f55e67559488f53c40a978e Mon Sep 17 00:00:00 2001 From: danieljvickers Date: Fri, 5 Jun 2026 00:02:19 -0400 Subject: [PATCH 4/4] Added python prohibit check for particle cloud packing method --- toolchain/mfc/case_validator.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/toolchain/mfc/case_validator.py b/toolchain/mfc/case_validator.py index 4848f77892..a9033a9fdb 100644 --- a/toolchain/mfc/case_validator.py +++ b/toolchain/mfc/case_validator.py @@ -595,6 +595,17 @@ def check_ibm(self): self.prohibit(not ib and num_ibs > 0, "num_ibs is set, but ib is not enabled") self.prohibit(ib_state_wrt and not ib, "ib_state_wrt requires ib to be enabled") + for i in range(1, num_particle_clouds + 1): + packing_method = self.get(f"particle_cloud({i})%packing_method", None) + self.prohibit( + packing_method is None, + f"particle_cloud({i})%packing_method must be specified (1 = rejection sampling)", + ) + self.prohibit( + packing_method is not None and packing_method not in [1], + f"particle_cloud({i})%packing_method must be 1 (rejection sampling is the only supported method)", + ) + num_ib_airfoils_max = get_fortran_constants().get("num_ib_airfoils_max", 5) num_stl_models_max = get_fortran_constants().get("num_stl_models_max", 10) num_stl_models = self.get("num_stl_models", 0)