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 use m_collisions
332
333 implicit none
334
335 private :: s_compute_image_points, s_compute_interpolation_coeffs, s_interpolate_image_point, s_find_ghost_points, &
336 & s_find_num_ghost_points
337 ; public :: ib_gbl_idx_lookup, s_initialize_ibm_module, s_ibm_setup, s_ibm_correct_state, s_finalize_ibm_module
338
339 type(integer_field), public :: ib_markers
340
341# 31 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
342#if defined(MFC_OpenACC)
343# 31 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
344!$acc declare create(ib_markers)
345# 31 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
346#elif defined(MFC_OpenMP)
347# 31 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
348!$omp declare target (ib_markers)
349# 31 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
350#endif
351
352 type(ghost_point), dimension(:), allocatable :: ghost_points
353
354# 34 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
355#if defined(MFC_OpenACC)
356# 34 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
357!$acc declare create(ghost_points)
358# 34 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
359#elif defined(MFC_OpenMP)
360# 34 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
361!$omp declare target (ghost_points)
362# 34 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
363#endif
364
365 integer, dimension(:), allocatable :: ib_gbl_idx_lookup
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(ib_gbl_idx_lookup)
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 (ib_gbl_idx_lookup)
375# 37 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
376#endif
377
378 integer :: num_gps !< Number of ghost points
379#if defined(MFC_OpenACC)
380
381# 41 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
382#if defined(MFC_OpenACC)
383# 41 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
384!$acc declare create(gp_layers, num_gps)
385# 41 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
386#elif defined(MFC_OpenMP)
387# 41 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
388!$omp declare target (gp_layers, num_gps)
389# 41 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
390#endif
391#elif defined(MFC_OpenMP)
392
393# 43 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
394#if defined(MFC_OpenACC)
395# 43 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
396!$acc declare create(num_gps)
397# 43 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
398#elif defined(MFC_OpenMP)
399# 43 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
400!$omp declare target (num_gps)
401# 43 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
402#endif
403#endif
404 logical :: moving_immersed_boundary_flag
405
406 ! IB MPI buffers
407 integer, allocatable :: send_ids(:), recv_ids(:)
408 real(wp), allocatable :: send_ft(:,:), recv_ft(:,:)
409 real(wp), allocatable :: recv_forces_snap(:,:), recv_torques_snap(:,:)
410
411contains
412
413 !> Allocates memory for the variables in the IBM module
414 impure subroutine s_initialize_ibm_module()
415
416 if (p > 0) then
417#ifdef MFC_DEBUG
418# 58 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
419 block
420# 58 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
421 use iso_fortran_env, only: output_unit
422# 58 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
423
424# 58 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
425 print *, 'm_ibm.fpp:58: ', '@:ALLOCATE(ib_markers%sf(-buff_size:m+buff_size, -buff_size:n+buff_size, -buff_size:p+buff_size))'
426# 58 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
427
428# 58 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
429 call flush (output_unit)
430# 58 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
431 end block
432# 58 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
433#endif
434# 58 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
435 allocate (ib_markers%sf(-buff_size:m+buff_size, -buff_size:n+buff_size, -buff_size:p+buff_size))
436# 58 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
437
438# 58 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
439
440# 58 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
441#if defined(MFC_OpenACC)
442# 58 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
443!$acc enter data create(ib_markers%sf)
444# 58 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
445#elif defined(MFC_OpenMP)
446# 58 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
447!$omp target enter data map(always,alloc:ib_markers%sf)
448# 58 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
449#endif
450 else
451#ifdef MFC_DEBUG
452# 60 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
453 block
454# 60 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
455 use iso_fortran_env, only: output_unit
456# 60 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
457
458# 60 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
459 print *, 'm_ibm.fpp:60: ', '@:ALLOCATE(ib_markers%sf(-buff_size:m+buff_size, -buff_size:n+buff_size, 0:0))'
460# 60 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
461
462# 60 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
463 call flush (output_unit)
464# 60 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
465 end block
466# 60 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
467#endif
468# 60 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
469 allocate (ib_markers%sf(-buff_size:m+buff_size, -buff_size:n+buff_size, 0:0))
470# 60 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
471
472# 60 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
473
474# 60 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
475#if defined(MFC_OpenACC)
476# 60 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
477!$acc enter data create(ib_markers%sf)
478# 60 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
479#elif defined(MFC_OpenMP)
480# 60 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
481!$omp target enter data map(always,alloc:ib_markers%sf)
482# 60 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
483#endif
484 end if
485
486#ifdef _CRAYFTN
487# 63 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
488 block
489# 63 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
490#ifdef MFC_DEBUG
491# 63 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
492 block
493# 63 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
494 use iso_fortran_env, only: output_unit
495# 63 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
496
497# 63 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
498 print *, 'm_ibm.fpp:63: ', '@:ACC_SETUP_SFs(ib_markers)'
499# 63 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
500
501# 63 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
502 call flush (output_unit)
503# 63 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
504 end block
505# 63 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
506#endif
507# 63 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
508
509# 63 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
510
511# 63 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
512#if defined(MFC_OpenACC)
513# 63 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
514!$acc enter data copyin(ib_markers)
515# 63 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
516#elif defined(MFC_OpenMP)
517# 63 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
518!$omp target enter data map(to:ib_markers)
519# 63 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
520#endif
521# 63 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
522 if (associated(ib_markers%sf)) then
523# 63 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
524
525# 63 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
526#if defined(MFC_OpenACC)
527# 63 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
528!$acc enter data copyin(ib_markers%sf)
529# 63 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
530#elif defined(MFC_OpenMP)
531# 63 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
532!$omp target enter data map(to:ib_markers%sf)
533# 63 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
534#endif
535# 63 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
536 end if
537# 63 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
538 end block
539# 63 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
540#endif
541
542
543# 65 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
544#if defined(MFC_OpenACC)
545# 65 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
546!$acc enter data copyin(num_gps)
547# 65 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
548#elif defined(MFC_OpenMP)
549# 65 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
550!$omp target enter data map(to:num_gps)
551# 65 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
552#endif
553
554 if (collision_model > 0) call s_initialize_collisions_module()
555
556 end subroutine s_initialize_ibm_module
557
558 !> Initializes the values of various IBM variables, such as ghost points and image points.
559 impure subroutine s_ibm_setup()
560
561 integer :: i, j, k
562 integer :: max_num_gps
563
564 call nvtxstartrange("SETUP-IBM-MODULE")
565
566 ! GPU routines require updated cell centers
567
568# 80 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
569#if defined(MFC_OpenACC)
570# 80 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
571!$acc update device(num_ibs, num_gbl_ibs, x_cc, y_cc, dx, dy, x_domain, y_domain, ib_bc_x%beg, ib_bc_y%beg)
572# 80 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
573#elif defined(MFC_OpenMP)
574# 80 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
575!$omp target update to(num_ibs, num_gbl_ibs, x_cc, y_cc, dx, dy, x_domain, y_domain, ib_bc_x%beg, ib_bc_y%beg)
576# 80 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
577#endif
578 if (p /= 0) then
579
580# 82 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
581#if defined(MFC_OpenACC)
582# 82 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
583!$acc update device(z_cc, dz, z_domain, ib_bc_z%beg)
584# 82 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
585#elif defined(MFC_OpenMP)
586# 82 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
587!$omp target update to(z_cc, dz, z_domain, ib_bc_z%beg)
588# 82 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
589#endif
590 end if
591
592# 84 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
593#if defined(MFC_OpenACC)
594# 84 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
595!$acc update device(patch_ib(1:num_ibs))
596# 84 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
597#elif defined(MFC_OpenMP)
598# 84 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
599!$omp target update to(patch_ib(1:num_ibs))
600# 84 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
601#endif
602
603 ! do all set up for moving immersed boundaries
604
605# 87 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
606
607# 87 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
608#if defined(MFC_OpenACC)
609# 87 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
610!$acc parallel loop gang vector default(present) private(i)
611# 87 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
612#elif defined(MFC_OpenMP)
613# 87 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
614
615# 87 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
616
617# 87 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
618
619# 87 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
620!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i)
621# 87 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
622#endif
623 do i = 1, num_ibs
624 if (patch_ib(i)%moving_ibm /= 0) then
625 call s_compute_moment_of_inertia(i, patch_ib(i)%angular_vel)
626 end if
627 call s_update_ib_rotation_matrix(i)
628 end do
629
630# 94 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
631#if defined(MFC_OpenACC)
632# 94 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
633!$acc end parallel loop
634# 94 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
635#elif defined(MFC_OpenMP)
636# 94 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
637
638# 94 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
639!$omp end target teams loop
640# 94 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
641#endif
642
643# 95 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
644#if defined(MFC_OpenACC)
645# 95 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
646!$acc update host(patch_ib(1:num_ibs))
647# 95 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
648#elif defined(MFC_OpenMP)
649# 95 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
650!$omp target update from(patch_ib(1:num_ibs))
651# 95 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
652#endif
653
654 ! allocate some arrays for MPI communication, if required by this simulation
655#ifdef MFC_MPI
656 if (num_procs > 1) then
657#ifdef MFC_DEBUG
658# 100 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
659 block
660# 100 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
661 use iso_fortran_env, only: output_unit
662# 100 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
663
664# 100 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
665 print *, 'm_ibm.fpp:100: ', '@:ALLOCATE(send_ids(size(patch_ib)), send_ft(6, size(patch_ib)))'
666# 100 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
667
668# 100 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
669 call flush (output_unit)
670# 100 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
671 end block
672# 100 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
673#endif
674# 100 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
675 allocate (send_ids(size(patch_ib)), send_ft(6, size(patch_ib)))
676# 100 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
677
678# 100 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
679
680# 100 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
681
682# 100 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
683#if defined(MFC_OpenACC)
684# 100 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
685!$acc enter data create(send_ids, send_ft)
686# 100 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
687#elif defined(MFC_OpenMP)
688# 100 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
689!$omp target enter data map(always,alloc:send_ids, send_ft)
690# 100 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
691#endif
692 allocate (recv_forces_snap(size(patch_ib), 3), recv_torques_snap(size(patch_ib), 3), recv_ids(size(patch_ib)), &
693 & recv_ft(6, size(patch_ib)))
694 end if
695#endif
696
697 call s_update_ib_lookup()
698
699 ! recompute the new ib_patch locations
700 ib_markers%sf = 0._wp
701
702# 110 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
703#if defined(MFC_OpenACC)
704# 110 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
705!$acc update device(ib_markers%sf)
706# 110 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
707#elif defined(MFC_OpenMP)
708# 110 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
709!$omp target update to(ib_markers%sf)
710# 110 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
711#endif
712 call s_apply_ib_patches(ib_markers)
713
714# 112 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
715#if defined(MFC_OpenACC)
716# 112 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
717!$acc update host(ib_markers%sf)
718# 112 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
719#elif defined(MFC_OpenMP)
720# 112 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
721!$omp target update from(ib_markers%sf)
722# 112 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
723#endif
724 do i = 1, num_ibs
725 if (patch_ib(i)%moving_ibm /= 0) call s_compute_centroid_offset(i) ! offsets are computed after IB markers are generated
726
727# 115 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
728#if defined(MFC_OpenACC)
729# 115 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
730!$acc update device(patch_ib(i))
731# 115 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
732#elif defined(MFC_OpenMP)
733# 115 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
734!$omp target update to(patch_ib(i))
735# 115 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
736#endif
737 end do
738
739 ! find the number of ghost points and set them to be the maximum total across ranks
740 call s_find_num_ghost_points(num_gps)
741 if (moving_immersed_boundary_flag) then
742 call s_mpi_allreduce_integer_sum(num_gps, max_num_gps)
743 max_num_gps = min(max_num_gps*2, (m + 1)*(n + 1)*(p + 1))
744 else
745 max_num_gps = num_gps
746 end if
747
748 ! set the size of the ghost point arrays to be the amount of points total, plus a factor of 2 buffer
749
750# 128 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
751#if defined(MFC_OpenACC)
752# 128 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
753!$acc update device(num_gps)
754# 128 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
755#elif defined(MFC_OpenMP)
756# 128 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
757!$omp target update to(num_gps)
758# 128 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
759#endif
760#ifdef MFC_DEBUG
761# 129 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
762 block
763# 129 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
764 use iso_fortran_env, only: output_unit
765# 129 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
766
767# 129 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
768 print *, 'm_ibm.fpp:129: ', '@:ALLOCATE(ghost_points(1:max_num_gps))'
769# 129 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
770
771# 129 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
772 call flush (output_unit)
773# 129 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
774 end block
775# 129 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
776#endif
777# 129 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
778 allocate (ghost_points(1:max_num_gps))
779# 129 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
780
781# 129 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
782
783# 129 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
784#if defined(MFC_OpenACC)
785# 129 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
786!$acc enter data create(ghost_points)
787# 129 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
788#elif defined(MFC_OpenMP)
789# 129 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
790!$omp target enter data map(always,alloc:ghost_points)
791# 129 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
792#endif
793
794
795# 131 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
796#if defined(MFC_OpenACC)
797# 131 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
798!$acc enter data copyin(ghost_points)
799# 131 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
800#elif defined(MFC_OpenMP)
801# 131 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
802!$omp target enter data map(to:ghost_points)
803# 131 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
804#endif
805 ! Ghost-cell IBM, Tseng & Ferziger JCP (2003), Mittal & Iaccarino ARFM (2005)
806 call s_find_ghost_points(ghost_points)
807 call s_apply_levelset(ghost_points, num_gps)
808
809 call s_compute_image_points(ghost_points)
810 call s_compute_interpolation_coeffs(ghost_points)
811
812 call nvtxendrange
813
814 end subroutine s_ibm_setup
815
816 !> Update the conservative variables at the ghost points
817 subroutine s_ibm_correct_state(q_cons_vf, q_prim_vf, pb_in, mv_in)
818
819 type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf !< Primitive Variables
820 type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf !< Primitive Variables
821 real(stp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:,1:), optional, intent(inout) :: pb_in, mv_in
822 integer :: i, j, k, l, q, r !< Iterator variables
823 integer :: patch_id, patch_id_temp !< Patch ID of ghost point
824 real(wp) :: rho, gamma, pi_inf, dyn_pres !< Mixture variables
825 real(wp), dimension(2) :: re_k
826 real(wp) :: g_k
827 real(wp) :: qv_k
828 real(wp) :: pres_ip
829 real(wp), dimension(3) :: vel_ip, vel_norm_ip
830 real(wp) :: c_ip
831
832# 166 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
833 real(wp), dimension(num_fluids) :: gs
834 real(wp), dimension(num_fluids) :: alpha_rho_ip, alpha_ip
835 real(wp), dimension(nb) :: r_ip, v_ip, pb_ip, mv_ip
836 real(wp), dimension(nb*nmom) :: nmom_ip
837 real(wp), dimension(nb*nnode) :: presb_ip, massv_ip
838# 172 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
839 ! Primitive variables at the image point associated with a ghost point, interpolated from surrounding fluid cells.
840
841 real(wp), dimension(3) :: norm !< Normal vector from GP to IP
842 real(wp), dimension(3) :: physical_loc !< Physical loc of GP
843 real(wp), dimension(3) :: vel_g !< Velocity of GP
844 real(wp), dimension(3) :: radial_vector !< vector from centroid to ghost point
845 real(wp), dimension(3) :: rotation_velocity !< speed of the ghost point due to rotation
846 real(wp) :: nbub
847 real(wp) :: buf
848 type(ghost_point) :: gp
849 type(ghost_point) :: innerp
850
851 ! set the Moving IBM interior conservative variables
852
853# 185 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
854
855# 185 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
856#if defined(MFC_OpenACC)
857# 185 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
858!$acc parallel loop collapse(3) gang vector default(present) private(i, j, k, patch_id, rho)
859# 185 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
860#elif defined(MFC_OpenMP)
861# 185 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
862
863# 185 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
864
865# 185 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
866
867# 185 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
868!$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)
869# 185 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
870#endif
871 do l = 0, p
872 do k = 0, n
873 do j = 0, m
874 patch_id = ib_markers%sf(j, k, l)
875 if (patch_id /= 0) then
876 call s_decode_patch_periodicity(patch_id, patch_id_temp)
877 call s_get_neighborhood_idx(patch_id_temp, patch_id)
878 if (patch_id > 0) then
879 q_prim_vf(eqn_idx%E)%sf(j, k, l) = 1._wp
880 rho = 0._wp
881 do i = 1, num_fluids
882 rho = rho + q_prim_vf(eqn_idx%cont%beg + i - 1)%sf(j, k, l)
883 end do
884
885 ! Sets the momentum
886 do i = 1, num_dims
887 q_cons_vf(eqn_idx%mom%beg + i - 1)%sf(j, k, l) = patch_ib(patch_id)%vel(i)*rho
888 q_prim_vf(eqn_idx%mom%beg + i - 1)%sf(j, k, l) = patch_ib(patch_id)%vel(i)
889 end do
890 end if ! patch_id > 0
891 end if
892 end do
893 end do
894 end do
895
896# 210 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
897#if defined(MFC_OpenACC)
898# 210 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
899!$acc end parallel loop
900# 210 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
901#elif defined(MFC_OpenMP)
902# 210 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
903
904# 210 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
905!$omp end target teams loop
906# 210 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
907#endif
908
909 if (num_gps > 0) then
910
911# 213 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
912
913# 213 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
914#if defined(MFC_OpenACC)
915# 213 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
916!$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)
917# 213 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
918#elif defined(MFC_OpenMP)
919# 213 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
920
921# 213 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
922
923# 213 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
924
925# 213 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
926!$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)
927# 213 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
928#endif
929# 216 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
930 do i = 1, num_gps
931 gp = ghost_points(i)
932 j = gp%loc(1)
933 k = gp%loc(2)
934 l = gp%loc(3)
935 patch_id = ghost_points(i)%ib_patch_id
936
937 ! Calculate physical location of GP
938 if (p > 0) then
939 physical_loc = [x_cc(j), y_cc(k), z_cc(l)]
940 else
941 physical_loc = [x_cc(j), y_cc(k), 0._wp]
942 end if
943
944 ! Interpolate primitive variables at image point associated w/ GP
945 if (bubbles_euler .and. .not. qbmm) then
946 call s_interpolate_image_point(q_prim_vf, gp, alpha_rho_ip, alpha_ip, pres_ip, vel_ip, c_ip, r_ip, v_ip, &
947 & pb_ip, mv_ip)
948 else if (qbmm .and. polytropic) then
949 call s_interpolate_image_point(q_prim_vf, gp, alpha_rho_ip, alpha_ip, pres_ip, vel_ip, c_ip, r_ip, v_ip, &
950 & pb_ip, mv_ip, nmom_ip)
951 else if (qbmm .and. .not. polytropic) then
952 call s_interpolate_image_point(q_prim_vf, gp, alpha_rho_ip, alpha_ip, pres_ip, vel_ip, c_ip, r_ip, v_ip, &
953 & pb_ip, mv_ip, nmom_ip, pb_in, mv_in, presb_ip, massv_ip)
954 else
955 call s_interpolate_image_point(q_prim_vf, gp, alpha_rho_ip, alpha_ip, pres_ip, vel_ip, c_ip)
956 end if
957
958 dyn_pres = 0._wp
959
960 ! Set q_prim_vf params at GP so that mixture vars calculated properly
961
962# 247 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
963#if defined(MFC_OpenACC)
964# 247 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
965!$acc loop seq
966# 247 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
967#elif defined(MFC_OpenMP)
968# 247 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
969
970# 247 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
971#endif
972 do q = 1, num_fluids
973 q_prim_vf(q)%sf(j, k, l) = alpha_rho_ip(q)
974 q_prim_vf(eqn_idx%adv%beg + q - 1)%sf(j, k, l) = alpha_ip(q)
975 end do
976
977 if (surface_tension) then
978 q_prim_vf(eqn_idx%c)%sf(j, k, l) = c_ip
979 end if
980
981 ! set the pressure
982 if (patch_ib(patch_id)%moving_ibm <= 1) then
983 q_prim_vf(eqn_idx%E)%sf(j, k, l) = pres_ip
984 else
985 q_prim_vf(eqn_idx%E)%sf(j, k, l) = 0._wp
986
987# 262 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
988#if defined(MFC_OpenACC)
989# 262 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
990!$acc loop seq
991# 262 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
992#elif defined(MFC_OpenMP)
993# 262 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
994
995# 262 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
996#endif
997 do q = 1, num_fluids
998 ! Pressure correction for moving IB: accounts for acceleration of IB surface
999 q_prim_vf(eqn_idx%E)%sf(j, k, l) = q_prim_vf(eqn_idx%E)%sf(j, k, &
1000 & l) + pres_ip/(1._wp - 2._wp*abs(gp%levelset*alpha_rho_ip(q)/pres_ip) &
1001 & *dot_product(patch_ib(patch_id)%force/patch_ib(patch_id)%mass, gp%levelset_norm))
1002 end do
1003 end if
1004
1005 if (model_eqns /= 4) then
1006 ! If in simulation, use acc mixture subroutines
1007 if (elasticity) then
1008 call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv_k, alpha_ip, alpha_rho_ip, re_k, &
1009 & g_k, gs)
1010 else
1011 call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv_k, alpha_ip, alpha_rho_ip, re_k)
1012 end if
1013 end if
1014
1015 ! Calculate velocity of ghost cell
1016 if (gp%slip) then
1017 norm(1:3) = gp%levelset_norm
1018 buf = sqrt(sum(norm**2))
1019 norm = norm/buf
1020 vel_norm_ip = sum(vel_ip*norm)*norm
1021 vel_g = vel_ip - vel_norm_ip
1022 if (patch_ib(patch_id)%moving_ibm /= 0) then
1023 ! compute the linear velocity of the ghost point due to rotation
1024 radial_vector = physical_loc - [patch_ib(patch_id)%x_centroid, patch_ib(patch_id)%y_centroid, &
1025 & patch_ib(patch_id)%z_centroid]
1026 call s_cross_product(patch_ib(patch_id)%angular_vel, radial_vector, rotation_velocity)
1027
1028 ! add only the component of the IB's motion that is normal to the surface
1029 vel_g = vel_g + sum((patch_ib(patch_id)%vel + rotation_velocity)*norm)*norm
1030 end if
1031 else
1032 if (patch_ib(patch_id)%moving_ibm == 0) then
1033 ! we know the object is not moving if moving_ibm is 0 (false)
1034 vel_g = 0._wp
1035 else
1036 ! get the vector that points from the centroid to the ghost
1037 radial_vector = physical_loc - [patch_ib(patch_id)%x_centroid, patch_ib(patch_id)%y_centroid, &
1038 & patch_ib(patch_id)%z_centroid]
1039 ! convert the angular velocity from the inertial reference frame to the fluids frame, then convert to linear
1040 ! velocity
1041 call s_cross_product(patch_ib(patch_id)%angular_vel, radial_vector, rotation_velocity)
1042 do q = 1, 3
1043 ! if mibm is 1 or 2, then the boundary may be moving
1044 vel_g(q) = patch_ib(patch_id)%vel(q) ! add the linear velocity
1045 vel_g(q) = vel_g(q) + rotation_velocity(q) ! add the rotational velocity
1046 end do
1047 end if
1048 end if
1049
1050 ! Set momentum
1051
1052# 317 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1053#if defined(MFC_OpenACC)
1054# 317 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1055!$acc loop seq
1056# 317 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1057#elif defined(MFC_OpenMP)
1058# 317 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1059
1060# 317 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1061#endif
1062 do q = eqn_idx%mom%beg, eqn_idx%mom%end
1063 q_cons_vf(q)%sf(j, k, l) = rho*vel_g(q - eqn_idx%mom%beg + 1)
1064 dyn_pres = dyn_pres + q_cons_vf(q)%sf(j, k, l)*vel_g(q - eqn_idx%mom%beg + 1)/2._wp
1065 end do
1066
1067 ! Set continuity and adv vars
1068
1069# 324 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1070#if defined(MFC_OpenACC)
1071# 324 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1072!$acc loop seq
1073# 324 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1074#elif defined(MFC_OpenMP)
1075# 324 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1076
1077# 324 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1078#endif
1079 do q = 1, num_fluids
1080 q_cons_vf(q)%sf(j, k, l) = alpha_rho_ip(q)
1081 q_cons_vf(eqn_idx%adv%beg + q - 1)%sf(j, k, l) = alpha_ip(q)
1082 end do
1083
1084 ! Set color function
1085 if (surface_tension) then
1086 q_cons_vf(eqn_idx%c)%sf(j, k, l) = c_ip
1087 end if
1088
1089 ! Set Energy
1090 if (bubbles_euler) then
1091 q_cons_vf(eqn_idx%E)%sf(j, k, l) = (1 - alpha_ip(1))*(gamma*pres_ip + pi_inf + dyn_pres)
1092 else
1093 q_cons_vf(eqn_idx%E)%sf(j, k, l) = gamma*pres_ip + pi_inf + dyn_pres
1094 end if
1095 ! Set bubble vars
1096 if (bubbles_euler .and. .not. qbmm) then
1097 call s_comp_n_from_prim(alpha_ip(1), r_ip, nbub, weight)
1098
1099# 344 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1100#if defined(MFC_OpenACC)
1101# 344 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1102!$acc loop seq
1103# 344 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1104#elif defined(MFC_OpenMP)
1105# 344 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1106
1107# 344 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1108#endif
1109 do q = 1, nb
1110 q_cons_vf(eqn_idx%bub%beg + (q - 1)*2)%sf(j, k, l) = nbub*r_ip(q)
1111 q_cons_vf(eqn_idx%bub%beg + (q - 1)*2 + 1)%sf(j, k, l) = nbub*v_ip(q)
1112 if (.not. polytropic) then
1113 q_cons_vf(eqn_idx%bub%beg + (q - 1)*4)%sf(j, k, l) = nbub*r_ip(q)
1114 q_cons_vf(eqn_idx%bub%beg + (q - 1)*4 + 1)%sf(j, k, l) = nbub*v_ip(q)
1115 q_cons_vf(eqn_idx%bub%beg + (q - 1)*4 + 2)%sf(j, k, l) = nbub*pb_ip(q)
1116 q_cons_vf(eqn_idx%bub%beg + (q - 1)*4 + 3)%sf(j, k, l) = nbub*mv_ip(q)
1117 end if
1118 end do
1119 end if
1120
1121 if (qbmm) then
1122 nbub = nmom_ip(1)
1123
1124# 359 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1125#if defined(MFC_OpenACC)
1126# 359 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1127!$acc loop seq
1128# 359 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1129#elif defined(MFC_OpenMP)
1130# 359 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1131
1132# 359 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1133#endif
1134 do q = 1, nb*nmom
1135 q_cons_vf(eqn_idx%bub%beg + q - 1)%sf(j, k, l) = nbub*nmom_ip(q)
1136 end do
1137
1138
1139# 364 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1140#if defined(MFC_OpenACC)
1141# 364 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1142!$acc loop seq
1143# 364 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1144#elif defined(MFC_OpenMP)
1145# 364 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1146
1147# 364 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1148#endif
1149 do q = 1, nb
1150 q_cons_vf(eqn_idx%bub%beg + (q - 1)*nmom)%sf(j, k, l) = nbub
1151 end do
1152
1153 if (.not. polytropic) then
1154
1155# 370 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1156#if defined(MFC_OpenACC)
1157# 370 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1158!$acc loop seq
1159# 370 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1160#elif defined(MFC_OpenMP)
1161# 370 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1162
1163# 370 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1164#endif
1165 do q = 1, nb
1166
1167# 372 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1168#if defined(MFC_OpenACC)
1169# 372 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1170!$acc loop seq
1171# 372 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1172#elif defined(MFC_OpenMP)
1173# 372 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1174
1175# 372 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1176#endif
1177 do r = 1, nnode
1178 pb_in(j, k, l, r, q) = presb_ip((q - 1)*nnode + r)
1179 mv_in(j, k, l, r, q) = massv_ip((q - 1)*nnode + r)
1180 end do
1181 end do
1182 end if
1183 end if
1184
1185 if (model_eqns == 3) then
1186
1187# 382 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1188#if defined(MFC_OpenACC)
1189# 382 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1190!$acc loop seq
1191# 382 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1192#elif defined(MFC_OpenMP)
1193# 382 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1194
1195# 382 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1196#endif
1197 do q = eqn_idx%int_en%beg, eqn_idx%int_en%end
1198 q_cons_vf(q)%sf(j, k, &
1199 & l) = alpha_ip(q - eqn_idx%int_en%beg + 1)*(gammas(q - eqn_idx%int_en%beg + 1)*pres_ip &
1200 & + pi_infs(q - eqn_idx%int_en%beg + 1))
1201 end do
1202 end if
1203 end do
1204
1205# 390 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1206#if defined(MFC_OpenACC)
1207# 390 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1208!$acc end parallel loop
1209# 390 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1210#elif defined(MFC_OpenMP)
1211# 390 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1212
1213# 390 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1214!$omp end target teams loop
1215# 390 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1216#endif
1217 end if
1218
1219 end subroutine s_ibm_correct_state
1220
1221 !> Compute the image points for each ghost point
1222 impure subroutine s_compute_image_points(ghost_points_in)
1223
1224 type(ghost_point), dimension(num_gps), intent(inout) :: ghost_points_in
1225 real(wp) :: dist
1226 real(wp), dimension(3) :: norm
1227 real(wp), dimension(3) :: physical_loc
1228 real(wp) :: temp_loc
1229 real(wp), pointer, dimension(:) :: s_cc => null()
1230 integer :: bound
1231 type(ghost_point) :: gp
1232 integer :: q, dim !< Iterator variables
1233 integer :: i, j, k, l !< Location indexes
1234 integer :: patch_id !< IB Patch ID
1235 integer :: dir
1236 integer :: index
1237 logical :: bounds_error
1238
1239 bounds_error = .false.
1240
1241
1242# 415 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1243
1244# 415 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1245#if defined(MFC_OpenACC)
1246# 415 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1247!$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)
1248# 415 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1249#elif defined(MFC_OpenMP)
1250# 415 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1251
1252# 415 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1253
1254# 415 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1255
1256# 415 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1257!$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)
1258# 415 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1259#endif
1260# 417 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1261 do q = 1, num_gps
1262 gp = ghost_points_in(q)
1263 i = gp%loc(1)
1264 j = gp%loc(2)
1265 k = gp%loc(3)
1266
1267 ! Calculate physical location of ghost point
1268 if (p > 0) then
1269 physical_loc = [x_cc(i), y_cc(j), z_cc(k)]
1270 else
1271 physical_loc = [x_cc(i), y_cc(j), 0._wp]
1272 end if
1273
1274 ! Calculate and store the precise location of the image point
1275 patch_id = gp%ib_patch_id
1276 dist = abs(real(gp%levelset, kind=wp))
1277 norm(:) = gp%levelset_norm
1278 ghost_points_in(q)%ip_loc(:) = physical_loc(:) + 2*dist*norm(:)
1279
1280 ! Find the closest grid point to the image point
1281 do dim = 1, num_dims
1282 ! s_cc points to the dim array we need
1283 if (dim == 1) then
1284 s_cc => x_cc
1285 bound = m + buff_size - 1
1286 else if (dim == 2) then
1287 s_cc => y_cc
1288 bound = n + buff_size - 1
1289 else
1290 s_cc => z_cc
1291 bound = p + buff_size - 1
1292 end if
1293
1294 if (f_approx_equal(norm(dim), 0._wp)) then
1295 ! if the ghost point is almost equal to a cell location, we set it equal and continue
1296 ghost_points_in(q)%ip_grid(dim) = ghost_points_in(q)%loc(dim)
1297 else
1298 if (norm(dim) > 0) then
1299 dir = 1
1300 else
1301 dir = -1
1302 end if
1303
1304 index = ghost_points_in(q)%loc(dim)
1305 temp_loc = ghost_points_in(q)%ip_loc(dim)
1306 do while ((temp_loc < s_cc(index) .or. temp_loc > s_cc(index + 1)) .and. (.not. bounds_error))
1307 index = index + dir
1308 if (index < -buff_size .or. index > bound) then
1309#if !defined(MFC_OpenACC) && !defined(MFC_OpenMP)
1310 print *, "A required image point is not located in this computational domain."
1311 print *, "Ghost Point is located at :"
1312 if (p == 0) then
1313 print *, [x_cc(i), y_cc(j)]
1314 else
1315 print *, [x_cc(i), y_cc(j), z_cc(k)]
1316 end if
1317 print *, "We are searching in dimension ", dim, " for image point at ", ghost_points_in(q)%ip_loc(:)
1318 print *, "Domain size: ", [x_cc(-buff_size), y_cc(-buff_size), z_cc(-buff_size)]
1319 print *, "x: ", x_cc(-buff_size), " to: ", x_cc(m + buff_size - 1)
1320 print *, "y: ", y_cc(-buff_size), " to: ", y_cc(n + buff_size - 1)
1321 if (p /= 0) print *, "z: ", z_cc(-buff_size), " to: ", z_cc(p + buff_size - 1)
1322 print *, "Image point is located approximately ", &
1323 & (ghost_points_in(q)%loc(dim) - ghost_points_in(q) %ip_loc(dim))/(s_cc(1) - s_cc(0)), &
1324 & " grid cells away"
1325 print *, "Levelset ", dist, " and Norm: ", norm(:)
1326 print *, &
1327 & "A short term fix may include increasing buff_size further in m_helper_basic (currently set to a minimum of 10)"
1328#endif
1329 bounds_error = .true.
1330 end if
1331 end do
1332
1333 ghost_points_in(q)%ip_grid(dim) = index
1334 if (ghost_points_in(q)%DB(dim) == -1) then
1335 ghost_points_in(q)%ip_grid(dim) = ghost_points_in(q)%loc(dim) + 1
1336 else if (ghost_points_in(q)%DB(dim) == 1) then
1337 ghost_points_in(q)%ip_grid(dim) = ghost_points_in(q)%loc(dim) - 1
1338 end if
1339 end if
1340 end do
1341 end do
1342
1343# 498 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1344#if defined(MFC_OpenACC)
1345# 498 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1346!$acc end parallel loop
1347# 498 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1348#elif defined(MFC_OpenMP)
1349# 498 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1350
1351# 498 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1352!$omp end target teams loop
1353# 498 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1354#endif
1355
1356 if (bounds_error) then
1357# 500 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1358 call s_prohibit_abort("bounds_error", "Ghost Point and Image Point on Different Processors. Exiting")
1359# 500 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1360 end if
1361
1362 end subroutine s_compute_image_points
1363
1364 !> Count the number of ghost points for memory allocation
1365 subroutine s_find_num_ghost_points(num_gps_out)
1366
1367 integer, intent(out) :: num_gps_out
1368 integer :: i, j, k, ii, jj, kk, gp_layers_z !< Iterator variables
1369 integer :: num_gps_local !< local copies of the gp count to support GPU compute
1370 logical :: is_gp
1371
1372 num_gps_local = 0
1373 gp_layers_z = gp_layers
1374 if (p == 0) gp_layers_z = 0
1375
1376
1377# 516 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1378
1379# 516 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1380#if defined(MFC_OpenACC)
1381# 516 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1382!$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)
1383# 516 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1384#elif defined(MFC_OpenMP)
1385# 516 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1386
1387# 516 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1388
1389# 516 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1390
1391# 516 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1392!$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)
1393# 516 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1394#endif
1395# 518 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1396 do i = 0, m
1397 do j = 0, n
1398 do k = 0, p
1399 if (ib_markers%sf(i, j, k) /= 0) then
1400 is_gp = .false.
1401 marker_search: do ii = i - gp_layers, i + gp_layers
1402 do jj = j - gp_layers, j + gp_layers
1403 do kk = k - gp_layers_z, k + gp_layers_z
1404 if (ib_markers%sf(ii, jj, kk) == 0) then
1405 ! if any neighbors are not in the IB, it is a ghost point
1406 is_gp = .true.
1407 exit marker_search
1408 end if
1409 end do
1410 end do
1411 end do marker_search
1412
1413 if (is_gp) then
1414
1415# 536 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1416#if defined(MFC_OpenACC)
1417# 536 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1418!$acc atomic update
1419# 536 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1420#elif defined(MFC_OpenMP)
1421# 536 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1422!$omp atomic update
1423# 536 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1424#endif
1425 num_gps_local = num_gps_local + 1
1426 end if
1427 end if
1428 end do
1429 end do
1430 end do
1431
1432# 543 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1433#if defined(MFC_OpenACC)
1434# 543 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1435!$acc end parallel loop
1436# 543 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1437#elif defined(MFC_OpenMP)
1438# 543 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1439
1440# 543 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1441!$omp end target teams loop
1442# 543 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1443#endif
1444
1445 num_gps_out = num_gps_local
1446
1447 end subroutine s_find_num_ghost_points
1448
1449 !> Locate all ghost points in the domain
1450 subroutine s_find_ghost_points(ghost_points_in)
1451
1452 type(ghost_point), dimension(num_gps), intent(inout) :: ghost_points_in
1453 integer :: i, j, k, ii, jj, kk, gp_layers_z !< Iterator variables
1454 integer :: xp, yp, zp !< periodicities
1455 integer :: count, count_i, local_idx
1456 integer :: patch_id, encoded_patch_id, neighborhood_patch_id
1457 logical :: is_gp
1458
1459 count = 0
1460 count_i = 0
1461 gp_layers_z = gp_layers
1462 if (p == 0) gp_layers_z = 0
1463
1464
1465# 564 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1466
1467# 564 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1468#if defined(MFC_OpenACC)
1469# 564 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1470!$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, neighborhood_patch_id, xp, yp, zp) firstprivate(gp_layers, gp_layers_z) copyin(count, count_i, x_domain, y_domain, z_domain)
1471# 564 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1472#elif defined(MFC_OpenMP)
1473# 564 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1474
1475# 564 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1476
1477# 564 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1478
1479# 564 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1480!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, j, k, ii, jj, kk, is_gp, local_idx, patch_id, encoded_patch_id, neighborhood_patch_id, xp, yp, zp) firstprivate(gp_layers, gp_layers_z) map(to:count, count_i, x_domain, y_domain, z_domain)
1481# 564 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1482#endif
1483# 567 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1484 do i = 0, m
1485 do j = 0, n
1486 do k = 0, p
1487 if (ib_markers%sf(i, j, k) /= 0) then
1488 is_gp = .false.
1489 marker_search: do ii = i - gp_layers, i + gp_layers
1490 do jj = j - gp_layers, j + gp_layers
1491 do kk = k - gp_layers_z, k + gp_layers_z
1492 if (ib_markers%sf(ii, jj, kk) == 0) then
1493 ! if any neighbors are not in the IB, it is a ghost point
1494 is_gp = .true.
1495 exit marker_search
1496 end if
1497 end do
1498 end do
1499 end do marker_search
1500
1501 if (is_gp) then
1502
1503# 585 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1504#if defined(MFC_OpenACC)
1505# 585 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1506!$acc atomic capture
1507# 585 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1508#elif defined(MFC_OpenMP)
1509# 585 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1510!$omp atomic capture
1511# 585 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1512#endif
1513 count = count + 1
1514 local_idx = count
1515#if defined(MFC_OpenACC)
1516# 588 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1517!$acc end atomic
1518# 588 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1519#elif defined(MFC_OpenMP)
1520# 588 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1521!$omp end atomic
1522# 588 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1523#endif
1524
1525 ghost_points_in(local_idx)%loc = [i, j, k]
1526 encoded_patch_id = ib_markers%sf(i, j, k)
1527 call s_decode_patch_periodicity(encoded_patch_id, patch_id, xp, yp, zp)
1528 call s_get_neighborhood_idx(patch_id, neighborhood_patch_id)
1529 ghost_points_in(local_idx)%ib_patch_id = neighborhood_patch_id
1530 ghost_points_in(local_idx)%x_periodicity = xp
1531 ghost_points_in(local_idx)%y_periodicity = yp
1532 ghost_points_in(local_idx)%z_periodicity = zp
1533 ghost_points_in(local_idx)%slip = patch_ib(neighborhood_patch_id)%slip
1534
1535 if ((x_cc(i) - dx(i)) < x_domain%beg) then
1536 ghost_points_in(local_idx)%DB(1) = -1
1537 else if ((x_cc(i) + dx(i)) > x_domain%end) then
1538 ghost_points_in(local_idx)%DB(1) = 1
1539 else
1540 ghost_points_in(local_idx)%DB(1) = 0
1541 end if
1542
1543 if ((y_cc(j) - dy(j)) < y_domain%beg) then
1544 ghost_points_in(local_idx)%DB(2) = -1
1545 else if ((y_cc(j) + dy(j)) > y_domain%end) then
1546 ghost_points_in(local_idx)%DB(2) = 1
1547 else
1548 ghost_points_in(local_idx)%DB(2) = 0
1549 end if
1550
1551 if (p /= 0) then
1552 if ((z_cc(k) - dz(k)) < z_domain%beg) then
1553 ghost_points_in(local_idx)%DB(3) = -1
1554 else if ((z_cc(k) + dz(k)) > z_domain%end) then
1555 ghost_points_in(local_idx)%DB(3) = 1
1556 else
1557 ghost_points_in(local_idx)%DB(3) = 0
1558 end if
1559 end if
1560 end if
1561 end if
1562 end do
1563 end do
1564 end do
1565
1566# 630 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1567#if defined(MFC_OpenACC)
1568# 630 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1569!$acc end parallel loop
1570# 630 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1571#elif defined(MFC_OpenMP)
1572# 630 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1573
1574# 630 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1575!$omp end target teams loop
1576# 630 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1577#endif
1578
1579 end subroutine s_find_ghost_points
1580
1581 !> Compute the interpolation coefficients for image points
1582 subroutine s_compute_interpolation_coeffs(ghost_points_in)
1583
1584 type(ghost_point), dimension(num_gps), intent(inout) :: ghost_points_in
1585 real(wp), dimension(2, 2, 2) :: dist
1586 real(wp), dimension(2, 2, 2) :: alpha
1587 real(wp), dimension(2, 2, 2) :: interp_coeffs
1588 real(wp) :: buf
1589 real(wp), dimension(2, 2, 2) :: eta
1590 type(ghost_point) :: gp
1591 integer :: q, i, j, k, ii, jj, kk !< Grid indexes and iterators
1592 integer :: patch_id
1593 logical :: is_cell_center
1594
1595
1596# 648 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1597
1598# 648 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1599#if defined(MFC_OpenACC)
1600# 648 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1601!$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)
1602# 648 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1603#elif defined(MFC_OpenMP)
1604# 648 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1605
1606# 648 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1607
1608# 648 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1609
1610# 648 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1611!$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)
1612# 648 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1613#endif
1614 do q = 1, num_gps
1615 gp = ghost_points_in(q)
1616 ! Get the interpolation points
1617 i = gp%ip_grid(1)
1618 j = gp%ip_grid(2)
1619 if (p /= 0) then
1620 k = gp%ip_grid(3)
1621 else
1622 k = 0
1623 end if
1624
1625 ! get the distance to a cell in each direction
1626 dist = 0._wp
1627 buf = 1._wp
1628 do ii = 0, 1
1629 do jj = 0, 1
1630 if (p == 0) then
1631 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)
1632 else
1633 do kk = 0, 1
1634 dist(1 + ii, 1 + jj, &
1635 & 1 + kk) = sqrt((x_cc(i + ii) - gp%ip_loc(1))**2 + (y_cc(j + jj) - gp%ip_loc(2))**2 + (z_cc(k &
1636 & + kk) - gp%ip_loc(3))**2)
1637 end do
1638 end if
1639 end do
1640 end do
1641
1642 ! check if we are arbitrarily close to a cell center
1643 interp_coeffs = 0._wp
1644 is_cell_center = .false.
1645 check_is_cell_center: do ii = 0, 1
1646 do jj = 0, 1
1647 if (dist(ii + 1, jj + 1, 1) <= 1.e-16_wp) then
1648 interp_coeffs(ii + 1, jj + 1, 1) = 1._wp
1649 is_cell_center = .true.
1650 exit check_is_cell_center
1651 else
1652 if (p /= 0) then
1653 if (dist(ii + 1, jj + 1, 2) <= 1.e-16_wp) then
1654 interp_coeffs(ii + 1, jj + 1, 2) = 1._wp
1655 is_cell_center = .true.
1656 exit check_is_cell_center
1657 end if
1658 end if
1659 end if
1660 end do
1661 end do check_is_cell_center
1662
1663 if (.not. is_cell_center) then
1664 ! if we are not arbitrarily close, interpolate
1665 alpha = 1._wp
1666 patch_id = gp%ib_patch_id
1667 if (ib_markers%sf(i, j, k) /= 0) alpha(1, 1, 1) = 0._wp
1668 if (ib_markers%sf(i + 1, j, k) /= 0) alpha(2, 1, 1) = 0._wp
1669 if (ib_markers%sf(i, j + 1, k) /= 0) alpha(1, 2, 1) = 0._wp
1670 if (ib_markers%sf(i + 1, j + 1, k) /= 0) alpha(2, 2, 1) = 0._wp
1671
1672 if (p == 0) then
1673 eta(:,:,1) = 1._wp/dist(:,:,1)**2
1674 buf = sum(alpha(:,:,1)*eta(:,:,1))
1675 if (buf > 0._wp) then
1676 interp_coeffs(:,:,1) = alpha(:,:,1)*eta(:,:,1)/buf
1677 else
1678 buf = sum(eta(:,:,1))
1679 interp_coeffs(:,:,1) = eta(:,:,1)/buf
1680 end if
1681 else
1682 if (ib_markers%sf(i, j, k + 1) /= 0) alpha(1, 1, 2) = 0._wp
1683 if (ib_markers%sf(i + 1, j, k + 1) /= 0) alpha(2, 1, 2) = 0._wp
1684 if (ib_markers%sf(i, j + 1, k + 1) /= 0) alpha(1, 2, 2) = 0._wp
1685 if (ib_markers%sf(i + 1, j + 1, k + 1) /= 0) alpha(2, 2, 2) = 0._wp
1686 eta = 1._wp/dist**2
1687 buf = sum(alpha*eta)
1688
1689 if (buf > 0._wp) then
1690 interp_coeffs = alpha*eta/buf
1691 else
1692 buf = sum(eta)
1693 interp_coeffs = eta/buf
1694 end if
1695 end if
1696 end if
1697
1698 ghost_points_in(q)%interp_coeffs = interp_coeffs
1699 end do
1700
1701# 735 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1702#if defined(MFC_OpenACC)
1703# 735 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1704!$acc end parallel loop
1705# 735 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1706#elif defined(MFC_OpenMP)
1707# 735 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1708
1709# 735 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1710!$omp end target teams loop
1711# 735 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1712#endif
1713
1714 end subroutine s_compute_interpolation_coeffs
1715
1716 !> Interpolate primitive variables to a ghost point's image point using bilinear or trilinear interpolation
1717 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, &
1718 & nmom_IP, pb_in, mv_in, presb_IP, massv_IP)
1719
1720
1721# 743 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1722#if MFC_OpenACC
1723# 743 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1724!$acc routine seq
1725# 743 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1726#elif MFC_OpenMP
1727# 743 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1728
1729# 743 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1730
1731# 743 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1732!$omp declare target device_type(any)
1733# 743 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1734#endif
1735
1736 type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf !< Primitive Variables
1737 real(stp), optional, dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:,1:), intent(in) :: pb_in, mv_in
1738 type(ghost_point), intent(in) :: gp
1739 real(wp), intent(inout) :: pres_ip
1740 real(wp), dimension(3), intent(inout) :: vel_ip
1741 real(wp), intent(inout) :: c_ip
1742# 754 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1743 real(wp), dimension(num_fluids), intent(inout) :: alpha_ip, alpha_rho_ip
1744# 756 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1745 real(wp), optional, dimension(:), intent(inout) :: r_ip, v_ip, pb_ip, mv_ip
1746 real(wp), optional, dimension(:), intent(inout) :: nmom_ip
1747 real(wp), optional, dimension(:), intent(inout) :: presb_ip, massv_ip
1748 integer :: i, j, k, l, q !< Iterator variables
1749 integer :: i1, i2, j1, j2, k1, k2 !< Iterator variables
1750 real(wp) :: coeff
1751
1752 i1 = gp%ip_grid(1); i2 = i1 + 1
1753 j1 = gp%ip_grid(2); j2 = j1 + 1
1754 k1 = gp%ip_grid(3); k2 = k1 + 1
1755
1756 if (p == 0) then
1757 k1 = 0
1758 k2 = 0
1759 end if
1760
1761 alpha_rho_ip = 0._wp
1762 alpha_ip = 0._wp
1763 pres_ip = 0._wp
1764 vel_ip = 0._wp
1765
1766 if (surface_tension) c_ip = 0._wp
1767
1768 if (bubbles_euler) then
1769 r_ip = 0._wp
1770 v_ip = 0._wp
1771 if (.not. polytropic) then
1772 mv_ip = 0._wp
1773 pb_ip = 0._wp
1774 end if
1775 end if
1776
1777 if (qbmm) then
1778 nmom_ip = 0._wp
1779 if (.not. polytropic) then
1780 presb_ip = 0._wp
1781 massv_ip = 0._wp
1782 end if
1783 end if
1784
1785
1786# 796 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1787#if defined(MFC_OpenACC)
1788# 796 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1789!$acc loop seq
1790# 796 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1791#elif defined(MFC_OpenMP)
1792# 796 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1793
1794# 796 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1795#endif
1796 do i = i1, i2
1797
1798# 798 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1799#if defined(MFC_OpenACC)
1800# 798 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1801!$acc loop seq
1802# 798 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1803#elif defined(MFC_OpenMP)
1804# 798 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1805
1806# 798 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1807#endif
1808 do j = j1, j2
1809
1810# 800 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1811#if defined(MFC_OpenACC)
1812# 800 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1813!$acc loop seq
1814# 800 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1815#elif defined(MFC_OpenMP)
1816# 800 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1817
1818# 800 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1819#endif
1820 do k = k1, k2
1821 coeff = gp%interp_coeffs(i - i1 + 1, j - j1 + 1, k - k1 + 1)
1822
1823 pres_ip = pres_ip + coeff*q_prim_vf(eqn_idx%E)%sf(i, j, k)
1824
1825
1826# 806 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1827#if defined(MFC_OpenACC)
1828# 806 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1829!$acc loop seq
1830# 806 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1831#elif defined(MFC_OpenMP)
1832# 806 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1833
1834# 806 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1835#endif
1836 do q = eqn_idx%mom%beg, eqn_idx%mom%end
1837 vel_ip(q + 1 - eqn_idx%mom%beg) = vel_ip(q + 1 - eqn_idx%mom%beg) + coeff*q_prim_vf(q)%sf(i, j, k)
1838 end do
1839
1840
1841# 811 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1842#if defined(MFC_OpenACC)
1843# 811 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1844!$acc loop seq
1845# 811 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1846#elif defined(MFC_OpenMP)
1847# 811 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1848
1849# 811 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1850#endif
1851 do l = eqn_idx%cont%beg, eqn_idx%cont%end
1852 alpha_rho_ip(l) = alpha_rho_ip(l) + coeff*q_prim_vf(l)%sf(i, j, k)
1853 alpha_ip(l) = alpha_ip(l) + coeff*q_prim_vf(eqn_idx%adv%beg + l - 1)%sf(i, j, k)
1854 end do
1855
1856 if (surface_tension) then
1857 c_ip = c_ip + coeff*q_prim_vf(eqn_idx%c)%sf(i, j, k)
1858 end if
1859
1860 if (bubbles_euler .and. .not. qbmm) then
1861
1862# 822 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1863#if defined(MFC_OpenACC)
1864# 822 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1865!$acc loop seq
1866# 822 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1867#elif defined(MFC_OpenMP)
1868# 822 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1869
1870# 822 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1871#endif
1872 do l = 1, nb
1873 if (polytropic) then
1874 r_ip(l) = r_ip(l) + coeff*q_prim_vf(eqn_idx%bub%beg + (l - 1)*2)%sf(i, j, k)
1875 v_ip(l) = v_ip(l) + coeff*q_prim_vf(eqn_idx%bub%beg + 1 + (l - 1)*2)%sf(i, j, k)
1876 else
1877 r_ip(l) = r_ip(l) + coeff*q_prim_vf(eqn_idx%bub%beg + (l - 1)*4)%sf(i, j, k)
1878 v_ip(l) = v_ip(l) + coeff*q_prim_vf(eqn_idx%bub%beg + 1 + (l - 1)*4)%sf(i, j, k)
1879 pb_ip(l) = pb_ip(l) + coeff*q_prim_vf(eqn_idx%bub%beg + 2 + (l - 1)*4)%sf(i, j, k)
1880 mv_ip(l) = mv_ip(l) + coeff*q_prim_vf(eqn_idx%bub%beg + 3 + (l - 1)*4)%sf(i, j, k)
1881 end if
1882 end do
1883 end if
1884
1885 if (qbmm) then
1886 do l = 1, nb*nmom
1887 nmom_ip(l) = nmom_ip(l) + coeff*q_prim_vf(eqn_idx%bub%beg - 1 + l)%sf(i, j, k)
1888 end do
1889 if (.not. polytropic) then
1890 do q = 1, nb
1891 do l = 1, nnode
1892 presb_ip((q - 1)*nnode + l) = presb_ip((q - 1)*nnode + l) + coeff*real(pb_in(i, j, k, l, q), &
1893 & kind=wp)
1894 massv_ip((q - 1)*nnode + l) = massv_ip((q - 1)*nnode + l) + coeff*real(mv_in(i, j, k, l, q), &
1895 & kind=wp)
1896 end do
1897 end do
1898 end if
1899 end if
1900 end do
1901 end do
1902 end do
1903
1904 end subroutine s_interpolate_image_point
1905
1906 !> Resets the current indexes of immersed boundaries and replaces them after updating
1907 !> the position of each moving immersed boundary
1908 impure subroutine s_update_mib(num_ibs)
1909
1910 integer, intent(in) :: num_ibs
1911 integer :: i, j, k, z_gp_layers
1912
1913 call nvtxstartrange("UPDATE-MIBM")
1914
1915 ! Clears the existing immersed boundary indices
1916 z_gp_layers = 0; if (p /= 0) z_gp_layers = gp_layers + 1
1917
1918# 868 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1919
1920# 868 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1921#if defined(MFC_OpenACC)
1922# 868 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1923!$acc parallel loop gang vector default(present) private(i, j, k)
1924# 868 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1925#elif defined(MFC_OpenMP)
1926# 868 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1927
1928# 868 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1929
1930# 868 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1931
1932# 868 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1933!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, j, k)
1934# 868 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1935#endif
1936 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
1937 ib_markers%sf(i, j, k) = 0._wp
1938 end do; end do; end do
1939
1940# 872 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1941#if defined(MFC_OpenACC)
1942# 872 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1943!$acc end parallel loop
1944# 872 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1945#elif defined(MFC_OpenMP)
1946# 872 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1947
1948# 872 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1949!$omp end target teams loop
1950# 872 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1951#endif
1952
1953 ! recalulcate the rotation matrix based upon the new angles
1954
1955# 875 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1956
1957# 875 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1958#if defined(MFC_OpenACC)
1959# 875 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1960!$acc parallel loop gang vector default(present) private(i)
1961# 875 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1962#elif defined(MFC_OpenMP)
1963# 875 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1964
1965# 875 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1966
1967# 875 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1968
1969# 875 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1970!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i)
1971# 875 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1972#endif
1973 do i = 1, num_ibs
1974 if (patch_ib(i)%moving_ibm /= 0) then
1975 call s_update_ib_rotation_matrix(i)
1976 end if
1977 end do
1978
1979# 881 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1980#if defined(MFC_OpenACC)
1981# 881 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1982!$acc end parallel loop
1983# 881 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1984#elif defined(MFC_OpenMP)
1985# 881 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1986
1987# 881 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1988!$omp end target teams loop
1989# 881 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1990#endif
1991
1992 ! recompute the new ib_patch locations and broadcast them.
1993 call nvtxstartrange("APPLY-IB-PATCHES")
1994 call s_apply_ib_patches(ib_markers)
1995 call nvtxendrange
1996
1997 call nvtxstartrange("COMPUTE-GHOST-POINTS")
1998 ! recalculate the ghost point locations and coefficients
1999 call s_find_num_ghost_points(num_gps)
2000 call s_find_ghost_points(ghost_points)
2001 call nvtxendrange
2002
2003 call nvtxstartrange("COMPUTE-IMAGE-POINTS")
2004 call s_apply_levelset(ghost_points, num_gps)
2005 call s_compute_image_points(ghost_points)
2006 call s_compute_interpolation_coeffs(ghost_points)
2007 call nvtxendrange
2008
2009 call nvtxendrange
2010
2011 end subroutine s_update_mib
2012
2013 !> Compute pressure and viscous forces and torques on immersed bodies via volume integration
2014 subroutine s_compute_ib_forces(q_prim_vf, fluid_pp)
2015
2016 type(scalar_field), dimension(1:sys_size), intent(in) :: q_prim_vf
2017 type(physical_parameters), dimension(1:num_fluids), intent(in) :: fluid_pp
2018 integer :: i, j, k, l, encoded_ib_idx, ib_idx, ib_idx_temp, fluid_idx
2019 real(wp), dimension(num_ibs, 3) :: forces, torques
2020 ! viscous stress tensor with temp vectors to hold divergence calculations
2021 real(wp), dimension(1:3,1:3) :: viscous_stress_div, viscous_stress_div_1, viscous_stress_div_2
2022 real(wp), dimension(1:3) :: local_force_contribution, radial_vector, local_torque_contribution
2023 real(wp) :: cell_volume, dx, dy, dz, dynamic_viscosity
2024
2025# 919 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2026 real(wp), dimension(num_fluids) :: dynamic_viscosities
2027# 921 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2028
2029 call nvtxstartrange("COMPUTE-IB-FORCES")
2030
2031 forces = 0._wp
2032 torques = 0._wp
2033
2034 if (viscous) then
2035 do fluid_idx = 1, num_fluids
2036 if (fluid_pp(fluid_idx)%Re(1) > 0._wp) then
2037 dynamic_viscosities(fluid_idx) = 1._wp/fluid_pp(fluid_idx)%Re(1)
2038 else
2039 dynamic_viscosities(fluid_idx) = 0._wp
2040 end if
2041 end do
2042 end if
2043
2044
2045# 937 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2046
2047# 937 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2048#if defined(MFC_OpenACC)
2049# 937 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2050!$acc parallel loop collapse(3) gang vector default(present) private(ib_idx, ib_idx_temp, encoded_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(dynamic_viscosities)
2051# 937 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2052#elif defined(MFC_OpenMP)
2053# 937 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2054
2055# 937 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2056
2057# 937 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2058
2059# 937 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2060!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(ib_idx, ib_idx_temp, encoded_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:dynamic_viscosities)
2061# 937 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2062#endif
2063# 941 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2064 do i = 0, m
2065 do j = 0, n
2066 do k = 0, p
2067 encoded_ib_idx = ib_markers%sf(i, j, k)
2068 if (encoded_ib_idx /= 0) then
2069 call s_decode_patch_periodicity(encoded_ib_idx, ib_idx_temp)
2070 call s_get_neighborhood_idx(ib_idx_temp, ib_idx) ! global patch ID -> local index
2071 if (ib_idx > 0) then
2072 ! get the vector pointing to the grid cell from the IB centroid
2073 if (num_dims == 3) then
2074 radial_vector = [x_cc(i), y_cc(j), z_cc(k)] - [patch_ib(ib_idx)%x_centroid, &
2075 & patch_ib(ib_idx)%y_centroid, patch_ib(ib_idx)%z_centroid]
2076 else
2077 radial_vector = [x_cc(i), y_cc(j), 0._wp] - [patch_ib(ib_idx)%x_centroid, &
2078 & patch_ib(ib_idx)%y_centroid, 0._wp]
2079 end if
2080 dx = x_cc(i + 1) - x_cc(i)
2081 dy = y_cc(j + 1) - y_cc(j)
2082
2083 local_force_contribution(:) = 0._wp
2084 do fluid_idx = 0, num_fluids - 1
2085 ! Get the pressure contribution to force via a finite difference to compute the 2D components of the
2086 ! gradient of the pressure and cell volume
2087 local_force_contribution(1) = local_force_contribution(1) - (q_prim_vf(eqn_idx%E &
2088 & + fluid_idx)%sf(i + 1, j, &
2089 & k) - q_prim_vf(eqn_idx%E + fluid_idx)%sf(i - 1, j, k))/(2._wp*dx) ! force is the negative pressure gradient
2090 local_force_contribution(2) = local_force_contribution(2) - (q_prim_vf(eqn_idx%E &
2091 & + fluid_idx)%sf(i, j + 1, k) - q_prim_vf(eqn_idx%E + fluid_idx)%sf(i, &
2092 & j - 1, k))/(2._wp*dy)
2093 cell_volume = abs(dx*dy)
2094 ! add the 3D component of the pressure gradient, if we are working in 3 dimensions
2095 if (num_dims == 3) then
2096 dz = z_cc(k + 1) - z_cc(k)
2097 local_force_contribution(3) = local_force_contribution(3) - (q_prim_vf(eqn_idx%E &
2098 & + fluid_idx)%sf(i, j, &
2099 & k + 1) - q_prim_vf(eqn_idx%E + fluid_idx)%sf(i, j, k - 1))/(2._wp*dz)
2100 cell_volume = abs(cell_volume*dz)
2101 end if
2102 end do
2103
2104 ! get the viscous stress and add its contribution if that is considered
2105 if (viscous) then
2106 ! compute the volume-weighted local dynamic viscosity
2107 dynamic_viscosity = 0._wp
2108 do fluid_idx = 1, num_fluids
2109 ! local dynamic viscosity is the dynamic viscosity of the fluid times alpha of the fluid
2110 dynamic_viscosity = dynamic_viscosity + (q_prim_vf(fluid_idx + eqn_idx%adv%beg - 1)%sf(i, j, &
2111 & k)*dynamic_viscosities(fluid_idx))
2112 end do
2113
2114 ! get the linear force components first
2115 call s_compute_viscous_stress_tensor(viscous_stress_div_1, q_prim_vf, dynamic_viscosity, i - 1, &
2116 & j, k)
2117 call s_compute_viscous_stress_tensor(viscous_stress_div_2, q_prim_vf, dynamic_viscosity, i + 1, &
2118 & j, k)
2119 ! get x derivative of the first-row of viscous stress tensor
2120 viscous_stress_div(1,1:3) = (viscous_stress_div_2(1,1:3) - viscous_stress_div_1(1,1:3))/(2._wp*dx)
2121 ! add the x components of the divergence to the force
2122 local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(1,1:3)
2123
2124 call s_compute_viscous_stress_tensor(viscous_stress_div_1, q_prim_vf, dynamic_viscosity, i, &
2125 & j - 1, k)
2126 call s_compute_viscous_stress_tensor(viscous_stress_div_2, q_prim_vf, dynamic_viscosity, i, &
2127 & j + 1, k)
2128 ! get y derivative of the second-row of viscous stress tensor
2129 viscous_stress_div(2,1:3) = (viscous_stress_div_2(2,1:3) - viscous_stress_div_1(2,1:3))/(2._wp*dy)
2130 ! add the y components of the divergence to the force
2131 local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(2,1:3)
2132
2133 if (num_dims == 3) then
2134 call s_compute_viscous_stress_tensor(viscous_stress_div_1, q_prim_vf, dynamic_viscosity, i, &
2135 & j, k - 1)
2136 call s_compute_viscous_stress_tensor(viscous_stress_div_2, q_prim_vf, dynamic_viscosity, i, &
2137 & j, k + 1)
2138 viscous_stress_div(3,1:3) = (viscous_stress_div_2(3,1:3) - viscous_stress_div_1(3, &
2139 & 1:3))/(2._wp*dz)
2140 ! add the z components of the divergence to the force
2141 local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(3,1:3)
2142 end if
2143 end if
2144
2145 call s_cross_product(radial_vector, local_force_contribution, local_torque_contribution)
2146
2147 ! Update the force and torque values atomically to prevent race conditions
2148 do l = 1, 3
2149
2150# 1026 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2151#if defined(MFC_OpenACC)
2152# 1026 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2153!$acc atomic update
2154# 1026 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2155#elif defined(MFC_OpenMP)
2156# 1026 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2157!$omp atomic update
2158# 1026 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2159#endif
2160 forces(ib_idx, l) = forces(ib_idx, l) + (local_force_contribution(l)*cell_volume)
2161
2162# 1028 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2163#if defined(MFC_OpenACC)
2164# 1028 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2165!$acc atomic update
2166# 1028 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2167#elif defined(MFC_OpenMP)
2168# 1028 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2169!$omp atomic update
2170# 1028 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2171#endif
2172 torques(ib_idx, l) = torques(ib_idx, l) + local_torque_contribution(l)*cell_volume
2173 end do
2174 end if ! ib_idx > 0
2175 end if
2176 end do
2177 end do
2178 end do
2179
2180# 1036 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2181#if defined(MFC_OpenACC)
2182# 1036 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2183!$acc end parallel loop
2184# 1036 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2185#elif defined(MFC_OpenMP)
2186# 1036 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2187
2188# 1036 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2189!$omp end target teams loop
2190# 1036 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2191#endif
2192
2193 call s_apply_collision_forces(ghost_points, num_gps, ib_markers, forces, torques)
2194
2195 ! reduce the forces across local neighborhood ranks
2196 call s_communicate_ib_forces(forces, torques)
2197
2198 ! consider body forces after reducing to avoid double counting
2199 do i = 1, num_ibs
2200 if (bf_x) then
2201 forces(i, 1) = forces(i, 1) + accel_bf(1)*patch_ib(i)%mass
2202 end if
2203 if (bf_y) then
2204 forces(i, 2) = forces(i, 2) + accel_bf(2)*patch_ib(i)%mass
2205 end if
2206 if (bf_z) then
2207 forces(i, 3) = forces(i, 3) + accel_bf(3)*patch_ib(i)%mass
2208 end if
2209 end do
2210
2211 ! apply the summed forces
2212
2213# 1057 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2214
2215# 1057 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2216#if defined(MFC_OpenACC)
2217# 1057 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2218!$acc parallel loop gang vector default(present) private(i) copyin(forces, torques)
2219# 1057 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2220#elif defined(MFC_OpenMP)
2221# 1057 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2222
2223# 1057 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2224
2225# 1057 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2226
2227# 1057 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2228!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i) map(to:forces, torques)
2229# 1057 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2230#endif
2231 do i = 1, num_ibs
2232 patch_ib(i)%force(:) = forces(i,:)
2233 patch_ib(i)%torque(:) = torques(i,:)
2234 end do
2235
2236# 1062 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2237#if defined(MFC_OpenACC)
2238# 1062 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2239!$acc end parallel loop
2240# 1062 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2241#elif defined(MFC_OpenMP)
2242# 1062 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2243
2244# 1062 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2245!$omp end target teams loop
2246# 1062 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2247#endif
2248
2249 call nvtxendrange
2250
2251 end subroutine s_compute_ib_forces
2252
2253 !> Computes the center of mass for IB patch types where we are unable to determine their center of mass analytically.
2254 !> These patches include things like NACA airfoils and STL models
2255 subroutine s_compute_centroid_offset(ib_marker)
2256
2257 integer, intent(in) :: ib_marker
2258 integer :: i, j, k, num_cells, num_cells_local, decoded_gbl_id
2259 real(wp), dimension(1:3) :: center_of_mass, center_of_mass_local
2260
2261 ! Offset only needs to be computes for specific geometries
2262
2263 if (patch_ib(ib_marker)%geometry == 4 .or. patch_ib(ib_marker)%geometry == 5 .or. patch_ib(ib_marker)%geometry == 11 &
2264 & .or. patch_ib(ib_marker)%geometry == 12) then
2265 center_of_mass_local = [0._wp, 0._wp, 0._wp]
2266 num_cells_local = 0
2267
2268 ! get the summed mass distribution and number of cells to divide by
2269 do i = 0, m
2270 do j = 0, n
2271 do k = 0, p
2272 if (ib_markers%sf(i, j, k) /= 0) then
2273 call s_decode_patch_periodicity(ib_markers%sf(i, j, k), decoded_gbl_id)
2274 if (decoded_gbl_id == patch_ib(ib_marker)%gbl_patch_id) then
2275 num_cells_local = num_cells_local + 1
2276 center_of_mass_local = center_of_mass_local + [x_cc(i), y_cc(j), 0._wp]
2277 if (num_dims == 3) center_of_mass_local(3) = center_of_mass_local(3) + z_cc(k)
2278 end if
2279 end if
2280 end do
2281 end do
2282 end do
2283
2284 ! reduce the mass contribution over all MPI ranks and compute COM
2285 call s_mpi_allreduce_integer_sum(num_cells_local, num_cells)
2286 if (num_cells /= 0) then
2287 call s_mpi_allreduce_sum(center_of_mass_local(1), center_of_mass(1))
2288 call s_mpi_allreduce_sum(center_of_mass_local(2), center_of_mass(2))
2289 call s_mpi_allreduce_sum(center_of_mass_local(3), center_of_mass(3))
2290 center_of_mass = center_of_mass/real(num_cells, wp)
2291 else
2292 patch_ib(ib_marker)%centroid_offset = [0._wp, 0._wp, 0._wp]
2293 return
2294 end if
2295
2296 ! assign the centroid offset as a vector pointing from the true COM to the "centroid" in the input file and replace the
2297 ! current centroid
2298 patch_ib(ib_marker)%centroid_offset = [patch_ib(ib_marker)%x_centroid, patch_ib(ib_marker)%y_centroid, &
2299 & patch_ib(ib_marker)%z_centroid] - center_of_mass
2300 patch_ib(ib_marker)%x_centroid = center_of_mass(1)
2301 patch_ib(ib_marker)%y_centroid = center_of_mass(2)
2302 patch_ib(ib_marker)%z_centroid = center_of_mass(3)
2303
2304 ! rotate the centroid offset back into the local coords of the IB
2305 patch_ib(ib_marker)%centroid_offset = matmul(patch_ib(ib_marker)%rotation_matrix_inverse, &
2306 & patch_ib(ib_marker)%centroid_offset)
2307 else
2308 patch_ib(ib_marker)%centroid_offset(:) = [0._wp, 0._wp, 0._wp]
2309 end if
2310
2311 end subroutine s_compute_centroid_offset
2312
2313 !> Computes the moment of inertia for an immersed boundary
2314 subroutine s_compute_moment_of_inertia(ib_idx, axis)
2315
2316
2317# 1131 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2318#if MFC_OpenACC
2319# 1131 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2320!$acc routine seq
2321# 1131 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2322#elif MFC_OpenMP
2323# 1131 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2324
2325# 1131 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2326
2327# 1131 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2328!$omp declare target device_type(any)
2329# 1131 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2330#endif
2331
2332 real(wp), dimension(3), intent(in) :: axis !< the axis about which we compute the moment. Only required in 3D.
2333 integer, intent(in) :: ib_idx
2334 real(wp) :: moment, distance_to_axis, cell_volume
2335 real(wp), dimension(3) :: position, closest_point_along_axis, vector_to_axis, normal_axis
2336 integer :: i, j, k, count, ib_marker
2337
2338 if (p == 0) then
2339 normal_axis = [0, 0, 1]
2340 else if (sqrt(sum(axis**2)) < sgm_eps) then
2341 ! if the object is not actually rotating at this time, return a dummy value and exit
2342 patch_ib(ib_idx)%moment = 1._wp
2343 return
2344 else
2345 normal_axis = axis/sqrt(sum(axis))
2346 end if
2347
2348 ! if the IB is in 2D or a 3D sphere, we can compute this exactly
2349 if (patch_ib(ib_idx)%geometry == 2) then ! circle
2350 patch_ib(ib_idx)%moment = 0.5_wp*patch_ib(ib_idx)%mass*(patch_ib(ib_idx)%radius)**2
2351 else if (patch_ib(ib_idx)%geometry == 3) then ! rectangle
2352 patch_ib(ib_idx)%moment = patch_ib(ib_idx)%mass*(patch_ib(ib_idx)%length_x**2 + patch_ib(ib_idx) %length_y**2)/6._wp
2353 else if (patch_ib(ib_idx)%geometry == 6) then ! ellipse
2354 patch_ib(ib_idx)%moment = 0.0625_wp*patch_ib(ib_idx)%mass*(patch_ib(ib_idx)%length_x**2 + patch_ib(ib_idx) %length_y**2)
2355 else if (patch_ib(ib_idx)%geometry == 8) then ! sphere
2356 patch_ib(ib_idx)%moment = 0.4*patch_ib(ib_idx)%mass*(patch_ib(ib_idx)%radius)**2
2357 else ! we do not have an analytic moment of inertia calculation and need to approximate it directly via a sum
2358 count = 0
2359 moment = 0._wp
2360 cell_volume = (x_cc(1) - x_cc(0))*(y_cc(1) - y_cc(0))
2361 ! computed without grid stretching. Update in the loop to perform with stretching
2362 if (p /= 0) then
2363 cell_volume = cell_volume*(z_cc(1) - z_cc(0))
2364 end if
2365
2366 ib_marker = patch_ib(ib_idx)%gbl_patch_id
2367
2368 do i = 0, m
2369 do j = 0, n
2370 do k = 0, p
2371 if (ib_markers%sf(i, j, k) == ib_marker) then
2372 count = count + 1 ! increment the count of total cells in the boundary
2373
2374 ! get the position in local coordinates so that the axis passes through 0, 0, 0
2375 if (num_dims < 3) then
2376 position = [x_cc(i), y_cc(j), 0._wp] - [patch_ib(ib_idx)%x_centroid, patch_ib(ib_idx)%y_centroid, &
2377 & 0._wp]
2378 else
2379 position = [x_cc(i), y_cc(j), z_cc(k)] - [patch_ib(ib_idx)%x_centroid, &
2380 & patch_ib(ib_idx)%y_centroid, patch_ib(ib_idx)%z_centroid]
2381 end if
2382
2383 ! project the position along the axis to find the closest distance to the rotation axis
2384 closest_point_along_axis = normal_axis*dot_product(normal_axis, position)
2385 vector_to_axis = position - closest_point_along_axis
2386 distance_to_axis = dot_product(vector_to_axis, vector_to_axis) ! saves the distance to the axis squared
2387
2388 ! compute the position component of the moment
2389 moment = moment + distance_to_axis
2390 end if
2391 end do
2392 end do
2393 end do
2394
2395 ! write the final moment assuming the points are all uniform density
2396 patch_ib(ib_idx)%moment = moment*patch_ib(ib_idx)%mass/(count*cell_volume)
2397 end if
2398
2399 end subroutine s_compute_moment_of_inertia
2400
2401 !> Wrap immersed boundary positions across periodic domain boundaries
2402 subroutine s_wrap_periodic_ibs()
2403
2404 integer :: patch_id
2405
2406 do patch_id = 1, num_ibs
2407 ! check domain wraps in x, y,
2408# 1210 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2409 ! check for periodicity
2410 if (ib_bc_x%beg == bc_periodic) then
2411 ! check if the boundary has left the domain, and then correct
2412 if (patch_ib(patch_id)%x_centroid < x_domain%beg) then
2413 ! if the boundary exited "left", wrap it back around to the "right"
2414 patch_ib(patch_id)%x_centroid = patch_ib(patch_id)%x_centroid + (x_domain%end &
2415 & - x_domain%beg)
2416 else if (patch_ib(patch_id)%x_centroid > x_domain%end) then
2417 ! if the boundary exited "right", wrap it back around to the "left"
2418 patch_ib(patch_id)%x_centroid = patch_ib(patch_id)%x_centroid - (x_domain%end &
2419 & - x_domain%beg)
2420 end if
2421 end if
2422# 1210 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2423 ! check for periodicity
2424 if (ib_bc_y%beg == bc_periodic) then
2425 ! check if the boundary has left the domain, and then correct
2426 if (patch_ib(patch_id)%y_centroid < y_domain%beg) then
2427 ! if the boundary exited "left", wrap it back around to the "right"
2428 patch_ib(patch_id)%y_centroid = patch_ib(patch_id)%y_centroid + (y_domain%end &
2429 & - y_domain%beg)
2430 else if (patch_ib(patch_id)%y_centroid > y_domain%end) then
2431 ! if the boundary exited "right", wrap it back around to the "left"
2432 patch_ib(patch_id)%y_centroid = patch_ib(patch_id)%y_centroid - (y_domain%end &
2433 & - y_domain%beg)
2434 end if
2435 end if
2436# 1224 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2437
2438 if (p /= 0) then
2439 ! check for periodicity
2440 if (ib_bc_z%beg == bc_periodic) then
2441 ! check if the boundary has left the domain, and then correct
2442 if (patch_ib(patch_id)%z_centroid < z_domain%beg) then
2443 ! if the boundary exited "left", wrap it back around to the "right"
2444 patch_ib(patch_id)%z_centroid = patch_ib(patch_id)%z_centroid + (z_domain%end - z_domain%beg)
2445 else if (patch_ib(patch_id)%z_centroid > z_domain%end) then
2446 ! if the boundary exited "right", wrap it back around to the "left"
2447 patch_ib(patch_id)%z_centroid = patch_ib(patch_id)%z_centroid - (z_domain%end - z_domain%beg)
2448 end if
2449 end if
2450 end if
2451 end do
2452
2453 end subroutine s_wrap_periodic_ibs
2454
2455 !> @brief Swaps ownership of IBs and passes ownership of IBs to neighbor processors
2456 !> Reduces forces and torques across the local neighborhood without a global allreduce. Accumulation phase: 2 passes per
2457 !! dimension receiving from the low-index (-X) neighbor. Pass 1: add received values; save what was received as recv_snap. Pass
2458 !! 2: send current (post-pass-1) values; add received; subtract recv_snap to remove double-counting of the direct contribution
2459 !! already added in pass 1. Back-propagation phase: 2 passes per dimension receiving from the high-index (+X) neighbor, each
2460 !! overwriting local forces with the neighbor's accumulated total.
2461 subroutine s_communicate_ib_forces(forces, torques)
2462
2463 real(wp), dimension(num_ibs, 3), intent(inout) :: forces, torques
2464
2465#ifdef MFC_MPI
2466 integer :: i, j, k, pack_pos, unpack_pos, buf_size, ierr
2467 integer :: send_neighbor, recv_neighbor, recv_count, tag
2468 character(len=1), allocatable :: ib_force_send_buf(:), ib_force_recv_buf(:)
2469
2470 if (num_procs == 1) return
2471
2472 buf_size = storage_size(0)/8 + (storage_size(0)/8 + 6*storage_size(0._wp)/8)*size(patch_ib)
2473 allocate (ib_force_send_buf(buf_size), ib_force_recv_buf(buf_size))
2474
2475 ! Accumulation phase: propagate contributions toward the high-index corner.
2476# 1264 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2477 if (num_dims >= 1) then
2478 send_neighbor = merge(bc_x%end, mpi_proc_null, bc_x%end >= 0)
2479 recv_neighbor = merge(bc_x%beg, mpi_proc_null, bc_x%beg >= 0)
2480
2481 recv_forces_snap = 0._wp
2482 recv_torques_snap = 0._wp
2483 tag = 300
2484
2485 do k = 1, 2*ib_neighborhood_radius
2486 ! send forces to +x neighbor; receive from -x neighbor. Add received values then
2487 pack_pos = 0
2488
2489# 1275 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2490
2491# 1275 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2492#if defined(MFC_OpenACC)
2493# 1275 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2494!$acc parallel loop gang vector default(present) private(i) copyin(forces, torques)
2495# 1275 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2496#elif defined(MFC_OpenMP)
2497# 1275 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2498
2499# 1275 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2500
2501# 1275 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2502
2503# 1275 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2504!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i) map(to:forces, torques)
2505# 1275 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2506#endif
2507 do i = 1, num_ibs
2508 send_ids(i) = patch_ib(i)%gbl_patch_id
2509 send_ft(1:3,i) = forces(i,:)
2510 send_ft(4:6,i) = torques(i,:)
2511 end do
2512
2513# 1281 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2514#if defined(MFC_OpenACC)
2515# 1281 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2516!$acc end parallel loop
2517# 1281 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2518#elif defined(MFC_OpenMP)
2519# 1281 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2520
2521# 1281 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2522!$omp end target teams loop
2523# 1281 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2524#endif
2525
2526# 1282 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2527#if defined(MFC_OpenACC)
2528# 1282 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2529!$acc update host(send_ids, send_ft)
2530# 1282 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2531#elif defined(MFC_OpenMP)
2532# 1282 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2533!$omp target update from(send_ids, send_ft)
2534# 1282 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2535#endif
2536 call mpi_pack(num_ibs, 1, mpi_integer, ib_force_send_buf, buf_size, pack_pos, mpi_comm_world, ierr)
2537 call mpi_pack(send_ids, num_ibs, mpi_integer, ib_force_send_buf, buf_size, pack_pos, mpi_comm_world, ierr)
2538 call mpi_pack(send_ft, 6*num_ibs, mpi_p, ib_force_send_buf, buf_size, pack_pos, mpi_comm_world, ierr)
2539 call mpi_sendrecv(ib_force_send_buf, pack_pos, mpi_packed, send_neighbor, tag, ib_force_recv_buf, buf_size, &
2540 & mpi_packed, recv_neighbor, tag, mpi_comm_world, mpi_status_ignore, ierr)
2541
2542 if (recv_neighbor /= mpi_proc_null) then
2543 unpack_pos = 0
2544 call mpi_unpack(ib_force_recv_buf, buf_size, unpack_pos, recv_count, 1, mpi_integer, mpi_comm_world, ierr)
2545 call mpi_unpack(ib_force_recv_buf, buf_size, unpack_pos, recv_ids, recv_count, mpi_integer, &
2546 & mpi_comm_world, ierr)
2547 call mpi_unpack(ib_force_recv_buf, buf_size, unpack_pos, recv_ft, 6*recv_count, mpi_p, mpi_comm_world, ierr)
2548
2549# 1295 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2550
2551# 1295 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2552#if defined(MFC_OpenACC)
2553# 1295 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2554!$acc parallel loop gang vector default(present) private(i, j) copy(forces, torques, recv_forces_snap, recv_torques_snap) copyin(recv_ft, recv_ids)
2555# 1295 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2556#elif defined(MFC_OpenMP)
2557# 1295 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2558
2559# 1295 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2560
2561# 1295 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2562
2563# 1295 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2564!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, j) map(tofrom:forces, torques, recv_forces_snap, recv_torques_snap) map(to:recv_ft, recv_ids)
2565# 1295 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2566#endif
2567# 1297 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2568 do i = 1, recv_count
2569 call s_get_neighborhood_idx(recv_ids(i), j)
2570 if (j > 0) then
2571 ! add forces and subtract recv_snap prevent double-counting
2572 forces(j,:) = forces(j,:) + recv_ft(1:3,i) - recv_forces_snap(j,:)
2573 torques(j,:) = torques(j,:) + recv_ft(4:6,i) - recv_torques_snap(j,:)
2574 recv_forces_snap(j,:) = recv_ft(1:3,i)
2575 recv_torques_snap(j,:) = recv_ft(4:6,i)
2576 end if
2577 end do
2578
2579# 1307 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2580#if defined(MFC_OpenACC)
2581# 1307 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2582!$acc end parallel loop
2583# 1307 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2584#elif defined(MFC_OpenMP)
2585# 1307 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2586
2587# 1307 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2588!$omp end target teams loop
2589# 1307 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2590#endif
2591 end if
2592 tag = tag + 2
2593 end do
2594 end if
2595# 1264 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2596 if (num_dims >= 2) then
2597 send_neighbor = merge(bc_y%end, mpi_proc_null, bc_y%end >= 0)
2598 recv_neighbor = merge(bc_y%beg, mpi_proc_null, bc_y%beg >= 0)
2599
2600 recv_forces_snap = 0._wp
2601 recv_torques_snap = 0._wp
2602 tag = 300
2603
2604 do k = 1, 2*ib_neighborhood_radius
2605 ! send forces to +y neighbor; receive from -y neighbor. Add received values then
2606 pack_pos = 0
2607
2608# 1275 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2609
2610# 1275 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2611#if defined(MFC_OpenACC)
2612# 1275 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2613!$acc parallel loop gang vector default(present) private(i) copyin(forces, torques)
2614# 1275 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2615#elif defined(MFC_OpenMP)
2616# 1275 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2617
2618# 1275 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2619
2620# 1275 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2621
2622# 1275 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2623!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i) map(to:forces, torques)
2624# 1275 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2625#endif
2626 do i = 1, num_ibs
2627 send_ids(i) = patch_ib(i)%gbl_patch_id
2628 send_ft(1:3,i) = forces(i,:)
2629 send_ft(4:6,i) = torques(i,:)
2630 end do
2631
2632# 1281 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2633#if defined(MFC_OpenACC)
2634# 1281 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2635!$acc end parallel loop
2636# 1281 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2637#elif defined(MFC_OpenMP)
2638# 1281 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2639
2640# 1281 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2641!$omp end target teams loop
2642# 1281 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2643#endif
2644
2645# 1282 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2646#if defined(MFC_OpenACC)
2647# 1282 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2648!$acc update host(send_ids, send_ft)
2649# 1282 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2650#elif defined(MFC_OpenMP)
2651# 1282 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2652!$omp target update from(send_ids, send_ft)
2653# 1282 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2654#endif
2655 call mpi_pack(num_ibs, 1, mpi_integer, ib_force_send_buf, buf_size, pack_pos, mpi_comm_world, ierr)
2656 call mpi_pack(send_ids, num_ibs, mpi_integer, ib_force_send_buf, buf_size, pack_pos, mpi_comm_world, ierr)
2657 call mpi_pack(send_ft, 6*num_ibs, mpi_p, ib_force_send_buf, buf_size, pack_pos, mpi_comm_world, ierr)
2658 call mpi_sendrecv(ib_force_send_buf, pack_pos, mpi_packed, send_neighbor, tag, ib_force_recv_buf, buf_size, &
2659 & mpi_packed, recv_neighbor, tag, mpi_comm_world, mpi_status_ignore, ierr)
2660
2661 if (recv_neighbor /= mpi_proc_null) then
2662 unpack_pos = 0
2663 call mpi_unpack(ib_force_recv_buf, buf_size, unpack_pos, recv_count, 1, mpi_integer, mpi_comm_world, ierr)
2664 call mpi_unpack(ib_force_recv_buf, buf_size, unpack_pos, recv_ids, recv_count, mpi_integer, &
2665 & mpi_comm_world, ierr)
2666 call mpi_unpack(ib_force_recv_buf, buf_size, unpack_pos, recv_ft, 6*recv_count, mpi_p, mpi_comm_world, ierr)
2667
2668# 1295 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2669
2670# 1295 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2671#if defined(MFC_OpenACC)
2672# 1295 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2673!$acc parallel loop gang vector default(present) private(i, j) copy(forces, torques, recv_forces_snap, recv_torques_snap) copyin(recv_ft, recv_ids)
2674# 1295 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2675#elif defined(MFC_OpenMP)
2676# 1295 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2677
2678# 1295 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2679
2680# 1295 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2681
2682# 1295 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2683!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, j) map(tofrom:forces, torques, recv_forces_snap, recv_torques_snap) map(to:recv_ft, recv_ids)
2684# 1295 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2685#endif
2686# 1297 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2687 do i = 1, recv_count
2688 call s_get_neighborhood_idx(recv_ids(i), j)
2689 if (j > 0) then
2690 ! add forces and subtract recv_snap prevent double-counting
2691 forces(j,:) = forces(j,:) + recv_ft(1:3,i) - recv_forces_snap(j,:)
2692 torques(j,:) = torques(j,:) + recv_ft(4:6,i) - recv_torques_snap(j,:)
2693 recv_forces_snap(j,:) = recv_ft(1:3,i)
2694 recv_torques_snap(j,:) = recv_ft(4:6,i)
2695 end if
2696 end do
2697
2698# 1307 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2699#if defined(MFC_OpenACC)
2700# 1307 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2701!$acc end parallel loop
2702# 1307 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2703#elif defined(MFC_OpenMP)
2704# 1307 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2705
2706# 1307 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2707!$omp end target teams loop
2708# 1307 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2709#endif
2710 end if
2711 tag = tag + 2
2712 end do
2713 end if
2714# 1264 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2715 if (num_dims >= 3) then
2716 send_neighbor = merge(bc_z%end, mpi_proc_null, bc_z%end >= 0)
2717 recv_neighbor = merge(bc_z%beg, mpi_proc_null, bc_z%beg >= 0)
2718
2719 recv_forces_snap = 0._wp
2720 recv_torques_snap = 0._wp
2721 tag = 300
2722
2723 do k = 1, 2*ib_neighborhood_radius
2724 ! send forces to +z neighbor; receive from -z neighbor. Add received values then
2725 pack_pos = 0
2726
2727# 1275 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2728
2729# 1275 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2730#if defined(MFC_OpenACC)
2731# 1275 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2732!$acc parallel loop gang vector default(present) private(i) copyin(forces, torques)
2733# 1275 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2734#elif defined(MFC_OpenMP)
2735# 1275 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2736
2737# 1275 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2738
2739# 1275 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2740
2741# 1275 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2742!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i) map(to:forces, torques)
2743# 1275 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2744#endif
2745 do i = 1, num_ibs
2746 send_ids(i) = patch_ib(i)%gbl_patch_id
2747 send_ft(1:3,i) = forces(i,:)
2748 send_ft(4:6,i) = torques(i,:)
2749 end do
2750
2751# 1281 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2752#if defined(MFC_OpenACC)
2753# 1281 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2754!$acc end parallel loop
2755# 1281 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2756#elif defined(MFC_OpenMP)
2757# 1281 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2758
2759# 1281 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2760!$omp end target teams loop
2761# 1281 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2762#endif
2763
2764# 1282 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2765#if defined(MFC_OpenACC)
2766# 1282 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2767!$acc update host(send_ids, send_ft)
2768# 1282 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2769#elif defined(MFC_OpenMP)
2770# 1282 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2771!$omp target update from(send_ids, send_ft)
2772# 1282 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2773#endif
2774 call mpi_pack(num_ibs, 1, mpi_integer, ib_force_send_buf, buf_size, pack_pos, mpi_comm_world, ierr)
2775 call mpi_pack(send_ids, num_ibs, mpi_integer, ib_force_send_buf, buf_size, pack_pos, mpi_comm_world, ierr)
2776 call mpi_pack(send_ft, 6*num_ibs, mpi_p, ib_force_send_buf, buf_size, pack_pos, mpi_comm_world, ierr)
2777 call mpi_sendrecv(ib_force_send_buf, pack_pos, mpi_packed, send_neighbor, tag, ib_force_recv_buf, buf_size, &
2778 & mpi_packed, recv_neighbor, tag, mpi_comm_world, mpi_status_ignore, ierr)
2779
2780 if (recv_neighbor /= mpi_proc_null) then
2781 unpack_pos = 0
2782 call mpi_unpack(ib_force_recv_buf, buf_size, unpack_pos, recv_count, 1, mpi_integer, mpi_comm_world, ierr)
2783 call mpi_unpack(ib_force_recv_buf, buf_size, unpack_pos, recv_ids, recv_count, mpi_integer, &
2784 & mpi_comm_world, ierr)
2785 call mpi_unpack(ib_force_recv_buf, buf_size, unpack_pos, recv_ft, 6*recv_count, mpi_p, mpi_comm_world, ierr)
2786
2787# 1295 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2788
2789# 1295 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2790#if defined(MFC_OpenACC)
2791# 1295 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2792!$acc parallel loop gang vector default(present) private(i, j) copy(forces, torques, recv_forces_snap, recv_torques_snap) copyin(recv_ft, recv_ids)
2793# 1295 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2794#elif defined(MFC_OpenMP)
2795# 1295 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2796
2797# 1295 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2798
2799# 1295 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2800
2801# 1295 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2802!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, j) map(tofrom:forces, torques, recv_forces_snap, recv_torques_snap) map(to:recv_ft, recv_ids)
2803# 1295 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2804#endif
2805# 1297 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2806 do i = 1, recv_count
2807 call s_get_neighborhood_idx(recv_ids(i), j)
2808 if (j > 0) then
2809 ! add forces and subtract recv_snap prevent double-counting
2810 forces(j,:) = forces(j,:) + recv_ft(1:3,i) - recv_forces_snap(j,:)
2811 torques(j,:) = torques(j,:) + recv_ft(4:6,i) - recv_torques_snap(j,:)
2812 recv_forces_snap(j,:) = recv_ft(1:3,i)
2813 recv_torques_snap(j,:) = recv_ft(4:6,i)
2814 end if
2815 end do
2816
2817# 1307 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2818#if defined(MFC_OpenACC)
2819# 1307 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2820!$acc end parallel loop
2821# 1307 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2822#elif defined(MFC_OpenMP)
2823# 1307 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2824
2825# 1307 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2826!$omp end target teams loop
2827# 1307 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2828#endif
2829 end if
2830 tag = tag + 2
2831 end do
2832 end if
2833# 1313 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2834
2835 ! Send final sums back to neighbors in -X direction
2836# 1316 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2837 if (num_dims >= 1) then
2838 send_neighbor = merge(bc_x%beg, mpi_proc_null, bc_x%beg >= 0)
2839 recv_neighbor = merge(bc_x%end, mpi_proc_null, bc_x%end >= 0)
2840
2841 do k = 1, 2*ib_neighborhood_radius
2842 pack_pos = 0
2843
2844# 1322 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2845
2846# 1322 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2847#if defined(MFC_OpenACC)
2848# 1322 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2849!$acc parallel loop gang vector default(present) private(i) copyin(forces, torques)
2850# 1322 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2851#elif defined(MFC_OpenMP)
2852# 1322 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2853
2854# 1322 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2855
2856# 1322 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2857
2858# 1322 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2859!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i) map(to:forces, torques)
2860# 1322 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2861#endif
2862 do i = 1, num_ibs
2863 send_ids(i) = patch_ib(i)%gbl_patch_id
2864 send_ft(1:3,i) = forces(i,:)
2865 send_ft(4:6,i) = torques(i,:)
2866 end do
2867
2868# 1328 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2869#if defined(MFC_OpenACC)
2870# 1328 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2871!$acc end parallel loop
2872# 1328 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2873#elif defined(MFC_OpenMP)
2874# 1328 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2875
2876# 1328 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2877!$omp end target teams loop
2878# 1328 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2879#endif
2880
2881# 1329 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2882#if defined(MFC_OpenACC)
2883# 1329 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2884!$acc update host(send_ids, send_ft)
2885# 1329 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2886#elif defined(MFC_OpenMP)
2887# 1329 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2888!$omp target update from(send_ids, send_ft)
2889# 1329 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2890#endif
2891 call mpi_pack(num_ibs, 1, mpi_integer, ib_force_send_buf, buf_size, pack_pos, mpi_comm_world, ierr)
2892 call mpi_pack(send_ids, num_ibs, mpi_integer, ib_force_send_buf, buf_size, pack_pos, mpi_comm_world, ierr)
2893 call mpi_pack(send_ft, 6*num_ibs, mpi_p, ib_force_send_buf, buf_size, pack_pos, mpi_comm_world, ierr)
2894 call mpi_sendrecv(ib_force_send_buf, pack_pos, mpi_packed, send_neighbor, tag, ib_force_recv_buf, buf_size, &
2895 & mpi_packed, recv_neighbor, tag, mpi_comm_world, mpi_status_ignore, ierr)
2896 if (recv_neighbor /= mpi_proc_null) then
2897 unpack_pos = 0
2898 call mpi_unpack(ib_force_recv_buf, buf_size, unpack_pos, recv_count, 1, mpi_integer, mpi_comm_world, ierr)
2899 call mpi_unpack(ib_force_recv_buf, buf_size, unpack_pos, recv_ids, recv_count, mpi_integer, &
2900 & mpi_comm_world, ierr)
2901 call mpi_unpack(ib_force_recv_buf, buf_size, unpack_pos, recv_ft, 6*recv_count, mpi_p, mpi_comm_world, ierr)
2902
2903# 1341 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2904
2905# 1341 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2906#if defined(MFC_OpenACC)
2907# 1341 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2908!$acc parallel loop gang vector default(present) private(i, j) copy(forces, torques) copyin(recv_ft, recv_ids)
2909# 1341 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2910#elif defined(MFC_OpenMP)
2911# 1341 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2912
2913# 1341 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2914
2915# 1341 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2916
2917# 1341 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2918!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, j) map(tofrom:forces, torques) map(to:recv_ft, recv_ids)
2919# 1341 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2920#endif
2921 do i = 1, recv_count
2922 call s_get_neighborhood_idx(recv_ids(i), j)
2923 if (j > 0) then
2924 forces(j,:) = recv_ft(1:3,i)
2925 torques(j,:) = recv_ft(4:6,i)
2926 end if
2927 end do
2928
2929# 1349 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2930#if defined(MFC_OpenACC)
2931# 1349 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2932!$acc end parallel loop
2933# 1349 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2934#elif defined(MFC_OpenMP)
2935# 1349 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2936
2937# 1349 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2938!$omp end target teams loop
2939# 1349 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2940#endif
2941 end if
2942 tag = tag + 2
2943 end do
2944 end if
2945# 1316 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2946 if (num_dims >= 2) then
2947 send_neighbor = merge(bc_y%beg, mpi_proc_null, bc_y%beg >= 0)
2948 recv_neighbor = merge(bc_y%end, mpi_proc_null, bc_y%end >= 0)
2949
2950 do k = 1, 2*ib_neighborhood_radius
2951 pack_pos = 0
2952
2953# 1322 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2954
2955# 1322 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2956#if defined(MFC_OpenACC)
2957# 1322 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2958!$acc parallel loop gang vector default(present) private(i) copyin(forces, torques)
2959# 1322 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2960#elif defined(MFC_OpenMP)
2961# 1322 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2962
2963# 1322 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2964
2965# 1322 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2966
2967# 1322 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2968!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i) map(to:forces, torques)
2969# 1322 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2970#endif
2971 do i = 1, num_ibs
2972 send_ids(i) = patch_ib(i)%gbl_patch_id
2973 send_ft(1:3,i) = forces(i,:)
2974 send_ft(4:6,i) = torques(i,:)
2975 end do
2976
2977# 1328 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2978#if defined(MFC_OpenACC)
2979# 1328 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2980!$acc end parallel loop
2981# 1328 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2982#elif defined(MFC_OpenMP)
2983# 1328 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2984
2985# 1328 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2986!$omp end target teams loop
2987# 1328 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2988#endif
2989
2990# 1329 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2991#if defined(MFC_OpenACC)
2992# 1329 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2993!$acc update host(send_ids, send_ft)
2994# 1329 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2995#elif defined(MFC_OpenMP)
2996# 1329 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2997!$omp target update from(send_ids, send_ft)
2998# 1329 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2999#endif
3000 call mpi_pack(num_ibs, 1, mpi_integer, ib_force_send_buf, buf_size, pack_pos, mpi_comm_world, ierr)
3001 call mpi_pack(send_ids, num_ibs, mpi_integer, ib_force_send_buf, buf_size, pack_pos, mpi_comm_world, ierr)
3002 call mpi_pack(send_ft, 6*num_ibs, mpi_p, ib_force_send_buf, buf_size, pack_pos, mpi_comm_world, ierr)
3003 call mpi_sendrecv(ib_force_send_buf, pack_pos, mpi_packed, send_neighbor, tag, ib_force_recv_buf, buf_size, &
3004 & mpi_packed, recv_neighbor, tag, mpi_comm_world, mpi_status_ignore, ierr)
3005 if (recv_neighbor /= mpi_proc_null) then
3006 unpack_pos = 0
3007 call mpi_unpack(ib_force_recv_buf, buf_size, unpack_pos, recv_count, 1, mpi_integer, mpi_comm_world, ierr)
3008 call mpi_unpack(ib_force_recv_buf, buf_size, unpack_pos, recv_ids, recv_count, mpi_integer, &
3009 & mpi_comm_world, ierr)
3010 call mpi_unpack(ib_force_recv_buf, buf_size, unpack_pos, recv_ft, 6*recv_count, mpi_p, mpi_comm_world, ierr)
3011
3012# 1341 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3013
3014# 1341 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3015#if defined(MFC_OpenACC)
3016# 1341 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3017!$acc parallel loop gang vector default(present) private(i, j) copy(forces, torques) copyin(recv_ft, recv_ids)
3018# 1341 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3019#elif defined(MFC_OpenMP)
3020# 1341 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3021
3022# 1341 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3023
3024# 1341 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3025
3026# 1341 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3027!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, j) map(tofrom:forces, torques) map(to:recv_ft, recv_ids)
3028# 1341 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3029#endif
3030 do i = 1, recv_count
3031 call s_get_neighborhood_idx(recv_ids(i), j)
3032 if (j > 0) then
3033 forces(j,:) = recv_ft(1:3,i)
3034 torques(j,:) = recv_ft(4:6,i)
3035 end if
3036 end do
3037
3038# 1349 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3039#if defined(MFC_OpenACC)
3040# 1349 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3041!$acc end parallel loop
3042# 1349 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3043#elif defined(MFC_OpenMP)
3044# 1349 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3045
3046# 1349 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3047!$omp end target teams loop
3048# 1349 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3049#endif
3050 end if
3051 tag = tag + 2
3052 end do
3053 end if
3054# 1316 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3055 if (num_dims >= 3) then
3056 send_neighbor = merge(bc_z%beg, mpi_proc_null, bc_z%beg >= 0)
3057 recv_neighbor = merge(bc_z%end, mpi_proc_null, bc_z%end >= 0)
3058
3059 do k = 1, 2*ib_neighborhood_radius
3060 pack_pos = 0
3061
3062# 1322 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3063
3064# 1322 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3065#if defined(MFC_OpenACC)
3066# 1322 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3067!$acc parallel loop gang vector default(present) private(i) copyin(forces, torques)
3068# 1322 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3069#elif defined(MFC_OpenMP)
3070# 1322 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3071
3072# 1322 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3073
3074# 1322 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3075
3076# 1322 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3077!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i) map(to:forces, torques)
3078# 1322 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3079#endif
3080 do i = 1, num_ibs
3081 send_ids(i) = patch_ib(i)%gbl_patch_id
3082 send_ft(1:3,i) = forces(i,:)
3083 send_ft(4:6,i) = torques(i,:)
3084 end do
3085
3086# 1328 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3087#if defined(MFC_OpenACC)
3088# 1328 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3089!$acc end parallel loop
3090# 1328 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3091#elif defined(MFC_OpenMP)
3092# 1328 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3093
3094# 1328 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3095!$omp end target teams loop
3096# 1328 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3097#endif
3098
3099# 1329 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3100#if defined(MFC_OpenACC)
3101# 1329 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3102!$acc update host(send_ids, send_ft)
3103# 1329 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3104#elif defined(MFC_OpenMP)
3105# 1329 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3106!$omp target update from(send_ids, send_ft)
3107# 1329 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3108#endif
3109 call mpi_pack(num_ibs, 1, mpi_integer, ib_force_send_buf, buf_size, pack_pos, mpi_comm_world, ierr)
3110 call mpi_pack(send_ids, num_ibs, mpi_integer, ib_force_send_buf, buf_size, pack_pos, mpi_comm_world, ierr)
3111 call mpi_pack(send_ft, 6*num_ibs, mpi_p, ib_force_send_buf, buf_size, pack_pos, mpi_comm_world, ierr)
3112 call mpi_sendrecv(ib_force_send_buf, pack_pos, mpi_packed, send_neighbor, tag, ib_force_recv_buf, buf_size, &
3113 & mpi_packed, recv_neighbor, tag, mpi_comm_world, mpi_status_ignore, ierr)
3114 if (recv_neighbor /= mpi_proc_null) then
3115 unpack_pos = 0
3116 call mpi_unpack(ib_force_recv_buf, buf_size, unpack_pos, recv_count, 1, mpi_integer, mpi_comm_world, ierr)
3117 call mpi_unpack(ib_force_recv_buf, buf_size, unpack_pos, recv_ids, recv_count, mpi_integer, &
3118 & mpi_comm_world, ierr)
3119 call mpi_unpack(ib_force_recv_buf, buf_size, unpack_pos, recv_ft, 6*recv_count, mpi_p, mpi_comm_world, ierr)
3120
3121# 1341 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3122
3123# 1341 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3124#if defined(MFC_OpenACC)
3125# 1341 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3126!$acc parallel loop gang vector default(present) private(i, j) copy(forces, torques) copyin(recv_ft, recv_ids)
3127# 1341 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3128#elif defined(MFC_OpenMP)
3129# 1341 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3130
3131# 1341 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3132
3133# 1341 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3134
3135# 1341 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3136!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, j) map(tofrom:forces, torques) map(to:recv_ft, recv_ids)
3137# 1341 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3138#endif
3139 do i = 1, recv_count
3140 call s_get_neighborhood_idx(recv_ids(i), j)
3141 if (j > 0) then
3142 forces(j,:) = recv_ft(1:3,i)
3143 torques(j,:) = recv_ft(4:6,i)
3144 end if
3145 end do
3146
3147# 1349 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3148#if defined(MFC_OpenACC)
3149# 1349 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3150!$acc end parallel loop
3151# 1349 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3152#elif defined(MFC_OpenMP)
3153# 1349 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3154
3155# 1349 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3156!$omp end target teams loop
3157# 1349 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3158#endif
3159 end if
3160 tag = tag + 2
3161 end do
3162 end if
3163# 1355 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3164#endif
3165
3166 end subroutine s_communicate_ib_forces
3167
3168 subroutine s_handoff_ib_ownership()
3169
3170 integer :: i, j, k, output_idx, local_output_idx
3171 integer :: old_num_local_ibs
3172 integer :: new_count, recv_count
3173 integer :: pack_pos, unpack_pos, buf_size, patch_bytes
3174 integer :: send_neighbor, recv_neighbor, ierr
3175 integer :: dx, dy, dz, tag, nbr_idx, nreqs
3176 real(wp), dimension(3) :: centroid
3177 logical :: is_new
3178 type(ib_patch_parameters) :: tmp_patch
3179 integer, dimension(num_local_ibs_max) :: local_ib_idx_old
3180 ! 26 neighbors max in 3D (8 in 2D); each gets its own recv buffer
3181 integer, parameter :: max_nbrs = 26
3182 character(len=1), allocatable :: send_buf(:), recv_bufs(:,:)
3183 integer, dimension(2*max_nbrs) :: requests
3184 integer, dimension(max_nbrs) :: recv_neighbor_list
3185
3186#ifdef MFC_MPI
3187 if (num_procs > 1) then
3188 ! save a copy of the local IB's global indices to cross-reference for later.
3189 local_ib_idx_old = 0
3190 old_num_local_ibs = num_local_ibs
3191 do i = 1, num_local_ibs
3192 local_ib_idx_old(i) = patch_ib(local_ib_patch_ids(i))%gbl_patch_id
3193 end do
3194
3195 ! delete any particles that no longer need to be tracked and coalesce the array
3196 output_idx = 0
3197 local_output_idx = 0
3198 do i = 1, num_ibs
3199 centroid = [patch_ib(i)%x_centroid, patch_ib(i)%y_centroid, 0._wp]
3200 if (num_dims == 3) centroid(3) = patch_ib(i)%z_centroid
3201
3202 ! delete if not in neighborhood
3203 if (f_neighborhood_ranks_own_location(centroid)) then
3204 output_idx = output_idx + 1
3205 if (i /= output_idx) then
3206 patch_ib(output_idx) = patch_ib(i)
3207 end if
3208
3209 ! check if in local domain
3210 if (f_local_rank_owns_location(centroid)) then
3211 local_output_idx = local_output_idx + 1
3212 local_ib_patch_ids(local_output_idx) = output_idx
3213 end if
3214 end if
3215 end do
3216 num_ibs = output_idx
3217 num_local_ibs = local_output_idx
3218
3219# 1409 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3220#if defined(MFC_OpenACC)
3221# 1409 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3222!$acc update device(patch_ib)
3223# 1409 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3224#elif defined(MFC_OpenMP)
3225# 1409 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3226!$omp target update to(patch_ib)
3227# 1409 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3228#endif
3229 call s_update_ib_lookup()
3230
3231 ! Broadcast newly-owned patches to all neighborhood neighbors
3232 patch_bytes = storage_size(tmp_patch)/8
3233 buf_size = storage_size(0)/8 + patch_bytes*num_local_ibs_max
3234 allocate (send_buf(buf_size), recv_bufs(buf_size, max_nbrs))
3235
3236 ! Write placeholder count at position 0
3237 pack_pos = 0
3238 call mpi_pack(0, 1, mpi_integer, send_buf, buf_size, pack_pos, mpi_comm_world, ierr)
3239
3240 ! pack new patches and count them
3241 new_count = 0
3242 do i = 1, num_local_ibs
3243 k = local_ib_patch_ids(i)
3244 is_new = .true.
3245 do j = 1, old_num_local_ibs
3246 if (patch_ib(k)%gbl_patch_id == local_ib_idx_old(j)) then
3247 is_new = .false.
3248 exit
3249 end if
3250 end do
3251 if (is_new) then
3252 call mpi_pack(patch_ib(k), patch_bytes, mpi_byte, send_buf, buf_size, pack_pos, mpi_comm_world, ierr)
3253 new_count = new_count + 1
3254 end if
3255 end do
3256
3257 ! Overwrite the placeholder with the real count
3258 pack_pos = 0
3259 call mpi_pack(new_count, 1, mpi_integer, send_buf, buf_size, pack_pos, mpi_comm_world, ierr)
3260 pack_pos = storage_size(0)/8 + new_count*patch_bytes
3261
3262 ! Post all receives first, then sends
3263 nreqs = 0
3264 nbr_idx = 0
3265 do dz = merge(-1, 0, num_dims == 3), merge(1, 0, num_dims == 3)
3266 do dy = -1, 1
3267 do dx = -1, 1
3268 if (dx == 0 .and. dy == 0 .and. dz == 0) cycle
3269 nbr_idx = nbr_idx + 1
3270 tag = 200 + (dx + 1)*9 + (dy + 1)*3 + (dz + 1)
3271 recv_neighbor = ib_neighbor_ranks(-dx, -dy, -dz)
3272 recv_neighbor_list(nbr_idx) = mpi_proc_null
3273 if (recv_neighbor < 0) cycle
3274 recv_neighbor_list(nbr_idx) = recv_neighbor
3275 nreqs = nreqs + 1
3276 call mpi_irecv(recv_bufs(:,nbr_idx), buf_size, mpi_packed, recv_neighbor, tag, mpi_comm_world, &
3277 & requests(nreqs), ierr)
3278 end do
3279 end do
3280 end do
3281
3282 do dz = merge(-1, 0, num_dims == 3), merge(1, 0, num_dims == 3)
3283 do dy = -1, 1
3284 do dx = -1, 1
3285 if (dx == 0 .and. dy == 0 .and. dz == 0) cycle
3286 tag = 200 + (dx + 1)*9 + (dy + 1)*3 + (dz + 1)
3287 send_neighbor = ib_neighbor_ranks(dx, dy, dz)
3288 if (send_neighbor < 0) cycle
3289 nreqs = nreqs + 1
3290 call mpi_isend(send_buf, pack_pos, mpi_packed, send_neighbor, tag, mpi_comm_world, requests(nreqs), ierr)
3291 end do
3292 end do
3293 end do
3294
3295 call mpi_waitall(nreqs, requests, mpi_statuses_ignore, ierr)
3296
3297 ! Unpack all received buffers
3298 do nbr_idx = 1, merge(26, 8, num_dims == 3)
3299 if (recv_neighbor_list(nbr_idx) == mpi_proc_null) cycle
3300 unpack_pos = 0
3301 call mpi_unpack(recv_bufs(:,nbr_idx), buf_size, unpack_pos, recv_count, 1, mpi_integer, mpi_comm_world, ierr)
3302 do i = 1, recv_count
3303 call mpi_unpack(recv_bufs(:,nbr_idx), buf_size, unpack_pos, tmp_patch, patch_bytes, mpi_byte, mpi_comm_world, &
3304 & ierr)
3305 call s_get_neighborhood_idx(tmp_patch%gbl_patch_id, j)
3306 if (j < 0) then
3307 num_ibs = num_ibs + 1
3308 if (.not. (num_ibs <= size(patch_ib))) then
3309# 1489 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3310 call s_mpi_abort("m_ibm.fpp:1489: " // "Assertion failed: num_ibs <= size(patch_ib). " &
3311# 1489 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3312 & // 'patch_ib overflow in neighborhood handoff')
3313# 1489 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3314 end if
3315 patch_ib(num_ibs) = tmp_patch
3316 end if
3317 end do
3318 end do
3319
3320 deallocate (send_buf, recv_bufs)
3321
3322# 1496 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3323#if defined(MFC_OpenACC)
3324# 1496 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3325!$acc update device(patch_ib)
3326# 1496 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3327#elif defined(MFC_OpenMP)
3328# 1496 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3329!$omp target update to(patch_ib)
3330# 1496 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3331#endif
3332 call s_update_ib_lookup()
3333 end if
3334#endif
3335
3336 end subroutine s_handoff_ib_ownership
3337
3338 subroutine s_get_neighborhood_idx(gbl_idx, neighborhood_idx)
3339
3340
3341# 1505 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3342#if MFC_OpenACC
3343# 1505 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3344!$acc routine seq
3345# 1505 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3346#elif MFC_OpenMP
3347# 1505 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3348
3349# 1505 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3350
3351# 1505 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3352!$omp declare target device_type(any)
3353# 1505 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3354#endif
3355
3356 integer, intent(in) :: gbl_idx
3357 integer, intent(out) :: neighborhood_idx
3358 integer :: i
3359
3360 neighborhood_idx = ib_gbl_idx_lookup(gbl_idx)
3361
3362 end subroutine s_get_neighborhood_idx
3363
3364 subroutine s_update_ib_lookup()
3365
3366 integer :: i
3367
3368 ib_gbl_idx_lookup = -1
3369
3370# 1520 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3371#if defined(MFC_OpenACC)
3372# 1520 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3373!$acc update device(ib_gbl_idx_lookup)
3374# 1520 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3375#elif defined(MFC_OpenMP)
3376# 1520 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3377!$omp target update to(ib_gbl_idx_lookup)
3378# 1520 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3379#endif
3380
3381
3382# 1522 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3383
3384# 1522 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3385#if defined(MFC_OpenACC)
3386# 1522 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3387!$acc parallel loop gang vector default(present) private(i)
3388# 1522 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3389#elif defined(MFC_OpenMP)
3390# 1522 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3391
3392# 1522 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3393
3394# 1522 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3395
3396# 1522 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3397!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i)
3398# 1522 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3399#endif
3400 do i = 1, num_ibs
3401 ib_gbl_idx_lookup(patch_ib(i)%gbl_patch_id) = i
3402 end do
3403
3404# 1526 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3405#if defined(MFC_OpenACC)
3406# 1526 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3407!$acc end parallel loop
3408# 1526 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3409#elif defined(MFC_OpenMP)
3410# 1526 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3411
3412# 1526 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3413!$omp end target teams loop
3414# 1526 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3415#endif
3416
3417
3418# 1528 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3419#if defined(MFC_OpenACC)
3420# 1528 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3421!$acc update host(ib_gbl_idx_lookup)
3422# 1528 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3423#elif defined(MFC_OpenMP)
3424# 1528 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3425!$omp target update from(ib_gbl_idx_lookup)
3426# 1528 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3427#endif
3428
3429 end subroutine s_update_ib_lookup
3430
3431 !> Finalize the IBM module
3432 impure subroutine s_finalize_ibm_module()
3433
3434 integer :: i
3435
3436#ifdef MFC_DEBUG
3437# 1537 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3438 block
3439# 1537 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3440 use iso_fortran_env, only: output_unit
3441# 1537 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3442
3443# 1537 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3444 print *, 'm_ibm.fpp:1537: ', '@:DEALLOCATE(ib_markers%sf)'
3445# 1537 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3446
3447# 1537 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3448 call flush (output_unit)
3449# 1537 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3450 end block
3451# 1537 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3452#endif
3453# 1537 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3454
3455# 1537 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3456#if defined(MFC_OpenACC)
3457# 1537 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3458!$acc exit data delete(ib_markers%sf)
3459# 1537 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3460#elif defined(MFC_OpenMP)
3461# 1537 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3462!$omp target exit data map(release:ib_markers%sf)
3463# 1537 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3464#endif
3465# 1537 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3466 deallocate (ib_markers%sf)
3467#ifdef MFC_DEBUG
3468# 1538 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3469 block
3470# 1538 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3471 use iso_fortran_env, only: output_unit
3472# 1538 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3473
3474# 1538 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3475 print *, 'm_ibm.fpp:1538: ', '@:DEALLOCATE(ib_gbl_idx_lookup)'
3476# 1538 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3477
3478# 1538 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3479 call flush (output_unit)
3480# 1538 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3481 end block
3482# 1538 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3483#endif
3484# 1538 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3485
3486# 1538 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3487#if defined(MFC_OpenACC)
3488# 1538 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3489!$acc exit data delete(ib_gbl_idx_lookup)
3490# 1538 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3491#elif defined(MFC_OpenMP)
3492# 1538 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3493!$omp target exit data map(release:ib_gbl_idx_lookup)
3494# 1538 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3495#endif
3496# 1538 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3497 deallocate (ib_gbl_idx_lookup)
3498 do i = 1, num_ib_airfoils_max
3499 if (allocated(ib_airfoil_grids(i)%upper)) then
3500#ifdef MFC_DEBUG
3501# 1541 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3502 block
3503# 1541 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3504 use iso_fortran_env, only: output_unit
3505# 1541 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3506
3507# 1541 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3508 print *, 'm_ibm.fpp:1541: ', '@:DEALLOCATE(ib_airfoil_grids(i)%upper)'
3509# 1541 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3510
3511# 1541 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3512 call flush (output_unit)
3513# 1541 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3514 end block
3515# 1541 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3516#endif
3517# 1541 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3518
3519# 1541 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3520#if defined(MFC_OpenACC)
3521# 1541 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3522!$acc exit data delete(ib_airfoil_grids(i)%upper)
3523# 1541 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3524#elif defined(MFC_OpenMP)
3525# 1541 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3526!$omp target exit data map(release:ib_airfoil_grids(i)%upper)
3527# 1541 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3528#endif
3529# 1541 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3530 deallocate (ib_airfoil_grids(i)%upper)
3531#ifdef MFC_DEBUG
3532# 1542 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3533 block
3534# 1542 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3535 use iso_fortran_env, only: output_unit
3536# 1542 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3537
3538# 1542 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3539 print *, 'm_ibm.fpp:1542: ', '@:DEALLOCATE(ib_airfoil_grids(i)%lower)'
3540# 1542 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3541
3542# 1542 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3543 call flush (output_unit)
3544# 1542 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3545 end block
3546# 1542 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3547#endif
3548# 1542 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3549
3550# 1542 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3551#if defined(MFC_OpenACC)
3552# 1542 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3553!$acc exit data delete(ib_airfoil_grids(i)%lower)
3554# 1542 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3555#elif defined(MFC_OpenMP)
3556# 1542 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3557!$omp target exit data map(release:ib_airfoil_grids(i)%lower)
3558# 1542 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3559#endif
3560# 1542 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3561 deallocate (ib_airfoil_grids(i)%lower)
3562 end if
3563 end do
3564
3565 if (allocated(models)) then
3566 if (size(models) > 0) then
3567#ifdef MFC_DEBUG
3568# 1548 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3569 block
3570# 1548 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3571 use iso_fortran_env, only: output_unit
3572# 1548 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3573
3574# 1548 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3575 print *, 'm_ibm.fpp:1548: ', '@:DEALLOCATE(models)'
3576# 1548 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3577
3578# 1548 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3579 call flush (output_unit)
3580# 1548 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3581 end block
3582# 1548 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3583#endif
3584# 1548 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3585
3586# 1548 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3587#if defined(MFC_OpenACC)
3588# 1548 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3589!$acc exit data delete(models)
3590# 1548 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3591#elif defined(MFC_OpenMP)
3592# 1548 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3593!$omp target exit data map(release:models)
3594# 1548 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3595#endif
3596# 1548 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3597 deallocate (models)
3598 else
3599 deallocate (models)
3600 end if
3601 end if
3602
3603 if (collision_model > 0) call s_finalize_collisions_module()
3604
3605#ifdef MFC_MPI
3606 if (num_procs > 1) then
3607#ifdef MFC_DEBUG
3608# 1558 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3609 block
3610# 1558 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3611 use iso_fortran_env, only: output_unit
3612# 1558 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3613
3614# 1558 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3615 print *, 'm_ibm.fpp:1558: ', '@:DEALLOCATE(send_ids, send_ft)'
3616# 1558 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3617
3618# 1558 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3619 call flush (output_unit)
3620# 1558 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3621 end block
3622# 1558 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3623#endif
3624# 1558 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3625
3626# 1558 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3627#if defined(MFC_OpenACC)
3628# 1558 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3629!$acc exit data delete(send_ids, send_ft)
3630# 1558 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3631#elif defined(MFC_OpenMP)
3632# 1558 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3633!$omp target exit data map(release:send_ids, send_ft)
3634# 1558 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3635#endif
3636# 1558 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3637 deallocate (send_ids, send_ft)
3638 deallocate (recv_forces_snap, recv_torques_snap, recv_ids, recv_ft)
3639 end if
3640#endif
3641
3642 end subroutine s_finalize_ibm_module
3643
3644end module m_ibm
q_cons_vf
type(scalar_field), dimension(sys_size), intent(inout) q_cons_vf
Definition m_phase_change.fpp.f90:1065
k
integer, intent(in) k
Definition m_phase_change.fpp.f90:1067
j
integer, intent(in) j
Definition m_phase_change.fpp.f90:1067
l
integer, intent(in) l
Definition m_phase_change.fpp.f90:1067
m_collisions
Ghost-node immersed boundary method: locates ghost/image points, computes interpolation coefficients,...
Definition m_collisions.fpp.f90:318
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:967
m_global_parameters::p
integer p
Definition m_global_parameters.fpp.f90:419
m_global_parameters::n
integer n
Definition m_global_parameters.fpp.f90:404
m_global_parameters::m
integer m
Definition m_global_parameters.fpp.f90:398
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:1909
m_ibm::send_ft
real(wp), dimension(:,:), allocatable send_ft
Definition m_ibm.fpp.f90:408
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:2015
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:560
m_ibm::s_update_ib_lookup
subroutine s_update_ib_lookup()
Definition m_ibm.fpp.f90:3365
m_ibm::s_communicate_ib_forces
subroutine s_communicate_ib_forces(forces, torques)
Swaps ownership of IBs and passes ownership of IBs to neighbor processors Reduces forces and torques ...
Definition m_ibm.fpp.f90:2462
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:1583
m_ibm::s_handoff_ib_ownership
subroutine s_handoff_ib_ownership()
Definition m_ibm.fpp.f90:3169
m_ibm::s_finalize_ibm_module
impure subroutine, public s_finalize_ibm_module()
Finalize the IBM module.
Definition m_ibm.fpp.f90:3433
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:1719
m_ibm::ghost_points
type(ghost_point), dimension(:), allocatable ghost_points
Definition m_ibm.fpp.f90:352
m_ibm::num_gps
integer num_gps
Number of ghost points.
Definition m_ibm.fpp.f90:378
m_ibm::recv_forces_snap
real(wp), dimension(:,:), allocatable recv_forces_snap
Definition m_ibm.fpp.f90:409
m_ibm::moving_immersed_boundary_flag
logical moving_immersed_boundary_flag
Definition m_ibm.fpp.f90:404
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:1451
m_ibm::s_compute_moment_of_inertia
subroutine s_compute_moment_of_inertia(ib_idx, axis)
Computes the moment of inertia for an immersed boundary.
Definition m_ibm.fpp.f90:2315
m_ibm::s_get_neighborhood_idx
subroutine s_get_neighborhood_idx(gbl_idx, neighborhood_idx)
Definition m_ibm.fpp.f90:3339
m_ibm::recv_ft
real(wp), dimension(:,:), allocatable recv_ft
Definition m_ibm.fpp.f90:408
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:2256
m_ibm::ib_markers
type(integer_field), public ib_markers
Definition m_ibm.fpp.f90:339
m_ibm::recv_ids
integer, dimension(:), allocatable recv_ids
Definition m_ibm.fpp.f90:407
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:415
m_ibm::s_wrap_periodic_ibs
subroutine s_wrap_periodic_ibs()
Wrap immersed boundary positions across periodic domain boundaries.
Definition m_ibm.fpp.f90:2403
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:818
m_ibm::send_ids
integer, dimension(:), allocatable send_ids
Definition m_ibm.fpp.f90:407
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:1223
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:1366
m_ibm::recv_torques_snap
real(wp), dimension(:,:), allocatable recv_torques_snap
Definition m_ibm.fpp.f90:409
m_ibm::ib_gbl_idx_lookup
integer, dimension(:), allocatable, public ib_gbl_idx_lookup
Definition m_ibm.fpp.f90:365
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:733
m_derived_types::integer_field
Derived type annexing an integer scalar field (SF).
Definition m_derived_types.fpp.f90:341