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