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