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