MFC
Exascale flow solver
Loading...
Searching...
No Matches
m_particle_cloud.fpp.f90
Go to the documentation of this file.
1# 1 "/home/runner/work/MFC/MFC/src/simulation/m_particle_cloud.fpp"
2!>
3!! @file m_particle_cloud.fpp
4!! @brief Generates particle beds: converts particle_cloud specifications into
5!! individual sphere/circle particle_cloud_ibs entries before reduction.
6
7# 1 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 1
8# 1 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 1
9# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
10# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
11# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
12# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
13# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
14# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
15
16# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
17# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
18# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
19
20# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
21
22# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
23
24# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
25
26# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
27
28# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
29
30# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
31
32# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
33
34# 167 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
35! New line at end of file is required for FYPP
36# 2 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
37# 1 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 1
38# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
39# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
40# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
41# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
42# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
43# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
44
45# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
46# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
47# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
48
49# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
50
51# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
52
53# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
54
55# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
56
57# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
58
59# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
60
61# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
62
63# 167 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
64! New line at end of file is required for FYPP
65# 2 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 2
66
67# 4 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
68# 5 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
69# 6 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
70# 7 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
71# 8 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
72
73# 20 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
74
75# 43 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
76
77# 48 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
78
79# 53 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
80
81# 58 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
82
83# 63 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
84
85# 68 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
86
87# 76 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
88
89# 81 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
90
91# 86 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
92
93# 91 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
94
95# 96 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
96
97# 101 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
98
99# 106 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
100
101# 111 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
102
103# 116 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
104
105# 121 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
106
107# 151 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
108
109# 192 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
110
111# 206 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
112
113# 231 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
114
115# 242 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
116
117# 244 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
118# 255 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
119
120# 284 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
121
122# 294 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
123
124# 304 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
125
126# 313 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
127
128# 330 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
129
130# 340 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
131
132# 347 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
133
134# 353 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
135
136# 359 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
137
138# 365 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
139
140# 371 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
141
142# 377 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
143! New line at end of file is required for FYPP
144# 3 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
145# 1 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 1
146# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
147# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
148# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
149# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
150# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
151# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
152
153# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
154# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
155# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
156
157# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
158
159# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
160
161# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
162
163# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
164
165# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
166
167# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
168
169# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
170
171# 167 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
172! New line at end of file is required for FYPP
173# 2 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 2
174
175# 7 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
176
177# 17 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
178
179# 22 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
180
181# 27 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
182
183# 32 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
184
185# 37 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
186
187# 42 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
188
189# 47 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
190
191# 52 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
192
193# 57 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
194
195# 62 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
196
197# 73 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
198
199# 78 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
200
201# 83 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
202
203# 88 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
204
205# 103 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
206
207# 131 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
208
209# 160 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
210
211# 175 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
212
213# 193 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
214
215# 215 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
216
217# 244 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
218
219# 259 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
220
221# 269 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
222
223# 278 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
224
225# 294 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
226
227# 304 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
228
229# 311 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
230! New line at end of file is required for FYPP
231# 4 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
232
233! GPU parallel region (scalar reductions, maxval/minval)
234# 23 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
235
236! GPU parallel loop over threads (most common GPU macro)
237# 43 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
238
239! Required closing for GPU_PARALLEL_LOOP
240# 55 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
241
242! Mark routine for device compilation
243# 112 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
244
245! Declare device-resident data
246# 130 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
247
248! Inner loop within a GPU parallel region
249# 145 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
250
251! Scoped GPU data region
252# 164 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
253
254! Host code with device pointers (for MPI with GPU buffers)
255# 193 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
256
257! Allocate device memory (unscoped)
258# 207 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
259
260! Free device memory
261# 219 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
262
263! Atomic operation on device
264# 231 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
265
266! End atomic capture block
267# 242 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
268
269! Copy data between host and device
270# 254 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
271
272! Synchronization barrier
273# 266 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
274
275! Import GPU library module (openacc or omp_lib)
276# 275 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
277
278! Emit code only for AMD compiler
279# 282 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
280
281! Emit code for non-Cray compilers
282# 289 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
283
284! Emit code only for Cray compiler
285# 296 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
286
287! Emit code for non-NVIDIA compilers
288# 303 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
289
290# 305 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
291# 306 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
292! New line at end of file is required for FYPP
293# 2 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 2
294
295# 14 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
296
297! Caution: This macro requires the use of a binding script to set CUDA_VISIBLE_DEVICES, such that we have one GPU device per MPI
298! rank. That's because for both cudaMemAdvise (preferred location) and cudaMemPrefetchAsync we use location = device_id = 0. For an
299! example see misc/nvidia_uvm/bind.sh.
300# 55 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
301
302! Allocate and create GPU device memory
303# 75 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
304
305! Free GPU device memory and deallocate
306# 83 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
307
308! Cray-specific GPU pointer setup for vector fields
309# 107 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
310
311! Cray-specific GPU pointer setup for scalar fields
312# 123 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
313
314! Cray-specific GPU pointer setup for acoustic source spatials
315# 148 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
316
317# 154 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
318
319# 161 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
320! New line at end of file is required for FYPP
321# 7 "/home/runner/work/MFC/MFC/src/simulation/m_particle_cloud.fpp" 2
322
323!> @brief Generates particle beds by converting particle_cloud patch specifications into individual immersed boundary patches before
324!! domain reduction. Each rank runs the same deterministic placement so no MPI broadcast of particle positions is needed.
325module m_particle_cloud
326
327 use m_global_parameters
328 use m_constants
329 use m_mpi_common
330 use m_collisions
331
332 implicit none
333
334 private
335
336 public :: s_generate_particle_clouds
337
338contains
339
340 !> Generate all particle beds and fill particle_cloud_ibs. Called on all ranks before s_reduce_ib_patch_array. Each packing
341 !! method owns and allocates its own per-cloud working array (see s_particle_cloud_lattice / s_particle_cloud_rejection_pack)
342 !! and hands back only the entries that fall within this rank's IB neighborhood. Only the first num_particle_cloud_ibs of them
343 !! are actually written - callers must use that count, not size(particle_cloud_ibs), since the remainder of the array is left
344 !! uninitialized.
345 impure subroutine s_generate_particle_clouds(particle_cloud_ibs, num_particle_cloud_ibs)
346
347 type(ib_patch_parameters), allocatable, intent(out), dimension(:) :: particle_cloud_ibs
348 integer, intent(out) :: num_particle_cloud_ibs
349 type(ib_patch_parameters), allocatable :: cloud_ibs(:)
350 integer :: cloud_idx, glbl_idx, num_cloud_ibs, n_total_particles
351 real(wp) :: t_start, t_end
352
353 if (num_particle_clouds == 0) then
354 allocate (particle_cloud_ibs(0))
355 num_particle_cloud_ibs = 0
356 return
357 end if
358
359 call cpu_time(t_start)
360
361 n_total_particles = 0
362 do cloud_idx = 1, num_particle_clouds
363 n_total_particles = n_total_particles + particle_cloud(cloud_idx)%num_particles
364 end do
365 allocate (particle_cloud_ibs(min(num_ib_patches_max_namelist, n_total_particles)))
366
367 num_particle_cloud_ibs = 0
368 glbl_idx = num_ibs
369
370 do cloud_idx = 1, num_particle_clouds
371 ! Dispatch on packing method only: rejection sampling handles both box and hemisphere-shell
372 ! geometries (the per-candidate geometry sampling lives in s_sample_cloud_candidate), while lattice
373 ! packing is box-only - the hemisphere-shell + lattice combination is rejected in case_validator.py.
374 select case (particle_cloud(cloud_idx)%packing_method)
375 case (1) ! rejection (random) packing method
376 call s_particle_cloud_rejection_pack(cloud_idx, glbl_idx, cloud_ibs, num_cloud_ibs)
377 case (2) ! lattice packing method
378 call s_particle_cloud_lattice(cloud_idx, glbl_idx, cloud_ibs, num_cloud_ibs)
379 case default
380 call s_mpi_abort("Particle cloud packing method is not a known packing method of MFC. Exiting.")
381 end select
382
383 if (num_particle_cloud_ibs + num_cloud_ibs > num_ib_patches_max_namelist) then
384# 68 "/home/runner/work/MFC/MFC/src/simulation/m_particle_cloud.fpp"
385 call s_prohibit_abort("num_particle_cloud_ibs + num_cloud_ibs > num_ib_patches_max_namelist", "Too many particle-cloud IBs in one rank's neighborhood. Modify case file or increase num_ib_patches_max_namelist.")
386# 68 "/home/runner/work/MFC/MFC/src/simulation/m_particle_cloud.fpp"
387 end if
388# 70 "/home/runner/work/MFC/MFC/src/simulation/m_particle_cloud.fpp"
389 particle_cloud_ibs(num_particle_cloud_ibs + 1:num_particle_cloud_ibs + num_cloud_ibs) = cloud_ibs(1:num_cloud_ibs)
390 num_particle_cloud_ibs = num_particle_cloud_ibs + num_cloud_ibs
391 deallocate (cloud_ibs)
392 end do
393
394 call cpu_time(t_end)
395 if (proc_rank == 0) print '(a,i0,a,f0.3,a)', 'Particle beds placed ', glbl_idx - num_ibs, ' particles in ', &
396 & t_end - t_start, ' seconds.'
397
398 end subroutine s_generate_particle_clouds
399
400 !> Rejection-samples particle centres into a box or hemisphere-shell region with a minimum centre-to-centre spacing. Rejection
401 !! sampling needs every placed particle tracked (regardless of which rank's neighborhood it falls in) to detect overlaps
402 !! deterministically, so cloud_ibs is allocated here to the cloud's full requested particle count and only pared down to this
403 !! rank's neighborhood afterwards, via s_reduce_particle_cloud_ibs. Only the per-candidate geometry sampling differs between box
404 !! and hemisphere shell; it is delegated to s_sample_cloud_candidate, and every other step (overlap rejection via the spatial
405 !! hash, acceptance, reduction) is geometry-independent.
406 subroutine s_particle_cloud_rejection_pack(cloud_idx, glbl_idx, cloud_ibs, num_cloud_ibs)
407
408 integer, intent(in) :: cloud_idx
409 integer, intent(inout) :: glbl_idx
410 type(ib_patch_parameters), allocatable, intent(out), dimension(:) :: cloud_ibs
411 integer, intent(out) :: num_cloud_ibs
412 integer :: ib_idx, n_placed, geom, seed, alloc_stat
413 integer(8) :: n_attempts, max_attempts
414 real(wp) :: min_dist, rx, ry, rz
415 logical :: overlaps, reject, periodic_pack
416 real(wp), allocatable :: placed(:,:)
417 integer :: hash_size, slot
418 integer :: bx, by, bz, nx_bins, ny_bins, nz_bins
419 integer, allocatable :: hash_head(:), chain_next(:)
420 real(wp) :: xmin, ymin, zmin, length_x, length_y, length_z
421
422 allocate (cloud_ibs(particle_cloud(cloud_idx)%num_particles), stat=alloc_stat)
423 if (alloc_stat /= 0) then
424 call s_mpi_abort("Error :: Ran out of CPU memory trying to allocate particle cloud IB array. " &
425 & // "Current system resources cannot perform rejection packing with the specified number of particles.")
426 end if
427 ib_idx = 0
428
429 min_dist = 2._wp*particle_cloud(cloud_idx)%radius + particle_cloud(cloud_idx)%min_spacing
430 periodic_pack = particle_cloud(cloud_idx)%cloud_geometry == 1 .and. particle_cloud(cloud_idx)%periodic == 1
431 length_x = particle_cloud(cloud_idx)%length_x
432 length_y = particle_cloud(cloud_idx)%length_y
433 length_z = particle_cloud(cloud_idx)%length_z
434 xmin = particle_cloud(cloud_idx)%x_centroid - 0.5_wp*length_x
435 ymin = particle_cloud(cloud_idx)%y_centroid - 0.5_wp*length_y
436 zmin = particle_cloud(cloud_idx)%z_centroid - 0.5_wp*length_z
437 nx_bins = max(1, ceiling(length_x/min_dist))
438 ny_bins = max(1, ceiling(length_y/min_dist))
439 nz_bins = max(1, ceiling(length_z/min_dist))
440 if (num_dims < 3) nz_bins = 1
441
442 if (num_dims < 3) then
443 geom = 2 ! circle for 2D
444 else
445 geom = 8 ! sphere for 3D
446 end if
447
448 max_attempts = int(particle_cloud(cloud_idx)%num_particles, 8)*1000_8
449 n_placed = 0
450 n_attempts = 0
451 seed = particle_cloud(cloud_idx)%seed
452 if (seed == 0) seed = 1 + cloud_idx*1013904223
453
454 allocate (placed(3, particle_cloud(cloud_idx)%num_particles))
455
456 ! Hash table: 4x overprovisioned for ~25% load factor, minimum 16 buckets. chain_next(i) links placed particle i to the
457 ! previous occupant of its bucket.
458 hash_size = max(16, 4*particle_cloud(cloud_idx)%num_particles)
459 allocate (hash_head(hash_size))
460 allocate (chain_next(particle_cloud(cloud_idx)%num_particles))
461 hash_head = -1
462 chain_next = -1
463
464 do while (n_placed < particle_cloud(cloud_idx)%num_particles .and. n_attempts < max_attempts)
465 n_attempts = n_attempts + 1
466
467 call s_sample_cloud_candidate(cloud_idx, seed, rx, ry, rz, reject)
468 if (reject) cycle
469
470 call s_check_cloud_particle_overlap(rx, ry, rz, placed, hash_head, chain_next, hash_size, min_dist, periodic_pack, &
471 & xmin, ymin, zmin, length_x, length_y, length_z, nx_bins, ny_bins, nz_bins, &
472 & overlaps, bx, by, bz)
473
474 if (.not. overlaps) then
475 n_placed = n_placed + 1
476 placed(1, n_placed) = rx
477 placed(2, n_placed) = ry
478 placed(3, n_placed) = rz
479
480 ! Insert into hash grid as head of bucket chain
481 slot = f_bin_hash(bx, by, bz, hash_size)
482 chain_next(n_placed) = hash_head(slot)
483 hash_head(slot) = n_placed
484
485 glbl_idx = glbl_idx + 1
486 call s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, rx, ry, rz, cloud_ibs)
487 end if
488 end do
489
490 if (n_placed < particle_cloud(cloud_idx)%num_particles) then
491 call s_mpi_abort("Error :: Failed to place all particles in particle bed")
492 end if
493
494 deallocate (placed, hash_head, chain_next)
495
496 call s_reduce_particle_cloud_ibs(cloud_ibs, ib_idx)
497 num_cloud_ibs = ib_idx
498
499 end subroutine s_particle_cloud_rejection_pack
500
501 !> Draws one rejection-sampling candidate centre (rx, ry, rz) for cloud_idx, advancing seed in place. For box geometry the
502 !! candidate is uniform in the box and never rejected. For a hemisphere shell the candidate is uniform in the shell volume - 2D
503 !! uses theta uniform on [0, pi] with the sqrt radial CDF; 3D uses uniform phi, uniform cos(polar) on [0, 1], and the cube-root
504 !! radial CDF - and reject is set when it lands within one particle radius of the flat face (the plane at y_centroid in 2D,
505 !! z_centroid in 3D), a hard geometric cut applied after sampling that preserves uniformity over the remaining region.
506 subroutine s_sample_cloud_candidate(cloud_idx, seed, rx, ry, rz, reject)
507
508 integer, intent(in) :: cloud_idx
509 integer, intent(inout) :: seed
510 real(wp), intent(out) :: rx, ry, rz
511 logical, intent(out) :: reject
512 real(wp) :: xmin, xmax, ymin, ymax, zmin, zmax
513 real(wp) :: theta, phi, r_shell, rho, u, zdir, r_inner, r_outer
514
515 reject = .false.
516
517 select case (particle_cloud(cloud_idx)%cloud_geometry)
518 case (1) ! box
519 xmin = particle_cloud(cloud_idx)%x_centroid - 0.5_wp*particle_cloud(cloud_idx)%length_x
520 xmax = particle_cloud(cloud_idx)%x_centroid + 0.5_wp*particle_cloud(cloud_idx)%length_x
521 ymin = particle_cloud(cloud_idx)%y_centroid - 0.5_wp*particle_cloud(cloud_idx)%length_y
522 ymax = particle_cloud(cloud_idx)%y_centroid + 0.5_wp*particle_cloud(cloud_idx)%length_y
523 zmin = particle_cloud(cloud_idx)%z_centroid - 0.5_wp*particle_cloud(cloud_idx)%length_z
524 zmax = particle_cloud(cloud_idx)%z_centroid + 0.5_wp*particle_cloud(cloud_idx)%length_z
525
526 rx = xmin + f_xorshift(seed)*(xmax - xmin)
527 ry = ymin + f_xorshift(seed)*(ymax - ymin)
528 if (num_dims < 3) then
529 rz = particle_cloud(cloud_idx)%z_centroid
530 else
531 rz = zmin + f_xorshift(seed)*(zmax - zmin)
532 end if
533 case (2) ! hemisphere shell
534 r_inner = particle_cloud(cloud_idx)%shell_inner_radius + particle_cloud(cloud_idx)%radius
535 r_outer = particle_cloud(cloud_idx)%shell_outer_radius - particle_cloud(cloud_idx)%radius
536
537 if (num_dims < 3) then
538 theta = pi*f_xorshift(seed)
539 u = f_xorshift(seed)
540 r_shell = sqrt((r_outer**2 - r_inner**2)*u + r_inner**2)
541 rx = particle_cloud(cloud_idx)%x_centroid + r_shell*cos(theta)
542 ry = particle_cloud(cloud_idx)%y_centroid + r_shell*sin(theta)
543 rz = particle_cloud(cloud_idx)%z_centroid
544 if (ry < particle_cloud(cloud_idx)%y_centroid + particle_cloud(cloud_idx)%radius) reject = .true.
545 else
546 phi = 2._wp*pi*f_xorshift(seed)
547 zdir = f_xorshift(seed)
548 rho = sqrt(max(0._wp, 1._wp - zdir**2))
549 u = f_xorshift(seed)
550 r_shell = ((r_outer**3 - r_inner**3)*u + r_inner**3)**(1._wp/3._wp)
551 rx = particle_cloud(cloud_idx)%x_centroid + r_shell*rho*cos(phi)
552 ry = particle_cloud(cloud_idx)%y_centroid + r_shell*rho*sin(phi)
553 rz = particle_cloud(cloud_idx)%z_centroid + r_shell*zdir
554 if (rz < particle_cloud(cloud_idx)%z_centroid + particle_cloud(cloud_idx)%radius) reject = .true.
555 end if
556 case default
557 call s_mpi_abort("Particle cloud geometry is not a known cloud geometry of MFC. Exiting.")
558 end select
559
560 end subroutine s_sample_cloud_candidate
561
562 !> Places particles on the optimally dense lattice for the cloud region: a triangular lattice in 2D, a face-centered cubic
563 !! lattice in 3D. The lattice spacing is set by the particle density (num_particles over the region area/volume); if that
564 !! spacing falls below the required centre-to-centre distance (2*radius + min_spacing), the region is too dense and the run is
565 !! aborted. No two lattice sites can overlap, so unlike rejection packing each site's IB neighborhood membership
566 !! (get_neighbor_bounds() must already have run) is checked as it is generated and only in-neighborhood sites are stored;
567 !! cloud_ibs is therefore allocated to the neighborhood-sized cap rather than the cloud's full particle count.
568 subroutine s_particle_cloud_lattice(cloud_idx, glbl_idx, cloud_ibs, num_cloud_ibs)
569
570 integer, intent(in) :: cloud_idx
571 integer, intent(inout) :: glbl_idx
572 type(ib_patch_parameters), allocatable, intent(out), dimension(:) :: cloud_ibs
573 integer, intent(out) :: num_cloud_ibs
574 integer :: ib_idx, n_placed, n_target, geom
575 integer :: row, col, ncx, ncy, ix, jy, kz, b
576 real(wp) :: xmin, xmax, ymin, ymax, zmin, zmax, min_dist
577 real(wp) :: spacing, row_dy, cell, x0, px, py
578 real(wp), dimension(4) :: bx_off, by_off, bz_off
579 real(wp), dimension(3) :: centroid
580
581 allocate (cloud_ibs(min(num_ib_patches_max_namelist, particle_cloud(cloud_idx)%num_particles)))
582 ib_idx = 0
583
584 xmin = particle_cloud(cloud_idx)%x_centroid - 0.5_wp*particle_cloud(cloud_idx)%length_x
585 xmax = particle_cloud(cloud_idx)%x_centroid + 0.5_wp*particle_cloud(cloud_idx)%length_x
586 ymin = particle_cloud(cloud_idx)%y_centroid - 0.5_wp*particle_cloud(cloud_idx)%length_y
587 ymax = particle_cloud(cloud_idx)%y_centroid + 0.5_wp*particle_cloud(cloud_idx)%length_y
588 zmin = particle_cloud(cloud_idx)%z_centroid - 0.5_wp*particle_cloud(cloud_idx)%length_z
589 zmax = particle_cloud(cloud_idx)%z_centroid + 0.5_wp*particle_cloud(cloud_idx)%length_z
590
591 min_dist = 2._wp*particle_cloud(cloud_idx)%radius + particle_cloud(cloud_idx)%min_spacing
592 n_target = particle_cloud(cloud_idx)%num_particles
593 n_placed = 0
594
595 if (num_dims < 3) then
596 geom = 2 ! circle for 2D
597 ! Triangular lattice: area per particle = (sqrt(3)/2)*spacing**2.
598 spacing = sqrt(2._wp*(xmax - xmin)*(ymax - ymin)/(sqrt(3._wp)*real(n_target, wp)))
599 else
600 geom = 8 ! sphere for 3D
601 ! Face-centered cubic lattice: volume per particle = spacing**3/sqrt(2).
602 spacing = (sqrt(2._wp)*(xmax - xmin)*(ymax - ymin)*(zmax - zmin)/real(n_target, wp))**(1._wp/3._wp)
603 end if
604
605 if (spacing < min_dist) then
606 call s_mpi_abort("Error :: Particle cloud is too dense for lattice packing; " &
607 & // "reduce num_particles or min_spacing, or enlarge the cloud region")
608 end if
609
610 if (num_dims < 3) then
611 ! Triangular lattice: rows pitched by spacing*sqrt(3)/2, odd rows shifted by half a spacing.
612 row_dy = spacing*sqrt(3._wp)/2._wp
613 row = 0
614 do while (n_placed < n_target)
615 py = ymin + real(row, wp)*row_dy
616 x0 = xmin
617 if (mod(row, 2) == 1) x0 = xmin + 0.5_wp*spacing
618 col = 0
619 px = x0
620 do while (px <= xmax .and. n_placed < n_target)
621 glbl_idx = glbl_idx + 1
622 centroid = [px, py, particle_cloud(cloud_idx)%z_centroid]
623 if (f_neighborhood_ranks_own_location(centroid)) then
624 call s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, centroid(1), centroid(2), centroid(3), &
625 & cloud_ibs)
626 end if
627 n_placed = n_placed + 1
628 col = col + 1
629 px = x0 + real(col, wp)*spacing
630 end do
631 row = row + 1
632 end do
633 else
634 ! Face-centered cubic lattice via the conventional cubic cell (side = spacing*sqrt(2)) and its four basis points.
635 cell = spacing*sqrt(2._wp)
636 bx_off = [0._wp, 0.5_wp, 0.5_wp, 0._wp]*cell
637 by_off = [0._wp, 0.5_wp, 0._wp, 0.5_wp]*cell
638 bz_off = [0._wp, 0._wp, 0.5_wp, 0.5_wp]*cell
639 ncx = max(1, ceiling((xmax - xmin)/cell))
640 ncy = max(1, ceiling((ymax - ymin)/cell))
641 kz = 0
642 do while (n_placed < n_target)
643 do jy = 0, ncy - 1
644 do ix = 0, ncx - 1
645 do b = 1, 4
646 if (n_placed >= n_target) exit
647 centroid = [xmin + real(ix, wp)*cell + bx_off(b), ymin + real(jy, wp)*cell + by_off(b), &
648 & zmin + real(kz, wp)*cell + bz_off(b)]
649 glbl_idx = glbl_idx + 1
650 if (f_neighborhood_ranks_own_location(centroid)) then
651 call s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, centroid(1), centroid(2), &
652 & centroid(3), cloud_ibs)
653 end if
654 n_placed = n_placed + 1
655 end do
656 end do
657 end do
658 kz = kz + 1
659 end do
660 end if
661
662 num_cloud_ibs = ib_idx
663
664 end subroutine s_particle_cloud_lattice
665
666 !> Writes a single placed particle into particle_cloud_ibs at the next free slot, advancing ib_idx. The caller decides whether
667 !! this particle belongs in the array (neighborhood membership, for lattice packing, or unconditionally for rejection packing -
668 !! see s_particle_cloud_lattice / s_particle_cloud_rejection_pack) and supplies its already-assigned, absolute global patch id
669 !! via glbl_idx - s_reduce_ib_patch_array copies gbl_patch_id as-is. Shared by all packing methods so the per-particle
670 !! ib_patch_parameters setup stays in one place.
671 subroutine s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, px, py, pz, particle_cloud_ibs)
672
673 integer, intent(in) :: cloud_idx, glbl_idx, geom
674 integer, intent(inout) :: ib_idx
675 real(wp), intent(in) :: px, py, pz
676 type(ib_patch_parameters), intent(inout), dimension(:) :: particle_cloud_ibs
677
678 ib_idx = ib_idx + 1
679 if (ib_idx > size(particle_cloud_ibs)) then
680# 360 "/home/runner/work/MFC/MFC/src/simulation/m_particle_cloud.fpp"
681 call s_prohibit_abort("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.")
682# 360 "/home/runner/work/MFC/MFC/src/simulation/m_particle_cloud.fpp"
683 end if
684# 362 "/home/runner/work/MFC/MFC/src/simulation/m_particle_cloud.fpp"
685
686 particle_cloud_ibs(ib_idx)%gbl_patch_id = glbl_idx
687 particle_cloud_ibs(ib_idx)%geometry = geom
688 particle_cloud_ibs(ib_idx)%x_centroid = px
689 particle_cloud_ibs(ib_idx)%y_centroid = py
690 particle_cloud_ibs(ib_idx)%z_centroid = pz
691 particle_cloud_ibs(ib_idx)%step_x_centroid = px
692 particle_cloud_ibs(ib_idx)%step_y_centroid = py
693 particle_cloud_ibs(ib_idx)%step_z_centroid = pz
694 particle_cloud_ibs(ib_idx)%angles(:) = 0._wp
695 particle_cloud_ibs(ib_idx)%step_angles(:) = 0._wp
696 particle_cloud_ibs(ib_idx)%vel(:) = 0._wp
697 particle_cloud_ibs(ib_idx)%step_vel(:) = 0._wp
698 particle_cloud_ibs(ib_idx)%angular_vel(:) = 0._wp
699 particle_cloud_ibs(ib_idx)%step_angular_vel(:) = 0._wp
700 particle_cloud_ibs(ib_idx)%force(:) = 0._wp
701 particle_cloud_ibs(ib_idx)%torque(:) = 0._wp
702 particle_cloud_ibs(ib_idx)%centroid_offset(:) = 0._wp
703 particle_cloud_ibs(ib_idx)%rotation_matrix = 0._wp
704 particle_cloud_ibs(ib_idx)%rotation_matrix(1, 1) = 1._wp
705 particle_cloud_ibs(ib_idx)%rotation_matrix(2, 2) = 1._wp
706 particle_cloud_ibs(ib_idx)%rotation_matrix(3, 3) = 1._wp
707 particle_cloud_ibs(ib_idx)%rotation_matrix_inverse = particle_cloud_ibs(ib_idx)%rotation_matrix
708 particle_cloud_ibs(ib_idx)%radius = particle_cloud(cloud_idx)%radius
709 particle_cloud_ibs(ib_idx)%mass = particle_cloud(cloud_idx)%mass
710 particle_cloud_ibs(ib_idx)%moment = dflt_real
711 particle_cloud_ibs(ib_idx)%moving_ibm = particle_cloud(cloud_idx)%moving_ibm
712 particle_cloud_ibs(ib_idx)%slip = .false.
713
714 ! Particles are inert surfaces. These must be set explicitly: particle_cloud_ibs is
715 ! allocated (not default-initialized) and s_reduce_ib_patch_array copies the whole
716 ! struct into patch_ib, overwriting the defaults from
717 ! s_assign_default_values_to_user_inputs -- so anything left unset here reaches the
718 ! solver as uninitialized memory (a nonzero v_blow injects a garbage wall-normal
719 ! velocity and NaNs the field).
720 particle_cloud_ibs(ib_idx)%v_blow = 0._wp
721 particle_cloud_ibs(ib_idx)%inj_species = 0
722 particle_cloud_ibs(ib_idx)%burn_rate_exp = 0._wp
723 particle_cloud_ibs(ib_idx)%burn_rate_pref = 0._wp
724
725 end subroutine s_add_cloud_particle
726
727 !> Convert a candidate particle centre to spatial-hash bin coordinates.
728 subroutine s_get_cloud_bin(px, py, pz, min_dist, periodic_pack, xmin, ymin, zmin, nx_bins, ny_bins, nz_bins, bx, by, bz)
729
730 real(wp), intent(in) :: px, py, pz, min_dist
731 logical, intent(in) :: periodic_pack
732 real(wp), intent(in) :: xmin, ymin, zmin
733 integer, intent(in) :: nx_bins, ny_bins, nz_bins
734 integer, intent(out) :: bx, by, bz
735
736 if (periodic_pack) then
737 bx = modulo(int(floor((px - xmin)/min_dist)), nx_bins)
738 by = modulo(int(floor((py - ymin)/min_dist)), ny_bins)
739 if (num_dims < 3) then
740 bz = 0
741 else
742 bz = modulo(int(floor((pz - zmin)/min_dist)), nz_bins)
743 end if
744 else
745 bx = int(floor(px/min_dist))
746 by = int(floor(py/min_dist))
747 if (num_dims < 3) then
748 bz = 0
749 else
750 bz = int(floor(pz/min_dist))
751 end if
752 end if
753
754 end subroutine s_get_cloud_bin
755
756 !> Check whether a candidate particle centre overlaps any already-placed particle in neighbouring spatial-hash bins. Scans the
757 !! 3x3(x3) bin neighborhood - O(1) average via hash lookup - and also hands back the candidate's own bin so the caller can
758 !! insert it without recomputing.
759 subroutine s_check_cloud_particle_overlap(px, py, pz, placed, hash_head, chain_next, hash_size, min_dist, periodic_pack, &
760 & xmin, ymin, zmin, length_x, length_y, length_z, nx_bins, ny_bins, nz_bins, overlaps, bx, by, bz)
761
762 real(wp), intent(in) :: px, py, pz, min_dist
763 real(wp), intent(in), dimension(:,:) :: placed
764 integer, intent(in), dimension(:) :: hash_head, chain_next
765 integer, intent(in) :: hash_size
766 logical, intent(in) :: periodic_pack
767 real(wp), intent(in) :: xmin, ymin, zmin, length_x, length_y, length_z
768 integer, intent(in) :: nx_bins, ny_bins, nz_bins
769 logical, intent(out) :: overlaps
770 integer, intent(out) :: bx, by, bz
771 integer :: nbx, nby, nbz, slot
772 integer :: dx_b, dy_b, dz_b, dz_lo, dz_hi, j
773 real(wp) :: dist_sq, min_dist_sq, dx, dy, dz
774
775 call s_get_cloud_bin(px, py, pz, min_dist, periodic_pack, xmin, ymin, zmin, nx_bins, ny_bins, nz_bins, bx, by, bz)
776
777 dz_lo = -1
778 dz_hi = 1
779 if (num_dims < 3) then
780 dz_lo = 0
781 dz_hi = 0
782 end if
783
784 min_dist_sq = min_dist**2
785 overlaps = .false.
786
787 do dx_b = -1, 1
788 do dy_b = -1, 1
789 do dz_b = dz_lo, dz_hi
790 nbx = bx + dx_b
791 nby = by + dy_b
792 nbz = bz + dz_b
793 if (periodic_pack) then
794 nbx = modulo(nbx, nx_bins)
795 nby = modulo(nby, ny_bins)
796 if (num_dims == 3) nbz = modulo(nbz, nz_bins)
797 end if
798 slot = f_bin_hash(nbx, nby, nbz, hash_size)
799 j = hash_head(slot)
800 do while (j > 0)
801 dx = abs(px - placed(1, j))
802 dy = abs(py - placed(2, j))
803 if (periodic_pack) then
804 dx = min(dx, length_x - dx)
805 dy = min(dy, length_y - dy)
806 end if
807 if (num_dims < 3) then
808 dist_sq = dx**2 + dy**2
809 else
810 dz = abs(pz - placed(3, j))
811 if (periodic_pack) dz = min(dz, length_z - dz)
812 dist_sq = dx**2 + dy**2 + dz**2
813 end if
814 if (dist_sq < min_dist_sq) then
815 overlaps = .true.
816 return
817 end if
818 j = chain_next(j)
819 end do
820 end do
821 end do
822 end do
823
824 end subroutine s_check_cloud_particle_overlap
825
826 !> Compacts cloud_ibs(1:num_ibs) in place, discarding entries outside this rank's IB neighborhood (get_neighbor_bounds() must
827 !! already have run) and updating num_ibs to the retained count. Used by rejection packing, which cannot filter as it places
828 !! particles (see s_particle_cloud_rejection_pack), to pare its full, unfiltered placement down to this rank's neighborhood.
829 subroutine s_reduce_particle_cloud_ibs(cloud_ibs, num_cloud_ibs)
830
831 type(ib_patch_parameters), intent(inout), dimension(:) :: cloud_ibs
832 integer, intent(inout) :: num_cloud_ibs
833 integer :: i, write_idx
834 real(wp), dimension(3) :: centroid
835
836 write_idx = 0
837 do i = 1, num_cloud_ibs
838 centroid = [cloud_ibs(i)%x_centroid, cloud_ibs(i)%y_centroid, 0._wp]
839 if (num_dims == 3) centroid(3) = cloud_ibs(i)%z_centroid
840 if (f_neighborhood_ranks_own_location(centroid)) then
841 write_idx = write_idx + 1
842 if (write_idx /= i) cloud_ibs(write_idx) = cloud_ibs(i)
843 end if
844 end do
845 num_cloud_ibs = write_idx
846
847 end subroutine s_reduce_particle_cloud_ibs
848
849 !> Xorshift PRNG. Advances seed in-place and returns a value in [0, 1).
850 function f_xorshift(seed) result(rval)
851
852 integer, intent(inout) :: seed
853 real(wp) :: rval
854
855 seed = ieor(seed, ishft(seed, 13))
856 seed = ieor(seed, ishft(seed, -17))
857 seed = ieor(seed, ishft(seed, 5))
858
859 ! Mask off the sign bit rather than abs(): at seed = -huge-1, abs() overflows and returns the value
860 ! unchanged, yielding rval = -1 and, in the shell sampler, a NaN centre that then passes every overlap
861 ! comparison (all comparisons against NaN are false) and gets placed.
862 rval = real(iand(seed, huge(seed)), wp)/real(huge(seed), wp)
863
864 end function f_xorshift
865
866 !> Hash bin coordinates to a 1-indexed slot in [1, hash_size]. Uses large prime multipliers to spread bins across buckets. Hash
867 !! collisions are benign: the distance check catches false neighbours.
868 function f_bin_hash(bx, by, bz, hash_size) result(slot)
869
870 integer, intent(in) :: bx, by, bz, hash_size
871 integer :: slot
872 integer(8) :: key
873
874 key = ieor(ieor(int(bx, 8)*73856093_8, int(by, 8)*19349663_8), int(bz, 8)*83492791_8)
875 slot = int(mod(abs(key), int(hash_size, 8))) + 1
876
877 end function f_bin_hash
878
879end module m_particle_cloud
m_collisions
Ghost-node immersed boundary method: locates ghost/image points, computes interpolation coefficients,...
Definition m_collisions.fpp.f90:324
m_collisions::f_neighborhood_ranks_own_location
logical function, public f_neighborhood_ranks_own_location(location)
function checks if this local MPI processor owns this specific collision
Definition m_collisions.fpp.f90:1171
m_constants
Compile-time constant parameters: default values, tolerances, and physical constants.
Definition m_constants.fpp.f90:7
m_constants::num_ib_patches_max_namelist
integer, parameter num_ib_patches_max_namelist
Definition m_constants.fpp.f90:33
m_constants::dflt_real
real(wp), parameter dflt_real
Default real value.
Definition m_constants.fpp.f90:12
m_constants::pi
real(wp), parameter pi
Pi.
Definition m_constants.fpp.f90:16
m_global_parameters
Global parameters for the computational domain, fluid properties, and simulation algorithm configurat...
Definition m_global_parameters.fpp.f90:334
m_global_parameters::proc_rank
integer proc_rank
Rank of the local processor.
Definition m_global_parameters.fpp.f90:354
m_mpi_common
MPI communication layer: domain decomposition, halo exchange, reductions, and parallel I/O setup.
Definition m_mpi_common.fpp.f90:334
m_mpi_common::s_mpi_abort
impure subroutine s_mpi_abort(prnt, code)
The subroutine terminates the MPI execution environment.
Definition m_mpi_common.fpp.f90:945
m_mpi_common::s_prohibit_abort
impure subroutine s_prohibit_abort(condition, message)
Print a case file error with the prohibited condition and message, then abort execution.
Definition m_mpi_common.fpp.f90:737
m_particle_cloud
Generates particle beds by converting particle_cloud patch specifications into individual immersed bo...
Definition m_particle_cloud.fpp.f90:325
m_particle_cloud::s_add_cloud_particle
subroutine s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, px, py, pz, particle_cloud_ibs)
Writes a single placed particle into particle_cloud_ibs at the next free slot, advancing ib_idx....
Definition m_particle_cloud.fpp.f90:672
m_particle_cloud::s_reduce_particle_cloud_ibs
subroutine s_reduce_particle_cloud_ibs(cloud_ibs, num_cloud_ibs)
Compacts cloud_ibs(1:num_ibs) in place, discarding entries outside this rank's IB neighborhood (get_n...
Definition m_particle_cloud.fpp.f90:830
m_particle_cloud::f_xorshift
real(wp) function f_xorshift(seed)
Xorshift PRNG. Advances seed in-place and returns a value in [0, 1).
Definition m_particle_cloud.fpp.f90:851
m_particle_cloud::s_check_cloud_particle_overlap
subroutine s_check_cloud_particle_overlap(px, py, pz, placed, hash_head, chain_next, hash_size, min_dist, periodic_pack, xmin, ymin, zmin, length_x, length_y, length_z, nx_bins, ny_bins, nz_bins, overlaps, bx, by, bz)
Check whether a candidate particle centre overlaps any already-placed particle in neighbouring spatia...
Definition m_particle_cloud.fpp.f90:761
m_particle_cloud::f_bin_hash
integer function f_bin_hash(bx, by, bz, hash_size)
Hash bin coordinates to a 1-indexed slot in [1, hash_size]. Uses large prime multipliers to spread bi...
Definition m_particle_cloud.fpp.f90:869
m_particle_cloud::s_generate_particle_clouds
impure subroutine, public s_generate_particle_clouds(particle_cloud_ibs, num_particle_cloud_ibs)
Generate all particle beds and fill particle_cloud_ibs. Called on all ranks before s_reduce_ib_patch_...
Definition m_particle_cloud.fpp.f90:346
m_particle_cloud::s_particle_cloud_lattice
subroutine s_particle_cloud_lattice(cloud_idx, glbl_idx, cloud_ibs, num_cloud_ibs)
Places particles on the optimally dense lattice for the cloud region: a triangular lattice in 2D,...
Definition m_particle_cloud.fpp.f90:569
m_particle_cloud::s_get_cloud_bin
subroutine s_get_cloud_bin(px, py, pz, min_dist, periodic_pack, xmin, ymin, zmin, nx_bins, ny_bins, nz_bins, bx, by, bz)
Convert a candidate particle centre to spatial-hash bin coordinates.
Definition m_particle_cloud.fpp.f90:729
m_particle_cloud::s_particle_cloud_rejection_pack
subroutine s_particle_cloud_rejection_pack(cloud_idx, glbl_idx, cloud_ibs, num_cloud_ibs)
Rejection-samples particle centres into a box or hemisphere-shell region with a minimum centre-to-cen...
Definition m_particle_cloud.fpp.f90:407
m_particle_cloud::s_sample_cloud_candidate
subroutine s_sample_cloud_candidate(cloud_idx, seed, rx, ry, rz, reject)
Draws one rejection-sampling candidate centre (rx, ry, rz) for cloud_idx, advancing seed in place....
Definition m_particle_cloud.fpp.f90:507