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# 15 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
90# 16 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
91# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
92
93# 24 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
94
95# 53 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
96
97# 65 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
98
99# 75 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
100
101# 105 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
102
103# 117 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
104
105# 127 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
106
107# 174 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
108! New line at end of file is required for FYPP
109# 2 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
110# 1 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 1
111# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
112# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
113# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
114# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
115# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
116# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
117
118# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
119# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
120# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
121
122# 15 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
123# 16 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
124# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
125
126# 24 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
127
128# 53 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
129
130# 65 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
131
132# 75 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
133
134# 105 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
135
136# 117 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
137
138# 127 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
139
140# 174 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
141! New line at end of file is required for FYPP
142# 2 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 2
143
144# 4 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
145# 5 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
146# 6 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
147# 7 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
148# 8 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
149
150# 20 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
151
152# 43 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
153
154# 48 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
155
156# 53 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
157
158# 58 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
159
160# 63 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
161
162# 68 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
163
164# 81 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
165
166# 86 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
167
168# 91 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
169
170# 96 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
171
172# 101 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
173
174# 106 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
175
176# 111 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
177
178# 116 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
179
180# 121 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
181
182# 126 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
183
184# 156 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
185
186# 197 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
187
188# 211 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
189
190# 236 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
191
192# 247 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
193
194# 249 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
195# 260 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
196
197# 310 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
198
199# 320 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
200
201# 330 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
202
203# 339 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
204
205# 356 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
206
207# 366 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
208
209# 373 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
210
211# 379 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
212
213# 385 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
214
215# 391 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
216
217# 397 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
218
219# 403 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
220! New line at end of file is required for FYPP
221# 3 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
222# 1 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 1
223# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
224# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
225# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
226# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
227# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
228# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
229
230# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
231# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
232# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
233
234# 15 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
235# 16 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
236# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
237
238# 24 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
239
240# 53 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
241
242# 65 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
243
244# 75 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
245
246# 105 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
247
248# 117 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
249
250# 127 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
251
252# 174 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
253! New line at end of file is required for FYPP
254# 2 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 2
255
256# 7 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
257
258# 17 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
259
260# 22 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
261
262# 27 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
263
264# 32 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
265
266# 37 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
267
268# 42 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
269
270# 47 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
271
272# 52 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
273
274# 57 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
275
276# 62 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
277
278# 73 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
279
280# 78 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
281
282# 83 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
283
284# 88 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
285
286# 103 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
287
288# 131 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
289
290# 160 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
291
292# 175 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
293
294# 193 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
295
296# 215 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
297
298# 244 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
299
300# 259 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
301
302# 269 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
303
304# 278 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
305
306# 294 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
307
308# 304 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
309
310# 311 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
311! New line at end of file is required for FYPP
312# 4 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
313
314! GPU parallel region (scalar reductions, maxval/minval)
315# 23 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
316
317! GPU parallel loop over threads (most common GPU macro)
318# 43 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
319
320! Required closing for GPU_PARALLEL_LOOP
321# 55 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
322
323! Mark routine for device compilation
324# 112 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
325
326! Declare device-resident data
327# 130 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
328
329! Inner loop within a GPU parallel region
330# 145 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
331
332! Scoped GPU data region
333# 164 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
334
335! Host code with device pointers (for MPI with GPU buffers)
336# 193 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
337
338! Allocate device memory (unscoped)
339# 207 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
340
341! Free device memory
342# 219 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
343
344! Atomic operation on device
345# 231 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
346
347! End atomic capture block
348# 242 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
349
350! Copy data between host and device
351# 254 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
352
353! Synchronization barrier
354# 266 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
355
356! Import GPU library module (openacc or omp_lib)
357# 275 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
358
359! Emit code only for AMD compiler
360# 282 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
361
362! Emit code for non-Cray compilers
363# 289 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
364
365! Emit code only for Cray compiler
366# 296 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
367
368! Emit code for non-NVIDIA compilers
369# 303 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
370
371# 305 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
372# 306 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
373! New line at end of file is required for FYPP
374# 2 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 2
375
376# 14 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
377
378! Caution: This macro requires the use of a binding script to set CUDA_VISIBLE_DEVICES, such that we have one GPU device per MPI
379! rank. That's because for both cudaMemAdvise (preferred location) and cudaMemPrefetchAsync we use location = device_id = 0. For an
380! example see misc/nvidia_uvm/bind.sh.
381# 52 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
382
383! Allocate and create GPU device memory
384# 72 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
385
386! Free GPU device memory and deallocate
387# 80 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
388
389! Cray-specific GPU pointer setup for vector fields
390# 104 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
391
392! Cray-specific GPU pointer setup for scalar fields
393# 120 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
394
395! Cray-specific GPU pointer setup for acoustic source spatials
396# 145 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
397
398# 151 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
399
400# 158 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
401! New line at end of file is required for FYPP
402# 11 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp" 2
403
404!> @brief Immersed boundary patch geometry constructors for 2D and 3D shapes
405module m_ib_patches
406
407 use m_patch_geometries
408 use m_model ! Subroutine(s) related to STL files
409 use m_derived_types ! Definitions of the derived types
410 use m_global_parameters
411 use m_helper_basic
412 use m_helper
413 use m_mpi_common
414
415 implicit none
416
417 private; public :: s_apply_ib_patches, s_update_ib_rotation_matrix, s_instantiate_stl_models, s_decode_patch_periodicity, &
418 & s_encode_patch_periodicity, s_initialize_ib_airfoils, s_get_periodicities, s_get_ib_bound
419
420contains
421
422 !> Apply all immersed boundary patch geometries to mark interior cells in the IB marker array
423 impure subroutine s_apply_ib_patches(ib_markers)
424
425 type(integer_field), intent(inout) :: ib_markers
426
427 if (many_ib_patch_parallelism) then
428 call s_apply_ib_patches_ib_parallelism(ib_markers)
429 else
430 call s_apply_ib_patches_grid_cell_parallelism(ib_markers)
431 end if
432
433 end subroutine s_apply_ib_patches
434
435 subroutine s_apply_ib_patches_grid_cell_parallelism(ib_markers)
436
437 type(integer_field), intent(inout) :: ib_markers
438 integer :: patch_id, airfoil_id, model_id, encoded_patch_id, i, j, k, il, ir, jl, jr, kl, kr, xp, yp, zp !< iterators
439 integer :: xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper !< periodic bounds
440 real(wp), dimension(3) :: center, xyz_local, length
441 real(wp) :: radius, eta
442
443
444# 51 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
445#if defined(MFC_OpenACC)
446# 51 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
447!$acc update host(patch_ib(1:num_ibs))
448# 51 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
449#elif defined(MFC_OpenMP)
450# 51 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
451!$omp target update from(patch_ib(1:num_ibs))
452# 51 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
453#endif
454
455 ! 3D Patch Geometries
456 if (num_dims == 3) then
457 call s_get_periodicities(xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper)
458 do xp = xp_lower, xp_upper
459 do yp = yp_lower, yp_upper
460 do zp = zp_lower, zp_upper
461 do patch_id = 1, num_ibs
462 center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(glb_bounds(1)%end - glb_bounds(1)%beg)
463 center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(glb_bounds(2)%end - glb_bounds(2)%beg)
464 center(3) = patch_ib(patch_id)%z_centroid + real(zp, wp)*(glb_bounds(3)%end - glb_bounds(3)%beg)
465
466 ! encode the periodicity information into the patch_id
467 call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, zp, encoded_patch_id)
468
469 ! find the indices to the left and right of the IB in i, j, k
470 call s_get_bounding_indices(patch_ib(patch_id), center, il, ir, jl, jr, kl, kr)
471
472 ! skip patches whose bounding box does not overlap this rank's domain
473 if (ir < il .or. jr < jl .or. kr < kl) cycle
474
475
476# 73 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
477
478# 73 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
479#if defined(MFC_OpenACC)
480# 73 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
481!$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)
482# 73 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
483#elif defined(MFC_OpenMP)
484# 73 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
485
486# 73 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
487
488# 73 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
489
490# 73 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
491!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) &
492# 73 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
493!$omp& private(i, j, k, xyz_local, length, radius, airfoil_id, eta) map(to:patch_id, encoded_patch_id, center)
494# 73 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
495#endif
496# 75 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
497 do k = kl, kr
498 do j = jl, jr
499 do i = il, ir
500 ! get coordinate frame centered on IB
501 xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), z_cc(k) - center(3)]
502 ! rotate the frame into the IB's coordinates
503 xyz_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse, xyz_local)
504
505 ! perform the interior check for the patch geometry of this IB
506 if (patch_ib(patch_id)%geometry == 8) then
507 ! sphere geometry
508 radius = patch_ib(patch_id)%radius
509
510 if (f_is_inside_sphere(xyz_local(1), xyz_local(2), xyz_local(3), &
511 & radius)) ib_markers%sf(i, j, k) = encoded_patch_id
512 else if (patch_ib(patch_id)%geometry == 9) then
513 ! cuboid geometry
514 length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, &
515 & patch_ib(patch_id)%length_z]
516 if (f_is_inside_cuboid(xyz_local(1), xyz_local(2), xyz_local(3), &
517 & length)) ib_markers%sf(i, j, k) = encoded_patch_id
518 else if (patch_ib(patch_id)%geometry == 10) then
519 ! cylinder geometry
520 radius = patch_ib(patch_id)%radius
521 if (f_is_inside_cylinder(xyz_local(2), xyz_local(3), xyz_local(1), radius, &
522 & patch_ib(patch_id)%length_x)) ib_markers%sf(i, j, k) = encoded_patch_id
523 else if (patch_ib(patch_id)%geometry == 11) then
524 ! 3D airfoil geometry
525 airfoil_id = patch_ib(patch_id)%airfoil_id
526 xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset
527 if (f_is_inside_airfoil(xyz_local(1), xyz_local(2), xyz_local(3), &
528 & patch_ib(patch_id)%length_z, airfoil_id)) ib_markers%sf(i, j, &
529 & k) = encoded_patch_id
530 else if (patch_ib(patch_id)%geometry == 12) then
531 ! STL model geometry
532 xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset
533 model_id = patch_ib(patch_id)%model_id
534 eta = f_model_is_inside(gpu_ntrs(model_id), model_id, xyz_local)
535 if (eta > stl_models(model_id)%model_threshold) then
536 ib_markers%sf(i, j, k) = encoded_patch_id
537 end if
538 end if
539 end do
540 end do
541 end do
542
543# 120 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
544#if defined(MFC_OpenACC)
545# 120 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
546!$acc end parallel loop
547# 120 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
548#elif defined(MFC_OpenMP)
549# 120 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
550
551# 120 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
552!$omp end target teams loop
553# 120 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
554#endif
555 end do
556 end do
557 end do
558 end do
559
560 ! 2D Patch Geometries
561 else if (num_dims == 2) then
562 call s_get_periodicities(xp_lower, xp_upper, yp_lower, yp_upper)
563 do xp = xp_lower, xp_upper
564 do yp = yp_lower, yp_upper
565 do patch_id = 1, num_ibs
566 center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(glb_bounds(1)%end - glb_bounds(1)%beg)
567 center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(glb_bounds(2)%end - glb_bounds(2)%beg)
568 center(3) = 0._wp
569
570 ! encode the periodicity information into the patch_id
571 call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, 0, encoded_patch_id)
572
573 ! find the indices to the left and right of the IB in i, j, k
574 call s_get_bounding_indices(patch_ib(patch_id), center, il, ir, jl, jr, kl, kr)
575
576 ! skip patches whose bounding box does not overlap this rank's domain
577 if (ir < il .or. jr < jl) cycle
578
579
580# 145 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
581
582# 145 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
583#if defined(MFC_OpenACC)
584# 145 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
585!$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)
586# 145 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
587#elif defined(MFC_OpenMP)
588# 145 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
589
590# 145 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
591
592# 145 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
593
594# 145 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
595!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(2) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) &
596# 145 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
597!$omp& private(i, j, xyz_local, airfoil_id, eta, length, radius) map(to:patch_id, encoded_patch_id, center)
598# 145 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
599#endif
600# 147 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
601 do j = jl, jr
602 do i = il, ir
603 ! get coordinate frame centered on IB
604 xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp]
605 ! rotate the frame into the IB's coordinates
606 xyz_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse, xyz_local)
607
608 ! perform the interior check for the patch geometry of this IB
609 if (patch_ib(patch_id)%geometry == 2) then
610 ! circular geometries
611 radius = patch_ib(patch_id)%radius
612 if (f_is_inside_cylinder(xyz_local(1), xyz_local(2), 0._wp, radius, 0._wp)) ib_markers%sf(i, &
613 & j, 0) = encoded_patch_id
614 else if (patch_ib(patch_id)%geometry == 3) then
615 ! rectangular geometries
616 length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, 0._wp]
617 if (f_is_inside_cuboid(xyz_local(1), xyz_local(2), xyz_local(3), length)) ib_markers%sf(i, j, &
618 & 0) = encoded_patch_id
619 else if (patch_ib(patch_id)%geometry == 4) then
620 ! 2D airfoil geometry
621 airfoil_id = patch_ib(patch_id)%airfoil_id
622 xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset
623 if (f_is_inside_airfoil(xyz_local(1), xyz_local(2), 0._wp, 0._wp, &
624 & airfoil_id)) ib_markers%sf(i, j, 0) = encoded_patch_id
625 else if (patch_ib(patch_id)%geometry == 5) then
626 ! STL model geometry
627 xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset
628 model_id = patch_ib(patch_id)%model_id
629 eta = f_model_is_inside(gpu_ntrs(model_id), model_id, xyz_local)
630 if (eta > stl_models(model_id)%model_threshold) then
631 ib_markers%sf(i, j, 0) = encoded_patch_id
632 end if
633 else if (patch_ib(patch_id)%geometry == 6) then
634 ! ellipse geometry
635 length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, 0._wp]
636 if (f_is_inside_ellipse(xyz_local(1), xyz_local(2), length)) ib_markers%sf(i, j, &
637 & 0) = encoded_patch_id
638 end if
639 end do
640 end do
641
642# 187 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
643#if defined(MFC_OpenACC)
644# 187 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
645!$acc end parallel loop
646# 187 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
647#elif defined(MFC_OpenMP)
648# 187 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
649
650# 187 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
651!$omp end target teams loop
652# 187 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
653#endif
654 end do
655 end do
656 end do
657 end if
658
659 end subroutine s_apply_ib_patches_grid_cell_parallelism
660
661 subroutine s_apply_ib_patches_ib_parallelism(ib_markers)
662
663 type(integer_field), intent(inout) :: ib_markers
664 integer :: patch_id, airfoil_id, model_id, encoded_patch_id, i, j, k, il, ir, jl, jr, kl, kr, xp, yp, zp !< iterators
665 integer :: xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper !< periodic bounds
666 real(wp), dimension(3) :: center, xyz_local, length
667 real(wp) :: radius, eta
668
669 if (num_dims == 3) then
670 ! get the periodicities
671 call s_get_periodicities(xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper)
672
673 do xp = xp_lower, xp_upper
674 do yp = yp_lower, yp_upper
675 do zp = zp_lower, zp_upper
676
677# 210 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
678
679# 210 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
680#if defined(MFC_OpenACC)
681# 210 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
682!$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)
683# 210 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
684#elif defined(MFC_OpenMP)
685# 210 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
686
687# 210 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
688
689# 210 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
690
691# 210 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
692!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) &
693# 210 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
694!$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)
695# 210 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
696#endif
697# 212 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
698 do patch_id = 1, num_ibs
699 center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(glb_bounds(1)%end - glb_bounds(1)%beg)
700 center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(glb_bounds(2)%end - glb_bounds(2)%beg)
701 center(3) = patch_ib(patch_id)%z_centroid + real(zp, wp)*(glb_bounds(3)%end - glb_bounds(3)%beg)
702
703 ! encode the periodicity information into the patch_id
704 call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, zp, encoded_patch_id)
705
706 ! find the indices to the left and right of the IB in i, j, k
707 call s_get_bounding_indices(patch_ib(patch_id), center, il, ir, jl, jr, kl, kr)
708
709 do k = kl, kr
710 do j = jl, jr
711 do i = il, ir
712 ! get coordinate frame centered on IB
713 xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), z_cc(k) - center(3)]
714 ! rotate the frame into the IB's coordinates
715 xyz_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse, xyz_local)
716
717 ! perform the interior check for the patch geometry of this IB
718 if (patch_ib(patch_id)%geometry == 8) then
719 ! sphere geometry
720 radius = patch_ib(patch_id)%radius
721
722 if (f_is_inside_sphere(xyz_local(1), xyz_local(2), xyz_local(3), &
723 & radius)) ib_markers%sf(i, j, k) = encoded_patch_id
724 else if (patch_ib(patch_id)%geometry == 9) then
725 ! cuboid geometry
726 length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, &
727 & patch_ib(patch_id)%length_z]
728 if (f_is_inside_cuboid(xyz_local(1), xyz_local(2), xyz_local(3), &
729 & length)) ib_markers%sf(i, j, k) = encoded_patch_id
730 else if (patch_ib(patch_id)%geometry == 10) then
731 ! cylinder geometry
732 radius = patch_ib(patch_id)%radius
733 if (f_is_inside_cylinder(xyz_local(2), xyz_local(3), xyz_local(1), radius, &
734 & patch_ib(patch_id)%length_x)) ib_markers%sf(i, j, k) = encoded_patch_id
735 else if (patch_ib(patch_id)%geometry == 11) then
736 ! 3D airfoil geometry
737 airfoil_id = patch_ib(patch_id)%airfoil_id
738 xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset
739 if (f_is_inside_airfoil(xyz_local(1), xyz_local(2), xyz_local(3), &
740 & patch_ib(patch_id)%length_z, airfoil_id)) ib_markers%sf(i, j, &
741 & k) = encoded_patch_id
742 else if (patch_ib(patch_id)%geometry == 12) then
743 ! STL model geometry
744 xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset
745 model_id = patch_ib(patch_id)%model_id
746 eta = f_model_is_inside(gpu_ntrs(model_id), model_id, xyz_local)
747 if (eta > stl_models(model_id)%model_threshold) then
748 ib_markers%sf(i, j, k) = encoded_patch_id
749 end if
750 end if
751 end do
752 end do
753 end do
754 end do
755
756# 269 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
757#if defined(MFC_OpenACC)
758# 269 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
759!$acc end parallel loop
760# 269 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
761#elif defined(MFC_OpenMP)
762# 269 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
763
764# 269 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
765!$omp end target teams loop
766# 269 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
767#endif
768 end do
769 end do
770 end do
771 else if (num_dims == 2) then
772 ! get the periodicities
773 call s_get_periodicities(xp_lower, xp_upper, yp_lower, yp_upper)
774
775 do xp = xp_lower, xp_upper
776 do yp = yp_lower, yp_upper
777
778# 279 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
779
780# 279 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
781#if defined(MFC_OpenACC)
782# 279 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
783!$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)
784# 279 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
785#elif defined(MFC_OpenMP)
786# 279 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
787
788# 279 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
789
790# 279 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
791
792# 279 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
793!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) &
794# 279 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
795!$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)
796# 279 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
797#endif
798# 281 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
799 do patch_id = 1, num_ibs
800 center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(glb_bounds(1)%end - glb_bounds(1)%beg)
801 center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(glb_bounds(2)%end - glb_bounds(2)%beg)
802 center(3) = 0._wp
803
804 ! encode the periodicity information into the patch_id
805 call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, 0, encoded_patch_id)
806
807 ! find the indices to the left and right of the IB in i, j, k
808 call s_get_bounding_indices(patch_ib(patch_id), center, il, ir, jl, jr, kl, kr)
809
810 do j = jl, jr
811 do i = il, ir
812 ! get coordinate frame centered on IB
813 xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp]
814 ! rotate the frame into the IB's coordinates
815 xyz_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse, xyz_local)
816
817 ! perform the interior check for the patch geometry of this IB
818 if (patch_ib(patch_id)%geometry == 2) then
819 ! circular geometries
820 radius = patch_ib(patch_id)%radius
821 if (f_is_inside_cylinder(xyz_local(1), xyz_local(2), 0._wp, radius, 0._wp)) ib_markers%sf(i, &
822 & j, 0) = encoded_patch_id
823 else if (patch_ib(patch_id)%geometry == 3) then
824 ! rectangular geometries
825 length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, 0._wp]
826 if (f_is_inside_cuboid(xyz_local(1), xyz_local(2), xyz_local(3), length)) ib_markers%sf(i, j, &
827 & 0) = encoded_patch_id
828 else if (patch_ib(patch_id)%geometry == 4) then
829 ! 2D airfoil geometry
830 airfoil_id = patch_ib(patch_id)%airfoil_id
831 xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset
832 if (f_is_inside_airfoil(xyz_local(1), xyz_local(2), 0._wp, 0._wp, &
833 & airfoil_id)) ib_markers%sf(i, j, 0) = encoded_patch_id
834 else if (patch_ib(patch_id)%geometry == 5) then
835 ! STL model geometry
836 xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset
837 model_id = patch_ib(patch_id)%model_id
838 eta = f_model_is_inside(gpu_ntrs(model_id), model_id, xyz_local)
839 if (eta > stl_models(model_id)%model_threshold) then
840 ib_markers%sf(i, j, 0) = encoded_patch_id
841 end if
842 else if (patch_ib(patch_id)%geometry == 6) then
843 ! ellipse geometry
844 length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, 0._wp]
845 if (f_is_inside_ellipse(xyz_local(1), xyz_local(2), length)) ib_markers%sf(i, j, &
846 & 0) = encoded_patch_id
847 end if
848 end do
849 end do
850 end do
851
852# 333 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
853#if defined(MFC_OpenACC)
854# 333 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
855!$acc end parallel loop
856# 333 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
857#elif defined(MFC_OpenMP)
858# 333 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
859
860# 333 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
861!$omp end target teams loop
862# 333 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
863#endif
864 end do
865 end do
866 end if
867
868 end subroutine s_apply_ib_patches_ib_parallelism
869
870 !> Initialize the NACA surface grids for all airfoil IB patches. Must be called after the grid is established (so dx is valid)
871 !! and before s_apply_ib_patches or s_apply_levelset.
872 subroutine s_initialize_ib_airfoils()
873
874 integer :: i, j, airfoil_id
875 integer :: np, np1, np2
876 real(wp) :: ca_in, pa, ma, ta
877 real(wp) :: xc, xa, yc, dycdxc, yt, xu, yu, xl, yl, sin_c, cos_c
878
879 do i = 1, num_ibs
880 if (patch_ib(i)%geometry /= 4 .and. patch_ib(i)%geometry /= 11) cycle
881
882 airfoil_id = patch_ib(i)%airfoil_id
883 ca_in = ib_airfoil(airfoil_id)%c
884 pa = ib_airfoil(airfoil_id)%p
885 ma = ib_airfoil(airfoil_id)%m
886 ta = ib_airfoil(airfoil_id)%t
887
888 np1 = int((pa*ca_in/dx(0))*20)
889 np2 = int(((ca_in - pa*ca_in)/dx(0))*20)
890 np = np1 + np2 + 1
891 ib_airfoil_grids(airfoil_id)%Np = np
892
893# 362 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
894#if defined(MFC_OpenACC)
895# 362 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
896!$acc update device(ib_airfoil_grids(airfoil_id)%Np)
897# 362 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
898#elif defined(MFC_OpenMP)
899# 362 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
900!$omp target update to(ib_airfoil_grids(airfoil_id)%Np)
901# 362 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
902#endif
903
904 if (.not. allocated(ib_airfoil_grids(airfoil_id)%upper)) then
905#ifdef MFC_DEBUG
906# 365 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
907 block
908# 365 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
909 use iso_fortran_env, only: output_unit
910# 365 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
911
912# 365 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
913 print *, 'm_ib_patches.fpp:365: ', '@:ALLOCATE(ib_airfoil_grids(airfoil_id)%upper(1:Np))'
914# 365 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
915
916# 365 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
917 call flush (output_unit)
918# 365 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
919 end block
920# 365 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
921#endif
922# 365 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
923 allocate (ib_airfoil_grids(airfoil_id)%upper(1:np))
924# 365 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
925
926# 365 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
927
928# 365 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
929#if defined(MFC_OpenACC)
930# 365 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
931!$acc enter data create(ib_airfoil_grids(airfoil_id)%upper)
932# 365 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
933#elif defined(MFC_OpenMP)
934# 365 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
935!$omp target enter data map(always,alloc:ib_airfoil_grids(airfoil_id)%upper)
936# 365 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
937#endif
938#ifdef MFC_DEBUG
939# 366 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
940 block
941# 366 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
942 use iso_fortran_env, only: output_unit
943# 366 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
944
945# 366 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
946 print *, 'm_ib_patches.fpp:366: ', '@:ALLOCATE(ib_airfoil_grids(airfoil_id)%lower(1:Np))'
947# 366 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
948
949# 366 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
950 call flush (output_unit)
951# 366 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
952 end block
953# 366 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
954#endif
955# 366 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
956 allocate (ib_airfoil_grids(airfoil_id)%lower(1:np))
957# 366 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
958
959# 366 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
960
961# 366 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
962#if defined(MFC_OpenACC)
963# 366 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
964!$acc enter data create(ib_airfoil_grids(airfoil_id)%lower)
965# 366 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
966#elif defined(MFC_OpenMP)
967# 366 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
968!$omp target enter data map(always,alloc:ib_airfoil_grids(airfoil_id)%lower)
969# 366 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
970#endif
971
972 ib_airfoil_grids(airfoil_id)%upper(1)%x = 0._wp
973 ib_airfoil_grids(airfoil_id)%upper(1)%y = 0._wp
974 ib_airfoil_grids(airfoil_id)%lower(1)%x = 0._wp
975 ib_airfoil_grids(airfoil_id)%lower(1)%y = 0._wp
976
977 do j = 1, np1 + np2 - 1
978 if (j <= np1) then
979 xc = j*(pa*ca_in/np1)
980 xa = xc/ca_in
981 yc = (ma/pa**2)*(2*pa*xa - xa**2)
982 dycdxc = (2*ma/pa**2)*(pa - xa)
983 else
984 xc = pa*ca_in + (j - np1)*((ca_in - pa*ca_in)/np2)
985 xa = xc/ca_in
986 yc = (ma/(1 - pa)**2)*(1 - 2*pa + 2*pa*xa - xa**2)
987 dycdxc = (2*ma/(1 - pa)**2)*(pa - xa)
988 end if
989
990 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)
991 sin_c = dycdxc/(1 + dycdxc**2)**0.5_wp
992 cos_c = 1/(1 + dycdxc**2)**0.5_wp
993
994 xu = (xa - yt*sin_c)*ca_in
995 yu = (yc + yt*cos_c)*ca_in
996 xl = (xa + yt*sin_c)*ca_in
997 yl = (yc - yt*cos_c)*ca_in
998
999 ib_airfoil_grids(airfoil_id)%upper(j + 1)%x = xu
1000 ib_airfoil_grids(airfoil_id)%upper(j + 1)%y = yu
1001 ib_airfoil_grids(airfoil_id)%lower(j + 1)%x = xl
1002 ib_airfoil_grids(airfoil_id)%lower(j + 1)%y = yl
1003 end do
1004
1005 ib_airfoil_grids(airfoil_id)%upper(np)%x = ca_in
1006 ib_airfoil_grids(airfoil_id)%upper(np)%y = 0._wp
1007 ib_airfoil_grids(airfoil_id)%lower(np)%x = ca_in
1008 ib_airfoil_grids(airfoil_id)%lower(np)%y = 0._wp
1009
1010
1011# 406 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1012#if defined(MFC_OpenACC)
1013# 406 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1014!$acc update device(ib_airfoil_grids(airfoil_id)%upper, ib_airfoil_grids(airfoil_id)%lower)
1015# 406 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1016#elif defined(MFC_OpenMP)
1017# 406 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1018!$omp target update to(ib_airfoil_grids(airfoil_id)%upper, ib_airfoil_grids(airfoil_id)%lower)
1019# 406 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1020#endif
1021 end if
1022
1023
1024# 409 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1025#if defined(MFC_OpenACC)
1026# 409 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1027!$acc update device(ib_airfoil(airfoil_id))
1028# 409 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1029#elif defined(MFC_OpenMP)
1030# 409 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1031!$omp target update to(ib_airfoil(airfoil_id))
1032# 409 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1033#endif
1034 end do
1035
1036 end subroutine s_initialize_ib_airfoils
1037
1038 !> Compute a rotation matrix for converting to the rotating frame of the boundary
1039 subroutine s_update_ib_rotation_matrix(patch_id)
1040
1041
1042# 417 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1043#if MFC_OpenACC
1044# 417 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1045!$acc routine seq
1046# 417 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1047#elif MFC_OpenMP
1048# 417 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1049
1050# 417 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1051
1052# 417 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1053!$omp declare target device_type(any)
1054# 417 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1055#endif
1056
1057 integer, intent(in) :: patch_id
1058 real(wp), dimension(3, 3, 3) :: rotation
1059 real(wp) :: angle
1060
1061 ! construct the x, y, and z rotation matrices
1062
1063 if (num_dims == 3) then
1064 ! also compute the x and y axes in 3D
1065 angle = patch_ib(patch_id)%angles(1)
1066 rotation(1, 1,:) = [1._wp, 0._wp, 0._wp]
1067 rotation(1, 2,:) = [0._wp, cos(angle), -sin(angle)]
1068 rotation(1, 3,:) = [0._wp, sin(angle), cos(angle)]
1069
1070 angle = patch_ib(patch_id)%angles(2)
1071 rotation(2, 1,:) = [cos(angle), 0._wp, sin(angle)]
1072 rotation(2, 2,:) = [0._wp, 1._wp, 0._wp]
1073 rotation(2, 3,:) = [-sin(angle), 0._wp, cos(angle)]
1074
1075 ! apply the y rotation to the x rotation
1076 patch_ib(patch_id)%rotation_matrix(:,:) = matmul(rotation(1,:,:), rotation(2,:,:))
1077 patch_ib(patch_id)%rotation_matrix_inverse(:,:) = matmul(transpose(rotation(2,:,:)), transpose(rotation(1,:,:)))
1078 end if
1079
1080 ! z component first, since it applies in 2D and 3D
1081 angle = patch_ib(patch_id)%angles(3)
1082 rotation(3, 1,:) = [cos(angle), -sin(angle), 0._wp]
1083 rotation(3, 2,:) = [sin(angle), cos(angle), 0._wp]
1084 rotation(3, 3,:) = [0._wp, 0._wp, 1._wp]
1085
1086 if (num_dims == 3) then
1087 ! apply the z rotation to the xy rotation in 3D
1088 patch_ib(patch_id)%rotation_matrix(:,:) = matmul(patch_ib(patch_id)%rotation_matrix(:,:), rotation(3,:,:))
1089 patch_ib(patch_id)%rotation_matrix_inverse(:,:) = matmul(transpose(rotation(3,:,:)), &
1090 & patch_ib(patch_id)%rotation_matrix_inverse(:,:))
1091 else
1092 ! write out only the z rotation in 2D
1093 patch_ib(patch_id)%rotation_matrix(:,:) = rotation(3,:,:)
1094 patch_ib(patch_id)%rotation_matrix_inverse(:,:) = transpose(rotation(3,:,:))
1095 end if
1096
1097 end subroutine s_update_ib_rotation_matrix
1098
1099 subroutine s_get_ib_bound(patch, bound)
1100
1101
1102# 463 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1103#if MFC_OpenACC
1104# 463 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1105!$acc routine seq
1106# 463 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1107#elif MFC_OpenMP
1108# 463 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1109
1110# 463 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1111
1112# 463 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1113!$omp declare target device_type(any)
1114# 463 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1115#endif
1116
1117 type(ib_patch_parameters), intent(in) :: patch
1118 real(wp), intent(out) :: bound
1119 real(wp), dimension(2) :: lx, ly, lz
1120
1121 if (patch%geometry == 2 .or. patch%geometry == 8) then
1122 ! circle and sphere geometries
1123 bound = patch%radius
1124 else if (patch%geometry == 3) then
1125 bound = 0.5_wp*sqrt(patch%length_x**2 + patch%length_y**2)
1126 else if (patch%geometry == 4 .or. patch%geometry == 11) then
1127 ! rectangular geometries
1128 bound = ib_airfoil(patch%airfoil_id)%c
1129 else if (patch%geometry == 5) then
1130 ! STL model geometry
1131 lx(1) = stl_bounding_boxes(patch%model_id, 1, 1)
1132 lx(2) = stl_bounding_boxes(patch%model_id, 1, 3)
1133 ly(1) = stl_bounding_boxes(patch%model_id, 2, 1)
1134 ly(2) = stl_bounding_boxes(patch%model_id, 2, 3)
1135
1136 bound = 0.5_wp*sqrt((lx(2) - lx(1))**2 + (ly(2) - ly(1))**2)
1137 else if (patch%geometry == 6) then
1138 ! ellipse geometry
1139 bound = 0.5_wp*max(patch%length_x, patch%length_y)
1140 else if (patch%geometry == 9) then
1141 ! cuboid geometries
1142 bound = 0.5_wp*sqrt(patch%length_x**2 + patch%length_y**2 + patch%length_z**2)
1143 else if (patch%geometry == 10) then
1144 ! cylinder geometry
1145 bound = sqrt(patch%radius**2 + patch%length_x**2)
1146 else if (patch%geometry == 12) then
1147 ! Local-space bounding box extents (min=1, max=2 in the third index)
1148 lx(1) = stl_bounding_boxes(patch%model_id, 1, 1) + patch%centroid_offset(1)
1149 lx(2) = stl_bounding_boxes(patch%model_id, 1, 3) + patch%centroid_offset(1)
1150 ly(1) = stl_bounding_boxes(patch%model_id, 2, 1) + patch%centroid_offset(2)
1151 ly(2) = stl_bounding_boxes(patch%model_id, 2, 3) + patch%centroid_offset(2)
1152 lz(1) = stl_bounding_boxes(patch%model_id, 3, 1) + patch%centroid_offset(3)
1153 lz(2) = stl_bounding_boxes(patch%model_id, 3, 3) + patch%centroid_offset(3)
1154
1155 bound = 0.5_wp*sqrt((lx(2) - lx(1))**2 + (ly(2) - ly(1))**2 + (lz(2) - lz(1))**2)
1156 end if
1157
1158 end subroutine s_get_ib_bound
1159
1160 subroutine s_get_bounding_indices(patch, center, il, ir, jl, jr, kl, kr)
1161
1162
1163# 510 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1164#if MFC_OpenACC
1165# 510 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1166!$acc routine seq
1167# 510 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1168#elif MFC_OpenMP
1169# 510 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1170
1171# 510 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1172
1173# 510 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1174!$omp declare target device_type(any)
1175# 510 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1176#endif
1177
1178 type(ib_patch_parameters), intent(in) :: patch
1179 real(wp), dimension(3), intent(in) :: center
1180 integer, intent(out) :: il, ir, jl, jr, kl, kr
1181 real(wp), dimension(3) :: bbox_min, bbox_max, local_corner, world_corner
1182 real(wp), dimension(2) :: lx, ly, lz
1183 real(wp) :: bound
1184 integer :: cx, cy, cz
1185 logical :: outside_domain
1186
1187 if (patch%geometry == 5) then
1188 ! STL model geometry
1189 lx(1) = stl_bounding_boxes(patch%model_id, 1, 1) + patch%centroid_offset(1)
1190 lx(2) = stl_bounding_boxes(patch%model_id, 1, 3) + patch%centroid_offset(1)
1191 ly(1) = stl_bounding_boxes(patch%model_id, 2, 1) + patch%centroid_offset(2)
1192 ly(2) = stl_bounding_boxes(patch%model_id, 2, 3) + patch%centroid_offset(2)
1193
1194 bbox_min = 1e12
1195 bbox_max = -1e12
1196 ! Enumerate all 4 corners of the local bounding box, rotate to world space, track world-space AABB
1197 do cx = 1, 2
1198 do cy = 1, 2
1199 local_corner = [lx(cx), ly(cy), 0._wp]
1200 world_corner = matmul(patch%rotation_matrix, local_corner) + center
1201 bbox_min(1) = min(bbox_min(1), world_corner(1))
1202 bbox_min(2) = min(bbox_min(2), world_corner(2))
1203 bbox_max(1) = max(bbox_max(1), world_corner(1))
1204 bbox_max(2) = max(bbox_max(2), world_corner(2))
1205 end do
1206 end do
1207 else if (patch%geometry == 12) then
1208 ! Local-space bounding box extents (min=1, max=2 in the third index)
1209 lx(1) = stl_bounding_boxes(patch%model_id, 1, 1) + patch%centroid_offset(1)
1210 lx(2) = stl_bounding_boxes(patch%model_id, 1, 3) + patch%centroid_offset(1)
1211 ly(1) = stl_bounding_boxes(patch%model_id, 2, 1) + patch%centroid_offset(2)
1212 ly(2) = stl_bounding_boxes(patch%model_id, 2, 3) + patch%centroid_offset(2)
1213 lz(1) = stl_bounding_boxes(patch%model_id, 3, 1) + patch%centroid_offset(3)
1214 lz(2) = stl_bounding_boxes(patch%model_id, 3, 3) + patch%centroid_offset(3)
1215
1216 bbox_min = 1e12
1217 bbox_max = -1e12
1218 ! Enumerate all 8 corners of the local bounding box, rotate to world space, track world-space AABB
1219 do cx = 1, 2
1220 do cy = 1, 2
1221 do cz = 1, 2
1222 local_corner = [lx(cx), ly(cy), lz(cz)]
1223 world_corner = matmul(patch%rotation_matrix, local_corner) + center
1224 bbox_min(1) = min(bbox_min(1), world_corner(1))
1225 bbox_min(2) = min(bbox_min(2), world_corner(2))
1226 bbox_min(3) = min(bbox_min(3), world_corner(3))
1227 bbox_max(1) = max(bbox_max(1), world_corner(1))
1228 bbox_max(2) = max(bbox_max(2), world_corner(2))
1229 bbox_max(3) = max(bbox_max(3), world_corner(3))
1230 end do
1231 end do
1232 end do
1233 else
1234 ! All other IBs
1235 call s_get_ib_bound(patch, bound)
1236 bbox_min = center - bound
1237 bbox_max = center + bound
1238 end if
1239
1240 ! completely skip patches whose bounding box does not overlap this rank's domain
1241 outside_domain = bbox_min(1) > x_cc(m + gp_layers + 1) .or. bbox_max(1) < x_cc(-gp_layers - 1) .or. bbox_min(2) > y_cc(n &
1242 & + gp_layers + 1) .or. bbox_max(2) < y_cc(-gp_layers - 1)
1243 if (num_dims == 3) then
1244 outside_domain = outside_domain .or. bbox_min(3) > z_cc(p + gp_layers + 1) .or. bbox_max(3) < z_cc(-gp_layers - 1)
1245 end if
1246
1247 if (outside_domain) then
1248 il = 1; ir = 0
1249 jl = 1; jr = 0
1250 kl = 1; kr = 0
1251 return
1252 end if
1253
1254 il = -gp_layers - 1
1255 jl = -gp_layers - 1
1256 kl = -gp_layers - 1
1257 ir = m + gp_layers + 1
1258 jr = n + gp_layers + 1
1259 kr = p + gp_layers + 1
1260 call get_indices_from_bounds(bbox_min(1), bbox_max(1), x_cc, il, ir)
1261 call get_indices_from_bounds(bbox_min(2), bbox_max(2), y_cc, jl, jr)
1262 if (num_dims == 3) call get_indices_from_bounds(bbox_min(3), bbox_max(3), z_cc, kl, kr)
1263
1264 end subroutine s_get_bounding_indices
1265
1266 subroutine get_indices_from_bounds(left_bound, right_bound, cell_centers, left_index, right_index)
1267
1268
1269# 602 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1270#if MFC_OpenACC
1271# 602 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1272!$acc routine seq
1273# 602 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1274#elif MFC_OpenMP
1275# 602 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1276
1277# 602 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1278
1279# 602 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1280!$omp declare target device_type(any)
1281# 602 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1282#endif
1283
1284 real(wp), intent(in) :: left_bound, right_bound
1285 integer, intent(inout) :: left_index, right_index
1286 real(wp), dimension(-buff_size:), intent(in) :: cell_centers
1287 integer :: itr_left, itr_middle, itr_right
1288
1289 itr_left = left_index
1290 itr_right = right_index
1291
1292 do while (itr_left + 1 < itr_right)
1293 itr_middle = (itr_left + itr_right)/2
1294 if (cell_centers(itr_middle) < left_bound) then
1295 itr_left = itr_middle
1296 else if (cell_centers(itr_middle) > left_bound) then
1297 itr_right = itr_middle
1298 else
1299 itr_left = itr_middle
1300 exit
1301 end if
1302 end do
1303 left_index = itr_left
1304
1305 itr_right = right_index
1306 do while (itr_left + 1 < itr_right)
1307 itr_middle = (itr_left + itr_right)/2
1308 if (cell_centers(itr_middle) < right_bound) then
1309 itr_left = itr_middle
1310 else if (cell_centers(itr_middle) > right_bound) then
1311 itr_right = itr_middle
1312 else
1313 itr_right = itr_middle
1314 exit
1315 end if
1316 end do
1317 right_index = itr_right
1318
1319 end subroutine get_indices_from_bounds
1320
1321 !> Encode the patch ID with a unique offset containing periodicity information
1322 subroutine s_encode_patch_periodicity(patch_id, x_periodicity, y_periodicity, z_periodicity, encoded_patch_id)
1323
1324
1325# 644 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1326#if MFC_OpenACC
1327# 644 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1328!$acc routine seq
1329# 644 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1330#elif MFC_OpenMP
1331# 644 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1332
1333# 644 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1334
1335# 644 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1336!$omp declare target device_type(any)
1337# 644 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1338#endif
1339
1340 integer, intent(in) :: patch_id, x_periodicity, y_periodicity, z_periodicity
1341 integer, intent(out) :: encoded_patch_id
1342 integer :: temp_x_per, temp_y_per, temp_z_per, offset
1343
1344 encoded_patch_id = patch_id
1345
1346 temp_x_per = x_periodicity; if (x_periodicity == -1) temp_x_per = 2
1347 temp_y_per = y_periodicity; if (y_periodicity == -1) temp_y_per = 2
1348 temp_z_per = z_periodicity; if (z_periodicity == -1) temp_z_per = 2
1349
1350 offset = (num_gbl_ibs + 1)*temp_x_per + 3*(num_gbl_ibs + 1)*temp_y_per + 9*(num_gbl_ibs + 1)*temp_z_per
1351 encoded_patch_id = patch_id + offset
1352
1353 end subroutine s_encode_patch_periodicity
1354
1355 !> Decode the encoded ID to recover the original patch ID and periodicity
1356 subroutine s_decode_patch_periodicity(encoded_patch_id, patch_id, x_periodicity, y_periodicity, z_periodicity)
1357
1358
1359# 664 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1360#if MFC_OpenACC
1361# 664 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1362!$acc routine seq
1363# 664 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1364#elif MFC_OpenMP
1365# 664 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1366
1367# 664 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1368
1369# 664 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1370!$omp declare target device_type(any)
1371# 664 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1372#endif
1373
1374 integer, intent(in) :: encoded_patch_id
1375 integer, intent(out) :: patch_id
1376 integer, intent(out), optional :: x_periodicity, y_periodicity, z_periodicity
1377 integer :: offset, remainder, xp, yp, zp, base
1378
1379 base = num_gbl_ibs + 1
1380
1381 patch_id = mod(encoded_patch_id - 1, base) + 1
1382 offset = (encoded_patch_id - patch_id)/base
1383
1384 xp = mod(offset, 3)
1385 remainder = offset/3
1386 yp = mod(remainder, 3)
1387 zp = remainder/3
1388
1389 ! Reverse map: 2 -> -1, 0 -> 0, 1 -> 1
1390 if (present(x_periodicity) .and. present(y_periodicity) .and. present(z_periodicity)) then
1391 x_periodicity = xp; if (xp == 2) x_periodicity = -1
1392 y_periodicity = yp; if (yp == 2) y_periodicity = -1
1393 z_periodicity = zp; if (zp == 2) z_periodicity = -1
1394 end if
1395
1396 end subroutine s_decode_patch_periodicity
1397
1398 !> Determine the periodic wrapping bounds in each direction
1399 subroutine s_get_periodicities(xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper)
1400
1401 integer, intent(out) :: xp_lower, xp_upper, yp_lower, yp_upper
1402 integer, intent(out), optional :: zp_lower, zp_upper
1403
1404 ! check domain wraps in x, y
1405
1406# 699 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1407 ! check for periodicity
1408 if (ib_bc_x%beg == bc_periodic) then
1409 xp_lower = -1
1410 xp_upper = 1
1411 else
1412 ! if it is not periodic, then both elements are 0
1413 xp_lower = 0
1414 xp_upper = 0
1415 end if
1416# 699 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1417 ! check for periodicity
1418 if (ib_bc_y%beg == bc_periodic) then
1419 yp_lower = -1
1420 yp_upper = 1
1421 else
1422 ! if it is not periodic, then both elements are 0
1423 yp_lower = 0
1424 yp_upper = 0
1425 end if
1426# 709 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1427
1428 ! z only if 3D
1429 if (present(zp_lower) .and. num_dims == 3) then
1430 if (ib_bc_z%beg == bc_periodic) then
1431 zp_lower = -1
1432 zp_upper = 1
1433 else
1434 zp_lower = 0
1435 zp_upper = 0
1436 end if
1437 else if (present(zp_lower)) then
1438 zp_lower = 0
1439 zp_upper = 0
1440 end if
1441
1442 end subroutine s_get_periodicities
1443
1444 !> Archimedes spiral function
1445 pure elemental function f_r(myth, offset, a)
1446
1447
1448# 729 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1449#if MFC_OpenACC
1450# 729 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1451!$acc routine seq
1452# 729 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1453#elif MFC_OpenMP
1454# 729 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1455
1456# 729 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1457
1458# 729 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1459!$omp declare target device_type(any)
1460# 729 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1461#endif
1462 real(wp), intent(in) :: myth, offset, a
1463 real(wp) :: b
1464 real(wp) :: f_r
1465
1466 ! r(th) = a + b*th
1467
1468 b = 2._wp*a/(2._wp*pi)
1469 f_r = a + b*myth + offset
1470
1471 end function f_r
1472
1473end module m_ib_patches
j
integer, intent(in) j
Definition m_phase_change.fpp.f90:1008
m_derived_types
Shared derived types for field data, patch geometry, bubble dynamics, and MPI I/O structures.
Definition m_derived_types.fpp.f90:335
m_global_parameters
Global parameters for the computational domain, fluid properties, and simulation algorithm configurat...
Definition m_global_parameters.fpp.f90:346
m_global_parameters::glb_bounds
type(bounds_info), dimension(3) glb_bounds
Definition m_global_parameters.fpp.f90:405
m_global_parameters::y_cc
real(wp), dimension(:), allocatable, target y_cc
Definition m_global_parameters.fpp.f90:410
m_global_parameters::z_cc
real(wp), dimension(:), allocatable, target z_cc
Definition m_global_parameters.fpp.f90:410
m_global_parameters::x_cc
real(wp), dimension(:), allocatable, target x_cc
Definition m_global_parameters.fpp.f90:410
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:864
m_global_parameters::dx
real(wp), dimension(:), allocatable, target dx
Definition m_global_parameters.fpp.f90:414
m_helper_basic
Basic floating-point utilities: approximate equality, default detection, and coordinate bounds.
Definition m_helper_basic.fpp.f90:335
m_helper
Utility routines for bubble model setup, coordinate transforms, array sampling, and special functions...
Definition m_helper.fpp.f90:346
m_ib_patches
Allocate memory and read initial condition data for IC extrusion.
Definition m_ib_patches.fpp.f90:405
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:1040
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:1267
m_ib_patches::s_get_ib_bound
subroutine, public s_get_ib_bound(patch, bound)
Definition m_ib_patches.fpp.f90:1100
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:1161
m_ib_patches::s_apply_ib_patches_ib_parallelism
subroutine s_apply_ib_patches_ib_parallelism(ib_markers)
Definition m_ib_patches.fpp.f90:662
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:436
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:424
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:1400
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:1357
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:873
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:1323
m_ib_patches::f_r
pure elemental real(wp) function f_r(myth, offset, a)
Archimedes spiral function.
Definition m_ib_patches.fpp.f90:1446
m_model
Binary STL file reader and processor for immersed boundary geometry.
Definition m_model.fpp.f90:336
m_model::gpu_ntrs
integer, dimension(:), allocatable, public gpu_ntrs
GPU-friendly flat arrays for STL model data.
Definition m_model.fpp.f90:357
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:804
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:1275
m_mpi_common
MPI communication layer: domain decomposition, halo exchange, reductions, and parallel I/O setup.
Definition m_mpi_common.fpp.f90:346
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:336
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:354
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:411
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:381
m_patch_geometries::f_is_inside_ellipse
logical function, public f_is_inside_ellipse(x, y, length)
Definition m_patch_geometries.fpp.f90:514
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:443
m_derived_types::integer_field
Derived type annexing an integer scalar field (SF).
Definition m_derived_types.fpp.f90:359