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# 76 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
303
304# 91 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
305
306# 102 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
307
308# 115 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
309
310# 143 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
311
312# 154 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
313
314# 165 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
315
316# 176 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
317
318# 187 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
319
320# 198 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
321
322# 208 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
323
324# 214 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
325
326# 220 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
327
328# 226 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
329
330# 232 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
331
332# 234 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
333# 235 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
334! New line at end of file is required for FYPP
335# 2 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 2
336
337# 14 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
338
339! Caution:
340! This macro requires the use of a binding script to set CUDA_VISIBLE_DEVICES, such that we have one GPU device per MPI rank.
341! That's because for both cudaMemAdvise (preferred location) and cudaMemPrefetchAsync we use location = device_id = 0.
342! For an example see misc/nvidia_uvm/bind.sh.
343# 63 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
344
345# 81 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
346
347# 88 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
348
349# 111 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
350
351# 127 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
352
353# 153 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
354
355# 159 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
356
357# 167 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
358! New line at end of file is required for FYPP
359# 11 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp" 2
360
361!> @brief Immersed boundary patch geometry constructors for 2D and 3D shapes
362module m_ib_patches
363
364 use m_model ! Subroutine(s) related to STL files
365
366 use m_derived_types ! Definitions of the derived types
367
368 use m_global_parameters !< definitions of the global parameters
369
370 use m_helper_basic !< functions to compare floating point numbers
371
372 use m_helper
373
374 use m_mpi_common
375
376 implicit none
377
378 private; public :: s_apply_ib_patches, s_update_ib_rotation_matrix, s_instantiate_stl_models, models, f_convert_cyl_to_cart
379
380 real(wp) :: x_centroid, y_centroid, z_centroid
381 real(wp) :: length_x, length_y, length_z
382
383# 33 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
384#if defined(MFC_OpenACC)
385# 33 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
386!$acc declare create(x_centroid, y_centroid, z_centroid)
387# 33 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
388#elif defined(MFC_OpenMP)
389# 33 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
390!$omp declare target (x_centroid, y_centroid, z_centroid)
391# 33 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
392#endif
393
394# 34 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
395#if defined(MFC_OpenACC)
396# 34 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
397!$acc declare create(length_x, length_y, length_z)
398# 34 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
399#elif defined(MFC_OpenMP)
400# 34 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
401!$omp declare target (length_x, length_y, length_z)
402# 34 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
403#endif
404
405 integer :: smooth_patch_id
406 real(wp) :: smooth_coeff
407
408# 38 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
409#if defined(MFC_OpenACC)
410# 38 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
411!$acc declare create(smooth_patch_id, smooth_coeff)
412# 38 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
413#elif defined(MFC_OpenMP)
414# 38 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
415!$omp declare target (smooth_patch_id, smooth_coeff)
416# 38 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
417#endif
418 !! These variables are analogous in both meaning and use to the similarly
419 !! named components in the ic_patch_parameters type (see m_derived_types.f90
420 !! for additional details). They are employed as a means to more concisely
421 !! perform the actions necessary to lay out a particular patch on the grid.
422
423 real(wp) :: cart_x, cart_y, cart_z
424 real(wp) :: sph_phi !<
425
426# 46 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
427#if defined(MFC_OpenACC)
428# 46 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
429!$acc declare create(cart_x, cart_y, cart_z, sph_phi)
430# 46 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
431#elif defined(MFC_OpenMP)
432# 46 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
433!$omp declare target (cart_x, cart_y, cart_z, sph_phi)
434# 46 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
435#endif
436 !! Variables to be used to hold cell locations in Cartesian coordinates if
437 !! 3D simulation is using cylindrical coordinates
438
439 type(bounds_info) :: x_boundary, y_boundary, z_boundary
440
441# 51 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
442#if defined(MFC_OpenACC)
443# 51 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
444!$acc declare create(x_boundary, y_boundary, z_boundary)
445# 51 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
446#elif defined(MFC_OpenMP)
447# 51 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
448!$omp declare target (x_boundary, y_boundary, z_boundary)
449# 51 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
450#endif
451 !! These variables combine the centroid and length parameters associated with
452 !! a particular patch to yield the locations of the patch boundaries in the
453 !! x-, y- and z-coordinate directions. They are used as a means to concisely
454 !! perform the actions necessary to lay out a particular patch on the grid.
455
456 character(len=5) :: istr ! string to store int to string result for error checking
457
458 type(t_model_array), allocatable, target :: models(:)
459
460# 60 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
461#if defined(MFC_OpenACC)
462# 60 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
463!$acc declare create(models)
464# 60 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
465#elif defined(MFC_OpenMP)
466# 60 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
467!$omp declare target (models)
468# 60 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
469#endif
470 !! array of STL models that can be allocated and then used in IB marker and levelset compute
471
472contains
473
474 !> @brief Applies all immersed boundary patch geometries to mark interior cells in the IB marker array.
475 impure subroutine s_apply_ib_patches(ib_markers_sf)
476
477 integer, dimension(:, :, :), intent(inout), optional :: ib_markers_sf
478
479 integer :: i
480
481 ! 3D Patch Geometries
482 if (p > 0) then
483
484 !> IB Patches
485 !> @{
486 ! Spherical patch
487 do i = 1, num_ibs
488
489 if (patch_ib(i)%geometry == 8) then
490 call s_ib_sphere(i, ib_markers_sf)
491 elseif (patch_ib(i)%geometry == 9) then
492 call s_ib_cuboid(i, ib_markers_sf)
493 elseif (patch_ib(i)%geometry == 10) then
494 call s_ib_cylinder(i, ib_markers_sf)
495 elseif (patch_ib(i)%geometry == 11) then
496 call s_ib_3d_airfoil(i, ib_markers_sf)
497 ! STL+IBM patch
498 elseif (patch_ib(i)%geometry == 12) then
499 call s_ib_model(i, ib_markers_sf)
500 end if
501 end do
502 !> @}
503
504 ! 2D Patch Geometries
505 elseif (n > 0) then
506
507 !> IB Patches
508 !> @{
509 do i = 1, num_ibs
510 if (patch_ib(i)%geometry == 2) then
511 call s_ib_circle(i, ib_markers_sf)
512 elseif (patch_ib(i)%geometry == 3) then
513 call s_ib_rectangle(i, ib_markers_sf)
514 elseif (patch_ib(i)%geometry == 4) then
515 call s_ib_airfoil(i, ib_markers_sf)
516 ! STL+IBM patch
517 elseif (patch_ib(i)%geometry == 5) then
518 call s_ib_model(i, ib_markers_sf)
519 elseif (patch_ib(i)%geometry == 6) then
520 call s_ib_ellipse(i, ib_markers_sf)
521 end if
522 end do
523 !> @}
524
525 end if
526
527 end subroutine s_apply_ib_patches
528
529 subroutine s_instantiate_stl_models()
530
531 ! Variables for IBM+STL
532 real(wp) :: normals(1:3) !< Boundary normal buffer
533 integer :: boundary_vertex_count, boundary_edge_count, total_vertices !< Boundary vertex
534 real(wp), allocatable, dimension(:, :, :) :: boundary_v !< Boundary vertex buffer
535 real(wp), allocatable, dimension(:, :) :: interpolated_boundary_v !< Interpolated vertex buffer
536 real(wp) :: distance !< Levelset distance buffer
537 logical :: interpolate !< Logical variable to determine whether or not the model should be interpolated
538
539 integer :: i, j, k !< Generic loop iterators
540 integer :: patch_id
541
542 type(t_bbox) :: bbox, bbox_old
543 type(t_model) :: model
544 type(ic_model_parameters) :: params
545
546 real(wp) :: eta
547 real(wp), dimension(1:3) :: point, model_center
548 real(wp) :: grid_mm(1:3, 1:2)
549
550 real(wp), dimension(1:4, 1:4) :: transform, transform_n
551
552 call random_seed()
553
554 do patch_id = 1, num_ibs
555 if (patch_ib(patch_id)%geometry == 5 .or. patch_ib(patch_id)%geometry == 12) then
556 allocate (models(patch_id)%model)
557 print *, " * Reading model: "//trim(patch_ib(patch_id)%model_filepath)
558
559 model = f_model_read(patch_ib(patch_id)%model_filepath)
560 params%scale(:) = patch_ib(patch_id)%model_scale(:)
561 params%translate(:) = patch_ib(patch_id)%model_translate(:)
562 params%rotate(:) = patch_ib(patch_id)%model_rotate(:)
563 params%spc = patch_ib(patch_id)%model_spc
564 params%threshold = patch_ib(patch_id)%model_threshold
565
566 if (f_approx_equal(dot_product(params%scale, params%scale), 0._wp)) then
567 params%scale(:) = 1._wp
568 end if
569
570 if (proc_rank == 0) then
571 print *, " * Transforming model."
572 end if
573
574 ! Get the model center before transforming the model
575 bbox_old = f_create_bbox(model)
576 model_center(1:3) = (bbox_old%min(1:3) + bbox_old%max(1:3))/2._wp
577
578 ! Compute the transform matrices for vertices and normals
579 transform = f_create_transform_matrix(params, model_center)
580 transform_n = f_create_transform_matrix(params)
581
582 call s_transform_model(model, transform, transform_n)
583
584 ! Recreate the bounding box after transformation
585 bbox = f_create_bbox(model)
586
587 ! Show the number of vertices in the original STL model
588 if (proc_rank == 0) then
589 print *, ' * Number of input model vertices:', 3*model%ntrs
590 end if
591
592 call f_check_boundary(model, boundary_v, boundary_vertex_count, boundary_edge_count)
593
594 ! Check if the model needs interpolation
595 if (p > 0) then
596 call f_check_interpolation_3d(model, (/dx, dy, dz/), interpolate)
597 else
598 call f_check_interpolation_2d(boundary_v, boundary_edge_count, (/minval(dx), minval(dy), 0._wp/), interpolate)
599 end if
600
601 ! Show the number of edges and boundary edges in 2D STL models
602 if (proc_rank == 0 .and. p == 0) then
603 print *, ' * Number of 2D model boundary edges:', boundary_edge_count
604 end if
605
606 ! Interpolate the STL model along the edges (2D) and on triangle facets (3D)
607 if (interpolate) then
608 if (proc_rank == 0) then
609 print *, ' * Interpolating STL vertices.'
610 end if
611
612 if (p > 0) then
613 call f_interpolate_3d(model, (/dx, dy, dz/), interpolated_boundary_v, total_vertices)
614 else
615 call f_interpolate_2d(boundary_v, boundary_edge_count, (/dx, dy, dz/), interpolated_boundary_v, total_vertices)
616 end if
617
618 if (proc_rank == 0) then
619 print *, ' * Total number of interpolated boundary vertices:', total_vertices
620 end if
621 end if
622
623 if (proc_rank == 0) then
624 write (*, "(A, 3(2X, F20.10))") " > Model: Min:", bbox%min(1:3)
625 write (*, "(A, 3(2X, F20.10))") " > Cen:", (bbox%min(1:3) + bbox%max(1:3))/2._wp
626 write (*, "(A, 3(2X, F20.10))") " > Max:", bbox%max(1:3)
627
628 grid_mm(1, :) = (/minval(x_cc(0:m)) - 0.e5_wp*dx(0), maxval(x_cc(0:m)) + 0.e5_wp*dx(m)/)
629 grid_mm(2, :) = (/minval(y_cc(0:n)) - 0.e5_wp*dy(0), maxval(y_cc(0:n)) + 0.e5_wp*dy(n)/)
630
631 if (p > 0) then
632 grid_mm(3, :) = (/minval(z_cc(0:p)) - 0.e5_wp*dz(0), maxval(z_cc(0:p)) + 0.e5_wp*dz(p)/)
633 else
634 grid_mm(3, :) = (/0._wp, 0._wp/)
635 end if
636
637 write (*, "(A, 3(2X, F20.10))") " > Domain: Min:", grid_mm(:, 1)
638 write (*, "(A, 3(2X, F20.10))") " > Cen:", (grid_mm(:, 1) + grid_mm(:, 2))/2._wp
639 write (*, "(A, 3(2X, F20.10))") " > Max:", grid_mm(:, 2)
640 end if
641
642 models(patch_id)%model = model
643 models(patch_id)%boundary_v = boundary_v
644 models(patch_id)%boundary_edge_count = boundary_edge_count
645 models(patch_id)%interpolate = interpolate
646 if (interpolate) then
647 models(patch_id)%interpolated_boundary_v = interpolated_boundary_v
648 models(patch_id)%total_vertices = total_vertices
649 end if
650
651 end if
652 end do
653
654 end subroutine s_instantiate_stl_models
655
656 !> The circular patch is a 2D geometry that may be used, for
657 !! example, in creating a bubble or a droplet. The geometry
658 !! of the patch is well-defined when its centroid and radius
659 !! are provided. Note that the circular patch DOES allow for
660 !! the smoothing of its boundary.
661 !! @param patch_id is the patch identifier
662 !! @param ib_markers_sf Array to track patch ids
663 subroutine s_ib_circle(patch_id, ib_markers_sf)
664
665 integer, intent(in) :: patch_id
666 integer, dimension(0:m, 0:n, 0:p), intent(inout) :: ib_markers_sf
667
668 real(wp), dimension(1:2) :: center
669 real(wp) :: radius
670
671 integer :: i, j, k !< Generic loop iterators
672
673 ! Transferring the circular patch's radius, centroid, smearing patch
674 ! identity and smearing coefficient information
675
676 center(1) = patch_ib(patch_id)%x_centroid
677 center(2) = patch_ib(patch_id)%y_centroid
678 radius = patch_ib(patch_id)%radius
679
680 ! Checking whether the circle covers a particular cell in the domain
681 ! and verifying whether the current patch has permission to write to
682 ! that cell. If both queries check out, the primitive variables of
683 ! the current patch are assigned to this cell.
684
685
686# 276 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
687
688# 276 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
689#if defined(MFC_OpenACC)
690# 276 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
691!$acc parallel loop collapse(2) gang vector default(present) private(i, j) copy(ib_markers_sf) copyin(patch_id, center, radius)
692# 276 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
693#elif defined(MFC_OpenMP)
694# 276 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
695
696# 276 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
697
698# 276 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
699
700# 276 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
701!$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(tofrom:ib_markers_sf) map(to:patch_id, center, radius)
702# 276 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
703#endif
704# 276 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
705
706# 278 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
707 do j = 0, n
708 do i = 0, m
709 if ((x_cc(i) - center(1))**2 &
710 + (y_cc(j) - center(2))**2 <= radius**2) &
711 then
712 ib_markers_sf(i, j, 0) = patch_id
713 end if
714 end do
715 end do
716
717# 287 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
718
719# 287 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
720#if defined(MFC_OpenACC)
721# 287 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
722!$acc end parallel loop
723# 287 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
724#elif defined(MFC_OpenMP)
725# 287 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
726
727# 287 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
728
729# 287 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
730!$omp end target teams loop
731# 287 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
732#endif
733# 287 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
734
735
736 end subroutine s_ib_circle
737
738 !> @brief Marks cells inside a 2D NACA 4-digit airfoil immersed boundary using upper and lower surface grids.
739 !! @param patch_id is the patch identifier
740 !! @param ib_markers_sf Array to track patch ids
741 subroutine s_ib_airfoil(patch_id, ib_markers_sf)
742
743 integer, intent(in) :: patch_id
744 integer, dimension(0:m, 0:n, 0:p), intent(inout) :: ib_markers_sf
745
746 real(wp) :: f, ca_in, pa, ma, ta
747 real(wp) :: xa, yt, xu, yu, xl, yl, xc, yc, dycdxc, sin_c, cos_c
748 integer :: i, j, k
749 integer :: Np1, Np2
750
751 real(wp), dimension(1:3) :: xy_local, offset !< x and y coordinates in local IB frame
752 real(wp), dimension(1:2) :: center !< x and y coordinates in local IB frame
753 real(wp), dimension(1:3, 1:3) :: inverse_rotation
754
755 center(1) = patch_ib(patch_id)%x_centroid
756 center(2) = patch_ib(patch_id)%y_centroid
757 ca_in = patch_ib(patch_id)%c
758 pa = patch_ib(patch_id)%p
759 ma = patch_ib(patch_id)%m
760 ta = patch_ib(patch_id)%t
761 inverse_rotation(:, :) = patch_ib(patch_id)%rotation_matrix_inverse(:, :)
762 offset(:) = patch_ib(patch_id)%centroid_offset(:)
763
764 np1 = int((pa*ca_in/dx(0))*20)
765 np2 = int(((ca_in - pa*ca_in)/dx(0))*20)
766 np = np1 + np2 + 1
767
768# 320 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
769#if defined(MFC_OpenACC)
770# 320 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
771!$acc update device(Np)
772# 320 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
773#elif defined(MFC_OpenMP)
774# 320 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
775!$omp target update to(Np)
776# 320 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
777#endif
778
779 if (.not. allocated(airfoil_grid_u)) then
780#ifdef MFC_DEBUG
781# 323 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
782 block
783# 323 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
784 use iso_fortran_env, only: output_unit
785# 323 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
786
787# 323 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
788 print *, 'm_ib_patches.fpp:323: ', '@:ALLOCATE(airfoil_grid_u(1:Np))'
789# 323 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
790
791# 323 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
792 call flush (output_unit)
793# 323 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
794 end block
795# 323 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
796#endif
797# 323 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
798 allocate (airfoil_grid_u(1:np))
799# 323 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
800
801# 323 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
802
803# 323 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
804#if defined(MFC_OpenACC)
805# 323 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
806!$acc enter data create(airfoil_grid_u)
807# 323 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
808#elif defined(MFC_OpenMP)
809# 323 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
810!$omp target enter data map(always,alloc:airfoil_grid_u)
811# 323 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
812#endif
813#ifdef MFC_DEBUG
814# 324 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
815 block
816# 324 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
817 use iso_fortran_env, only: output_unit
818# 324 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
819
820# 324 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
821 print *, 'm_ib_patches.fpp:324: ', '@:ALLOCATE(airfoil_grid_l(1:Np))'
822# 324 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
823
824# 324 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
825 call flush (output_unit)
826# 324 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
827 end block
828# 324 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
829#endif
830# 324 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
831 allocate (airfoil_grid_l(1:np))
832# 324 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
833
834# 324 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
835
836# 324 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
837#if defined(MFC_OpenACC)
838# 324 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
839!$acc enter data create(airfoil_grid_l)
840# 324 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
841#elif defined(MFC_OpenMP)
842# 324 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
843!$omp target enter data map(always,alloc:airfoil_grid_l)
844# 324 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
845#endif
846
847 ! TODO :: The below instantiations are already handles by the loop below
848 airfoil_grid_u(1)%x = 0._wp
849 airfoil_grid_u(1)%y = 0._wp
850
851 airfoil_grid_l(1)%x = 0._wp
852 airfoil_grid_l(1)%y = 0._wp
853
854 do i = 1, np1 + np2 - 1
855 ! 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.
856 if (i <= np1) then
857 xc = i*(pa*ca_in/np1)
858 xa = xc/ca_in
859 yc = (ma/pa**2)*(2*pa*xa - xa**2)
860 dycdxc = (2*ma/pa**2)*(pa - xa)
861 else
862 xc = pa*ca_in + (i - np1)*((ca_in - pa*ca_in)/np2)
863 xa = xc/ca_in
864 yc = (ma/(1 - pa)**2)*(1 - 2*pa + 2*pa*xa - xa**2)
865 dycdxc = (2*ma/(1 - pa)**2)*(pa - xa)
866 end if
867
868 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)
869 sin_c = dycdxc/(1 + dycdxc**2)**0.5_wp
870 cos_c = 1/(1 + dycdxc**2)**0.5_wp
871
872 xu = xa - yt*sin_c
873 yu = yc + yt*cos_c
874
875 xl = xa + yt*sin_c
876 yl = yc - yt*cos_c
877
878 xu = xu*ca_in
879 yu = yu*ca_in
880
881 xl = xl*ca_in
882 yl = yl*ca_in
883
884 airfoil_grid_u(i + 1)%x = xu
885 airfoil_grid_u(i + 1)%y = yu
886
887 airfoil_grid_l(i + 1)%x = xl
888 airfoil_grid_l(i + 1)%y = yl
889
890 end do
891
892 airfoil_grid_u(np)%x = ca_in
893 airfoil_grid_u(np)%y = 0._wp
894
895 airfoil_grid_l(np)%x = ca_in
896 airfoil_grid_l(np)%y = 0._wp
897
898
899# 377 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
900#if defined(MFC_OpenACC)
901# 377 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
902!$acc update device(airfoil_grid_l, airfoil_grid_u)
903# 377 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
904#elif defined(MFC_OpenMP)
905# 377 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
906!$omp target update to(airfoil_grid_l, airfoil_grid_u)
907# 377 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
908#endif
909
910 end if
911
912
913# 381 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
914
915# 381 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
916#if defined(MFC_OpenACC)
917# 381 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
918!$acc parallel loop collapse(2) gang vector default(present) private(i, j, xy_local, k, f) copy(ib_markers_sf) copyin(patch_id, center, inverse_rotation, offset, ma, ca_in, airfoil_grid_u, airfoil_grid_l)
919# 381 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
920#elif defined(MFC_OpenMP)
921# 381 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
922
923# 381 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
924
925# 381 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
926
927# 381 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
928!$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(tofrom:ib_markers_sf) map(to:patch_id, center, inverse_rotation, offset, ma, ca_in, airfoil_grid_u, airfoil_grid_l)
929# 381 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
930#endif
931# 381 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
932
933# 383 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
934 do j = 0, n
935 do i = 0, m
936 xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp] ! get coordinate frame centered on IB
937 xy_local = matmul(inverse_rotation, xy_local) ! rotate the frame into the IB's coordinates
938 xy_local = xy_local - offset ! airfoils are a patch that require a centroid offset
939
940 if (xy_local(1) >= 0._wp .and. xy_local(1) <= ca_in) then
941 xa = xy_local(1)/ca_in
942 if (xa <= pa) then
943 yc = (ma/pa**2)*(2*pa*xa - xa**2)
944 dycdxc = (2*ma/pa**2)*(pa - xa)
945 else
946 yc = (ma/(1 - pa)**2)*(1 - 2*pa + 2*pa*xa - xa**2)
947 dycdxc = (2*ma/(1 - pa)**2)*(pa - xa)
948 end if
949 if (xy_local(2) >= 0._wp) then
950 k = 1
951 do while (airfoil_grid_u(k)%x < xy_local(1) .and. k <= np)
952 k = k + 1
953 end do
954 if (f_approx_equal(airfoil_grid_u(k)%x, xy_local(1))) then
955 if (xy_local(2) <= airfoil_grid_u(k)%y) then
956 !!IB
957 ib_markers_sf(i, j, 0) = patch_id
958 end if
959 else
960 f = (airfoil_grid_u(k)%x - xy_local(1))/(airfoil_grid_u(k)%x - airfoil_grid_u(k - 1)%x)
961 if (xy_local(2) <= ((1._wp - f)*airfoil_grid_u(k)%y + f*airfoil_grid_u(k - 1)%y)) then
962 !!IB
963 ib_markers_sf(i, j, 0) = patch_id
964 end if
965 end if
966 else
967 k = 1
968 do while (airfoil_grid_l(k)%x < xy_local(1))
969 k = k + 1
970 end do
971 if (f_approx_equal(airfoil_grid_l(k)%x, xy_local(1))) then
972 if (xy_local(2) >= airfoil_grid_l(k)%y) then
973 !!IB
974 ib_markers_sf(i, j, 0) = patch_id
975 end if
976 else
977 f = (airfoil_grid_l(k)%x - xy_local(1))/(airfoil_grid_l(k)%x - airfoil_grid_l(k - 1)%x)
978
979 if (xy_local(2) >= ((1._wp - f)*airfoil_grid_l(k)%y + f*airfoil_grid_l(k - 1)%y)) then
980 !!IB
981 ib_markers_sf(i, j, 0) = patch_id
982 end if
983 end if
984 end if
985 end if
986 end do
987 end do
988
989# 437 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
990
991# 437 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
992#if defined(MFC_OpenACC)
993# 437 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
994!$acc end parallel loop
995# 437 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
996#elif defined(MFC_OpenMP)
997# 437 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
998
999# 437 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1000
1001# 437 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1002!$omp end target teams loop
1003# 437 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1004#endif
1005# 437 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1006
1007
1008 end subroutine s_ib_airfoil
1009
1010 !> @brief Marks cells inside a 3D extruded NACA 4-digit airfoil immersed boundary with finite span.
1011 !! @param patch_id is the patch identifier
1012 !! @param ib_markers_sf Array to track patch ids
1013 subroutine s_ib_3d_airfoil(patch_id, ib_markers_sf)
1014
1015 integer, intent(in) :: patch_id
1016 integer, dimension(0:m, 0:n, 0:p), intent(inout) :: ib_markers_sf
1017
1018 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
1019 integer :: i, j, k, l
1020 integer :: Np1, Np2
1021
1022 real(wp), dimension(1:3) :: xyz_local, center, offset !< x, y, z coordinates in local IB frame
1023 real(wp), dimension(1:3, 1:3) :: inverse_rotation
1024
1025 center(1) = patch_ib(patch_id)%x_centroid
1026 center(2) = patch_ib(patch_id)%y_centroid
1027 center(3) = patch_ib(patch_id)%z_centroid
1028 lz = patch_ib(patch_id)%length_z
1029 ca_in = patch_ib(patch_id)%c
1030 pa = patch_ib(patch_id)%p
1031 ma = patch_ib(patch_id)%m
1032 ta = patch_ib(patch_id)%t
1033 inverse_rotation(:, :) = patch_ib(patch_id)%rotation_matrix_inverse(:, :)
1034 offset(:) = patch_ib(patch_id)%centroid_offset(:)
1035
1036 z_max = lz/2
1037 z_min = -lz/2
1038
1039 np1 = int((pa*ca_in/dx(0))*20)
1040 np2 = int(((ca_in - pa*ca_in)/dx(0))*20)
1041 np = np1 + np2 + 1
1042
1043# 473 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1044#if defined(MFC_OpenACC)
1045# 473 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1046!$acc update device(Np)
1047# 473 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1048#elif defined(MFC_OpenMP)
1049# 473 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1050!$omp target update to(Np)
1051# 473 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1052#endif
1053
1054 if (.not. allocated(airfoil_grid_u)) then
1055
1056#ifdef MFC_DEBUG
1057# 477 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1058 block
1059# 477 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1060 use iso_fortran_env, only: output_unit
1061# 477 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1062
1063# 477 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1064 print *, 'm_ib_patches.fpp:477: ', '@:ALLOCATE(airfoil_grid_u(1:Np))'
1065# 477 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1066
1067# 477 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1068 call flush (output_unit)
1069# 477 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1070 end block
1071# 477 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1072#endif
1073# 477 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1074 allocate (airfoil_grid_u(1:np))
1075# 477 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1076
1077# 477 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1078
1079# 477 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1080#if defined(MFC_OpenACC)
1081# 477 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1082!$acc enter data create(airfoil_grid_u)
1083# 477 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1084#elif defined(MFC_OpenMP)
1085# 477 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1086!$omp target enter data map(always,alloc:airfoil_grid_u)
1087# 477 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1088#endif
1089#ifdef MFC_DEBUG
1090# 478 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1091 block
1092# 478 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1093 use iso_fortran_env, only: output_unit
1094# 478 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1095
1096# 478 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1097 print *, 'm_ib_patches.fpp:478: ', '@:ALLOCATE(airfoil_grid_l(1:Np))'
1098# 478 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1099
1100# 478 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1101 call flush (output_unit)
1102# 478 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1103 end block
1104# 478 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1105#endif
1106# 478 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1107 allocate (airfoil_grid_l(1:np))
1108# 478 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1109
1110# 478 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1111
1112# 478 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1113#if defined(MFC_OpenACC)
1114# 478 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1115!$acc enter data create(airfoil_grid_l)
1116# 478 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1117#elif defined(MFC_OpenMP)
1118# 478 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1119!$omp target enter data map(always,alloc:airfoil_grid_l)
1120# 478 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1121#endif
1122
1123 airfoil_grid_u(1)%x = 0._wp
1124 airfoil_grid_u(1)%y = 0._wp
1125
1126 airfoil_grid_l(1)%x = 0._wp
1127 airfoil_grid_l(1)%y = 0._wp
1128
1129 do i = 1, np1 + np2 - 1
1130 if (i <= np1) then
1131 xc = i*(pa*ca_in/np1)
1132 xa = xc/ca_in
1133 yc = (ma/pa**2)*(2*pa*xa - xa**2)
1134 dycdxc = (2*ma/pa**2)*(pa - xa)
1135 else
1136 xc = pa*ca_in + (i - np1)*((ca_in - pa*ca_in)/np2)
1137 xa = xc/ca_in
1138 yc = (ma/(1 - pa)**2)*(1 - 2*pa + 2*pa*xa - xa**2)
1139 dycdxc = (2*ma/(1 - pa)**2)*(pa - xa)
1140 end if
1141
1142 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)
1143 sin_c = dycdxc/(1 + dycdxc**2)**0.5_wp
1144 cos_c = 1/(1 + dycdxc**2)**0.5_wp
1145
1146 xu = xa - yt*sin_c
1147 yu = yc + yt*cos_c
1148
1149 xl = xa + yt*sin_c
1150 yl = yc - yt*cos_c
1151
1152 xu = xu*ca_in
1153 yu = yu*ca_in
1154
1155 xl = xl*ca_in
1156 yl = yl*ca_in
1157
1158 airfoil_grid_u(i + 1)%x = xu
1159 airfoil_grid_u(i + 1)%y = yu
1160
1161 airfoil_grid_l(i + 1)%x = xl
1162 airfoil_grid_l(i + 1)%y = yl
1163
1164 end do
1165
1166 airfoil_grid_u(np)%x = ca_in
1167 airfoil_grid_u(np)%y = 0._wp
1168
1169 airfoil_grid_l(np)%x = ca_in
1170 airfoil_grid_l(np)%y = 0._wp
1171
1172
1173# 529 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1174#if defined(MFC_OpenACC)
1175# 529 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1176!$acc update device(airfoil_grid_l, airfoil_grid_u)
1177# 529 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1178#elif defined(MFC_OpenMP)
1179# 529 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1180!$omp target update to(airfoil_grid_l, airfoil_grid_u)
1181# 529 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1182#endif
1183 end if
1184
1185
1186# 532 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1187
1188# 532 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1189#if defined(MFC_OpenACC)
1190# 532 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1191!$acc parallel loop collapse(3) gang vector default(present) private(i, j, l, xyz_local, k, f) copy(ib_markers_sf) copyin(patch_id, center, inverse_rotation, offset, ma, ca_in, airfoil_grid_u, airfoil_grid_l, z_min, z_max)
1192# 532 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1193#elif defined(MFC_OpenMP)
1194# 532 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1195
1196# 532 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1197
1198# 532 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1199
1200# 532 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1201!$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(tofrom:ib_markers_sf) map(to:patch_id, center, inverse_rotation, offset, ma, ca_in, airfoil_grid_u, airfoil_grid_l, z_min, z_max)
1202# 532 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1203#endif
1204# 532 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1205
1206# 534 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1207 do l = 0, p
1208 do j = 0, n
1209 do i = 0, m
1210 xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), z_cc(l) - center(3)] ! get coordinate frame centered on IB
1211 xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinates
1212 xyz_local = xyz_local - offset ! airfoils are a patch that require a centroid offset
1213
1214 if (xyz_local(3) >= z_min .and. xyz_local(3) <= z_max) then
1215
1216 if (xyz_local(1) >= 0._wp .and. xyz_local(1) <= ca_in) then
1217 if (xyz_local(2) >= 0._wp) then
1218 k = 1
1219 do while (airfoil_grid_u(k)%x < xyz_local(1))
1220 k = k + 1
1221 end do
1222 if (f_approx_equal(airfoil_grid_u(k)%x, xyz_local(1))) then
1223 if (xyz_local(2) <= airfoil_grid_u(k)%y) then
1224 !IB
1225 ib_markers_sf(i, j, l) = patch_id
1226 end if
1227 else
1228 f = (airfoil_grid_u(k)%x - xyz_local(1))/(airfoil_grid_u(k)%x - airfoil_grid_u(k - 1)%x)
1229 if (xyz_local(2) <= ((1._wp - f)*airfoil_grid_u(k)%y + f*airfoil_grid_u(k - 1)%y)) then
1230 !!IB
1231 ib_markers_sf(i, j, l) = patch_id
1232 end if
1233 end if
1234 else
1235 k = 1
1236 do while (airfoil_grid_l(k)%x < xyz_local(1))
1237 k = k + 1
1238 end do
1239 if (f_approx_equal(airfoil_grid_l(k)%x, xyz_local(1))) then
1240 if (xyz_local(2) >= airfoil_grid_l(k)%y) then
1241 !!IB
1242 ib_markers_sf(i, j, l) = patch_id
1243 end if
1244 else
1245 f = (airfoil_grid_l(k)%x - xyz_local(1))/(airfoil_grid_l(k)%x - airfoil_grid_l(k - 1)%x)
1246
1247 if (xyz_local(2) >= ((1._wp - f)*airfoil_grid_l(k)%y + f*airfoil_grid_l(k - 1)%y)) then
1248 !!IB
1249 ib_markers_sf(i, j, l) = patch_id
1250 end if
1251 end if
1252 end if
1253 end if
1254 end if
1255 end do
1256 end do
1257 end do
1258
1259# 585 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1260
1261# 585 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1262#if defined(MFC_OpenACC)
1263# 585 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1264!$acc end parallel loop
1265# 585 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1266#elif defined(MFC_OpenMP)
1267# 585 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1268
1269# 585 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1270
1271# 585 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1272!$omp end target teams loop
1273# 585 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1274#endif
1275# 585 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1276
1277
1278 end subroutine s_ib_3d_airfoil
1279
1280 !> The rectangular patch is a 2D geometry that may be used,
1281 !! for example, in creating a solid boundary, or pre-/post-
1282 !! shock region, in alignment with the axes of the Cartesian
1283 !! coordinate system. The geometry of such a patch is well-
1284 !! defined when its centroid and lengths in the x- and y-
1285 !! coordinate directions are provided. Please note that the
1286 !! rectangular patch DOES NOT allow for the smoothing of its
1287 !! boundaries.
1288 !! @param patch_id is the patch identifier
1289 !! @param ib_markers_sf Array to track patch ids
1290 subroutine s_ib_rectangle(patch_id, ib_markers_sf)
1291
1292 integer, intent(in) :: patch_id
1293 integer, dimension(0:m, 0:n, 0:p), intent(inout) :: ib_markers_sf
1294
1295 integer :: i, j, k !< generic loop iterators
1296 real(wp) :: pi_inf, gamma, lit_gamma !< Equation of state parameters
1297 real(wp), dimension(1:3) :: xy_local !< x and y coordinates in local IB frame
1298 real(wp), dimension(1:2) :: length, center !< x and y coordinates in local IB frame
1299 real(wp), dimension(1:3, 1:3) :: inverse_rotation
1300
1301 pi_inf = fluid_pp(1)%pi_inf
1302 gamma = fluid_pp(1)%gamma
1303 lit_gamma = (1._wp + gamma)/gamma
1304
1305 ! Transferring the rectangle's centroid and length information
1306 center(1) = patch_ib(patch_id)%x_centroid
1307 center(2) = patch_ib(patch_id)%y_centroid
1308 length(1) = patch_ib(patch_id)%length_x
1309 length(2) = patch_ib(patch_id)%length_y
1310 inverse_rotation(:, :) = patch_ib(patch_id)%rotation_matrix_inverse(:, :)
1311
1312 ! Checking whether the rectangle covers a particular cell in the
1313 ! domain and verifying whether the current patch has the permission
1314 ! to write to that cell. If both queries check out, the primitive
1315 ! variables of the current patch are assigned to this cell.
1316
1317# 625 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1318
1319# 625 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1320#if defined(MFC_OpenACC)
1321# 625 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1322!$acc parallel loop collapse(2) gang vector default(present) private(i, j, xy_local) copy(ib_markers_sf) copyin(patch_id, center, length, inverse_rotation, x_cc, y_cc)
1323# 625 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1324#elif defined(MFC_OpenMP)
1325# 625 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1326
1327# 625 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1328
1329# 625 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1330
1331# 625 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1332!$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(tofrom:ib_markers_sf) map(to:patch_id, center, length, inverse_rotation, x_cc, y_cc)
1333# 625 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1334#endif
1335# 625 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1336
1337# 627 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1338 do j = 0, n
1339 do i = 0, m
1340 ! get the x and y coordinates in the local IB frame
1341 xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp]
1342 xy_local = matmul(inverse_rotation, xy_local)
1343
1344 if (-0.5_wp*length(1) <= xy_local(1) .and. &
1345 0.5_wp*length(1) >= xy_local(1) .and. &
1346 -0.5_wp*length(2) <= xy_local(2) .and. &
1347 0.5_wp*length(2) >= xy_local(2)) then
1348
1349 ! Updating the patch identities bookkeeping variable
1350 ib_markers_sf(i, j, 0) = patch_id
1351
1352 end if
1353 end do
1354 end do
1355
1356# 644 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1357
1358# 644 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1359#if defined(MFC_OpenACC)
1360# 644 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1361!$acc end parallel loop
1362# 644 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1363#elif defined(MFC_OpenMP)
1364# 644 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1365
1366# 644 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1367
1368# 644 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1369!$omp end target teams loop
1370# 644 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1371#endif
1372# 644 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1373
1374
1375 end subroutine s_ib_rectangle
1376
1377 !> The spherical patch is a 3D geometry that may be used,
1378 !! for example, in creating a bubble or a droplet. The patch
1379 !! geometry is well-defined when its centroid and radius are
1380 !! provided. Please note that the spherical patch DOES allow
1381 !! for the smoothing of its boundary.
1382 !! @param patch_id is the patch identifier
1383 !! @param ib_markers_sf Array to track patch ids
1384 subroutine s_ib_sphere(patch_id, ib_markers_sf)
1385
1386 integer, intent(in) :: patch_id
1387 integer, dimension(0:m, 0:n, 0:p), intent(inout) :: ib_markers_sf
1388
1389 ! Generic loop iterators
1390 integer :: i, j, k
1391 real(wp) :: radius
1392 real(wp), dimension(1:3) :: center
1393
1394 !! Variables to initialize the pressure field that corresponds to the
1395 !! bubble-collapse test case found in Tiwari et al. (2013)
1396
1397 ! Transferring spherical patch's radius, centroid, smoothing patch
1398 ! identity and smoothing coefficient information
1399 center(1) = patch_ib(patch_id)%x_centroid
1400 center(2) = patch_ib(patch_id)%y_centroid
1401 center(3) = patch_ib(patch_id)%z_centroid
1402 radius = patch_ib(patch_id)%radius
1403
1404 ! Checking whether the sphere covers a particular cell in the domain
1405 ! and verifying whether the current patch has permission to write to
1406 ! that cell. If both queries check out, the primitive variables of
1407 ! the current patch are assigned to this cell.
1408
1409# 679 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1410
1411# 679 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1412#if defined(MFC_OpenACC)
1413# 679 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1414!$acc parallel loop collapse(3) gang vector default(present) private(i, j, k, cart_y, cart_z) copy(ib_markers_sf) copyin(patch_id, center, radius)
1415# 679 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1416#elif defined(MFC_OpenMP)
1417# 679 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1418
1419# 679 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1420
1421# 679 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1422
1423# 679 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1424!$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(tofrom:ib_markers_sf) map(to:patch_id, center, radius)
1425# 679 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1426#endif
1427# 679 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1428
1429# 681 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1430 do k = 0, p
1431 do j = 0, n
1432 do i = 0, m
1433 if (grid_geometry == 3) then
1434 call s_convert_cylindrical_to_cartesian_coord(y_cc(j), z_cc(k))
1435 else
1436 cart_y = y_cc(j)
1437 cart_z = z_cc(k)
1438 end if
1439 ! Updating the patch identities bookkeeping variable
1440 if (((x_cc(i) - center(1))**2 &
1441 + (cart_y - center(2))**2 &
1442 + (cart_z - center(3))**2 <= radius**2)) then
1443 ib_markers_sf(i, j, k) = patch_id
1444 end if
1445 end do
1446 end do
1447 end do
1448
1449# 699 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1450
1451# 699 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1452#if defined(MFC_OpenACC)
1453# 699 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1454!$acc end parallel loop
1455# 699 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1456#elif defined(MFC_OpenMP)
1457# 699 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1458
1459# 699 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1460
1461# 699 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1462!$omp end target teams loop
1463# 699 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1464#endif
1465# 699 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1466
1467
1468 end subroutine s_ib_sphere
1469
1470 !> The cuboidal patch is a 3D geometry that may be used, for
1471 !! example, in creating a solid boundary, or pre-/post-shock
1472 !! region, which is aligned with the axes of the Cartesian
1473 !! coordinate system. The geometry of such a patch is well-
1474 !! defined when its centroid and lengths in the x-, y- and
1475 !! z-coordinate directions are provided. Please notice that
1476 !! the cuboidal patch DOES NOT allow for the smearing of its
1477 !! boundaries.
1478 !! @param patch_id is the patch identifier
1479 !! @param ib_markers_sf Array to track patch ids
1480 subroutine s_ib_cuboid(patch_id, ib_markers_sf)
1481
1482 integer, intent(in) :: patch_id
1483 integer, dimension(0:m, 0:n, 0:p), intent(inout) :: ib_markers_sf
1484
1485 integer :: i, j, k !< Generic loop iterators
1486 real(wp), dimension(1:3) :: xyz_local, center, length !< x and y coordinates in local IB frame
1487 real(wp), dimension(1:3, 1:3) :: inverse_rotation
1488
1489 ! Transferring the cuboid's centroid and length information
1490 center(1) = patch_ib(patch_id)%x_centroid
1491 center(2) = patch_ib(patch_id)%y_centroid
1492 center(3) = patch_ib(patch_id)%z_centroid
1493 length(1) = patch_ib(patch_id)%length_x
1494 length(2) = patch_ib(patch_id)%length_y
1495 length(3) = patch_ib(patch_id)%length_z
1496 inverse_rotation(:, :) = patch_ib(patch_id)%rotation_matrix_inverse(:, :)
1497
1498 ! Checking whether the cuboid covers a particular cell in the domain
1499 ! and verifying whether the current patch has permission to write to
1500 ! to that cell. If both queries check out, the primitive variables
1501 ! of the current patch are assigned to this cell.
1502
1503# 735 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1504
1505# 735 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1506#if defined(MFC_OpenACC)
1507# 735 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1508!$acc parallel loop collapse(3) gang vector default(present) private(i, j, k, xyz_local, cart_y, cart_z) copy(ib_markers_sf) copyin(patch_id, center, length, inverse_rotation)
1509# 735 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1510#elif defined(MFC_OpenMP)
1511# 735 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1512
1513# 735 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1514
1515# 735 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1516
1517# 735 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1518!$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(tofrom:ib_markers_sf) map(to:patch_id, center, length, inverse_rotation)
1519# 735 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1520#endif
1521# 735 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1522
1523# 737 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1524 do k = 0, p
1525 do j = 0, n
1526 do i = 0, m
1527
1528 if (grid_geometry == 3) then
1529 ! TODO :: This does not work and is not covered by any tests. This should be fixed
1530 call s_convert_cylindrical_to_cartesian_coord(y_cc(j), z_cc(k))
1531 else
1532 cart_y = y_cc(j)
1533 cart_z = z_cc(k)
1534 end if
1535 xyz_local = [x_cc(i), cart_y, cart_z] - center ! get coordinate frame centered on IB
1536 xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinates
1537
1538 if (-0.5*length(1) <= xyz_local(1) .and. &
1539 0.5*length(1) >= xyz_local(1) .and. &
1540 -0.5*length(2) <= xyz_local(2) .and. &
1541 0.5*length(2) >= xyz_local(2) .and. &
1542 -0.5*length(3) <= xyz_local(3) .and. &
1543 0.5*length(3) >= xyz_local(3)) then
1544
1545 ! Updating the patch identities bookkeeping variable
1546 ib_markers_sf(i, j, k) = patch_id
1547 end if
1548 end do
1549 end do
1550 end do
1551
1552# 764 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1553
1554# 764 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1555#if defined(MFC_OpenACC)
1556# 764 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1557!$acc end parallel loop
1558# 764 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1559#elif defined(MFC_OpenMP)
1560# 764 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1561
1562# 764 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1563
1564# 764 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1565!$omp end target teams loop
1566# 764 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1567#endif
1568# 764 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1569
1570
1571 end subroutine s_ib_cuboid
1572
1573 !> The cylindrical patch is a 3D geometry that may be used,
1574 !! for example, in setting up a cylindrical solid boundary
1575 !! confinement, like a blood vessel. The geometry of this
1576 !! patch is well-defined when the centroid, the radius and
1577 !! the length along the cylinder's axis, parallel to the x-,
1578 !! y- or z-coordinate direction, are provided. Please note
1579 !! that the cylindrical patch DOES allow for the smoothing
1580 !! of its lateral boundary.
1581 !! @param patch_id is the patch identifier
1582 !! @param ib_markers_sf Array to track patch ids
1583 subroutine s_ib_cylinder(patch_id, ib_markers_sf)
1584
1585 integer, intent(in) :: patch_id
1586 integer, dimension(0:m, 0:n, 0:p), intent(inout) :: ib_markers_sf
1587
1588 integer :: i, j, k !< Generic loop iterators
1589 real(wp) :: radius
1590 real(wp), dimension(1:3) :: xyz_local, center, length !< x and y coordinates in local IB frame
1591 real(wp), dimension(1:3, 1:3) :: inverse_rotation
1592
1593 ! Transferring the cylindrical patch's centroid, length, radius,
1594 ! smoothing patch identity and smoothing coefficient information
1595
1596 center(1) = patch_ib(patch_id)%x_centroid
1597 center(2) = patch_ib(patch_id)%y_centroid
1598 center(3) = patch_ib(patch_id)%z_centroid
1599 length(1) = patch_ib(patch_id)%length_x
1600 length(2) = patch_ib(patch_id)%length_y
1601 length(3) = patch_ib(patch_id)%length_z
1602 radius = patch_ib(patch_id)%radius
1603 inverse_rotation(:, :) = patch_ib(patch_id)%rotation_matrix_inverse(:, :)
1604
1605 ! Checking whether the cylinder covers a particular cell in the
1606 ! domain and verifying whether the current patch has the permission
1607 ! to write to that cell. If both queries check out, the primitive
1608 ! variables of the current patch are assigned to this cell.
1609
1610# 804 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1611
1612# 804 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1613#if defined(MFC_OpenACC)
1614# 804 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1615!$acc parallel loop collapse(3) gang vector default(present) private(i, j, k, xyz_local, cart_y, cart_z) copy(ib_markers_sf) copyin(patch_id, center, length, radius, inverse_rotation)
1616# 804 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1617#elif defined(MFC_OpenMP)
1618# 804 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1619
1620# 804 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1621
1622# 804 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1623
1624# 804 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1625!$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(tofrom:ib_markers_sf) map(to:patch_id, center, length, radius, inverse_rotation)
1626# 804 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1627#endif
1628# 804 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1629
1630# 806 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1631 do k = 0, p
1632 do j = 0, n
1633 do i = 0, m
1634
1635 if (grid_geometry == 3) then
1636 call s_convert_cylindrical_to_cartesian_coord(y_cc(j), z_cc(k))
1637 else
1638 cart_y = y_cc(j)
1639 cart_z = z_cc(k)
1640 end if
1641 xyz_local = [x_cc(i), cart_y, cart_z] - center ! get coordinate frame centered on IB
1642 xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinates
1643
1644 if (((.not. f_is_default(length(1)) .and. &
1645 xyz_local(2)**2 &
1646 + xyz_local(3)**2 <= radius**2 .and. &
1647 -0.5_wp*length(1) <= xyz_local(1) .and. &
1648 0.5_wp*length(1) >= xyz_local(1)) &
1649 .or. &
1650 (.not. f_is_default(length(2)) .and. &
1651 xyz_local(1)**2 &
1652 + xyz_local(3)**2 <= radius**2 .and. &
1653 -0.5_wp*length(2) <= xyz_local(2) .and. &
1654 0.5_wp*length(2) >= xyz_local(2)) &
1655 .or. &
1656 (.not. f_is_default(length(3)) .and. &
1657 xyz_local(1)**2 &
1658 + xyz_local(2)**2 <= radius**2 .and. &
1659 -0.5_wp*length(3) <= xyz_local(3) .and. &
1660 0.5_wp*length(3) >= xyz_local(3)))) then
1661
1662 ! Updating the patch identities bookkeeping variable
1663 ib_markers_sf(i, j, k) = patch_id
1664 end if
1665 end do
1666 end do
1667 end do
1668
1669# 843 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1670
1671# 843 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1672#if defined(MFC_OpenACC)
1673# 843 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1674!$acc end parallel loop
1675# 843 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1676#elif defined(MFC_OpenMP)
1677# 843 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1678
1679# 843 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1680
1681# 843 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1682!$omp end target teams loop
1683# 843 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1684#endif
1685# 843 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1686
1687
1688 end subroutine s_ib_cylinder
1689
1690 !> @brief Marks cells inside a 2D elliptical immersed boundary defined by semi-axis lengths and rotation.
1691 subroutine s_ib_ellipse(patch_id, ib_markers_sf)
1692
1693 integer, intent(in) :: patch_id
1694 integer, dimension(0:m, 0:n, 0:p), intent(inout) :: ib_markers_sf
1695
1696 integer :: i, j, k !< generic loop iterators
1697 real(wp), dimension(1:3) :: xy_local !< x and y coordinates in local IB frame
1698 real(wp), dimension(1:2) :: ellipse_coeffs !< a and b in the ellipse coefficients
1699 real(wp), dimension(1:2) :: center !< x and y coordinates in local IB frame
1700 real(wp), dimension(1:3, 1:3) :: inverse_rotation
1701
1702 ! Transferring the ellipse's centroid and length information
1703 center(1) = patch_ib(patch_id)%x_centroid
1704 center(2) = patch_ib(patch_id)%y_centroid
1705 ellipse_coeffs(1) = 0.5_wp*patch_ib(patch_id)%length_x
1706 ellipse_coeffs(2) = 0.5_wp*patch_ib(patch_id)%length_y
1707 inverse_rotation(:, :) = patch_ib(patch_id)%rotation_matrix_inverse(:, :)
1708
1709 ! Checking whether the ellipse covers a particular cell in the
1710 ! domain
1711
1712# 868 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1713
1714# 868 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1715#if defined(MFC_OpenACC)
1716# 868 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1717!$acc parallel loop collapse(2) gang vector default(present) private(i, j, xy_local) copy(ib_markers_sf) copyin(patch_id, center, ellipse_coeffs, inverse_rotation, x_cc, y_cc)
1718# 868 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1719#elif defined(MFC_OpenMP)
1720# 868 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1721
1722# 868 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1723
1724# 868 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1725
1726# 868 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1727!$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(tofrom:ib_markers_sf) map(to:patch_id, center, ellipse_coeffs, inverse_rotation, x_cc, y_cc)
1728# 868 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1729#endif
1730# 868 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1731
1732# 870 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1733 do j = 0, n
1734 do i = 0, m
1735 ! get the x and y coordinates in the local IB frame
1736 xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp]
1737 xy_local = matmul(inverse_rotation, xy_local)
1738
1739 ! Ellipse condition (x/a)^2 + (y/b)^2 <= 1
1740 if ((xy_local(1)/ellipse_coeffs(1))**2 + (xy_local(2)/ellipse_coeffs(2))**2 <= 1._wp) then
1741 ! Updating the patch identities bookkeeping variable
1742 ib_markers_sf(i, j, 0) = patch_id
1743 end if
1744 end do
1745 end do
1746
1747# 883 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1748
1749# 883 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1750#if defined(MFC_OpenACC)
1751# 883 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1752!$acc end parallel loop
1753# 883 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1754#elif defined(MFC_OpenMP)
1755# 883 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1756
1757# 883 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1758
1759# 883 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1760!$omp end target teams loop
1761# 883 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1762#endif
1763# 883 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1764
1765
1766 end subroutine s_ib_ellipse
1767
1768 !> The STL patch is a 2/3D geometry that is imported from an STL file.
1769 !! @param patch_id is the patch identifier
1770 !! @param ib_markers_sf Array to track patch ids
1771 subroutine s_ib_model(patch_id, ib_markers_sf)
1772
1773 integer, intent(in) :: patch_id
1774 integer, dimension(0:m, 0:n, 0:p), intent(inout) :: ib_markers_sf
1775
1776 integer :: i, j, k !< Generic loop iterators
1777
1778 type(t_model), pointer :: model
1779
1780 real(wp) :: eta
1781 real(wp), dimension(1:3) :: point, local_point, offset
1782 real(wp), dimension(1:3) :: center, xyz_local
1783 real(wp), dimension(1:3, 1:3) :: inverse_rotation
1784
1785 model => models(patch_id)%model
1786 center = 0._wp
1787 center(1) = patch_ib(patch_id)%x_centroid
1788 center(2) = patch_ib(patch_id)%y_centroid
1789 if (p > 0) center(3) = patch_ib(patch_id)%z_centroid
1790 inverse_rotation(:, :) = patch_ib(patch_id)%rotation_matrix_inverse(:, :)
1791 offset(:) = patch_ib(patch_id)%centroid_offset(:)
1792
1793 do i = 0, m
1794 do j = 0, n
1795 do k = 0, p
1796
1797 xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp]
1798 if (p > 0) then
1799 xyz_local(3) = z_cc(k) - center(3)
1800 end if
1801 xyz_local = matmul(inverse_rotation, xyz_local)
1802 xyz_local = xyz_local - offset
1803
1804 if (grid_geometry == 3) then
1805 xyz_local = f_convert_cyl_to_cart(xyz_local)
1806 end if
1807
1808 if (p == 0) then
1809 eta = f_model_is_inside(model, xyz_local, (/dx(i), dy(j), 0._wp/), patch_ib(patch_id)%model_spc)
1810 else
1811 eta = f_model_is_inside(model, xyz_local, (/dx(i), dy(j), dz(k)/), patch_ib(patch_id)%model_spc)
1812 end if
1813
1814 ! Reading STL boundary vertices and compute the levelset and levelset_norm
1815 if (eta > patch_ib(patch_id)%model_threshold) then
1816 ib_markers_sf(i, j, k) = patch_id
1817 end if
1818
1819 end do
1820 end do
1821 end do
1822
1823 end subroutine s_ib_model
1824
1825 !> Subroutine that computes a rotation matrix for converting to the rotating frame of the boundary
1826 subroutine s_update_ib_rotation_matrix(patch_id)
1827
1828 integer, intent(in) :: patch_id
1829 integer :: i
1830
1831 real(wp), dimension(3, 3, 3) :: rotation
1832 real(wp) :: angle
1833
1834 ! construct the x, y, and z rotation matrices
1835 if (num_dims == 3) then
1836 ! also compute the x and y axes in 3D
1837 angle = patch_ib(patch_id)%angles(1)
1838 rotation(1, 1, :) = [1._wp, 0._wp, 0._wp]
1839 rotation(1, 2, :) = [0._wp, cos(angle), -sin(angle)]
1840 rotation(1, 3, :) = [0._wp, sin(angle), cos(angle)]
1841
1842 angle = patch_ib(patch_id)%angles(2)
1843 rotation(2, 1, :) = [cos(angle), 0._wp, sin(angle)]
1844 rotation(2, 2, :) = [0._wp, 1._wp, 0._wp]
1845 rotation(2, 3, :) = [-sin(angle), 0._wp, cos(angle)]
1846
1847 ! apply the y rotation to the x rotation
1848 patch_ib(patch_id)%rotation_matrix(:, :) = matmul(rotation(1, :, :), rotation(2, :, :))
1849 patch_ib(patch_id)%rotation_matrix_inverse(:, :) = matmul(transpose(rotation(2, :, :)), transpose(rotation(1, :, :)))
1850 end if
1851
1852 ! z component first, since it applies in 2D and 3D
1853 angle = patch_ib(patch_id)%angles(3)
1854 rotation(3, 1, :) = [cos(angle), -sin(angle), 0._wp]
1855 rotation(3, 2, :) = [sin(angle), cos(angle), 0._wp]
1856 rotation(3, 3, :) = [0._wp, 0._wp, 1._wp]
1857
1858 if (num_dims == 3) then
1859 ! apply the z rotation to the xy rotation in 3D
1860 patch_ib(patch_id)%rotation_matrix(:, :) = matmul(patch_ib(patch_id)%rotation_matrix(:, :), rotation(3, :, :))
1861 patch_ib(patch_id)%rotation_matrix_inverse(:, :) = matmul(transpose(rotation(3, :, :)), patch_ib(patch_id)%rotation_matrix_inverse(:, :))
1862 else
1863 ! write out only the z rotation in 2D
1864 patch_ib(patch_id)%rotation_matrix(:, :) = rotation(3, :, :)
1865 patch_ib(patch_id)%rotation_matrix_inverse(:, :) = transpose(rotation(3, :, :))
1866 end if
1867
1868 end subroutine s_update_ib_rotation_matrix
1869
1870 !> @brief Converts cylindrical (r, theta) coordinates to Cartesian (y, z) and stores in module variables.
1871 subroutine s_convert_cylindrical_to_cartesian_coord(cyl_y, cyl_z)
1872
1873# 991 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1874#if MFC_OpenACC
1875# 991 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1876!$acc routine seq
1877# 991 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1878#elif MFC_OpenMP
1879# 991 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1880
1881# 991 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1882
1883# 991 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1884!$omp declare target device_type(any)
1885# 991 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1886#endif
1887
1888 real(wp), intent(in) :: cyl_y, cyl_z
1889
1890 cart_y = cyl_y*sin(cyl_z)
1891 cart_z = cyl_y*cos(cyl_z)
1892
1893 end subroutine s_convert_cylindrical_to_cartesian_coord
1894
1895 !> @brief Converts a 3D cylindrical coordinate vector (x, r, theta) to Cartesian (x, y, z).
1896 pure function f_convert_cyl_to_cart(cyl) result(cart)
1897
1898
1899# 1003 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1900#if MFC_OpenACC
1901# 1003 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1902!$acc routine seq
1903# 1003 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1904#elif MFC_OpenMP
1905# 1003 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1906
1907# 1003 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1908
1909# 1003 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1910!$omp declare target device_type(any)
1911# 1003 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1912#endif
1913
1914 real(wp), dimension(1:3), intent(in) :: cyl
1915 real(wp), dimension(1:3) :: cart
1916
1917 cart = (/cyl(1), &
1918 cyl(2)*sin(cyl(3)), &
1919 cyl(2)*cos(cyl(3))/)
1920
1921 end function f_convert_cyl_to_cart
1922
1923 !> @brief Converts cylindrical coordinates (x, r) to the spherical azimuthal angle phi and stores in a module variable.
1924 subroutine s_convert_cylindrical_to_spherical_coord(cyl_x, cyl_y)
1925
1926# 1016 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1927#if MFC_OpenACC
1928# 1016 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1929!$acc routine seq
1930# 1016 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1931#elif MFC_OpenMP
1932# 1016 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1933
1934# 1016 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1935
1936# 1016 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1937!$omp declare target device_type(any)
1938# 1016 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1939#endif
1940
1941 real(wp), intent(IN) :: cyl_x, cyl_y
1942
1943 sph_phi = atan(cyl_y/cyl_x)
1944
1945 end subroutine s_convert_cylindrical_to_spherical_coord
1946
1947 !> Archimedes spiral function
1948 !! @param myth Angle
1949 !! @param offset Thickness
1950 !! @param a Starting position
1951 pure elemental function f_r(myth, offset, a)
1952
1953# 1029 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1954#if MFC_OpenACC
1955# 1029 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1956!$acc routine seq
1957# 1029 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1958#elif MFC_OpenMP
1959# 1029 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1960
1961# 1029 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1962
1963# 1029 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1964!$omp declare target device_type(any)
1965# 1029 "/home/runner/work/MFC/MFC/src/simulation/m_ib_patches.fpp"
1966#endif
1967 real(wp), intent(in) :: myth, offset, a
1968 real(wp) :: b
1969 real(wp) :: f_r
1970
1971 !r(th) = a + b*th
1972
1973 b = 2._wp*a/(2._wp*pi)
1974 f_r = a + b*myth + offset
1975 end function f_r
1976
1977end module m_ib_patches
k
integer, intent(in) k
Definition m_phase_change.fpp.f90:1121
j
integer, intent(in) j
Definition m_phase_change.fpp.f90:1121
l
integer, intent(in) l
Definition m_phase_change.fpp.f90:1121
m_derived_types
Shared derived types for field data, patch geometry, bubble dynamics, and MPI I/O structures.
Definition m_derived_types.fpp.f90:292
m_global_parameters
Global parameters for the computational domain, fluid properties, and simulation algorithm configurat...
Definition m_global_parameters.fpp.f90:303
m_global_parameters::patch_ib
type(ib_patch_parameters), dimension(num_patches_max) patch_ib
Definition m_global_parameters.fpp.f90:1005
m_global_parameters::proc_rank
integer proc_rank
Rank of the local processor.
Definition m_global_parameters.fpp.f90:327
m_global_parameters::airfoil_grid_u
type(vec3_dt), dimension(:), allocatable airfoil_grid_u
Definition m_global_parameters.fpp.f90:1006
m_global_parameters::y_cc
real(wp), dimension(:), allocatable, target y_cc
Definition m_global_parameters.fpp.f90:368
m_global_parameters::z_cc
real(wp), dimension(:), allocatable, target z_cc
Definition m_global_parameters.fpp.f90:368
m_global_parameters::x_cc
real(wp), dimension(:), allocatable, target x_cc
Definition m_global_parameters.fpp.f90:368
m_global_parameters::p
integer p
Definition m_global_parameters.fpp.f90:331
m_global_parameters::n
integer n
Definition m_global_parameters.fpp.f90:331
m_global_parameters::np
integer np
Definition m_global_parameters.fpp.f90:1007
m_global_parameters::dy
real(wp), dimension(:), allocatable, target dy
Definition m_global_parameters.fpp.f90:375
m_global_parameters::m
integer m
Definition m_global_parameters.fpp.f90:331
m_global_parameters::dz
real(wp), dimension(:), allocatable, target dz
Definition m_global_parameters.fpp.f90:375
m_global_parameters::dx
real(wp), dimension(:), allocatable, target dx
Definition m_global_parameters.fpp.f90:375
m_global_parameters::num_ibs
integer num_ibs
Definition m_global_parameters.fpp.f90:1003
m_helper_basic
Basic floating-point utilities: approximate equality, default detection, and coordinate bounds.
Definition m_helper_basic.fpp.f90:292
m_helper_basic::f_approx_equal
logical elemental function, public f_approx_equal(a, b, tol_input)
This procedure checks if two floating point numbers of wp are within tolerance.
Definition m_helper_basic.fpp.f90:315
m_helper
Utility routines for bubble model setup, coordinate transforms, array sampling, and special functions...
Definition m_helper.fpp.f90:303
m_helper::s_transform_model
subroutine, public s_transform_model(model, matrix, matrix_n)
This procedure transforms a model by a matrix, one triangle at a time.
Definition m_helper.fpp.f90:926
m_helper::f_create_bbox
type(t_bbox) function, public f_create_bbox(model)
This procedure creates a bounding box for a model.
Definition m_helper.fpp.f90:942
m_helper::f_create_transform_matrix
real(wp) function, dimension(1:4, 1:4), public f_create_transform_matrix(param, center)
This procedure creates a transformation matrix.
Definition m_helper.fpp.f90:830
m_ib_patches
Allocate memory and read initial condition data for IC extrusion.
Definition m_ib_patches.fpp.f90:362
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:1827
m_ib_patches::smooth_coeff
real(wp) smooth_coeff
Definition m_ib_patches.fpp.f90:406
m_ib_patches::cart_x
real(wp) cart_x
Definition m_ib_patches.fpp.f90:423
m_ib_patches::models
type(t_model_array), dimension(:), allocatable, target, public models
Definition m_ib_patches.fpp.f90:458
m_ib_patches::s_ib_rectangle
subroutine s_ib_rectangle(patch_id, ib_markers_sf)
The rectangular patch is a 2D geometry that may be used, for example, in creating a solid boundary,...
Definition m_ib_patches.fpp.f90:1291
m_ib_patches::z_boundary
type(bounds_info) z_boundary
Definition m_ib_patches.fpp.f90:439
m_ib_patches::x_centroid
real(wp) x_centroid
Definition m_ib_patches.fpp.f90:380
m_ib_patches::length_z
real(wp) length_z
Definition m_ib_patches.fpp.f90:381
m_ib_patches::smooth_patch_id
integer smooth_patch_id
Definition m_ib_patches.fpp.f90:405
m_ib_patches::cart_z
real(wp) cart_z
Definition m_ib_patches.fpp.f90:423
m_ib_patches::s_apply_ib_patches
impure subroutine, public s_apply_ib_patches(ib_markers_sf)
Applies all immersed boundary patch geometries to mark interior cells in the IB marker array.
Definition m_ib_patches.fpp.f90:476
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:1872
m_ib_patches::s_ib_ellipse
subroutine s_ib_ellipse(patch_id, ib_markers_sf)
Marks cells inside a 2D elliptical immersed boundary defined by semi-axis lengths and rotation.
Definition m_ib_patches.fpp.f90:1692
m_ib_patches::y_boundary
type(bounds_info) y_boundary
Definition m_ib_patches.fpp.f90:439
m_ib_patches::y_centroid
real(wp) y_centroid
Definition m_ib_patches.fpp.f90:380
m_ib_patches::sph_phi
real(wp) sph_phi
Definition m_ib_patches.fpp.f90:424
m_ib_patches::istr
character(len=5) istr
Definition m_ib_patches.fpp.f90:456
m_ib_patches::length_y
real(wp) length_y
Definition m_ib_patches.fpp.f90:381
m_ib_patches::z_centroid
real(wp) z_centroid
Definition m_ib_patches.fpp.f90:380
m_ib_patches::x_boundary
type(bounds_info) x_boundary
Definition m_ib_patches.fpp.f90:439
m_ib_patches::s_ib_circle
subroutine s_ib_circle(patch_id, ib_markers_sf)
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:664
m_ib_patches::s_ib_cuboid
subroutine s_ib_cuboid(patch_id, ib_markers_sf)
The cuboidal patch is a 3D geometry that may be used, for example, in creating a solid boundary,...
Definition m_ib_patches.fpp.f90:1481
m_ib_patches::s_ib_sphere
subroutine s_ib_sphere(patch_id, ib_markers_sf)
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:1385
m_ib_patches::s_ib_3d_airfoil
subroutine s_ib_3d_airfoil(patch_id, ib_markers_sf)
Marks cells inside a 3D extruded NACA 4-digit airfoil immersed boundary with finite span.
Definition m_ib_patches.fpp.f90:1014
m_ib_patches::length_x
real(wp) length_x
Definition m_ib_patches.fpp.f90:381
m_ib_patches::s_ib_airfoil
subroutine s_ib_airfoil(patch_id, ib_markers_sf)
Marks cells inside a 2D NACA 4-digit airfoil immersed boundary using upper and lower surface grids.
Definition m_ib_patches.fpp.f90:742
m_ib_patches::s_instantiate_stl_models
subroutine, public s_instantiate_stl_models()
Definition m_ib_patches.fpp.f90:530
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:1897
m_ib_patches::s_ib_model
subroutine s_ib_model(patch_id, ib_markers_sf)
The STL patch is a 2/3D geometry that is imported from an STL file.
Definition m_ib_patches.fpp.f90:1772
m_ib_patches::cart_y
real(wp) cart_y
Definition m_ib_patches.fpp.f90:423
m_ib_patches::s_ib_cylinder
subroutine s_ib_cylinder(patch_id, ib_markers_sf)
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:1584
m_ib_patches::f_r
pure elemental real(wp) function f_r(myth, offset, a)
Archimedes spiral function.
Definition m_ib_patches.fpp.f90:1952
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:1925
m_model
Binary STL file reader and processor for immersed boundary geometry.
Definition m_model.fpp.f90:293
m_model::f_check_interpolation_2d
subroutine, public f_check_interpolation_2d(boundary_v, boundary_edge_count, spacing, interpolate)
This procedure check if interpolates is needed for 2D models.
Definition m_model.fpp.f90:1018
m_model::f_check_interpolation_3d
subroutine, public f_check_interpolation_3d(model, spacing, interpolate)
This procedure check if interpolates is needed for 3D models.
Definition m_model.fpp.f90:1048
m_model::f_check_boundary
subroutine, public f_check_boundary(model, boundary_v, boundary_vertex_count, boundary_edge_count)
This procedure checks and labels edges shared by two or more triangles facets of the 2D STL model.
Definition m_model.fpp.f90:891
m_model::f_interpolate_2d
subroutine, public f_interpolate_2d(boundary_v, boundary_edge_count, spacing, interpolated_boundary_v, total_vertices)
This procedure interpolates 2D models.
Definition m_model.fpp.f90:1095
m_model::f_model_read
impure type(t_model) function, public f_model_read(filepath)
This procedure reads a mesh from a file.
Definition m_model.fpp.f90:585
m_model::f_interpolate_3d
impure subroutine, public f_interpolate_3d(model, spacing, interpolated_boundary_v, total_vertices)
This procedure interpolates 3D models.
Definition m_model.fpp.f90:1189
m_mpi_common
MPI communication layer: domain decomposition, halo exchange, reductions, and parallel I/O setup.
Definition m_mpi_common.fpp.f90:303
m_derived_types::bounds_info
Derived type adding beginning (beg) and end bounds info as attributes.
Definition m_derived_types.fpp.f90:409
m_derived_types::ic_model_parameters
Defines parameters for a Model Patch.
Definition m_derived_types.fpp.f90:427
m_derived_types::t_bbox
Definition m_derived_types.fpp.f90:458
m_derived_types::t_model_array
Definition m_derived_types.fpp.f90:469
m_derived_types::t_model
Definition m_derived_types.fpp.f90:463