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# 9 "/home/runner/work/MFC/MFC/src/common/include/case.fpp"
13
14! For moving immersed boundaries in simulation
15# 14 "/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# 74 "/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# 26 "/home/runner/work/MFC/MFC/src/common/include/2dHardcodedIC.fpp"
68
69# 349 "/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! New line at end of file is required for FYPP
103# 2 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
104# 1 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 1
105# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
106# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
107# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
108# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
109# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
110# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
111
112# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
113# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
114# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
115
116# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
117
118# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
119
120# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
121
122# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
123
124# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
125
126# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
127
128# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
129! New line at end of file is required for FYPP
130# 2 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 2
131
132# 4 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
133# 5 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
134# 6 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
135# 7 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
136# 8 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
137
138# 20 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
139
140# 43 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
141
142# 48 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
143
144# 53 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
145
146# 58 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
147
148# 63 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
149
150# 68 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
151
152# 76 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
153
154# 81 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
155
156# 86 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
157
158# 91 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
159
160# 96 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
161
162# 101 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
163
164# 106 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
165
166# 111 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
167
168# 116 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
169
170# 121 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
171
172# 151 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
173
174# 192 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
175
176# 206 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
177
178# 231 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
179
180# 242 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
181
182# 244 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
183# 255 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
184
185# 284 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
186
187# 294 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
188
189# 304 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
190
191# 313 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
192
193# 330 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
194
195# 340 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
196
197# 347 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
198
199# 353 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
200
201# 359 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
202
203# 365 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
204
205# 371 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
206
207# 377 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
208! New line at end of file is required for FYPP
209# 3 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
210# 1 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 1
211# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
212# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
213# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
214# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
215# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
216# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
217
218# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
219# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
220# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
221
222# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
223
224# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
225
226# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
227
228# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
229
230# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
231
232# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
233
234# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
235! New line at end of file is required for FYPP
236# 2 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 2
237
238# 7 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
239
240# 17 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
241
242# 22 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
243
244# 27 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
245
246# 32 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
247
248# 37 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
249
250# 42 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
251
252# 47 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
253
254# 52 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
255
256# 57 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
257
258# 62 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
259
260# 73 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
261
262# 78 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
263
264# 83 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
265
266# 88 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
267
268# 103 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
269
270# 131 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
271
272# 160 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
273
274# 175 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
275
276# 193 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
277
278# 215 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
279
280# 244 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
281
282# 259 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
283
284# 269 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
285
286# 278 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
287
288# 294 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
289
290# 304 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
291
292# 311 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
293! New line at end of file is required for FYPP
294# 4 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
295
296! GPU parallel region (scalar reductions, maxval/minval)
297# 23 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
298
299! GPU parallel loop over threads (most common GPU macro)
300# 43 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
301
302! Required closing for GPU_PARALLEL_LOOP
303# 55 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
304
305! Mark routine for device compilation
306# 112 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
307
308! Declare device-resident data
309# 130 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
310
311! Inner loop within a GPU parallel region
312# 145 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
313
314! Scoped GPU data region
315# 164 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
316
317! Host code with device pointers (for MPI with GPU buffers)
318# 193 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
319
320! Allocate device memory (unscoped)
321# 207 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
322
323! Free device memory
324# 219 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
325
326! Atomic operation on device
327# 231 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
328
329! End atomic capture block
330# 242 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
331
332! Copy data between host and device
333# 254 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
334
335! Synchronization barrier
336# 266 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
337
338! Import GPU library module (openacc or omp_lib)
339# 275 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
340
341! Emit code only for AMD compiler
342# 282 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
343
344! Emit code for non-Cray compilers
345# 289 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
346
347! Emit code only for Cray compiler
348# 296 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
349
350! Emit code for non-NVIDIA compilers
351# 303 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
352
353# 305 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
354# 306 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
355! New line at end of file is required for FYPP
356# 2 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 2
357
358# 14 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
359
360! Caution: This macro requires the use of a binding script to set CUDA_VISIBLE_DEVICES, such that we have one GPU device per MPI
361! rank. That's because for both cudaMemAdvise (preferred location) and cudaMemPrefetchAsync we use location = device_id = 0. For an
362! example see misc/nvidia_uvm/bind.sh. NVIDIA unified memory page placement hint
363# 57 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
364
365! Allocate and create GPU device memory
366# 77 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
367
368! Free GPU device memory and deallocate
369# 85 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
370
371! Cray-specific GPU pointer setup for vector fields
372# 109 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
373
374! Cray-specific GPU pointer setup for scalar fields
375# 125 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
376
377! Cray-specific GPU pointer setup for acoustic source spatials
378# 150 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
379
380# 156 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
381
382# 163 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
383! New line at end of file is required for FYPP
384# 11 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp" 2
385
386!> @brief Immersed boundary patch geometry constructors for 2D and 3D shapes
387module m_ib_patches
388
389 use m_model ! Subroutine(s) related to STL files
390 use m_derived_types ! Definitions of the derived types
391 use m_global_parameters
392 use m_helper_basic
393 use m_helper
394 use m_mpi_common
395
396 implicit none
397
398 private; public :: s_apply_ib_patches, s_update_ib_rotation_matrix, f_convert_cyl_to_cart, s_instantiate_stl_models, &
399 & s_decode_patch_periodicity
400
401 real(wp) :: cart_y, cart_z
402
403# 28 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
404#if defined(MFC_OpenACC)
405# 28 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
406!$acc declare create(cart_y, cart_z)
407# 28 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
408#elif defined(MFC_OpenMP)
409# 28 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
410!$omp declare target (cart_y, cart_z)
411# 28 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
412#endif
413 ! Variables to be used to hold cell locations in Cartesian coordinates if 3D simulation is using cylindrical coordinates
414
415contains
416
417 !> Apply all immersed boundary patch geometries to mark interior cells in the IB marker array
418 impure subroutine s_apply_ib_patches(ib_markers)
419
420 type(integer_field), intent(inout) :: ib_markers
421 integer :: i, xp, yp, zp !< iterators
422 integer :: xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper !< periodic bounds
423
424 ! 3D Patch Geometries
425
426 if (p > 0) then
427 !> IB Patches
428 !> @{
429 call s_get_periodicities(xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper)
430 do xp = xp_lower, xp_upper
431 do yp = yp_lower, yp_upper
432 do zp = zp_lower, zp_upper
433 do i = 1, num_ibs
434 if (patch_ib(i)%geometry == 8) then
435 call s_ib_sphere(i, ib_markers, xp, yp, zp)
436 else if (patch_ib(i)%geometry == 9) then
437 call s_ib_cuboid(i, ib_markers, xp, yp, zp)
438 else if (patch_ib(i)%geometry == 10) then
439 call s_ib_cylinder(i, ib_markers, xp, yp, zp)
440 else if (patch_ib(i)%geometry == 11) then
441 call s_ib_3d_airfoil(i, ib_markers, xp, yp, zp)
442 else if (patch_ib(i)%geometry == 12) then
443 call s_ib_3d_model(i, ib_markers, xp, yp, zp)
444 end if
445 end do
446 end do
447 end do
448 end do
449 !> @}
450
451 ! 2D Patch Geometries
452 else if (n > 0) then
453 !> IB Patches
454 !> @{
455 call s_get_periodicities(xp_lower, xp_upper, yp_lower, yp_upper)
456 do xp = xp_lower, xp_upper
457 do yp = yp_lower, yp_upper
458 do i = 1, num_ibs
459 if (patch_ib(i)%geometry == 2) then
460 call s_ib_circle(i, ib_markers, xp, yp)
461 else if (patch_ib(i)%geometry == 3) then
462 call s_ib_rectangle(i, ib_markers, xp, yp)
463 else if (patch_ib(i)%geometry == 4) then
464 call s_ib_airfoil(i, ib_markers, xp, yp)
465 else if (patch_ib(i)%geometry == 5) then
466 call s_ib_model(i, ib_markers, xp, yp)
467 else if (patch_ib(i)%geometry == 6) then
468 call s_ib_ellipse(i, ib_markers, xp, yp)
469 end if
470 end do
471 end do
472 end do
473 !> @}
474 end if
475
476 end subroutine s_apply_ib_patches
477
478 !> Mark cells inside a circular immersed boundary
479 subroutine s_ib_circle(patch_id, ib_markers, xp, yp)
480
481 integer, intent(in) :: patch_id
482 integer, intent(in) :: xp, yp !< integers containing the periodicity projection information
483 type(integer_field), intent(inout) :: ib_markers
484 real(wp), dimension(1:2) :: center
485 real(wp) :: radius
486 integer :: i, j, il, ir, jl, jr !< Generic loop iterators
487 integer :: encoded_patch_id
488
489 ! Transferring the circular patch's radius, centroid, smearing patch identity and smearing coefficient information
490
491 center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg)
492 center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg)
493 radius = patch_ib(patch_id)%radius
494
495 ! encode the periodicity information into the patch_id
496 call s_encode_patch_periodicity(patch_id, xp, yp, 0, encoded_patch_id)
497
498 ! find the indices to the left and right of the IB in i, j, k
499 il = -gp_layers - 1
500 jl = -gp_layers - 1
501 ir = m + gp_layers + 1
502 jr = n + gp_layers + 1
503 call get_bounding_indices(center(1) - radius, center(1) + radius, x_cc, il, ir)
504 call get_bounding_indices(center(2) - radius, center(2) + radius, y_cc, jl, jr)
505
506 ! Assign primitive variables if circle covers cell and patch has write permission
507
508
509# 124 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
510
511# 124 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
512#if defined(MFC_OpenACC)
513# 124 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
514!$acc parallel loop collapse(2) gang vector default(present) private(i, j) copyin(encoded_patch_id, center, radius)
515# 124 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
516#elif defined(MFC_OpenMP)
517# 124 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
518
519# 124 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
520
521# 124 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
522
523# 124 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
524!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(2) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, j) map(to:encoded_patch_id, center, radius)
525# 124 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
526#endif
527 do j = jl, jr
528 do i = il, ir
529 if ((x_cc(i) - center(1))**2 + (y_cc(j) - center(2))**2 <= radius**2) then
530 ib_markers%sf(i, j, 0) = encoded_patch_id
531 end if
532 end do
533 end do
534
535# 132 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
536#if defined(MFC_OpenACC)
537# 132 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
538!$acc end parallel loop
539# 132 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
540#elif defined(MFC_OpenMP)
541# 132 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
542
543# 132 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
544!$omp end target teams loop
545# 132 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
546#endif
547
548 end subroutine s_ib_circle
549
550 !> Mark cells inside a 2D NACA 4-digit airfoil immersed boundary
551 subroutine s_ib_airfoil(patch_id, ib_markers, xp, yp)
552
553 integer, intent(in) :: patch_id
554 type(integer_field), intent(inout) :: ib_markers
555 integer, intent(in) :: xp, yp !< integers containing the periodicity projection information
556 real(wp) :: f, ca_in, pa, ma, ta
557 real(wp) :: xa, yt, xu, yu, xl, yl, xc, yc, dycdxc, sin_c, cos_c
558 integer :: i, j, k, il, ir, jl, jr
559 integer :: Np1, Np2
560 integer :: encoded_patch_id
561 real(wp), dimension(1:3) :: xy_local, offset !< x and y coordinates in local IB frame
562 real(wp), dimension(1:2) :: center !< x and y coordinates in local IB frame
563 real(wp), dimension(1:3,1:3) :: inverse_rotation
564
565 center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg)
566 center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg)
567 ca_in = patch_ib(patch_id)%c
568 pa = patch_ib(patch_id)%p
569 ma = patch_ib(patch_id)%m
570 ta = patch_ib(patch_id)%t
571 inverse_rotation(:,:) = patch_ib(patch_id)%rotation_matrix_inverse(:,:)
572 offset(:) = patch_ib(patch_id)%centroid_offset(:)
573
574 np1 = int((pa*ca_in/dx(0))*20)
575 np2 = int(((ca_in - pa*ca_in)/dx(0))*20)
576 np = np1 + np2 + 1
577
578# 163 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
579#if defined(MFC_OpenACC)
580# 163 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
581!$acc update device(Np)
582# 163 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
583#elif defined(MFC_OpenMP)
584# 163 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
585!$omp target update to(Np)
586# 163 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
587#endif
588
589 if (.not. allocated(airfoil_grid_u)) then
590#ifdef MFC_DEBUG
591# 166 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
592 block
593# 166 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
594 use iso_fortran_env, only: output_unit
595# 166 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
596
597# 166 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
598 print *, 'm_ib_patches.fpp:166: ', '@:ALLOCATE(airfoil_grid_u(1:Np))'
599# 166 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
600
601# 166 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
602 call flush (output_unit)
603# 166 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
604 end block
605# 166 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
606#endif
607# 166 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
608 allocate (airfoil_grid_u(1:np))
609# 166 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
610
611# 166 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
612
613# 166 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
614#if defined(MFC_OpenACC)
615# 166 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
616!$acc enter data create(airfoil_grid_u)
617# 166 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
618#elif defined(MFC_OpenMP)
619# 166 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
620!$omp target enter data map(always,alloc:airfoil_grid_u)
621# 166 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
622#endif
623#ifdef MFC_DEBUG
624# 167 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
625 block
626# 167 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
627 use iso_fortran_env, only: output_unit
628# 167 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
629
630# 167 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
631 print *, 'm_ib_patches.fpp:167: ', '@:ALLOCATE(airfoil_grid_l(1:Np))'
632# 167 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
633
634# 167 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
635 call flush (output_unit)
636# 167 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
637 end block
638# 167 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
639#endif
640# 167 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
641 allocate (airfoil_grid_l(1:np))
642# 167 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
643
644# 167 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
645
646# 167 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
647#if defined(MFC_OpenACC)
648# 167 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
649!$acc enter data create(airfoil_grid_l)
650# 167 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
651#elif defined(MFC_OpenMP)
652# 167 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
653!$omp target enter data map(always,alloc:airfoil_grid_l)
654# 167 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
655#endif
656
657 ! TODO :: The below instantiations are already handled by the loop below
658 airfoil_grid_u(1)%x = 0._wp
659 airfoil_grid_u(1)%y = 0._wp
660
661 airfoil_grid_l(1)%x = 0._wp
662 airfoil_grid_l(1)%y = 0._wp
663
664 do i = 1, np1 + np2 - 1
665 ! TODO :: This allocates the upper and lower airfoil arrays, and does not need to be performed each time the IB
666 ! markers are updated. Place this as a separate subroutine.
667 if (i <= np1) then
668 xc = i*(pa*ca_in/np1)
669 xa = xc/ca_in
670 yc = (ma/pa**2)*(2*pa*xa - xa**2)
671 dycdxc = (2*ma/pa**2)*(pa - xa)
672 else
673 xc = pa*ca_in + (i - np1)*((ca_in - pa*ca_in)/np2)
674 xa = xc/ca_in
675 yc = (ma/(1 - pa)**2)*(1 - 2*pa + 2*pa*xa - xa**2)
676 dycdxc = (2*ma/(1 - pa)**2)*(pa - xa)
677 end if
678
679 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)
680 sin_c = dycdxc/(1 + dycdxc**2)**0.5_wp
681 cos_c = 1/(1 + dycdxc**2)**0.5_wp
682
683 xu = xa - yt*sin_c
684 yu = yc + yt*cos_c
685
686 xl = xa + yt*sin_c
687 yl = yc - yt*cos_c
688
689 xu = xu*ca_in
690 yu = yu*ca_in
691
692 xl = xl*ca_in
693 yl = yl*ca_in
694
695 airfoil_grid_u(i + 1)%x = xu
696 airfoil_grid_u(i + 1)%y = yu
697
698 airfoil_grid_l(i + 1)%x = xl
699 airfoil_grid_l(i + 1)%y = yl
700 end do
701
702 airfoil_grid_u(np)%x = ca_in
703 airfoil_grid_u(np)%y = 0._wp
704
705 airfoil_grid_l(np)%x = ca_in
706 airfoil_grid_l(np)%y = 0._wp
707
708
709# 220 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
710#if defined(MFC_OpenACC)
711# 220 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
712!$acc update device(airfoil_grid_l, airfoil_grid_u)
713# 220 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
714#elif defined(MFC_OpenMP)
715# 220 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
716!$omp target update to(airfoil_grid_l, airfoil_grid_u)
717# 220 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
718#endif
719 end if
720
721 ! encode the periodicity information into the patch_id
722 call s_encode_patch_periodicity(patch_id, xp, yp, 0, encoded_patch_id)
723
724 ! find the indices to the left and right of the IB in i, j, k
725 il = -gp_layers - 1
726 jl = -gp_layers - 1
727 ir = m + gp_layers + 1
728 jr = n + gp_layers + 1
729 ! maximum distance any marker can be from the center is the length of the airfoil
730 call get_bounding_indices(center(1) - ca_in, center(1) + ca_in, x_cc, il, ir)
731 call get_bounding_indices(center(2) - ca_in, center(2) + ca_in, y_cc, jl, jr)
732
733
734# 235 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
735
736# 235 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
737#if defined(MFC_OpenACC)
738# 235 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
739!$acc parallel loop collapse(2) gang vector default(present) private(i, j, xy_local, k, f) copyin(encoded_patch_id, center, inverse_rotation, offset, ma, ca_in, airfoil_grid_u, airfoil_grid_l)
740# 235 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
741#elif defined(MFC_OpenMP)
742# 235 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
743
744# 235 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
745
746# 235 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
747
748# 235 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
749!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(2) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, j, xy_local, k, f) map(to:encoded_patch_id, center, inverse_rotation, offset, ma, ca_in, airfoil_grid_u, airfoil_grid_l)
750# 235 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
751#endif
752# 237 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
753 do j = jl, jr
754 do i = il, ir
755 xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp] ! get coordinate frame centered on IB
756 xy_local = matmul(inverse_rotation, xy_local) ! rotate the frame into the IB's coordinates
757 xy_local = xy_local - offset ! airfoils are a patch that require a centroid offset
758
759 if (xy_local(1) >= 0._wp .and. xy_local(1) <= ca_in) then
760 xa = xy_local(1)/ca_in
761 if (xa <= pa) then
762 yc = (ma/pa**2)*(2*pa*xa - xa**2)
763 dycdxc = (2*ma/pa**2)*(pa - xa)
764 else
765 yc = (ma/(1 - pa)**2)*(1 - 2*pa + 2*pa*xa - xa**2)
766 dycdxc = (2*ma/(1 - pa)**2)*(pa - xa)
767 end if
768 if (xy_local(2) >= 0._wp) then
769 k = 1
770 do while (airfoil_grid_u(k)%x < xy_local(1) .and. k <= np)
771 k = k + 1
772 end do
773 if (f_approx_equal(airfoil_grid_u(k)%x, xy_local(1))) then
774 if (xy_local(2) <= airfoil_grid_u(k)%y) then
775 ib_markers%sf(i, j, 0) = encoded_patch_id
776 end if
777 else
778 f = (airfoil_grid_u(k)%x - xy_local(1))/(airfoil_grid_u(k)%x - airfoil_grid_u(k - 1)%x)
779 if (xy_local(2) <= ((1._wp - f)*airfoil_grid_u(k)%y + f*airfoil_grid_u(k - 1)%y)) then
780 ib_markers%sf(i, j, 0) = encoded_patch_id
781 end if
782 end if
783 else
784 k = 1
785 do while (airfoil_grid_l(k)%x < xy_local(1))
786 k = k + 1
787 end do
788 if (f_approx_equal(airfoil_grid_l(k)%x, xy_local(1))) then
789 if (xy_local(2) >= airfoil_grid_l(k)%y) then
790 ib_markers%sf(i, j, 0) = encoded_patch_id
791 end if
792 else
793 f = (airfoil_grid_l(k)%x - xy_local(1))/(airfoil_grid_l(k)%x - airfoil_grid_l(k - 1)%x)
794
795 if (xy_local(2) >= ((1._wp - f)*airfoil_grid_l(k)%y + f*airfoil_grid_l(k - 1)%y)) then
796 ib_markers%sf(i, j, 0) = encoded_patch_id
797 end if
798 end if
799 end if
800 end if
801 end do
802 end do
803
804# 287 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
805#if defined(MFC_OpenACC)
806# 287 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
807!$acc end parallel loop
808# 287 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
809#elif defined(MFC_OpenMP)
810# 287 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
811
812# 287 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
813!$omp end target teams loop
814# 287 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
815#endif
816
817 end subroutine s_ib_airfoil
818
819 !> Mark cells inside a 3D extruded NACA 4-digit airfoil immersed boundary with finite span
820 subroutine s_ib_3d_airfoil(patch_id, ib_markers, xp, yp, zp)
821
822 integer, intent(in) :: patch_id
823 type(integer_field), intent(inout) :: ib_markers
824 integer, intent(in) :: xp, yp, zp !< integers containing the periodicity projection information
825 real(wp) :: lz, z_max, z_min, f, ca_in, pa, ma, ta, xa, yt, xu, yu, xl, yl, xc, yc, dycdxc, sin_c, cos_c
826 integer :: i, j, k, l, il, ir, jl, jr, ll, lr
827 integer :: Np1, Np2
828 integer :: encoded_patch_id
829 real(wp), dimension(1:3) :: xyz_local, center, offset !< x, y, z coordinates in local IB frame
830 real(wp), dimension(1:3,1:3) :: inverse_rotation
831
832 center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg)
833 center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg)
834 center(3) = patch_ib(patch_id)%z_centroid + real(zp, wp)*(z_domain%end - z_domain%beg)
835 lz = patch_ib(patch_id)%length_z
836 ca_in = patch_ib(patch_id)%c
837 pa = patch_ib(patch_id)%p
838 ma = patch_ib(patch_id)%m
839 ta = patch_ib(patch_id)%t
840 inverse_rotation(:,:) = patch_ib(patch_id)%rotation_matrix_inverse(:,:)
841 offset(:) = patch_ib(patch_id)%centroid_offset(:)
842
843 z_max = lz/2
844 z_min = -lz/2
845
846 np1 = int((pa*ca_in/dx(0))*20)
847 np2 = int(((ca_in - pa*ca_in)/dx(0))*20)
848 np = np1 + np2 + 1
849
850# 321 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
851#if defined(MFC_OpenACC)
852# 321 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
853!$acc update device(Np)
854# 321 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
855#elif defined(MFC_OpenMP)
856# 321 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
857!$omp target update to(Np)
858# 321 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
859#endif
860
861 if (.not. allocated(airfoil_grid_u)) then
862#ifdef MFC_DEBUG
863# 324 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
864 block
865# 324 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
866 use iso_fortran_env, only: output_unit
867# 324 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
868
869# 324 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
870 print *, 'm_ib_patches.fpp:324: ', '@:ALLOCATE(airfoil_grid_u(1:Np))'
871# 324 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
872
873# 324 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
874 call flush (output_unit)
875# 324 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
876 end block
877# 324 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
878#endif
879# 324 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
880 allocate (airfoil_grid_u(1:np))
881# 324 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
882
883# 324 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
884
885# 324 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
886#if defined(MFC_OpenACC)
887# 324 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
888!$acc enter data create(airfoil_grid_u)
889# 324 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
890#elif defined(MFC_OpenMP)
891# 324 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
892!$omp target enter data map(always,alloc:airfoil_grid_u)
893# 324 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
894#endif
895#ifdef MFC_DEBUG
896# 325 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
897 block
898# 325 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
899 use iso_fortran_env, only: output_unit
900# 325 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
901
902# 325 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
903 print *, 'm_ib_patches.fpp:325: ', '@:ALLOCATE(airfoil_grid_l(1:Np))'
904# 325 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
905
906# 325 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
907 call flush (output_unit)
908# 325 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
909 end block
910# 325 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
911#endif
912# 325 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
913 allocate (airfoil_grid_l(1:np))
914# 325 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
915
916# 325 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
917
918# 325 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
919#if defined(MFC_OpenACC)
920# 325 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
921!$acc enter data create(airfoil_grid_l)
922# 325 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
923#elif defined(MFC_OpenMP)
924# 325 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
925!$omp target enter data map(always,alloc:airfoil_grid_l)
926# 325 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
927#endif
928
929 airfoil_grid_u(1)%x = 0._wp
930 airfoil_grid_u(1)%y = 0._wp
931
932 airfoil_grid_l(1)%x = 0._wp
933 airfoil_grid_l(1)%y = 0._wp
934
935 do i = 1, np1 + np2 - 1
936 if (i <= np1) then
937 xc = i*(pa*ca_in/np1)
938 xa = xc/ca_in
939 yc = (ma/pa**2)*(2*pa*xa - xa**2)
940 dycdxc = (2*ma/pa**2)*(pa - xa)
941 else
942 xc = pa*ca_in + (i - np1)*((ca_in - pa*ca_in)/np2)
943 xa = xc/ca_in
944 yc = (ma/(1 - pa)**2)*(1 - 2*pa + 2*pa*xa - xa**2)
945 dycdxc = (2*ma/(1 - pa)**2)*(pa - xa)
946 end if
947
948 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)
949 sin_c = dycdxc/(1 + dycdxc**2)**0.5_wp
950 cos_c = 1/(1 + dycdxc**2)**0.5_wp
951
952 xu = xa - yt*sin_c
953 yu = yc + yt*cos_c
954
955 xl = xa + yt*sin_c
956 yl = yc - yt*cos_c
957
958 xu = xu*ca_in
959 yu = yu*ca_in
960
961 xl = xl*ca_in
962 yl = yl*ca_in
963
964 airfoil_grid_u(i + 1)%x = xu
965 airfoil_grid_u(i + 1)%y = yu
966
967 airfoil_grid_l(i + 1)%x = xl
968 airfoil_grid_l(i + 1)%y = yl
969 end do
970
971 airfoil_grid_u(np)%x = ca_in
972 airfoil_grid_u(np)%y = 0._wp
973
974 airfoil_grid_l(np)%x = ca_in
975 airfoil_grid_l(np)%y = 0._wp
976
977
978# 375 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
979#if defined(MFC_OpenACC)
980# 375 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
981!$acc update device(airfoil_grid_l, airfoil_grid_u)
982# 375 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
983#elif defined(MFC_OpenMP)
984# 375 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
985!$omp target update to(airfoil_grid_l, airfoil_grid_u)
986# 375 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
987#endif
988 end if
989
990 ! encode the periodicity information into the patch_id
991 call s_encode_patch_periodicity(patch_id, xp, yp, zp, encoded_patch_id)
992
993 ! find the indices to the left and right of the IB in i, j, k
994 il = -gp_layers - 1
995 jl = -gp_layers - 1
996 ll = -gp_layers - 1
997 ir = m + gp_layers + 1
998 jr = n + gp_layers + 1
999 lr = p + gp_layers + 1
1000 ! maximum distance any marker can be from the center is the length of the airfoil
1001 call get_bounding_indices(center(1) - ca_in, center(1) + ca_in, x_cc, il, ir)
1002 call get_bounding_indices(center(2) - ca_in, center(2) + ca_in, y_cc, jl, jr)
1003 call get_bounding_indices(center(3) - ca_in, center(3) + ca_in, z_cc, ll, lr)
1004
1005
1006# 393 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1007
1008# 393 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1009#if defined(MFC_OpenACC)
1010# 393 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1011!$acc parallel loop collapse(3) gang vector default(present) private(i, j, l, xyz_local, k, f) copyin(encoded_patch_id, center, inverse_rotation, offset, ma, ca_in, airfoil_grid_u, airfoil_grid_l, z_min, z_max)
1012# 393 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1013#elif defined(MFC_OpenMP)
1014# 393 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1015
1016# 393 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1017
1018# 393 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1019
1020# 393 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1021!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, j, l, xyz_local, k, f) map(to:encoded_patch_id, center, inverse_rotation, offset, ma, ca_in, airfoil_grid_u, airfoil_grid_l, z_min, z_max)
1022# 393 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1023#endif
1024# 395 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1025 do l = ll, lr
1026 do j = jl, jr
1027 do i = il, ir
1028 xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), &
1029 & z_cc(l) - center(3)] ! get coordinate frame centered on IB
1030 xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinates
1031 xyz_local = xyz_local - offset ! airfoils are a patch that require a centroid offset
1032
1033 if (xyz_local(3) >= z_min .and. xyz_local(3) <= z_max) then
1034 if (xyz_local(1) >= 0._wp .and. xyz_local(1) <= ca_in) then
1035 if (xyz_local(2) >= 0._wp) then
1036 k = 1
1037 do while (airfoil_grid_u(k)%x < xyz_local(1))
1038 k = k + 1
1039 end do
1040 if (f_approx_equal(airfoil_grid_u(k)%x, xyz_local(1))) then
1041 if (xyz_local(2) <= airfoil_grid_u(k)%y) then
1042 ! IB
1043 ib_markers%sf(i, j, l) = encoded_patch_id
1044 end if
1045 else
1046 f = (airfoil_grid_u(k)%x - xyz_local(1))/(airfoil_grid_u(k)%x - airfoil_grid_u(k - 1)%x)
1047 if (xyz_local(2) <= ((1._wp - f)*airfoil_grid_u(k)%y + f*airfoil_grid_u(k - 1)%y)) then
1048 ib_markers%sf(i, j, l) = encoded_patch_id
1049 end if
1050 end if
1051 else
1052 k = 1
1053 do while (airfoil_grid_l(k)%x < xyz_local(1))
1054 k = k + 1
1055 end do
1056 if (f_approx_equal(airfoil_grid_l(k)%x, xyz_local(1))) then
1057 if (xyz_local(2) >= airfoil_grid_l(k)%y) then
1058 ib_markers%sf(i, j, l) = encoded_patch_id
1059 end if
1060 else
1061 f = (airfoil_grid_l(k)%x - xyz_local(1))/(airfoil_grid_l(k)%x - airfoil_grid_l(k - 1)%x)
1062
1063 if (xyz_local(2) >= ((1._wp - f)*airfoil_grid_l(k)%y + f*airfoil_grid_l(k - 1)%y)) then
1064 ib_markers%sf(i, j, l) = encoded_patch_id
1065 end if
1066 end if
1067 end if
1068 end if
1069 end if
1070 end do
1071 end do
1072 end do
1073
1074# 443 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1075#if defined(MFC_OpenACC)
1076# 443 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1077!$acc end parallel loop
1078# 443 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1079#elif defined(MFC_OpenMP)
1080# 443 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1081
1082# 443 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1083!$omp end target teams loop
1084# 443 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1085#endif
1086
1087 end subroutine s_ib_3d_airfoil
1088
1089 !> Mark cells inside a rectangular immersed boundary
1090 subroutine s_ib_rectangle(patch_id, ib_markers, xp, yp)
1091
1092 integer, intent(in) :: patch_id
1093 type(integer_field), intent(inout) :: ib_markers
1094 integer, intent(in) :: xp, yp !< integers containing the periodicity projection information
1095 integer :: i, j, il, ir, jl, jr !< generic loop iterators
1096 integer :: encoded_patch_id
1097 real(wp) :: corner_distance !< Equation of state parameters
1098 real(wp), dimension(1:3) :: xy_local !< x and y coordinates in local IB frame
1099 real(wp), dimension(1:2) :: length, center !< x and y coordinates in local IB frame
1100 real(wp), dimension(1:3,1:3) :: inverse_rotation
1101
1102 ! Transferring the rectangle's centroid and length information
1103
1104 center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg)
1105 center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg)
1106 length(1) = patch_ib(patch_id)%length_x
1107 length(2) = patch_ib(patch_id)%length_y
1108 inverse_rotation(:,:) = patch_ib(patch_id)%rotation_matrix_inverse(:,:)
1109
1110 ! encode the periodicity information into the patch_id
1111 call s_encode_patch_periodicity(patch_id, xp, yp, 0, encoded_patch_id)
1112
1113 ! find the indices to the left and right of the IB in i, j, k
1114 il = -gp_layers - 1
1115 jl = -gp_layers - 1
1116 ir = m + gp_layers + 1
1117 jr = n + gp_layers + 1
1118 corner_distance = sqrt(dot_product(length, length))/2._wp ! maximum distance any marker can be from the center
1119 call get_bounding_indices(center(1) - corner_distance, center(1) + corner_distance, x_cc, il, ir)
1120 call get_bounding_indices(center(2) - corner_distance, center(2) + corner_distance, y_cc, jl, jr)
1121
1122 ! Assign primitive variables if rectangle covers cell and patch has write permission
1123
1124# 481 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1125
1126# 481 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1127#if defined(MFC_OpenACC)
1128# 481 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1129!$acc parallel loop collapse(2) gang vector default(present) private(i, j, xy_local) copyin(encoded_patch_id, center, length, inverse_rotation, x_cc, y_cc)
1130# 481 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1131#elif defined(MFC_OpenMP)
1132# 481 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1133
1134# 481 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1135
1136# 481 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1137
1138# 481 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1139!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(2) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, j, xy_local) map(to:encoded_patch_id, center, length, inverse_rotation, x_cc, y_cc)
1140# 481 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1141#endif
1142# 483 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1143 do j = jl, jr
1144 do i = il, ir
1145 ! get the x and y coordinates in the local IB frame
1146 xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp]
1147 xy_local = matmul(inverse_rotation, xy_local)
1148
1149 if (-0.5_wp*length(1) <= xy_local(1) .and. 0.5_wp*length(1) >= xy_local(1) .and. -0.5_wp*length(2) <= xy_local(2) &
1150 & .and. 0.5_wp*length(2) >= xy_local(2)) then
1151 ! Updating the patch identities bookkeeping variable
1152 ib_markers%sf(i, j, 0) = encoded_patch_id
1153 end if
1154 end do
1155 end do
1156
1157# 496 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1158#if defined(MFC_OpenACC)
1159# 496 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1160!$acc end parallel loop
1161# 496 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1162#elif defined(MFC_OpenMP)
1163# 496 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1164
1165# 496 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1166!$omp end target teams loop
1167# 496 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1168#endif
1169
1170 end subroutine s_ib_rectangle
1171
1172 !> Mark cells inside a spherical immersed boundary
1173 subroutine s_ib_sphere(patch_id, ib_markers, xp, yp, zp)
1174
1175 integer, intent(in) :: patch_id
1176 type(integer_field), intent(inout) :: ib_markers
1177 integer, intent(in) :: xp, yp, zp !< integers containing the periodicity projection information
1178 ! Generic loop iterators
1179 integer :: i, j, k
1180 integer :: il, ir, jl, jr, kl, kr
1181 integer :: encoded_patch_id
1182 real(wp) :: radius
1183 real(wp), dimension(1:3) :: center
1184
1185 ! Variables to initialize the pressure field that corresponds to the bubble-collapse test case found in Tiwari et al. (2013)
1186
1187 ! Transferring spherical patch's radius, centroid, smoothing patch identity and smoothing coefficient information
1188
1189 center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg)
1190 center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg)
1191 center(3) = patch_ib(patch_id)%z_centroid + real(zp, wp)*(z_domain%end - z_domain%beg)
1192 radius = patch_ib(patch_id)%radius
1193
1194 ! completely skip particles no in the domain
1195 if (center(1) - radius > x_cc(m + gp_layers + 1) .or. center(1) + radius < x_cc(-gp_layers - 1) .or. center(2) &
1196 & - radius > y_cc(n + gp_layers + 1) .or. center(2) + radius < y_cc(-gp_layers - 1) .or. center(3) - radius > z_cc(p &
1197 & + gp_layers + 1) .or. center(3) + radius < z_cc(-gp_layers - 1)) then
1198 return
1199 end if
1200
1201 ! encode the periodicity information into the patch_id
1202 call s_encode_patch_periodicity(patch_id, xp, yp, zp, encoded_patch_id)
1203
1204 ! find the indices to the left and right of the IB in i, j, k
1205 il = -gp_layers - 1
1206 jl = -gp_layers - 1
1207 kl = -gp_layers - 1
1208 ir = m + gp_layers + 1
1209 jr = n + gp_layers + 1
1210 kr = p + gp_layers + 1
1211 call get_bounding_indices(center(1) - radius, center(1) + radius, x_cc, il, ir)
1212 call get_bounding_indices(center(2) - radius, center(2) + radius, y_cc, jl, jr)
1213 call get_bounding_indices(center(3) - radius, center(3) + radius, z_cc, kl, kr)
1214
1215 ! Checking whether the sphere covers a particular cell in the domain and verifying whether the current patch has permission
1216 ! to write to that cell. If both queries check out, the primitive variables of the current patch are assigned to this cell.
1217
1218# 545 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1219
1220# 545 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1221#if defined(MFC_OpenACC)
1222# 545 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1223!$acc parallel loop collapse(3) gang vector default(present) private(i, j, k, cart_y, cart_z) copyin(encoded_patch_id, center, radius)
1224# 545 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1225#elif defined(MFC_OpenMP)
1226# 545 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1227
1228# 545 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1229
1230# 545 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1231
1232# 545 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1233!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, j, k, cart_y, cart_z) map(to:encoded_patch_id, center, radius)
1234# 545 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1235#endif
1236 do k = kl, kr
1237 do j = jl, jr
1238 do i = il, ir
1239 if (grid_geometry == 3) then
1240 call s_convert_cylindrical_to_cartesian_coord(y_cc(j), z_cc(k))
1241 else
1242 cart_y = y_cc(j)
1243 cart_z = z_cc(k)
1244 end if
1245 ! Updating the patch identities bookkeeping variable
1246 if (((x_cc(i) - center(1))**2 + (cart_y - center(2))**2 + (cart_z - center(3))**2 <= radius**2)) then
1247 ib_markers%sf(i, j, k) = encoded_patch_id
1248 end if
1249 end do
1250 end do
1251 end do
1252
1253# 562 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1254#if defined(MFC_OpenACC)
1255# 562 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1256!$acc end parallel loop
1257# 562 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1258#elif defined(MFC_OpenMP)
1259# 562 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1260
1261# 562 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1262!$omp end target teams loop
1263# 562 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1264#endif
1265
1266 end subroutine s_ib_sphere
1267
1268 !> Mark cells inside a cuboidal immersed boundary
1269 subroutine s_ib_cuboid(patch_id, ib_markers, xp, yp, zp)
1270
1271 integer, intent(in) :: patch_id
1272 type(integer_field), intent(inout) :: ib_markers
1273 integer, intent(in) :: xp, yp, zp !< integers containing the periodicity projection information
1274 integer :: i, j, k, ir, il, jr, jl, kr, kl !< Generic loop iterators
1275 integer :: encoded_patch_id
1276 real(wp), dimension(1:3) :: xyz_local, center, length !< x and y coordinates in local IB frame
1277 real(wp), dimension(1:3,1:3) :: inverse_rotation
1278 real(wp) :: corner_distance
1279
1280 ! Transferring the cuboid's centroid and length information
1281
1282 center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg)
1283 center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg)
1284 center(3) = patch_ib(patch_id)%z_centroid + real(zp, wp)*(z_domain%end - z_domain%beg)
1285 length(1) = patch_ib(patch_id)%length_x
1286 length(2) = patch_ib(patch_id)%length_y
1287 length(3) = patch_ib(patch_id)%length_z
1288 inverse_rotation(:,:) = patch_ib(patch_id)%rotation_matrix_inverse(:,:)
1289
1290 ! encode the periodicity information into the patch_id
1291 call s_encode_patch_periodicity(patch_id, xp, yp, zp, encoded_patch_id)
1292
1293 ! find the indices to the left and right of the IB in i, j, k
1294 il = -gp_layers - 1
1295 jl = -gp_layers - 1
1296 kl = -gp_layers - 1
1297 ir = m + gp_layers + 1
1298 jr = n + gp_layers + 1
1299 kr = p + gp_layers + 1
1300 corner_distance = sqrt(dot_product(length, length))/2._wp ! maximum distance any marker can be from the center
1301 call get_bounding_indices(center(1) - corner_distance, center(1) + corner_distance, x_cc, il, ir)
1302 call get_bounding_indices(center(2) - corner_distance, center(2) + corner_distance, y_cc, jl, jr)
1303 call get_bounding_indices(center(3) - corner_distance, center(3) + corner_distance, z_cc, kl, kr)
1304
1305 ! Checking whether the cuboid covers a particular cell in the domain and verifying whether the current patch has permission
1306 ! to write to to that cell. If both queries check out, the primitive variables of the current patch are assigned to this
1307 ! cell.
1308
1309# 606 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1310
1311# 606 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1312#if defined(MFC_OpenACC)
1313# 606 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1314!$acc parallel loop collapse(3) gang vector default(present) private(i, j, k, xyz_local, cart_y, cart_z) copyin(encoded_patch_id, center, length, inverse_rotation)
1315# 606 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1316#elif defined(MFC_OpenMP)
1317# 606 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1318
1319# 606 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1320
1321# 606 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1322
1323# 606 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1324!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, j, k, xyz_local, cart_y, cart_z) map(to:encoded_patch_id, center, length, inverse_rotation)
1325# 606 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1326#endif
1327# 608 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1328 do k = kl, kr
1329 do j = jl, jr
1330 do i = il, ir
1331 if (grid_geometry == 3) then
1332 ! TODO :: This does not work and is not covered by any tests. This should be fixed
1333 call s_convert_cylindrical_to_cartesian_coord(y_cc(j), z_cc(k))
1334 else
1335 cart_y = y_cc(j)
1336 cart_z = z_cc(k)
1337 end if
1338 xyz_local = [x_cc(i), cart_y, cart_z] - center ! get coordinate frame centered on IB
1339 xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinates
1340
1341 if (-0.5*length(1) <= xyz_local(1) .and. 0.5*length(1) >= xyz_local(1) .and. -0.5*length(2) <= xyz_local(2) &
1342 & .and. 0.5*length(2) >= xyz_local(2) .and. -0.5*length(3) <= xyz_local(3) .and. 0.5*length(3) &
1343 & >= xyz_local(3)) then
1344 ! Updating the patch identities bookkeeping variable
1345 ib_markers%sf(i, j, k) = encoded_patch_id
1346 end if
1347 end do
1348 end do
1349 end do
1350
1351# 630 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1352#if defined(MFC_OpenACC)
1353# 630 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1354!$acc end parallel loop
1355# 630 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1356#elif defined(MFC_OpenMP)
1357# 630 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1358
1359# 630 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1360!$omp end target teams loop
1361# 630 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1362#endif
1363
1364 end subroutine s_ib_cuboid
1365
1366 !> Mark cells inside a cylindrical immersed boundary
1367 subroutine s_ib_cylinder(patch_id, ib_markers, xp, yp, zp)
1368
1369 integer, intent(in) :: patch_id
1370 type(integer_field), intent(inout) :: ib_markers
1371 integer, intent(in) :: xp, yp, zp !< integers containing the periodicity projection information
1372 integer :: i, j, k, il, ir, jl, jr, kl, kr !< Generic loop iterators
1373 integer :: encoded_patch_id
1374 real(wp) :: radius
1375 real(wp), dimension(1:3) :: xyz_local, center, length !< x and y coordinates in local IB frame
1376 real(wp), dimension(1:3,1:3) :: inverse_rotation
1377 real(wp) :: corner_distance
1378
1379 ! Transferring the cylindrical patch's centroid, length, radius,
1380
1381 center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg)
1382 center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg)
1383 center(3) = patch_ib(patch_id)%z_centroid + real(zp, wp)*(z_domain%end - z_domain%beg)
1384 length(1) = patch_ib(patch_id)%length_x
1385 length(2) = patch_ib(patch_id)%length_y
1386 length(3) = patch_ib(patch_id)%length_z
1387 radius = patch_ib(patch_id)%radius
1388 inverse_rotation(:,:) = patch_ib(patch_id)%rotation_matrix_inverse(:,:)
1389
1390 ! encode the periodicity information into the patch_id
1391 call s_encode_patch_periodicity(patch_id, xp, yp, zp, encoded_patch_id)
1392
1393 il = -gp_layers - 1
1394 jl = -gp_layers - 1
1395 kl = -gp_layers - 1
1396 ir = m + gp_layers + 1
1397 jr = n + gp_layers + 1
1398 kr = p + gp_layers + 1
1399 corner_distance = sqrt(radius**2 + maxval(length)**2) ! distance to rim of cylinder
1400 call get_bounding_indices(center(1) - corner_distance, center(1) + corner_distance, x_cc, il, ir)
1401 call get_bounding_indices(center(2) - corner_distance, center(2) + corner_distance, y_cc, jl, jr)
1402 call get_bounding_indices(center(3) - corner_distance, center(3) + corner_distance, z_cc, kl, kr)
1403
1404 ! Checking whether the cylinder covers a particular cell in the domain and verifying whether the current patch has the
1405 ! permission to write to that cell. If both queries check out, the primitive variables of the current patch are assigned to
1406 ! this cell.
1407
1408# 675 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1409
1410# 675 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1411#if defined(MFC_OpenACC)
1412# 675 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1413!$acc parallel loop collapse(3) gang vector default(present) private(i, j, k, xyz_local, cart_y, cart_z) copyin(encoded_patch_id, center, length, radius, inverse_rotation)
1414# 675 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1415#elif defined(MFC_OpenMP)
1416# 675 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1417
1418# 675 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1419
1420# 675 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1421
1422# 675 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1423!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, j, k, xyz_local, cart_y, cart_z) map(to:encoded_patch_id, center, length, radius, inverse_rotation)
1424# 675 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1425#endif
1426# 677 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1427 do k = kl, kr
1428 do j = jl, jr
1429 do i = il, ir
1430 if (grid_geometry == 3) then
1431 call s_convert_cylindrical_to_cartesian_coord(y_cc(j), z_cc(k))
1432 else
1433 cart_y = y_cc(j)
1434 cart_z = z_cc(k)
1435 end if
1436 xyz_local = [x_cc(i), cart_y, cart_z] - center ! get coordinate frame centered on IB
1437 xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinates
1438
1439 if (((.not. f_is_default(length(1)) .and. xyz_local(2)**2 + xyz_local(3)**2 <= radius**2 .and. &
1440 & -0.5_wp*length(1) <= xyz_local(1) .and. 0.5_wp*length(1) >= xyz_local(1)) &
1441 & .or. (.not. f_is_default(length(2)) .and. xyz_local(1)**2 + xyz_local(3)**2 <= radius**2 .and. &
1442 & -0.5_wp*length(2) <= xyz_local(2) .and. 0.5_wp*length(2) >= xyz_local(2)) &
1443 & .or. (.not. f_is_default(length(3)) .and. xyz_local(1)**2 + xyz_local(2)**2 <= radius**2 .and. &
1444 & -0.5_wp*length(3) <= xyz_local(3) .and. 0.5_wp*length(3) >= xyz_local(3)))) then
1445 ! Updating the patch identities bookkeeping variable
1446 ib_markers%sf(i, j, k) = encoded_patch_id
1447 end if
1448 end do
1449 end do
1450 end do
1451
1452# 701 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1453#if defined(MFC_OpenACC)
1454# 701 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1455!$acc end parallel loop
1456# 701 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1457#elif defined(MFC_OpenMP)
1458# 701 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1459
1460# 701 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1461!$omp end target teams loop
1462# 701 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1463#endif
1464
1465 end subroutine s_ib_cylinder
1466
1467 !> Mark cells inside a 2D elliptical immersed boundary
1468 subroutine s_ib_ellipse(patch_id, ib_markers, xp, yp)
1469
1470 integer, intent(in) :: patch_id
1471 type(integer_field), intent(inout) :: ib_markers
1472 integer, intent(in) :: xp, yp !< integers containing the periodicity projection information
1473 integer :: i, j, il, ir, jl, jr !< Generic loop iterators
1474 integer :: encoded_patch_id
1475 real(wp), dimension(1:3) :: xy_local !< x and y coordinates in local IB frame
1476 real(wp), dimension(1:2) :: ellipse_coeffs !< a and b in the ellipse coefficients
1477 real(wp), dimension(1:2) :: center !< x and y coordinates in local IB frame
1478 real(wp), dimension(1:3,1:3) :: inverse_rotation
1479
1480 ! Transferring the ellipse's centroid and length information
1481
1482 center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg)
1483 center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg)
1484 ellipse_coeffs(1) = 0.5_wp*patch_ib(patch_id)%length_x
1485 ellipse_coeffs(2) = 0.5_wp*patch_ib(patch_id)%length_y
1486 inverse_rotation(:,:) = patch_ib(patch_id)%rotation_matrix_inverse(:,:)
1487
1488 ! encode the periodicity information into the patch_id
1489 call s_encode_patch_periodicity(patch_id, xp, yp, 0, encoded_patch_id)
1490
1491 ! find the indices to the left and right of the IB in i, j, k
1492 il = -gp_layers - 1
1493 jl = -gp_layers - 1
1494 ir = m + gp_layers + 1
1495 jr = n + gp_layers + 1
1496 call get_bounding_indices(center(1) - maxval(ellipse_coeffs)*2._wp, center(1) + maxval(ellipse_coeffs)*2._wp, x_cc, il, ir)
1497 call get_bounding_indices(center(2) - maxval(ellipse_coeffs)*2._wp, center(2) + maxval(ellipse_coeffs)*2._wp, y_cc, jl, jr)
1498
1499 ! Checking whether the ellipse covers a particular cell in the domain
1500
1501# 738 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1502
1503# 738 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1504#if defined(MFC_OpenACC)
1505# 738 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1506!$acc parallel loop collapse(2) gang vector default(present) private(i, j, xy_local) copyin(encoded_patch_id, center, ellipse_coeffs, inverse_rotation, x_cc, y_cc)
1507# 738 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1508#elif defined(MFC_OpenMP)
1509# 738 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1510
1511# 738 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1512
1513# 738 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1514
1515# 738 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1516!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(2) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, j, xy_local) map(to:encoded_patch_id, center, ellipse_coeffs, inverse_rotation, x_cc, y_cc)
1517# 738 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1518#endif
1519# 740 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1520 do j = jl, jr
1521 do i = il, ir
1522 ! get the x and y coordinates in the local IB frame
1523 xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp]
1524 xy_local = matmul(inverse_rotation, xy_local)
1525
1526 ! Ellipse condition (x/a)^2 + (y/b)^2 <= 1
1527 if ((xy_local(1)/ellipse_coeffs(1))**2 + (xy_local(2)/ellipse_coeffs(2))**2 <= 1._wp) then
1528 ! Updating the patch identities bookkeeping variable
1529 ib_markers%sf(i, j, 0) = encoded_patch_id
1530 end if
1531 end do
1532 end do
1533
1534# 753 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1535#if defined(MFC_OpenACC)
1536# 753 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1537!$acc end parallel loop
1538# 753 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1539#elif defined(MFC_OpenMP)
1540# 753 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1541
1542# 753 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1543!$omp end target teams loop
1544# 753 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1545#endif
1546
1547 end subroutine s_ib_ellipse
1548
1549 !> The STL patch is a 2D geometry that is imported from an STL file.
1550 subroutine s_ib_model(patch_id, ib_markers, xp, yp)
1551
1552 integer, intent(in) :: patch_id
1553 type(integer_field), intent(inout) :: ib_markers
1554 integer, intent(in) :: xp, yp !< integers containing the periodicity projection information
1555 integer :: i, j, il, ir, jl, jr !< Generic loop iterators
1556 integer :: spc, encoded_patch_id
1557 integer :: cx, cy
1558 real(wp) :: lx(2), ly(2)
1559 real(wp), dimension(1:2) :: bbox_min, bbox_max
1560 real(wp), dimension(1:3) :: local_corner, world_corner
1561 real(wp) :: eta, threshold
1562 real(wp), dimension(1:3) :: offset
1563 real(wp), dimension(1:3) :: center, xy_local
1564 real(wp), dimension(1:3,1:3) :: inverse_rotation, rotation
1565
1566 center = 0._wp
1567 center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg)
1568 center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg)
1569 inverse_rotation(:,:) = patch_ib(patch_id)%rotation_matrix_inverse(:,:)
1570 rotation(:,:) = patch_ib(patch_id)%rotation_matrix(:,:)
1571 offset(:) = patch_ib(patch_id)%centroid_offset(:)
1572 spc = patch_ib(patch_id)%model_spc
1573 threshold = patch_ib(patch_id)%model_threshold
1574
1575 ! encode the periodicity information into the patch_id
1576 call s_encode_patch_periodicity(patch_id, xp, yp, 0, encoded_patch_id)
1577
1578 il = -gp_layers - 1
1579 jl = -gp_layers - 1
1580 ir = m + gp_layers + 1
1581 jr = n + gp_layers + 1
1582
1583 ! Local-space bounding box extents (min=1, max=2 in the third index)
1584 lx(1) = stl_bounding_boxes(patch_id, 1, 1) + offset(1)
1585 lx(2) = stl_bounding_boxes(patch_id, 1, 3) + offset(1)
1586 ly(1) = stl_bounding_boxes(patch_id, 2, 1) + offset(2)
1587 ly(2) = stl_bounding_boxes(patch_id, 2, 3) + offset(2)
1588
1589 bbox_min = 1e12
1590 bbox_max = -1e12
1591 ! Enumerate all 8 corners of the local bounding box, rotate to world space, track world-space AABB
1592 do cx = 1, 2
1593 do cy = 1, 2
1594 local_corner = [lx(cx), ly(cy), 0._wp]
1595 world_corner = matmul(rotation, local_corner) + center
1596 bbox_min(1) = min(bbox_min(1), world_corner(1))
1597 bbox_min(2) = min(bbox_min(2), world_corner(2))
1598 bbox_max(1) = max(bbox_max(1), world_corner(1))
1599 bbox_max(2) = max(bbox_max(2), world_corner(2))
1600 end do
1601 end do
1602
1603 call get_bounding_indices(bbox_min(1), bbox_max(1), x_cc, il, ir)
1604 call get_bounding_indices(bbox_min(2), bbox_max(2), y_cc, jl, jr)
1605
1606
1607# 814 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1608
1609# 814 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1610#if defined(MFC_OpenACC)
1611# 814 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1612!$acc parallel loop collapse(2) gang vector default(present) private(i, j, xy_local, eta) copyin(patch_id, encoded_patch_id, center, inverse_rotation, offset, spc, threshold)
1613# 814 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1614#elif defined(MFC_OpenMP)
1615# 814 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1616
1617# 814 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1618
1619# 814 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1620
1621# 814 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1622!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(2) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, j, xy_local, eta) map(to:patch_id, encoded_patch_id, center, inverse_rotation, offset, spc, threshold)
1623# 814 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1624#endif
1625# 816 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1626 do i = il, ir
1627 do j = jl, jr
1628 xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp]
1629 xy_local = matmul(inverse_rotation, xy_local)
1630 xy_local = xy_local - offset
1631
1632 eta = f_model_is_inside_flat(gpu_ntrs(patch_id), patch_id, xy_local)
1633
1634 ! Reading STL boundary vertices and compute the levelset and levelset_norm
1635 if (eta > threshold) then
1636 ib_markers%sf(i, j, 0) = encoded_patch_id
1637 end if
1638 end do
1639 end do
1640
1641# 830 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1642#if defined(MFC_OpenACC)
1643# 830 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1644!$acc end parallel loop
1645# 830 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1646#elif defined(MFC_OpenMP)
1647# 830 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1648
1649# 830 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1650!$omp end target teams loop
1651# 830 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1652#endif
1653
1654 end subroutine s_ib_model
1655
1656 !> The STL patch is a 3D geometry that is imported from an STL file.
1657 subroutine s_ib_3d_model(patch_id, ib_markers, xp, yp, zp)
1658
1659 integer, intent(in) :: patch_id
1660 type(integer_field), intent(inout) :: ib_markers
1661 integer, intent(in) :: xp, yp, zp !< integers containing the periodicity projection information
1662 integer :: i, j, k, il, ir, jl, jr, kl, kr !< Generic loop iterators
1663 integer :: spc, encoded_patch_id
1664 real(wp) :: eta, threshold
1665 real(wp), dimension(1:3) :: offset
1666 real(wp), dimension(1:3) :: center, xyz_local
1667 real(wp), dimension(1:3,1:3) :: inverse_rotation, rotation
1668 integer :: cx, cy, cz
1669 real(wp) :: lx(2), ly(2), lz(2)
1670 real(wp), dimension(1:3) :: bbox_min, bbox_max, local_corner, world_corner
1671
1672 center = 0._wp
1673 center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg)
1674 center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg)
1675 center(3) = patch_ib(patch_id)%z_centroid + real(zp, wp)*(z_domain%end - z_domain%beg)
1676 inverse_rotation(:,:) = patch_ib(patch_id)%rotation_matrix_inverse(:,:)
1677 offset(:) = patch_ib(patch_id)%centroid_offset(:)
1678 spc = patch_ib(patch_id)%model_spc
1679 threshold = patch_ib(patch_id)%model_threshold
1680 rotation(:,:) = patch_ib(patch_id)%rotation_matrix(:,:)
1681
1682 ! encode the periodicity information into the patch_id
1683 call s_encode_patch_periodicity(patch_id, xp, yp, zp, encoded_patch_id)
1684
1685 il = -gp_layers - 1
1686 jl = -gp_layers - 1
1687 kl = -gp_layers - 1
1688 ir = m + gp_layers + 1
1689 jr = n + gp_layers + 1
1690 kr = p + gp_layers + 1
1691
1692 ! Local-space bounding box extents (min=1, max=2 in the third index)
1693 lx(1) = stl_bounding_boxes(patch_id, 1, 1) + offset(1)
1694 lx(2) = stl_bounding_boxes(patch_id, 1, 3) + offset(1)
1695 ly(1) = stl_bounding_boxes(patch_id, 2, 1) + offset(2)
1696 ly(2) = stl_bounding_boxes(patch_id, 2, 3) + offset(2)
1697 lz(1) = stl_bounding_boxes(patch_id, 3, 1) + offset(3)
1698 lz(2) = stl_bounding_boxes(patch_id, 3, 3) + offset(3)
1699
1700 bbox_min = 1e12
1701 bbox_max = -1e12
1702 ! Enumerate all 8 corners of the local bounding box, rotate to world space, track world-space AABB
1703 do cx = 1, 2
1704 do cy = 1, 2
1705 do cz = 1, 2
1706 local_corner = [lx(cx), ly(cy), lz(cz)]
1707 world_corner = matmul(rotation, local_corner) + center
1708 bbox_min(1) = min(bbox_min(1), world_corner(1))
1709 bbox_min(2) = min(bbox_min(2), world_corner(2))
1710 bbox_min(3) = min(bbox_min(3), world_corner(3))
1711 bbox_max(1) = max(bbox_max(1), world_corner(1))
1712 bbox_max(2) = max(bbox_max(2), world_corner(2))
1713 bbox_max(3) = max(bbox_max(3), world_corner(3))
1714 end do
1715 end do
1716 end do
1717
1718 call get_bounding_indices(bbox_min(1), bbox_max(1), x_cc, il, ir)
1719 call get_bounding_indices(bbox_min(2), bbox_max(2), y_cc, jl, jr)
1720 call get_bounding_indices(bbox_min(3), bbox_max(3), z_cc, kl, kr)
1721
1722
1723# 900 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1724
1725# 900 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1726#if defined(MFC_OpenACC)
1727# 900 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1728!$acc parallel loop collapse(3) gang vector default(present) private(i, j, k, xyz_local, eta) copyin(patch_id, encoded_patch_id, center, inverse_rotation, offset, spc, threshold)
1729# 900 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1730#elif defined(MFC_OpenMP)
1731# 900 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1732
1733# 900 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1734
1735# 900 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1736
1737# 900 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1738!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, j, k, xyz_local, eta) map(to:patch_id, encoded_patch_id, center, inverse_rotation, offset, spc, threshold)
1739# 900 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1740#endif
1741# 902 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1742 do i = il, ir
1743 do j = jl, jr
1744 do k = kl, kr
1745 xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), z_cc(k) - center(3)]
1746 xyz_local = matmul(inverse_rotation, xyz_local)
1747 xyz_local = xyz_local - offset
1748
1749 eta = f_model_is_inside_flat(gpu_ntrs(patch_id), patch_id, xyz_local)
1750
1751 if (eta > patch_ib(patch_id)%model_threshold) then
1752 ib_markers%sf(i, j, k) = encoded_patch_id
1753 end if
1754 end do
1755 end do
1756 end do
1757
1758# 917 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1759#if defined(MFC_OpenACC)
1760# 917 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1761!$acc end parallel loop
1762# 917 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1763#elif defined(MFC_OpenMP)
1764# 917 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1765
1766# 917 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1767!$omp end target teams loop
1768# 917 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1769#endif
1770
1771 end subroutine s_ib_3d_model
1772
1773 !> Compute a rotation matrix for converting to the rotating frame of the boundary
1774 subroutine s_update_ib_rotation_matrix(patch_id)
1775
1776 integer, intent(in) :: patch_id
1777 real(wp), dimension(3, 3, 3) :: rotation
1778 real(wp) :: angle
1779
1780 ! construct the x, y, and z rotation matrices
1781
1782 if (num_dims == 3) then
1783 ! also compute the x and y axes in 3D
1784 angle = patch_ib(patch_id)%angles(1)
1785 rotation(1, 1,:) = [1._wp, 0._wp, 0._wp]
1786 rotation(1, 2,:) = [0._wp, cos(angle), -sin(angle)]
1787 rotation(1, 3,:) = [0._wp, sin(angle), cos(angle)]
1788
1789 angle = patch_ib(patch_id)%angles(2)
1790 rotation(2, 1,:) = [cos(angle), 0._wp, sin(angle)]
1791 rotation(2, 2,:) = [0._wp, 1._wp, 0._wp]
1792 rotation(2, 3,:) = [-sin(angle), 0._wp, cos(angle)]
1793
1794 ! apply the y rotation to the x rotation
1795 patch_ib(patch_id)%rotation_matrix(:,:) = matmul(rotation(1,:,:), rotation(2,:,:))
1796 patch_ib(patch_id)%rotation_matrix_inverse(:,:) = matmul(transpose(rotation(2,:,:)), transpose(rotation(1,:,:)))
1797 end if
1798
1799 ! z component first, since it applies in 2D and 3D
1800 angle = patch_ib(patch_id)%angles(3)
1801 rotation(3, 1,:) = [cos(angle), -sin(angle), 0._wp]
1802 rotation(3, 2,:) = [sin(angle), cos(angle), 0._wp]
1803 rotation(3, 3,:) = [0._wp, 0._wp, 1._wp]
1804
1805 if (num_dims == 3) then
1806 ! apply the z rotation to the xy rotation in 3D
1807 patch_ib(patch_id)%rotation_matrix(:,:) = matmul(patch_ib(patch_id)%rotation_matrix(:,:), rotation(3,:,:))
1808 patch_ib(patch_id)%rotation_matrix_inverse(:,:) = matmul(transpose(rotation(3,:,:)), &
1809 & patch_ib(patch_id)%rotation_matrix_inverse(:,:))
1810 else
1811 ! write out only the z rotation in 2D
1812 patch_ib(patch_id)%rotation_matrix(:,:) = rotation(3,:,:)
1813 patch_ib(patch_id)%rotation_matrix_inverse(:,:) = transpose(rotation(3,:,:))
1814 end if
1815
1816 end subroutine s_update_ib_rotation_matrix
1817
1818 !> Convert cylindrical (r, theta) coordinates to Cartesian (y, z)
1819 subroutine s_convert_cylindrical_to_cartesian_coord(cyl_y, cyl_z)
1820
1821
1822# 969 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1823#if MFC_OpenACC
1824# 969 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1825!$acc routine seq
1826# 969 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1827#elif MFC_OpenMP
1828# 969 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1829
1830# 969 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1831
1832# 969 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1833!$omp declare target device_type(any)
1834# 969 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1835#endif
1836
1837 real(wp), intent(in) :: cyl_y, cyl_z
1838
1839 cart_y = cyl_y*sin(cyl_z)
1840 cart_z = cyl_y*cos(cyl_z)
1841
1842 end subroutine s_convert_cylindrical_to_cartesian_coord
1843
1844 !> Convert a 3D cylindrical coordinate vector (x, r, theta) to Cartesian (x, y, z)
1845 pure function f_convert_cyl_to_cart(cyl) result(cart)
1846
1847
1848# 981 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1849#if MFC_OpenACC
1850# 981 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1851!$acc routine seq
1852# 981 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1853#elif MFC_OpenMP
1854# 981 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1855
1856# 981 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1857
1858# 981 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1859!$omp declare target device_type(any)
1860# 981 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1861#endif
1862
1863 real(wp), dimension(1:3), intent(in) :: cyl
1864 real(wp), dimension(1:3) :: cart
1865
1866 cart = (/cyl(1), cyl(2)*sin(cyl(3)), cyl(2)*cos(cyl(3))/)
1867
1868 end function f_convert_cyl_to_cart
1869
1870 subroutine get_bounding_indices(left_bound, right_bound, cell_centers, left_index, right_index)
1871
1872 real(wp), intent(in) :: left_bound, right_bound
1873 integer, intent(inout) :: left_index, right_index
1874 real(wp), dimension(-buff_size:), intent(in) :: cell_centers
1875 integer :: itr_left, itr_middle, itr_right
1876
1877 itr_left = left_index
1878 itr_right = right_index
1879
1880 do while (itr_left + 1 < itr_right)
1881 itr_middle = (itr_left + itr_right)/2
1882 if (cell_centers(itr_middle) < left_bound) then
1883 itr_left = itr_middle
1884 else if (cell_centers(itr_middle) > left_bound) then
1885 itr_right = itr_middle
1886 else
1887 itr_left = itr_middle
1888 exit
1889 end if
1890 end do
1891 left_index = itr_left
1892
1893 itr_right = right_index
1894 do while (itr_left + 1 < itr_right)
1895 itr_middle = (itr_left + itr_right)/2
1896 if (cell_centers(itr_middle) < right_bound) then
1897 itr_left = itr_middle
1898 else if (cell_centers(itr_middle) > right_bound) then
1899 itr_right = itr_middle
1900 else
1901 itr_right = itr_middle
1902 exit
1903 end if
1904 end do
1905 right_index = itr_right
1906
1907 end subroutine get_bounding_indices
1908
1909 !> Encode the patch ID with a unique offset containing periodicity information
1910 subroutine s_encode_patch_periodicity(patch_id, x_periodicity, y_periodicity, z_periodicity, encoded_patch_id)
1911
1912 integer, intent(in) :: patch_id, x_periodicity, y_periodicity, z_periodicity
1913 integer, intent(out) :: encoded_patch_id
1914 integer :: temp_x_per, temp_y_per, temp_z_per, offset
1915
1916 encoded_patch_id = patch_id
1917
1918 temp_x_per = x_periodicity; if (x_periodicity == -1) temp_x_per = 2
1919 temp_y_per = y_periodicity; if (y_periodicity == -1) temp_y_per = 2
1920 temp_z_per = z_periodicity; if (z_periodicity == -1) temp_z_per = 2
1921
1922 offset = (num_ibs + 1)*temp_x_per + 3*(num_ibs + 1)*temp_y_per + 9*(num_ibs + 1)*temp_z_per
1923 encoded_patch_id = patch_id + offset
1924
1925 end subroutine s_encode_patch_periodicity
1926
1927 !> Decode the encoded ID to recover the original patch ID and periodicity
1928 subroutine s_decode_patch_periodicity(encoded_patch_id, patch_id, x_periodicity, y_periodicity, z_periodicity)
1929
1930
1931# 1050 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1932#if MFC_OpenACC
1933# 1050 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1934!$acc routine seq
1935# 1050 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1936#elif MFC_OpenMP
1937# 1050 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1938
1939# 1050 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1940
1941# 1050 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1942!$omp declare target device_type(any)
1943# 1050 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1944#endif
1945
1946 integer, intent(in) :: encoded_patch_id
1947 integer, intent(out) :: patch_id
1948 integer, intent(out), optional :: x_periodicity, y_periodicity, z_periodicity
1949 integer :: offset, remainder, xp, yp, zp, base
1950
1951 base = num_ibs + 1
1952
1953 patch_id = mod(encoded_patch_id - 1, base) + 1
1954 offset = (encoded_patch_id - patch_id)/base
1955
1956 xp = mod(offset, 3)
1957 remainder = offset/3
1958 yp = mod(remainder, 3)
1959 zp = remainder/3
1960
1961 ! Reverse map: 2 -> -1, 0 -> 0, 1 -> 1
1962 if (present(x_periodicity) .and. present(y_periodicity) .and. present(z_periodicity)) then
1963 x_periodicity = xp; if (xp == 2) x_periodicity = -1
1964 y_periodicity = yp; if (yp == 2) y_periodicity = -1
1965 z_periodicity = zp; if (zp == 2) z_periodicity = -1
1966 end if
1967
1968 end subroutine s_decode_patch_periodicity
1969
1970 !> Determine the periodic wrapping bounds in each direction
1971 subroutine s_get_periodicities(xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper)
1972
1973 integer, intent(out) :: xp_lower, xp_upper, yp_lower, yp_upper
1974 integer, intent(out), optional :: zp_lower, zp_upper
1975
1976 ! check domain wraps in x, y
1977
1978# 1085 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1979 ! check for periodicity
1980 if (ib_bc_x%beg == bc_periodic) then
1981 xp_lower = -1
1982 xp_upper = 1
1983 else
1984 ! if it is not periodic, then both elements are 0
1985 xp_lower = 0
1986 xp_upper = 0
1987 end if
1988# 1085 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1989 ! check for periodicity
1990 if (ib_bc_y%beg == bc_periodic) then
1991 yp_lower = -1
1992 yp_upper = 1
1993 else
1994 ! if it is not periodic, then both elements are 0
1995 yp_lower = 0
1996 yp_upper = 0
1997 end if
1998# 1095 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1999
2000 ! z only if 3D
2001 if (present(zp_lower) .and. p /= 0) then
2002 if (ib_bc_z%beg == bc_periodic) then
2003 zp_lower = -1
2004 zp_upper = 1
2005 else
2006 zp_lower = 0
2007 zp_upper = 0
2008 end if
2009 end if
2010
2011 end subroutine s_get_periodicities
2012
2013 !> Archimedes spiral function
2014 pure elemental function f_r(myth, offset, a)
2015
2016
2017# 1112 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
2018#if MFC_OpenACC
2019# 1112 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
2020!$acc routine seq
2021# 1112 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
2022#elif MFC_OpenMP
2023# 1112 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
2024
2025# 1112 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
2026
2027# 1112 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
2028!$omp declare target device_type(any)
2029# 1112 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
2030#endif
2031 real(wp), intent(in) :: myth, offset, a
2032 real(wp) :: b
2033 real(wp) :: f_r
2034
2035 ! r(th) = a + b*th
2036
2037 b = 2._wp*a/(2._wp*pi)
2038 f_r = a + b*myth + offset
2039
2040 end function f_r
2041
2042end module m_ib_patches
k
integer, intent(in) k
Definition m_phase_change.fpp.f90:1067
j
integer, intent(in) j
Definition m_phase_change.fpp.f90:1067
l
integer, intent(in) l
Definition m_phase_change.fpp.f90:1067
m_derived_types
Shared derived types for field data, patch geometry, bubble dynamics, and MPI I/O structures.
Definition m_derived_types.fpp.f90:317
m_global_parameters
Global parameters for the computational domain, fluid properties, and simulation algorithm configurat...
Definition m_global_parameters.fpp.f90:328
m_global_parameters::x_domain
type(bounds_info) x_domain
Definition m_global_parameters.fpp.f90:751
m_global_parameters::airfoil_grid_u
type(vec3_dt), dimension(:), allocatable airfoil_grid_u
Definition m_global_parameters.fpp.f90:987
m_global_parameters::y_cc
real(wp), dimension(:), allocatable, target y_cc
Definition m_global_parameters.fpp.f90:387
m_global_parameters::y_domain
type(bounds_info) y_domain
Definition m_global_parameters.fpp.f90:751
m_global_parameters::x_cc
real(wp), dimension(:), allocatable, target x_cc
Definition m_global_parameters.fpp.f90:387
m_global_parameters::p
integer p
Definition m_global_parameters.fpp.f90:352
m_global_parameters::n
integer n
Definition m_global_parameters.fpp.f90:352
m_global_parameters::np
integer np
Definition m_global_parameters.fpp.f90:988
m_global_parameters::m
integer m
Definition m_global_parameters.fpp.f90:352
m_global_parameters::dx
real(wp), dimension(:), allocatable, target dx
Definition m_global_parameters.fpp.f90:392
m_global_parameters::num_ibs
integer num_ibs
Definition m_global_parameters.fpp.f90:980
m_global_parameters::patch_ib
type(ib_patch_parameters), dimension(num_ib_patches_max) patch_ib
Immersed boundary patch parameters.
Definition m_global_parameters.fpp.f90:986
m_helper_basic
Basic floating-point utilities: approximate equality, default detection, and coordinate bounds.
Definition m_helper_basic.fpp.f90:317
m_helper
Utility routines for bubble model setup, coordinate transforms, array sampling, and special functions...
Definition m_helper.fpp.f90:328
m_ib_patches
Allocate memory and read initial condition data for IC extrusion.
Definition m_ib_patches.fpp.f90:387
m_ib_patches::s_ib_cylinder
subroutine s_ib_cylinder(patch_id, ib_markers, xp, yp, zp)
Mark cells inside a cylindrical immersed boundary.
Definition m_ib_patches.fpp.f90:1368
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:1775
m_ib_patches::s_ib_airfoil
subroutine s_ib_airfoil(patch_id, ib_markers, xp, yp)
Mark cells inside a 2D NACA 4-digit airfoil immersed boundary.
Definition m_ib_patches.fpp.f90:552
m_ib_patches::s_ib_3d_model
subroutine s_ib_3d_model(patch_id, ib_markers, xp, yp, zp)
The STL patch is a 3D geometry that is imported from an STL file.
Definition m_ib_patches.fpp.f90:1658
m_ib_patches::s_ib_circle
subroutine s_ib_circle(patch_id, ib_markers, xp, yp)
Mark cells inside a circular immersed boundary.
Definition m_ib_patches.fpp.f90:480
m_ib_patches::s_get_periodicities
subroutine 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:1972
m_ib_patches::s_ib_ellipse
subroutine s_ib_ellipse(patch_id, ib_markers, xp, yp)
Mark cells inside a 2D elliptical immersed boundary.
Definition m_ib_patches.fpp.f90:1469
m_ib_patches::get_bounding_indices
subroutine get_bounding_indices(left_bound, right_bound, cell_centers, left_index, right_index)
Definition m_ib_patches.fpp.f90:1871
m_ib_patches::cart_z
real(wp) cart_z
Definition m_ib_patches.fpp.f90:401
m_ib_patches::s_ib_3d_airfoil
subroutine s_ib_3d_airfoil(patch_id, ib_markers, xp, yp, zp)
Mark cells inside a 3D extruded NACA 4-digit airfoil immersed boundary with finite span.
Definition m_ib_patches.fpp.f90:821
m_ib_patches::s_convert_cylindrical_to_cartesian_coord
subroutine s_convert_cylindrical_to_cartesian_coord(cyl_y, cyl_z)
Convert cylindrical (r, theta) coordinates to Cartesian (y, z).
Definition m_ib_patches.fpp.f90:1820
m_ib_patches::s_ib_model
subroutine s_ib_model(patch_id, ib_markers, xp, yp)
The STL patch is a 2D geometry that is imported from an STL file.
Definition m_ib_patches.fpp.f90:1551
m_ib_patches::s_encode_patch_periodicity
subroutine 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:1911
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:419
m_ib_patches::s_ib_rectangle
subroutine s_ib_rectangle(patch_id, ib_markers, xp, yp)
Mark cells inside a rectangular immersed boundary.
Definition m_ib_patches.fpp.f90:1091
m_ib_patches::f_convert_cyl_to_cart
pure real(wp) function, dimension(1:3), public f_convert_cyl_to_cart(cyl)
Convert a 3D cylindrical coordinate vector (x, r, theta) to Cartesian (x, y, z).
Definition m_ib_patches.fpp.f90:1846
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:1929
m_ib_patches::s_ib_sphere
subroutine s_ib_sphere(patch_id, ib_markers, xp, yp, zp)
Mark cells inside a spherical immersed boundary.
Definition m_ib_patches.fpp.f90:1174
m_ib_patches::cart_y
real(wp) cart_y
Definition m_ib_patches.fpp.f90:401
m_ib_patches::f_r
pure elemental real(wp) function f_r(myth, offset, a)
Archimedes spiral function.
Definition m_ib_patches.fpp.f90:2015
m_ib_patches::s_ib_cuboid
subroutine s_ib_cuboid(patch_id, ib_markers, xp, yp, zp)
Mark cells inside a cuboidal immersed boundary.
Definition m_ib_patches.fpp.f90:1270
m_model
Binary STL file reader and processor for immersed boundary geometry.
Definition m_model.fpp.f90:318
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:1383
m_mpi_common
MPI communication layer: domain decomposition, halo exchange, reductions, and parallel I/O setup.
Definition m_mpi_common.fpp.f90:328
m_derived_types::integer_field
Derived type annexing an integer scalar field (SF).
Definition m_derived_types.fpp.f90:341