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