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