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