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# 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
400contains
401
402 !> Apply all immersed boundary patch geometries to mark interior cells in the IB marker array
403 impure subroutine s_apply_ib_patches(ib_markers)
404
405 type(integer_field), intent(inout) :: ib_markers
406 integer :: i, xp, yp, zp !< iterators
407 integer :: xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper !< periodic bounds
408
409 ! 3D Patch Geometries
410
411 if (p > 0) then
412 !> IB Patches
413 !> @{
414 call s_get_periodicities(xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper)
415 do xp = xp_lower, xp_upper
416 do yp = yp_lower, yp_upper
417 do zp = zp_lower, zp_upper
418 do i = 1, num_ibs
419 if (patch_ib(i)%geometry == 8) then
420 call s_ib_sphere(i, ib_markers, xp, yp, zp)
421 else if (patch_ib(i)%geometry == 9) then
422 call s_ib_cuboid(i, ib_markers, xp, yp, zp)
423 else if (patch_ib(i)%geometry == 10) then
424 call s_ib_cylinder(i, ib_markers, xp, yp, zp)
425 else if (patch_ib(i)%geometry == 11) then
426 call s_ib_3d_airfoil(i, ib_markers, xp, yp, zp)
427 else if (patch_ib(i)%geometry == 12) then
428 call s_ib_3d_model(i, ib_markers, xp, yp, zp)
429 end if
430 end do
431 end do
432 end do
433 end do
434 !> @}
435
436 ! 2D Patch Geometries
437 else if (n > 0) then
438 !> IB Patches
439 !> @{
440 call s_get_periodicities(xp_lower, xp_upper, yp_lower, yp_upper)
441 do xp = xp_lower, xp_upper
442 do yp = yp_lower, yp_upper
443 do i = 1, num_ibs
444 if (patch_ib(i)%geometry == 2) then
445 call s_ib_circle(i, ib_markers, xp, yp)
446 else if (patch_ib(i)%geometry == 3) then
447 call s_ib_rectangle(i, ib_markers, xp, yp)
448 else if (patch_ib(i)%geometry == 4) then
449 call s_ib_airfoil(i, ib_markers, xp, yp)
450 else if (patch_ib(i)%geometry == 5) then
451 call s_ib_model(i, ib_markers, xp, yp)
452 else if (patch_ib(i)%geometry == 6) then
453 call s_ib_ellipse(i, ib_markers, xp, yp)
454 end if
455 end do
456 end do
457 end do
458 !> @}
459 end if
460
461 end subroutine s_apply_ib_patches
462
463 !> Mark cells inside a circular immersed boundary
464 subroutine s_ib_circle(patch_id, ib_markers, xp, yp)
465
466 integer, intent(in) :: patch_id
467 integer, intent(in) :: xp, yp !< integers containing the periodicity projection information
468 type(integer_field), intent(inout) :: ib_markers
469 real(wp), dimension(1:2) :: center
470 real(wp) :: radius
471 integer :: i, j, il, ir, jl, jr !< Generic loop iterators
472 integer :: encoded_patch_id
473
474 ! Transferring the circular patch's radius, centroid, smearing patch identity and smearing coefficient information
475
476 center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg)
477 center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg)
478 radius = patch_ib(patch_id)%radius
479
480 ! encode the periodicity information into the patch_id
481 call s_encode_patch_periodicity(patch_id, xp, yp, 0, encoded_patch_id)
482
483 ! find the indices to the left and right of the IB in i, j, k
484 il = -gp_layers - 1
485 jl = -gp_layers - 1
486 ir = m + gp_layers + 1
487 jr = n + gp_layers + 1
488 call get_bounding_indices(center(1) - radius, center(1) + radius, x_cc, il, ir)
489 call get_bounding_indices(center(2) - radius, center(2) + radius, y_cc, jl, jr)
490
491 ! Assign primitive variables if circle covers cell and patch has write permission
492
493
494# 119 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
495
496# 119 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
497#if defined(MFC_OpenACC)
498# 119 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
499!$acc parallel loop collapse(2) gang vector default(present) private(i, j) copyin(encoded_patch_id, center, radius)
500# 119 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
501#elif defined(MFC_OpenMP)
502# 119 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
503
504# 119 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
505
506# 119 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
507
508# 119 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
509!$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)
510# 119 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
511#endif
512 do j = jl, jr
513 do i = il, ir
514 if ((x_cc(i) - center(1))**2 + (y_cc(j) - center(2))**2 <= radius**2) then
515 ib_markers%sf(i, j, 0) = encoded_patch_id
516 end if
517 end do
518 end do
519
520# 127 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
521#if defined(MFC_OpenACC)
522# 127 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
523!$acc end parallel loop
524# 127 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
525#elif defined(MFC_OpenMP)
526# 127 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
527
528# 127 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
529!$omp end target teams loop
530# 127 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
531#endif
532
533 end subroutine s_ib_circle
534
535 !> Mark cells inside a 2D NACA 4-digit airfoil immersed boundary
536 subroutine s_ib_airfoil(patch_id, ib_markers, xp, yp)
537
538 integer, intent(in) :: patch_id
539 type(integer_field), intent(inout) :: ib_markers
540 integer, intent(in) :: xp, yp !< integers containing the periodicity projection information
541 real(wp) :: f, ca_in, pa, ma, ta
542 real(wp) :: xa, yt, xu, yu, xl, yl, xc, yc, dycdxc, sin_c, cos_c
543 integer :: i, j, k, il, ir, jl, jr
544 integer :: Np1, Np2
545 integer :: encoded_patch_id
546 real(wp), dimension(1:3) :: xy_local, offset !< x and y coordinates in local IB frame
547 real(wp), dimension(1:2) :: center !< x and y coordinates in local IB frame
548 real(wp), dimension(1:3,1:3) :: inverse_rotation
549
550 center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg)
551 center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg)
552 ca_in = patch_ib(patch_id)%c
553 pa = patch_ib(patch_id)%p
554 ma = patch_ib(patch_id)%m
555 ta = patch_ib(patch_id)%t
556 inverse_rotation(:,:) = patch_ib(patch_id)%rotation_matrix_inverse(:,:)
557 offset(:) = patch_ib(patch_id)%centroid_offset(:)
558
559 np1 = int((pa*ca_in/dx(0))*20)
560 np2 = int(((ca_in - pa*ca_in)/dx(0))*20)
561 np = np1 + np2 + 1
562
563# 158 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
564#if defined(MFC_OpenACC)
565# 158 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
566!$acc update device(Np)
567# 158 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
568#elif defined(MFC_OpenMP)
569# 158 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
570!$omp target update to(Np)
571# 158 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
572#endif
573
574 if (.not. allocated(airfoil_grid_u)) then
575#ifdef MFC_DEBUG
576# 161 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
577 block
578# 161 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
579 use iso_fortran_env, only: output_unit
580# 161 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
581
582# 161 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
583 print *, 'm_ib_patches.fpp:161: ', '@:ALLOCATE(airfoil_grid_u(1:Np))'
584# 161 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
585
586# 161 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
587 call flush (output_unit)
588# 161 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
589 end block
590# 161 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
591#endif
592# 161 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
593 allocate (airfoil_grid_u(1:np))
594# 161 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
595
596# 161 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
597
598# 161 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
599#if defined(MFC_OpenACC)
600# 161 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
601!$acc enter data create(airfoil_grid_u)
602# 161 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
603#elif defined(MFC_OpenMP)
604# 161 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
605!$omp target enter data map(always,alloc:airfoil_grid_u)
606# 161 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
607#endif
608#ifdef MFC_DEBUG
609# 162 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
610 block
611# 162 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
612 use iso_fortran_env, only: output_unit
613# 162 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
614
615# 162 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
616 print *, 'm_ib_patches.fpp:162: ', '@:ALLOCATE(airfoil_grid_l(1:Np))'
617# 162 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
618
619# 162 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
620 call flush (output_unit)
621# 162 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
622 end block
623# 162 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
624#endif
625# 162 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
626 allocate (airfoil_grid_l(1:np))
627# 162 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
628
629# 162 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
630
631# 162 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
632#if defined(MFC_OpenACC)
633# 162 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
634!$acc enter data create(airfoil_grid_l)
635# 162 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
636#elif defined(MFC_OpenMP)
637# 162 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
638!$omp target enter data map(always,alloc:airfoil_grid_l)
639# 162 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
640#endif
641
642 ! TODO :: The below instantiations are already handled by the loop below
643 airfoil_grid_u(1)%x = 0._wp
644 airfoil_grid_u(1)%y = 0._wp
645
646 airfoil_grid_l(1)%x = 0._wp
647 airfoil_grid_l(1)%y = 0._wp
648
649 do i = 1, np1 + np2 - 1
650 ! TODO :: This allocates the upper and lower airfoil arrays, and does not need to be performed each time the IB
651 ! markers are updated. Place this as a separate subroutine.
652 if (i <= np1) then
653 xc = i*(pa*ca_in/np1)
654 xa = xc/ca_in
655 yc = (ma/pa**2)*(2*pa*xa - xa**2)
656 dycdxc = (2*ma/pa**2)*(pa - xa)
657 else
658 xc = pa*ca_in + (i - np1)*((ca_in - pa*ca_in)/np2)
659 xa = xc/ca_in
660 yc = (ma/(1 - pa)**2)*(1 - 2*pa + 2*pa*xa - xa**2)
661 dycdxc = (2*ma/(1 - pa)**2)*(pa - xa)
662 end if
663
664 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)
665 sin_c = dycdxc/(1 + dycdxc**2)**0.5_wp
666 cos_c = 1/(1 + dycdxc**2)**0.5_wp
667
668 xu = xa - yt*sin_c
669 yu = yc + yt*cos_c
670
671 xl = xa + yt*sin_c
672 yl = yc - yt*cos_c
673
674 xu = xu*ca_in
675 yu = yu*ca_in
676
677 xl = xl*ca_in
678 yl = yl*ca_in
679
680 airfoil_grid_u(i + 1)%x = xu
681 airfoil_grid_u(i + 1)%y = yu
682
683 airfoil_grid_l(i + 1)%x = xl
684 airfoil_grid_l(i + 1)%y = yl
685 end do
686
687 airfoil_grid_u(np)%x = ca_in
688 airfoil_grid_u(np)%y = 0._wp
689
690 airfoil_grid_l(np)%x = ca_in
691 airfoil_grid_l(np)%y = 0._wp
692
693
694# 215 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
695#if defined(MFC_OpenACC)
696# 215 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
697!$acc update device(airfoil_grid_l, airfoil_grid_u)
698# 215 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
699#elif defined(MFC_OpenMP)
700# 215 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
701!$omp target update to(airfoil_grid_l, airfoil_grid_u)
702# 215 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
703#endif
704 end if
705
706 ! encode the periodicity information into the patch_id
707 call s_encode_patch_periodicity(patch_id, xp, yp, 0, encoded_patch_id)
708
709 ! find the indices to the left and right of the IB in i, j, k
710 il = -gp_layers - 1
711 jl = -gp_layers - 1
712 ir = m + gp_layers + 1
713 jr = n + gp_layers + 1
714 ! maximum distance any marker can be from the center is the length of the airfoil
715 call get_bounding_indices(center(1) - ca_in, center(1) + ca_in, x_cc, il, ir)
716 call get_bounding_indices(center(2) - ca_in, center(2) + ca_in, y_cc, jl, jr)
717
718
719# 230 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
720
721# 230 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
722#if defined(MFC_OpenACC)
723# 230 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
724!$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)
725# 230 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
726#elif defined(MFC_OpenMP)
727# 230 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
728
729# 230 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
730
731# 230 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
732
733# 230 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
734!$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)
735# 230 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
736#endif
737# 232 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
738 do j = jl, jr
739 do i = il, ir
740 xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp] ! get coordinate frame centered on IB
741 xy_local = matmul(inverse_rotation, xy_local) ! rotate the frame into the IB's coordinates
742 xy_local = xy_local - offset ! airfoils are a patch that require a centroid offset
743
744 if (xy_local(1) >= 0._wp .and. xy_local(1) <= ca_in) then
745 xa = xy_local(1)/ca_in
746 if (xa <= pa) then
747 yc = (ma/pa**2)*(2*pa*xa - xa**2)
748 dycdxc = (2*ma/pa**2)*(pa - xa)
749 else
750 yc = (ma/(1 - pa)**2)*(1 - 2*pa + 2*pa*xa - xa**2)
751 dycdxc = (2*ma/(1 - pa)**2)*(pa - xa)
752 end if
753 if (xy_local(2) >= 0._wp) then
754 k = 1
755 do while (airfoil_grid_u(k)%x < xy_local(1) .and. k <= np)
756 k = k + 1
757 end do
758 if (f_approx_equal(airfoil_grid_u(k)%x, xy_local(1))) then
759 if (xy_local(2) <= airfoil_grid_u(k)%y) then
760 ib_markers%sf(i, j, 0) = encoded_patch_id
761 end if
762 else
763 f = (airfoil_grid_u(k)%x - xy_local(1))/(airfoil_grid_u(k)%x - airfoil_grid_u(k - 1)%x)
764 if (xy_local(2) <= ((1._wp - f)*airfoil_grid_u(k)%y + f*airfoil_grid_u(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 (airfoil_grid_l(k)%x < xy_local(1))
771 k = k + 1
772 end do
773 if (f_approx_equal(airfoil_grid_l(k)%x, xy_local(1))) then
774 if (xy_local(2) >= airfoil_grid_l(k)%y) then
775 ib_markers%sf(i, j, 0) = encoded_patch_id
776 end if
777 else
778 f = (airfoil_grid_l(k)%x - xy_local(1))/(airfoil_grid_l(k)%x - airfoil_grid_l(k - 1)%x)
779
780 if (xy_local(2) >= ((1._wp - f)*airfoil_grid_l(k)%y + f*airfoil_grid_l(k - 1)%y)) then
781 ib_markers%sf(i, j, 0) = encoded_patch_id
782 end if
783 end if
784 end if
785 end if
786 end do
787 end do
788
789# 282 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
790#if defined(MFC_OpenACC)
791# 282 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
792!$acc end parallel loop
793# 282 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
794#elif defined(MFC_OpenMP)
795# 282 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
796
797# 282 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
798!$omp end target teams loop
799# 282 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
800#endif
801
802 end subroutine s_ib_airfoil
803
804 !> Mark cells inside a 3D extruded NACA 4-digit airfoil immersed boundary with finite span
805 subroutine s_ib_3d_airfoil(patch_id, ib_markers, xp, yp, zp)
806
807 integer, intent(in) :: patch_id
808 type(integer_field), intent(inout) :: ib_markers
809 integer, intent(in) :: xp, yp, zp !< integers containing the periodicity projection information
810 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
811 integer :: i, j, k, l, il, ir, jl, jr, ll, lr
812 integer :: Np1, Np2
813 integer :: encoded_patch_id
814 real(wp), dimension(1:3) :: xyz_local, center, offset !< x, y, z coordinates in local IB frame
815 real(wp), dimension(1:3,1:3) :: inverse_rotation
816
817 center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg)
818 center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg)
819 center(3) = patch_ib(patch_id)%z_centroid + real(zp, wp)*(z_domain%end - z_domain%beg)
820 lz = patch_ib(patch_id)%length_z
821 ca_in = patch_ib(patch_id)%c
822 pa = patch_ib(patch_id)%p
823 ma = patch_ib(patch_id)%m
824 ta = patch_ib(patch_id)%t
825 inverse_rotation(:,:) = patch_ib(patch_id)%rotation_matrix_inverse(:,:)
826 offset(:) = patch_ib(patch_id)%centroid_offset(:)
827
828 z_max = lz/2
829 z_min = -lz/2
830
831 np1 = int((pa*ca_in/dx(0))*20)
832 np2 = int(((ca_in - pa*ca_in)/dx(0))*20)
833 np = np1 + np2 + 1
834
835# 316 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
836#if defined(MFC_OpenACC)
837# 316 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
838!$acc update device(Np)
839# 316 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
840#elif defined(MFC_OpenMP)
841# 316 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
842!$omp target update to(Np)
843# 316 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
844#endif
845
846 if (.not. allocated(airfoil_grid_u)) then
847#ifdef MFC_DEBUG
848# 319 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
849 block
850# 319 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
851 use iso_fortran_env, only: output_unit
852# 319 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
853
854# 319 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
855 print *, 'm_ib_patches.fpp:319: ', '@:ALLOCATE(airfoil_grid_u(1:Np))'
856# 319 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
857
858# 319 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
859 call flush (output_unit)
860# 319 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
861 end block
862# 319 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
863#endif
864# 319 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
865 allocate (airfoil_grid_u(1:np))
866# 319 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
867
868# 319 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
869
870# 319 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
871#if defined(MFC_OpenACC)
872# 319 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
873!$acc enter data create(airfoil_grid_u)
874# 319 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
875#elif defined(MFC_OpenMP)
876# 319 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
877!$omp target enter data map(always,alloc:airfoil_grid_u)
878# 319 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
879#endif
880#ifdef MFC_DEBUG
881# 320 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
882 block
883# 320 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
884 use iso_fortran_env, only: output_unit
885# 320 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
886
887# 320 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
888 print *, 'm_ib_patches.fpp:320: ', '@:ALLOCATE(airfoil_grid_l(1:Np))'
889# 320 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
890
891# 320 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
892 call flush (output_unit)
893# 320 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
894 end block
895# 320 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
896#endif
897# 320 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
898 allocate (airfoil_grid_l(1:np))
899# 320 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
900
901# 320 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
902
903# 320 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
904#if defined(MFC_OpenACC)
905# 320 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
906!$acc enter data create(airfoil_grid_l)
907# 320 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
908#elif defined(MFC_OpenMP)
909# 320 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
910!$omp target enter data map(always,alloc:airfoil_grid_l)
911# 320 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
912#endif
913
914 airfoil_grid_u(1)%x = 0._wp
915 airfoil_grid_u(1)%y = 0._wp
916
917 airfoil_grid_l(1)%x = 0._wp
918 airfoil_grid_l(1)%y = 0._wp
919
920 do i = 1, np1 + np2 - 1
921 if (i <= np1) then
922 xc = i*(pa*ca_in/np1)
923 xa = xc/ca_in
924 yc = (ma/pa**2)*(2*pa*xa - xa**2)
925 dycdxc = (2*ma/pa**2)*(pa - xa)
926 else
927 xc = pa*ca_in + (i - np1)*((ca_in - pa*ca_in)/np2)
928 xa = xc/ca_in
929 yc = (ma/(1 - pa)**2)*(1 - 2*pa + 2*pa*xa - xa**2)
930 dycdxc = (2*ma/(1 - pa)**2)*(pa - xa)
931 end if
932
933 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)
934 sin_c = dycdxc/(1 + dycdxc**2)**0.5_wp
935 cos_c = 1/(1 + dycdxc**2)**0.5_wp
936
937 xu = xa - yt*sin_c
938 yu = yc + yt*cos_c
939
940 xl = xa + yt*sin_c
941 yl = yc - yt*cos_c
942
943 xu = xu*ca_in
944 yu = yu*ca_in
945
946 xl = xl*ca_in
947 yl = yl*ca_in
948
949 airfoil_grid_u(i + 1)%x = xu
950 airfoil_grid_u(i + 1)%y = yu
951
952 airfoil_grid_l(i + 1)%x = xl
953 airfoil_grid_l(i + 1)%y = yl
954 end do
955
956 airfoil_grid_u(np)%x = ca_in
957 airfoil_grid_u(np)%y = 0._wp
958
959 airfoil_grid_l(np)%x = ca_in
960 airfoil_grid_l(np)%y = 0._wp
961
962
963# 370 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
964#if defined(MFC_OpenACC)
965# 370 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
966!$acc update device(airfoil_grid_l, airfoil_grid_u)
967# 370 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
968#elif defined(MFC_OpenMP)
969# 370 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
970!$omp target update to(airfoil_grid_l, airfoil_grid_u)
971# 370 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
972#endif
973 end if
974
975 ! encode the periodicity information into the patch_id
976 call s_encode_patch_periodicity(patch_id, xp, yp, zp, encoded_patch_id)
977
978 ! find the indices to the left and right of the IB in i, j, k
979 il = -gp_layers - 1
980 jl = -gp_layers - 1
981 ll = -gp_layers - 1
982 ir = m + gp_layers + 1
983 jr = n + gp_layers + 1
984 lr = p + gp_layers + 1
985 ! maximum distance any marker can be from the center is the length of the airfoil
986 call get_bounding_indices(center(1) - ca_in, center(1) + ca_in, x_cc, il, ir)
987 call get_bounding_indices(center(2) - ca_in, center(2) + ca_in, y_cc, jl, jr)
988 call get_bounding_indices(center(3) - ca_in, center(3) + ca_in, z_cc, ll, lr)
989
990
991# 388 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
992
993# 388 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
994#if defined(MFC_OpenACC)
995# 388 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
996!$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)
997# 388 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
998#elif defined(MFC_OpenMP)
999# 388 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1000
1001# 388 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1002
1003# 388 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1004
1005# 388 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1006!$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)
1007# 388 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1008#endif
1009# 390 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1010 do l = ll, lr
1011 do j = jl, jr
1012 do i = il, ir
1013 ! get coordinate frame centered on IB
1014 xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), z_cc(l) - center(3)]
1015 xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinates
1016 xyz_local = xyz_local - offset ! airfoils are a patch that require a centroid offset
1017
1018 if (xyz_local(3) >= z_min .and. xyz_local(3) <= z_max) then
1019 if (xyz_local(1) >= 0._wp .and. xyz_local(1) <= ca_in) then
1020 if (xyz_local(2) >= 0._wp) then
1021 k = 1
1022 do while (airfoil_grid_u(k)%x < xyz_local(1))
1023 k = k + 1
1024 end do
1025 if (f_approx_equal(airfoil_grid_u(k)%x, xyz_local(1))) then
1026 if (xyz_local(2) <= airfoil_grid_u(k)%y) then
1027 ! IB
1028 ib_markers%sf(i, j, l) = encoded_patch_id
1029 end if
1030 else
1031 f = (airfoil_grid_u(k)%x - xyz_local(1))/(airfoil_grid_u(k)%x - airfoil_grid_u(k - 1)%x)
1032 if (xyz_local(2) <= ((1._wp - f)*airfoil_grid_u(k)%y + f*airfoil_grid_u(k - 1)%y)) then
1033 ib_markers%sf(i, j, l) = encoded_patch_id
1034 end if
1035 end if
1036 else
1037 k = 1
1038 do while (airfoil_grid_l(k)%x < xyz_local(1))
1039 k = k + 1
1040 end do
1041 if (f_approx_equal(airfoil_grid_l(k)%x, xyz_local(1))) then
1042 if (xyz_local(2) >= airfoil_grid_l(k)%y) then
1043 ib_markers%sf(i, j, l) = encoded_patch_id
1044 end if
1045 else
1046 f = (airfoil_grid_l(k)%x - xyz_local(1))/(airfoil_grid_l(k)%x - airfoil_grid_l(k - 1)%x)
1047
1048 if (xyz_local(2) >= ((1._wp - f)*airfoil_grid_l(k)%y + f*airfoil_grid_l(k - 1)%y)) then
1049 ib_markers%sf(i, j, l) = encoded_patch_id
1050 end if
1051 end if
1052 end if
1053 end if
1054 end if
1055 end do
1056 end do
1057 end do
1058
1059# 438 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1060#if defined(MFC_OpenACC)
1061# 438 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1062!$acc end parallel loop
1063# 438 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1064#elif defined(MFC_OpenMP)
1065# 438 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1066
1067# 438 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1068!$omp end target teams loop
1069# 438 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1070#endif
1071
1072 end subroutine s_ib_3d_airfoil
1073
1074 !> Mark cells inside a rectangular immersed boundary
1075 subroutine s_ib_rectangle(patch_id, ib_markers, xp, yp)
1076
1077 integer, intent(in) :: patch_id
1078 type(integer_field), intent(inout) :: ib_markers
1079 integer, intent(in) :: xp, yp !< integers containing the periodicity projection information
1080 integer :: i, j, il, ir, jl, jr !< generic loop iterators
1081 integer :: encoded_patch_id
1082 real(wp) :: corner_distance !< Equation of state parameters
1083 real(wp), dimension(1:3) :: xy_local !< x and y coordinates in local IB frame
1084 real(wp), dimension(1:2) :: length, center !< x and y coordinates in local IB frame
1085 real(wp), dimension(1:3,1:3) :: inverse_rotation
1086
1087 ! Transferring the rectangle's centroid and length information
1088
1089 center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg)
1090 center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg)
1091 length(1) = patch_ib(patch_id)%length_x
1092 length(2) = patch_ib(patch_id)%length_y
1093 inverse_rotation(:,:) = patch_ib(patch_id)%rotation_matrix_inverse(:,:)
1094
1095 ! encode the periodicity information into the patch_id
1096 call s_encode_patch_periodicity(patch_id, xp, yp, 0, encoded_patch_id)
1097
1098 ! find the indices to the left and right of the IB in i, j, k
1099 il = -gp_layers - 1
1100 jl = -gp_layers - 1
1101 ir = m + gp_layers + 1
1102 jr = n + gp_layers + 1
1103 corner_distance = sqrt(dot_product(length, length))/2._wp ! maximum distance any marker can be from the center
1104 call get_bounding_indices(center(1) - corner_distance, center(1) + corner_distance, x_cc, il, ir)
1105 call get_bounding_indices(center(2) - corner_distance, center(2) + corner_distance, y_cc, jl, jr)
1106
1107 ! Assign primitive variables if rectangle covers cell and patch has write permission
1108
1109# 476 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1110
1111# 476 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1112#if defined(MFC_OpenACC)
1113# 476 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1114!$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)
1115# 476 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1116#elif defined(MFC_OpenMP)
1117# 476 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1118
1119# 476 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1120
1121# 476 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1122
1123# 476 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1124!$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)
1125# 476 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1126#endif
1127# 478 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1128 do j = jl, jr
1129 do i = il, ir
1130 ! get the x and y coordinates in the local IB frame
1131 xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp]
1132 xy_local = matmul(inverse_rotation, xy_local)
1133
1134 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) &
1135 & .and. 0.5_wp*length(2) >= xy_local(2)) then
1136 ! Updating the patch identities bookkeeping variable
1137 ib_markers%sf(i, j, 0) = encoded_patch_id
1138 end if
1139 end do
1140 end do
1141
1142# 491 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1143#if defined(MFC_OpenACC)
1144# 491 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1145!$acc end parallel loop
1146# 491 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1147#elif defined(MFC_OpenMP)
1148# 491 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1149
1150# 491 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1151!$omp end target teams loop
1152# 491 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1153#endif
1154
1155 end subroutine s_ib_rectangle
1156
1157 !> Mark cells inside a spherical immersed boundary
1158 subroutine s_ib_sphere(patch_id, ib_markers, xp, yp, zp)
1159
1160 integer, intent(in) :: patch_id
1161 type(integer_field), intent(inout) :: ib_markers
1162 integer, intent(in) :: xp, yp, zp !< integers containing the periodicity projection information
1163 ! Generic loop iterators
1164 integer :: i, j, k
1165 integer :: il, ir, jl, jr, kl, kr
1166 integer :: encoded_patch_id
1167 real(wp) :: radius
1168 real(wp), dimension(1:3) :: center
1169
1170 ! Variables to initialize the pressure field that corresponds to the bubble-collapse test case found in Tiwari et al. (2013)
1171
1172 ! Transferring spherical patch's radius, centroid, smoothing patch identity and smoothing coefficient information
1173
1174 center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg)
1175 center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg)
1176 center(3) = patch_ib(patch_id)%z_centroid + real(zp, wp)*(z_domain%end - z_domain%beg)
1177 radius = patch_ib(patch_id)%radius
1178
1179 ! completely skip particles no in the domain
1180 if (center(1) - radius > x_cc(m + gp_layers + 1) .or. center(1) + radius < x_cc(-gp_layers - 1) .or. center(2) &
1181 & - radius > y_cc(n + gp_layers + 1) .or. center(2) + radius < y_cc(-gp_layers - 1) .or. center(3) - radius > z_cc(p &
1182 & + gp_layers + 1) .or. center(3) + radius < z_cc(-gp_layers - 1)) then
1183 return
1184 end if
1185
1186 ! encode the periodicity information into the patch_id
1187 call s_encode_patch_periodicity(patch_id, xp, yp, zp, encoded_patch_id)
1188
1189 ! find the indices to the left and right of the IB in i, j, k
1190 il = -gp_layers - 1
1191 jl = -gp_layers - 1
1192 kl = -gp_layers - 1
1193 ir = m + gp_layers + 1
1194 jr = n + gp_layers + 1
1195 kr = p + gp_layers + 1
1196 call get_bounding_indices(center(1) - radius, center(1) + radius, x_cc, il, ir)
1197 call get_bounding_indices(center(2) - radius, center(2) + radius, y_cc, jl, jr)
1198 call get_bounding_indices(center(3) - radius, center(3) + radius, z_cc, kl, kr)
1199
1200 ! Checking whether the sphere covers a particular cell in the domain and verifying whether the current patch has permission
1201 ! to write to that cell. If both queries check out, the primitive variables of the current patch are assigned to this cell.
1202
1203# 540 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1204
1205# 540 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1206#if defined(MFC_OpenACC)
1207# 540 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1208!$acc parallel loop collapse(3) gang vector default(present) private(i, j, k) copyin(encoded_patch_id, center, radius)
1209# 540 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1210#elif defined(MFC_OpenMP)
1211# 540 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1212
1213# 540 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1214
1215# 540 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1216
1217# 540 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1218!$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)
1219# 540 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1220#endif
1221 do k = kl, kr
1222 do j = jl, jr
1223 do i = il, ir
1224 ! Updating the patch identities bookkeeping variable
1225 if (((x_cc(i) - center(1))**2 + (y_cc(j) - center(2))**2 + (z_cc(k) - center(3))**2 <= radius**2)) then
1226 ib_markers%sf(i, j, k) = encoded_patch_id
1227 end if
1228 end do
1229 end do
1230 end do
1231
1232# 551 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1233#if defined(MFC_OpenACC)
1234# 551 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1235!$acc end parallel loop
1236# 551 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1237#elif defined(MFC_OpenMP)
1238# 551 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1239
1240# 551 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1241!$omp end target teams loop
1242# 551 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1243#endif
1244
1245 end subroutine s_ib_sphere
1246
1247 !> Mark cells inside a cuboidal immersed boundary
1248 subroutine s_ib_cuboid(patch_id, ib_markers, xp, yp, zp)
1249
1250 integer, intent(in) :: patch_id
1251 type(integer_field), intent(inout) :: ib_markers
1252 integer, intent(in) :: xp, yp, zp !< integers containing the periodicity projection information
1253 integer :: i, j, k, ir, il, jr, jl, kr, kl !< Generic loop iterators
1254 integer :: encoded_patch_id
1255 real(wp), dimension(1:3) :: xyz_local, center, length !< x and y coordinates in local IB frame
1256 real(wp), dimension(1:3,1:3) :: inverse_rotation
1257 real(wp) :: corner_distance
1258
1259 ! Transferring the cuboid's centroid and length information
1260
1261 center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg)
1262 center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg)
1263 center(3) = patch_ib(patch_id)%z_centroid + real(zp, wp)*(z_domain%end - z_domain%beg)
1264 length(1) = patch_ib(patch_id)%length_x
1265 length(2) = patch_ib(patch_id)%length_y
1266 length(3) = patch_ib(patch_id)%length_z
1267 inverse_rotation(:,:) = patch_ib(patch_id)%rotation_matrix_inverse(:,:)
1268
1269 ! encode the periodicity information into the patch_id
1270 call s_encode_patch_periodicity(patch_id, xp, yp, zp, encoded_patch_id)
1271
1272 ! find the indices to the left and right of the IB in i, j, k
1273 il = -gp_layers - 1
1274 jl = -gp_layers - 1
1275 kl = -gp_layers - 1
1276 ir = m + gp_layers + 1
1277 jr = n + gp_layers + 1
1278 kr = p + gp_layers + 1
1279 corner_distance = sqrt(dot_product(length, length))/2._wp ! maximum distance any marker can be from the center
1280 call get_bounding_indices(center(1) - corner_distance, center(1) + corner_distance, x_cc, il, ir)
1281 call get_bounding_indices(center(2) - corner_distance, center(2) + corner_distance, y_cc, jl, jr)
1282 call get_bounding_indices(center(3) - corner_distance, center(3) + corner_distance, z_cc, kl, kr)
1283
1284 ! Checking whether the cuboid covers a particular cell in the domain and verifying whether the current patch has permission
1285 ! to write to to that cell. If both queries check out, the primitive variables of the current patch are assigned to this
1286 ! cell.
1287
1288# 595 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1289
1290# 595 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1291#if defined(MFC_OpenACC)
1292# 595 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1293!$acc parallel loop collapse(3) gang vector default(present) private(i, j, k, xyz_local) copyin(encoded_patch_id, center, length, inverse_rotation)
1294# 595 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1295#elif defined(MFC_OpenMP)
1296# 595 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1297
1298# 595 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1299
1300# 595 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1301
1302# 595 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1303!$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, length, inverse_rotation)
1304# 595 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1305#endif
1306# 597 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1307 do k = kl, kr
1308 do j = jl, jr
1309 do i = il, ir
1310 xyz_local = [x_cc(i), y_cc(j), z_cc(k)] - center ! get coordinate frame centered on IB
1311 xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinates
1312
1313 if (-0.5*length(1) <= xyz_local(1) .and. 0.5*length(1) >= xyz_local(1) .and. -0.5*length(2) <= xyz_local(2) &
1314 & .and. 0.5*length(2) >= xyz_local(2) .and. -0.5*length(3) <= xyz_local(3) .and. 0.5*length(3) &
1315 & >= xyz_local(3)) then
1316 ! Updating the patch identities bookkeeping variable
1317 ib_markers%sf(i, j, k) = encoded_patch_id
1318 end if
1319 end do
1320 end do
1321 end do
1322
1323# 612 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1324#if defined(MFC_OpenACC)
1325# 612 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1326!$acc end parallel loop
1327# 612 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1328#elif defined(MFC_OpenMP)
1329# 612 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1330
1331# 612 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1332!$omp end target teams loop
1333# 612 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1334#endif
1335
1336 end subroutine s_ib_cuboid
1337
1338 !> Mark cells inside a cylindrical immersed boundary
1339 subroutine s_ib_cylinder(patch_id, ib_markers, xp, yp, zp)
1340
1341 integer, intent(in) :: patch_id
1342 type(integer_field), intent(inout) :: ib_markers
1343 integer, intent(in) :: xp, yp, zp !< integers containing the periodicity projection information
1344 integer :: i, j, k, il, ir, jl, jr, kl, kr !< Generic loop iterators
1345 integer :: encoded_patch_id
1346 real(wp) :: radius
1347 real(wp), dimension(1:3) :: xyz_local, center, length !< x and y coordinates in local IB frame
1348 real(wp), dimension(1:3,1:3) :: inverse_rotation
1349 real(wp) :: corner_distance
1350
1351 ! Transferring the cylindrical patch's centroid, length, radius,
1352
1353 center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg)
1354 center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg)
1355 center(3) = patch_ib(patch_id)%z_centroid + real(zp, wp)*(z_domain%end - z_domain%beg)
1356 length(1) = patch_ib(patch_id)%length_x
1357 length(2) = patch_ib(patch_id)%length_y
1358 length(3) = patch_ib(patch_id)%length_z
1359 radius = patch_ib(patch_id)%radius
1360 inverse_rotation(:,:) = patch_ib(patch_id)%rotation_matrix_inverse(:,:)
1361
1362 ! encode the periodicity information into the patch_id
1363 call s_encode_patch_periodicity(patch_id, xp, yp, zp, encoded_patch_id)
1364
1365 il = -gp_layers - 1
1366 jl = -gp_layers - 1
1367 kl = -gp_layers - 1
1368 ir = m + gp_layers + 1
1369 jr = n + gp_layers + 1
1370 kr = p + gp_layers + 1
1371 corner_distance = sqrt(radius**2 + maxval(length)**2) ! distance to rim of cylinder
1372 call get_bounding_indices(center(1) - corner_distance, center(1) + corner_distance, x_cc, il, ir)
1373 call get_bounding_indices(center(2) - corner_distance, center(2) + corner_distance, y_cc, jl, jr)
1374 call get_bounding_indices(center(3) - corner_distance, center(3) + corner_distance, z_cc, kl, kr)
1375
1376 ! Checking whether the cylinder covers a particular cell in the domain and verifying whether the current patch has the
1377 ! permission to write to that cell. If both queries check out, the primitive variables of the current patch are assigned to
1378 ! this cell.
1379
1380# 657 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1381
1382# 657 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1383#if defined(MFC_OpenACC)
1384# 657 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1385!$acc parallel loop collapse(3) gang vector default(present) private(i, j, k, xyz_local) copyin(encoded_patch_id, center, length, radius, inverse_rotation)
1386# 657 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1387#elif defined(MFC_OpenMP)
1388# 657 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1389
1390# 657 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1391
1392# 657 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1393
1394# 657 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1395!$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, length, radius, inverse_rotation)
1396# 657 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1397#endif
1398# 659 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1399 do k = kl, kr
1400 do j = jl, jr
1401 do i = il, ir
1402 xyz_local = [x_cc(i), y_cc(j), z_cc(k)] - center ! get coordinate frame centered on IB
1403 xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinates
1404
1405 if (((.not. f_is_default(length(1)) .and. xyz_local(2)**2 + xyz_local(3)**2 <= radius**2 .and. &
1406 & -0.5_wp*length(1) <= xyz_local(1) .and. 0.5_wp*length(1) >= xyz_local(1)) &
1407 & .or. (.not. f_is_default(length(2)) .and. xyz_local(1)**2 + xyz_local(3)**2 <= radius**2 .and. &
1408 & -0.5_wp*length(2) <= xyz_local(2) .and. 0.5_wp*length(2) >= xyz_local(2)) &
1409 & .or. (.not. f_is_default(length(3)) .and. xyz_local(1)**2 + xyz_local(2)**2 <= radius**2 .and. &
1410 & -0.5_wp*length(3) <= xyz_local(3) .and. 0.5_wp*length(3) >= xyz_local(3)))) then
1411 ! Updating the patch identities bookkeeping variable
1412 ib_markers%sf(i, j, k) = encoded_patch_id
1413 end if
1414 end do
1415 end do
1416 end do
1417
1418# 677 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1419#if defined(MFC_OpenACC)
1420# 677 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1421!$acc end parallel loop
1422# 677 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1423#elif defined(MFC_OpenMP)
1424# 677 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1425
1426# 677 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1427!$omp end target teams loop
1428# 677 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1429#endif
1430
1431 end subroutine s_ib_cylinder
1432
1433 !> Mark cells inside a 2D elliptical immersed boundary
1434 subroutine s_ib_ellipse(patch_id, ib_markers, xp, yp)
1435
1436 integer, intent(in) :: patch_id
1437 type(integer_field), intent(inout) :: ib_markers
1438 integer, intent(in) :: xp, yp !< integers containing the periodicity projection information
1439 integer :: i, j, il, ir, jl, jr !< Generic loop iterators
1440 integer :: encoded_patch_id
1441 real(wp), dimension(1:3) :: xy_local !< x and y coordinates in local IB frame
1442 real(wp), dimension(1:2) :: ellipse_coeffs !< a and b in the ellipse coefficients
1443 real(wp), dimension(1:2) :: center !< x and y coordinates in local IB frame
1444 real(wp), dimension(1:3,1:3) :: inverse_rotation
1445
1446 ! Transferring the ellipse's centroid and length information
1447
1448 center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg)
1449 center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg)
1450 ellipse_coeffs(1) = 0.5_wp*patch_ib(patch_id)%length_x
1451 ellipse_coeffs(2) = 0.5_wp*patch_ib(patch_id)%length_y
1452 inverse_rotation(:,:) = patch_ib(patch_id)%rotation_matrix_inverse(:,:)
1453
1454 ! encode the periodicity information into the patch_id
1455 call s_encode_patch_periodicity(patch_id, xp, yp, 0, encoded_patch_id)
1456
1457 ! find the indices to the left and right of the IB in i, j, k
1458 il = -gp_layers - 1
1459 jl = -gp_layers - 1
1460 ir = m + gp_layers + 1
1461 jr = n + gp_layers + 1
1462 call get_bounding_indices(center(1) - maxval(ellipse_coeffs)*2._wp, center(1) + maxval(ellipse_coeffs)*2._wp, x_cc, il, ir)
1463 call get_bounding_indices(center(2) - maxval(ellipse_coeffs)*2._wp, center(2) + maxval(ellipse_coeffs)*2._wp, y_cc, jl, jr)
1464
1465 ! Checking whether the ellipse covers a particular cell in the domain
1466
1467# 714 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1468
1469# 714 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1470#if defined(MFC_OpenACC)
1471# 714 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1472!$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)
1473# 714 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1474#elif defined(MFC_OpenMP)
1475# 714 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1476
1477# 714 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1478
1479# 714 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1480
1481# 714 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1482!$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)
1483# 714 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1484#endif
1485# 716 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1486 do j = jl, jr
1487 do i = il, ir
1488 ! get the x and y coordinates in the local IB frame
1489 xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp]
1490 xy_local = matmul(inverse_rotation, xy_local)
1491
1492 ! Ellipse condition (x/a)^2 + (y/b)^2 <= 1
1493 if ((xy_local(1)/ellipse_coeffs(1))**2 + (xy_local(2)/ellipse_coeffs(2))**2 <= 1._wp) then
1494 ! Updating the patch identities bookkeeping variable
1495 ib_markers%sf(i, j, 0) = encoded_patch_id
1496 end if
1497 end do
1498 end do
1499
1500# 729 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1501#if defined(MFC_OpenACC)
1502# 729 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1503!$acc end parallel loop
1504# 729 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1505#elif defined(MFC_OpenMP)
1506# 729 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1507
1508# 729 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1509!$omp end target teams loop
1510# 729 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1511#endif
1512
1513 end subroutine s_ib_ellipse
1514
1515 !> The STL patch is a 2D geometry that is imported from an STL file.
1516 subroutine s_ib_model(patch_id, ib_markers, xp, yp)
1517
1518 integer, intent(in) :: patch_id
1519 type(integer_field), intent(inout) :: ib_markers
1520 integer, intent(in) :: xp, yp !< integers containing the periodicity projection information
1521 integer :: i, j, il, ir, jl, jr !< Generic loop iterators
1522 integer :: spc, encoded_patch_id
1523 integer :: cx, cy
1524 real(wp) :: lx(2), ly(2)
1525 real(wp), dimension(1:2) :: bbox_min, bbox_max
1526 real(wp), dimension(1:3) :: local_corner, world_corner
1527 real(wp) :: eta, threshold
1528 real(wp), dimension(1:3) :: offset
1529 real(wp), dimension(1:3) :: center, xy_local
1530 real(wp), dimension(1:3,1:3) :: inverse_rotation, rotation
1531
1532 center = 0._wp
1533 center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg)
1534 center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg)
1535 inverse_rotation(:,:) = patch_ib(patch_id)%rotation_matrix_inverse(:,:)
1536 rotation(:,:) = patch_ib(patch_id)%rotation_matrix(:,:)
1537 offset(:) = patch_ib(patch_id)%centroid_offset(:)
1538 spc = patch_ib(patch_id)%model_spc
1539 threshold = patch_ib(patch_id)%model_threshold
1540
1541 ! encode the periodicity information into the patch_id
1542 call s_encode_patch_periodicity(patch_id, xp, yp, 0, encoded_patch_id)
1543
1544 il = -gp_layers - 1
1545 jl = -gp_layers - 1
1546 ir = m + gp_layers + 1
1547 jr = n + gp_layers + 1
1548
1549 ! Local-space bounding box extents (min=1, max=2 in the third index)
1550 lx(1) = stl_bounding_boxes(patch_id, 1, 1) + offset(1)
1551 lx(2) = stl_bounding_boxes(patch_id, 1, 3) + offset(1)
1552 ly(1) = stl_bounding_boxes(patch_id, 2, 1) + offset(2)
1553 ly(2) = stl_bounding_boxes(patch_id, 2, 3) + offset(2)
1554
1555 bbox_min = 1e12
1556 bbox_max = -1e12
1557 ! Enumerate all 8 corners of the local bounding box, rotate to world space, track world-space AABB
1558 do cx = 1, 2
1559 do cy = 1, 2
1560 local_corner = [lx(cx), ly(cy), 0._wp]
1561 world_corner = matmul(rotation, local_corner) + center
1562 bbox_min(1) = min(bbox_min(1), world_corner(1))
1563 bbox_min(2) = min(bbox_min(2), world_corner(2))
1564 bbox_max(1) = max(bbox_max(1), world_corner(1))
1565 bbox_max(2) = max(bbox_max(2), world_corner(2))
1566 end do
1567 end do
1568
1569 call get_bounding_indices(bbox_min(1), bbox_max(1), x_cc, il, ir)
1570 call get_bounding_indices(bbox_min(2), bbox_max(2), y_cc, jl, jr)
1571
1572
1573# 790 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1574
1575# 790 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1576#if defined(MFC_OpenACC)
1577# 790 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1578!$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)
1579# 790 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1580#elif defined(MFC_OpenMP)
1581# 790 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1582
1583# 790 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1584
1585# 790 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1586
1587# 790 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1588!$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)
1589# 790 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1590#endif
1591# 792 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1592 do i = il, ir
1593 do j = jl, jr
1594 xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp]
1595 xy_local = matmul(inverse_rotation, xy_local)
1596 xy_local = xy_local - offset
1597
1598 eta = f_model_is_inside_flat(gpu_ntrs(patch_id), patch_id, xy_local)
1599
1600 ! Reading STL boundary vertices and compute the levelset and levelset_norm
1601 if (eta > threshold) then
1602 ib_markers%sf(i, j, 0) = encoded_patch_id
1603 end if
1604 end do
1605 end do
1606
1607# 806 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1608#if defined(MFC_OpenACC)
1609# 806 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1610!$acc end parallel loop
1611# 806 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1612#elif defined(MFC_OpenMP)
1613# 806 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1614
1615# 806 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1616!$omp end target teams loop
1617# 806 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1618#endif
1619
1620 end subroutine s_ib_model
1621
1622 !> The STL patch is a 3D geometry that is imported from an STL file.
1623 subroutine s_ib_3d_model(patch_id, ib_markers, xp, yp, zp)
1624
1625 integer, intent(in) :: patch_id
1626 type(integer_field), intent(inout) :: ib_markers
1627 integer, intent(in) :: xp, yp, zp !< integers containing the periodicity projection information
1628 integer :: i, j, k, il, ir, jl, jr, kl, kr !< Generic loop iterators
1629 integer :: spc, encoded_patch_id
1630 real(wp) :: eta, threshold
1631 real(wp), dimension(1:3) :: offset
1632 real(wp), dimension(1:3) :: center, xyz_local
1633 real(wp), dimension(1:3,1:3) :: inverse_rotation, rotation
1634 integer :: cx, cy, cz
1635 real(wp) :: lx(2), ly(2), lz(2)
1636 real(wp), dimension(1:3) :: bbox_min, bbox_max, local_corner, world_corner
1637
1638 center = 0._wp
1639 center(1) = patch_ib(patch_id)%x_centroid + real(xp, wp)*(x_domain%end - x_domain%beg)
1640 center(2) = patch_ib(patch_id)%y_centroid + real(yp, wp)*(y_domain%end - y_domain%beg)
1641 center(3) = patch_ib(patch_id)%z_centroid + real(zp, wp)*(z_domain%end - z_domain%beg)
1642 inverse_rotation(:,:) = patch_ib(patch_id)%rotation_matrix_inverse(:,:)
1643 offset(:) = patch_ib(patch_id)%centroid_offset(:)
1644 spc = patch_ib(patch_id)%model_spc
1645 threshold = patch_ib(patch_id)%model_threshold
1646 rotation(:,:) = patch_ib(patch_id)%rotation_matrix(:,:)
1647
1648 ! encode the periodicity information into the patch_id
1649 call s_encode_patch_periodicity(patch_id, xp, yp, zp, encoded_patch_id)
1650
1651 il = -gp_layers - 1
1652 jl = -gp_layers - 1
1653 kl = -gp_layers - 1
1654 ir = m + gp_layers + 1
1655 jr = n + gp_layers + 1
1656 kr = p + gp_layers + 1
1657
1658 ! Local-space bounding box extents (min=1, max=2 in the third index)
1659 lx(1) = stl_bounding_boxes(patch_id, 1, 1) + offset(1)
1660 lx(2) = stl_bounding_boxes(patch_id, 1, 3) + offset(1)
1661 ly(1) = stl_bounding_boxes(patch_id, 2, 1) + offset(2)
1662 ly(2) = stl_bounding_boxes(patch_id, 2, 3) + offset(2)
1663 lz(1) = stl_bounding_boxes(patch_id, 3, 1) + offset(3)
1664 lz(2) = stl_bounding_boxes(patch_id, 3, 3) + offset(3)
1665
1666 bbox_min = 1e12
1667 bbox_max = -1e12
1668 ! Enumerate all 8 corners of the local bounding box, rotate to world space, track world-space AABB
1669 do cx = 1, 2
1670 do cy = 1, 2
1671 do cz = 1, 2
1672 local_corner = [lx(cx), ly(cy), lz(cz)]
1673 world_corner = matmul(rotation, local_corner) + center
1674 bbox_min(1) = min(bbox_min(1), world_corner(1))
1675 bbox_min(2) = min(bbox_min(2), world_corner(2))
1676 bbox_min(3) = min(bbox_min(3), world_corner(3))
1677 bbox_max(1) = max(bbox_max(1), world_corner(1))
1678 bbox_max(2) = max(bbox_max(2), world_corner(2))
1679 bbox_max(3) = max(bbox_max(3), world_corner(3))
1680 end do
1681 end do
1682 end do
1683
1684 call get_bounding_indices(bbox_min(1), bbox_max(1), x_cc, il, ir)
1685 call get_bounding_indices(bbox_min(2), bbox_max(2), y_cc, jl, jr)
1686 call get_bounding_indices(bbox_min(3), bbox_max(3), z_cc, kl, kr)
1687
1688
1689# 876 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1690
1691# 876 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1692#if defined(MFC_OpenACC)
1693# 876 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1694!$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)
1695# 876 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1696#elif defined(MFC_OpenMP)
1697# 876 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1698
1699# 876 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1700
1701# 876 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1702
1703# 876 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1704!$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)
1705# 876 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1706#endif
1707# 878 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1708 do i = il, ir
1709 do j = jl, jr
1710 do k = kl, kr
1711 xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), z_cc(k) - center(3)]
1712 xyz_local = matmul(inverse_rotation, xyz_local)
1713 xyz_local = xyz_local - offset
1714
1715 eta = f_model_is_inside_flat(gpu_ntrs(patch_id), patch_id, xyz_local)
1716
1717 if (eta > patch_ib(patch_id)%model_threshold) then
1718 ib_markers%sf(i, j, k) = encoded_patch_id
1719 end if
1720 end do
1721 end do
1722 end do
1723
1724# 893 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1725#if defined(MFC_OpenACC)
1726# 893 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1727!$acc end parallel loop
1728# 893 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1729#elif defined(MFC_OpenMP)
1730# 893 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1731
1732# 893 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1733!$omp end target teams loop
1734# 893 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1735#endif
1736
1737 end subroutine s_ib_3d_model
1738
1739 !> Compute a rotation matrix for converting to the rotating frame of the boundary
1740 subroutine s_update_ib_rotation_matrix(patch_id)
1741
1742 integer, intent(in) :: patch_id
1743 real(wp), dimension(3, 3, 3) :: rotation
1744 real(wp) :: angle
1745
1746 ! construct the x, y, and z rotation matrices
1747
1748 if (num_dims == 3) then
1749 ! also compute the x and y axes in 3D
1750 angle = patch_ib(patch_id)%angles(1)
1751 rotation(1, 1,:) = [1._wp, 0._wp, 0._wp]
1752 rotation(1, 2,:) = [0._wp, cos(angle), -sin(angle)]
1753 rotation(1, 3,:) = [0._wp, sin(angle), cos(angle)]
1754
1755 angle = patch_ib(patch_id)%angles(2)
1756 rotation(2, 1,:) = [cos(angle), 0._wp, sin(angle)]
1757 rotation(2, 2,:) = [0._wp, 1._wp, 0._wp]
1758 rotation(2, 3,:) = [-sin(angle), 0._wp, cos(angle)]
1759
1760 ! apply the y rotation to the x rotation
1761 patch_ib(patch_id)%rotation_matrix(:,:) = matmul(rotation(1,:,:), rotation(2,:,:))
1762 patch_ib(patch_id)%rotation_matrix_inverse(:,:) = matmul(transpose(rotation(2,:,:)), transpose(rotation(1,:,:)))
1763 end if
1764
1765 ! z component first, since it applies in 2D and 3D
1766 angle = patch_ib(patch_id)%angles(3)
1767 rotation(3, 1,:) = [cos(angle), -sin(angle), 0._wp]
1768 rotation(3, 2,:) = [sin(angle), cos(angle), 0._wp]
1769 rotation(3, 3,:) = [0._wp, 0._wp, 1._wp]
1770
1771 if (num_dims == 3) then
1772 ! apply the z rotation to the xy rotation in 3D
1773 patch_ib(patch_id)%rotation_matrix(:,:) = matmul(patch_ib(patch_id)%rotation_matrix(:,:), rotation(3,:,:))
1774 patch_ib(patch_id)%rotation_matrix_inverse(:,:) = matmul(transpose(rotation(3,:,:)), &
1775 & patch_ib(patch_id)%rotation_matrix_inverse(:,:))
1776 else
1777 ! write out only the z rotation in 2D
1778 patch_ib(patch_id)%rotation_matrix(:,:) = rotation(3,:,:)
1779 patch_ib(patch_id)%rotation_matrix_inverse(:,:) = transpose(rotation(3,:,:))
1780 end if
1781
1782 end subroutine s_update_ib_rotation_matrix
1783
1784 subroutine get_bounding_indices(left_bound, right_bound, cell_centers, left_index, right_index)
1785
1786 real(wp), intent(in) :: left_bound, right_bound
1787 integer, intent(inout) :: left_index, right_index
1788 real(wp), dimension(-buff_size:), intent(in) :: cell_centers
1789 integer :: itr_left, itr_middle, itr_right
1790
1791 itr_left = left_index
1792 itr_right = right_index
1793
1794 do while (itr_left + 1 < itr_right)
1795 itr_middle = (itr_left + itr_right)/2
1796 if (cell_centers(itr_middle) < left_bound) then
1797 itr_left = itr_middle
1798 else if (cell_centers(itr_middle) > left_bound) then
1799 itr_right = itr_middle
1800 else
1801 itr_left = itr_middle
1802 exit
1803 end if
1804 end do
1805 left_index = itr_left
1806
1807 itr_right = right_index
1808 do while (itr_left + 1 < itr_right)
1809 itr_middle = (itr_left + itr_right)/2
1810 if (cell_centers(itr_middle) < right_bound) then
1811 itr_left = itr_middle
1812 else if (cell_centers(itr_middle) > right_bound) then
1813 itr_right = itr_middle
1814 else
1815 itr_right = itr_middle
1816 exit
1817 end if
1818 end do
1819 right_index = itr_right
1820
1821 end subroutine get_bounding_indices
1822
1823 !> Encode the patch ID with a unique offset containing periodicity information
1824 subroutine s_encode_patch_periodicity(patch_id, x_periodicity, y_periodicity, z_periodicity, encoded_patch_id)
1825
1826 integer, intent(in) :: patch_id, x_periodicity, y_periodicity, z_periodicity
1827 integer, intent(out) :: encoded_patch_id
1828 integer :: temp_x_per, temp_y_per, temp_z_per, offset
1829
1830 encoded_patch_id = patch_id
1831
1832 temp_x_per = x_periodicity; if (x_periodicity == -1) temp_x_per = 2
1833 temp_y_per = y_periodicity; if (y_periodicity == -1) temp_y_per = 2
1834 temp_z_per = z_periodicity; if (z_periodicity == -1) temp_z_per = 2
1835
1836 offset = (num_ibs + 1)*temp_x_per + 3*(num_ibs + 1)*temp_y_per + 9*(num_ibs + 1)*temp_z_per
1837 encoded_patch_id = patch_id + offset
1838
1839 end subroutine s_encode_patch_periodicity
1840
1841 !> Decode the encoded ID to recover the original patch ID and periodicity
1842 subroutine s_decode_patch_periodicity(encoded_patch_id, patch_id, x_periodicity, y_periodicity, z_periodicity)
1843
1844
1845# 1002 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1846#if MFC_OpenACC
1847# 1002 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1848!$acc routine seq
1849# 1002 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1850#elif MFC_OpenMP
1851# 1002 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1852
1853# 1002 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1854
1855# 1002 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1856!$omp declare target device_type(any)
1857# 1002 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1858#endif
1859
1860 integer, intent(in) :: encoded_patch_id
1861 integer, intent(out) :: patch_id
1862 integer, intent(out), optional :: x_periodicity, y_periodicity, z_periodicity
1863 integer :: offset, remainder, xp, yp, zp, base
1864
1865 base = num_ibs + 1
1866
1867 patch_id = mod(encoded_patch_id - 1, base) + 1
1868 offset = (encoded_patch_id - patch_id)/base
1869
1870 xp = mod(offset, 3)
1871 remainder = offset/3
1872 yp = mod(remainder, 3)
1873 zp = remainder/3
1874
1875 ! Reverse map: 2 -> -1, 0 -> 0, 1 -> 1
1876 if (present(x_periodicity) .and. present(y_periodicity) .and. present(z_periodicity)) then
1877 x_periodicity = xp; if (xp == 2) x_periodicity = -1
1878 y_periodicity = yp; if (yp == 2) y_periodicity = -1
1879 z_periodicity = zp; if (zp == 2) z_periodicity = -1
1880 end if
1881
1882 end subroutine s_decode_patch_periodicity
1883
1884 !> Determine the periodic wrapping bounds in each direction
1885 subroutine s_get_periodicities(xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper)
1886
1887 integer, intent(out) :: xp_lower, xp_upper, yp_lower, yp_upper
1888 integer, intent(out), optional :: zp_lower, zp_upper
1889
1890 ! check domain wraps in x, y
1891
1892# 1037 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1893 ! check for periodicity
1894 if (ib_bc_x%beg == bc_periodic) then
1895 xp_lower = -1
1896 xp_upper = 1
1897 else
1898 ! if it is not periodic, then both elements are 0
1899 xp_lower = 0
1900 xp_upper = 0
1901 end if
1902# 1037 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1903 ! check for periodicity
1904 if (ib_bc_y%beg == bc_periodic) then
1905 yp_lower = -1
1906 yp_upper = 1
1907 else
1908 ! if it is not periodic, then both elements are 0
1909 yp_lower = 0
1910 yp_upper = 0
1911 end if
1912# 1047 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1913
1914 ! z only if 3D
1915 if (present(zp_lower) .and. p /= 0) then
1916 if (ib_bc_z%beg == bc_periodic) then
1917 zp_lower = -1
1918 zp_upper = 1
1919 else
1920 zp_lower = 0
1921 zp_upper = 0
1922 end if
1923 end if
1924
1925 end subroutine s_get_periodicities
1926
1927 !> Archimedes spiral function
1928 pure elemental function f_r(myth, offset, a)
1929
1930
1931# 1064 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1932#if MFC_OpenACC
1933# 1064 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1934!$acc routine seq
1935# 1064 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1936#elif MFC_OpenMP
1937# 1064 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1938
1939# 1064 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1940
1941# 1064 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1942!$omp declare target device_type(any)
1943# 1064 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1944#endif
1945 real(wp), intent(in) :: myth, offset, a
1946 real(wp) :: b
1947 real(wp) :: f_r
1948
1949 ! r(th) = a + b*th
1950
1951 b = 2._wp*a/(2._wp*pi)
1952 f_r = a + b*myth + offset
1953
1954 end function f_r
1955
1956end 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:774
m_global_parameters::airfoil_grid_u
type(vec3_dt), dimension(:), allocatable airfoil_grid_u
Definition m_global_parameters.fpp.f90:1010
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:774
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:1011
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:1003
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:1009
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:1340
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:1741
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:537
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:1624
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:465
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:1886
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:1435
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:1785
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:806
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:1517
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:1825
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:404
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:1076
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:1843
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:1159
m_ib_patches::f_r
pure elemental real(wp) function f_r(myth, offset, a)
Archimedes spiral function.
Definition m_ib_patches.fpp.f90:1929
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:1249
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