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