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