MFC
Exascale flow solver
Loading...
Searching...
No Matches
m_ib_patches.fpp.f90
Go to the documentation of this file.
1# 1 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
2!>
3!! @file
4!! @brief Contains module m_ib_patches
5
6# 1 "/home/runner/work/MFC/MFC/src/common/include/case.fpp" 1
7! This file exists so that Fypp can be run without generating case.fpp files for
8! each target. This is useful when generating documentation, for example. This
9! should also let MFC be built with CMake directly, without invoking mfc.sh.
10
11! For pre-process.
12# 8 "/home/runner/work/MFC/MFC/src/common/include/case.fpp"
13
14! For moving immersed boundaries in simulation
15# 12 "/home/runner/work/MFC/MFC/src/common/include/case.fpp"
16# 6 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp" 2
17# 1 "/home/runner/work/MFC/MFC/src/common/include/ExtrusionHardcodedIC.fpp" 1
18!> Allocate memory and read initial condition data for IC extrusion.
19!>
20!> @details
21!> This macro handles the complete initialization process for IC extrusion by:
22!>
23!> **Memory Allocation:**
24!> - stored_values(xRows, yRows, sys_size) - stores primitive variable data from files
25!> - x_coords(nrows) - stores x-coordinates from input files
26!> - y_coords(nrows) - stores y-coordinates from input files (3D case only)
27!>
28!> **File Reading Operations:**
29!> - Reads primitive variable data from multiple files with pattern:
30!> `prim.<file_number>.00.<timestep>.dat` where timestep uses `zeros_default` padding
31!> - Files are read from directory specified by `init_dir` parameter
32!> - Supports 1D, 2D, and 3D computational domains
33!>
34!> **Grid Structure Detection:**
35!> - 1D/2D: Counts lines in first file to determine xRows
36!> - 3D: Analyzes coordinate patterns to determine xRows and yRows structure
37!>
38!> **MPI Domain Mapping:**
39!> - Calculates global_offset_x and global_offset_y for MPI subdomain positioning
40!> - Maps file coordinates to local computational grid coordinates
41!>
42!> **Data Assignment:**
43!> - Populates q_prim_vf primitive variable arrays with file data
44!> - Handles momentum component indexing with special treatment for eqn_idx%mom%end
45!> - Sets eqn_idx%mom%end component to zero for 2D/3D cases
46!>
47!> **State Management:**
48!> - Uses files_loaded flag to prevent redundant file operations
49!> - Preserves data across multiple macro calls within same simulation
50!>
51!> @note File pattern uses `zeros_default` parameter (default: "000000") for timestep padding
52!> @note Directory path is hardcoded in `init_dir` parameter - modify as needed
53!> @warning Aborts execution if file reading errors occur.
54
55# 56 "/home/runner/work/MFC/MFC/src/common/include/ExtrusionHardcodedIC.fpp"
56
57# 194 "/home/runner/work/MFC/MFC/src/common/include/ExtrusionHardcodedIC.fpp"
58
59# 205 "/home/runner/work/MFC/MFC/src/common/include/ExtrusionHardcodedIC.fpp"
60# 7 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp" 2
61# 1 "/home/runner/work/MFC/MFC/src/common/include/1dHardcodedIC.fpp" 1
62# 5 "/home/runner/work/MFC/MFC/src/common/include/1dHardcodedIC.fpp"
63
64# 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# 32 "/home/runner/work/MFC/MFC/src/common/include/2dHardcodedIC.fpp"
68
69# 395 "/home/runner/work/MFC/MFC/src/common/include/2dHardcodedIC.fpp"
70# 9 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp" 2
71# 1 "/home/runner/work/MFC/MFC/src/common/include/3dHardcodedIC.fpp" 1
72# 66 "/home/runner/work/MFC/MFC/src/common/include/3dHardcodedIC.fpp"
73
74# 186 "/home/runner/work/MFC/MFC/src/common/include/3dHardcodedIC.fpp"
75# 10 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp" 2
76# 1 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 1
77# 1 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 1
78# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
79# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
80# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
81# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
82# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
83# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
84
85# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
86# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
87# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
88
89# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
90
91# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
92
93# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
94
95# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
96
97# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
98
99# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
100
101# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
102! 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, s_instantiate_stl_models, s_decode_patch_periodicity, &
399 & s_initialize_ib_airfoils
400
401contains
402
403 !> Initialize the NACA surface grids for all airfoil IB patches. Must be called after the grid is established (so dx is valid)
404 !! and before s_apply_ib_patches or s_apply_levelset.
405 subroutine s_initialize_ib_airfoils()
406
407 integer :: i, j, airfoil_id
408 integer :: np, np1, np2
409 real(wp) :: ca_in, pa, ma, ta
410 real(wp) :: xc, xa, yc, dycdxc, yt, xu, yu, xl, yl, sin_c, cos_c
411
412 do i = 1, num_ibs
413 if (patch_ib(i)%geometry /= 4 .and. patch_ib(i)%geometry /= 11) cycle
414
415 airfoil_id = patch_ib(i)%airfoil_id
416 ca_in = ib_airfoil(airfoil_id)%c
417 pa = ib_airfoil(airfoil_id)%p
418 ma = ib_airfoil(airfoil_id)%m
419 ta = ib_airfoil(airfoil_id)%t
420
421 np1 = int((pa*ca_in/dx(0))*20)
422 np2 = int(((ca_in - pa*ca_in)/dx(0))*20)
423 np = np1 + np2 + 1
424 ib_airfoil_grids(airfoil_id)%Np = np
425
426# 51 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
427#if defined(MFC_OpenACC)
428# 51 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
429!$acc update device(ib_airfoil_grids(airfoil_id)%Np)
430# 51 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
431#elif defined(MFC_OpenMP)
432# 51 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
433!$omp target update to(ib_airfoil_grids(airfoil_id)%Np)
434# 51 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
435#endif
436
437 if (.not. allocated(ib_airfoil_grids(airfoil_id)%upper)) then
438#ifdef MFC_DEBUG
439# 54 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
440 block
441# 54 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
442 use iso_fortran_env, only: output_unit
443# 54 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
444
445# 54 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
446 print *, 'm_ib_patches.fpp:54: ', '@:ALLOCATE(ib_airfoil_grids(airfoil_id)%upper(1:Np))'
447# 54 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
448
449# 54 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
450 call flush (output_unit)
451# 54 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
452 end block
453# 54 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
454#endif
455# 54 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
456 allocate (ib_airfoil_grids(airfoil_id)%upper(1:np))
457# 54 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
458
459# 54 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
460
461# 54 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
462#if defined(MFC_OpenACC)
463# 54 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
464!$acc enter data create(ib_airfoil_grids(airfoil_id)%upper)
465# 54 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
466#elif defined(MFC_OpenMP)
467# 54 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
468!$omp target enter data map(always,alloc:ib_airfoil_grids(airfoil_id)%upper)
469# 54 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
470#endif
471#ifdef MFC_DEBUG
472# 55 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
473 block
474# 55 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
475 use iso_fortran_env, only: output_unit
476# 55 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
477
478# 55 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
479 print *, 'm_ib_patches.fpp:55: ', '@:ALLOCATE(ib_airfoil_grids(airfoil_id)%lower(1:Np))'
480# 55 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
481
482# 55 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
483 call flush (output_unit)
484# 55 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
485 end block
486# 55 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
487#endif
488# 55 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
489 allocate (ib_airfoil_grids(airfoil_id)%lower(1:np))
490# 55 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
491
492# 55 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
493
494# 55 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
495#if defined(MFC_OpenACC)
496# 55 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
497!$acc enter data create(ib_airfoil_grids(airfoil_id)%lower)
498# 55 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
499#elif defined(MFC_OpenMP)
500# 55 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
501!$omp target enter data map(always,alloc:ib_airfoil_grids(airfoil_id)%lower)
502# 55 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
503#endif
504
505 ib_airfoil_grids(airfoil_id)%upper(1)%x = 0._wp
506 ib_airfoil_grids(airfoil_id)%upper(1)%y = 0._wp
507 ib_airfoil_grids(airfoil_id)%lower(1)%x = 0._wp
508 ib_airfoil_grids(airfoil_id)%lower(1)%y = 0._wp
509
510 do j = 1, np1 + np2 - 1
511 if (j <= np1) then
512 xc = j*(pa*ca_in/np1)
513 xa = xc/ca_in
514 yc = (ma/pa**2)*(2*pa*xa - xa**2)
515 dycdxc = (2*ma/pa**2)*(pa - xa)
516 else
517 xc = pa*ca_in + (j - np1)*((ca_in - pa*ca_in)/np2)
518 xa = xc/ca_in
519 yc = (ma/(1 - pa)**2)*(1 - 2*pa + 2*pa*xa - xa**2)
520 dycdxc = (2*ma/(1 - pa)**2)*(pa - xa)
521 end if
522
523 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)
524 sin_c = dycdxc/(1 + dycdxc**2)**0.5_wp
525 cos_c = 1/(1 + dycdxc**2)**0.5_wp
526
527 xu = (xa - yt*sin_c)*ca_in
528 yu = (yc + yt*cos_c)*ca_in
529 xl = (xa + yt*sin_c)*ca_in
530 yl = (yc - yt*cos_c)*ca_in
531
532 ib_airfoil_grids(airfoil_id)%upper(j + 1)%x = xu
533 ib_airfoil_grids(airfoil_id)%upper(j + 1)%y = yu
534 ib_airfoil_grids(airfoil_id)%lower(j + 1)%x = xl
535 ib_airfoil_grids(airfoil_id)%lower(j + 1)%y = yl
536 end do
537
538 ib_airfoil_grids(airfoil_id)%upper(np)%x = ca_in
539 ib_airfoil_grids(airfoil_id)%upper(np)%y = 0._wp
540 ib_airfoil_grids(airfoil_id)%lower(np)%x = ca_in
541 ib_airfoil_grids(airfoil_id)%lower(np)%y = 0._wp
542
543
544# 95 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
545#if defined(MFC_OpenACC)
546# 95 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
547!$acc update device(ib_airfoil_grids(airfoil_id)%upper, ib_airfoil_grids(airfoil_id)%lower)
548# 95 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
549#elif defined(MFC_OpenMP)
550# 95 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
551!$omp target update to(ib_airfoil_grids(airfoil_id)%upper, ib_airfoil_grids(airfoil_id)%lower)
552# 95 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
553#endif
554 end if
555 end do
556
557 end subroutine s_initialize_ib_airfoils
558
559 !> Apply all immersed boundary patch geometries to mark interior cells in the IB marker array
560 impure subroutine s_apply_ib_patches(ib_markers)
561
562 type(integer_field), intent(inout) :: ib_markers
563 integer :: i, xp, yp, zp !< iterators
564 integer :: xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper !< periodic bounds
565
566 ! 3D Patch Geometries
567
568 if (p > 0) then
569 !> IB Patches
570 !> @{
571 call s_get_periodicities(xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper)
572 do xp = xp_lower, xp_upper
573 do yp = yp_lower, yp_upper
574 do zp = zp_lower, zp_upper
575 do i = 1, num_ibs
576 if (patch_ib(i)%geometry == 8) then
577 call s_ib_sphere(i, ib_markers, xp, yp, zp)
578 else if (patch_ib(i)%geometry == 9) then
579 call s_ib_cuboid(i, ib_markers, xp, yp, zp)
580 else if (patch_ib(i)%geometry == 10) then
581 call s_ib_cylinder(i, ib_markers, xp, yp, zp)
582 else if (patch_ib(i)%geometry == 11) then
583 call s_ib_3d_airfoil(i, ib_markers, xp, yp, zp)
584 else if (patch_ib(i)%geometry == 12) then
585 call s_ib_3d_model(i, ib_markers, xp, yp, zp)
586 end if
587 end do
588 end do
589 end do
590 end do
591 !> @}
592
593 ! 2D Patch Geometries
594 else if (n > 0) then
595 !> IB Patches
596 !> @{
597 call s_get_periodicities(xp_lower, xp_upper, yp_lower, yp_upper)
598 do xp = xp_lower, xp_upper
599 do yp = yp_lower, yp_upper
600 do i = 1, num_ibs
601 if (patch_ib(i)%geometry == 2) then
602 call s_ib_circle(i, ib_markers, xp, yp)
603 else if (patch_ib(i)%geometry == 3) then
604 call s_ib_rectangle(i, ib_markers, xp, yp)
605 else if (patch_ib(i)%geometry == 4) then
606 call s_ib_airfoil(i, ib_markers, xp, yp)
607 else if (patch_ib(i)%geometry == 5) then
608 call s_ib_model(i, ib_markers, xp, yp)
609 else if (patch_ib(i)%geometry == 6) then
610 call s_ib_ellipse(i, ib_markers, xp, yp)
611 end if
612 end do
613 end do
614 end do
615 !> @}
616 end if
617
618 end subroutine s_apply_ib_patches
619
620 !> Mark cells inside a circular immersed boundary
621 subroutine s_ib_circle(patch_id, ib_markers, xp, yp)
622
623 integer, intent(in) :: patch_id
624 integer, intent(in) :: xp, yp !< integers containing the periodicity projection information
625 type(integer_field), intent(inout) :: ib_markers
626 real(wp), dimension(1:2) :: center
627 real(wp) :: radius
628 integer :: i, j, il, ir, jl, jr !< Generic loop iterators
629 integer :: encoded_patch_id
630
631 ! Transferring the circular patch's radius, centroid, smearing patch identity and smearing coefficient information
632
633 center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg)
634 center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg)
635 radius = patch_ib(patch_id)%radius
636
637 ! encode the periodicity information into the patch_id
638 call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, 0, encoded_patch_id)
639
640 ! find the indices to the left and right of the IB in i, j, k
641 il = -gp_layers - 1
642 jl = -gp_layers - 1
643 ir = m + gp_layers + 1
644 jr = n + gp_layers + 1
645 call get_bounding_indices(center(1) - radius, center(1) + radius, x_cc, il, ir)
646 call get_bounding_indices(center(2) - radius, center(2) + radius, y_cc, jl, jr)
647
648 ! Assign primitive variables if circle covers cell and patch has write permission
649
650
651# 192 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
652
653# 192 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
654#if defined(MFC_OpenACC)
655# 192 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
656!$acc parallel loop collapse(2) gang vector default(present) private(i, j) copyin(encoded_patch_id, center, radius)
657# 192 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
658#elif defined(MFC_OpenMP)
659# 192 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
660
661# 192 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
662
663# 192 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
664
665# 192 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
666!$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)
667# 192 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
668#endif
669 do j = jl, jr
670 do i = il, ir
671 if ((x_cc(i) - center(1))**2 + (y_cc(j) - center(2))**2 <= radius**2) then
672 ib_markers%sf(i, j, 0) = encoded_patch_id
673 end if
674 end do
675 end do
676
677# 200 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
678#if defined(MFC_OpenACC)
679# 200 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
680!$acc end parallel loop
681# 200 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
682#elif defined(MFC_OpenMP)
683# 200 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
684
685# 200 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
686!$omp end target teams loop
687# 200 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
688#endif
689
690 end subroutine s_ib_circle
691
692 !> Mark cells inside a 2D NACA 4-digit airfoil immersed boundary
693 subroutine s_ib_airfoil(patch_id, ib_markers, xp, yp)
694
695 integer, intent(in) :: patch_id
696 type(integer_field), intent(inout) :: ib_markers
697 integer, intent(in) :: xp, yp !< integers containing the periodicity projection information
698 real(wp) :: f, ca_in
699 integer :: i, j, k, il, ir, jl, jr
700 integer :: Np_local, airfoil_id
701 integer :: encoded_patch_id
702 real(wp), dimension(1:3) :: xy_local, offset !< x and y coordinates in local IB frame
703 real(wp), dimension(1:2) :: center !< x and y coordinates in local IB frame
704
705 airfoil_id = patch_ib(patch_id)%airfoil_id
706 center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg)
707 center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg)
708 ca_in = ib_airfoil(airfoil_id)%c
709 np_local = ib_airfoil_grids(airfoil_id)%Np
710 offset(:) = patch_ib(patch_id)%centroid_offset(:)
711
712 ! encode the periodicity information into the patch_id
713 call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, 0, encoded_patch_id)
714
715 ! find the indices to the left and right of the IB in i, j, k
716 il = -gp_layers - 1
717 jl = -gp_layers - 1
718 ir = m + gp_layers + 1
719 jr = n + gp_layers + 1
720 ! maximum distance any marker can be from the center is the length of the airfoil
721 call get_bounding_indices(center(1) - ca_in, center(1) + ca_in, x_cc, il, ir)
722 call get_bounding_indices(center(2) - ca_in, center(2) + ca_in, y_cc, jl, jr)
723
724
725# 236 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
726
727# 236 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
728#if defined(MFC_OpenACC)
729# 236 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
730!$acc parallel loop collapse(2) gang vector default(present) private(i, j, xy_local, k, f) copyin(encoded_patch_id, center, offset, ca_in, airfoil_id, Np_local, ib_airfoil_grids(airfoil_id)%upper, ib_airfoil_grids(airfoil_id)%lower)
731# 236 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
732#elif defined(MFC_OpenMP)
733# 236 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
734
735# 236 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
736
737# 236 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
738
739# 236 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
740!$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, offset, ca_in, airfoil_id, Np_local, ib_airfoil_grids(airfoil_id)%upper, ib_airfoil_grids(airfoil_id)%lower)
741# 236 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
742#endif
743# 238 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
744 do j = jl, jr
745 do i = il, ir
746 xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp] ! get coordinate frame centered on IB
747 xy_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse, xy_local) ! rotate the frame into the IB's coordinates
748 xy_local = xy_local - offset ! airfoils are a patch that require a centroid offset
749
750 if (xy_local(1) >= 0._wp .and. xy_local(1) <= ca_in) then
751 if (xy_local(2) >= 0._wp) then
752 k = 1
753 do while (ib_airfoil_grids(airfoil_id)%upper(k)%x < xy_local(1) .and. k <= np_local)
754 k = k + 1
755 end do
756 if (f_approx_equal(ib_airfoil_grids(airfoil_id)%upper(k)%x, xy_local(1))) then
757 if (xy_local(2) <= ib_airfoil_grids(airfoil_id)%upper(k)%y) then
758 ib_markers%sf(i, j, 0) = encoded_patch_id
759 end if
760 else
761 f = (ib_airfoil_grids(airfoil_id)%upper(k)%x - xy_local(1))/(ib_airfoil_grids(airfoil_id)%upper(k)%x &
762 & - ib_airfoil_grids(airfoil_id)%upper(k - 1)%x)
763 if (xy_local(2) <= ((1._wp - f)*ib_airfoil_grids(airfoil_id)%upper(k)%y &
764 & + f*ib_airfoil_grids(airfoil_id)%upper(k - 1)%y)) then
765 ib_markers%sf(i, j, 0) = encoded_patch_id
766 end if
767 end if
768 else
769 k = 1
770 do while (ib_airfoil_grids(airfoil_id)%lower(k)%x < xy_local(1))
771 k = k + 1
772 end do
773 if (f_approx_equal(ib_airfoil_grids(airfoil_id)%lower(k)%x, xy_local(1))) then
774 if (xy_local(2) >= ib_airfoil_grids(airfoil_id)%lower(k)%y) then
775 ib_markers%sf(i, j, 0) = encoded_patch_id
776 end if
777 else
778 f = (ib_airfoil_grids(airfoil_id)%lower(k)%x - xy_local(1))/(ib_airfoil_grids(airfoil_id)%lower(k)%x &
779 & - ib_airfoil_grids(airfoil_id)%lower(k - 1)%x)
780 if (xy_local(2) >= ((1._wp - f)*ib_airfoil_grids(airfoil_id)%lower(k)%y &
781 & + f*ib_airfoil_grids(airfoil_id)%lower(k - 1)%y)) then
782 ib_markers%sf(i, j, 0) = encoded_patch_id
783 end if
784 end if
785 end if
786 end if
787 end do
788 end do
789
790# 283 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
791#if defined(MFC_OpenACC)
792# 283 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
793!$acc end parallel loop
794# 283 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
795#elif defined(MFC_OpenMP)
796# 283 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
797
798# 283 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
799!$omp end target teams loop
800# 283 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
801#endif
802
803 end subroutine s_ib_airfoil
804
805 !> Mark cells inside a 3D extruded NACA 4-digit airfoil immersed boundary with finite span
806 subroutine s_ib_3d_airfoil(patch_id, ib_markers, xp, yp, zp)
807
808 integer, intent(in) :: patch_id
809 type(integer_field), intent(inout) :: ib_markers
810 integer, intent(in) :: xp, yp, zp !< integers containing the periodicity projection information
811 real(wp) :: lz, z_max, z_min, f, ca_in
812 integer :: i, j, k, l, il, ir, jl, jr, ll, lr
813 integer :: airfoil_id
814 integer :: encoded_patch_id
815 real(wp), dimension(1:3) :: xyz_local, center, offset !< x, y, z coordinates in local IB frame
816
817 airfoil_id = patch_ib(patch_id)%airfoil_id
818 center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg)
819 center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg)
820 center(3) = patch_ib(patch_id)%z_centroid + real(zp, wp)*(z_domain%end - z_domain%beg)
821 lz = patch_ib(patch_id)%length_z
822 ca_in = ib_airfoil(airfoil_id)%c
823 offset(:) = patch_ib(patch_id)%centroid_offset(:)
824
825 z_max = lz/2
826 z_min = -lz/2
827
828 ! encode the periodicity information into the patch_id
829 call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, zp, encoded_patch_id)
830
831 ! find the indices to the left and right of the IB in i, j, k
832 il = -gp_layers - 1
833 jl = -gp_layers - 1
834 ll = -gp_layers - 1
835 ir = m + gp_layers + 1
836 jr = n + gp_layers + 1
837 lr = p + gp_layers + 1
838 ! maximum distance any marker can be from the center is the length of the airfoil
839 call get_bounding_indices(center(1) - ca_in, center(1) + ca_in, x_cc, il, ir)
840 call get_bounding_indices(center(2) - ca_in, center(2) + ca_in, y_cc, jl, jr)
841 call get_bounding_indices(center(3) - ca_in, center(3) + ca_in, z_cc, ll, lr)
842
843
844# 325 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
845
846# 325 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
847#if defined(MFC_OpenACC)
848# 325 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
849!$acc parallel loop collapse(3) gang vector default(present) private(i, j, l, xyz_local, k, f) copyin(encoded_patch_id, center, offset, ca_in, airfoil_id, ib_airfoil_grids(airfoil_id)%upper, ib_airfoil_grids(airfoil_id)%lower, z_min, z_max)
850# 325 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
851#elif defined(MFC_OpenMP)
852# 325 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
853
854# 325 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
855
856# 325 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
857
858# 325 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
859!$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, offset, ca_in, airfoil_id, ib_airfoil_grids(airfoil_id)%upper, ib_airfoil_grids(airfoil_id)%lower, z_min, z_max)
860# 325 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
861#endif
862# 327 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
863 do l = ll, lr
864 do j = jl, jr
865 do i = il, ir
866 ! get coordinate frame centered on IB
867 xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), z_cc(l) - center(3)]
868 ! rotate the frame into the IB's coordinates
869 xyz_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse, xyz_local)
870 xyz_local = xyz_local - offset ! airfoils are a patch that require a centroid offset
871
872 if (xyz_local(3) >= z_min .and. xyz_local(3) <= z_max) then
873 if (xyz_local(1) >= 0._wp .and. xyz_local(1) <= ca_in) then
874 if (xyz_local(2) >= 0._wp) then
875 k = 1
876 do while (ib_airfoil_grids(airfoil_id)%upper(k)%x < xyz_local(1))
877 k = k + 1
878 end do
879 if (f_approx_equal(ib_airfoil_grids(airfoil_id)%upper(k)%x, xyz_local(1))) then
880 if (xyz_local(2) <= ib_airfoil_grids(airfoil_id)%upper(k)%y) then
881 ! IB
882 ib_markers%sf(i, j, l) = encoded_patch_id
883 end if
884 else
885 f = (ib_airfoil_grids(airfoil_id)%upper(k)%x - xyz_local(1)) &
886 & /(ib_airfoil_grids(airfoil_id)%upper(k)%x - ib_airfoil_grids(airfoil_id)%upper(k - 1)%x)
887 if (xyz_local(2) <= ((1._wp - f)*ib_airfoil_grids(airfoil_id)%upper(k)%y &
888 & + f*ib_airfoil_grids(airfoil_id)%upper(k - 1)%y)) then
889 ib_markers%sf(i, j, l) = encoded_patch_id
890 end if
891 end if
892 else
893 k = 1
894 do while (ib_airfoil_grids(airfoil_id)%lower(k)%x < xyz_local(1))
895 k = k + 1
896 end do
897 if (f_approx_equal(ib_airfoil_grids(airfoil_id)%lower(k)%x, xyz_local(1))) then
898 if (xyz_local(2) >= ib_airfoil_grids(airfoil_id)%lower(k)%y) then
899 ib_markers%sf(i, j, l) = encoded_patch_id
900 end if
901 else
902 f = (ib_airfoil_grids(airfoil_id)%lower(k)%x - xyz_local(1)) &
903 & /(ib_airfoil_grids(airfoil_id)%lower(k)%x - ib_airfoil_grids(airfoil_id)%lower(k - 1)%x)
904 if (xyz_local(2) >= ((1._wp - f)*ib_airfoil_grids(airfoil_id)%lower(k)%y &
905 & + f*ib_airfoil_grids(airfoil_id)%lower(k - 1)%y)) then
906 ib_markers%sf(i, j, l) = encoded_patch_id
907 end if
908 end if
909 end if
910 end if
911 end if
912 end do
913 end do
914 end do
915
916# 379 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
917#if defined(MFC_OpenACC)
918# 379 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
919!$acc end parallel loop
920# 379 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
921#elif defined(MFC_OpenMP)
922# 379 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
923
924# 379 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
925!$omp end target teams loop
926# 379 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
927#endif
928
929 end subroutine s_ib_3d_airfoil
930
931 !> Mark cells inside a rectangular immersed boundary
932 subroutine s_ib_rectangle(patch_id, ib_markers, xp, yp)
933
934 integer, intent(in) :: patch_id
935 type(integer_field), intent(inout) :: ib_markers
936 integer, intent(in) :: xp, yp !< integers containing the periodicity projection information
937 integer :: i, j, il, ir, jl, jr !< generic loop iterators
938 integer :: encoded_patch_id
939 real(wp) :: corner_distance !< Equation of state parameters
940 real(wp), dimension(1:3) :: xy_local !< x and y coordinates in local IB frame
941 real(wp), dimension(1:2) :: length, center !< x and y coordinates in local IB frame
942
943 ! Transferring the rectangle's centroid and length information
944
945 center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg)
946 center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg)
947
948 ! encode the periodicity information into the patch_id
949 call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, 0, encoded_patch_id)
950
951 ! find the indices to the left and right of the IB in i, j, k
952 il = -gp_layers - 1
953 jl = -gp_layers - 1
954 ir = m + gp_layers + 1
955 jr = n + gp_layers + 1
956 ! maximum distance any marker can be from the center
957 corner_distance = 0.5_wp*sqrt(patch_ib(patch_id)%length_x**2 + patch_ib(patch_id)%length_y**2)
958 call get_bounding_indices(center(1) - corner_distance, center(1) + corner_distance, x_cc, il, ir)
959 call get_bounding_indices(center(2) - corner_distance, center(2) + corner_distance, y_cc, jl, jr)
960
961 ! Assign primitive variables if rectangle covers cell and patch has write permission
962
963# 414 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
964
965# 414 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
966#if defined(MFC_OpenACC)
967# 414 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
968!$acc parallel loop collapse(2) gang vector default(present) private(i, j, xy_local) copyin(encoded_patch_id, center)
969# 414 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
970#elif defined(MFC_OpenMP)
971# 414 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
972
973# 414 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
974
975# 414 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
976
977# 414 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
978!$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)
979# 414 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
980#endif
981 do j = jl, jr
982 do i = il, ir
983 ! get the x and y coordinates in the local IB frame
984 xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp]
985 xy_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse, xy_local)
986
987 if (-0.5_wp*patch_ib(patch_id)%length_x <= xy_local(1) .and. 0.5_wp*patch_ib(patch_id)%length_x >= xy_local(1) &
988 & .and. -0.5_wp*patch_ib(patch_id)%length_y <= xy_local(2) &
989 & .and. 0.5_wp*patch_ib(patch_id)%length_y >= xy_local(2)) then
990 ! Updating the patch identities bookkeeping variable
991 ib_markers%sf(i, j, 0) = encoded_patch_id
992 end if
993 end do
994 end do
995
996# 429 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
997#if defined(MFC_OpenACC)
998# 429 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
999!$acc end parallel loop
1000# 429 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1001#elif defined(MFC_OpenMP)
1002# 429 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1003
1004# 429 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1005!$omp end target teams loop
1006# 429 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1007#endif
1008
1009 end subroutine s_ib_rectangle
1010
1011 !> Mark cells inside a spherical immersed boundary
1012 subroutine s_ib_sphere(patch_id, ib_markers, xp, yp, zp)
1013
1014 integer, intent(in) :: patch_id
1015 type(integer_field), intent(inout) :: ib_markers
1016 integer, intent(in) :: xp, yp, zp !< integers containing the periodicity projection information
1017 ! Generic loop iterators
1018 integer :: i, j, k
1019 integer :: il, ir, jl, jr, kl, kr
1020 integer :: encoded_patch_id
1021 real(wp) :: radius
1022 real(wp), dimension(1:3) :: center
1023
1024 center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg)
1025 center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg)
1026 center(3) = patch_ib(patch_id)%z_centroid + real(zp, wp)*(z_domain%end - z_domain%beg)
1027 radius = patch_ib(patch_id)%radius
1028
1029 ! completely skip particles not in the domain
1030 if (center(1) - radius > x_cc(m + gp_layers + 1) .or. center(1) + radius < x_cc(-gp_layers - 1) .or. center(2) &
1031 & - radius > y_cc(n + gp_layers + 1) .or. center(2) + radius < y_cc(-gp_layers - 1) .or. center(3) - radius > z_cc(p &
1032 & + gp_layers + 1) .or. center(3) + radius < z_cc(-gp_layers - 1)) then
1033 return
1034 end if
1035
1036 ! encode the periodicity information into the patch_id
1037 call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, zp, encoded_patch_id)
1038
1039 ! find the indices to the left and right of the IB in i, j, k
1040 il = -gp_layers - 1
1041 jl = -gp_layers - 1
1042 kl = -gp_layers - 1
1043 ir = m + gp_layers + 1
1044 jr = n + gp_layers + 1
1045 kr = p + gp_layers + 1
1046 call get_bounding_indices(center(1) - radius, center(1) + radius, x_cc, il, ir)
1047 call get_bounding_indices(center(2) - radius, center(2) + radius, y_cc, jl, jr)
1048 call get_bounding_indices(center(3) - radius, center(3) + radius, z_cc, kl, kr)
1049
1050 ! Checking whether the sphere covers a particular cell in the domain and verifying whether the current patch has permission
1051 ! to write to that cell. If both queries check out, the primitive variables of the current patch are assigned to this cell.
1052
1053# 474 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1054
1055# 474 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1056#if defined(MFC_OpenACC)
1057# 474 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1058!$acc parallel loop collapse(3) gang vector default(present) private(i, j, k) copyin(encoded_patch_id, center, radius)
1059# 474 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1060#elif defined(MFC_OpenMP)
1061# 474 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1062
1063# 474 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1064
1065# 474 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1066
1067# 474 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1068!$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) map(to:encoded_patch_id, center, radius)
1069# 474 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1070#endif
1071 do k = kl, kr
1072 do j = jl, jr
1073 do i = il, ir
1074 ! Updating the patch identities bookkeeping variable
1075 if (((x_cc(i) - center(1))**2 + (y_cc(j) - center(2))**2 + (z_cc(k) - center(3))**2 <= radius**2)) then
1076 ib_markers%sf(i, j, k) = encoded_patch_id
1077 end if
1078 end do
1079 end do
1080 end do
1081
1082# 485 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1083#if defined(MFC_OpenACC)
1084# 485 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1085!$acc end parallel loop
1086# 485 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1087#elif defined(MFC_OpenMP)
1088# 485 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1089
1090# 485 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1091!$omp end target teams loop
1092# 485 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1093#endif
1094
1095 end subroutine s_ib_sphere
1096
1097 !> Mark cells inside a cuboidal immersed boundary
1098 subroutine s_ib_cuboid(patch_id, ib_markers, xp, yp, zp)
1099
1100 integer, intent(in) :: patch_id
1101 type(integer_field), intent(inout) :: ib_markers
1102 integer, intent(in) :: xp, yp, zp !< integers containing the periodicity projection information
1103 integer :: i, j, k, ir, il, jr, jl, kr, kl !< Generic loop iterators
1104 integer :: encoded_patch_id
1105 real(wp), dimension(1:3) :: xyz_local, center !< x and y coordinates in local IB frame
1106 real(wp) :: corner_distance
1107
1108 ! Transferring the cuboid's centroid
1109
1110 center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg)
1111 center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg)
1112 center(3) = patch_ib(patch_id)%z_centroid + real(zp, wp)*(z_domain%end - z_domain%beg)
1113
1114 ! encode the periodicity information into the patch_id
1115 call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, zp, encoded_patch_id)
1116
1117 ! find the indices to the left and right of the IB in i, j, k
1118 il = -gp_layers - 1
1119 jl = -gp_layers - 1
1120 kl = -gp_layers - 1
1121 ir = m + gp_layers + 1
1122 jr = n + gp_layers + 1
1123 kr = p + gp_layers + 1
1124 corner_distance = 0.5_wp*sqrt(patch_ib(patch_id)%length_x**2 + patch_ib(patch_id)%length_y**2 &
1125 & + patch_ib(patch_id)%length_z**2) ! maximum distance any marker can be from the center
1126 call get_bounding_indices(center(1) - corner_distance, center(1) + corner_distance, x_cc, il, ir)
1127 call get_bounding_indices(center(2) - corner_distance, center(2) + corner_distance, y_cc, jl, jr)
1128 call get_bounding_indices(center(3) - corner_distance, center(3) + corner_distance, z_cc, kl, kr)
1129
1130 ! Checking whether the cuboid covers a particular cell in the domain and verifying whether the current patch has permission
1131 ! to write to to that cell. If both queries check out, the primitive variables of the current patch are assigned to this
1132 ! cell.
1133
1134# 525 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1135
1136# 525 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1137#if defined(MFC_OpenACC)
1138# 525 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1139!$acc parallel loop collapse(3) gang vector default(present) private(i, j, k, xyz_local) copyin(encoded_patch_id, center)
1140# 525 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1141#elif defined(MFC_OpenMP)
1142# 525 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1143
1144# 525 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1145
1146# 525 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1147
1148# 525 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1149!$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) map(to:encoded_patch_id, center)
1150# 525 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1151#endif
1152 do k = kl, kr
1153 do j = jl, jr
1154 do i = il, ir
1155 xyz_local = [x_cc(i), y_cc(j), z_cc(k)] - center ! get coordinate frame centered on IB
1156 ! rotate the frame into the IB's coordinates
1157 xyz_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse, xyz_local)
1158
1159 if (-0.5_wp*patch_ib(patch_id)%length_x <= xyz_local(1) &
1160 & .and. 0.5_wp*patch_ib(patch_id)%length_x >= xyz_local(1) .and. &
1161 & -0.5_wp*patch_ib(patch_id)%length_y <= xyz_local(2) &
1162 & .and. 0.5_wp*patch_ib(patch_id)%length_y >= xyz_local(2) .and. &
1163 & -0.5_wp*patch_ib(patch_id)%length_z <= xyz_local(3) &
1164 & .and. 0.5_wp*patch_ib(patch_id)%length_z >= xyz_local(3)) then
1165 ! Updating the patch identities bookkeeping variable
1166 ib_markers%sf(i, j, k) = encoded_patch_id
1167 end if
1168 end do
1169 end do
1170 end do
1171
1172# 545 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1173#if defined(MFC_OpenACC)
1174# 545 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1175!$acc end parallel loop
1176# 545 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1177#elif defined(MFC_OpenMP)
1178# 545 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1179
1180# 545 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1181!$omp end target teams loop
1182# 545 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1183#endif
1184
1185 end subroutine s_ib_cuboid
1186
1187 !> Mark cells inside a cylindrical immersed boundary
1188 subroutine s_ib_cylinder(patch_id, ib_markers, xp, yp, zp)
1189
1190 integer, intent(in) :: patch_id
1191 type(integer_field), intent(inout) :: ib_markers
1192 integer, intent(in) :: xp, yp, zp !< integers containing the periodicity projection information
1193 integer :: i, j, k, il, ir, jl, jr, kl, kr !< Generic loop iterators
1194 integer :: encoded_patch_id
1195 real(wp), dimension(1:3) :: xyz_local, center, length !< x and y coordinates in local IB frame
1196 real(wp), dimension(1:3,1:3) :: inverse_rotation
1197 real(wp) :: corner_distance
1198
1199 ! Transferring the cylindrical patch's centroid
1200
1201 center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg)
1202 center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg)
1203 center(3) = patch_ib(patch_id)%z_centroid + real(zp, wp)*(z_domain%end - z_domain%beg)
1204 length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, patch_ib(patch_id)%length_z]
1205
1206 ! encode the periodicity information into the patch_id
1207 call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, zp, encoded_patch_id)
1208
1209 il = -gp_layers - 1
1210 jl = -gp_layers - 1
1211 kl = -gp_layers - 1
1212 ir = m + gp_layers + 1
1213 jr = n + gp_layers + 1
1214 kr = p + gp_layers + 1
1215 corner_distance = sqrt(patch_ib(patch_id)%radius**2 + maxval(length)**2) ! distance to rim of cylinder
1216 call get_bounding_indices(center(1) - corner_distance, center(1) + corner_distance, x_cc, il, ir)
1217 call get_bounding_indices(center(2) - corner_distance, center(2) + corner_distance, y_cc, jl, jr)
1218 call get_bounding_indices(center(3) - corner_distance, center(3) + corner_distance, z_cc, kl, kr)
1219
1220 ! Checking whether the cylinder covers a particular cell in the domain and verifying whether the current patch has the
1221 ! permission to write to that cell. If both queries check out, the primitive variables of the current patch are assigned to
1222 ! this cell.
1223
1224# 585 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1225
1226# 585 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1227#if defined(MFC_OpenACC)
1228# 585 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1229!$acc parallel loop collapse(3) gang vector default(present) private(i, j, k, xyz_local) copyin(encoded_patch_id, center)
1230# 585 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1231#elif defined(MFC_OpenMP)
1232# 585 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1233
1234# 585 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1235
1236# 585 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1237
1238# 585 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1239!$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) map(to:encoded_patch_id, center)
1240# 585 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1241#endif
1242 do k = kl, kr
1243 do j = jl, jr
1244 do i = il, ir
1245 xyz_local = [x_cc(i), y_cc(j), z_cc(k)] - center ! get coordinate frame centered on IB
1246 ! rotate the frame into the IB's coordinates
1247 xyz_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse, xyz_local)
1248
1249 if (((.not. f_is_default(patch_ib(patch_id)%length_x) .and. xyz_local(2)**2 + xyz_local(3) &
1250 & **2 <= patch_ib(patch_id)%radius**2 .and. -0.5_wp*patch_ib(patch_id)%length_x <= xyz_local(1) &
1251 & .and. 0.5_wp*patch_ib(patch_id)%length_x >= xyz_local(1)) &
1252 & .or. (.not. f_is_default(patch_ib(patch_id)%length_y) .and. xyz_local(1)**2 + xyz_local(3) &
1253 & **2 <= patch_ib(patch_id)%radius**2 .and. -0.5_wp*patch_ib(patch_id)%length_y <= xyz_local(2) &
1254 & .and. 0.5_wp*patch_ib(patch_id)%length_y >= xyz_local(2)) &
1255 & .or. (.not. f_is_default(patch_ib(patch_id)%length_z) .and. xyz_local(1)**2 + xyz_local(2) &
1256 & **2 <= patch_ib(patch_id)%radius**2 .and. -0.5_wp*patch_ib(patch_id)%length_z <= xyz_local(3) &
1257 & .and. 0.5_wp*patch_ib(patch_id)%length_z >= xyz_local(3)))) then
1258 ! Updating the patch identities bookkeeping variable
1259 ib_markers%sf(i, j, k) = encoded_patch_id
1260 end if
1261 end do
1262 end do
1263 end do
1264
1265# 608 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1266#if defined(MFC_OpenACC)
1267# 608 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1268!$acc end parallel loop
1269# 608 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1270#elif defined(MFC_OpenMP)
1271# 608 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1272
1273# 608 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1274!$omp end target teams loop
1275# 608 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1276#endif
1277
1278 end subroutine s_ib_cylinder
1279
1280 !> Mark cells inside a 2D elliptical immersed boundary
1281 subroutine s_ib_ellipse(patch_id, ib_markers, xp, yp)
1282
1283 integer, intent(in) :: patch_id
1284 type(integer_field), intent(inout) :: ib_markers
1285 integer, intent(in) :: xp, yp !< integers containing the periodicity projection information
1286 integer :: i, j, il, ir, jl, jr !< Generic loop iterators
1287 integer :: encoded_patch_id
1288 real(wp), dimension(1:3) :: xy_local !< x and y coordinates in local IB frame
1289 real(wp), dimension(1:2) :: center !< x and y coordinates in local IB frame
1290 real(wp) :: bounding_radius
1291
1292 ! Transferring the ellipse's centroid
1293
1294 center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg)
1295 center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg)
1296
1297 ! encode the periodicity information into the patch_id
1298 call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, 0, encoded_patch_id)
1299
1300 ! find the indices to the left and right of the IB in i, j, k
1301 bounding_radius = 0.5_wp*max(patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y)
1302 il = -gp_layers - 1
1303 jl = -gp_layers - 1
1304 ir = m + gp_layers + 1
1305 jr = n + gp_layers + 1
1306 call get_bounding_indices(center(1) - bounding_radius*2._wp, center(1) + bounding_radius*2._wp, x_cc, il, ir)
1307 call get_bounding_indices(center(2) - bounding_radius*2._wp, center(2) + bounding_radius*2._wp, y_cc, jl, jr)
1308
1309 ! Checking whether the ellipse covers a particular cell in the domain
1310
1311# 642 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1312
1313# 642 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1314#if defined(MFC_OpenACC)
1315# 642 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1316!$acc parallel loop collapse(2) gang vector default(present) private(i, j, xy_local) copyin(encoded_patch_id, center)
1317# 642 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1318#elif defined(MFC_OpenMP)
1319# 642 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1320
1321# 642 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1322
1323# 642 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1324
1325# 642 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1326!$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)
1327# 642 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1328#endif
1329 do j = jl, jr
1330 do i = il, ir
1331 ! get the x and y coordinates in the local IB frame
1332 xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp]
1333 xy_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse(:,:), xy_local)
1334
1335 ! Ellipse condition (x/a)^2 + (y/b)^2 <= 1
1336 if ((xy_local(1)/(0.5_wp*patch_ib(patch_id)%length_x))**2 + (xy_local(2)/(0.5_wp*patch_ib(patch_id)%length_y)) &
1337 & **2 <= 1._wp) then
1338 ! Updating the patch identities bookkeeping variable
1339 ib_markers%sf(i, j, 0) = encoded_patch_id
1340 end if
1341 end do
1342 end do
1343
1344# 657 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1345#if defined(MFC_OpenACC)
1346# 657 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1347!$acc end parallel loop
1348# 657 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1349#elif defined(MFC_OpenMP)
1350# 657 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1351
1352# 657 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1353!$omp end target teams loop
1354# 657 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1355#endif
1356
1357 end subroutine s_ib_ellipse
1358
1359 !> The STL patch is a 2D geometry that is imported from an STL file.
1360 subroutine s_ib_model(patch_id, ib_markers, xp, yp)
1361
1362 integer, intent(in) :: patch_id
1363 type(integer_field), intent(inout) :: ib_markers
1364 integer, intent(in) :: xp, yp !< integers containing the periodicity projection information
1365 integer :: i, j, il, ir, jl, jr !< Generic loop iterators
1366 integer :: model_id, encoded_patch_id
1367 integer :: cx, cy
1368 real(wp) :: lx(2), ly(2)
1369 real(wp), dimension(1:2) :: bbox_min, bbox_max
1370 real(wp), dimension(1:3) :: local_corner, world_corner
1371 real(wp) :: eta, threshold
1372 real(wp), dimension(1:3) :: offset
1373 real(wp), dimension(1:3) :: center, xy_local
1374 real(wp), dimension(1:3,1:3) :: inverse_rotation, rotation
1375
1376 center = 0._wp
1377 center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg)
1378 center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg)
1379 inverse_rotation(:,:) = patch_ib(patch_id)%rotation_matrix_inverse(:,:)
1380 rotation(:,:) = patch_ib(patch_id)%rotation_matrix(:,:)
1381 offset(:) = patch_ib(patch_id)%centroid_offset(:)
1382 model_id = patch_ib(patch_id)%model_id
1383 threshold = stl_models(model_id)%model_threshold
1384
1385 ! encode the periodicity information into the patch_id
1386 call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, 0, encoded_patch_id)
1387
1388 il = -gp_layers - 1
1389 jl = -gp_layers - 1
1390 ir = m + gp_layers + 1
1391 jr = n + gp_layers + 1
1392
1393 ! Local-space bounding box extents (min=1, max=2 in the third index)
1394 lx(1) = stl_bounding_boxes(model_id, 1, 1) + offset(1)
1395 lx(2) = stl_bounding_boxes(model_id, 1, 3) + offset(1)
1396 ly(1) = stl_bounding_boxes(model_id, 2, 1) + offset(2)
1397 ly(2) = stl_bounding_boxes(model_id, 2, 3) + offset(2)
1398
1399 bbox_min = 1e12
1400 bbox_max = -1e12
1401 ! Enumerate all 8 corners of the local bounding box, rotate to world space, track world-space AABB
1402 do cx = 1, 2
1403 do cy = 1, 2
1404 local_corner = [lx(cx), ly(cy), 0._wp]
1405 world_corner = matmul(rotation, local_corner) + center
1406 bbox_min(1) = min(bbox_min(1), world_corner(1))
1407 bbox_min(2) = min(bbox_min(2), world_corner(2))
1408 bbox_max(1) = max(bbox_max(1), world_corner(1))
1409 bbox_max(2) = max(bbox_max(2), world_corner(2))
1410 end do
1411 end do
1412
1413 call get_bounding_indices(bbox_min(1), bbox_max(1), x_cc, il, ir)
1414 call get_bounding_indices(bbox_min(2), bbox_max(2), y_cc, jl, jr)
1415
1416
1417# 718 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1418
1419# 718 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1420#if defined(MFC_OpenACC)
1421# 718 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1422!$acc parallel loop collapse(2) gang vector default(present) private(i, j, xy_local, eta) copyin(patch_id, model_id, encoded_patch_id, center, inverse_rotation, offset, threshold)
1423# 718 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1424#elif defined(MFC_OpenMP)
1425# 718 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1426
1427# 718 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1428
1429# 718 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1430
1431# 718 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1432!$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, model_id, encoded_patch_id, center, inverse_rotation, offset, threshold)
1433# 718 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1434#endif
1435# 720 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1436 do i = il, ir
1437 do j = jl, jr
1438 xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp]
1439 xy_local = matmul(inverse_rotation, xy_local)
1440 xy_local = xy_local - offset
1441
1442 eta = f_model_is_inside_flat(gpu_ntrs(model_id), model_id, xy_local)
1443
1444 ! Reading STL boundary vertices and compute the levelset and levelset_norm
1445 if (eta > threshold) then
1446 ib_markers%sf(i, j, 0) = encoded_patch_id
1447 end if
1448 end do
1449 end do
1450
1451# 734 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1452#if defined(MFC_OpenACC)
1453# 734 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1454!$acc end parallel loop
1455# 734 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1456#elif defined(MFC_OpenMP)
1457# 734 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1458
1459# 734 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1460!$omp end target teams loop
1461# 734 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1462#endif
1463
1464 end subroutine s_ib_model
1465
1466 !> The STL patch is a 3D geometry that is imported from an STL file.
1467 subroutine s_ib_3d_model(patch_id, ib_markers, xp, yp, zp)
1468
1469 integer, intent(in) :: patch_id
1470 type(integer_field), intent(inout) :: ib_markers
1471 integer, intent(in) :: xp, yp, zp !< integers containing the periodicity projection information
1472 integer :: i, j, k, il, ir, jl, jr, kl, kr !< Generic loop iterators
1473 integer :: model_id, encoded_patch_id
1474 real(wp) :: eta, threshold
1475 real(wp), dimension(1:3) :: offset
1476 real(wp), dimension(1:3) :: center, xyz_local
1477 real(wp), dimension(1:3,1:3) :: inverse_rotation, rotation
1478 integer :: cx, cy, cz
1479 real(wp) :: lx(2), ly(2), lz(2)
1480 real(wp), dimension(1:3) :: bbox_min, bbox_max, local_corner, world_corner
1481
1482 center = 0._wp
1483 center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg)
1484 center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg)
1485 center(3) = patch_ib(patch_id)%z_centroid + real(zp, wp)*(z_domain%end - z_domain%beg)
1486 inverse_rotation(:,:) = patch_ib(patch_id)%rotation_matrix_inverse(:,:)
1487 offset(:) = patch_ib(patch_id)%centroid_offset(:)
1488 model_id = patch_ib(patch_id)%model_id
1489 threshold = stl_models(model_id)%model_threshold
1490 rotation(:,:) = patch_ib(patch_id)%rotation_matrix(:,:)
1491
1492 ! encode the periodicity information into the patch_id
1493 call s_encode_patch_periodicity(patch_ib(patch_id)%gbl_patch_id, xp, yp, zp, encoded_patch_id)
1494
1495 il = -gp_layers - 1
1496 jl = -gp_layers - 1
1497 kl = -gp_layers - 1
1498 ir = m + gp_layers + 1
1499 jr = n + gp_layers + 1
1500 kr = p + gp_layers + 1
1501
1502 ! Local-space bounding box extents (min=1, max=2 in the third index)
1503 lx(1) = stl_bounding_boxes(model_id, 1, 1) + offset(1)
1504 lx(2) = stl_bounding_boxes(model_id, 1, 3) + offset(1)
1505 ly(1) = stl_bounding_boxes(model_id, 2, 1) + offset(2)
1506 ly(2) = stl_bounding_boxes(model_id, 2, 3) + offset(2)
1507 lz(1) = stl_bounding_boxes(model_id, 3, 1) + offset(3)
1508 lz(2) = stl_bounding_boxes(model_id, 3, 3) + offset(3)
1509
1510 bbox_min = 1e12
1511 bbox_max = -1e12
1512 ! Enumerate all 8 corners of the local bounding box, rotate to world space, track world-space AABB
1513 do cx = 1, 2
1514 do cy = 1, 2
1515 do cz = 1, 2
1516 local_corner = [lx(cx), ly(cy), lz(cz)]
1517 world_corner = matmul(rotation, local_corner) + center
1518 bbox_min(1) = min(bbox_min(1), world_corner(1))
1519 bbox_min(2) = min(bbox_min(2), world_corner(2))
1520 bbox_min(3) = min(bbox_min(3), world_corner(3))
1521 bbox_max(1) = max(bbox_max(1), world_corner(1))
1522 bbox_max(2) = max(bbox_max(2), world_corner(2))
1523 bbox_max(3) = max(bbox_max(3), world_corner(3))
1524 end do
1525 end do
1526 end do
1527
1528 call get_bounding_indices(bbox_min(1), bbox_max(1), x_cc, il, ir)
1529 call get_bounding_indices(bbox_min(2), bbox_max(2), y_cc, jl, jr)
1530 call get_bounding_indices(bbox_min(3), bbox_max(3), z_cc, kl, kr)
1531
1532
1533# 804 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1534
1535# 804 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1536#if defined(MFC_OpenACC)
1537# 804 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1538!$acc parallel loop collapse(3) gang vector default(present) private(i, j, k, xyz_local, eta) copyin(patch_id, model_id, encoded_patch_id, center, inverse_rotation, offset, threshold)
1539# 804 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1540#elif defined(MFC_OpenMP)
1541# 804 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1542
1543# 804 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1544
1545# 804 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1546
1547# 804 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1548!$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, model_id, encoded_patch_id, center, inverse_rotation, offset, threshold)
1549# 804 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1550#endif
1551# 806 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1552 do i = il, ir
1553 do j = jl, jr
1554 do k = kl, kr
1555 xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), z_cc(k) - center(3)]
1556 xyz_local = matmul(inverse_rotation, xyz_local)
1557 xyz_local = xyz_local - offset
1558
1559 eta = f_model_is_inside_flat(gpu_ntrs(model_id), model_id, xyz_local)
1560
1561 if (eta > threshold) then
1562 ib_markers%sf(i, j, k) = encoded_patch_id
1563 end if
1564 end do
1565 end do
1566 end do
1567
1568# 821 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1569#if defined(MFC_OpenACC)
1570# 821 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1571!$acc end parallel loop
1572# 821 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1573#elif defined(MFC_OpenMP)
1574# 821 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1575
1576# 821 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1577!$omp end target teams loop
1578# 821 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1579#endif
1580
1581 end subroutine s_ib_3d_model
1582
1583 !> Compute a rotation matrix for converting to the rotating frame of the boundary
1584 subroutine s_update_ib_rotation_matrix(patch_id)
1585
1586
1587# 828 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1588#if MFC_OpenACC
1589# 828 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1590!$acc routine seq
1591# 828 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1592#elif MFC_OpenMP
1593# 828 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1594
1595# 828 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1596
1597# 828 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1598!$omp declare target device_type(any)
1599# 828 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1600#endif
1601
1602 integer, intent(in) :: patch_id
1603 real(wp), dimension(3, 3, 3) :: rotation
1604 real(wp) :: angle
1605
1606 ! construct the x, y, and z rotation matrices
1607
1608 if (num_dims == 3) then
1609 ! also compute the x and y axes in 3D
1610 angle = patch_ib(patch_id)%angles(1)
1611 rotation(1, 1,:) = [1._wp, 0._wp, 0._wp]
1612 rotation(1, 2,:) = [0._wp, cos(angle), -sin(angle)]
1613 rotation(1, 3,:) = [0._wp, sin(angle), cos(angle)]
1614
1615 angle = patch_ib(patch_id)%angles(2)
1616 rotation(2, 1,:) = [cos(angle), 0._wp, sin(angle)]
1617 rotation(2, 2,:) = [0._wp, 1._wp, 0._wp]
1618 rotation(2, 3,:) = [-sin(angle), 0._wp, cos(angle)]
1619
1620 ! apply the y rotation to the x rotation
1621 patch_ib(patch_id)%rotation_matrix(:,:) = matmul(rotation(1,:,:), rotation(2,:,:))
1622 patch_ib(patch_id)%rotation_matrix_inverse(:,:) = matmul(transpose(rotation(2,:,:)), transpose(rotation(1,:,:)))
1623 end if
1624
1625 ! z component first, since it applies in 2D and 3D
1626 angle = patch_ib(patch_id)%angles(3)
1627 rotation(3, 1,:) = [cos(angle), -sin(angle), 0._wp]
1628 rotation(3, 2,:) = [sin(angle), cos(angle), 0._wp]
1629 rotation(3, 3,:) = [0._wp, 0._wp, 1._wp]
1630
1631 if (num_dims == 3) then
1632 ! apply the z rotation to the xy rotation in 3D
1633 patch_ib(patch_id)%rotation_matrix(:,:) = matmul(patch_ib(patch_id)%rotation_matrix(:,:), rotation(3,:,:))
1634 patch_ib(patch_id)%rotation_matrix_inverse(:,:) = matmul(transpose(rotation(3,:,:)), &
1635 & patch_ib(patch_id)%rotation_matrix_inverse(:,:))
1636 else
1637 ! write out only the z rotation in 2D
1638 patch_ib(patch_id)%rotation_matrix(:,:) = rotation(3,:,:)
1639 patch_ib(patch_id)%rotation_matrix_inverse(:,:) = transpose(rotation(3,:,:))
1640 end if
1641
1642 end subroutine s_update_ib_rotation_matrix
1643
1644 subroutine get_bounding_indices(left_bound, right_bound, cell_centers, left_index, right_index)
1645
1646 real(wp), intent(in) :: left_bound, right_bound
1647 integer, intent(inout) :: left_index, right_index
1648 real(wp), dimension(-buff_size:), intent(in) :: cell_centers
1649 integer :: itr_left, itr_middle, itr_right
1650
1651 itr_left = left_index
1652 itr_right = right_index
1653
1654 do while (itr_left + 1 < itr_right)
1655 itr_middle = (itr_left + itr_right)/2
1656 if (cell_centers(itr_middle) < left_bound) then
1657 itr_left = itr_middle
1658 else if (cell_centers(itr_middle) > left_bound) then
1659 itr_right = itr_middle
1660 else
1661 itr_left = itr_middle
1662 exit
1663 end if
1664 end do
1665 left_index = itr_left
1666
1667 itr_right = right_index
1668 do while (itr_left + 1 < itr_right)
1669 itr_middle = (itr_left + itr_right)/2
1670 if (cell_centers(itr_middle) < right_bound) then
1671 itr_left = itr_middle
1672 else if (cell_centers(itr_middle) > right_bound) then
1673 itr_right = itr_middle
1674 else
1675 itr_right = itr_middle
1676 exit
1677 end if
1678 end do
1679 right_index = itr_right
1680
1681 end subroutine get_bounding_indices
1682
1683 !> Encode the patch ID with a unique offset containing periodicity information
1684 subroutine s_encode_patch_periodicity(patch_id, x_periodicity, y_periodicity, z_periodicity, encoded_patch_id)
1685
1686 integer, intent(in) :: patch_id, x_periodicity, y_periodicity, z_periodicity
1687 integer, intent(out) :: encoded_patch_id
1688 integer :: temp_x_per, temp_y_per, temp_z_per, offset
1689
1690 encoded_patch_id = patch_id
1691
1692 temp_x_per = x_periodicity; if (x_periodicity == -1) temp_x_per = 2
1693 temp_y_per = y_periodicity; if (y_periodicity == -1) temp_y_per = 2
1694 temp_z_per = z_periodicity; if (z_periodicity == -1) temp_z_per = 2
1695
1696 offset = (num_gbl_ibs + 1)*temp_x_per + 3*(num_gbl_ibs + 1)*temp_y_per + 9*(num_gbl_ibs + 1)*temp_z_per
1697 encoded_patch_id = patch_id + offset
1698
1699 end subroutine s_encode_patch_periodicity
1700
1701 !> Decode the encoded ID to recover the original patch ID and periodicity
1702 subroutine s_decode_patch_periodicity(encoded_patch_id, patch_id, x_periodicity, y_periodicity, z_periodicity)
1703
1704
1705# 932 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1706#if MFC_OpenACC
1707# 932 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1708!$acc routine seq
1709# 932 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1710#elif MFC_OpenMP
1711# 932 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1712
1713# 932 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1714
1715# 932 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1716!$omp declare target device_type(any)
1717# 932 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1718#endif
1719
1720 integer, intent(in) :: encoded_patch_id
1721 integer, intent(out) :: patch_id
1722 integer, intent(out), optional :: x_periodicity, y_periodicity, z_periodicity
1723 integer :: offset, remainder, xp, yp, zp, base
1724
1725 base = num_gbl_ibs + 1
1726
1727 patch_id = mod(encoded_patch_id - 1, base) + 1
1728 offset = (encoded_patch_id - patch_id)/base
1729
1730 xp = mod(offset, 3)
1731 remainder = offset/3
1732 yp = mod(remainder, 3)
1733 zp = remainder/3
1734
1735 ! Reverse map: 2 -> -1, 0 -> 0, 1 -> 1
1736 if (present(x_periodicity) .and. present(y_periodicity) .and. present(z_periodicity)) then
1737 x_periodicity = xp; if (xp == 2) x_periodicity = -1
1738 y_periodicity = yp; if (yp == 2) y_periodicity = -1
1739 z_periodicity = zp; if (zp == 2) z_periodicity = -1
1740 end if
1741
1742 end subroutine s_decode_patch_periodicity
1743
1744 !> Determine the periodic wrapping bounds in each direction
1745 subroutine s_get_periodicities(xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper)
1746
1747 integer, intent(out) :: xp_lower, xp_upper, yp_lower, yp_upper
1748 integer, intent(out), optional :: zp_lower, zp_upper
1749
1750 ! check domain wraps in x, y
1751
1752# 967 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1753 ! check for periodicity
1754 if (ib_bc_x%beg == bc_periodic) then
1755 xp_lower = -1
1756 xp_upper = 1
1757 else
1758 ! if it is not periodic, then both elements are 0
1759 xp_lower = 0
1760 xp_upper = 0
1761 end if
1762# 967 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1763 ! check for periodicity
1764 if (ib_bc_y%beg == bc_periodic) then
1765 yp_lower = -1
1766 yp_upper = 1
1767 else
1768 ! if it is not periodic, then both elements are 0
1769 yp_lower = 0
1770 yp_upper = 0
1771 end if
1772# 977 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1773
1774 ! z only if 3D
1775 if (present(zp_lower) .and. p /= 0) then
1776 if (ib_bc_z%beg == bc_periodic) then
1777 zp_lower = -1
1778 zp_upper = 1
1779 else
1780 zp_lower = 0
1781 zp_upper = 0
1782 end if
1783 end if
1784
1785 end subroutine s_get_periodicities
1786
1787 !> Archimedes spiral function
1788 pure elemental function f_r(myth, offset, a)
1789
1790
1791# 994 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1792#if MFC_OpenACC
1793# 994 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1794!$acc routine seq
1795# 994 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1796#elif MFC_OpenMP
1797# 994 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1798
1799# 994 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1800
1801# 994 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1802!$omp declare target device_type(any)
1803# 994 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1804#endif
1805 real(wp), intent(in) :: myth, offset, a
1806 real(wp) :: b
1807 real(wp) :: f_r
1808
1809 ! r(th) = a + b*th
1810
1811 b = 2._wp*a/(2._wp*pi)
1812 f_r = a + b*myth + offset
1813
1814 end function f_r
1815
1816end module m_ib_patches
j
integer, intent(in) j
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::ib_airfoil
type(ib_airfoil_parameters), dimension(num_ib_airfoils_max) ib_airfoil
Per-airfoil NACA user inputs (namelist).
Definition m_global_parameters.fpp.f90:1039
m_global_parameters::patch_ib
type(ib_patch_parameters), dimension(num_ib_patches_max_namelist) patch_ib
Immersed boundary patch parameters.
Definition m_global_parameters.fpp.f90:1035
m_global_parameters::ib_airfoil_grids
type(ib_airfoil_grid), dimension(num_ib_airfoils_max) ib_airfoil_grids
Per-airfoil computed surface grids.
Definition m_global_parameters.fpp.f90:1040
m_global_parameters::dx
real(wp), dimension(:), allocatable, target dx
Definition m_global_parameters.fpp.f90:524
m_global_parameters::num_ibs
integer num_ibs
Definition m_global_parameters.fpp.f90:408
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:1189
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:1585
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:694
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:1468
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:622
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:1746
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:1282
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:1645
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:807
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:1361
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:1685
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:561
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:933
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:1703
m_ib_patches::s_initialize_ib_airfoils
subroutine, public s_initialize_ib_airfoils()
Initialize the NACA surface grids for all airfoil IB patches. Must be called after the grid is establ...
Definition m_ib_patches.fpp.f90:406
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:1013
m_ib_patches::f_r
pure elemental real(wp) function f_r(myth, offset, a)
Archimedes spiral function.
Definition m_ib_patches.fpp.f90:1789
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:1099
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