MFC
Exascale flow solver
Loading...
Searching...
No Matches
m_ibm.fpp.f90
Go to the documentation of this file.
1# 1 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2!>
3!! @file
4!! @brief Contains module m_ibm
5
6# 1 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 1
7# 1 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 1
8# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
9# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
10# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
11# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
12# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
13# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
14
15# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
16# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
17# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
18
19# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
20
21# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
22
23# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
24
25# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
26
27# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
28
29# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
30
31# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
32! New line at end of file is required for FYPP
33# 2 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
34# 1 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 1
35# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
36# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
37# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
38# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
39# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
40# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
41
42# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
43# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
44# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
45
46# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
47
48# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
49
50# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
51
52# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
53
54# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
55
56# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
57
58# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
59! New line at end of file is required for FYPP
60# 2 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 2
61
62# 4 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
63# 5 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
64# 6 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
65# 7 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
66# 8 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
67
68# 20 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
69
70# 43 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
71
72# 48 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
73
74# 53 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
75
76# 58 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
77
78# 63 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
79
80# 68 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
81
82# 76 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
83
84# 81 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
85
86# 86 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
87
88# 91 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
89
90# 96 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
91
92# 101 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
93
94# 106 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
95
96# 111 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
97
98# 116 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
99
100# 121 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
101
102# 151 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
103
104# 192 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
105
106# 207 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
107
108# 232 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
109
110# 243 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
111
112# 245 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
113# 255 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
114
115# 283 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
116
117# 293 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
118
119# 303 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
120
121# 312 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
122
123# 329 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
124
125# 339 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
126
127# 346 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
128
129# 352 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
130
131# 358 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
132
133# 364 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
134
135# 370 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
136
137# 376 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
138! New line at end of file is required for FYPP
139# 3 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
140# 1 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 1
141# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
142# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
143# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
144# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
145# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
146# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
147
148# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
149# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
150# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
151
152# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
153
154# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
155
156# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
157
158# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
159
160# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
161
162# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
163
164# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
165! New line at end of file is required for FYPP
166# 2 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 2
167
168# 7 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
169
170# 17 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
171
172# 22 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
173
174# 27 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
175
176# 32 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
177
178# 37 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
179
180# 42 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
181
182# 47 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
183
184# 52 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
185
186# 57 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
187
188# 62 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
189
190# 73 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
191
192# 78 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
193
194# 83 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
195
196# 88 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
197
198# 103 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
199
200# 131 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
201
202# 160 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
203
204# 175 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
205
206# 192 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
207
208# 213 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
209
210# 241 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
211
212# 256 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
213
214# 266 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
215
216# 275 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
217
218# 291 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
219
220# 301 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
221
222# 308 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
223! New line at end of file is required for FYPP
224# 4 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
225
226# 21 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
227
228# 37 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
229
230# 50 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
231
232# 76 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
233
234# 91 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
235
236# 102 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
237
238# 115 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
239
240# 143 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
241
242# 154 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
243
244# 165 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
245
246# 176 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
247
248# 187 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
249
250# 198 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
251
252# 208 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
253
254# 214 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
255
256# 220 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
257
258# 226 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
259
260# 232 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
261
262# 234 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
263# 235 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
264! New line at end of file is required for FYPP
265# 2 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 2
266
267# 14 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
268
269! Caution:
270! This macro requires the use of a binding script to set CUDA_VISIBLE_DEVICES, such that we have one GPU device per MPI rank.
271! That's because for both cudaMemAdvise (preferred location) and cudaMemPrefetchAsync we use location = device_id = 0.
272! For an example see misc/nvidia_uvm/bind.sh.
273# 63 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
274
275# 81 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
276
277# 88 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
278
279# 111 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
280
281# 127 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
282
283# 153 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
284
285# 159 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
286
287# 167 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
288! New line at end of file is required for FYPP
289# 6 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp" 2
290
291!> @brief Ghost-node immersed boundary method: locates ghost/image points, computes interpolation coefficients, and corrects the flow state
292module m_ibm
293
294 use m_derived_types !< definitions of the derived types
295
296 use m_global_parameters !< definitions of the global parameters
297
298 use m_mpi_proxy !< message passing interface (mpi) module proxy
299
300 use m_variables_conversion !< state variables type conversion procedures
301
302 use m_helper
303
304 use m_helper_basic !< functions to compare floating point numbers
305
306 use m_constants
307
308 use m_compute_levelset
309
310 use m_ib_patches
311
312 use m_viscous
313
314 implicit none
315
316 private :: s_compute_image_points, &
317 s_compute_interpolation_coeffs, &
318 s_interpolate_image_point, &
319 s_find_ghost_points, &
320 s_find_num_ghost_points
321 ; public :: s_initialize_ibm_module, &
322 s_ibm_setup, &
323 s_ibm_correct_state, &
324 s_finalize_ibm_module
325
326 type(integer_field), public :: ib_markers
327
328# 43 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
329#if defined(MFC_OpenACC)
330# 43 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
331!$acc declare create(ib_markers)
332# 43 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
333#elif defined(MFC_OpenMP)
334# 43 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
335!$omp declare target (ib_markers)
336# 43 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
337#endif
338
339 type(ghost_point), dimension(:), allocatable :: ghost_points
340 type(ghost_point), dimension(:), allocatable :: inner_points
341
342# 47 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
343#if defined(MFC_OpenACC)
344# 47 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
345!$acc declare create(ghost_points, inner_points)
346# 47 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
347#elif defined(MFC_OpenMP)
348# 47 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
349!$omp declare target (ghost_points, inner_points)
350# 47 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
351#endif
352
353 integer :: num_gps !< Number of ghost points
354 integer :: num_inner_gps !< Number of ghost points
355#if defined(MFC_OpenACC)
356
357# 52 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
358#if defined(MFC_OpenACC)
359# 52 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
360!$acc declare create(gp_layers, num_gps, num_inner_gps)
361# 52 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
362#elif defined(MFC_OpenMP)
363# 52 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
364!$omp declare target (gp_layers, num_gps, num_inner_gps)
365# 52 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
366#endif
367#elif defined(MFC_OpenMP)
368
369# 54 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
370#if defined(MFC_OpenACC)
371# 54 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
372!$acc declare create(num_gps, num_inner_gps)
373# 54 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
374#elif defined(MFC_OpenMP)
375# 54 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
376!$omp declare target (num_gps, num_inner_gps)
377# 54 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
378#endif
379#endif
380 logical :: moving_immersed_boundary_flag
381
382contains
383
384 !> Allocates memory for the variables in the IBM module
385 impure subroutine s_initialize_ibm_module()
386
387 if (p > 0) then
388#ifdef MFC_DEBUG
389# 64 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
390 block
391# 64 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
392 use iso_fortran_env, only: output_unit
393# 64 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
394
395# 64 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
396 print *, 'm_ibm.fpp:64: ', '@:ALLOCATE(ib_markers%sf(-buff_size:m+buff_size, -buff_size:n+buff_size, -buff_size:p+buff_size))'
397# 64 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
398
399# 64 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
400 call flush (output_unit)
401# 64 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
402 end block
403# 64 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
404#endif
405# 64 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
406 allocate (ib_markers%sf(-buff_size:m+buff_size, -buff_size:n+buff_size, -buff_size:p+buff_size))
407# 64 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
408
409# 64 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
410
411# 64 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
412#if defined(MFC_OpenACC)
413# 64 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
414!$acc enter data create(ib_markers%sf)
415# 64 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
416#elif defined(MFC_OpenMP)
417# 64 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
418!$omp target enter data map(always,alloc:ib_markers%sf)
419# 64 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
420#endif
421# 66 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
422 else
423#ifdef MFC_DEBUG
424# 67 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
425 block
426# 67 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
427 use iso_fortran_env, only: output_unit
428# 67 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
429
430# 67 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
431 print *, 'm_ibm.fpp:67: ', '@:ALLOCATE(ib_markers%sf(-buff_size:m+buff_size, -buff_size:n+buff_size, 0:0))'
432# 67 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
433
434# 67 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
435 call flush (output_unit)
436# 67 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
437 end block
438# 67 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
439#endif
440# 67 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
441 allocate (ib_markers%sf(-buff_size:m+buff_size, -buff_size:n+buff_size, 0:0))
442# 67 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
443
444# 67 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
445
446# 67 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
447#if defined(MFC_OpenACC)
448# 67 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
449!$acc enter data create(ib_markers%sf)
450# 67 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
451#elif defined(MFC_OpenMP)
452# 67 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
453!$omp target enter data map(always,alloc:ib_markers%sf)
454# 67 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
455#endif
456# 69 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
457 end if
458
459#ifdef MFC_DEBUG
460# 71 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
461 block
462# 71 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
463 use iso_fortran_env, only: output_unit
464# 71 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
465
466# 71 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
467 print *, 'm_ibm.fpp:71: ', '@:ALLOCATE(models(num_ibs))'
468# 71 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
469
470# 71 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
471 call flush (output_unit)
472# 71 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
473 end block
474# 71 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
475#endif
476# 71 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
477 allocate (models(num_ibs))
478# 71 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
479
480# 71 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
481
482# 71 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
483#if defined(MFC_OpenACC)
484# 71 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
485!$acc enter data create(models)
486# 71 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
487#elif defined(MFC_OpenMP)
488# 71 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
489!$omp target enter data map(always,alloc:models)
490# 71 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
491#endif
492
493#ifdef _CRAYFTN
494# 73 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
495 block
496# 73 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
497
498# 73 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
499#ifdef MFC_DEBUG
500# 73 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
501 block
502# 73 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
503 use iso_fortran_env, only: output_unit
504# 73 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
505
506# 73 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
507 print *, 'm_ibm.fpp:73: ', '@:ACC_SETUP_SFs(ib_markers)'
508# 73 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
509
510# 73 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
511 call flush (output_unit)
512# 73 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
513 end block
514# 73 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
515#endif
516# 73 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
517
518# 73 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
519
520# 73 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
521#if defined(MFC_OpenACC)
522# 73 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
523!$acc enter data copyin(ib_markers)
524# 73 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
525#elif defined(MFC_OpenMP)
526# 73 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
527!$omp target enter data map(to:ib_markers)
528# 73 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
529#endif
530# 73 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
531 if (associated(ib_markers%sf)) then
532# 73 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
533
534# 73 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
535#if defined(MFC_OpenACC)
536# 73 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
537!$acc enter data copyin(ib_markers%sf)
538# 73 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
539#elif defined(MFC_OpenMP)
540# 73 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
541!$omp target enter data map(to:ib_markers%sf)
542# 73 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
543#endif
544# 73 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
545 end if
546# 73 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
547 end block
548# 73 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
549#endif
550
551
552# 75 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
553#if defined(MFC_OpenACC)
554# 75 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
555!$acc enter data copyin(num_gps, num_inner_gps)
556# 75 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
557#elif defined(MFC_OpenMP)
558# 75 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
559!$omp target enter data map(to:num_gps, num_inner_gps)
560# 75 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
561#endif
562
563 end subroutine s_initialize_ibm_module
564
565 !> Initializes the values of various IBM variables, such as ghost points and
566 !! image points.
567 impure subroutine s_ibm_setup()
568
569 integer :: i, j, k
570 integer :: max_num_gps, max_num_inner_gps
571
572 ! do all set up for moving immersed boundaries
573 moving_immersed_boundary_flag = .false.
574 do i = 1, num_ibs
575 if (patch_ib(i)%moving_ibm /= 0) then
576 call s_compute_moment_of_inertia(i, patch_ib(i)%angular_vel)
577 moving_immersed_boundary_flag = .true.
578 end if
579 call s_update_ib_rotation_matrix(i)
580 end do
581
582# 95 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
583#if defined(MFC_OpenACC)
584# 95 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
585!$acc update device(patch_ib(1:num_ibs))
586# 95 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
587#elif defined(MFC_OpenMP)
588# 95 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
589!$omp target update to(patch_ib(1:num_ibs))
590# 95 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
591#endif
592
593 ! GPU routines require updated cell centers
594
595# 98 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
596#if defined(MFC_OpenACC)
597# 98 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
598!$acc update device(x_cc, y_cc)
599# 98 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
600#elif defined(MFC_OpenMP)
601# 98 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
602!$omp target update to(x_cc, y_cc)
603# 98 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
604#endif
605 if (p /= 0) then
606
607# 100 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
608#if defined(MFC_OpenACC)
609# 100 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
610!$acc update device(z_cc)
611# 100 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
612#elif defined(MFC_OpenMP)
613# 100 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
614!$omp target update to(z_cc)
615# 100 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
616#endif
617 end if
618
619 ! allocate STL models
620 call s_instantiate_stl_models()
621
622 ! recompute the new ib_patch locations and broadcast them.
623 ib_markers%sf = 0._wp
624 call s_apply_ib_patches(ib_markers%sf(0:m, 0:n, 0:p))
625
626# 109 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
627#if defined(MFC_OpenACC)
628# 109 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
629!$acc update device(ib_markers%sf)
630# 109 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
631#elif defined(MFC_OpenMP)
632# 109 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
633!$omp target update to(ib_markers%sf)
634# 109 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
635#endif
636 call s_populate_ib_buffers()
637
638# 111 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
639#if defined(MFC_OpenACC)
640# 111 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
641!$acc update host(ib_markers%sf)
642# 111 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
643#elif defined(MFC_OpenMP)
644# 111 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
645!$omp target update from(ib_markers%sf)
646# 111 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
647#endif
648 do i = 1, num_ibs
649 if (patch_ib(i)%moving_ibm /= 0) call s_compute_centroid_offset(i) ! offsets are computed after IB markers are generated
650
651# 114 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
652#if defined(MFC_OpenACC)
653# 114 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
654!$acc update device(patch_ib(i))
655# 114 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
656#elif defined(MFC_OpenMP)
657# 114 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
658!$omp target update to(patch_ib(i))
659# 114 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
660#endif
661 end do
662
663 ! find the number of ghost points and set them to be the maximum total across ranks
664 call s_find_num_ghost_points(num_gps, num_inner_gps)
665 call s_mpi_allreduce_integer_sum(num_gps, max_num_gps)
666 call s_mpi_allreduce_integer_sum(num_inner_gps, max_num_inner_gps)
667 max_num_gps = min(max_num_gps, (m + 1)*(n + 1)*(p + 1)/2)
668 max_num_inner_gps = min(max_num_inner_gps, (m + 1)*(n + 1)*(p + 1)/2)
669
670 ! set the size of the ghost point arrays to be the amount of points total, plus a factor of 2 buffer
671
672# 125 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
673#if defined(MFC_OpenACC)
674# 125 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
675!$acc update device(num_gps, num_inner_gps)
676# 125 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
677#elif defined(MFC_OpenMP)
678# 125 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
679!$omp target update to(num_gps, num_inner_gps)
680# 125 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
681#endif
682#ifdef MFC_DEBUG
683# 126 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
684 block
685# 126 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
686 use iso_fortran_env, only: output_unit
687# 126 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
688
689# 126 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
690 print *, 'm_ibm.fpp:126: ', '@:ALLOCATE(ghost_points(1:int((max_num_gps + max_num_inner_gps) * 2.0)))'
691# 126 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
692
693# 126 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
694 call flush (output_unit)
695# 126 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
696 end block
697# 126 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
698#endif
699# 126 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
700 allocate (ghost_points(1:int((max_num_gps + max_num_inner_gps) * 2.0)))
701# 126 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
702
703# 126 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
704
705# 126 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
706#if defined(MFC_OpenACC)
707# 126 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
708!$acc enter data create(ghost_points)
709# 126 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
710#elif defined(MFC_OpenMP)
711# 126 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
712!$omp target enter data map(always,alloc:ghost_points)
713# 126 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
714#endif
715#ifdef MFC_DEBUG
716# 127 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
717 block
718# 127 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
719 use iso_fortran_env, only: output_unit
720# 127 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
721
722# 127 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
723 print *, 'm_ibm.fpp:127: ', '@:ALLOCATE(inner_points(1:int((max_num_gps + max_num_inner_gps) * 2.0)))'
724# 127 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
725
726# 127 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
727 call flush (output_unit)
728# 127 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
729 end block
730# 127 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
731#endif
732# 127 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
733 allocate (inner_points(1:int((max_num_gps + max_num_inner_gps) * 2.0)))
734# 127 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
735
736# 127 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
737
738# 127 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
739#if defined(MFC_OpenACC)
740# 127 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
741!$acc enter data create(inner_points)
742# 127 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
743#elif defined(MFC_OpenMP)
744# 127 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
745!$omp target enter data map(always,alloc:inner_points)
746# 127 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
747#endif
748
749
750# 129 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
751#if defined(MFC_OpenACC)
752# 129 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
753!$acc enter data copyin(ghost_points, inner_points)
754# 129 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
755#elif defined(MFC_OpenMP)
756# 129 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
757!$omp target enter data map(to:ghost_points, inner_points)
758# 129 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
759#endif
760
761 call s_find_ghost_points(ghost_points, inner_points)
762
763# 132 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
764#if defined(MFC_OpenACC)
765# 132 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
766!$acc update device(ghost_points, inner_points)
767# 132 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
768#elif defined(MFC_OpenMP)
769# 132 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
770!$omp target update to(ghost_points, inner_points)
771# 132 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
772#endif
773
774 call s_apply_levelset(ghost_points, num_gps)
775
776# 135 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
777#if defined(MFC_OpenACC)
778# 135 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
779!$acc update device(ghost_points)
780# 135 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
781#elif defined(MFC_OpenMP)
782# 135 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
783!$omp target update to(ghost_points)
784# 135 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
785#endif
786
787 call s_compute_image_points(ghost_points)
788
789# 138 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
790#if defined(MFC_OpenACC)
791# 138 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
792!$acc update device(ghost_points)
793# 138 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
794#elif defined(MFC_OpenMP)
795# 138 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
796!$omp target update to(ghost_points)
797# 138 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
798#endif
799
800 call s_compute_interpolation_coeffs(ghost_points)
801
802# 141 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
803#if defined(MFC_OpenACC)
804# 141 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
805!$acc update device(ghost_points)
806# 141 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
807#elif defined(MFC_OpenMP)
808# 141 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
809!$omp target update to(ghost_points)
810# 141 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
811#endif
812
813 end subroutine s_ibm_setup
814
815 !> @brief Exchanges immersed boundary marker data across MPI subdomain boundaries.
816 subroutine s_populate_ib_buffers()
817
818# 149 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
819# 150 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
820 if (bc_x%beg >= 0) then
821 call s_mpi_sendrecv_ib_buffers(ib_markers, 1, -1)
822 end if
823# 150 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
824 if (bc_x%end >= 0) then
825 call s_mpi_sendrecv_ib_buffers(ib_markers, 1, 1)
826 end if
827# 154 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
828# 149 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
829# 150 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
830 if (bc_y%beg >= 0) then
831 call s_mpi_sendrecv_ib_buffers(ib_markers, 2, -1)
832 end if
833# 150 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
834 if (bc_y%end >= 0) then
835 call s_mpi_sendrecv_ib_buffers(ib_markers, 2, 1)
836 end if
837# 154 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
838# 149 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
839# 150 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
840 if (bc_z%beg >= 0) then
841 call s_mpi_sendrecv_ib_buffers(ib_markers, 3, -1)
842 end if
843# 150 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
844 if (bc_z%end >= 0) then
845 call s_mpi_sendrecv_ib_buffers(ib_markers, 3, 1)
846 end if
847# 154 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
848# 155 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
849
850 end subroutine s_populate_ib_buffers
851
852 !> Subroutine that updates the conservative variables at the ghost points
853 !! @param pb_in Internal bubble pressure
854 !! @param mv_in Mass of vapor in bubble
855 subroutine s_ibm_correct_state(q_cons_vf, q_prim_vf, pb_in, mv_in)
856
857 type(scalar_field), &
858 dimension(sys_size), &
859 intent(INOUT) :: q_cons_vf !< Primitive Variables
860
861 type(scalar_field), &
862 dimension(sys_size), &
863 intent(INOUT) :: q_prim_vf !< Primitive Variables
864
865 real(stp), dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), optional, intent(INOUT) :: pb_in, mv_in
866
867 integer :: i, j, k, l, q, r!< Iterator variables
868 integer :: patch_id !< Patch ID of ghost point
869 real(wp) :: rho, gamma, pi_inf, dyn_pres !< Mixture variables
870 real(wp), dimension(2) :: re_k
871 real(wp) :: g_k
872 real(wp) :: qv_k
873
874 real(wp) :: pres_ip
875 real(wp), dimension(3) :: vel_ip, vel_norm_ip
876 real(wp) :: c_ip
877# 190 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
878 real(wp), dimension(num_fluids) :: gs
879 real(wp), dimension(num_fluids) :: alpha_rho_ip, alpha_ip
880 real(wp), dimension(nb) :: r_ip, v_ip, pb_ip, mv_ip
881 real(wp), dimension(nb*nmom) :: nmom_ip
882 real(wp), dimension(nb*nnode) :: presb_ip, massv_ip
883# 196 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
884 !! Primitive variables at the image point associated with a ghost point,
885 !! interpolated from surrounding fluid cells.
886
887 real(wp), dimension(3) :: norm !< Normal vector from GP to IP
888 real(wp), dimension(3) :: physical_loc !< Physical loc of GP
889 real(wp), dimension(3) :: vel_g !< Velocity of GP
890 real(wp), dimension(3) :: radial_vector !< vector from centroid to ghost point
891 real(wp), dimension(3) :: rotation_velocity !< speed of the ghost point due to rotation
892
893 real(wp) :: nbub
894 real(wp) :: buf
895 type(ghost_point) :: gp
896 type(ghost_point) :: innerp
897
898 ! set the Moving IBM interior Pressure Values
899
900# 211 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
901
902# 211 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
903#if defined(MFC_OpenACC)
904# 211 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
905!$acc parallel loop collapse(3) gang vector default(present) private(i, j, k, patch_id, rho) copyin(E_idx, momxb)
906# 211 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
907#elif defined(MFC_OpenMP)
908# 211 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
909
910# 211 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
911
912# 211 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
913
914# 211 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
915!$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, patch_id, rho) map(to:E_idx, momxb)
916# 211 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
917#endif
918# 211 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
919
920 do l = 0, p
921 do k = 0, n
922 do j = 0, m
923 patch_id = ib_markers%sf(j, k, l)
924 if (patch_id /= 0) then
925 q_prim_vf(e_idx)%sf(j, k, l) = 1._wp
926 if (patch_ib(patch_id)%moving_ibm > 0) then
927 rho = 0._wp
928 do i = 1, num_fluids
929 rho = rho + q_prim_vf(contxb + i - 1)%sf(j, k, l)
930 end do
931
932 ! Sets the momentum
933 do i = 1, num_dims
934 q_cons_vf(momxb + i - 1)%sf(j, k, l) = patch_ib(patch_id)%vel(i)*rho
935 q_prim_vf(momxb + i - 1)%sf(j, k, l) = patch_ib(patch_id)%vel(i)
936 end do
937 end if
938 end if
939 end do
940 end do
941 end do
942
943# 234 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
944
945# 234 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
946#if defined(MFC_OpenACC)
947# 234 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
948!$acc end parallel loop
949# 234 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
950#elif defined(MFC_OpenMP)
951# 234 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
952
953# 234 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
954
955# 234 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
956!$omp end target teams loop
957# 234 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
958#endif
959# 234 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
960
961
962 if (num_gps > 0) then
963
964# 237 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
965
966# 237 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
967#if defined(MFC_OpenACC)
968# 237 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
969!$acc parallel loop gang vector default(present) private(i, physical_loc, dyn_pres, alpha_rho_IP, alpha_IP, pres_IP, vel_IP, vel_g, vel_norm_IP, r_IP, v_IP, pb_IP, mv_IP, nmom_IP, presb_IP, massv_IP, rho, gamma, pi_inf, Re_K, G_K, Gs, gp, innerp, norm, buf, radial_vector, rotation_velocity, j, k, l, q, qv_K, c_IP, nbub, patch_id)
970# 237 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
971#elif defined(MFC_OpenMP)
972# 237 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
973
974# 237 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
975
976# 237 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
977
978# 237 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
979!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, physical_loc, dyn_pres, alpha_rho_IP, alpha_IP, pres_IP, vel_IP, vel_g, vel_norm_IP, r_IP, v_IP, pb_IP, mv_IP, nmom_IP, presb_IP, massv_IP, rho, gamma, pi_inf, Re_K, G_K, Gs, gp, innerp, norm, buf, radial_vector, rotation_velocity, j, k, l, q, qv_K, c_IP, nbub, patch_id)
980# 237 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
981#endif
982# 237 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
983
984 do i = 1, num_gps
985
986 gp = ghost_points(i)
987 j = gp%loc(1)
988 k = gp%loc(2)
989 l = gp%loc(3)
990 patch_id = ghost_points(i)%ib_patch_id
991
992 ! Calculate physical location of GP
993 if (p > 0) then
994 physical_loc = [x_cc(j), y_cc(k), z_cc(l)]
995 else
996 physical_loc = [x_cc(j), y_cc(k), 0._wp]
997 end if
998
999 !Interpolate primitive variables at image point associated w/ GP
1000 if (bubbles_euler .and. .not. qbmm) then
1001 call s_interpolate_image_point(q_prim_vf, gp, &
1002 alpha_rho_ip, alpha_ip, pres_ip, vel_ip, c_ip, &
1003 r_ip, v_ip, pb_ip, mv_ip)
1004 else if (qbmm .and. polytropic) then
1005 call s_interpolate_image_point(q_prim_vf, gp, &
1006 alpha_rho_ip, alpha_ip, pres_ip, vel_ip, c_ip, &
1007 r_ip, v_ip, pb_ip, mv_ip, nmom_ip)
1008 else if (qbmm .and. .not. polytropic) then
1009 call s_interpolate_image_point(q_prim_vf, gp, &
1010 alpha_rho_ip, alpha_ip, pres_ip, vel_ip, c_ip, &
1011 r_ip, v_ip, pb_ip, mv_ip, nmom_ip, pb_in, mv_in, presb_ip, massv_ip)
1012 else
1013 call s_interpolate_image_point(q_prim_vf, gp, &
1014 alpha_rho_ip, alpha_ip, pres_ip, vel_ip, c_ip)
1015 end if
1016
1017 dyn_pres = 0._wp
1018
1019 ! Set q_prim_vf params at GP so that mixture vars calculated properly
1020
1021# 274 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1022#if defined(MFC_OpenACC)
1023# 274 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1024!$acc loop seq
1025# 274 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1026#elif defined(MFC_OpenMP)
1027# 274 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1028
1029# 274 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1030#endif
1031 do q = 1, num_fluids
1032 q_prim_vf(q)%sf(j, k, l) = alpha_rho_ip(q)
1033 q_prim_vf(advxb + q - 1)%sf(j, k, l) = alpha_ip(q)
1034 end do
1035
1036 if (surface_tension) then
1037 q_prim_vf(c_idx)%sf(j, k, l) = c_ip
1038 end if
1039
1040 ! set the pressure
1041 if (patch_ib(patch_id)%moving_ibm <= 1) then
1042 q_prim_vf(e_idx)%sf(j, k, l) = pres_ip
1043 else
1044 q_prim_vf(e_idx)%sf(j, k, l) = 0._wp
1045
1046# 289 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1047#if defined(MFC_OpenACC)
1048# 289 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1049!$acc loop seq
1050# 289 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1051#elif defined(MFC_OpenMP)
1052# 289 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1053
1054# 289 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1055#endif
1056 do q = 1, num_fluids
1057 ! Se the pressure inside a moving immersed boundary based upon the pressure of the image point. acceleration, and normal vector direction
1058 q_prim_vf(e_idx)%sf(j, k, l) = q_prim_vf(e_idx)%sf(j, k, l) + pres_ip/(1._wp - 2._wp*abs(gp%levelset*alpha_rho_ip(q)/pres_ip)*dot_product(patch_ib(patch_id)%force/patch_ib(patch_id)%mass, gp%levelset_norm))
1059 end do
1060 end if
1061
1062 if (model_eqns /= 4) then
1063 ! If in simulation, use acc mixture subroutines
1064 if (elasticity) then
1065 call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv_k, alpha_ip, &
1066 alpha_rho_ip, re_k, g_k, gs)
1067 else
1068 call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv_k, alpha_ip, &
1069 alpha_rho_ip, re_k)
1070 end if
1071 end if
1072
1073 ! Calculate velocity of ghost cell
1074 if (gp%slip) then
1075 norm(1:3) = gp%levelset_norm
1076 buf = sqrt(sum(norm**2))
1077 norm = norm/buf
1078 vel_norm_ip = sum(vel_ip*norm)*norm
1079 vel_g = vel_ip - vel_norm_ip
1080 if (patch_ib(patch_id)%moving_ibm /= 0) then
1081 ! compute the linear velocity of the ghost point due to rotation
1082 radial_vector = physical_loc - [patch_ib(patch_id)%x_centroid, &
1083 patch_ib(patch_id)%y_centroid, patch_ib(patch_id)%z_centroid]
1084 call s_cross_product(matmul(patch_ib(patch_id)%rotation_matrix, patch_ib(patch_id)%angular_vel), radial_vector, rotation_velocity)
1085
1086 ! add only the component of the IB's motion that is normal to the surface
1087 vel_g = vel_g + sum((patch_ib(patch_id)%vel + rotation_velocity)*norm)*norm
1088 end if
1089 else
1090 if (patch_ib(patch_id)%moving_ibm == 0) then
1091 ! we know the object is not moving if moving_ibm is 0 (false)
1092 vel_g = 0._wp
1093 else
1094 ! get the vector that points from the centroid to the ghost
1095 radial_vector = physical_loc - [patch_ib(patch_id)%x_centroid, &
1096 patch_ib(patch_id)%y_centroid, patch_ib(patch_id)%z_centroid]
1097 ! convert the angular velocity from the inertial reference frame to the fluids frame, then convert to linear velocity
1098 call s_cross_product(matmul(patch_ib(patch_id)%rotation_matrix, patch_ib(patch_id)%angular_vel), radial_vector, rotation_velocity)
1099 do q = 1, 3
1100 ! if mibm is 1 or 2, then the boundary may be moving
1101 vel_g(q) = patch_ib(patch_id)%vel(q) ! add the linear velocity
1102 vel_g(q) = vel_g(q) + rotation_velocity(q) ! add the rotational velocity
1103 end do
1104 end if
1105 end if
1106
1107 ! Set momentum
1108
1109# 342 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1110#if defined(MFC_OpenACC)
1111# 342 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1112!$acc loop seq
1113# 342 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1114#elif defined(MFC_OpenMP)
1115# 342 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1116
1117# 342 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1118#endif
1119 do q = momxb, momxe
1120 q_cons_vf(q)%sf(j, k, l) = rho*vel_g(q - momxb + 1)
1121 dyn_pres = dyn_pres + q_cons_vf(q)%sf(j, k, l)* &
1122 vel_g(q - momxb + 1)/2._wp
1123 end do
1124
1125 ! Set continuity and adv vars
1126
1127# 350 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1128#if defined(MFC_OpenACC)
1129# 350 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1130!$acc loop seq
1131# 350 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1132#elif defined(MFC_OpenMP)
1133# 350 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1134
1135# 350 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1136#endif
1137 do q = 1, num_fluids
1138 q_cons_vf(q)%sf(j, k, l) = alpha_rho_ip(q)
1139 q_cons_vf(advxb + q - 1)%sf(j, k, l) = alpha_ip(q)
1140 end do
1141
1142 ! Set color function
1143 if (surface_tension) then
1144 q_cons_vf(c_idx)%sf(j, k, l) = c_ip
1145 end if
1146
1147 ! Set Energy
1148 if (bubbles_euler) then
1149 q_cons_vf(e_idx)%sf(j, k, l) = (1 - alpha_ip(1))*(gamma*pres_ip + pi_inf + dyn_pres)
1150 else
1151 q_cons_vf(e_idx)%sf(j, k, l) = gamma*pres_ip + pi_inf + dyn_pres
1152 end if
1153 ! Set bubble vars
1154 if (bubbles_euler .and. .not. qbmm) then
1155 call s_comp_n_from_prim(alpha_ip(1), r_ip, nbub, weight)
1156
1157# 370 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1158#if defined(MFC_OpenACC)
1159# 370 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1160!$acc loop seq
1161# 370 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1162#elif defined(MFC_OpenMP)
1163# 370 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1164
1165# 370 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1166#endif
1167 do q = 1, nb
1168 q_cons_vf(bubxb + (q - 1)*2)%sf(j, k, l) = nbub*r_ip(q)
1169 q_cons_vf(bubxb + (q - 1)*2 + 1)%sf(j, k, l) = nbub*v_ip(q)
1170 if (.not. polytropic) then
1171 q_cons_vf(bubxb + (q - 1)*4)%sf(j, k, l) = nbub*r_ip(q)
1172 q_cons_vf(bubxb + (q - 1)*4 + 1)%sf(j, k, l) = nbub*v_ip(q)
1173 q_cons_vf(bubxb + (q - 1)*4 + 2)%sf(j, k, l) = nbub*pb_ip(q)
1174 q_cons_vf(bubxb + (q - 1)*4 + 3)%sf(j, k, l) = nbub*mv_ip(q)
1175 end if
1176 end do
1177 end if
1178
1179 if (qbmm) then
1180
1181 nbub = nmom_ip(1)
1182
1183# 386 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1184#if defined(MFC_OpenACC)
1185# 386 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1186!$acc loop seq
1187# 386 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1188#elif defined(MFC_OpenMP)
1189# 386 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1190
1191# 386 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1192#endif
1193 do q = 1, nb*nmom
1194 q_cons_vf(bubxb + q - 1)%sf(j, k, l) = nbub*nmom_ip(q)
1195 end do
1196
1197
1198# 391 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1199#if defined(MFC_OpenACC)
1200# 391 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1201!$acc loop seq
1202# 391 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1203#elif defined(MFC_OpenMP)
1204# 391 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1205
1206# 391 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1207#endif
1208 do q = 1, nb
1209 q_cons_vf(bubxb + (q - 1)*nmom)%sf(j, k, l) = nbub
1210 end do
1211
1212 if (.not. polytropic) then
1213
1214# 397 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1215#if defined(MFC_OpenACC)
1216# 397 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1217!$acc loop seq
1218# 397 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1219#elif defined(MFC_OpenMP)
1220# 397 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1221
1222# 397 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1223#endif
1224 do q = 1, nb
1225
1226# 399 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1227#if defined(MFC_OpenACC)
1228# 399 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1229!$acc loop seq
1230# 399 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1231#elif defined(MFC_OpenMP)
1232# 399 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1233
1234# 399 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1235#endif
1236 do r = 1, nnode
1237 pb_in(j, k, l, r, q) = presb_ip((q - 1)*nnode + r)
1238 mv_in(j, k, l, r, q) = massv_ip((q - 1)*nnode + r)
1239 end do
1240 end do
1241 end if
1242 end if
1243
1244 if (model_eqns == 3) then
1245
1246# 409 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1247#if defined(MFC_OpenACC)
1248# 409 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1249!$acc loop seq
1250# 409 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1251#elif defined(MFC_OpenMP)
1252# 409 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1253
1254# 409 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1255#endif
1256 do q = intxb, intxe
1257 q_cons_vf(q)%sf(j, k, l) = alpha_ip(q - intxb + 1)*(gammas(q - intxb + 1)*pres_ip &
1258 + pi_infs(q - intxb + 1))
1259 end do
1260 end if
1261 end do
1262
1263# 416 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1264
1265# 416 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1266#if defined(MFC_OpenACC)
1267# 416 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1268!$acc end parallel loop
1269# 416 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1270#elif defined(MFC_OpenMP)
1271# 416 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1272
1273# 416 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1274
1275# 416 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1276!$omp end target teams loop
1277# 416 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1278#endif
1279# 416 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1280
1281 end if
1282
1283 !Correct the state of the inner points in IBs
1284 if (num_inner_gps > 0) then
1285
1286# 421 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1287
1288# 421 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1289#if defined(MFC_OpenACC)
1290# 421 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1291!$acc parallel loop gang vector default(present) private(i, physical_loc, dyn_pres, alpha_rho_IP, alpha_IP, vel_g, rho, gamma, pi_inf, Re_K, innerp, j, k, l, q)
1292# 421 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1293#elif defined(MFC_OpenMP)
1294# 421 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1295
1296# 421 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1297
1298# 421 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1299
1300# 421 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1301!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, physical_loc, dyn_pres, alpha_rho_IP, alpha_IP, vel_g, rho, gamma, pi_inf, Re_K, innerp, j, k, l, q)
1302# 421 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1303#endif
1304# 421 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1305
1306 do i = 1, num_inner_gps
1307
1308 innerp = inner_points(i)
1309 j = innerp%loc(1)
1310 k = innerp%loc(2)
1311 l = innerp%loc(3)
1312
1313
1314# 429 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1315#if defined(MFC_OpenACC)
1316# 429 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1317!$acc loop seq
1318# 429 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1319#elif defined(MFC_OpenMP)
1320# 429 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1321
1322# 429 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1323#endif
1324 do q = momxb, momxe
1325 q_cons_vf(q)%sf(j, k, l) = 0._wp
1326 end do
1327 end do
1328
1329# 434 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1330
1331# 434 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1332#if defined(MFC_OpenACC)
1333# 434 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1334!$acc end parallel loop
1335# 434 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1336#elif defined(MFC_OpenMP)
1337# 434 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1338
1339# 434 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1340
1341# 434 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1342!$omp end target teams loop
1343# 434 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1344#endif
1345# 434 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1346
1347 end if
1348
1349 end subroutine s_ibm_correct_state
1350
1351 !> Function that computes the image points for each ghost point
1352 !! @param ghost_points_in Ghost Points
1353 impure subroutine s_compute_image_points(ghost_points_in)
1354
1355 type(ghost_point), dimension(num_gps), intent(INOUT) :: ghost_points_in
1356
1357 real(wp) :: dist
1358 real(wp), dimension(3) :: norm
1359 real(wp), dimension(3) :: physical_loc
1360 real(wp) :: temp_loc
1361 real(wp), pointer, dimension(:) :: s_cc => null()
1362 integer :: bound
1363 type(ghost_point) :: gp
1364
1365 integer :: q, dim !< Iterator variables
1366 integer :: i, j, k, l !< Location indexes
1367 integer :: patch_id !< IB Patch ID
1368 integer :: dir
1369 integer :: index
1370
1371 do q = 1, num_gps
1372 gp = ghost_points_in(q)
1373 i = gp%loc(1)
1374 j = gp%loc(2)
1375 k = gp%loc(3)
1376
1377 ! Calculate physical location of ghost point
1378 if (p > 0) then
1379 physical_loc = [x_cc(i), y_cc(j), z_cc(k)]
1380 else
1381 physical_loc = [x_cc(i), y_cc(j), 0._wp]
1382 end if
1383
1384 ! Calculate and store the precise location of the image point
1385 patch_id = gp%ib_patch_id
1386 dist = abs(real(gp%levelset, kind=wp))
1387 norm(:) = gp%levelset_norm
1388 ghost_points_in(q)%ip_loc(:) = physical_loc(:) + 2*dist*norm(:)
1389
1390 ! Find the closest grid point to the image point
1391 do dim = 1, num_dims
1392
1393 ! s_cc points to the dim array we need
1394 if (dim == 1) then
1395 s_cc => x_cc
1396 bound = m + buff_size - 1
1397 elseif (dim == 2) then
1398 s_cc => y_cc
1399 bound = n + buff_size - 1
1400 else
1401 s_cc => z_cc
1402 bound = p + buff_size - 1
1403 end if
1404
1405 if (f_approx_equal(norm(dim), 0._wp)) then
1406 ! if the ghost point is almost equal to a cell location, we set it equal and continue
1407 ghost_points_in(q)%ip_grid(dim) = ghost_points_in(q)%loc(dim)
1408 else
1409 if (norm(dim) > 0) then
1410 dir = 1
1411 else
1412 dir = -1
1413 end if
1414
1415 index = ghost_points_in(q)%loc(dim)
1416 temp_loc = ghost_points_in(q)%ip_loc(dim)
1417 do while ((temp_loc < s_cc(index) &
1418 .or. temp_loc > s_cc(index + 1)))
1419 index = index + dir
1420 if (index < -buff_size .or. index > bound) then
1421 print *, "A required image point is not located in this computational domain."
1422 print *, "Ghost Point is located at ", [x_cc(i), y_cc(j), z_cc(k)], " while moving in dimension ", dim
1423 print *, "We are searching for image point at ", ghost_points_in(q)%ip_loc(:)
1424 print *, "We can only support points located inside the box from ", [x_cc(-buff_size), y_cc(-buff_size), z_cc(-buff_size)]
1425 print *, "To ", [x_cc(m + buff_size - 1), y_cc(n + buff_size - 1), z_cc(p + buff_size - 1)]
1426 print *, "Image point is located approximately ", (ghost_points_in(q)%loc(dim) - ghost_points_in(q)%ip_loc(dim))/(s_cc(1) - s_cc(0)), " grid cells away"
1427 print *, "Levelset ", dist, " and Norm: ", norm(:)
1428 print *, "A short term fix may include increasing buff_size further in m_helper_basic (currently set to a minimum of 10)"
1429 error stop "Ghost Point and Image Point on Different Processors"
1430 end if
1431 end do
1432 ghost_points_in(q)%ip_grid(dim) = index
1433 if (ghost_points_in(q)%DB(dim) == -1) then
1434 ghost_points_in(q)%ip_grid(dim) = ghost_points_in(q)%loc(dim) + 1
1435 else if (ghost_points_in(q)%DB(dim) == 1) then
1436 ghost_points_in(q)%ip_grid(dim) = ghost_points_in(q)%loc(dim) - 1
1437 end if
1438 end if
1439 end do
1440 end do
1441
1442 end subroutine s_compute_image_points
1443
1444 !> Function that finds the number of ghost points, used for allocating
1445 !! memory.
1446 subroutine s_find_num_ghost_points(num_gps_out, num_inner_gps_out)
1447
1448 integer, intent(out) :: num_gps_out
1449 integer, intent(out) :: num_inner_gps_out
1450
1451 integer, dimension(2*gp_layers + 1, 2*gp_layers + 1) &
1452 :: subsection_2d
1453 integer, dimension(2*gp_layers + 1, 2*gp_layers + 1, 2*gp_layers + 1) &
1454 :: subsection_3d
1455 integer :: i, j, k!< Iterator variables
1456
1457 num_gps_out = 0
1458 num_inner_gps_out = 0
1459
1460 do i = 0, m
1461 do j = 0, n
1462 if (p == 0) then
1463 if (ib_markers%sf(i, j, 0) /= 0) then
1464 subsection_2d = ib_markers%sf( &
1465 i - gp_layers:i + gp_layers, &
1466 j - gp_layers:j + gp_layers, 0)
1467 if (any(subsection_2d == 0)) then
1468 num_gps_out = num_gps_out + 1
1469 else
1470 num_inner_gps_out = num_inner_gps_out + 1
1471 end if
1472 end if
1473 else
1474 do k = 0, p
1475 if (ib_markers%sf(i, j, k) /= 0) then
1476 subsection_3d = ib_markers%sf( &
1477 i - gp_layers:i + gp_layers, &
1478 j - gp_layers:j + gp_layers, &
1479 k - gp_layers:k + gp_layers)
1480 if (any(subsection_3d == 0)) then
1481 num_gps_out = num_gps_out + 1
1482 else
1483 num_inner_gps_out = num_inner_gps_out + 1
1484 end if
1485 end if
1486 end do
1487 end if
1488 end do
1489 end do
1490
1491 end subroutine s_find_num_ghost_points
1492
1493 !> Function that finds the ghost points
1494 subroutine s_find_ghost_points(ghost_points_in, inner_points_in)
1495
1496 type(ghost_point), dimension(num_gps), intent(INOUT) :: ghost_points_in
1497 type(ghost_point), dimension(num_inner_gps), intent(INOUT) :: inner_points_in
1498 integer, dimension(2*gp_layers + 1, 2*gp_layers + 1) &
1499 :: subsection_2d
1500 integer, dimension(2*gp_layers + 1, 2*gp_layers + 1, 2*gp_layers + 1) &
1501 :: subsection_3d
1502 integer :: i, j, k !< Iterator variables
1503 integer :: count, count_i
1504 integer :: patch_id
1505
1506 count = 1
1507 count_i = 1
1508
1509 do i = 0, m
1510 do j = 0, n
1511 if (p == 0) then
1512 ! 2D
1513 if (ib_markers%sf(i, j, 0) /= 0) then
1514 subsection_2d = ib_markers%sf( &
1515 i - gp_layers:i + gp_layers, &
1516 j - gp_layers:j + gp_layers, 0)
1517 if (any(subsection_2d == 0)) then
1518 ghost_points_in(count)%loc = [i, j, 0]
1519 patch_id = ib_markers%sf(i, j, 0)
1520 ghost_points_in(count)%ib_patch_id = &
1521 patch_id
1522
1523 ghost_points_in(count)%slip = patch_ib(patch_id)%slip
1524 ! ghost_points(count)%rank = proc_rank
1525
1526 if ((x_cc(i) - dx(i)) < x_domain%beg) then
1527 ghost_points_in(count)%DB(1) = -1
1528 else if ((x_cc(i) + dx(i)) > x_domain%end) then
1529 ghost_points_in(count)%DB(1) = 1
1530 else
1531 ghost_points_in(count)%DB(1) = 0
1532 end if
1533
1534 if ((y_cc(j) - dy(j)) < y_domain%beg) then
1535 ghost_points_in(count)%DB(2) = -1
1536 else if ((y_cc(j) + dy(j)) > y_domain%end) then
1537 ghost_points_in(count)%DB(2) = 1
1538 else
1539 ghost_points_in(count)%DB(2) = 0
1540 end if
1541
1542 count = count + 1
1543
1544 else
1545 inner_points_in(count_i)%loc = [i, j, 0]
1546 patch_id = ib_markers%sf(i, j, 0)
1547 inner_points_in(count_i)%ib_patch_id = &
1548 patch_id
1549 inner_points_in(count_i)%slip = patch_ib(patch_id)%slip
1550 count_i = count_i + 1
1551
1552 end if
1553 end if
1554 else
1555 ! 3D
1556 do k = 0, p
1557 if (ib_markers%sf(i, j, k) /= 0) then
1558 subsection_3d = ib_markers%sf( &
1559 i - gp_layers:i + gp_layers, &
1560 j - gp_layers:j + gp_layers, &
1561 k - gp_layers:k + gp_layers)
1562 if (any(subsection_3d == 0)) then
1563 ghost_points_in(count)%loc = [i, j, k]
1564 patch_id = ib_markers%sf(i, j, k)
1565 ghost_points_in(count)%ib_patch_id = &
1566 ib_markers%sf(i, j, k)
1567 ghost_points_in(count)%slip = patch_ib(patch_id)%slip
1568
1569 if ((x_cc(i) - dx(i)) < x_domain%beg) then
1570 ghost_points_in(count)%DB(1) = -1
1571 else if ((x_cc(i) + dx(i)) > x_domain%end) then
1572 ghost_points_in(count)%DB(1) = 1
1573 else
1574 ghost_points_in(count)%DB(1) = 0
1575 end if
1576
1577 if ((y_cc(j) - dy(j)) < y_domain%beg) then
1578 ghost_points_in(count)%DB(2) = -1
1579 else if ((y_cc(j) + dy(j)) > y_domain%end) then
1580 ghost_points_in(count)%DB(2) = 1
1581 else
1582 ghost_points_in(count)%DB(2) = 0
1583 end if
1584
1585 if ((z_cc(k) - dz(k)) < z_domain%beg) then
1586 ghost_points_in(count)%DB(3) = -1
1587 else if ((z_cc(k) + dz(k)) > z_domain%end) then
1588 ghost_points_in(count)%DB(3) = 1
1589 else
1590 ghost_points_in(count)%DB(3) = 0
1591 end if
1592
1593 count = count + 1
1594 else
1595 inner_points_in(count_i)%loc = [i, j, k]
1596 patch_id = ib_markers%sf(i, j, k)
1597 inner_points_in(count_i)%ib_patch_id = &
1598 ib_markers%sf(i, j, k)
1599 inner_points_in(count_i)%slip = patch_ib(patch_id)%slip
1600
1601 count_i = count_i + 1
1602 end if
1603 end if
1604 end do
1605 end if
1606 end do
1607 end do
1608
1609 end subroutine s_find_ghost_points
1610
1611 !> Function that computes the interpolation coefficients of image points
1612 subroutine s_compute_interpolation_coeffs(ghost_points_in)
1613
1614 type(ghost_point), dimension(num_gps), intent(INOUT) :: ghost_points_in
1615
1616 real(wp), dimension(2, 2, 2) :: dist
1617 real(wp), dimension(2, 2, 2) :: alpha
1618 real(wp), dimension(2, 2, 2) :: interp_coeffs
1619 real(wp) :: buf
1620 real(wp), dimension(2, 2, 2) :: eta
1621 type(ghost_point) :: gp
1622 integer :: i !< Iterator variables
1623 integer :: i1, i2, j1, j2, k1, k2 !< Grid indexes
1624 integer :: patch_id
1625
1626 ! 2D
1627 if (p <= 0) then
1628 do i = 1, num_gps
1629 gp = ghost_points_in(i)
1630 ! Get the interpolation points
1631 i1 = gp%ip_grid(1); i2 = i1 + 1
1632 j1 = gp%ip_grid(2); j2 = j1 + 1
1633
1634 dist = 0._wp
1635 buf = 1._wp
1636 dist(1, 1, 1) = sqrt( &
1637 (x_cc(i1) - gp%ip_loc(1))**2 + &
1638 (y_cc(j1) - gp%ip_loc(2))**2)
1639 dist(2, 1, 1) = sqrt( &
1640 (x_cc(i2) - gp%ip_loc(1))**2 + &
1641 (y_cc(j1) - gp%ip_loc(2))**2)
1642 dist(1, 2, 1) = sqrt( &
1643 (x_cc(i1) - gp%ip_loc(1))**2 + &
1644 (y_cc(j2) - gp%ip_loc(2))**2)
1645 dist(2, 2, 1) = sqrt( &
1646 (x_cc(i2) - gp%ip_loc(1))**2 + &
1647 (y_cc(j2) - gp%ip_loc(2))**2)
1648
1649 interp_coeffs = 0._wp
1650
1651 if (dist(1, 1, 1) <= 1.e-16_wp) then
1652 interp_coeffs(1, 1, 1) = 1._wp
1653 else if (dist(2, 1, 1) <= 1.e-16_wp) then
1654 interp_coeffs(2, 1, 1) = 1._wp
1655 else if (dist(1, 2, 1) <= 1.e-16_wp) then
1656 interp_coeffs(1, 2, 1) = 1._wp
1657 else if (dist(2, 2, 1) <= 1.e-16_wp) then
1658 interp_coeffs(2, 2, 1) = 1._wp
1659 else
1660 eta(:, :, 1) = 1._wp/dist(:, :, 1)**2
1661 alpha = 1._wp
1662 patch_id = gp%ib_patch_id
1663 if (ib_markers%sf(i1, j1, 0) /= 0) alpha(1, 1, 1) = 0._wp
1664 if (ib_markers%sf(i2, j1, 0) /= 0) alpha(2, 1, 1) = 0._wp
1665 if (ib_markers%sf(i1, j2, 0) /= 0) alpha(1, 2, 1) = 0._wp
1666 if (ib_markers%sf(i2, j2, 0) /= 0) alpha(2, 2, 1) = 0._wp
1667 buf = sum(alpha(:, :, 1)*eta(:, :, 1))
1668 if (buf > 0._wp) then
1669 interp_coeffs(:, :, 1) = alpha(:, :, 1)*eta(:, :, 1)/buf
1670 else
1671 buf = sum(eta(:, :, 1))
1672 interp_coeffs(:, :, 1) = eta(:, :, 1)/buf
1673 end if
1674 end if
1675
1676 ghost_points_in(i)%interp_coeffs = interp_coeffs
1677 end do
1678
1679 else
1680 do i = 1, num_gps
1681 gp = ghost_points_in(i)
1682 ! Get the interpolation points
1683 i1 = gp%ip_grid(1); i2 = i1 + 1
1684 j1 = gp%ip_grid(2); j2 = j1 + 1
1685 k1 = gp%ip_grid(3); k2 = k1 + 1
1686
1687 ! Get interpolation weights (Chaudhuri et al. 2011, JCP)
1688 dist(1, 1, 1) = sqrt( &
1689 (x_cc(i1) - gp%ip_loc(1))**2 + &
1690 (y_cc(j1) - gp%ip_loc(2))**2 + &
1691 (z_cc(k1) - gp%ip_loc(3))**2)
1692 dist(2, 1, 1) = sqrt( &
1693 (x_cc(i2) - gp%ip_loc(1))**2 + &
1694 (y_cc(j1) - gp%ip_loc(2))**2 + &
1695 (z_cc(k1) - gp%ip_loc(3))**2)
1696 dist(1, 2, 1) = sqrt( &
1697 (x_cc(i1) - gp%ip_loc(1))**2 + &
1698 (y_cc(j2) - gp%ip_loc(2))**2 + &
1699 (z_cc(k1) - gp%ip_loc(3))**2)
1700 dist(2, 2, 1) = sqrt( &
1701 (x_cc(i2) - gp%ip_loc(1))**2 + &
1702 (y_cc(j2) - gp%ip_loc(2))**2 + &
1703 (z_cc(k1) - gp%ip_loc(3))**2)
1704 dist(1, 1, 2) = sqrt( &
1705 (x_cc(i1) - gp%ip_loc(1))**2 + &
1706 (y_cc(j1) - gp%ip_loc(2))**2 + &
1707 (z_cc(k2) - gp%ip_loc(3))**2)
1708 dist(2, 1, 2) = sqrt( &
1709 (x_cc(i2) - gp%ip_loc(1))**2 + &
1710 (y_cc(j1) - gp%ip_loc(2))**2 + &
1711 (z_cc(k2) - gp%ip_loc(3))**2)
1712 dist(1, 2, 2) = sqrt( &
1713 (x_cc(i1) - gp%ip_loc(1))**2 + &
1714 (y_cc(j2) - gp%ip_loc(2))**2 + &
1715 (z_cc(k2) - gp%ip_loc(3))**2)
1716 dist(2, 2, 2) = sqrt( &
1717 (x_cc(i2) - gp%ip_loc(1))**2 + &
1718 (y_cc(j2) - gp%ip_loc(2))**2 + &
1719 (z_cc(k2) - gp%ip_loc(3))**2)
1720 interp_coeffs = 0._wp
1721 buf = 1._wp
1722 if (dist(1, 1, 1) <= 1.e-16_wp) then
1723 interp_coeffs(1, 1, 1) = 1._wp
1724 else if (dist(2, 1, 1) <= 1.e-16_wp) then
1725 interp_coeffs(2, 1, 1) = 1._wp
1726 else if (dist(1, 2, 1) <= 1.e-16_wp) then
1727 interp_coeffs(1, 2, 1) = 1._wp
1728 else if (dist(2, 2, 1) <= 1.e-16_wp) then
1729 interp_coeffs(2, 2, 1) = 1._wp
1730 else if (dist(1, 1, 2) <= 1.e-16_wp) then
1731 interp_coeffs(1, 1, 2) = 1._wp
1732 else if (dist(2, 1, 2) <= 1.e-16_wp) then
1733 interp_coeffs(2, 1, 2) = 1._wp
1734 else if (dist(1, 2, 2) <= 1.e-16_wp) then
1735 interp_coeffs(1, 2, 2) = 1._wp
1736 else if (dist(2, 2, 2) <= 1.e-16_wp) then
1737 interp_coeffs(2, 2, 2) = 1._wp
1738 else
1739 eta = 1._wp/dist**2
1740 alpha = 1._wp
1741 if (ib_markers%sf(i1, j1, k1) /= 0) alpha(1, 1, 1) = 0._wp
1742 if (ib_markers%sf(i2, j1, k1) /= 0) alpha(2, 1, 1) = 0._wp
1743 if (ib_markers%sf(i1, j2, k1) /= 0) alpha(1, 2, 1) = 0._wp
1744 if (ib_markers%sf(i2, j2, k1) /= 0) alpha(2, 2, 1) = 0._wp
1745 if (ib_markers%sf(i1, j1, k2) /= 0) alpha(1, 1, 2) = 0._wp
1746 if (ib_markers%sf(i2, j1, k2) /= 0) alpha(2, 1, 2) = 0._wp
1747 if (ib_markers%sf(i1, j2, k2) /= 0) alpha(1, 2, 2) = 0._wp
1748 if (ib_markers%sf(i2, j2, k2) /= 0) alpha(2, 2, 2) = 0._wp
1749 buf = sum(alpha*eta)
1750 if (buf > 0._wp) then
1751 interp_coeffs = alpha*eta/buf
1752 else
1753 buf = sum(eta)
1754 interp_coeffs = eta/buf
1755 end if
1756 end if
1757
1758 ghost_points_in(i)%interp_coeffs = interp_coeffs
1759 end do
1760 end if
1761
1762 end subroutine s_compute_interpolation_coeffs
1763
1764 !> Function that uses the interpolation coefficients and the current state
1765 !! at the cell centers in order to estimate the state at the image point
1766 !! @param gp Ghost point data structure
1767 !> @brief Interpolates primitive variables from the fluid domain to a ghost point's image point using bilinear or trilinear interpolation.
1768 !! @param alpha_rho_IP Partial density at image point
1769 !! @param alpha_IP Volume fraction at image point
1770 !! @param pres_IP Pressure at image point
1771 !! @param vel_IP Velocity at image point
1772 !! @param c_IP Speed of sound at image point
1773 !! @param r_IP Bubble radius at image point
1774 !! @param v_IP Bubble radial velocity at image point
1775 !! @param pb_IP Bubble pressure at image point
1776 !! @param mv_IP Bubble vapor mass at image point
1777 !! @param nmom_IP Bubble moment at image point
1778 !! @param pb_in Internal bubble pressure array
1779 !! @param mv_in Mass of vapor in bubble array
1780 !! @param presb_IP Bubble node pressure at image point
1781 !! @param massv_IP Bubble node vapor mass at image point
1782 subroutine s_interpolate_image_point(q_prim_vf, gp, alpha_rho_IP, alpha_IP, &
1783 pres_IP, vel_IP, c_IP, r_IP, v_IP, pb_IP, &
1784 mv_IP, nmom_IP, pb_in, mv_in, presb_IP, massv_IP)
1785
1786# 873 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1787#if MFC_OpenACC
1788# 873 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1789!$acc routine seq
1790# 873 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1791#elif MFC_OpenMP
1792# 873 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1793
1794# 873 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1795
1796# 873 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1797!$omp declare target device_type(any)
1798# 873 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1799#endif
1800 type(scalar_field), &
1801 dimension(sys_size), &
1802 intent(IN) :: q_prim_vf !< Primitive Variables
1803
1804 real(stp), optional, dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(IN) :: pb_in, mv_in
1805
1806 type(ghost_point), intent(IN) :: gp
1807 real(wp), intent(INOUT) :: pres_ip
1808 real(wp), dimension(3), intent(INOUT) :: vel_ip
1809 real(wp), intent(INOUT) :: c_ip
1810# 887 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1811 real(wp), dimension(num_fluids), intent(INOUT) :: alpha_ip, alpha_rho_ip
1812# 889 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1813 real(wp), optional, dimension(:), intent(INOUT) :: r_ip, v_ip, pb_ip, mv_ip
1814 real(wp), optional, dimension(:), intent(INOUT) :: nmom_ip
1815 real(wp), optional, dimension(:), intent(INOUT) :: presb_ip, massv_ip
1816
1817 integer :: i, j, k, l, q !< Iterator variables
1818 integer :: i1, i2, j1, j2, k1, k2 !< Iterator variables
1819 real(wp) :: coeff
1820
1821 i1 = gp%ip_grid(1); i2 = i1 + 1
1822 j1 = gp%ip_grid(2); j2 = j1 + 1
1823 k1 = gp%ip_grid(3); k2 = k1 + 1
1824
1825 if (p == 0) then
1826 k1 = 0
1827 k2 = 0
1828 end if
1829
1830 alpha_rho_ip = 0._wp
1831 alpha_ip = 0._wp
1832 pres_ip = 0._wp
1833 vel_ip = 0._wp
1834
1835 if (surface_tension) c_ip = 0._wp
1836
1837 if (bubbles_euler) then
1838 r_ip = 0._wp
1839 v_ip = 0._wp
1840 if (.not. polytropic) then
1841 mv_ip = 0._wp
1842 pb_ip = 0._wp
1843 end if
1844 end if
1845
1846 if (qbmm) then
1847 nmom_ip = 0._wp
1848 if (.not. polytropic) then
1849 presb_ip = 0._wp
1850 massv_ip = 0._wp
1851 end if
1852 end if
1853
1854
1855# 930 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1856#if defined(MFC_OpenACC)
1857# 930 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1858!$acc loop seq
1859# 930 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1860#elif defined(MFC_OpenMP)
1861# 930 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1862
1863# 930 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1864#endif
1865 do i = i1, i2
1866
1867# 932 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1868#if defined(MFC_OpenACC)
1869# 932 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1870!$acc loop seq
1871# 932 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1872#elif defined(MFC_OpenMP)
1873# 932 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1874
1875# 932 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1876#endif
1877 do j = j1, j2
1878
1879# 934 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1880#if defined(MFC_OpenACC)
1881# 934 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1882!$acc loop seq
1883# 934 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1884#elif defined(MFC_OpenMP)
1885# 934 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1886
1887# 934 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1888#endif
1889 do k = k1, k2
1890
1891 coeff = gp%interp_coeffs(i - i1 + 1, j - j1 + 1, k - k1 + 1)
1892
1893 pres_ip = pres_ip + coeff* &
1894 q_prim_vf(e_idx)%sf(i, j, k)
1895
1896
1897# 942 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1898#if defined(MFC_OpenACC)
1899# 942 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1900!$acc loop seq
1901# 942 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1902#elif defined(MFC_OpenMP)
1903# 942 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1904
1905# 942 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1906#endif
1907 do q = momxb, momxe
1908 vel_ip(q + 1 - momxb) = vel_ip(q + 1 - momxb) + coeff* &
1909 q_prim_vf(q)%sf(i, j, k)
1910 end do
1911
1912
1913# 948 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1914#if defined(MFC_OpenACC)
1915# 948 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1916!$acc loop seq
1917# 948 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1918#elif defined(MFC_OpenMP)
1919# 948 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1920
1921# 948 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1922#endif
1923 do l = contxb, contxe
1924 alpha_rho_ip(l) = alpha_rho_ip(l) + coeff* &
1925 q_prim_vf(l)%sf(i, j, k)
1926 alpha_ip(l) = alpha_ip(l) + coeff* &
1927 q_prim_vf(advxb + l - 1)%sf(i, j, k)
1928 end do
1929
1930 if (surface_tension) then
1931 c_ip = c_ip + coeff*q_prim_vf(c_idx)%sf(i, j, k)
1932 end if
1933
1934 if (bubbles_euler .and. .not. qbmm) then
1935
1936# 961 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1937#if defined(MFC_OpenACC)
1938# 961 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1939!$acc loop seq
1940# 961 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1941#elif defined(MFC_OpenMP)
1942# 961 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1943
1944# 961 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1945#endif
1946 do l = 1, nb
1947 if (polytropic) then
1948 r_ip(l) = r_ip(l) + coeff*q_prim_vf(bubxb + (l - 1)*2)%sf(i, j, k)
1949 v_ip(l) = v_ip(l) + coeff*q_prim_vf(bubxb + 1 + (l - 1)*2)%sf(i, j, k)
1950 else
1951 r_ip(l) = r_ip(l) + coeff*q_prim_vf(bubxb + (l - 1)*4)%sf(i, j, k)
1952 v_ip(l) = v_ip(l) + coeff*q_prim_vf(bubxb + 1 + (l - 1)*4)%sf(i, j, k)
1953 pb_ip(l) = pb_ip(l) + coeff*q_prim_vf(bubxb + 2 + (l - 1)*4)%sf(i, j, k)
1954 mv_ip(l) = mv_ip(l) + coeff*q_prim_vf(bubxb + 3 + (l - 1)*4)%sf(i, j, k)
1955 end if
1956 end do
1957 end if
1958
1959 if (qbmm) then
1960 do l = 1, nb*nmom
1961 nmom_ip(l) = nmom_ip(l) + coeff*q_prim_vf(bubxb - 1 + l)%sf(i, j, k)
1962 end do
1963 if (.not. polytropic) then
1964 do q = 1, nb
1965 do l = 1, nnode
1966 presb_ip((q - 1)*nnode + l) = presb_ip((q - 1)*nnode + l) + &
1967 coeff*real(pb_in(i, j, k, l, q), kind=wp)
1968 massv_ip((q - 1)*nnode + l) = massv_ip((q - 1)*nnode + l) + &
1969 coeff*real(mv_in(i, j, k, l, q), kind=wp)
1970 end do
1971 end do
1972 end if
1973
1974 end if
1975
1976 end do
1977 end do
1978 end do
1979
1980 end subroutine s_interpolate_image_point
1981
1982 !> Resets the current indexes of immersed boundaries and replaces them after updating
1983 !> the position of each moving immersed boundary
1984 impure subroutine s_update_mib(num_ibs)
1985
1986 integer, intent(in) :: num_ibs
1987
1988 integer :: i, ierr
1989
1990 ! Clears the existing immersed boundary indices
1991 ib_markers%sf = 0._wp
1992
1993 ! recalulcate the rotation matrix based upon the new angles
1994 do i = 1, num_ibs
1995 if (patch_ib(i)%moving_ibm /= 0) then
1996 call s_update_ib_rotation_matrix(i)
1997 end if
1998 end do
1999
2000
2001# 1016 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2002#if defined(MFC_OpenACC)
2003# 1016 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2004!$acc update device(patch_ib)
2005# 1016 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2006#elif defined(MFC_OpenMP)
2007# 1016 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2008!$omp target update to(patch_ib)
2009# 1016 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2010#endif
2011
2012 ! recompute the new ib_patch locations and broadcast them.
2013 call s_apply_ib_patches(ib_markers%sf(0:m, 0:n, 0:p))
2014
2015# 1020 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2016#if defined(MFC_OpenACC)
2017# 1020 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2018!$acc update device(ib_markers%sf)
2019# 1020 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2020#elif defined(MFC_OpenMP)
2021# 1020 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2022!$omp target update to(ib_markers%sf)
2023# 1020 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2024#endif
2025 call s_populate_ib_buffers()
2026
2027# 1022 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2028#if defined(MFC_OpenACC)
2029# 1022 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2030!$acc update host(ib_markers%sf)
2031# 1022 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2032#elif defined(MFC_OpenMP)
2033# 1022 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2034!$omp target update from(ib_markers%sf)
2035# 1022 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2036#endif
2037
2038 ! recalculate the ghost point locations and coefficients
2039 call s_find_num_ghost_points(num_gps, num_inner_gps)
2040
2041# 1026 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2042#if defined(MFC_OpenACC)
2043# 1026 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2044!$acc update device(num_gps, num_inner_gps)
2045# 1026 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2046#elif defined(MFC_OpenMP)
2047# 1026 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2048!$omp target update to(num_gps, num_inner_gps)
2049# 1026 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2050#endif
2051
2052 call s_find_ghost_points(ghost_points, inner_points)
2053
2054# 1029 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2055#if defined(MFC_OpenACC)
2056# 1029 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2057!$acc update device(ghost_points, inner_points)
2058# 1029 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2059#elif defined(MFC_OpenMP)
2060# 1029 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2061!$omp target update to(ghost_points, inner_points)
2062# 1029 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2063#endif
2064
2065 call s_apply_levelset(ghost_points, num_gps)
2066
2067# 1032 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2068#if defined(MFC_OpenACC)
2069# 1032 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2070!$acc update device(ghost_points)
2071# 1032 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2072#elif defined(MFC_OpenMP)
2073# 1032 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2074!$omp target update to(ghost_points)
2075# 1032 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2076#endif
2077
2078 call s_compute_image_points(ghost_points)
2079
2080# 1035 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2081#if defined(MFC_OpenACC)
2082# 1035 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2083!$acc update device(ghost_points)
2084# 1035 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2085#elif defined(MFC_OpenMP)
2086# 1035 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2087!$omp target update to(ghost_points)
2088# 1035 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2089#endif
2090
2091 call s_compute_interpolation_coeffs(ghost_points)
2092
2093# 1038 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2094#if defined(MFC_OpenACC)
2095# 1038 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2096!$acc update device(ghost_points)
2097# 1038 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2098#elif defined(MFC_OpenMP)
2099# 1038 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2100!$omp target update to(ghost_points)
2101# 1038 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2102#endif
2103
2104 end subroutine s_update_mib
2105
2106 !> @brief Computes pressure and viscous forces and torques on immersed bodies via a volume integration method.
2107 subroutine s_compute_ib_forces(q_prim_vf, fluid_pp)
2108
2109 ! real(wp), dimension(idwbuff(1)%beg:idwbuff(1)%end, &
2110 ! idwbuff(2)%beg:idwbuff(2)%end, &
2111 ! idwbuff(3)%beg:idwbuff(3)%end), intent(in) :: pressure
2112 type(scalar_field), dimension(1:sys_size), intent(in) :: q_prim_vf
2113 type(physical_parameters), dimension(1:num_fluids), intent(in) :: fluid_pp
2114
2115 integer :: gp_id, i, j, k, l, q, ib_idx, fluid_idx
2116 real(wp), dimension(num_ibs, 3) :: forces, torques
2117 real(wp), dimension(1:3, 1:3) :: viscous_stress_div, viscous_stress_div_1, viscous_stress_div_2, viscous_cross_1, viscous_cross_2 ! viscous stress tensor with temp vectors to hold divergence calculations
2118 real(wp), dimension(1:3) :: local_force_contribution, radial_vector, local_torque_contribution, vel
2119 real(wp) :: cell_volume, dx, dy, dz, dynamic_viscosity
2120# 1059 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2121 real(wp), dimension(num_fluids) :: dynamic_viscosities
2122# 1061 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2123 forces = 0._wp
2124 torques = 0._wp
2125
2126 if (viscous) then
2127 do fluid_idx = 1, num_fluids
2128 if (fluid_pp(fluid_idx)%Re(1) /= 0._wp) then
2129 dynamic_viscosities(fluid_idx) = 1._wp/fluid_pp(fluid_idx)%Re(1)
2130 else
2131 dynamic_viscosities(fluid_idx) = 0._wp
2132 end if
2133 end do
2134 end if
2135
2136
2137# 1074 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2138
2139# 1074 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2140#if defined(MFC_OpenACC)
2141# 1074 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2142!$acc parallel loop collapse(3) gang vector default(present) private(ib_idx, fluid_idx, radial_vector, local_force_contribution, cell_volume, local_torque_contribution, dynamic_viscosity, viscous_stress_div, viscous_stress_div_1, viscous_stress_div_2, viscous_cross_1, viscous_cross_2, dx, dy, dz) copy(forces, torques) copyin(ib_markers, patch_ib, dynamic_viscosities)
2143# 1074 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2144#elif defined(MFC_OpenMP)
2145# 1074 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2146
2147# 1074 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2148
2149# 1074 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2150
2151# 1074 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2152!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(ib_idx, fluid_idx, radial_vector, local_force_contribution, cell_volume, local_torque_contribution, dynamic_viscosity, viscous_stress_div, viscous_stress_div_1, viscous_stress_div_2, viscous_cross_1, viscous_cross_2, dx, dy, dz) map(tofrom:forces, torques) map(to:ib_markers, patch_ib, dynamic_viscosities)
2153# 1074 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2154#endif
2155# 1074 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2156
2157 do i = 0, m
2158 do j = 0, n
2159 do k = 0, p
2160 ib_idx = ib_markers%sf(i, j, k)
2161 if (ib_idx /= 0) then
2162 ! get the vector pointing to the grid cell from the IB centroid
2163 if (num_dims == 3) then
2164 radial_vector = [x_cc(i), y_cc(j), z_cc(k)] - [patch_ib(ib_idx)%x_centroid, patch_ib(ib_idx)%y_centroid, patch_ib(ib_idx)%z_centroid]
2165 else
2166 radial_vector = [x_cc(i), y_cc(j), 0._wp] - [patch_ib(ib_idx)%x_centroid, patch_ib(ib_idx)%y_centroid, 0._wp]
2167 end if
2168 dx = x_cc(i + 1) - x_cc(i)
2169 dy = y_cc(j + 1) - y_cc(j)
2170
2171 local_force_contribution(:) = 0._wp
2172 do fluid_idx = 0, num_fluids - 1
2173 ! Get the pressure contribution to force via a finite difference to compute the 2D components of the gradient of the pressure and cell volume
2174 local_force_contribution(1) = local_force_contribution(1) - (q_prim_vf(e_idx + fluid_idx)%sf(i + 1, j, k) - q_prim_vf(e_idx + fluid_idx)%sf(i - 1, j, k))/(2._wp*dx) ! force is the negative pressure gradient
2175 local_force_contribution(2) = local_force_contribution(2) - (q_prim_vf(e_idx + fluid_idx)%sf(i, j + 1, k) - q_prim_vf(e_idx + fluid_idx)%sf(i, j - 1, k))/(2._wp*dy)
2176 cell_volume = abs(dx*dy)
2177 ! add the 3D component of the pressure gradient, if we are working in 3 dimensions
2178 if (num_dims == 3) then
2179 dz = z_cc(k + 1) - z_cc(k)
2180 local_force_contribution(3) = local_force_contribution(3) - (q_prim_vf(e_idx + fluid_idx)%sf(i, j, k + 1) - q_prim_vf(e_idx + fluid_idx)%sf(i, j, k - 1))/(2._wp*dz)
2181 cell_volume = abs(cell_volume*dz)
2182 end if
2183 end do
2184
2185 ! Update the force values atomically to prevent race conditions
2186 call s_cross_product(radial_vector, local_force_contribution, local_torque_contribution)
2187
2188 ! get the viscous stress and add its contribution if that is considered
2189 ! TODO :: This is really bad code
2190 if (viscous) then
2191 ! compute the volume-weighted local dynamic viscosity
2192 dynamic_viscosity = 0._wp
2193 do fluid_idx = 1, num_fluids
2194 ! local dynamic viscosity is the dynamic viscosity of the fluid times alpha of the fluid
2195 dynamic_viscosity = dynamic_viscosity + (q_prim_vf(fluid_idx + advxb - 1)%sf(i, j, k)*dynamic_viscosities(fluid_idx))
2196 end do
2197
2198 ! get the linear force component first
2199 call s_compute_viscous_stress_tensor(viscous_stress_div_1, q_prim_vf, dynamic_viscosity, i - 1, j, k)
2200 call s_compute_viscous_stress_tensor(viscous_stress_div_2, q_prim_vf, dynamic_viscosity, i + 1, j, k)
2201 viscous_stress_div = (viscous_stress_div_2 - viscous_stress_div_1)/(2._wp*dx) ! get the x derivative of the viscous stress tensor
2202 local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(1, 1:3) ! add te x components of the derivative to the force
2203 do l = 1, 3
2204 ! take the cross products for the torque component
2205 call s_cross_product(radial_vector, viscous_stress_div_1(l, 1:3), viscous_cross_1(l, 1:3))
2206 call s_cross_product(radial_vector, viscous_stress_div_2(l, 1:3), viscous_cross_2(l, 1:3))
2207 end do
2208
2209 viscous_stress_div = (viscous_cross_2 - viscous_cross_1)/(2._wp*dx) ! get the x derivative of the cross product
2210 local_torque_contribution(1:3) = local_torque_contribution(1:3) + viscous_stress_div(1, 1:3) ! apply the cross product derivative to the torque
2211
2212 call s_compute_viscous_stress_tensor(viscous_stress_div_1, q_prim_vf, dynamic_viscosity, i, j - 1, k)
2213 call s_compute_viscous_stress_tensor(viscous_stress_div_2, q_prim_vf, dynamic_viscosity, i, j + 1, k)
2214 viscous_stress_div = (viscous_stress_div_2 - viscous_stress_div_1)/(2._wp*dy)
2215 local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(2, 1:3)
2216 do l = 1, 3
2217 call s_cross_product(radial_vector, viscous_stress_div_1(l, 1:3), viscous_cross_1(l, 1:3))
2218 call s_cross_product(radial_vector, viscous_stress_div_2(l, 1:3), viscous_cross_2(l, 1:3))
2219 end do
2220
2221 viscous_stress_div = (viscous_cross_2 - viscous_cross_1)/(2._wp*dy)
2222 local_torque_contribution(1:3) = local_torque_contribution(1:3) + viscous_stress_div(2, 1:3)
2223
2224 if (num_dims == 3) then
2225 call s_compute_viscous_stress_tensor(viscous_stress_div_1, q_prim_vf, dynamic_viscosity, i, j, k - 1)
2226 call s_compute_viscous_stress_tensor(viscous_stress_div_2, q_prim_vf, dynamic_viscosity, i, j, k + 1)
2227 viscous_stress_div = (viscous_stress_div_2 - viscous_stress_div_1)/(2._wp*dz)
2228 local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(3, 1:3)
2229 do l = 1, 3
2230 call s_cross_product(radial_vector, viscous_stress_div_1(l, 1:3), viscous_cross_1(l, 1:3))
2231 call s_cross_product(radial_vector, viscous_stress_div_2(l, 1:3), viscous_cross_2(l, 1:3))
2232 end do
2233 viscous_stress_div = (viscous_cross_2 - viscous_cross_1)/(2._wp*dz)
2234 local_torque_contribution(1:3) = local_torque_contribution(1:3) + viscous_stress_div(3, 1:3)
2235 end if
2236 end if
2237
2238 do l = 1, 3
2239
2240# 1157 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2241#if defined(MFC_OpenACC)
2242# 1157 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2243!$acc atomic update
2244# 1157 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2245#elif defined(MFC_OpenMP)
2246# 1157 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2247!$omp atomic update
2248# 1157 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2249#endif
2250 forces(ib_idx, l) = forces(ib_idx, l) + (local_force_contribution(l)*cell_volume)
2251
2252# 1159 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2253#if defined(MFC_OpenACC)
2254# 1159 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2255!$acc atomic update
2256# 1159 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2257#elif defined(MFC_OpenMP)
2258# 1159 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2259!$omp atomic update
2260# 1159 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2261#endif
2262 torques(ib_idx, l) = torques(ib_idx, l) + local_torque_contribution(l)*cell_volume
2263 end do
2264 end if
2265 end do
2266 end do
2267 end do
2268
2269# 1166 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2270
2271# 1166 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2272#if defined(MFC_OpenACC)
2273# 1166 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2274!$acc end parallel loop
2275# 1166 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2276#elif defined(MFC_OpenMP)
2277# 1166 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2278
2279# 1166 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2280
2281# 1166 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2282!$omp end target teams loop
2283# 1166 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2284#endif
2285# 1166 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2286
2287
2288 ! reduce the forces across all MPI ranks
2289 call s_mpi_allreduce_vectors_sum(forces, forces, num_ibs, 3)
2290 call s_mpi_allreduce_vectors_sum(torques, torques, num_ibs, 3)
2291
2292 ! consider body forces after reducing to avoid double counting
2293 do i = 1, num_ibs
2294 if (bf_x) then
2295 forces(i, 1) = forces(i, 1) + accel_bf(1)*patch_ib(i)%mass
2296 end if
2297 if (bf_y) then
2298 forces(i, 2) = forces(i, 2) + accel_bf(2)*patch_ib(i)%mass
2299 end if
2300 if (bf_z) then
2301 forces(i, 3) = forces(i, 3) + accel_bf(3)*patch_ib(i)%mass
2302 end if
2303 end do
2304
2305 ! apply the summed forces
2306 do i = 1, num_ibs
2307 patch_ib(i)%force(:) = forces(i, :)
2308 patch_ib(i)%torque(:) = matmul(patch_ib(i)%rotation_matrix_inverse, torques(i, :)) ! torques must be converted to the local coordinates of the IB
2309 end do
2310
2311 end subroutine s_compute_ib_forces
2312
2313 !> Subroutine to deallocate memory reserved for the IBM module
2314 impure subroutine s_finalize_ibm_module()
2315
2316#ifdef MFC_DEBUG
2317# 1196 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2318 block
2319# 1196 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2320 use iso_fortran_env, only: output_unit
2321# 1196 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2322
2323# 1196 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2324 print *, 'm_ibm.fpp:1196: ', '@:DEALLOCATE(ib_markers%sf)'
2325# 1196 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2326
2327# 1196 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2328 call flush (output_unit)
2329# 1196 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2330 end block
2331# 1196 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2332#endif
2333# 1196 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2334
2335# 1196 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2336#if defined(MFC_OpenACC)
2337# 1196 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2338!$acc exit data delete(ib_markers%sf)
2339# 1196 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2340#elif defined(MFC_OpenMP)
2341# 1196 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2342!$omp target exit data map(release:ib_markers%sf)
2343# 1196 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2344#endif
2345# 1196 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2346 deallocate (ib_markers%sf)
2347 if (allocated(airfoil_grid_u)) then
2348#ifdef MFC_DEBUG
2349# 1198 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2350 block
2351# 1198 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2352 use iso_fortran_env, only: output_unit
2353# 1198 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2354
2355# 1198 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2356 print *, 'm_ibm.fpp:1198: ', '@:DEALLOCATE(airfoil_grid_u)'
2357# 1198 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2358
2359# 1198 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2360 call flush (output_unit)
2361# 1198 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2362 end block
2363# 1198 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2364#endif
2365# 1198 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2366
2367# 1198 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2368#if defined(MFC_OpenACC)
2369# 1198 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2370!$acc exit data delete(airfoil_grid_u)
2371# 1198 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2372#elif defined(MFC_OpenMP)
2373# 1198 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2374!$omp target exit data map(release:airfoil_grid_u)
2375# 1198 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2376#endif
2377# 1198 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2378 deallocate (airfoil_grid_u)
2379#ifdef MFC_DEBUG
2380# 1199 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2381 block
2382# 1199 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2383 use iso_fortran_env, only: output_unit
2384# 1199 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2385
2386# 1199 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2387 print *, 'm_ibm.fpp:1199: ', '@:DEALLOCATE(airfoil_grid_l)'
2388# 1199 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2389
2390# 1199 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2391 call flush (output_unit)
2392# 1199 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2393 end block
2394# 1199 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2395#endif
2396# 1199 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2397
2398# 1199 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2399#if defined(MFC_OpenACC)
2400# 1199 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2401!$acc exit data delete(airfoil_grid_l)
2402# 1199 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2403#elif defined(MFC_OpenMP)
2404# 1199 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2405!$omp target exit data map(release:airfoil_grid_l)
2406# 1199 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2407#endif
2408# 1199 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2409 deallocate (airfoil_grid_l)
2410 end if
2411
2412 end subroutine s_finalize_ibm_module
2413
2414 !> Computes the center of mass for IB patch types where we are unable to determine their center of mass analytically.
2415 !> These patches include things like NACA airfoils and STL models
2416 subroutine s_compute_centroid_offset(ib_marker)
2417
2418 integer, intent(in) :: ib_marker
2419
2420 integer :: i, j, k, num_cells, num_cells_local
2421 real(wp), dimension(1:3) :: center_of_mass, center_of_mass_local
2422
2423 ! Offset only needs to be computes for specific geometries
2424 if (patch_ib(ib_marker)%geometry == 4 .or. &
2425 patch_ib(ib_marker)%geometry == 5 .or. &
2426 patch_ib(ib_marker)%geometry == 11 .or. &
2427 patch_ib(ib_marker)%geometry == 12) then
2428
2429 center_of_mass_local = [0._wp, 0._wp, 0._wp]
2430 num_cells_local = 0
2431
2432 ! get the summed mass distribution and number of cells to divide by
2433 do i = 0, m
2434 do j = 0, n
2435 do k = 0, p
2436 if (ib_markers%sf(i, j, k) == ib_marker) then
2437 num_cells_local = num_cells_local + 1
2438 center_of_mass_local = center_of_mass_local + [x_cc(i), y_cc(j), 0._wp]
2439 if (num_dims == 3) center_of_mass_local(3) = center_of_mass_local(3) + z_cc(k)
2440 end if
2441 end do
2442 end do
2443 end do
2444
2445 ! reduce the mass contribution over all MPI ranks and compute COM
2446 call s_mpi_allreduce_integer_sum(num_cells_local, num_cells)
2447 if (num_cells /= 0) then
2448 call s_mpi_allreduce_sum(center_of_mass_local(1), center_of_mass(1))
2449 call s_mpi_allreduce_sum(center_of_mass_local(2), center_of_mass(2))
2450 call s_mpi_allreduce_sum(center_of_mass_local(3), center_of_mass(3))
2451 center_of_mass = center_of_mass/real(num_cells, wp)
2452 else
2453 patch_ib(ib_marker)%centroid_offset = [0._wp, 0._wp, 0._wp]
2454 return
2455 end if
2456
2457 ! assign the centroid offset as a vector pointing from the true COM to the "centroid" in the input file and replace the current centroid
2458 patch_ib(ib_marker)%centroid_offset = [patch_ib(ib_marker)%x_centroid, patch_ib(ib_marker)%y_centroid, patch_ib(ib_marker)%z_centroid] &
2459 - center_of_mass
2460 patch_ib(ib_marker)%x_centroid = center_of_mass(1)
2461 patch_ib(ib_marker)%y_centroid = center_of_mass(2)
2462 patch_ib(ib_marker)%z_centroid = center_of_mass(3)
2463
2464 ! rotate the centroid offset back into the local coords of the IB
2465 patch_ib(ib_marker)%centroid_offset = matmul(patch_ib(ib_marker)%rotation_matrix_inverse, patch_ib(ib_marker)%centroid_offset)
2466 else
2467 patch_ib(ib_marker)%centroid_offset(:) = [0._wp, 0._wp, 0._wp]
2468 end if
2469
2470 end subroutine s_compute_centroid_offset
2471
2472 !> Computes the moment of inertia for an immersed boundary
2473 !! @param ib_marker Immersed boundary marker index
2474 subroutine s_compute_moment_of_inertia(ib_marker, axis)
2475
2476 real(wp), dimension(3), intent(in) :: axis !< the axis about which we compute the moment. Only required in 3D.
2477 integer, intent(in) :: ib_marker
2478
2479 real(wp) :: moment, distance_to_axis, cell_volume
2480 real(wp), dimension(3) :: position, closest_point_along_axis, vector_to_axis, normal_axis
2481 integer :: i, j, k, count
2482
2483 if (p == 0) then
2484 normal_axis = [0, 0, 1]
2485 else if (sqrt(sum(axis**2)) == 0) then
2486 ! if the object is not actually rotating at this time, return a dummy value and exit
2487 patch_ib(ib_marker)%moment = 1._wp
2488 return
2489 else
2490 normal_axis = axis/sqrt(sum(axis))
2491 end if
2492
2493 ! if the IB is in 2D or a 3D sphere, we can compute this exactly
2494 if (patch_ib(ib_marker)%geometry == 2) then ! circle
2495 patch_ib(ib_marker)%moment = 0.5_wp*patch_ib(ib_marker)%mass*(patch_ib(ib_marker)%radius)**2
2496 elseif (patch_ib(ib_marker)%geometry == 3) then ! rectangle
2497 patch_ib(ib_marker)%moment = patch_ib(ib_marker)%mass*(patch_ib(ib_marker)%length_x**2 + patch_ib(ib_marker)%length_y**2)/6._wp
2498 elseif (patch_ib(ib_marker)%geometry == 6) then ! ellipse
2499 patch_ib(ib_marker)%moment = 0.0625_wp*patch_ib(ib_marker)%mass*(patch_ib(ib_marker)%length_x**2 + patch_ib(ib_marker)%length_y**2)
2500 elseif (patch_ib(ib_marker)%geometry == 8) then ! sphere
2501 patch_ib(ib_marker)%moment = 0.4*patch_ib(ib_marker)%mass*(patch_ib(ib_marker)%radius)**2
2502
2503 else ! we do not have an analytic moment of inertia calculation and need to approximate it directly via a sum
2504 count = 0
2505 moment = 0._wp
2506 cell_volume = (x_cc(1) - x_cc(0))*(y_cc(1) - y_cc(0)) ! computed without grid stretching. Update in the loop to perform with stretching
2507 if (p /= 0) then
2508 cell_volume = cell_volume*(z_cc(1) - z_cc(0))
2509 end if
2510
2511
2512# 1301 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2513
2514# 1301 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2515#if defined(MFC_OpenACC)
2516# 1301 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2517!$acc parallel loop collapse(3) gang vector default(present) private(position, closest_point_along_axis, vector_to_axis, distance_to_axis) copy(moment, count) copyin(ib_marker, cell_volume, normal_axis)
2518# 1301 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2519#elif defined(MFC_OpenMP)
2520# 1301 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2521
2522# 1301 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2523
2524# 1301 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2525
2526# 1301 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2527!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(position, closest_point_along_axis, vector_to_axis, distance_to_axis) map(tofrom:moment, count) map(to:ib_marker, cell_volume, normal_axis)
2528# 1301 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2529#endif
2530# 1301 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2531
2532 do i = 0, m
2533 do j = 0, n
2534 do k = 0, p
2535 if (ib_markers%sf(i, j, k) == ib_marker) then
2536
2537# 1306 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2538#if defined(MFC_OpenACC)
2539# 1306 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2540!$acc atomic update
2541# 1306 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2542#elif defined(MFC_OpenMP)
2543# 1306 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2544!$omp atomic update
2545# 1306 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2546#endif
2547 count = count + 1 ! increment the count of total cells in the boundary
2548
2549 ! get the position in local coordinates so that the axis passes through 0, 0, 0
2550 if (p == 0) then
2551 position = [x_cc(i), y_cc(j), 0._wp] - [patch_ib(ib_marker)%x_centroid, patch_ib(ib_marker)%y_centroid, 0._wp]
2552 else
2553 position = [x_cc(i), y_cc(j), z_cc(k)] - [patch_ib(ib_marker)%x_centroid, patch_ib(ib_marker)%y_centroid, patch_ib(ib_marker)%z_centroid]
2554 end if
2555
2556 ! project the position along the axis to find the closest distance to the rotation axis
2557 closest_point_along_axis = normal_axis*dot_product(normal_axis, position)
2558 vector_to_axis = position - closest_point_along_axis
2559 distance_to_axis = dot_product(vector_to_axis, vector_to_axis) ! saves the distance to the axis squared
2560
2561 ! compute the position component of the moment
2562
2563# 1322 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2564#if defined(MFC_OpenACC)
2565# 1322 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2566!$acc atomic update
2567# 1322 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2568#elif defined(MFC_OpenMP)
2569# 1322 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2570!$omp atomic update
2571# 1322 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2572#endif
2573 moment = moment + distance_to_axis
2574 end if
2575 end do
2576 end do
2577 end do
2578
2579# 1328 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2580
2581# 1328 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2582#if defined(MFC_OpenACC)
2583# 1328 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2584!$acc end parallel loop
2585# 1328 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2586#elif defined(MFC_OpenMP)
2587# 1328 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2588
2589# 1328 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2590
2591# 1328 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2592!$omp end target teams loop
2593# 1328 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2594#endif
2595# 1328 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2596
2597
2598 ! write the final moment assuming the points are all uniform density
2599 patch_ib(ib_marker)%moment = moment*patch_ib(ib_marker)%mass/(count*cell_volume)
2600
2601# 1332 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2602#if defined(MFC_OpenACC)
2603# 1332 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2604!$acc update device(patch_ib(ib_marker)%moment)
2605# 1332 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2606#elif defined(MFC_OpenMP)
2607# 1332 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2608!$omp target update to(patch_ib(ib_marker)%moment)
2609# 1332 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2610#endif
2611 end if
2612
2613 end subroutine s_compute_moment_of_inertia
2614
2615 !> @brief Computes the cross product c = a x b of two 3D vectors.
2616 subroutine s_cross_product(a, b, c)
2617
2618# 1339 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2619#if MFC_OpenACC
2620# 1339 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2621!$acc routine seq
2622# 1339 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2623#elif MFC_OpenMP
2624# 1339 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2625
2626# 1339 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2627
2628# 1339 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2629!$omp declare target device_type(any)
2630# 1339 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2631#endif
2632 real(wp), intent(in) :: a(3), b(3)
2633 real(wp), intent(out) :: c(3)
2634
2635 c(1) = a(2)*b(3) - a(3)*b(2)
2636 c(2) = a(3)*b(1) - a(1)*b(3)
2637 c(3) = a(1)*b(2) - a(2)*b(1)
2638 end subroutine s_cross_product
2639
2640end module m_ibm
q_cons_vf
type(scalar_field), dimension(sys_size), intent(inout) q_cons_vf
Definition m_phase_change.fpp.f90:1119
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_compute_levelset
Computes signed-distance level-set fields and surface normals for immersed-boundary patch geometries.
Definition m_compute_levelset.fpp.f90:292
m_constants
Compile-time constant parameters: default values, tolerances, and physical constants.
Definition m_constants.fpp.f90:7
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::buff_size
integer buff_size
The number of cells that are necessary to be able to store enough boundary conditions data to march t...
Definition m_global_parameters.fpp.f90:906
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::m
integer m
Definition m_global_parameters.fpp.f90:331
m_helper_basic
Basic floating-point utilities: approximate equality, default detection, and coordinate bounds.
Definition m_helper_basic.fpp.f90:292
m_helper
Utility routines for bubble model setup, coordinate transforms, array sampling, and special functions...
Definition m_helper.fpp.f90:303
m_ib_patches
Allocate memory and read initial condition data for IC extrusion.
Definition m_ib_patches.fpp.f90:362
m_ibm
Ghost-node immersed boundary method: locates ghost/image points, computes interpolation coefficients,...
Definition m_ibm.fpp.f90:292
m_ibm::s_update_mib
impure subroutine s_update_mib(num_ibs)
Resets the current indexes of immersed boundaries and replaces them after updating the position of ea...
Definition m_ibm.fpp.f90:1985
m_ibm::s_find_ghost_points
subroutine, private s_find_ghost_points(ghost_points_in, inner_points_in)
Function that finds the ghost points.
Definition m_ibm.fpp.f90:1495
m_ibm::s_compute_ib_forces
subroutine s_compute_ib_forces(q_prim_vf, fluid_pp)
Computes pressure and viscous forces and torques on immersed bodies via a volume integration method.
Definition m_ibm.fpp.f90:2108
m_ibm::s_ibm_setup
impure subroutine, public s_ibm_setup()
Initializes the values of various IBM variables, such as ghost points and image points.
Definition m_ibm.fpp.f90:568
m_ibm::s_compute_interpolation_coeffs
subroutine, private s_compute_interpolation_coeffs(ghost_points_in)
Function that computes the interpolation coefficients of image points.
Definition m_ibm.fpp.f90:1613
m_ibm::s_finalize_ibm_module
impure subroutine, public s_finalize_ibm_module()
Subroutine to deallocate memory reserved for the IBM module.
Definition m_ibm.fpp.f90:2315
m_ibm::s_interpolate_image_point
subroutine, private s_interpolate_image_point(q_prim_vf, gp, alpha_rho_ip, alpha_ip, pres_ip, vel_ip, c_ip, r_ip, v_ip, pb_ip, mv_ip, nmom_ip, pb_in, mv_in, presb_ip, massv_ip)
Function that uses the interpolation coefficients and the current state at the cell centers in order ...
Definition m_ibm.fpp.f90:1785
m_ibm::s_populate_ib_buffers
subroutine s_populate_ib_buffers()
Exchanges immersed boundary marker data across MPI subdomain boundaries.
Definition m_ibm.fpp.f90:817
m_ibm::ghost_points
type(ghost_point), dimension(:), allocatable ghost_points
Definition m_ibm.fpp.f90:339
m_ibm::num_gps
integer num_gps
Number of ghost points.
Definition m_ibm.fpp.f90:353
m_ibm::inner_points
type(ghost_point), dimension(:), allocatable inner_points
Definition m_ibm.fpp.f90:340
m_ibm::moving_immersed_boundary_flag
logical moving_immersed_boundary_flag
Definition m_ibm.fpp.f90:380
m_ibm::s_compute_moment_of_inertia
subroutine s_compute_moment_of_inertia(ib_marker, axis)
Computes the moment of inertia for an immersed boundary.
Definition m_ibm.fpp.f90:2475
m_ibm::s_cross_product
subroutine s_cross_product(a, b, c)
Computes the cross product c = a x b of two 3D vectors.
Definition m_ibm.fpp.f90:2617
m_ibm::s_compute_centroid_offset
subroutine s_compute_centroid_offset(ib_marker)
Computes the center of mass for IB patch types where we are unable to determine their center of mass ...
Definition m_ibm.fpp.f90:2417
m_ibm::ib_markers
type(integer_field), public ib_markers
Definition m_ibm.fpp.f90:326
m_ibm::s_initialize_ibm_module
impure subroutine, public s_initialize_ibm_module()
Allocates memory for the variables in the IBM module.
Definition m_ibm.fpp.f90:386
m_ibm::s_ibm_correct_state
subroutine, public s_ibm_correct_state(q_cons_vf, q_prim_vf, pb_in, mv_in)
Subroutine that updates the conservative variables at the ghost points.
Definition m_ibm.fpp.f90:856
m_ibm::s_find_num_ghost_points
subroutine, private s_find_num_ghost_points(num_gps_out, num_inner_gps_out)
Function that finds the number of ghost points, used for allocating memory.
Definition m_ibm.fpp.f90:1447
m_ibm::s_compute_image_points
impure subroutine, private s_compute_image_points(ghost_points_in)
Function that computes the image points for each ghost point.
Definition m_ibm.fpp.f90:1354
m_ibm::num_inner_gps
integer num_inner_gps
Number of ghost points.
Definition m_ibm.fpp.f90:354
m_mpi_proxy
MPI halo exchange, domain decomposition, and buffer packing/unpacking for the simulation solver.
Definition m_mpi_proxy.fpp.f90:303
m_variables_conversion
Conservative-to-primitive variable conversion, mixture property evaluation, and pressure computation.
Definition m_variables_conversion.fpp.f90:303
m_viscous
Computes viscous stress tensors and diffusive flux contributions for the Navier–Stokes equations.
Definition m_viscous.fpp.f90:302
m_derived_types::ghost_point
Ghost Point for Immersed Boundaries.
Definition m_derived_types.fpp.f90:727
m_derived_types::integer_field
Derived type annexing an integer scalar field (SF).
Definition m_derived_types.fpp.f90:317