MFC
Exascale flow solver
Loading...
Searching...
No Matches
m_ib_patches.fpp.f90
Go to the documentation of this file.
1# 1 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
2!>
3!! @file
4!! @brief Contains module m_ib_patches
5
6# 1 "/home/runner/work/MFC/MFC/src/common/include/case.fpp" 1
7! This file exists so that Fypp can be run without generating case.fpp files for
8! each target. This is useful when generating documentation, for example. This
9! should also let MFC be built with CMake directly, without invoking mfc.sh.
10
11! For pre-process.
12# 8 "/home/runner/work/MFC/MFC/src/common/include/case.fpp"
13
14! For moving immersed boundaries in simulation
15# 12 "/home/runner/work/MFC/MFC/src/common/include/case.fpp"
16# 6 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp" 2
17# 1 "/home/runner/work/MFC/MFC/src/common/include/ExtrusionHardcodedIC.fpp" 1
18!> Allocate memory and read initial condition data for IC extrusion.
19!>
20!> @details
21!> This macro handles the complete initialization process for IC extrusion by:
22!>
23!> **Memory Allocation:**
24!> - stored_values(xRows, yRows, sys_size) - stores primitive variable data from files
25!> - x_coords(nrows) - stores x-coordinates from input files
26!> - y_coords(nrows) - stores y-coordinates from input files (3D case only)
27!>
28!> **File Reading Operations:**
29!> - Reads primitive variable data from multiple files with pattern:
30!> `prim.<file_number>.00.<file_extension>.dat`
31!> - Files are read from directory specified by `files_dir` parameter
32!> - Supports 1D, 2D, and 3D computational domains
33!>
34!> **Grid Structure Detection:**
35!> - 1D/2D: Counts lines in first file to determine xRows
36!> - 3D: Analyzes coordinate patterns to determine xRows and yRows structure
37!>
38!> **MPI Domain Mapping:**
39!> - Calculates global_offset_x and global_offset_y for MPI subdomain positioning
40!> - Maps file coordinates to local computational grid coordinates
41!>
42!> **Data Assignment:**
43!> - Populates q_prim_vf primitive variable arrays with file data
44!> - Handles momentum component indexing with special treatment for eqn_idx%mom%end
45!> - Sets eqn_idx%mom%end component to zero for 2D/3D cases
46!>
47!> **State Management:**
48!> - Uses files_loaded flag to prevent redundant file operations
49!> - Preserves data across multiple macro calls within same simulation
50!>
51!> @note File pattern timestep field is controlled by the `file_extension` parameter
52!> @note Directory path is set via the `files_dir` parameter
53!> @warning Aborts execution if file reading errors occur.
54
55# 67 "/home/runner/work/MFC/MFC/src/common/include/ExtrusionHardcodedIC.fpp"
56
57# 231 "/home/runner/work/MFC/MFC/src/common/include/ExtrusionHardcodedIC.fpp"
58
59# 250 "/home/runner/work/MFC/MFC/src/common/include/ExtrusionHardcodedIC.fpp"
60# 7 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp" 2
61# 1 "/home/runner/work/MFC/MFC/src/common/include/1dHardcodedIC.fpp" 1
62# 5 "/home/runner/work/MFC/MFC/src/common/include/1dHardcodedIC.fpp"
63
64# 72 "/home/runner/work/MFC/MFC/src/common/include/1dHardcodedIC.fpp"
65# 8 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp" 2
66# 1 "/home/runner/work/MFC/MFC/src/common/include/2dHardcodedIC.fpp" 1
67# 38 "/home/runner/work/MFC/MFC/src/common/include/2dHardcodedIC.fpp"
68
69# 569 "/home/runner/work/MFC/MFC/src/common/include/2dHardcodedIC.fpp"
70# 9 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp" 2
71# 1 "/home/runner/work/MFC/MFC/src/common/include/3dHardcodedIC.fpp" 1
72# 134 "/home/runner/work/MFC/MFC/src/common/include/3dHardcodedIC.fpp"
73
74# 296 "/home/runner/work/MFC/MFC/src/common/include/3dHardcodedIC.fpp"
75# 10 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp" 2
76# 1 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 1
77# 1 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 1
78# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
79# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
80# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
81# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
82# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
83# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
84
85# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
86# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
87# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
88
89# 12 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
90# 1 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp" 1
91! Default for Fypp runs outside the toolchain (e.g. documentation). Searched last, so the
92! build/include/<target>/thermochem.fpp the toolchain writes from m_thermochem takes precedence.
93# 4 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp"
94# 5 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp"
95# 13 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 2
96
97# 16 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
98# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
99
100# 24 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
101
102# 53 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
103
104# 65 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
105
106# 75 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
107
108# 105 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
109
110# 117 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
111
112# 127 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
113
114# 174 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
115! New line at end of file is required for FYPP
116# 2 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
117# 1 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 1
118# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
119# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
120# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
121# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
122# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
123# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
124
125# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
126# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
127# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
128
129# 12 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
130# 1 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp" 1
131! Default for Fypp runs outside the toolchain (e.g. documentation). Searched last, so the
132! build/include/<target>/thermochem.fpp the toolchain writes from m_thermochem takes precedence.
133# 4 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp"
134# 5 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp"
135# 13 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 2
136
137# 16 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
138# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
139
140# 24 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
141
142# 53 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
143
144# 65 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
145
146# 75 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
147
148# 105 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
149
150# 117 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
151
152# 127 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
153
154# 174 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
155! New line at end of file is required for FYPP
156# 2 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 2
157
158# 4 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
159# 5 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
160# 6 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
161# 7 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
162# 8 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
163
164# 20 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
165
166# 43 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
167
168# 48 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
169
170# 53 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
171
172# 58 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
173
174# 63 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
175
176# 68 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
177
178# 81 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
179
180# 95 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
181
182# 100 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
183
184# 105 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
185
186# 110 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
187
188# 115 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
189
190# 120 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
191
192# 125 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
193
194# 130 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
195
196# 135 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
197
198# 165 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
199
200# 219 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
201
202# 233 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
203
204# 258 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
205
206# 269 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
207
208# 271 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
209# 282 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
210
211# 332 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
212
213# 342 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
214
215# 352 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
216
217# 361 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
218
219# 378 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
220
221# 388 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
222
223# 395 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
224
225# 401 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
226
227# 407 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
228
229# 413 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
230
231# 419 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
232
233# 425 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
234! New line at end of file is required for FYPP
235# 3 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
236# 1 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 1
237# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
238# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
239# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
240# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
241# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
242# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
243
244# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
245# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
246# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
247
248# 12 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
249# 1 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp" 1
250! Default for Fypp runs outside the toolchain (e.g. documentation). Searched last, so the
251! build/include/<target>/thermochem.fpp the toolchain writes from m_thermochem takes precedence.
252# 4 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp"
253# 5 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp"
254# 13 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 2
255
256# 16 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
257# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
258
259# 24 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
260
261# 53 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
262
263# 65 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
264
265# 75 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
266
267# 105 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
268
269# 117 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
270
271# 127 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
272
273# 174 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
274! New line at end of file is required for FYPP
275# 2 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 2
276
277# 7 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
278
279# 17 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
280
281# 22 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
282
283# 27 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
284
285# 32 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
286
287# 37 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
288
289# 42 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
290
291# 47 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
292
293# 52 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
294
295# 57 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
296
297# 62 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
298
299# 73 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
300
301# 78 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
302
303# 83 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
304
305# 88 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
306
307# 103 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
308
309# 131 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
310
311# 160 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
312
313# 175 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
314
315# 193 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
316
317# 215 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
318
319# 244 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
320
321# 259 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
322
323# 269 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
324
325# 278 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
326
327# 294 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
328
329# 304 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
330
331# 311 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
332! New line at end of file is required for FYPP
333# 4 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
334
335! GPU parallel region (scalar reductions, maxval/minval)
336# 23 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
337
338! GPU parallel loop over threads (most common GPU macro)
339# 43 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
340
341! Required closing for GPU_PARALLEL_LOOP
342# 55 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
343
344! Mark routine for device compilation
345# 112 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
346
347! Declare device-resident data
348# 130 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
349
350! Inner loop within a GPU parallel region
351# 145 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
352
353! Scoped GPU data region
354# 164 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
355
356! Host code with device pointers (for MPI with GPU buffers)
357# 193 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
358
359! Allocate device memory (unscoped)
360# 207 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
361
362! Free device memory
363# 219 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
364
365! Atomic operation on device
366# 231 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
367
368! End atomic capture block
369# 242 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
370
371! Copy data between host and device
372# 254 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
373
374! Synchronization barrier
375# 266 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
376
377! Import GPU library module (openacc or omp_lib)
378# 275 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
379
380! Emit code only for AMD compiler
381# 282 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
382
383! Emit code for non-Cray compilers
384# 289 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
385
386! Emit code only for Cray compiler
387# 296 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
388
389! Emit code for non-NVIDIA compilers
390# 303 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
391
392# 305 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
393# 306 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
394! New line at end of file is required for FYPP
395# 2 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 2
396
397# 14 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
398
399! Caution: This macro requires the use of a binding script to set CUDA_VISIBLE_DEVICES, such that we have one GPU device per MPI
400! rank. That's because for both cudaMemAdvise (preferred location) and cudaMemPrefetchAsync we use location = device_id = 0. For an
401! example see misc/nvidia_uvm/bind.sh.
402# 52 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
403
404! Allocate and create GPU device memory
405# 72 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
406
407! Free GPU device memory and deallocate
408# 80 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
409
410! Cray-specific GPU pointer setup for vector fields
411# 104 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
412
413! Cray-specific GPU pointer setup for scalar fields
414# 120 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
415
416! Cray-specific GPU pointer setup for acoustic source spatials
417# 145 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
418
419# 151 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
420
421# 158 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
422! New line at end of file is required for FYPP
423# 11 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp" 2
424
425!> @brief Immersed boundary patch geometry constructors for 2D and 3D shapes
426module m_ib_patches
427
428 use m_patch_geometries
429 use m_model ! Subroutine(s) related to STL files
430 use m_derived_types ! Definitions of the derived types
431 use m_global_parameters
432 use m_helper_basic
433 use m_helper
434 use m_mpi_common
435 use m_constants
436
437 implicit none
438
439 private; public :: s_apply_ib_patches, s_update_ib_rotation_matrix, s_instantiate_stl_models, s_decode_patch_periodicity, &
440 & s_encode_patch_periodicity, s_initialize_ib_airfoils, s_get_periodicities, s_get_ib_bound, s_get_neighborhood_idx, &
441 & s_update_ib_lookup, s_compact_ib_lookup, s_merge_ib_lookup
442
443 !> lookup arrays for converting global IB indices to local indices
444 integer, dimension(num_ib_patches_max_namelist) :: ib_lookup_keys, ib_lookup_vals
445
446# 32 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
447#if defined(MFC_OpenACC)
448# 32 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
449!$acc declare create(ib_lookup_keys, ib_lookup_vals)
450# 32 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
451#elif defined(MFC_OpenMP)
452# 32 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
453!$omp declare target (ib_lookup_keys, ib_lookup_vals)
454# 32 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
455#endif
456# 33 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
457
458 !> Holds each step's arrivals while they are sorted and merged in. Host only.
459 integer, dimension(num_ib_patches_max_namelist) :: ib_new_keys, ib_new_vals
460
461contains
462
463 !> Apply all immersed boundary patch geometries to mark interior cells in the IB marker array
464 impure subroutine s_apply_ib_patches(ib_markers)
465
466 type(integer_field), intent(inout) :: ib_markers
467
468 if (many_ib_patch_parallelism) then
469 call s_apply_ib_patches_ib_parallelism(ib_markers)
470 else
471 call s_apply_ib_patches_grid_cell_parallelism(ib_markers)
472 end if
473
474 end subroutine s_apply_ib_patches
475
476 subroutine s_apply_ib_patches_grid_cell_parallelism(ib_markers)
477
478 type(integer_field), intent(inout) :: ib_markers
479 integer :: patch_id, airfoil_id, model_id, encoded_patch_id, i, j, k, il, ir, jl, jr, kl, kr, xp, yp, zp !< iterators
480 integer :: xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper !< periodic bounds
481 real(wp), dimension(3) :: center, xyz_local, length
482 real(wp) :: radius, eta
483
484
485# 60 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
486#if defined(MFC_OpenACC)
487# 60 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
488!$acc update host(patch_ib(1:num_ibs))
489# 60 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
490#elif defined(MFC_OpenMP)
491# 60 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
492!$omp target update from(patch_ib(1:num_ibs))
493# 60 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
494#endif
495# 61 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
496
497 ! 3D Patch Geometries
498 if (num_dims == 3) then
499 call s_get_periodicities(xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper)
500 do xp = xp_lower, xp_upper
501 do yp = yp_lower, yp_upper
502 do zp = zp_lower, zp_upper
503 do patch_id = 1, num_ibs
504 center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(glb_bounds(1)%end - glb_bounds(1)%beg)
505 center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(glb_bounds(2)%end - glb_bounds(2)%beg)
506 center(3) = patch_ib(patch_id)%z_centroid + real(zp, wp)*(glb_bounds(3)%end - glb_bounds(3)%beg)
507
508 ! encode the periodicity information into the patch_id
509 call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, zp, encoded_patch_id)
510
511 ! find the indices to the left and right of the IB in i, j, k
512 call s_get_bounding_indices(patch_ib(patch_id), center, il, ir, jl, jr, kl, kr)
513
514 ! skip patches whose bounding box does not overlap this rank's domain
515 if (ir < il .or. jr < jl .or. kr < kl) cycle
516
517
518# 82 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
519
520# 82 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
521#if defined(MFC_OpenACC)
522# 82 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
523!$acc parallel loop collapse(3) gang vector default(present) private(i, j, k, xyz_local, length, radius, airfoil_id, eta) copyin(patch_id, encoded_patch_id, center)
524# 82 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
525#elif defined(MFC_OpenMP)
526# 82 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
527
528# 82 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
529
530# 82 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
531
532# 82 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
533!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) &
534# 82 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
535!$omp& private(i, j, k, xyz_local, length, radius, airfoil_id, eta) map(to:patch_id, encoded_patch_id, center)
536# 82 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
537#endif
538# 84 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
539 do k = kl, kr
540 do j = jl, jr
541 do i = il, ir
542 ! get coordinate frame centered on IB
543 xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), z_cc(k) - center(3)]
544 ! rotate the frame into the IB's coordinates
545 xyz_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse, xyz_local)
546
547 ! perform the interior check for the patch geometry of this IB
548 if (patch_ib(patch_id)%geometry == 8) then
549 ! sphere geometry
550 radius = patch_ib(patch_id)%radius
551
552 if (f_is_inside_sphere(xyz_local(1), xyz_local(2), xyz_local(3), &
553 & radius)) ib_markers%sf(i, j, k) = encoded_patch_id
554 else if (patch_ib(patch_id)%geometry == 9) then
555 ! cuboid geometry
556 length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, &
557 & patch_ib(patch_id)%length_z]
558 if (f_is_inside_cuboid(xyz_local(1), xyz_local(2), xyz_local(3), &
559 & length)) ib_markers%sf(i, j, k) = encoded_patch_id
560 else if (patch_ib(patch_id)%geometry == 10) then
561 ! cylinder geometry
562 radius = patch_ib(patch_id)%radius
563 if (f_is_inside_cylinder(xyz_local(2), xyz_local(3), xyz_local(1), radius, &
564 & patch_ib(patch_id)%length_x)) ib_markers%sf(i, j, k) = encoded_patch_id
565 else if (patch_ib(patch_id)%geometry == 11) then
566 ! 3D airfoil geometry
567 airfoil_id = patch_ib(patch_id)%airfoil_id
568 xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset
569 if (f_is_inside_airfoil(xyz_local(1), xyz_local(2), xyz_local(3), &
570 & patch_ib(patch_id)%length_z, airfoil_id)) ib_markers%sf(i, j, &
571 & k) = encoded_patch_id
572 else if (patch_ib(patch_id)%geometry == 12) then
573 ! STL model geometry
574 xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset
575 model_id = patch_ib(patch_id)%model_id
576 eta = f_model_is_inside(gpu_ntrs(model_id), model_id, xyz_local)
577 if (eta > stl_models(model_id)%model_threshold) then
578 ib_markers%sf(i, j, k) = encoded_patch_id
579 end if
580 end if
581 end do
582 end do
583 end do
584
585# 129 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
586#if defined(MFC_OpenACC)
587# 129 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
588!$acc end parallel loop
589# 129 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
590#elif defined(MFC_OpenMP)
591# 129 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
592
593# 129 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
594!$omp end target teams loop
595# 129 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
596#endif
597# 130 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
598 end do
599 end do
600 end do
601 end do
602
603 ! 2D Patch Geometries
604 else if (num_dims == 2) then
605 call s_get_periodicities(xp_lower, xp_upper, yp_lower, yp_upper)
606 do xp = xp_lower, xp_upper
607 do yp = yp_lower, yp_upper
608 do patch_id = 1, num_ibs
609 center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(glb_bounds(1)%end - glb_bounds(1)%beg)
610 center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(glb_bounds(2)%end - glb_bounds(2)%beg)
611 center(3) = 0._wp
612
613 ! encode the periodicity information into the patch_id
614 call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, 0, encoded_patch_id)
615
616 ! find the indices to the left and right of the IB in i, j, k
617 call s_get_bounding_indices(patch_ib(patch_id), center, il, ir, jl, jr, kl, kr)
618
619 ! skip patches whose bounding box does not overlap this rank's domain
620 if (ir < il .or. jr < jl) cycle
621
622
623# 154 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
624
625# 154 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
626#if defined(MFC_OpenACC)
627# 154 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
628!$acc parallel loop collapse(2) gang vector default(present) private(i, j, xyz_local, airfoil_id, eta, length, radius) copyin(patch_id, encoded_patch_id, center)
629# 154 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
630#elif defined(MFC_OpenMP)
631# 154 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
632
633# 154 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
634
635# 154 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
636
637# 154 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
638!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(2) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) &
639# 154 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
640!$omp& private(i, j, xyz_local, airfoil_id, eta, length, radius) map(to:patch_id, encoded_patch_id, center)
641# 154 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
642#endif
643# 156 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
644 do j = jl, jr
645 do i = il, ir
646 ! get coordinate frame centered on IB
647 xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp]
648 ! rotate the frame into the IB's coordinates
649 xyz_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse, xyz_local)
650
651 ! perform the interior check for the patch geometry of this IB
652 if (patch_ib(patch_id)%geometry == 2) then
653 ! circular geometries
654 radius = patch_ib(patch_id)%radius
655 if (f_is_inside_cylinder(xyz_local(1), xyz_local(2), 0._wp, radius, 0._wp)) ib_markers%sf(i, &
656 & j, 0) = encoded_patch_id
657 else if (patch_ib(patch_id)%geometry == 3) then
658 ! rectangular geometries
659 length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, 0._wp]
660 if (f_is_inside_cuboid(xyz_local(1), xyz_local(2), xyz_local(3), length)) ib_markers%sf(i, j, &
661 & 0) = encoded_patch_id
662 else if (patch_ib(patch_id)%geometry == 4) then
663 ! 2D airfoil geometry
664 airfoil_id = patch_ib(patch_id)%airfoil_id
665 xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset
666 if (f_is_inside_airfoil(xyz_local(1), xyz_local(2), 0._wp, 0._wp, &
667 & airfoil_id)) ib_markers%sf(i, j, 0) = encoded_patch_id
668 else if (patch_ib(patch_id)%geometry == 5) then
669 ! STL model geometry
670 xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset
671 model_id = patch_ib(patch_id)%model_id
672 eta = f_model_is_inside(gpu_ntrs(model_id), model_id, xyz_local)
673 if (eta > stl_models(model_id)%model_threshold) then
674 ib_markers%sf(i, j, 0) = encoded_patch_id
675 end if
676 else if (patch_ib(patch_id)%geometry == 6) then
677 ! ellipse geometry
678 length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, 0._wp]
679 if (f_is_inside_ellipse(xyz_local(1), xyz_local(2), length)) ib_markers%sf(i, j, &
680 & 0) = encoded_patch_id
681 end if
682 end do
683 end do
684
685# 196 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
686#if defined(MFC_OpenACC)
687# 196 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
688!$acc end parallel loop
689# 196 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
690#elif defined(MFC_OpenMP)
691# 196 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
692
693# 196 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
694!$omp end target teams loop
695# 196 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
696#endif
697# 197 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
698 end do
699 end do
700 end do
701 end if
702
703 end subroutine s_apply_ib_patches_grid_cell_parallelism
704
705 subroutine s_apply_ib_patches_ib_parallelism(ib_markers)
706
707 type(integer_field), intent(inout) :: ib_markers
708 integer :: patch_id, airfoil_id, model_id, encoded_patch_id, i, j, k, il, ir, jl, jr, kl, kr, xp, yp, zp !< iterators
709 integer :: xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper !< periodic bounds
710 real(wp), dimension(3) :: center, xyz_local, length
711 real(wp) :: radius, eta
712
713 if (num_dims == 3) then
714 ! get the periodicities
715 call s_get_periodicities(xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper)
716
717 do xp = xp_lower, xp_upper
718 do yp = yp_lower, yp_upper
719 do zp = zp_lower, zp_upper
720
721# 219 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
722
723# 219 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
724#if defined(MFC_OpenACC)
725# 219 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
726!$acc parallel loop gang vector default(present) private(i, il, ir, j, jl, jr, k, kl, kr, xyz_local, length, radius, patch_id, airfoil_id, model_id, encoded_patch_id, center, eta) copyin(xp, yp, zp)
727# 219 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
728#elif defined(MFC_OpenMP)
729# 219 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
730
731# 219 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
732
733# 219 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
734
735# 219 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
736!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) &
737# 219 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
738!$omp& private(i, il, ir, j, jl, jr, k, kl, kr, xyz_local, length, radius, patch_id, airfoil_id, model_id, encoded_patch_id, center, eta) map(to:xp, yp, zp)
739# 219 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
740#endif
741# 221 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
742 do patch_id = 1, num_ibs
743 center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(glb_bounds(1)%end - glb_bounds(1)%beg)
744 center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(glb_bounds(2)%end - glb_bounds(2)%beg)
745 center(3) = patch_ib(patch_id)%z_centroid + real(zp, wp)*(glb_bounds(3)%end - glb_bounds(3)%beg)
746
747 ! encode the periodicity information into the patch_id
748 call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, zp, encoded_patch_id)
749
750 ! find the indices to the left and right of the IB in i, j, k
751 call s_get_bounding_indices(patch_ib(patch_id), center, il, ir, jl, jr, kl, kr)
752
753 do k = kl, kr
754 do j = jl, jr
755 do i = il, ir
756 ! get coordinate frame centered on IB
757 xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), z_cc(k) - center(3)]
758 ! rotate the frame into the IB's coordinates
759 xyz_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse, xyz_local)
760
761 ! perform the interior check for the patch geometry of this IB. Writes to ib_markers use
762 ! an atomic max (not a plain assignment) because this loop is parallel over patch_id: the
763 ! soft-sphere collision model allows particles to physically interpenetrate by design, so
764 ! two different patches can both claim the same overlapping cell here. A plain write would
765 ! be a data race with a nondeterministic winner; the atomic max makes the higher
766 ! encoded_patch_id win consistently every time, regardless of thread scheduling.
767 if (patch_ib(patch_id)%geometry == 8) then
768 ! sphere geometry
769 radius = patch_ib(patch_id)%radius
770
771 if (f_is_inside_sphere(xyz_local(1), xyz_local(2), xyz_local(3), radius)) then
772
773# 251 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
774#if defined(MFC_OpenACC)
775# 251 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
776!$acc atomic update
777# 251 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
778#elif defined(MFC_OpenMP)
779# 251 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
780!$omp atomic update
781# 251 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
782#endif
783# 252 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
784 ib_markers%sf(i, j, k) = max(ib_markers%sf(i, j, k), encoded_patch_id)
785 end if
786 else if (patch_ib(patch_id)%geometry == 9) then
787 ! cuboid geometry
788 length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, &
789 & patch_ib(patch_id)%length_z]
790 if (f_is_inside_cuboid(xyz_local(1), xyz_local(2), xyz_local(3), length)) then
791
792# 259 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
793#if defined(MFC_OpenACC)
794# 259 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
795!$acc atomic update
796# 259 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
797#elif defined(MFC_OpenMP)
798# 259 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
799!$omp atomic update
800# 259 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
801#endif
802# 260 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
803 ib_markers%sf(i, j, k) = max(ib_markers%sf(i, j, k), encoded_patch_id)
804 end if
805 else if (patch_ib(patch_id)%geometry == 10) then
806 ! cylinder geometry
807 radius = patch_ib(patch_id)%radius
808 if (f_is_inside_cylinder(xyz_local(2), xyz_local(3), xyz_local(1), radius, &
809 & patch_ib(patch_id)%length_x)) then
810
811# 267 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
812#if defined(MFC_OpenACC)
813# 267 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
814!$acc atomic update
815# 267 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
816#elif defined(MFC_OpenMP)
817# 267 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
818!$omp atomic update
819# 267 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
820#endif
821# 268 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
822 ib_markers%sf(i, j, k) = max(ib_markers%sf(i, j, k), encoded_patch_id)
823 end if
824 else if (patch_ib(patch_id)%geometry == 11) then
825 ! 3D airfoil geometry
826 airfoil_id = patch_ib(patch_id)%airfoil_id
827 xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset
828 if (f_is_inside_airfoil(xyz_local(1), xyz_local(2), xyz_local(3), &
829 & patch_ib(patch_id)%length_z, airfoil_id)) then
830
831# 276 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
832#if defined(MFC_OpenACC)
833# 276 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
834!$acc atomic update
835# 276 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
836#elif defined(MFC_OpenMP)
837# 276 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
838!$omp atomic update
839# 276 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
840#endif
841# 277 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
842 ib_markers%sf(i, j, k) = max(ib_markers%sf(i, j, k), encoded_patch_id)
843 end if
844 else if (patch_ib(patch_id)%geometry == 12) then
845 ! STL model geometry
846 xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset
847 model_id = patch_ib(patch_id)%model_id
848 eta = f_model_is_inside(gpu_ntrs(model_id), model_id, xyz_local)
849 if (eta > stl_models(model_id)%model_threshold) then
850
851# 285 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
852#if defined(MFC_OpenACC)
853# 285 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
854!$acc atomic update
855# 285 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
856#elif defined(MFC_OpenMP)
857# 285 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
858!$omp atomic update
859# 285 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
860#endif
861# 286 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
862 ib_markers%sf(i, j, k) = max(ib_markers%sf(i, j, k), encoded_patch_id)
863 end if
864 end if
865 end do
866 end do
867 end do
868 end do
869
870# 293 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
871#if defined(MFC_OpenACC)
872# 293 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
873!$acc end parallel loop
874# 293 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
875#elif defined(MFC_OpenMP)
876# 293 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
877
878# 293 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
879!$omp end target teams loop
880# 293 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
881#endif
882# 294 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
883 end do
884 end do
885 end do
886 else if (num_dims == 2) then
887 ! get the periodicities
888 call s_get_periodicities(xp_lower, xp_upper, yp_lower, yp_upper)
889
890 do xp = xp_lower, xp_upper
891 do yp = yp_lower, yp_upper
892
893# 303 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
894
895# 303 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
896#if defined(MFC_OpenACC)
897# 303 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
898!$acc parallel loop gang vector default(present) private(i, il, ir, j, jl, jr, xyz_local, length, radius, patch_id, airfoil_id, model_id, encoded_patch_id, center, eta) copyin(xp, yp)
899# 303 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
900#elif defined(MFC_OpenMP)
901# 303 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
902
903# 303 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
904
905# 303 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
906
907# 303 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
908!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) &
909# 303 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
910!$omp& private(i, il, ir, j, jl, jr, xyz_local, length, radius, patch_id, airfoil_id, model_id, encoded_patch_id, center, eta) map(to:xp, yp)
911# 303 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
912#endif
913# 305 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
914 do patch_id = 1, num_ibs
915 center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(glb_bounds(1)%end - glb_bounds(1)%beg)
916 center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(glb_bounds(2)%end - glb_bounds(2)%beg)
917 center(3) = 0._wp
918
919 ! encode the periodicity information into the patch_id
920 call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, 0, encoded_patch_id)
921
922 ! find the indices to the left and right of the IB in i, j, k
923 call s_get_bounding_indices(patch_ib(patch_id), center, il, ir, jl, jr, kl, kr)
924
925 do j = jl, jr
926 do i = il, ir
927 ! get coordinate frame centered on IB
928 xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp]
929 ! rotate the frame into the IB's coordinates
930 xyz_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse, xyz_local)
931
932 ! perform the interior check for the patch geometry of this IB. Writes to ib_markers use an
933 ! atomic max (not a plain assignment) because this loop is parallel over patch_id: the
934 ! soft-sphere collision model allows particles to physically interpenetrate by design, so two
935 ! different patches can both claim the same overlapping cell here. A plain write would be a data
936 ! race with a nondeterministic winner; the atomic max makes the higher encoded_patch_id win
937 ! consistently every time, regardless of thread scheduling.
938 if (patch_ib(patch_id)%geometry == 2) then
939 ! circular geometries
940 radius = patch_ib(patch_id)%radius
941 if (f_is_inside_cylinder(xyz_local(1), xyz_local(2), 0._wp, radius, 0._wp)) then
942
943# 333 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
944#if defined(MFC_OpenACC)
945# 333 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
946!$acc atomic update
947# 333 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
948#elif defined(MFC_OpenMP)
949# 333 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
950!$omp atomic update
951# 333 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
952#endif
953# 334 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
954 ib_markers%sf(i, j, 0) = max(ib_markers%sf(i, j, 0), encoded_patch_id)
955 end if
956 else if (patch_ib(patch_id)%geometry == 3) then
957 ! rectangular geometries
958 length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, 0._wp]
959 if (f_is_inside_cuboid(xyz_local(1), xyz_local(2), xyz_local(3), length)) then
960
961# 340 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
962#if defined(MFC_OpenACC)
963# 340 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
964!$acc atomic update
965# 340 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
966#elif defined(MFC_OpenMP)
967# 340 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
968!$omp atomic update
969# 340 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
970#endif
971# 341 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
972 ib_markers%sf(i, j, 0) = max(ib_markers%sf(i, j, 0), encoded_patch_id)
973 end if
974 else if (patch_ib(patch_id)%geometry == 4) then
975 ! 2D airfoil geometry
976 airfoil_id = patch_ib(patch_id)%airfoil_id
977 xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset
978 if (f_is_inside_airfoil(xyz_local(1), xyz_local(2), 0._wp, 0._wp, airfoil_id)) then
979
980# 348 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
981#if defined(MFC_OpenACC)
982# 348 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
983!$acc atomic update
984# 348 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
985#elif defined(MFC_OpenMP)
986# 348 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
987!$omp atomic update
988# 348 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
989#endif
990# 349 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
991 ib_markers%sf(i, j, 0) = max(ib_markers%sf(i, j, 0), encoded_patch_id)
992 end if
993 else if (patch_ib(patch_id)%geometry == 5) then
994 ! STL model geometry
995 xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset
996 model_id = patch_ib(patch_id)%model_id
997 eta = f_model_is_inside(gpu_ntrs(model_id), model_id, xyz_local)
998 if (eta > stl_models(model_id)%model_threshold) then
999
1000# 357 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1001#if defined(MFC_OpenACC)
1002# 357 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1003!$acc atomic update
1004# 357 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1005#elif defined(MFC_OpenMP)
1006# 357 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1007!$omp atomic update
1008# 357 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1009#endif
1010# 358 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1011 ib_markers%sf(i, j, 0) = max(ib_markers%sf(i, j, 0), encoded_patch_id)
1012 end if
1013 else if (patch_ib(patch_id)%geometry == 6) then
1014 ! ellipse geometry
1015 length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, 0._wp]
1016 if (f_is_inside_ellipse(xyz_local(1), xyz_local(2), length)) then
1017
1018# 364 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1019#if defined(MFC_OpenACC)
1020# 364 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1021!$acc atomic update
1022# 364 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1023#elif defined(MFC_OpenMP)
1024# 364 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1025!$omp atomic update
1026# 364 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1027#endif
1028# 365 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1029 ib_markers%sf(i, j, 0) = max(ib_markers%sf(i, j, 0), encoded_patch_id)
1030 end if
1031 end if
1032 end do
1033 end do
1034 end do
1035
1036# 371 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1037#if defined(MFC_OpenACC)
1038# 371 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1039!$acc end parallel loop
1040# 371 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1041#elif defined(MFC_OpenMP)
1042# 371 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1043
1044# 371 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1045!$omp end target teams loop
1046# 371 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1047#endif
1048# 372 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1049 end do
1050 end do
1051 end if
1052
1053 end subroutine s_apply_ib_patches_ib_parallelism
1054
1055 !> Initialize the NACA surface grids for all airfoil IB patches. Must be called after the grid is established (so dx is valid)
1056 !! and before s_apply_ib_patches or s_apply_levelset.
1057 subroutine s_initialize_ib_airfoils()
1058
1059 integer :: i, j, airfoil_id, np
1060 real(wp) :: ca_in, pa, ma, ta
1061 real(wp) :: xa, yc, dycdxc, yt, xu, yu, xl, yl, sin_c, cos_c
1062 real(wp) :: edge(1:2) !< segment vector from vertex 1 to vertex 2
1063 real(wp) :: min_grid_spacing !< smallest cell width across all ranks
1064
1065 if (.not. any(patch_ib(1:num_ibs)%geometry == 4 .or. patch_ib(1:num_ibs)%geometry == 11)) return
1066
1067 ! size the surface on the global minimum spacing so every rank builds the same airfoil
1068 min_grid_spacing = min(dx_min, dy_min)
1069 if (num_procs > 1) call s_mpi_allreduce_min(min(dx_min, dy_min), min_grid_spacing)
1070
1071 do i = 1, num_ibs
1072 if (patch_ib(i)%geometry /= 4 .and. patch_ib(i)%geometry /= 11) cycle
1073
1074 airfoil_id = patch_ib(i)%airfoil_id
1075 ca_in = ib_airfoil(airfoil_id)%c
1076 pa = ib_airfoil(airfoil_id)%p
1077 ma = ib_airfoil(airfoil_id)%m
1078 ta = ib_airfoil(airfoil_id)%t
1079
1080 np = int(10._wp*ca_in/min_grid_spacing) + 1 ! 20 surface points per smallest cell
1081 ib_airfoil_grids(airfoil_id)%Np = np
1082
1083# 405 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1084#if defined(MFC_OpenACC)
1085# 405 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1086!$acc update device(ib_airfoil_grids(airfoil_id)%Np)
1087# 405 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1088#elif defined(MFC_OpenMP)
1089# 405 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1090!$omp target update to(ib_airfoil_grids(airfoil_id)%Np)
1091# 405 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1092#endif
1093# 406 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1094
1095 if (.not. allocated(ib_airfoil_grids(airfoil_id)%upper)) then
1096#ifdef MFC_DEBUG
1097# 408 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1098 block
1099# 408 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1100 use iso_fortran_env, only: output_unit
1101# 408 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1102
1103# 408 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1104 print *, 'm_ib_patches.fpp:408: ', '@:ALLOCATE(ib_airfoil_grids(airfoil_id)%upper(1:Np - 1, 1:3, 1:2))'
1105# 408 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1106
1107# 408 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1108 call flush (output_unit)
1109# 408 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1110 end block
1111# 408 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1112#endif
1113# 408 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1114 allocate (ib_airfoil_grids(airfoil_id)%upper(1:np - 1, 1:3, 1:2))
1115# 408 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1116
1117# 408 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1118
1119# 408 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1120#if defined(MFC_OpenACC)
1121# 408 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1122!$acc enter data create(ib_airfoil_grids(airfoil_id)%upper)
1123# 408 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1124#elif defined(MFC_OpenMP)
1125# 408 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1126!$omp target enter data map(always,alloc:ib_airfoil_grids(airfoil_id)%upper)
1127# 408 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1128#endif
1129# 409 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1130#ifdef MFC_DEBUG
1131# 409 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1132 block
1133# 409 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1134 use iso_fortran_env, only: output_unit
1135# 409 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1136
1137# 409 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1138 print *, 'm_ib_patches.fpp:409: ', '@:ALLOCATE(ib_airfoil_grids(airfoil_id)%lower(1:Np - 1, 1:3, 1:2))'
1139# 409 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1140
1141# 409 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1142 call flush (output_unit)
1143# 409 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1144 end block
1145# 409 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1146#endif
1147# 409 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1148 allocate (ib_airfoil_grids(airfoil_id)%lower(1:np - 1, 1:3, 1:2))
1149# 409 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1150
1151# 409 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1152
1153# 409 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1154#if defined(MFC_OpenACC)
1155# 409 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1156!$acc enter data create(ib_airfoil_grids(airfoil_id)%lower)
1157# 409 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1158#elif defined(MFC_OpenMP)
1159# 409 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1160!$omp target enter data map(always,alloc:ib_airfoil_grids(airfoil_id)%lower)
1161# 409 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1162#endif
1163# 410 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1164
1165 ! segments run from the leading edge to the trailing edge
1166 ib_airfoil_grids(airfoil_id)%upper(1, 1,:) = 0._wp
1167 ib_airfoil_grids(airfoil_id)%lower(1, 1,:) = 0._wp
1168
1169 do j = 1, np - 2
1170 ! cosine spacing clusters points at the leading and trailing edges, Katz & Plotkin (2001)
1171 ! "Low-Speed Aerodynamics", 2nd ed., Cambridge University Press
1172 xa = 0.5_wp*(1._wp - cos(pi*real(j, wp)/real(np - 1, wp)))
1173 if (xa <= pa) then
1174 yc = (ma/pa**2)*(2*pa*xa - xa**2)
1175 dycdxc = (2*ma/pa**2)*(pa - xa)
1176 else
1177 yc = (ma/(1 - pa)**2)*(1 - 2*pa + 2*pa*xa - xa**2)
1178 dycdxc = (2*ma/(1 - pa)**2)*(pa - xa)
1179 end if
1180
1181 yt = (5._wp*ta)*(0.2969_wp*xa**0.5_wp - 0.126_wp*xa - 0.3516_wp*xa**2._wp + 0.2843_wp*xa**3 - 0.1015_wp*xa**4)
1182 sin_c = dycdxc/(1 + dycdxc**2)**0.5_wp
1183 cos_c = 1/(1 + dycdxc**2)**0.5_wp
1184
1185 xu = (xa - yt*sin_c)*ca_in
1186 yu = (yc + yt*cos_c)*ca_in
1187 xl = (xa + yt*sin_c)*ca_in
1188 yl = (yc - yt*cos_c)*ca_in
1189
1190 ! each surface point ends segment j and starts segment j + 1
1191 ib_airfoil_grids(airfoil_id)%upper(j, 2,:) = [xu, yu]
1192 ib_airfoil_grids(airfoil_id)%upper(j + 1, 1,:) = [xu, yu]
1193 ib_airfoil_grids(airfoil_id)%lower(j, 2,:) = [xl, yl]
1194 ib_airfoil_grids(airfoil_id)%lower(j + 1, 1,:) = [xl, yl]
1195 end do
1196
1197 ib_airfoil_grids(airfoil_id)%upper(np - 1, 2,:) = [ca_in, 0._wp]
1198 ib_airfoil_grids(airfoil_id)%lower(np - 1, 2,:) = [ca_in, 0._wp]
1199
1200 ! outward unit normals: left of the edge on the upper surface, right of it on the lower
1201 do j = 1, np - 1
1202 edge = ib_airfoil_grids(airfoil_id)%upper(j, 2,:) - ib_airfoil_grids(airfoil_id)%upper(j, 1,:)
1203 ib_airfoil_grids(airfoil_id)%upper(j, 3,:) = [-edge(2), edge(1)]/norm2(edge)
1204 edge = ib_airfoil_grids(airfoil_id)%lower(j, 2,:) - ib_airfoil_grids(airfoil_id)%lower(j, 1,:)
1205 ib_airfoil_grids(airfoil_id)%lower(j, 3,:) = [edge(2), -edge(1)]/norm2(edge)
1206 end do
1207
1208
1209# 454 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1210#if defined(MFC_OpenACC)
1211# 454 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1212!$acc update device(ib_airfoil_grids(airfoil_id)%upper, ib_airfoil_grids(airfoil_id)%lower)
1213# 454 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1214#elif defined(MFC_OpenMP)
1215# 454 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1216!$omp target update to(ib_airfoil_grids(airfoil_id)%upper, ib_airfoil_grids(airfoil_id)%lower)
1217# 454 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1218#endif
1219# 455 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1220 end if
1221
1222
1223# 457 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1224#if defined(MFC_OpenACC)
1225# 457 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1226!$acc update device(ib_airfoil(airfoil_id))
1227# 457 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1228#elif defined(MFC_OpenMP)
1229# 457 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1230!$omp target update to(ib_airfoil(airfoil_id))
1231# 457 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1232#endif
1233# 458 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1234 end do
1235
1236 end subroutine s_initialize_ib_airfoils
1237
1238 !> Compute a rotation matrix for converting to the rotating frame of the boundary
1239 subroutine s_update_ib_rotation_matrix(patch_id)
1240
1241
1242# 465 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1243#if MFC_OpenACC
1244# 465 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1245!$acc routine seq
1246# 465 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1247#elif MFC_OpenMP
1248# 465 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1249
1250# 465 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1251
1252# 465 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1253!$omp declare target device_type(any)
1254# 465 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1255#endif
1256# 466 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1257
1258 integer, intent(in) :: patch_id
1259 real(wp), dimension(3, 3, 3) :: rotation
1260 real(wp) :: angle
1261
1262 ! construct the x, y, and z rotation matrices
1263
1264 if (num_dims == 3) then
1265 ! also compute the x and y axes in 3D
1266 angle = patch_ib(patch_id)%angles(1)
1267 rotation(1, 1,:) = [1._wp, 0._wp, 0._wp]
1268 rotation(1, 2,:) = [0._wp, cos(angle), -sin(angle)]
1269 rotation(1, 3,:) = [0._wp, sin(angle), cos(angle)]
1270
1271 angle = patch_ib(patch_id)%angles(2)
1272 rotation(2, 1,:) = [cos(angle), 0._wp, sin(angle)]
1273 rotation(2, 2,:) = [0._wp, 1._wp, 0._wp]
1274 rotation(2, 3,:) = [-sin(angle), 0._wp, cos(angle)]
1275
1276 ! apply the y rotation to the x rotation
1277 patch_ib(patch_id)%rotation_matrix(:,:) = matmul(rotation(1,:,:), rotation(2,:,:))
1278 patch_ib(patch_id)%rotation_matrix_inverse(:,:) = matmul(transpose(rotation(2,:,:)), transpose(rotation(1,:,:)))
1279 end if
1280
1281 ! z component first, since it applies in 2D and 3D
1282 angle = patch_ib(patch_id)%angles(3)
1283 rotation(3, 1,:) = [cos(angle), -sin(angle), 0._wp]
1284 rotation(3, 2,:) = [sin(angle), cos(angle), 0._wp]
1285 rotation(3, 3,:) = [0._wp, 0._wp, 1._wp]
1286
1287 if (num_dims == 3) then
1288 ! apply the z rotation to the xy rotation in 3D
1289 patch_ib(patch_id)%rotation_matrix(:,:) = matmul(patch_ib(patch_id)%rotation_matrix(:,:), rotation(3,:,:))
1290 patch_ib(patch_id)%rotation_matrix_inverse(:,:) = matmul(transpose(rotation(3,:,:)), &
1291 & patch_ib(patch_id)%rotation_matrix_inverse(:,:))
1292 else
1293 ! write out only the z rotation in 2D
1294 patch_ib(patch_id)%rotation_matrix(:,:) = rotation(3,:,:)
1295 patch_ib(patch_id)%rotation_matrix_inverse(:,:) = transpose(rotation(3,:,:))
1296 end if
1297
1298 end subroutine s_update_ib_rotation_matrix
1299
1300 subroutine s_get_ib_bound(patch, bound)
1301
1302
1303# 511 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1304#if MFC_OpenACC
1305# 511 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1306!$acc routine seq
1307# 511 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1308#elif MFC_OpenMP
1309# 511 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1310
1311# 511 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1312
1313# 511 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1314!$omp declare target device_type(any)
1315# 511 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1316#endif
1317# 512 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1318
1319 type(ib_patch_parameters), intent(in) :: patch
1320 real(wp), intent(out) :: bound
1321 real(wp), dimension(2) :: lx, ly, lz
1322
1323 if (patch%geometry == 2 .or. patch%geometry == 8) then
1324 ! circle and sphere geometries
1325 bound = patch%radius
1326 else if (patch%geometry == 3) then
1327 bound = 0.5_wp*sqrt(patch%length_x**2 + patch%length_y**2)
1328 else if (patch%geometry == 4 .or. patch%geometry == 11) then
1329 ! rectangular geometries
1330 bound = ib_airfoil(patch%airfoil_id)%c
1331 else if (patch%geometry == 5) then
1332 ! STL model geometry
1333 lx(1) = stl_bounding_boxes(patch%model_id, 1, 1)
1334 lx(2) = stl_bounding_boxes(patch%model_id, 1, 3)
1335 ly(1) = stl_bounding_boxes(patch%model_id, 2, 1)
1336 ly(2) = stl_bounding_boxes(patch%model_id, 2, 3)
1337
1338 bound = 0.5_wp*sqrt((lx(2) - lx(1))**2 + (ly(2) - ly(1))**2)
1339 else if (patch%geometry == 6) then
1340 ! ellipse geometry
1341 bound = 0.5_wp*max(patch%length_x, patch%length_y)
1342 else if (patch%geometry == 9) then
1343 ! cuboid geometries
1344 bound = 0.5_wp*sqrt(patch%length_x**2 + patch%length_y**2 + patch%length_z**2)
1345 else if (patch%geometry == 10) then
1346 ! cylinder geometry
1347 bound = sqrt(patch%radius**2 + patch%length_x**2)
1348 else if (patch%geometry == 12) then
1349 ! Local-space bounding box extents (min=1, max=2 in the third index)
1350 lx(1) = stl_bounding_boxes(patch%model_id, 1, 1) + patch%centroid_offset(1)
1351 lx(2) = stl_bounding_boxes(patch%model_id, 1, 3) + patch%centroid_offset(1)
1352 ly(1) = stl_bounding_boxes(patch%model_id, 2, 1) + patch%centroid_offset(2)
1353 ly(2) = stl_bounding_boxes(patch%model_id, 2, 3) + patch%centroid_offset(2)
1354 lz(1) = stl_bounding_boxes(patch%model_id, 3, 1) + patch%centroid_offset(3)
1355 lz(2) = stl_bounding_boxes(patch%model_id, 3, 3) + patch%centroid_offset(3)
1356
1357 bound = 0.5_wp*sqrt((lx(2) - lx(1))**2 + (ly(2) - ly(1))**2 + (lz(2) - lz(1))**2)
1358 end if
1359
1360 end subroutine s_get_ib_bound
1361
1362 subroutine s_get_bounding_indices(patch, center, il, ir, jl, jr, kl, kr)
1363
1364
1365# 558 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1366#if MFC_OpenACC
1367# 558 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1368!$acc routine seq
1369# 558 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1370#elif MFC_OpenMP
1371# 558 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1372
1373# 558 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1374
1375# 558 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1376!$omp declare target device_type(any)
1377# 558 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1378#endif
1379# 559 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1380
1381 type(ib_patch_parameters), intent(in) :: patch
1382 real(wp), dimension(3), intent(in) :: center
1383 integer, intent(out) :: il, ir, jl, jr, kl, kr
1384 real(wp), dimension(3) :: bbox_min, bbox_max, local_corner, world_corner
1385 real(wp), dimension(2) :: lx, ly, lz
1386 real(wp) :: bound
1387 integer :: cx, cy, cz
1388 logical :: outside_domain
1389
1390 if (patch%geometry == 5) then
1391 ! STL model geometry
1392 lx(1) = stl_bounding_boxes(patch%model_id, 1, 1) + patch%centroid_offset(1)
1393 lx(2) = stl_bounding_boxes(patch%model_id, 1, 3) + patch%centroid_offset(1)
1394 ly(1) = stl_bounding_boxes(patch%model_id, 2, 1) + patch%centroid_offset(2)
1395 ly(2) = stl_bounding_boxes(patch%model_id, 2, 3) + patch%centroid_offset(2)
1396
1397 bbox_min = 1e12
1398 bbox_max = -1e12
1399 ! Enumerate all 4 corners of the local bounding box, rotate to world space, track world-space AABB
1400 do cx = 1, 2
1401 do cy = 1, 2
1402 local_corner = [lx(cx), ly(cy), 0._wp]
1403 world_corner = matmul(patch%rotation_matrix, local_corner) + center
1404 bbox_min(1) = min(bbox_min(1), world_corner(1))
1405 bbox_min(2) = min(bbox_min(2), world_corner(2))
1406 bbox_max(1) = max(bbox_max(1), world_corner(1))
1407 bbox_max(2) = max(bbox_max(2), world_corner(2))
1408 end do
1409 end do
1410 else if (patch%geometry == 12) then
1411 ! Local-space bounding box extents (min=1, max=2 in the third index)
1412 lx(1) = stl_bounding_boxes(patch%model_id, 1, 1) + patch%centroid_offset(1)
1413 lx(2) = stl_bounding_boxes(patch%model_id, 1, 3) + patch%centroid_offset(1)
1414 ly(1) = stl_bounding_boxes(patch%model_id, 2, 1) + patch%centroid_offset(2)
1415 ly(2) = stl_bounding_boxes(patch%model_id, 2, 3) + patch%centroid_offset(2)
1416 lz(1) = stl_bounding_boxes(patch%model_id, 3, 1) + patch%centroid_offset(3)
1417 lz(2) = stl_bounding_boxes(patch%model_id, 3, 3) + patch%centroid_offset(3)
1418
1419 bbox_min = 1e12
1420 bbox_max = -1e12
1421 ! Enumerate all 8 corners of the local bounding box, rotate to world space, track world-space AABB
1422 do cx = 1, 2
1423 do cy = 1, 2
1424 do cz = 1, 2
1425 local_corner = [lx(cx), ly(cy), lz(cz)]
1426 world_corner = matmul(patch%rotation_matrix, local_corner) + center
1427 bbox_min(1) = min(bbox_min(1), world_corner(1))
1428 bbox_min(2) = min(bbox_min(2), world_corner(2))
1429 bbox_min(3) = min(bbox_min(3), world_corner(3))
1430 bbox_max(1) = max(bbox_max(1), world_corner(1))
1431 bbox_max(2) = max(bbox_max(2), world_corner(2))
1432 bbox_max(3) = max(bbox_max(3), world_corner(3))
1433 end do
1434 end do
1435 end do
1436 else
1437 ! All other IBs
1438 call s_get_ib_bound(patch, bound)
1439 bbox_min = center - bound
1440 bbox_max = center + bound
1441 end if
1442
1443 ! completely skip patches whose bounding box does not overlap this rank's domain
1444 ! Markers cover the full halo: image-point stencils of ghost points near a rank boundary read markers up to
1445 ! 2*gp_layers+1 cells into it, and an undrawn cell there reads as fluid.
1446 outside_domain = bbox_min(1) > x_cc(m + buff_size) .or. bbox_max(1) < x_cc(-buff_size) .or. bbox_min(2) > y_cc(n &
1447 & + buff_size) .or. bbox_max(2) < y_cc(-buff_size)
1448 if (num_dims == 3) then
1449 outside_domain = outside_domain .or. bbox_min(3) > z_cc(p + buff_size) .or. bbox_max(3) < z_cc(-buff_size)
1450 end if
1451
1452 if (outside_domain) then
1453 il = 1; ir = 0
1454 jl = 1; jr = 0
1455 kl = 1; kr = 0
1456 return
1457 end if
1458
1459 il = -buff_size
1460 jl = -buff_size
1461 kl = -buff_size
1462 ir = m + buff_size
1463 jr = n + buff_size
1464 kr = p + buff_size
1465 call get_indices_from_bounds(bbox_min(1), bbox_max(1), x_cc, il, ir)
1466 call get_indices_from_bounds(bbox_min(2), bbox_max(2), y_cc, jl, jr)
1467 if (num_dims == 3) call get_indices_from_bounds(bbox_min(3), bbox_max(3), z_cc, kl, kr)
1468
1469 end subroutine s_get_bounding_indices
1470
1471 subroutine get_indices_from_bounds(left_bound, right_bound, cell_centers, left_index, right_index)
1472
1473
1474# 652 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1475#if MFC_OpenACC
1476# 652 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1477!$acc routine seq
1478# 652 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1479#elif MFC_OpenMP
1480# 652 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1481
1482# 652 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1483
1484# 652 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1485!$omp declare target device_type(any)
1486# 652 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1487#endif
1488# 653 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1489
1490 real(wp), intent(in) :: left_bound, right_bound
1491 integer, intent(inout) :: left_index, right_index
1492 real(wp), dimension(-buff_size:), intent(in) :: cell_centers
1493 integer :: itr_left, itr_middle, itr_right
1494
1495 itr_left = left_index
1496 itr_right = right_index
1497
1498 do while (itr_left + 1 < itr_right)
1499 itr_middle = (itr_left + itr_right)/2
1500 if (cell_centers(itr_middle) < left_bound) then
1501 itr_left = itr_middle
1502 else if (cell_centers(itr_middle) > left_bound) then
1503 itr_right = itr_middle
1504 else
1505 itr_left = itr_middle
1506 exit
1507 end if
1508 end do
1509 left_index = itr_left
1510
1511 itr_right = right_index
1512 do while (itr_left + 1 < itr_right)
1513 itr_middle = (itr_left + itr_right)/2
1514 if (cell_centers(itr_middle) < right_bound) then
1515 itr_left = itr_middle
1516 else if (cell_centers(itr_middle) > right_bound) then
1517 itr_right = itr_middle
1518 else
1519 itr_right = itr_middle
1520 exit
1521 end if
1522 end do
1523 right_index = itr_right
1524
1525 end subroutine get_indices_from_bounds
1526
1527 !> Encode the patch ID with a unique offset containing periodicity information
1528 subroutine s_encode_patch_periodicity(patch_id, x_periodicity, y_periodicity, z_periodicity, encoded_patch_id)
1529
1530
1531# 694 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1532#if MFC_OpenACC
1533# 694 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1534!$acc routine seq
1535# 694 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1536#elif MFC_OpenMP
1537# 694 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1538
1539# 694 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1540
1541# 694 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1542!$omp declare target device_type(any)
1543# 694 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1544#endif
1545# 695 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1546
1547 integer, intent(in) :: patch_id, x_periodicity, y_periodicity, z_periodicity
1548 integer, intent(out) :: encoded_patch_id
1549 integer :: temp_x_per, temp_y_per, temp_z_per, offset
1550
1551 encoded_patch_id = patch_id
1552
1553 temp_x_per = x_periodicity; if (x_periodicity == -1) temp_x_per = 2
1554 temp_y_per = y_periodicity; if (y_periodicity == -1) temp_y_per = 2
1555 temp_z_per = z_periodicity; if (z_periodicity == -1) temp_z_per = 2
1556
1557 offset = (num_gbl_ibs + 1)*temp_x_per + 3*(num_gbl_ibs + 1)*temp_y_per + 9*(num_gbl_ibs + 1)*temp_z_per
1558 encoded_patch_id = patch_id + offset
1559
1560 end subroutine s_encode_patch_periodicity
1561
1562 !> Decode the encoded ID to recover the original patch ID and periodicity
1563 subroutine s_decode_patch_periodicity(encoded_patch_id, patch_id, x_periodicity, y_periodicity, z_periodicity)
1564
1565
1566# 714 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1567#if MFC_OpenACC
1568# 714 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1569!$acc routine seq
1570# 714 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1571#elif MFC_OpenMP
1572# 714 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1573
1574# 714 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1575
1576# 714 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1577!$omp declare target device_type(any)
1578# 714 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1579#endif
1580# 715 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1581
1582 integer, intent(in) :: encoded_patch_id
1583 integer, intent(out) :: patch_id
1584 integer, intent(out), optional :: x_periodicity, y_periodicity, z_periodicity
1585 integer :: offset, remainder, xp, yp, zp, base
1586
1587 base = num_gbl_ibs + 1
1588
1589 patch_id = mod(encoded_patch_id - 1, base) + 1
1590 offset = (encoded_patch_id - patch_id)/base
1591
1592 xp = mod(offset, 3)
1593 remainder = offset/3
1594 yp = mod(remainder, 3)
1595 zp = remainder/3
1596
1597 ! Reverse map: 2 -> -1, 0 -> 0, 1 -> 1
1598 if (present(x_periodicity) .and. present(y_periodicity) .and. present(z_periodicity)) then
1599 x_periodicity = xp; if (xp == 2) x_periodicity = -1
1600 y_periodicity = yp; if (yp == 2) y_periodicity = -1
1601 z_periodicity = zp; if (zp == 2) z_periodicity = -1
1602 end if
1603
1604 end subroutine s_decode_patch_periodicity
1605
1606 !> binary search to retrieve the local IB patch index using the global index
1607 subroutine s_get_neighborhood_idx(gbl_idx, neighborhood_idx, num_entries)
1608
1609
1610# 743 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1611#if MFC_OpenACC
1612# 743 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1613!$acc routine seq
1614# 743 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1615#elif MFC_OpenMP
1616# 743 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1617
1618# 743 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1619
1620# 743 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1621!$omp declare target device_type(any)
1622# 743 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1623#endif
1624# 744 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1625
1626 integer, intent(in) :: gbl_idx
1627 integer, intent(out) :: neighborhood_idx
1628 integer, intent(in), optional :: num_entries
1629 integer :: lo, hi, mid
1630
1631 neighborhood_idx = -1
1632 lo = 1
1633 hi = num_ibs
1634 if (present(num_entries)) hi = num_entries
1635
1636 do while (lo <= hi)
1637 mid = lo + (hi - lo)/2
1638 if (ib_lookup_keys(mid) == gbl_idx) then
1639 neighborhood_idx = ib_lookup_vals(mid)
1640 return
1641 else if (ib_lookup_keys(mid) < gbl_idx) then
1642 lo = mid + 1
1643 else
1644 hi = mid - 1
1645 end if
1646 end do
1647
1648 end subroutine s_get_neighborhood_idx
1649
1650 !> Completely rebuilds the ib lookup map, used at startup
1651 subroutine s_update_ib_lookup()
1652
1653 integer :: i
1654
1655 if (num_ibs > num_ib_patches_max_namelist) then
1656# 774 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1657 call s_prohibit_abort("num_ibs > num_ib_patches_max_namelist", "num_ibs exceeds the IB lookup capacity. Increase num_ib_patches_max_namelist.")
1658# 774 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1659 end if
1660# 776 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1661
1662 do i = 1, num_ibs
1663 ib_lookup_keys(i) = patch_ib(i)%gbl_patch_id
1664 ib_lookup_vals(i) = i
1665 end do
1666 call s_sort_int_key_value(ib_lookup_keys, ib_lookup_vals, num_ibs)
1667
1668
1669# 783 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1670#if defined(MFC_OpenACC)
1671# 783 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1672!$acc update device(ib_lookup_keys(1:num_ibs), ib_lookup_vals(1:num_ibs))
1673# 783 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1674#elif defined(MFC_OpenMP)
1675# 783 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1676!$omp target update to(ib_lookup_keys(1:num_ibs), ib_lookup_vals(1:num_ibs))
1677# 783 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1678#endif
1679# 784 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1680
1681 end subroutine s_update_ib_lookup
1682
1683 !> Drop the entries whose patches left the neighborhood and renumber the survivors onto the patch_ib slots they were compacted
1684 !! into. Keys are never reordered, so the map stays sorted for free and only the values move: O(num_ibs_old) against re-sorting
1685 !! the whole map.
1686 subroutine s_compact_ib_lookup(old_to_new, num_ibs_old)
1687
1688 integer, dimension(:), intent(in) :: old_to_new !< old patch_ib slot -> new slot, -1 if dropped
1689 integer, intent(in) :: num_ibs_old
1690 integer :: i, k
1691
1692 k = 0
1693 do i = 1, num_ibs_old
1694 if (old_to_new(ib_lookup_vals(i)) < 0) cycle
1695 k = k + 1
1696 ib_lookup_keys(k) = ib_lookup_keys(i)
1697 ib_lookup_vals(k) = old_to_new(ib_lookup_vals(i))
1698 end do
1699 if (.not. (k == num_ibs)) then
1700# 803 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1701 call s_mpi_abort("m_ib_patches.fpp:803: " // "Assertion failed: k == num_ibs. " &
1702# 803 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1703 & // 'IB lookup and patch_ib disagree on the surviving patch count')
1704# 803 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1705 end if
1706# 804 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1707
1708
1709# 805 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1710#if defined(MFC_OpenACC)
1711# 805 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1712!$acc update device(ib_lookup_keys(1:num_ibs), ib_lookup_vals(1:num_ibs))
1713# 805 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1714#elif defined(MFC_OpenMP)
1715# 805 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1716!$omp target update to(ib_lookup_keys(1:num_ibs), ib_lookup_vals(1:num_ibs))
1717# 805 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1718#endif
1719# 806 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1720
1721 end subroutine s_compact_ib_lookup
1722
1723 !> Fold patch_ib(num_ibs_pre+1:num_ibs) into the map: sort just the arrivals, then merge the two sorted runs downward from the
1724 !! top. O(num_ibs) plus the sort of the few arrivals. They are copied out first because the runs share this array and merging in
1725 !! place would overwrite entries still to be read.
1726 subroutine s_merge_ib_lookup(num_ibs_pre)
1727
1728 integer, intent(in) :: num_ibs_pre
1729 integer :: i, j, k, r
1730 logical :: take_old
1731
1732 call nvtxstartrange("MERGE-IB-LOOKUP")
1733
1734 r = num_ibs - num_ibs_pre
1735 if (r <= 0) return
1736
1737 do i = 1, r
1738 ib_new_keys(i) = patch_ib(num_ibs_pre + i)%gbl_patch_id
1739 ib_new_vals(i) = num_ibs_pre + i
1740 end do
1741 call s_sort_int_key_value(ib_new_keys, ib_new_vals, r)
1742
1743 i = num_ibs_pre; j = r; k = num_ibs
1744 do while (j >= 1)
1745 ! Fortran does not short-circuit .and., so the exhausted-head test stands on its own
1746 take_old = .false.
1747 if (i >= 1) take_old = ib_lookup_keys(i) > ib_new_keys(j)
1748
1749 if (take_old) then
1750 ib_lookup_keys(k) = ib_lookup_keys(i); ib_lookup_vals(k) = ib_lookup_vals(i)
1751 i = i - 1
1752 else
1753 ib_lookup_keys(k) = ib_new_keys(j); ib_lookup_vals(k) = ib_new_vals(j)
1754 j = j - 1
1755 end if
1756 k = k - 1
1757 end do
1758
1759
1760# 845 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1761#if defined(MFC_OpenACC)
1762# 845 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1763!$acc update device(ib_lookup_keys(1:num_ibs), ib_lookup_vals(1:num_ibs))
1764# 845 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1765#elif defined(MFC_OpenMP)
1766# 845 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1767!$omp target update to(ib_lookup_keys(1:num_ibs), ib_lookup_vals(1:num_ibs))
1768# 845 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1769#endif
1770# 846 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1771
1772 call nvtxendrange()
1773
1774 end subroutine s_merge_ib_lookup
1775
1776 !> Determine the periodic wrapping bounds in each direction
1777 subroutine s_get_periodicities(xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper)
1778
1779 integer, intent(out) :: xp_lower, xp_upper, yp_lower, yp_upper
1780 integer, intent(out), optional :: zp_lower, zp_upper
1781
1782 ! check domain wraps in x, y
1783
1784# 860 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1785 ! check for periodicity
1786 if (ib_bc_x%beg == bc_periodic) then
1787# 862 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1788 xp_lower = -1
1789# 863 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1790 xp_upper = 1
1791# 864 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1792 else
1793 ! if it is not periodic, then both elements are 0
1794 xp_lower = 0
1795# 867 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1796 xp_upper = 0
1797# 868 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1798 end if
1799# 860 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1800 ! check for periodicity
1801 if (ib_bc_y%beg == bc_periodic) then
1802# 862 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1803 yp_lower = -1
1804# 863 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1805 yp_upper = 1
1806# 864 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1807 else
1808 ! if it is not periodic, then both elements are 0
1809 yp_lower = 0
1810# 867 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1811 yp_upper = 0
1812# 868 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1813 end if
1814# 870 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1815
1816 ! z only if 3D
1817 if (present(zp_lower) .and. num_dims == 3) then
1818 if (ib_bc_z%beg == bc_periodic) then
1819 zp_lower = -1
1820 zp_upper = 1
1821 else
1822 zp_lower = 0
1823 zp_upper = 0
1824 end if
1825 else if (present(zp_lower)) then
1826 zp_lower = 0
1827 zp_upper = 0
1828 end if
1829
1830 end subroutine s_get_periodicities
1831
1832 !> Archimedes spiral function
1833 pure elemental function f_r(myth, offset, a)
1834
1835
1836# 890 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1837#if MFC_OpenACC
1838# 890 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1839!$acc routine seq
1840# 890 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1841#elif MFC_OpenMP
1842# 890 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1843
1844# 890 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1845
1846# 890 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1847!$omp declare target device_type(any)
1848# 890 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1849#endif
1850# 891 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1851 real(wp), intent(in) :: myth, offset, a
1852 real(wp) :: b
1853 real(wp) :: f_r
1854
1855 ! r(th) = a + b*th
1856
1857 b = 2._wp*a/(2._wp*pi)
1858 f_r = a + b*myth + offset
1859
1860 end function f_r
1861
1862end module m_ib_patches
k
integer, intent(in) k
Definition m_phase_change.fpp.f90:1043
j
integer, intent(in) j
Definition m_phase_change.fpp.f90:1043
m_constants
Compile-time constant parameters: default values, tolerances, and physical constants.
Definition m_constants.fpp.f90:7
m_derived_types
Shared derived types for field data, patch geometry, bubble dynamics, and MPI I/O structures.
Definition m_derived_types.fpp.f90:356
m_global_parameters
Global parameters for the computational domain, fluid properties, and simulation algorithm configurat...
Definition m_global_parameters.fpp.f90:367
m_global_parameters::glb_bounds
type(bounds_info), dimension(3) glb_bounds
Definition m_global_parameters.fpp.f90:428
m_global_parameters::y_cc
real(wp), dimension(:), allocatable, target y_cc
Definition m_global_parameters.fpp.f90:433
m_global_parameters::z_cc
real(wp), dimension(:), allocatable, target z_cc
Definition m_global_parameters.fpp.f90:433
m_global_parameters::num_procs
integer num_procs
Number of processors.
Definition m_global_parameters.fpp.f90:384
m_global_parameters::x_cc
real(wp), dimension(:), allocatable, target x_cc
Definition m_global_parameters.fpp.f90:433
m_global_parameters::ib_airfoil_grids
type(ib_airfoil_grid), dimension(num_ib_airfoils_max) ib_airfoil_grids
Per-airfoil computed surface grids.
Definition m_global_parameters.fpp.f90:915
m_helper_basic
Basic floating-point utilities: approximate equality, default detection, and coordinate bounds.
Definition m_helper_basic.fpp.f90:356
m_helper
Utility routines for bubble model setup, coordinate transforms, array sampling, and special functions...
Definition m_helper.fpp.f90:367
m_ib_patches
Allocate memory and read initial condition data for IC extrusion.
Definition m_ib_patches.fpp.f90:426
m_ib_patches::s_update_ib_rotation_matrix
subroutine, public s_update_ib_rotation_matrix(patch_id)
Compute a rotation matrix for converting to the rotating frame of the boundary.
Definition m_ib_patches.fpp.f90:1240
m_ib_patches::get_indices_from_bounds
subroutine get_indices_from_bounds(left_bound, right_bound, cell_centers, left_index, right_index)
Definition m_ib_patches.fpp.f90:1472
m_ib_patches::s_get_neighborhood_idx
subroutine, public s_get_neighborhood_idx(gbl_idx, neighborhood_idx, num_entries)
binary search to retrieve the local IB patch index using the global index
Definition m_ib_patches.fpp.f90:1608
m_ib_patches::s_get_ib_bound
subroutine, public s_get_ib_bound(patch, bound)
Definition m_ib_patches.fpp.f90:1301
m_ib_patches::ib_new_keys
integer, dimension(num_ib_patches_max_namelist) ib_new_keys
Holds each step's arrivals while they are sorted and merged in. Host only.
Definition m_ib_patches.fpp.f90:459
m_ib_patches::ib_new_vals
integer, dimension(num_ib_patches_max_namelist) ib_new_vals
Definition m_ib_patches.fpp.f90:459
m_ib_patches::s_merge_ib_lookup
subroutine, public s_merge_ib_lookup(num_ibs_pre)
Fold patch_ib(num_ibs_pre+1:num_ibs) into the map: sort just the arrivals, then merge the two sorted ...
Definition m_ib_patches.fpp.f90:1727
m_ib_patches::s_get_bounding_indices
subroutine s_get_bounding_indices(patch, center, il, ir, jl, jr, kl, kr)
Definition m_ib_patches.fpp.f90:1363
m_ib_patches::ib_lookup_vals
integer, dimension(num_ib_patches_max_namelist) ib_lookup_vals
Definition m_ib_patches.fpp.f90:444
m_ib_patches::s_apply_ib_patches_ib_parallelism
subroutine s_apply_ib_patches_ib_parallelism(ib_markers)
Definition m_ib_patches.fpp.f90:706
m_ib_patches::s_apply_ib_patches_grid_cell_parallelism
subroutine s_apply_ib_patches_grid_cell_parallelism(ib_markers)
Definition m_ib_patches.fpp.f90:477
m_ib_patches::s_apply_ib_patches
impure subroutine, public s_apply_ib_patches(ib_markers)
Apply all immersed boundary patch geometries to mark interior cells in the IB marker array.
Definition m_ib_patches.fpp.f90:465
m_ib_patches::ib_lookup_keys
integer, dimension(num_ib_patches_max_namelist) ib_lookup_keys
lookup arrays for converting global IB indices to local indices
Definition m_ib_patches.fpp.f90:444
m_ib_patches::s_update_ib_lookup
subroutine, public s_update_ib_lookup()
Completely rebuilds the ib lookup map, used at startup.
Definition m_ib_patches.fpp.f90:1652
m_ib_patches::s_get_periodicities
subroutine, public s_get_periodicities(xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper)
Determine the periodic wrapping bounds in each direction.
Definition m_ib_patches.fpp.f90:1778
m_ib_patches::s_decode_patch_periodicity
subroutine, public s_decode_patch_periodicity(encoded_patch_id, patch_id, x_periodicity, y_periodicity, z_periodicity)
Decode the encoded ID to recover the original patch ID and periodicity.
Definition m_ib_patches.fpp.f90:1564
m_ib_patches::s_initialize_ib_airfoils
subroutine, public s_initialize_ib_airfoils()
Initialize the NACA surface grids for all airfoil IB patches. Must be called after the grid is establ...
Definition m_ib_patches.fpp.f90:1058
m_ib_patches::s_encode_patch_periodicity
subroutine, public s_encode_patch_periodicity(patch_id, x_periodicity, y_periodicity, z_periodicity, encoded_patch_id)
Encode the patch ID with a unique offset containing periodicity information.
Definition m_ib_patches.fpp.f90:1529
m_ib_patches::s_compact_ib_lookup
subroutine, public s_compact_ib_lookup(old_to_new, num_ibs_old)
Drop the entries whose patches left the neighborhood and renumber the survivors onto the patch_ib slo...
Definition m_ib_patches.fpp.f90:1687
m_ib_patches::f_r
pure elemental real(wp) function f_r(myth, offset, a)
Archimedes spiral function.
Definition m_ib_patches.fpp.f90:1834
m_model
Binary STL file reader and processor for immersed boundary geometry.
Definition m_model.fpp.f90:357
m_model::gpu_ntrs
integer, dimension(:), allocatable, public gpu_ntrs
GPU-friendly flat arrays for STL model data.
Definition m_model.fpp.f90:378
m_model::f_model_is_inside
real(wp) function, public f_model_is_inside(ntrs, pid, point)
Determine if a point is inside a surface using the generalized winding number (Jacobson et al....
Definition m_model.fpp.f90:826
m_model::s_instantiate_stl_models
subroutine, public s_instantiate_stl_models()
Load, transform, and register STL/OBJ immersed-boundary models onto the simulation grid.
Definition m_model.fpp.f90:1303
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_allreduce_min
impure subroutine s_mpi_allreduce_min(var_loc, var_glb)
Reduce a local real value to its global minimum across all MPI ranks.
Definition m_mpi_common.fpp.f90:927
m_patch_geometries
Contains helper functions specific to various patch gemoetries for determining if a grid cell lies in...
Definition m_patch_geometries.fpp.f90:357
m_patch_geometries::f_is_inside_sphere
logical function, public f_is_inside_sphere(x, y, z, radius)
Check if the x, y, and z coordinates would be located inside a sphere with the patch_id's radius.
Definition m_patch_geometries.fpp.f90:375
m_patch_geometries::f_is_inside_cuboid
logical function, public f_is_inside_cuboid(x, y, z, length)
Check if the x, y, and possibly z coordinates would be located inside a cuboid with the patch_id's le...
Definition m_patch_geometries.fpp.f90:434
m_patch_geometries::f_is_inside_cylinder
logical function, public f_is_inside_cylinder(polar_x, polar_y, height, radius, length)
Check which length of the cylinder is not default. Use that direction as the height and the other two...
Definition m_patch_geometries.fpp.f90:403
m_patch_geometries::f_is_inside_ellipse
logical function, public f_is_inside_ellipse(x, y, length)
Definition m_patch_geometries.fpp.f90:524
m_patch_geometries::f_is_inside_airfoil
logical function, public f_is_inside_airfoil(x, y, z, length, airfoil_id)
Check if the x, y, are bounded by a NACA airfoil. Check if the z coordinate is inside the left and ri...
Definition m_patch_geometries.fpp.f90:467
m_derived_types::integer_field
Derived type annexing an integer scalar field (SF).
Definition m_derived_types.fpp.f90:380