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# 12 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
20# 1 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp" 1
21! Default for Fypp runs outside the toolchain (e.g. documentation). Searched last, so the
22! build/include/<target>/thermochem.fpp the toolchain writes from m_thermochem takes precedence.
23# 4 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp"
24# 5 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp"
25# 13 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 2
26
27# 16 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
28# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
29
30# 24 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
31
32# 53 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
33
34# 65 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
35
36# 75 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
37
38# 105 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
39
40# 117 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
41
42# 127 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
43
44# 174 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
45! New line at end of file is required for FYPP
46# 2 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
47# 1 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 1
48# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
49# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
50# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
51# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
52# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
53# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
54
55# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
56# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
57# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
58
59# 12 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
60# 1 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp" 1
61! Default for Fypp runs outside the toolchain (e.g. documentation). Searched last, so the
62! build/include/<target>/thermochem.fpp the toolchain writes from m_thermochem takes precedence.
63# 4 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp"
64# 5 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp"
65# 13 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 2
66
67# 16 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
68# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
69
70# 24 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
71
72# 53 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
73
74# 65 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
75
76# 75 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
77
78# 105 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
79
80# 117 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
81
82# 127 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
83
84# 174 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
85! New line at end of file is required for FYPP
86# 2 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 2
87
88# 4 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
89# 5 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
90# 6 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
91# 7 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
92# 8 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
93
94# 20 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
95
96# 43 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
97
98# 48 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
99
100# 53 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
101
102# 58 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
103
104# 63 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
105
106# 68 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
107
108# 81 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
109
110# 95 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
111
112# 100 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
113
114# 105 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
115
116# 110 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
117
118# 115 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
119
120# 120 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
121
122# 125 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
123
124# 130 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
125
126# 135 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
127
128# 165 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
129
130# 219 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
131
132# 233 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
133
134# 258 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
135
136# 269 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
137
138# 271 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
139# 282 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
140
141# 332 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
142
143# 342 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
144
145# 352 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
146
147# 361 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
148
149# 378 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
150
151# 388 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
152
153# 395 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
154
155# 401 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
156
157# 407 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
158
159# 413 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
160
161# 419 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
162
163# 425 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
164! New line at end of file is required for FYPP
165# 3 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
166# 1 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 1
167# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
168# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
169# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
170# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
171# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
172# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
173
174# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
175# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
176# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
177
178# 12 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
179# 1 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp" 1
180! Default for Fypp runs outside the toolchain (e.g. documentation). Searched last, so the
181! build/include/<target>/thermochem.fpp the toolchain writes from m_thermochem takes precedence.
182# 4 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp"
183# 5 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp"
184# 13 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 2
185
186# 16 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
187# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
188
189# 24 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
190
191# 53 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
192
193# 65 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
194
195# 75 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
196
197# 105 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
198
199# 117 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
200
201# 127 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
202
203# 174 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
204! New line at end of file is required for FYPP
205# 2 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 2
206
207# 7 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
208
209# 17 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
210
211# 22 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
212
213# 27 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
214
215# 32 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
216
217# 37 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
218
219# 42 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
220
221# 47 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
222
223# 52 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
224
225# 57 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
226
227# 62 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
228
229# 73 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
230
231# 78 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
232
233# 83 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
234
235# 88 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
236
237# 103 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
238
239# 131 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
240
241# 160 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
242
243# 175 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
244
245# 193 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
246
247# 215 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
248
249# 244 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
250
251# 259 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
252
253# 269 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
254
255# 278 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
256
257# 294 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
258
259# 304 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
260
261# 311 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
262! New line at end of file is required for FYPP
263# 4 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
264
265! GPU parallel region (scalar reductions, maxval/minval)
266# 23 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
267
268! GPU parallel loop over threads (most common GPU macro)
269# 43 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
270
271! Required closing for GPU_PARALLEL_LOOP
272# 55 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
273
274! Mark routine for device compilation
275# 112 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
276
277! Declare device-resident data
278# 130 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
279
280! Inner loop within a GPU parallel region
281# 145 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
282
283! Scoped GPU data region
284# 164 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
285
286! Host code with device pointers (for MPI with GPU buffers)
287# 193 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
288
289! Allocate device memory (unscoped)
290# 207 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
291
292! Free device memory
293# 219 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
294
295! Atomic operation on device
296# 231 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
297
298! End atomic capture block
299# 242 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
300
301! Copy data between host and device
302# 254 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
303
304! Synchronization barrier
305# 266 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
306
307! Import GPU library module (openacc or omp_lib)
308# 275 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
309
310! Emit code only for AMD compiler
311# 282 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
312
313! Emit code for non-Cray compilers
314# 289 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
315
316! Emit code only for Cray compiler
317# 296 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
318
319! Emit code for non-NVIDIA compilers
320# 303 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
321
322# 305 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
323# 306 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
324! New line at end of file is required for FYPP
325# 2 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 2
326
327# 14 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
328
329! Caution: This macro requires the use of a binding script to set CUDA_VISIBLE_DEVICES, such that we have one GPU device per MPI
330! rank. That's because for both cudaMemAdvise (preferred location) and cudaMemPrefetchAsync we use location = device_id = 0. For an
331! example see misc/nvidia_uvm/bind.sh.
332# 52 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
333
334! Allocate and create GPU device memory
335# 72 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
336
337! Free GPU device memory and deallocate
338# 80 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
339
340! Cray-specific GPU pointer setup for vector fields
341# 104 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
342
343! Cray-specific GPU pointer setup for scalar fields
344# 120 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
345
346! Cray-specific GPU pointer setup for acoustic source spatials
347# 145 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
348
349# 151 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
350
351# 158 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
352! New line at end of file is required for FYPP
353# 6 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp" 2
354
355!> @brief Ghost-node immersed boundary method: locates ghost/image points, computes interpolation coefficients, and corrects the
356!! flow state
357module m_ibm
358
359 use m_derived_types
360 use m_global_parameters
361 use m_mpi_proxy
362 use m_variables_conversion
363 use m_eos
364 use m_helper
365 use m_helper_basic
366 use m_constants
367 use m_compute_levelset
368 use m_ib_patches
369 use m_viscous
370 use m_model
371 use m_patch_geometries
372 use m_collisions
373 use m_thermochem, only: num_species, gas_constant, molecular_weights, get_mixture_molecular_weight, get_mixture_energy_mass, &
374 & get_mixture_thermal_conductivity_mixavg, get_species_mass_diffusivities_mixavg
375
376 use m_surface_thermochem, only: get_surface_net_production_rates, get_surface_reaction_heat_flux
377
378 implicit none
379
380 private :: s_compute_image_points, s_compute_interpolation_coeffs, s_interpolate_image_point, s_find_ghost_points, &
381 & s_find_num_ghost_points, s_compute_ghost_point_pressure, s_compute_ghost_point_velocity, s_interpolation_weights
382 ; public :: s_initialize_ibm_module, s_ibm_setup, s_ibm_correct_state, s_finalize_ibm_module, s_report_ibm_surface
383
384 !> Ghost points at which the reacting-surface Newton solve did not reach its tolerance, so the point fell back to a chemically
385 !! inert wall (keeping a prescribed Twall). Counted because that fallback is otherwise indistinguishable from a surface
386 !! mechanism that simply does nothing, and a run can look converged while the surface chemistry never engaged.
387 integer :: n_surface_not_converged = 0
388
389 !> Ghost points whose levelset distance was not positive, leaving the surface gradient (X_IP - X_s)/d undefined. Separate from
390 !! the counter above because this is a grid/geometry degeneracy, not a kinetics failure.
391 integer :: n_surface_ill_posed = 0
392 type(integer_field), public :: ib_markers
393
394# 45 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
395#if defined(MFC_OpenACC)
396# 45 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
397!$acc declare create(ib_markers)
398# 45 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
399#elif defined(MFC_OpenMP)
400# 45 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
401!$omp declare target (ib_markers)
402# 45 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
403#endif
404# 46 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
405
406 !> 1 at ghost points whose image point reaches the fluid, 0 elsewhere. Ghost points with a buried image point average over these
407 !! neighbors, and never over each other.
408 type(integer_field) :: corrected_gps
409
410# 50 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
411#if defined(MFC_OpenACC)
412# 50 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
413!$acc declare create(corrected_gps)
414# 50 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
415#elif defined(MFC_OpenMP)
416# 50 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
417!$omp declare target (corrected_gps)
418# 50 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
419#endif
420# 51 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
421
422 type(ghost_point), dimension(:), allocatable :: ghost_points
423
424# 53 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
425#if defined(MFC_OpenACC)
426# 53 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
427!$acc declare create(ghost_points)
428# 53 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
429#elif defined(MFC_OpenMP)
430# 53 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
431!$omp declare target (ghost_points)
432# 53 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
433#endif
434# 54 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
435
436 integer :: num_gps !< Number of ghost points
437#if defined(MFC_OpenACC)
438
439# 57 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
440#if defined(MFC_OpenACC)
441# 57 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
442!$acc declare create(gp_layers, num_gps)
443# 57 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
444#elif defined(MFC_OpenMP)
445# 57 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
446!$omp declare target (gp_layers, num_gps)
447# 57 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
448#endif
449# 58 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
450#elif defined(MFC_OpenMP)
451
452# 59 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
453#if defined(MFC_OpenACC)
454# 59 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
455!$acc declare create(num_gps)
456# 59 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
457#elif defined(MFC_OpenMP)
458# 59 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
459!$omp declare target (num_gps)
460# 59 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
461#endif
462# 60 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
463#endif
464 logical :: moving_immersed_boundary_flag
465 logical :: centroid_offsets_active = .false. !< some patch is a moving airfoil or STL, which carries a centroid_offset
466
467 ! IB force reduction (s_communicate_ib_forces): distinct neighborhood ranks in ascending order, and buffers allocated once
468 integer, parameter :: ib_rec_len = 7 !< gbl_patch_id, force(1:3), torque(1:3)
469 integer :: n_ib_nbrs = 0
470 integer, allocatable :: ib_nbrs(:), ib_send_counts(:), ib_reqs(:), ib_stats(:,:)
471 real(wp), allocatable :: ib_send_bufs(:,:), ib_recv_bufs(:,:)
472 private :: ib_rec_len, n_ib_nbrs, ib_nbrs, ib_send_counts, ib_reqs, ib_stats, ib_send_bufs, ib_recv_bufs
473
474contains
475
476 !> Allocates memory for the variables in the IBM module
477 impure subroutine s_initialize_ibm_module()
478
479 if (p > 0) then
480#ifdef MFC_DEBUG
481# 77 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
482 block
483# 77 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
484 use iso_fortran_env, only: output_unit
485# 77 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
486
487# 77 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
488 print *, 'm_ibm.fpp:77: ', '@:ALLOCATE(ib_markers%sf(-buff_size:m+buff_size, -buff_size:n+buff_size, -buff_size:p+buff_size))'
489# 77 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
490
491# 77 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
492 call flush (output_unit)
493# 77 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
494 end block
495# 77 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
496#endif
497# 77 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
498 allocate (ib_markers%sf(-buff_size:m+buff_size, -buff_size:n+buff_size, -buff_size:p+buff_size))
499# 77 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
500
501# 77 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
502
503# 77 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
504#if defined(MFC_OpenACC)
505# 77 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
506!$acc enter data create(ib_markers%sf)
507# 77 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
508#elif defined(MFC_OpenMP)
509# 77 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
510!$omp target enter data map(always,alloc:ib_markers%sf)
511# 77 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
512#endif
513# 78 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
514#ifdef MFC_DEBUG
515# 78 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
516 block
517# 78 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
518 use iso_fortran_env, only: output_unit
519# 78 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
520
521# 78 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
522 print *, 'm_ibm.fpp:78: ', '@:ALLOCATE(corrected_gps%sf(-buff_size:m+buff_size, -buff_size:n+buff_size, -buff_size:p+buff_size))'
523# 78 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
524
525# 78 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
526 call flush (output_unit)
527# 78 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
528 end block
529# 78 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
530#endif
531# 78 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
532 allocate (corrected_gps%sf(-buff_size:m+buff_size, -buff_size:n+buff_size, -buff_size:p+buff_size))
533# 78 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
534
535# 78 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
536
537# 78 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
538#if defined(MFC_OpenACC)
539# 78 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
540!$acc enter data create(corrected_gps%sf)
541# 78 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
542#elif defined(MFC_OpenMP)
543# 78 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
544!$omp target enter data map(always,alloc:corrected_gps%sf)
545# 78 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
546#endif
547# 79 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
548 else
549#ifdef MFC_DEBUG
550# 80 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
551 block
552# 80 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
553 use iso_fortran_env, only: output_unit
554# 80 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
555
556# 80 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
557 print *, 'm_ibm.fpp:80: ', '@:ALLOCATE(ib_markers%sf(-buff_size:m+buff_size, -buff_size:n+buff_size, 0:0))'
558# 80 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
559
560# 80 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
561 call flush (output_unit)
562# 80 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
563 end block
564# 80 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
565#endif
566# 80 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
567 allocate (ib_markers%sf(-buff_size:m+buff_size, -buff_size:n+buff_size, 0:0))
568# 80 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
569
570# 80 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
571
572# 80 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
573#if defined(MFC_OpenACC)
574# 80 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
575!$acc enter data create(ib_markers%sf)
576# 80 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
577#elif defined(MFC_OpenMP)
578# 80 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
579!$omp target enter data map(always,alloc:ib_markers%sf)
580# 80 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
581#endif
582# 81 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
583#ifdef MFC_DEBUG
584# 81 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
585 block
586# 81 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
587 use iso_fortran_env, only: output_unit
588# 81 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
589
590# 81 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
591 print *, 'm_ibm.fpp:81: ', '@:ALLOCATE(corrected_gps%sf(-buff_size:m+buff_size, -buff_size:n+buff_size, 0:0))'
592# 81 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
593
594# 81 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
595 call flush (output_unit)
596# 81 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
597 end block
598# 81 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
599#endif
600# 81 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
601 allocate (corrected_gps%sf(-buff_size:m+buff_size, -buff_size:n+buff_size, 0:0))
602# 81 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
603
604# 81 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
605
606# 81 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
607#if defined(MFC_OpenACC)
608# 81 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
609!$acc enter data create(corrected_gps%sf)
610# 81 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
611#elif defined(MFC_OpenMP)
612# 81 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
613!$omp target enter data map(always,alloc:corrected_gps%sf)
614# 81 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
615#endif
616# 82 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
617 end if
618
619#ifdef _CRAYFTN
620# 84 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
621 block
622# 84 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
623#ifdef MFC_DEBUG
624# 84 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
625 block
626# 84 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
627 use iso_fortran_env, only: output_unit
628# 84 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
629
630# 84 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
631 print *, 'm_ibm.fpp:84: ', '@:ACC_SETUP_SFs(ib_markers)'
632# 84 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
633
634# 84 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
635 call flush (output_unit)
636# 84 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
637 end block
638# 84 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
639#endif
640# 84 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
641
642# 84 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
643
644# 84 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
645#if defined(MFC_OpenACC)
646# 84 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
647!$acc enter data copyin(ib_markers)
648# 84 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
649#elif defined(MFC_OpenMP)
650# 84 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
651!$omp target enter data map(to:ib_markers)
652# 84 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
653#endif
654# 84 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
655 if (associated(ib_markers%sf)) then
656# 84 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
657
658# 84 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
659#if defined(MFC_OpenACC)
660# 84 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
661!$acc enter data copyin(ib_markers%sf)
662# 84 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
663#elif defined(MFC_OpenMP)
664# 84 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
665!$omp target enter data map(to:ib_markers%sf)
666# 84 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
667#endif
668# 84 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
669 end if
670# 84 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
671 end block
672# 84 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
673#endif
674# 85 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
675#ifdef _CRAYFTN
676# 85 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
677 block
678# 85 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
679#ifdef MFC_DEBUG
680# 85 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
681 block
682# 85 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
683 use iso_fortran_env, only: output_unit
684# 85 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
685
686# 85 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
687 print *, 'm_ibm.fpp:85: ', '@:ACC_SETUP_SFs(corrected_gps)'
688# 85 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
689
690# 85 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
691 call flush (output_unit)
692# 85 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
693 end block
694# 85 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
695#endif
696# 85 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
697
698# 85 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
699
700# 85 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
701#if defined(MFC_OpenACC)
702# 85 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
703!$acc enter data copyin(corrected_gps)
704# 85 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
705#elif defined(MFC_OpenMP)
706# 85 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
707!$omp target enter data map(to:corrected_gps)
708# 85 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
709#endif
710# 85 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
711 if (associated(corrected_gps%sf)) then
712# 85 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
713
714# 85 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
715#if defined(MFC_OpenACC)
716# 85 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
717!$acc enter data copyin(corrected_gps%sf)
718# 85 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
719#elif defined(MFC_OpenMP)
720# 85 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
721!$omp target enter data map(to:corrected_gps%sf)
722# 85 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
723#endif
724# 85 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
725 end if
726# 85 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
727 end block
728# 85 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
729#endif
730# 86 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
731
732
733# 87 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
734#if defined(MFC_OpenACC)
735# 87 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
736!$acc enter data copyin(num_gps)
737# 87 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
738#elif defined(MFC_OpenMP)
739# 87 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
740!$omp target enter data map(to:num_gps)
741# 87 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
742#endif
743# 88 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
744
745# 88 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
746#if defined(MFC_OpenACC)
747# 88 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
748!$acc update device(ib_second_order_vel, ib_ip_min_dist)
749# 88 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
750#elif defined(MFC_OpenMP)
751# 88 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
752!$omp target update to(ib_second_order_vel, ib_ip_min_dist)
753# 88 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
754#endif
755# 89 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
756
757 if (collision_model > 0) call s_initialize_collisions_module()
758
759 end subroutine s_initialize_ibm_module
760
761 !> Initializes the values of various IBM variables, such as ghost points and image points.
762 impure subroutine s_ibm_setup()
763
764 integer :: i, j, k, gid
765 integer(kind=8) :: n_need
766 real(wp) :: t_init !< initial time for prescribed kinematics
767 real(wp) :: max_num_gps_rank
768 integer(kind=8) :: max_num_gps
769
770 call nvtxstartrange("SETUP-IBM-MODULE")
771
772 ! GPU routines require updated cell centers
773
774# 106 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
775#if defined(MFC_OpenACC)
776# 106 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
777!$acc update device(num_ibs, num_gbl_ibs, x_cc, y_cc, dx, dy, ib_bc_x%beg, ib_bc_x%end, ib_bc_y%beg, ib_bc_y%end)
778# 106 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
779#elif defined(MFC_OpenMP)
780# 106 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
781!$omp target update to(num_ibs, num_gbl_ibs, x_cc, y_cc, dx, dy, ib_bc_x%beg, ib_bc_x%end, ib_bc_y%beg, ib_bc_y%end)
782# 106 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
783#endif
784# 107 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
785 if (p /= 0) then
786
787# 108 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
788#if defined(MFC_OpenACC)
789# 108 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
790!$acc update device(z_cc, dz, ib_bc_z%beg, ib_bc_z%end)
791# 108 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
792#elif defined(MFC_OpenMP)
793# 108 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
794!$omp target update to(z_cc, dz, ib_bc_z%beg, ib_bc_z%end)
795# 108 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
796#endif
797# 109 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
798 end if
799
800# 110 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
801#if defined(MFC_OpenACC)
802# 110 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
803!$acc update device(patch_ib(1:num_ibs), glb_bounds)
804# 110 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
805#elif defined(MFC_OpenMP)
806# 110 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
807!$omp target update to(patch_ib(1:num_ibs), glb_bounds)
808# 110 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
809#endif
810# 111 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
811
812 ! do all set up for moving immersed boundaries; prescribed kinematics are evaluated at the initial time so the
813 ! first stage already sees the correct body state, velocity and angular velocity
814 t_init = t_step_start*dt
815 if (cfl_dt) t_init = t_save*n_start
816
817# 116 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
818
819# 116 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
820#if defined(MFC_OpenACC)
821# 116 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
822!$acc parallel loop gang vector default(present) private(i) copyin(t_init)
823# 116 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
824#elif defined(MFC_OpenMP)
825# 116 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
826
827# 116 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
828
829# 116 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
830
831# 116 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
832!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i) map(to:t_init)
833# 116 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
834#endif
835# 117 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
836 do i = 1, num_ibs
837 if (patch_ib(i)%moving_ibm /= 0) then
838 if (patch_ib(i)%kin_model > 0) call s_prescribed_kinematics(i, t_init)
839 call s_compute_moment_of_inertia(patch_ib(i), patch_ib(i)%angular_vel, patch_ib(i)%moment)
840 end if
841 call s_update_ib_rotation_matrix(i)
842 end do
843
844# 124 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
845#if defined(MFC_OpenACC)
846# 124 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
847!$acc end parallel loop
848# 124 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
849#elif defined(MFC_OpenMP)
850# 124 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
851
852# 124 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
853!$omp end target teams loop
854# 124 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
855#endif
856# 125 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
857
858# 125 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
859#if defined(MFC_OpenACC)
860# 125 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
861!$acc update host(patch_ib(1:num_ibs))
862# 125 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
863#elif defined(MFC_OpenMP)
864# 125 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
865!$omp target update from(patch_ib(1:num_ibs))
866# 125 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
867#endif
868# 126 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
869
870 ! allocate some arrays for MPI communication, if required by this simulation
871#ifdef MFC_MPI
872 if (num_procs > 1) call s_setup_ib_force_comm()
873#endif
874
875 call s_update_ib_lookup()
876
877 ! recompute the new ib_patch locations
878 ib_markers%sf = 0._wp
879 corrected_gps%sf = 0
880
881# 137 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
882#if defined(MFC_OpenACC)
883# 137 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
884!$acc update device(ib_markers%sf, corrected_gps%sf)
885# 137 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
886#elif defined(MFC_OpenMP)
887# 137 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
888!$omp target update to(ib_markers%sf, corrected_gps%sf)
889# 137 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
890#endif
891# 138 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
892 call s_apply_ib_patches(ib_markers)
893
894# 139 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
895#if defined(MFC_OpenACC)
896# 139 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
897!$acc update host(ib_markers%sf)
898# 139 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
899#elif defined(MFC_OpenMP)
900# 139 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
901!$omp target update from(ib_markers%sf)
902# 139 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
903#endif
904# 140 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
905 call s_check_every_patch_marked()
906 ! One reduction decides whether any patch needs an offset, so cases without one (particle clouds have thousands of
907 ! patches) skip the per-patch collectives. Then loop over global ids: ranks hold different patches.
908 call s_mpi_allreduce_integer_sum(int(count([(f_needs_centroid_offset(patch_ib(i)), i=1, num_ibs)]), 8), n_need)
909 centroid_offsets_active = n_need > 0_8
910 if (centroid_offsets_active) then
911 do gid = 1, num_gbl_ibs
912 call s_get_neighborhood_idx(gid, i)
913 call s_compute_centroid_offset(gid, i)
914 end do
915 call s_restore_centroid_offsets(t_init)
916 end if
917 do i = 1, num_ibs
918
919# 153 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
920#if defined(MFC_OpenACC)
921# 153 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
922!$acc update device(patch_ib(i))
923# 153 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
924#elif defined(MFC_OpenMP)
925# 153 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
926!$omp target update to(patch_ib(i))
927# 153 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
928#endif
929# 154 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
930 end do
931
932 ! find the number of ghost points and set them to be the maximum total across ranks
933 call s_find_num_ghost_points(num_gps)
934 if (moving_immersed_boundary_flag) then
935 ! Every rank allocates the same size, which must hold any one rank's peak as the body moves: twice the global
936 ! total, but at least eight times the largest per-rank count (a thin body crossing few of many ranks)
937 call s_mpi_allreduce_integer_sum(int(num_gps, 8), max_num_gps)
938 call s_mpi_allreduce_max(real(num_gps, wp), max_num_gps_rank)
939 max_num_gps = max(max_num_gps*2_8, 8_8*int(max_num_gps_rank, 8))
940 max_num_gps = min(max_num_gps, int(m + 1, 8)*int(n + 1, 8)*int(p + 1, 8))
941 else
942 max_num_gps = int(num_gps, 8)
943 end if
944
945 ! set the size of the ghost point arrays to be the amount of points total, plus a factor of 2 buffer
946
947# 170 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
948#if defined(MFC_OpenACC)
949# 170 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
950!$acc update device(num_gps)
951# 170 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
952#elif defined(MFC_OpenMP)
953# 170 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
954!$omp target update to(num_gps)
955# 170 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
956#endif
957# 171 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
958#ifdef MFC_DEBUG
959# 171 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
960 block
961# 171 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
962 use iso_fortran_env, only: output_unit
963# 171 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
964
965# 171 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
966 print *, 'm_ibm.fpp:171: ', '@:ALLOCATE(ghost_points(1:max_num_gps))'
967# 171 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
968
969# 171 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
970 call flush (output_unit)
971# 171 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
972 end block
973# 171 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
974#endif
975# 171 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
976 allocate (ghost_points(1:max_num_gps))
977# 171 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
978
979# 171 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
980
981# 171 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
982#if defined(MFC_OpenACC)
983# 171 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
984!$acc enter data create(ghost_points)
985# 171 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
986#elif defined(MFC_OpenMP)
987# 171 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
988!$omp target enter data map(always,alloc:ghost_points)
989# 171 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
990#endif
991# 172 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
992
993
994# 173 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
995#if defined(MFC_OpenACC)
996# 173 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
997!$acc enter data copyin(ghost_points)
998# 173 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
999#elif defined(MFC_OpenMP)
1000# 173 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1001!$omp target enter data map(to:ghost_points)
1002# 173 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1003#endif
1004# 174 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1005 ! Ghost-cell IBM, Tseng & Ferziger JCP (2003), Mittal & Iaccarino ARFM (2005)
1006 call s_find_ghost_points()
1007 call s_apply_levelset(ghost_points, num_gps)
1008
1009 call s_compute_image_points()
1010 call s_compute_interpolation_coeffs()
1011
1012 call nvtxendrange
1013
1014 end subroutine s_ibm_setup
1015
1016 !> Pressure correction for a moving IB, accounting for the acceleration of the boundary surface. Clamped both ways: the
1017 !! linearization it comes from holds only while the correction is order one, and an unbounded one drives the ghost state to
1018 !! vacuum.
1019 subroutine s_compute_ghost_point_pressure(gp, gp_patch_id, rho, pres_IP, pres_GP)
1020
1021
1022# 190 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1023#if MFC_OpenACC
1024# 190 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1025!$acc routine seq
1026# 190 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1027#elif MFC_OpenMP
1028# 190 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1029
1030# 190 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1031
1032# 190 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1033!$omp declare target device_type(any)
1034# 190 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1035#endif
1036# 191 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1037
1038 type(ghost_point), intent(in) :: gp
1039 integer, intent(in) :: gp_patch_id
1040 real(wp), intent(in) :: rho, pres_ip
1041 real(wp), intent(out) :: pres_gp
1042
1043 pres_gp = pres_ip/min(max(1._wp - (abs(gp%levelset) + gp%ip_dist) &
1044 & *rho/pres_ip*dot_product(patch_ib(gp_patch_id)%force/patch_ib(gp_patch_id)%mass, &
1045 & gp%levelset_norm), 5.e-1_wp), 2._wp)
1046
1047 end subroutine s_compute_ghost_point_pressure
1048
1049 subroutine s_compute_ghost_point_velocity(gp, gp_patch_id, radial_vector, vel_IP, pres_IP, vel_GP)
1050
1051
1052# 205 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1053#if MFC_OpenACC
1054# 205 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1055!$acc routine seq
1056# 205 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1057#elif MFC_OpenMP
1058# 205 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1059
1060# 205 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1061
1062# 205 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1063!$omp declare target device_type(any)
1064# 205 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1065#endif
1066# 206 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1067
1068 type(ghost_point), intent(in) :: gp
1069 integer, intent(in) :: gp_patch_id
1070 real(wp), dimension(3), intent(in) :: radial_vector, vel_ip
1071 real(wp), intent(in) :: pres_ip
1072 real(wp), dimension(3), intent(out) :: vel_gp
1073 real(wp), dimension(3) :: norm, vel_norm_ip, rotation_velocity
1074 real(wp) :: buf, v_blow_eff
1075 integer :: q !< Iterator variable
1076
1077 ! Calculate velocity of ghost cell
1078 if (gp%slip) then
1079 norm(1:3) = gp%levelset_norm
1080 buf = sqrt(sum(norm**2))
1081 norm = norm/buf
1082 vel_norm_ip = sum(vel_ip*norm)*norm
1083 vel_gp = vel_ip - vel_norm_ip
1084 if (patch_ib(gp_patch_id)%moving_ibm /= 0) then
1085 ! compute the linear velocity of the ghost point due to rotation
1086 call s_cross_product(patch_ib(gp_patch_id)%angular_vel, radial_vector, rotation_velocity)
1087
1088 ! add only the component of the IB's motion that is normal to the surface
1089 vel_gp = vel_gp + sum((patch_ib(gp_patch_id)%vel + rotation_velocity)*norm)*norm
1090 end if
1091 else
1092 if (patch_ib(gp_patch_id)%moving_ibm == 0) then
1093 ! we know the object is not moving if moving_ibm is 0 (false)
1094 vel_gp = 0._wp
1095 else
1096 ! convert the angular velocity from the inertial reference frame to the fluids frame, then convert to linear
1097 ! velocity
1098 call s_cross_product(patch_ib(gp_patch_id)%angular_vel, radial_vector, rotation_velocity)
1099 do q = 1, 3
1100 ! if mibm is 1 or 2, then the boundary may be moving
1101 vel_gp(q) = patch_ib(gp_patch_id)%vel(q) ! add the linear velocity
1102 vel_gp(q) = vel_gp(q) + rotation_velocity(q) ! add the rotational velocity
1103 end do
1104 end if
1105 end if
1106
1107 ! Burning/injecting surface: superimpose wall-normal (outward) blowing on the
1108 ! ghost velocity so the immersed surface transpires/injects gas into the flow.
1109 if (patch_ib(gp_patch_id)%v_blow > 0._wp) then
1110 v_blow_eff = patch_ib(gp_patch_id)%v_blow
1111 ! Pressure-coupled burn rate (Vieille's law r_dot ~ p^n)
1112 if (patch_ib(gp_patch_id)%burn_rate_pref > 0._wp) then
1113 ! max(pres_IP, 0) guards the fractional power against a transient negative
1114 v_blow_eff = v_blow_eff*(max(pres_ip, &
1115 & 0._wp)/patch_ib(gp_patch_id)%burn_rate_pref)**patch_ib(gp_patch_id)%burn_rate_exp
1116 end if
1117 norm(1:3) = gp%levelset_norm
1118 buf = sqrt(sum(norm**2))
1119 if (buf > 0._wp) vel_gp = vel_gp + v_blow_eff*norm/buf
1120 end if
1121
1122 ! extrapolate linearly through the boundary intercept, Mittal et al. (2008)
1123 if (ib_second_order_vel) vel_gp = vel_gp - abs(gp%levelset)/max(gp%ip_dist, sgm_eps)*(vel_ip - vel_gp)
1124
1125 end subroutine s_compute_ghost_point_velocity
1126
1127 !> Update the conservative variables at the ghost points
1128 subroutine s_ibm_correct_state(q_cons_vf, q_prim_vf, pb_in, mv_in)
1129
1130 type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf !< Primitive Variables
1131 type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf !< Primitive Variables
1132 real(stp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:,1:), intent(inout) :: pb_in, mv_in
1133 integer :: i, j, k, l, q, r !< Iterator variables
1134 integer :: jj, kk, ll !< Neighbor iterators
1135 integer :: rad, rad_z, num_nbrs !< Neighbor stencil radius and population
1136 integer :: patch_id, patch_id_temp !< Patch ID of ghost point
1137 real(wp) :: rho, gamma, pi_inf, dyn_pres !< Mixture variables
1138 real(wp) :: vel_sum_g, e_ghost !< Ghost-point velocity magnitude and energy
1139 real(wp), dimension(2) :: re_k
1140 real(wp) :: g_k
1141 real(wp) :: qv_k
1142 real(wp) :: pres_ip, pres_gp
1143 real(wp), dimension(3) :: vel_ip
1144 real(wp) :: c_ip
1145
1146# 292 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1147 real(wp), dimension(num_fluids) :: gs
1148 real(wp), dimension(num_fluids) :: alpha_rho_ip, alpha_ip, alpha_rho_gp
1149 real(wp), dimension(nb) :: r_ip, v_ip, pb_ip, mv_ip
1150 real(wp), dimension(nb*nmom) :: nmom_ip
1151 real(wp), dimension(nb*nnode) :: presb_ip, massv_ip
1152# 298 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1153 real(wp), dimension(1) :: ys_ip, ys_g, ys_s, w_species
1154# 299 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1155 real(wp) :: alpha_q, alpha_rho_q, e_q
1156 real(wp) :: t_ip, mw_ip, e_ip !< Image-point temperature, mixture MW, and mass-specific internal energy (chemistry)
1157 real(wp) :: t_s, t_g, mw_s, mw_g, rho_s, rho_g, mdot_s, v_stefan, d
1158 logical :: surface_converged
1159 integer :: n_not_converged, n_ill_posed !< Per-call reacting-surface failure tallies (see the module-level counters)
1160 real(wp), dimension(3) :: norm !< Levelset normal at the ghost point, normalized below with buf
1161 ! Primitive variables at the image point associated with a ghost point, interpolated from surrounding fluid cells.
1162
1163 real(wp), dimension(3) :: physical_loc !< Physical loc of GP
1164 real(wp), dimension(3) :: vel_g !< Velocity of GP
1165 real(wp), dimension(3) :: radial_vector !< vector from centroid to ghost point
1166 real(wp) :: nbub
1167 real(wp) :: buf, buf_prim
1168 type(ghost_point) :: gp
1169
1170 ! set the Moving IBM interior conservative variables
1171
1172# 315 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1173
1174# 315 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1175#if defined(MFC_OpenACC)
1176# 315 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1177!$acc parallel loop collapse(3) gang vector default(present) private(i, j, k, patch_id, rho, patch_id_temp)
1178# 315 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1179#elif defined(MFC_OpenMP)
1180# 315 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1181
1182# 315 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1183
1184# 315 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1185
1186# 315 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1187!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) &
1188# 315 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1189!$omp& private(i, j, k, patch_id, rho, patch_id_temp)
1190# 315 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1191#endif
1192# 316 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1193 do l = 0, p
1194 do k = 0, n
1195 do j = 0, m
1196 patch_id = ib_markers%sf(j, k, l)
1197 if (patch_id /= 0) then
1198 call s_decode_patch_periodicity(patch_id, patch_id_temp)
1199 call s_get_neighborhood_idx(patch_id_temp, patch_id)
1200 if (patch_id > 0) then
1201 ! skip pressure correction with chemistry to prevent unphysical pressure
1202 if (.not. chemistry) q_prim_vf(eqn_idx%E)%sf(j, k, l) = 1._wp
1203 rho = 0._wp
1204 do i = 1, num_fluids
1205 rho = rho + q_prim_vf(eqn_idx%cont%beg + i - 1)%sf(j, k, l)
1206 end do
1207
1208 ! Sets the momentum
1209 do i = 1, num_dims
1210 q_cons_vf(eqn_idx%mom%beg + i - 1)%sf(j, k, l) = patch_ib(patch_id)%vel(i)*rho
1211 q_prim_vf(eqn_idx%mom%beg + i - 1)%sf(j, k, l) = patch_ib(patch_id)%vel(i)
1212 end do
1213 end if ! patch_id > 0
1214 end if
1215 end do
1216 end do
1217 end do
1218
1219# 341 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1220#if defined(MFC_OpenACC)
1221# 341 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1222!$acc end parallel loop
1223# 341 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1224#elif defined(MFC_OpenMP)
1225# 341 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1226
1227# 341 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1228!$omp end target teams loop
1229# 341 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1230#endif
1231# 342 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1232
1233 n_not_converged = 0
1234 n_ill_posed = 0
1235
1236 if (num_gps > 0) then
1237
1238# 347 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1239
1240# 347 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1241#if defined(MFC_OpenACC)
1242# 347 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1243!$acc parallel loop gang vector default(present) private(i, physical_loc, dyn_pres, alpha_rho_IP, alpha_IP, alpha_rho_GP, pres_IP, pres_GP, vel_IP, vel_g, r_IP, v_IP, pb_IP, mv_IP, nmom_IP, presb_IP, &
1244# 347 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1245!$acc& massv_IP, rho, gamma, pi_inf, Re_K, G_K, Gs, gp, norm, buf, radial_vector, j, k, l, q, qv_K, c_IP, nbub, patch_id, Ys_IP, W_species, T_IP, mw_IP, e_IP, Ys_g, Ys_s, T_s, T_g, mw_s, mw_g, rho_s, &
1246# 347 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1247!$acc& rho_g, mdot_s, v_stefan, d, surface_converged, vel_sum_g, E_ghost, alpha_q, alpha_rho_q, e_q) reduction(+:n_not_converged, n_ill_posed) present(ghost_points)
1248# 347 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1249#elif defined(MFC_OpenMP)
1250# 347 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1251
1252# 347 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1253
1254# 347 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1255
1256# 347 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1257!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, physical_loc, dyn_pres, &
1258# 347 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1259!$omp& alpha_rho_IP, alpha_IP, alpha_rho_GP, pres_IP, pres_GP, vel_IP, vel_g, r_IP, v_IP, pb_IP, mv_IP, nmom_IP, presb_IP, massv_IP, rho, gamma, pi_inf, Re_K, G_K, Gs, gp, norm, buf, radial_vector, j, &
1260# 347 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1261!$omp& k, l, q, qv_K, c_IP, nbub, patch_id, Ys_IP, W_species, T_IP, mw_IP, e_IP, Ys_g, Ys_s, T_s, T_g, mw_s, mw_g, rho_s, rho_g, mdot_s, v_stefan, d, surface_converged, vel_sum_g, E_ghost, alpha_q, &
1262# 347 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1263!$omp& alpha_rho_q, e_q) reduction(+:n_not_converged, n_ill_posed)
1264# 347 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1265#endif
1266# 353 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1267 do i = 1, num_gps
1268 gp = ghost_points(i)
1269 if (.not. gp%interp_valid) cycle
1270 j = gp%loc(1)
1271 k = gp%loc(2)
1272 l = gp%loc(3)
1273 patch_id = ghost_points(i)%ib_patch_id
1274
1275 ! Calculate physical location of GP
1276 physical_loc = [x_cc(j), y_cc(k), 0._wp]
1277 if (num_dims == 3) physical_loc(3) = z_cc(l)
1278
1279 ! Interpolate primitive variables at image point associated w/ GP
1280 if (bubbles_euler .and. .not. qbmm) then
1281 call s_interpolate_image_point(q_prim_vf, gp, alpha_rho_ip, alpha_ip, pres_ip, vel_ip, c_ip, r_ip, v_ip, &
1282 & pb_ip, mv_ip)
1283 else if (qbmm .and. polytropic) then
1284 call s_interpolate_image_point(q_prim_vf, gp, alpha_rho_ip, alpha_ip, pres_ip, vel_ip, c_ip, r_ip, v_ip, &
1285 & pb_ip, mv_ip, nmom_ip)
1286 else if (qbmm .and. .not. polytropic) then
1287 call s_interpolate_image_point(q_prim_vf, gp, alpha_rho_ip, alpha_ip, pres_ip, vel_ip, c_ip, r_ip, v_ip, &
1288 & pb_ip, mv_ip, nmom_ip, pb_in, mv_in, presb_ip, massv_ip)
1289 else if (chemistry) then
1290 call s_interpolate_image_point(q_prim_vf, gp, alpha_rho_ip, alpha_ip, pres_ip, vel_ip, c_ip, ys_ip=ys_ip)
1291 else
1292 call s_interpolate_image_point(q_prim_vf, gp, alpha_rho_ip, alpha_ip, pres_ip, vel_ip, c_ip)
1293 end if
1294
1295 ! Injecting (burning) surface: replace the mirrored ghost composition with pure
1296 ! injected fuel at the local pressure and the ambient (image-point) temperature.
1297 ! Setting a consistent injected density here (rather than reusing the heavy ambient
1298 ! rho) keeps the light fuel at a physical temperature and feeds the surface flame.
1299 if (chemistry .and. patch_ib(patch_id)%inj_species > 0) then
1300 call get_mixture_molecular_weight(ys_ip, mw_ip)
1301 t_ip = pres_ip*mw_ip/(alpha_rho_ip(1)*gas_constant)
1302 ys_ip = 0._wp
1303 ys_ip(patch_ib(patch_id)%inj_species) = 1._wp
1304 call get_mixture_molecular_weight(ys_ip, mw_ip)
1305 alpha_rho_ip(1) = pres_ip*mw_ip/(t_ip*gas_constant)
1306 end if
1307
1308 ! Thermal and heterogeneous reacting-surface boundary conditions.
1309 v_stefan = 0._wp
1310 surface_converged = .false.
1311
1312 if (chemistry .and. patch_ib(patch_id)%inj_species == 0) then
1313 ! Intrinsic gas state at the image point: rho_IP = (alpha*rho)_IP / alpha_IP
1314 call get_mixture_molecular_weight(ys_ip, mw_ip)
1315 t_ip = pres_ip*mw_ip*alpha_ip(1)/(alpha_rho_ip(1)*gas_constant)
1316
1317 if (patch_ib(patch_id)%surface_reaction == 1) then
1318 ! Heterogeneous reacting surface.
1319 d = abs(real(gp%levelset, kind=wp))
1320
1321 w_species(1:num_species) = molecular_weights(:)
1322
1323 ! d is the gas-side gradient length for every surface flux, so a ghost point sitting exactly on the
1324 ! immersed surface (grid-aligned rectangle patch, coincident STL facet) leaves the whole balance
1325 ! undefined. Left alone the 1/d produces NaN residuals, and a NaN then defeats the pivot test in
1326 ! s_solve_surface_linear_system -- NaN <= epsilon is false -- so the solve reports success and writes
1327 ! NaN into the ghost state. Such a point is not solvable, so it is recorded and skipped rather than
1328 ! having its distance rescaled to something it is not.
1329 if (d > 0._wp) then
1330 call s_solve_surface(pres_ip, t_ip, patch_ib(patch_id)%Twall, d, ys_ip, w_species, &
1331 & patch_ib(patch_id)%thermal_bc, ys_s, t_s, mdot_s, surface_converged)
1332 if (.not. surface_converged) n_not_converged = n_not_converged + 1
1333 else
1334 n_ill_posed = n_ill_posed + 1
1335 end if
1336
1337 if (surface_converged) then
1338 call get_mixture_molecular_weight(ys_s, mw_s)
1339
1340 ! Intrinsic gas density at the reacting surface.
1341 rho_s = pres_ip*mw_s/(gas_constant*t_s)
1342 if (rho_s > 0._wp) v_stefan = mdot_s/rho_s
1343 end if
1344 end if
1345
1346 ! Inert surface, or a reacting surface whose solve failed: zero species flux, so the surface composition
1347 ! is the image-point composition. Only the chemistry is dropped on failure; a prescribed Twall is kept,
1348 ! and thermal_bc = 2, which has no wall temperature without the energy balance, takes a zero gradient.
1349 if (.not. surface_converged) then
1350 ys_s(1:num_species) = ys_ip(1:num_species)
1351 t_s = t_ip
1352 if (patch_ib(patch_id)%thermal_bc == 1) t_s = patch_ib(patch_id)%Twall
1353 end if
1354
1355 call s_blend_ghost_temperature(t_ip, t_s, t_surface_min, t_surface_max, t_g)
1356 call s_blend_ghost_species(ys_ip, ys_s, ys_g)
1357
1358 call get_mixture_molecular_weight(ys_g, mw_g)
1359 alpha_rho_ip(1) = alpha_ip(1)*pres_ip*mw_g/(gas_constant*t_g)
1360 end if
1361
1362 ! Isothermal wall without chemistry: keep p_IP, set the ghost density from the mirrored temperature.
1363 if (.not. chemistry .and. patch_ib(patch_id)%thermal_bc == 1) then
1364 call s_phase_temperature(alpha_rho_ip(1)/alpha_ip(1), pres_ip, 1, t_ip)
1365 call s_blend_ghost_temperature(t_ip, patch_ib(patch_id)%Twall, 0._wp, huge(1._wp), t_g)
1366 call s_phase_density(pres_ip, t_g, 1, rho_g)
1367 alpha_rho_ip(1) = alpha_ip(1)*rho_g
1368 end if
1369
1370 ! If in simulation, use acc mixture subroutines
1371 if (hypoelasticity) then
1372 call s_convert_species_to_mixture_variables_kernel(rho, gamma, pi_inf, qv_k, alpha_ip, alpha_rho_ip, re_k, &
1373 & g_k, gs)
1374 else
1375 call s_convert_species_to_mixture_variables_kernel(rho, gamma, pi_inf, qv_k, alpha_ip, alpha_rho_ip, re_k)
1376 end if
1377
1378 if (surface_tension) q_prim_vf(eqn_idx%c)%sf(j, k, l) = c_ip
1379
1380 ! set the pressure and density
1381 if (patch_ib(patch_id)%moving_ibm <= 1) then
1382 pres_gp = pres_ip
1383 alpha_rho_gp = alpha_rho_ip
1384 else
1385 call s_compute_ghost_point_pressure(gp, patch_id, rho, pres_ip, pres_gp)
1386
1387 ! The adiabatic wall condition T_GP = T_IP the correction is derived from also
1388 ! fixes the ghost density: p + B = (n - 1)*cv*rho*T at both points under the one
1389 ! temperature leaves each partial density carrying the pressure ratio. The volume
1390 ! fractions are untouched, so only rho and qv move with it.
1391
1392# 477 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1393#if defined(MFC_OpenACC)
1394# 477 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1395!$acc loop seq
1396# 477 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1397#elif defined(MFC_OpenMP)
1398# 477 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1399
1400# 477 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1401#endif
1402# 478 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1403 do q = 1, num_fluids
1404 alpha_rho_gp(q) = alpha_rho_ip(q)*(pres_gp + isentrope_b(q))/(pres_ip + isentrope_b(q))
1405 end do
1406 call s_compute_mixture_coefficients(alpha_rho_gp, alpha_ip, rho, gamma, pi_inf, qv_k)
1407 end if
1408 q_prim_vf(eqn_idx%E)%sf(j, k, l) = pres_gp
1409
1410 ! Set q_prim_vf params at GP
1411
1412# 486 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1413#if defined(MFC_OpenACC)
1414# 486 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1415!$acc loop seq
1416# 486 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1417#elif defined(MFC_OpenMP)
1418# 486 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1419
1420# 486 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1421#endif
1422# 487 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1423 do q = 1, num_fluids
1424 q_prim_vf(q)%sf(j, k, l) = alpha_rho_gp(q)
1425 q_prim_vf(eqn_idx%adv%beg + q - 1)%sf(j, k, l) = alpha_ip(q)
1426 end do
1427
1428 ! get the vector that points from the centroid to the ghost
1429 radial_vector(1) = physical_loc(1) - (patch_ib(patch_id)%x_centroid + real(ghost_points(i)%x_periodicity, &
1430 & wp)*(glb_bounds(1)%end - glb_bounds(1)%beg))
1431 radial_vector(2) = physical_loc(2) - (patch_ib(patch_id)%y_centroid + real(ghost_points(i)%y_periodicity, &
1432 & wp)*(glb_bounds(2)%end - glb_bounds(2)%beg))
1433 radial_vector(3) = 0._wp
1434 if (num_dims == 3) radial_vector(3) = physical_loc(3) - (patch_ib(patch_id)%z_centroid &
1435 & + real(ghost_points(i)%z_periodicity, wp)*(glb_bounds(3)%end - glb_bounds(3)%beg))
1436
1437 ! Calculate velocity of ghost cell
1438 call s_compute_ghost_point_velocity(gp, patch_id, radial_vector, vel_ip, pres_ip, vel_g)
1439
1440 if (chemistry .and. patch_ib(patch_id)%inj_species == 0 .and. patch_ib(patch_id)%surface_reaction == 1 &
1441 & .and. surface_converged) then
1442 norm(1:3) = gp%levelset_norm
1443 buf = sqrt(sum(norm**2))
1444 if (buf > 0._wp) vel_g = vel_g + v_stefan*norm/buf
1445 end if
1446
1447 ! Set momentum
1448 vel_sum_g = 0._wp
1449
1450# 513 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1451#if defined(MFC_OpenACC)
1452# 513 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1453!$acc loop seq
1454# 513 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1455#elif defined(MFC_OpenMP)
1456# 513 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1457
1458# 513 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1459#endif
1460# 514 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1461 do q = eqn_idx%mom%beg, eqn_idx%mom%end
1462 q_cons_vf(q)%sf(j, k, l) = rho*vel_g(q - eqn_idx%mom%beg + 1)
1463 q_prim_vf(q)%sf(j, k, l) = vel_g(q - eqn_idx%mom%beg + 1)
1464 vel_sum_g = vel_sum_g + vel_g(q - eqn_idx%mom%beg + 1)**2._wp
1465 end do
1466 dyn_pres = 5.e-1_wp*rho*vel_sum_g
1467
1468 ! Set continuity and adv vars
1469
1470# 522 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1471#if defined(MFC_OpenACC)
1472# 522 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1473!$acc loop seq
1474# 522 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1475#elif defined(MFC_OpenMP)
1476# 522 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1477
1478# 522 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1479#endif
1480# 523 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1481 do q = 1, num_fluids
1482 q_cons_vf(q)%sf(j, k, l) = alpha_rho_gp(q)
1483 q_cons_vf(eqn_idx%adv%beg + q - 1)%sf(j, k, l) = alpha_ip(q)
1484 end do
1485
1486 ! Set color function
1487 if (surface_tension) q_cons_vf(eqn_idx%c)%sf(j, k, l) = c_ip
1488
1489 ! Set Energy
1490 if (chemistry) then
1491 ! Use the reconstructed thermal/species ghost state for an inert
1492 ! thermal surface or a converged heterogeneous reacting surface.
1493 if (patch_ib(patch_id)%inj_species == 0 .and. (patch_ib(patch_id)%surface_reaction == 0 &
1494 & .or. surface_converged)) then
1495
1496 call get_mixture_energy_mass(t_g, ys_g, e_ip)
1497
1498# 539 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1499#if defined(MFC_OpenACC)
1500# 539 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1501!$acc loop seq
1502# 539 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1503#elif defined(MFC_OpenMP)
1504# 539 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1505
1506# 539 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1507#endif
1508# 540 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1509 do q = 1, num_species
1510 q_cons_vf(eqn_idx%species%beg + q - 1)%sf(j, k, l) = rho*ys_g(q)
1511 end do
1512 else
1513 ! Ordinary chemistry/injection, or fallback after a failed
1514 ! heterogeneous surface solve: retain the image-point state.
1515 mw_ip = 0._wp
1516 call get_mixture_molecular_weight(ys_ip, mw_ip)
1517 t_ip = pres_ip*mw_ip/(rho*gas_constant)
1518 call get_mixture_energy_mass(t_ip, ys_ip, e_ip)
1519
1520# 550 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1521#if defined(MFC_OpenACC)
1522# 550 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1523!$acc loop seq
1524# 550 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1525#elif defined(MFC_OpenMP)
1526# 550 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1527
1528# 550 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1529#endif
1530# 551 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1531 do q = 1, num_species
1532 q_cons_vf(eqn_idx%species%beg + q - 1)%sf(j, k, l) = rho*ys_ip(q)
1533 end do
1534 end if
1535
1536 q_cons_vf(eqn_idx%E)%sf(j, k, l) = rho*e_ip + dyn_pres
1537 else
1538 call s_compute_energy(pres_gp, alpha_rho_gp, alpha_ip, vel_sum_g, e_ghost)
1539 q_cons_vf(eqn_idx%E)%sf(j, k, l) = e_ghost
1540 end if
1541 ! Set bubble vars
1542 if (bubbles_euler .and. .not. qbmm) then
1543 call s_comp_n_from_prim(alpha_ip(1), r_ip, nbub, weight)
1544
1545# 564 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1546#if defined(MFC_OpenACC)
1547# 564 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1548!$acc loop seq
1549# 564 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1550#elif defined(MFC_OpenMP)
1551# 564 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1552
1553# 564 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1554#endif
1555# 565 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1556 do q = 1, nb
1557 q_cons_vf(eqn_idx%bub%beg + (q - 1)*2)%sf(j, k, l) = nbub*r_ip(q)
1558 q_cons_vf(eqn_idx%bub%beg + (q - 1)*2 + 1)%sf(j, k, l) = nbub*v_ip(q)
1559 if (.not. polytropic) then
1560 q_cons_vf(eqn_idx%bub%beg + (q - 1)*4)%sf(j, k, l) = nbub*r_ip(q)
1561 q_cons_vf(eqn_idx%bub%beg + (q - 1)*4 + 1)%sf(j, k, l) = nbub*v_ip(q)
1562 q_cons_vf(eqn_idx%bub%beg + (q - 1)*4 + 2)%sf(j, k, l) = nbub*pb_ip(q)
1563 q_cons_vf(eqn_idx%bub%beg + (q - 1)*4 + 3)%sf(j, k, l) = nbub*mv_ip(q)
1564 end if
1565 end do
1566 end if
1567
1568 if (qbmm) then
1569 nbub = nmom_ip(1)
1570
1571# 579 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1572#if defined(MFC_OpenACC)
1573# 579 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1574!$acc loop seq
1575# 579 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1576#elif defined(MFC_OpenMP)
1577# 579 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1578
1579# 579 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1580#endif
1581# 580 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1582 do q = 1, nb*nmom
1583 q_cons_vf(eqn_idx%bub%beg + q - 1)%sf(j, k, l) = nbub*nmom_ip(q)
1584 end do
1585
1586
1587# 584 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1588#if defined(MFC_OpenACC)
1589# 584 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1590!$acc loop seq
1591# 584 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1592#elif defined(MFC_OpenMP)
1593# 584 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1594
1595# 584 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1596#endif
1597# 585 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1598 do q = 1, nb
1599 q_cons_vf(eqn_idx%bub%beg + (q - 1)*nmom)%sf(j, k, l) = nbub
1600 end do
1601
1602 if (.not. polytropic) then
1603
1604# 590 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1605#if defined(MFC_OpenACC)
1606# 590 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1607!$acc loop seq
1608# 590 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1609#elif defined(MFC_OpenMP)
1610# 590 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1611
1612# 590 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1613#endif
1614# 591 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1615 do q = 1, nb
1616
1617# 592 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1618#if defined(MFC_OpenACC)
1619# 592 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1620!$acc loop seq
1621# 592 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1622#elif defined(MFC_OpenMP)
1623# 592 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1624
1625# 592 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1626#endif
1627# 593 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1628 do r = 1, nnode
1629 pb_in(j, k, l, r, q) = presb_ip((q - 1)*nnode + r)
1630 mv_in(j, k, l, r, q) = massv_ip((q - 1)*nnode + r)
1631 end do
1632 end do
1633 end if
1634 end if
1635
1636 if (model_eqns == model_eqns_6eq) then
1637
1638# 602 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1639#if defined(MFC_OpenACC)
1640# 602 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1641!$acc loop seq
1642# 602 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1643#elif defined(MFC_OpenMP)
1644# 602 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1645
1646# 602 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1647#endif
1648# 603 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1649 do q = eqn_idx%int_en%beg, eqn_idx%int_en%end
1650 alpha_q = alpha_ip(q - eqn_idx%int_en%beg + 1)
1651 alpha_rho_q = alpha_rho_gp(q - eqn_idx%int_en%beg + 1)
1652 call s_phase_internal_energy(pres_gp, alpha_q, alpha_rho_q, q - eqn_idx%int_en%beg + 1, e_q)
1653 q_cons_vf(q)%sf(j, k, l) = e_q
1654 end do
1655 end if
1656 end do
1657
1658# 611 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1659#if defined(MFC_OpenACC)
1660# 611 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1661!$acc end parallel loop
1662# 611 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1663#elif defined(MFC_OpenMP)
1664# 611 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1665
1666# 611 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1667!$omp end target teams loop
1668# 611 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1669#endif
1670# 612 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1671
1672 ! A ghost point whose image point is buried in a neighboring IB has no fluid to mirror, so
1673 ! it takes the average of the ghost points corrected above, growing the stencil until it
1674 ! reaches one. Those neighbors already carry the wall condition, so the average does too.
1675
1676# 616 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1677
1678# 616 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1679#if defined(MFC_OpenACC)
1680# 616 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1681!$acc parallel loop gang vector default(present) private(i, j, k, l, q, r, jj, kk, ll, gp, rad, rad_z, num_nbrs, buf, buf_prim) present(ghost_points)
1682# 616 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1683#elif defined(MFC_OpenMP)
1684# 616 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1685
1686# 616 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1687
1688# 616 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1689
1690# 616 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1691!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) &
1692# 616 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1693!$omp& private(i, j, k, l, q, r, jj, kk, ll, gp, rad, rad_z, num_nbrs, buf, buf_prim)
1694# 616 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1695#endif
1696# 618 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1697 do i = 1, num_gps
1698 gp = ghost_points(i)
1699 if (gp%interp_valid) cycle
1700 j = gp%loc(1)
1701 k = gp%loc(2)
1702 l = gp%loc(3)
1703
1704 num_nbrs = 0
1705 rad = 0
1706 rad_z = 0
1707 do while (num_nbrs == 0 .and. rad < gp_layers)
1708 rad = rad + 1
1709 if (p /= 0) rad_z = rad
1710 do jj = j - rad, j + rad
1711 do kk = k - rad, k + rad
1712 do ll = l - rad_z, l + rad_z
1713 num_nbrs = num_nbrs + corrected_gps%sf(jj, kk, ll)
1714 end do
1715 end do
1716 end do
1717 end do
1718
1719 ! Fully enclosed: nothing to average from, so leave the interior placeholder in place
1720 if (num_nbrs == 0) cycle
1721
1722
1723# 643 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1724#if defined(MFC_OpenACC)
1725# 643 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1726!$acc loop seq
1727# 643 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1728#elif defined(MFC_OpenMP)
1729# 643 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1730
1731# 643 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1732#endif
1733# 644 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1734 do q = 1, sys_size
1735 buf = 0._wp
1736 buf_prim = 0._wp
1737 do jj = j - rad, j + rad
1738 do kk = k - rad, k + rad
1739 do ll = l - rad_z, l + rad_z
1740 if (corrected_gps%sf(jj, kk, ll) == 1) then
1741 buf = buf + q_cons_vf(q)%sf(jj, kk, ll)
1742 buf_prim = buf_prim + q_prim_vf(q)%sf(jj, kk, ll)
1743 end if
1744 end do
1745 end do
1746 end do
1747 q_cons_vf(q)%sf(j, k, l) = buf/real(num_nbrs, wp)
1748 q_prim_vf(q)%sf(j, k, l) = buf_prim/real(num_nbrs, wp)
1749 end do
1750
1751 if (qbmm .and. .not. polytropic) then
1752
1753# 662 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1754#if defined(MFC_OpenACC)
1755# 662 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1756!$acc loop seq
1757# 662 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1758#elif defined(MFC_OpenMP)
1759# 662 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1760
1761# 662 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1762#endif
1763# 663 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1764 do q = 1, nb
1765
1766# 664 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1767#if defined(MFC_OpenACC)
1768# 664 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1769!$acc loop seq
1770# 664 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1771#elif defined(MFC_OpenMP)
1772# 664 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1773
1774# 664 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1775#endif
1776# 665 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1777 do r = 1, nnode
1778 buf = 0._wp
1779 buf_prim = 0._wp
1780 do jj = j - rad, j + rad
1781 do kk = k - rad, k + rad
1782 do ll = l - rad_z, l + rad_z
1783 if (corrected_gps%sf(jj, kk, ll) == 1) then
1784 buf = buf + pb_in(jj, kk, ll, r, q)
1785 buf_prim = buf_prim + mv_in(jj, kk, ll, r, q)
1786 end if
1787 end do
1788 end do
1789 end do
1790 pb_in(j, k, l, r, q) = buf/real(num_nbrs, wp)
1791 mv_in(j, k, l, r, q) = buf_prim/real(num_nbrs, wp)
1792 end do
1793 end do
1794 end if
1795 end do
1796
1797# 684 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1798#if defined(MFC_OpenACC)
1799# 684 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1800!$acc end parallel loop
1801# 684 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1802#elif defined(MFC_OpenMP)
1803# 684 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1804
1805# 684 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1806!$omp end target teams loop
1807# 684 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1808#endif
1809# 685 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1810 end if
1811
1812 n_surface_not_converged = n_surface_not_converged + n_not_converged
1813 n_surface_ill_posed = n_surface_ill_posed + n_ill_posed
1814
1815 end subroutine s_ibm_correct_state
1816
1817 !> One line at the end of a run if the reacting surface ever failed to solve. Silence means every ghost point on every step
1818 !! reached the surface-balance tolerance, i.e. the surface chemistry was actually applied everywhere it was asked for.
1819 impure subroutine s_report_ibm_surface
1820
1821 integer :: n_not_converged_glb, n_ill_posed_glb
1822
1823 call s_mpi_reduce_int_sum(n_surface_not_converged, n_not_converged_glb)
1824 call s_mpi_reduce_int_sum(n_surface_ill_posed, n_ill_posed_glb)
1825
1826 if (proc_rank /= 0) return
1827
1828 if (n_not_converged_glb > 0) then
1829 print '(A,I0,A)', ' Immersed-boundary surface chemistry: the Newton solve did not converge at ', n_not_converged_glb, &
1830 & ' ghost-point updates, which fell back to a chemically inert wall (Twall kept for thermal_bc = 1, ' &
1831 & // 'zero temperature gradient for thermal_bc = 2).'
1832 end if
1833
1834 if (n_ill_posed_glb > 0) then
1835 print '(A,I0,A)', ' Immersed-boundary surface chemistry: ', n_ill_posed_glb, &
1836 & ' ghost-point updates had a zero levelset distance, leaving the surface balance undefined; those points ' &
1837 & // 'fell back to a chemically inert wall (Twall kept for thermal_bc = 1, zero temperature gradient for ' &
1838 & // 'thermal_bc = 2). Check the immersed geometry against the grid.'
1839 end if
1840
1841 end subroutine s_report_ibm_surface
1842
1843 !> Compute the image points for each ghost point
1844 impure subroutine s_compute_image_points()
1845
1846 real(wp) :: dist
1847 real(wp) :: min_cell_width
1848 real(wp), dimension(3) :: norm
1849 real(wp), dimension(3) :: physical_loc
1850 real(wp) :: temp_loc
1851 real(wp), pointer, dimension(:) :: s_cc => null()
1852 integer :: bound
1853 type(ghost_point) :: gp
1854 integer :: q, dim !< Iterator variables
1855 integer :: i, j, k, l !< Location indexes
1856 integer :: patch_id !< IB Patch ID
1857 integer :: dir
1858 integer :: index
1859 logical :: bounds_error
1860
1861 bounds_error = .false.
1862
1863
1864# 738 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1865
1866# 738 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1867#if defined(MFC_OpenACC)
1868# 738 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1869!$acc parallel loop gang vector default(present) private(q, gp, i, j, k, physical_loc, patch_id, dist, min_cell_width, norm, dim, bound, dir, index, temp_loc, s_cc) copy(bounds_error) &
1870# 738 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1871!$acc& present(ghost_points)
1872# 738 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1873#elif defined(MFC_OpenMP)
1874# 738 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1875
1876# 738 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1877
1878# 738 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1879
1880# 738 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1881!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) &
1882# 738 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1883!$omp& private(q, gp, i, j, k, physical_loc, patch_id, dist, min_cell_width, norm, dim, bound, dir, index, temp_loc, s_cc) map(tofrom:bounds_error)
1884# 738 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1885#endif
1886# 740 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1887 do q = 1, num_gps
1888 gp = ghost_points(q)
1889 i = gp%loc(1)
1890 j = gp%loc(2)
1891 k = gp%loc(3)
1892
1893 ! Calculate physical location of ghost point
1894 if (p > 0) then
1895 physical_loc = [x_cc(i), y_cc(j), z_cc(k)]
1896 else
1897 physical_loc = [x_cc(i), y_cc(j), 0._wp]
1898 end if
1899
1900 ! Calculate and store the precise location of the image point
1901 patch_id = gp%ib_patch_id
1902 dist = abs(real(gp%levelset, kind=wp))
1903 norm(:) = gp%levelset_norm
1904 ghost_points(q)%ip_dist = dist
1905 ! keep the image point stencil off the wall
1906 if (ib_second_order_vel) then
1907 min_cell_width = min(dx(i), dy(j))
1908 if (p > 0) min_cell_width = min(min_cell_width, dz(k))
1909 ghost_points(q)%ip_dist = max(dist, ib_ip_min_dist*min_cell_width)
1910 end if
1911 ghost_points(q)%ip_loc(:) = physical_loc(:) + (dist + ghost_points(q)%ip_dist)*norm(:)
1912
1913 ! Find the closest grid point to the image point
1914 do dim = 1, num_dims
1915 ! s_cc points to the dim array we need
1916 if (dim == 1) then
1917 s_cc => x_cc
1918 bound = m + buff_size - 1
1919 else if (dim == 2) then
1920 s_cc => y_cc
1921 bound = n + buff_size - 1
1922 else
1923 s_cc => z_cc
1924 bound = p + buff_size - 1
1925 end if
1926
1927 if (f_approx_equal(norm(dim), 0._wp)) then
1928 ! if the ghost point is almost equal to a cell location, we set it equal and continue
1929 ghost_points(q)%ip_grid(dim) = ghost_points(q)%loc(dim)
1930 else
1931 if (norm(dim) > 0) then
1932 dir = 1
1933 else
1934 dir = -1
1935 end if
1936
1937 index = ghost_points(q)%loc(dim)
1938 temp_loc = ghost_points(q)%ip_loc(dim)
1939 do while ((temp_loc < s_cc(index) .or. temp_loc > s_cc(index + 1)) .and. (.not. bounds_error))
1940 index = index + dir
1941 if (index < -buff_size .or. index > bound) then
1942#if !defined(MFC_OpenACC) && !defined(MFC_OpenMP)
1943 print *, "A required image point is not located in this computational domain."
1944 print *, "Ghost Point is located at :"
1945 if (p == 0) then
1946 print *, [x_cc(i), y_cc(j)]
1947 else
1948 print *, [x_cc(i), y_cc(j), z_cc(k)]
1949 end if
1950 print *, "We are searching in dimension ", dim, " for image point at ", ghost_points(q)%ip_loc(:)
1951 print *, "Domain size: "
1952 print *, "x: ", x_cc(-buff_size), " to: ", x_cc(m + buff_size - 1)
1953 print *, "y: ", y_cc(-buff_size), " to: ", y_cc(n + buff_size - 1)
1954 if (p /= 0) print *, "z: ", z_cc(-buff_size), " to: ", z_cc(p + buff_size - 1)
1955 print *, "Image point is located approximately ", &
1956 & (ghost_points(q)%loc(dim) - ghost_points(q) %ip_loc(dim))/(s_cc(1) - s_cc(0)), " grid cells away"
1957 print *, "Levelset ", dist, " and Norm: ", norm(:)
1958 print *, &
1959 & "A short term fix may include increasing buff_size further in m_helper_basic (currently set to a minimum of 10)"
1960#endif
1961 bounds_error = .true.
1962 end if
1963 end do
1964
1965 ghost_points(q)%ip_grid(dim) = index
1966 if (ghost_points(q)%DB(dim) == -1) then
1967 ghost_points(q)%ip_grid(dim) = ghost_points(q)%loc(dim) + 1
1968 else if (ghost_points(q)%DB(dim) == 1) then
1969 ghost_points(q)%ip_grid(dim) = ghost_points(q)%loc(dim) - 1
1970 end if
1971 end if
1972 end do
1973 end do
1974
1975# 827 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1976#if defined(MFC_OpenACC)
1977# 827 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1978!$acc end parallel loop
1979# 827 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1980#elif defined(MFC_OpenMP)
1981# 827 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1982
1983# 827 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1984!$omp end target teams loop
1985# 827 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1986#endif
1987# 828 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1988
1989 if (bounds_error) then
1990# 829 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1991 call s_prohibit_abort("bounds_error", "Ghost Point and Image Point on Different Processors. Exiting")
1992# 829 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1993 end if
1994# 830 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
1995
1996 end subroutine s_compute_image_points
1997
1998 !> Count the number of ghost points for memory allocation
1999 subroutine s_find_num_ghost_points(num_gps_out)
2000
2001 integer, intent(out) :: num_gps_out
2002 integer :: i, j, k, ii, jj, kk, gp_layers_z !< Iterator variables
2003 integer :: num_gps_local !< local copies of the gp count to support GPU compute
2004 logical :: is_gp
2005
2006 num_gps_local = 0
2007 gp_layers_z = gp_layers
2008 if (p == 0) gp_layers_z = 0
2009
2010
2011# 845 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2012
2013# 845 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2014#if defined(MFC_OpenACC)
2015# 845 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2016!$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)
2017# 845 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2018#elif defined(MFC_OpenMP)
2019# 845 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2020
2021# 845 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2022
2023# 845 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2024
2025# 845 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2026!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) &
2027# 845 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2028!$omp& private(i, j, k, ii, jj, kk, is_gp) firstprivate(gp_layers, gp_layers_z) map(tofrom:num_gps_local)
2029# 845 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2030#endif
2031# 847 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2032 do i = 0, m
2033 do j = 0, n
2034 do k = 0, p
2035 if (ib_markers%sf(i, j, k) /= 0) then
2036 is_gp = .false.
2037 marker_search: do ii = i - gp_layers, i + gp_layers
2038 do jj = j - gp_layers, j + gp_layers
2039 do kk = k - gp_layers_z, k + gp_layers_z
2040 if (ib_markers%sf(ii, jj, kk) == 0) then
2041 ! if any neighbors are not in the IB, it is a ghost point
2042 is_gp = .true.
2043 exit marker_search
2044 end if
2045 end do
2046 end do
2047 end do marker_search
2048
2049 if (is_gp) then
2050
2051# 865 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2052#if defined(MFC_OpenACC)
2053# 865 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2054!$acc atomic update
2055# 865 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2056#elif defined(MFC_OpenMP)
2057# 865 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2058!$omp atomic update
2059# 865 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2060#endif
2061# 866 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2062 num_gps_local = num_gps_local + 1
2063 end if
2064 end if
2065 end do
2066 end do
2067 end do
2068
2069# 872 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2070#if defined(MFC_OpenACC)
2071# 872 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2072!$acc end parallel loop
2073# 872 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2074#elif defined(MFC_OpenMP)
2075# 872 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2076
2077# 872 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2078!$omp end target teams loop
2079# 872 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2080#endif
2081# 873 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2082
2083 num_gps_out = num_gps_local
2084
2085 end subroutine s_find_num_ghost_points
2086
2087 !> Locate all ghost points in the domain
2088 subroutine s_find_ghost_points()
2089
2090 integer :: i, j, k, ii, jj, kk, gp_layers_z !< Iterator variables
2091 integer :: xp, yp, zp !< periodicities
2092 integer :: count, count_i, local_idx
2093 integer :: patch_id, encoded_patch_id, neighborhood_patch_id
2094 logical :: is_gp
2095
2096 count = 0
2097 count_i = 0
2098 gp_layers_z = gp_layers
2099 if (p == 0) gp_layers_z = 0
2100
2101
2102# 892 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2103
2104# 892 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2105#if defined(MFC_OpenACC)
2106# 892 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2107!$acc parallel loop collapse(3) gang vector default(present) private(i, j, k, ii, jj, kk, is_gp, local_idx, patch_id, encoded_patch_id, neighborhood_patch_id, xp, yp, zp) &
2108# 892 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2109!$acc& firstprivate(gp_layers, gp_layers_z) copyin(count, count_i, glb_bounds) present(ghost_points)
2110# 892 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2111#elif defined(MFC_OpenMP)
2112# 892 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2113
2114# 892 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2115
2116# 892 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2117
2118# 892 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2119!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) &
2120# 892 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2121!$omp& private(i, j, k, ii, jj, kk, is_gp, local_idx, patch_id, encoded_patch_id, neighborhood_patch_id, xp, yp, zp) firstprivate(gp_layers, gp_layers_z) map(to:count, count_i, glb_bounds)
2122# 892 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2123#endif
2124# 895 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2125 do i = 0, m
2126 do j = 0, n
2127 do k = 0, p
2128 if (ib_markers%sf(i, j, k) /= 0) then
2129 is_gp = .false.
2130 marker_search: do ii = i - gp_layers, i + gp_layers
2131 do jj = j - gp_layers, j + gp_layers
2132 do kk = k - gp_layers_z, k + gp_layers_z
2133 if (ib_markers%sf(ii, jj, kk) == 0) then
2134 ! if any neighbors are not in the IB, it is a ghost point
2135 is_gp = .true.
2136 exit marker_search
2137 end if
2138 end do
2139 end do
2140 end do marker_search
2141
2142 if (is_gp) then
2143
2144# 913 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2145#if defined(MFC_OpenACC)
2146# 913 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2147!$acc atomic capture
2148# 913 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2149#elif defined(MFC_OpenMP)
2150# 913 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2151!$omp atomic capture
2152# 913 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2153#endif
2154# 914 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2155 count = count + 1
2156 local_idx = count
2157#if defined(MFC_OpenACC)
2158# 916 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2159!$acc end atomic
2160# 916 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2161#elif defined(MFC_OpenMP)
2162# 916 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2163!$omp end atomic
2164# 916 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2165#endif
2166# 917 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2167
2168 ghost_points(local_idx)%loc = [i, j, k]
2169 encoded_patch_id = ib_markers%sf(i, j, k)
2170 call s_decode_patch_periodicity(encoded_patch_id, patch_id, xp, yp, zp)
2171 call s_get_neighborhood_idx(patch_id, neighborhood_patch_id)
2172 ghost_points(local_idx)%ib_patch_id = neighborhood_patch_id
2173 ghost_points(local_idx)%x_periodicity = xp
2174 ghost_points(local_idx)%y_periodicity = yp
2175 ghost_points(local_idx)%z_periodicity = zp
2176 ghost_points(local_idx)%slip = patch_ib(neighborhood_patch_id)%slip
2177
2178# 929 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2179 ghost_points(local_idx)%DB(1) = 0
2180# 930 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2181 if (1 <= num_dims) then ! Onlyrun the numeher of dimensions present
2182# 931 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2183 if (ib_bc_x%beg /= bc_periodic) then
2184# 932 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2185 if ((x_cc(i) - dx(i)) < glb_bounds(1)%beg) then
2186# 933 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2187 ! if the grid cell is in a wall on the "left"
2188 ghost_points(local_idx)%DB(1) = -1
2189# 935 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2190 else if ((x_cc(i) + dx(i)) > glb_bounds(1)%end) then
2191# 936 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2192 ! if the grid cell is in a wall on the "right"
2193 ghost_points(local_idx)%DB(1) = 1
2194# 938 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2195 end if
2196 end if
2197 end if
2198# 929 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2199 ghost_points(local_idx)%DB(2) = 0
2200# 930 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2201 if (2 <= num_dims) then ! Onlyrun the numeher of dimensions present
2202# 931 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2203 if (ib_bc_y%beg /= bc_periodic) then
2204# 932 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2205 if ((y_cc(j) - dy(j)) < glb_bounds(2)%beg) then
2206# 933 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2207 ! if the grid cell is in a wall on the "left"
2208 ghost_points(local_idx)%DB(2) = -1
2209# 935 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2210 else if ((y_cc(j) + dy(j)) > glb_bounds(2)%end) then
2211# 936 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2212 ! if the grid cell is in a wall on the "right"
2213 ghost_points(local_idx)%DB(2) = 1
2214# 938 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2215 end if
2216 end if
2217 end if
2218# 929 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2219 ghost_points(local_idx)%DB(3) = 0
2220# 930 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2221 if (3 <= num_dims) then ! Onlyrun the numeher of dimensions present
2222# 931 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2223 if (ib_bc_z%beg /= bc_periodic) then
2224# 932 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2225 if ((z_cc(k) - dz(k)) < glb_bounds(3)%beg) then
2226# 933 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2227 ! if the grid cell is in a wall on the "left"
2228 ghost_points(local_idx)%DB(3) = -1
2229# 935 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2230 else if ((z_cc(k) + dz(k)) > glb_bounds(3)%end) then
2231# 936 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2232 ! if the grid cell is in a wall on the "right"
2233 ghost_points(local_idx)%DB(3) = 1
2234# 938 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2235 end if
2236 end if
2237 end if
2238# 942 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2239 end if
2240 end if
2241 end do
2242 end do
2243 end do
2244
2245# 947 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2246#if defined(MFC_OpenACC)
2247# 947 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2248!$acc end parallel loop
2249# 947 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2250#elif defined(MFC_OpenMP)
2251# 947 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2252
2253# 947 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2254!$omp end target teams loop
2255# 947 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2256#endif
2257# 948 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2258
2259 end subroutine s_find_ghost_points
2260
2261 !> Compute the interpolation coefficients for image points
2262 subroutine s_compute_interpolation_coeffs()
2263
2264 real(wp), dimension(2, 2, 2) :: dist
2265 real(wp), dimension(2, 2, 2) :: fluid_mask
2266 real(wp), dimension(2, 2, 2) :: interp_coeffs
2267 type(ghost_point) :: gp
2268 integer :: q, i, j, k, ii, jj, kk !< Grid indexes and iterators
2269
2270
2271# 960 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2272
2273# 960 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2274#if defined(MFC_OpenACC)
2275# 960 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2276!$acc parallel loop gang vector default(present) private(q, i, j, k, ii, jj, kk, dist, gp, interp_coeffs, fluid_mask) present(ghost_points)
2277# 960 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2278#elif defined(MFC_OpenMP)
2279# 960 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2280
2281# 960 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2282
2283# 960 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2284
2285# 960 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2286!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) &
2287# 960 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2288!$omp& private(q, i, j, k, ii, jj, kk, dist, gp, interp_coeffs, fluid_mask)
2289# 960 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2290#endif
2291# 961 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2292 do q = 1, num_gps
2293 gp = ghost_points(q)
2294 ! Get the interpolation points
2295 i = gp%ip_grid(1)
2296 j = gp%ip_grid(2)
2297 if (p /= 0) then
2298 k = gp%ip_grid(3)
2299 else
2300 k = 0
2301 end if
2302
2303 ! get the distance to a cell in each direction
2304 dist = 0._wp
2305 do ii = 0, 1
2306 do jj = 0, 1
2307 if (p == 0) then
2308 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)
2309 else
2310 do kk = 0, 1
2311 dist(1 + ii, 1 + jj, &
2312 & 1 + kk) = sqrt((x_cc(i + ii) - gp%ip_loc(1))**2 + (y_cc(j + jj) - gp%ip_loc(2))**2 + (z_cc(k &
2313 & + kk) - gp%ip_loc(3))**2)
2314 end do
2315 end if
2316 end do
2317 end do
2318
2319 ! flag the fluid cells of the stencil
2320 fluid_mask = 0._wp
2321 do kk = 0, min(p, 1)
2322 do jj = 0, 1
2323 do ii = 0, 1
2324 if (ib_markers%sf(i + ii, j + jj, k + kk) == 0) fluid_mask(ii + 1, jj + 1, kk + 1) = 1._wp
2325 end do
2326 end do
2327 end do
2328
2329 call s_interpolation_weights(gp, i, j, k, dist, fluid_mask, interp_coeffs)
2330
2331 ! An image point buried in a neighboring IB gets no weights at all: its ghost point is
2332 ! averaged from its own neighbors instead, in s_ibm_correct_state.
2333 ghost_points(q)%interp_coeffs = interp_coeffs
2334 ghost_points(q)%interp_valid = any(interp_coeffs > 0._wp)
2335 if (ghost_points(q)%interp_valid) corrected_gps%sf(gp%loc(1), gp%loc(2), gp%loc(3)) = 1
2336 end do
2337
2338# 1006 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2339#if defined(MFC_OpenACC)
2340# 1006 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2341!$acc end parallel loop
2342# 1006 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2343#elif defined(MFC_OpenMP)
2344# 1006 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2345
2346# 1006 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2347!$omp end target teams loop
2348# 1006 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2349#endif
2350# 1007 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2351
2352 end subroutine s_compute_interpolation_coeffs
2353
2354 !> Fluid-masked image-point weights: bilinear for second-order IBM, else inverse distance squared
2355 subroutine s_interpolation_weights(gp, i, j, k, dist, fluid_mask, interp_coeffs)
2356
2357
2358# 1013 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2359#if MFC_OpenACC
2360# 1013 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2361!$acc routine seq
2362# 1013 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2363#elif MFC_OpenMP
2364# 1013 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2365
2366# 1013 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2367
2368# 1013 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2369!$omp declare target device_type(any)
2370# 1013 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2371#endif
2372# 1014 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2373
2374 type(ghost_point), intent(in) :: gp
2375 integer, intent(in) :: i, j, k !< lower corner of the stencil
2376 real(wp), dimension(2, 2, 2), intent(in) :: dist, fluid_mask
2377 real(wp), dimension(2, 2, 2), intent(out) :: interp_coeffs
2378 real(wp), dimension(2, 2, 2) :: weights
2379 real(wp), dimension(2, 3) :: axis_weight !< (stencil side, dimension) linear weights
2380 real(wp) :: weight_sum
2381 integer :: ii, jj, kk, num_z_layers
2382
2383 num_z_layers = min(p, 1) + 1
2384
2385 if (ib_second_order_vel) then
2386 ! bilinear weights, Mittal et al. (2008)
2387 axis_weight(2, 1) = (gp%ip_loc(1) - x_cc(i))/(x_cc(i + 1) - x_cc(i))
2388 axis_weight(2, 2) = (gp%ip_loc(2) - y_cc(j))/(y_cc(j + 1) - y_cc(j))
2389 axis_weight(2, 3) = 0._wp
2390 if (p /= 0) axis_weight(2, 3) = (gp%ip_loc(3) - z_cc(k))/(z_cc(k + 1) - z_cc(k))
2391 axis_weight(2,:) = min(max(axis_weight(2,:), 0._wp), 1._wp) ! clamp to the cell edges
2392 axis_weight(1,:) = 1._wp - axis_weight(2,:)
2393 do kk = 1, num_z_layers
2394 do jj = 1, 2
2395 do ii = 1, 2
2396 weights(ii, jj, kk) = axis_weight(ii, 1)*axis_weight(jj, 2)*axis_weight(kk, 3)
2397 end do
2398 end do
2399 end do
2400 else
2401 ! inverse distance squared, Chaudhuri et al. JCP (2011)
2402 if (minval(dist(:,:,1:num_z_layers)) <= 1.e-16_wp) then
2403 weights(:,:,1:num_z_layers) = merge(1._wp, 0._wp, dist(:,:,1:num_z_layers) <= 1.e-16_wp)
2404 else
2405 weights(:,:,1:num_z_layers) = 1._wp/dist(:,:,1:num_z_layers)**2
2406 end if
2407 end if
2408
2409 interp_coeffs = 0._wp
2410 weights(:,:,1:num_z_layers) = fluid_mask(:,:,1:num_z_layers)*weights(:,:,1:num_z_layers)
2411 weight_sum = sum(weights(:,:,1:num_z_layers))
2412 if (weight_sum > 0._wp) interp_coeffs(:,:,1:num_z_layers) = weights(:,:,1:num_z_layers)/weight_sum
2413
2414 end subroutine s_interpolation_weights
2415
2416 !> Interpolate primitive variables to a ghost point's image point using bilinear or trilinear interpolation
2417 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, &
2418 & nmom_IP, pb_in, mv_in, presb_IP, massv_IP, Ys_IP)
2419
2420# 1060 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2421#if MFC_OpenACC
2422# 1060 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2423!$acc routine seq
2424# 1060 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2425#elif MFC_OpenMP
2426# 1060 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2427
2428# 1060 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2429
2430# 1060 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2431!$omp declare target device_type(any)
2432# 1060 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2433#endif
2434# 1061 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2435
2436 type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf !< Primitive Variables
2437 real(stp), optional, dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:,1:), intent(in) :: pb_in, mv_in
2438 type(ghost_point), intent(in) :: gp
2439 real(wp), intent(inout) :: pres_ip
2440 real(wp), dimension(3), intent(inout) :: vel_ip
2441 real(wp), intent(inout) :: c_ip
2442# 1071 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2443 real(wp), dimension(num_fluids), intent(inout) :: alpha_ip, alpha_rho_ip
2444# 1073 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2445 real(wp), optional, dimension(:), intent(inout) :: r_ip, v_ip, pb_ip, mv_ip
2446 real(wp), optional, dimension(:), intent(inout) :: nmom_ip
2447 real(wp), optional, dimension(:), intent(inout) :: presb_ip, massv_ip
2448 real(wp), optional, dimension(:), intent(inout) :: ys_ip !< Interpolated species mass fractions (chemistry)
2449 integer :: i, j, k, l, q !< Iterator variables
2450 integer :: i1, i2, j1, j2, k1, k2 !< Iterator variables
2451 real(wp) :: coeff
2452
2453 i1 = gp%ip_grid(1); i2 = i1 + 1
2454 j1 = gp%ip_grid(2); j2 = j1 + 1
2455 k1 = gp%ip_grid(3); k2 = k1 + 1
2456
2457 if (p == 0) then
2458 k1 = 0
2459 k2 = 0
2460 end if
2461
2462 alpha_rho_ip = 0._wp
2463 alpha_ip = 0._wp
2464 pres_ip = 0._wp
2465 vel_ip = 0._wp
2466
2467 if (chemistry) ys_ip = 0._wp
2468
2469 if (surface_tension) c_ip = 0._wp
2470
2471 if (bubbles_euler) then
2472 r_ip = 0._wp
2473 v_ip = 0._wp
2474 if (.not. polytropic) then
2475 mv_ip = 0._wp
2476 pb_ip = 0._wp
2477 end if
2478 end if
2479
2480 if (qbmm) then
2481 nmom_ip = 0._wp
2482 if (.not. polytropic) then
2483 presb_ip = 0._wp
2484 massv_ip = 0._wp
2485 end if
2486 end if
2487
2488
2489# 1116 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2490#if defined(MFC_OpenACC)
2491# 1116 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2492!$acc loop seq
2493# 1116 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2494#elif defined(MFC_OpenMP)
2495# 1116 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2496
2497# 1116 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2498#endif
2499# 1117 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2500 do i = i1, i2
2501
2502# 1118 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2503#if defined(MFC_OpenACC)
2504# 1118 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2505!$acc loop seq
2506# 1118 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2507#elif defined(MFC_OpenMP)
2508# 1118 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2509
2510# 1118 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2511#endif
2512# 1119 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2513 do j = j1, j2
2514
2515# 1120 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2516#if defined(MFC_OpenACC)
2517# 1120 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2518!$acc loop seq
2519# 1120 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2520#elif defined(MFC_OpenMP)
2521# 1120 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2522
2523# 1120 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2524#endif
2525# 1121 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2526 do k = k1, k2
2527 coeff = gp%interp_coeffs(i - i1 + 1, j - j1 + 1, k - k1 + 1)
2528
2529 pres_ip = pres_ip + coeff*q_prim_vf(eqn_idx%E)%sf(i, j, k)
2530
2531
2532# 1126 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2533#if defined(MFC_OpenACC)
2534# 1126 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2535!$acc loop seq
2536# 1126 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2537#elif defined(MFC_OpenMP)
2538# 1126 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2539
2540# 1126 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2541#endif
2542# 1127 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2543 do q = eqn_idx%mom%beg, eqn_idx%mom%end
2544 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)
2545 end do
2546
2547
2548# 1131 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2549#if defined(MFC_OpenACC)
2550# 1131 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2551!$acc loop seq
2552# 1131 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2553#elif defined(MFC_OpenMP)
2554# 1131 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2555
2556# 1131 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2557#endif
2558# 1132 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2559 do l = eqn_idx%cont%beg, eqn_idx%cont%end
2560 alpha_rho_ip(l) = alpha_rho_ip(l) + coeff*q_prim_vf(l)%sf(i, j, k)
2561 alpha_ip(l) = alpha_ip(l) + coeff*q_prim_vf(eqn_idx%adv%beg + l - 1)%sf(i, j, k)
2562 end do
2563
2564 if (surface_tension) then
2565 c_ip = c_ip + coeff*q_prim_vf(eqn_idx%c)%sf(i, j, k)
2566 end if
2567
2568 if (chemistry) then
2569
2570# 1142 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2571#if defined(MFC_OpenACC)
2572# 1142 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2573!$acc loop seq
2574# 1142 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2575#elif defined(MFC_OpenMP)
2576# 1142 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2577
2578# 1142 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2579#endif
2580# 1143 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2581 do q = 1, num_species
2582 ys_ip(q) = ys_ip(q) + coeff*q_prim_vf(eqn_idx%species%beg + q - 1)%sf(i, j, k)
2583 end do
2584 end if
2585
2586 if (bubbles_euler .and. .not. qbmm) then
2587
2588# 1149 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2589#if defined(MFC_OpenACC)
2590# 1149 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2591!$acc loop seq
2592# 1149 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2593#elif defined(MFC_OpenMP)
2594# 1149 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2595
2596# 1149 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2597#endif
2598# 1150 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2599 do l = 1, nb
2600 if (polytropic) then
2601 r_ip(l) = r_ip(l) + coeff*q_prim_vf(eqn_idx%bub%beg + (l - 1)*2)%sf(i, j, k)
2602 v_ip(l) = v_ip(l) + coeff*q_prim_vf(eqn_idx%bub%beg + 1 + (l - 1)*2)%sf(i, j, k)
2603 else
2604 r_ip(l) = r_ip(l) + coeff*q_prim_vf(eqn_idx%bub%beg + (l - 1)*4)%sf(i, j, k)
2605 v_ip(l) = v_ip(l) + coeff*q_prim_vf(eqn_idx%bub%beg + 1 + (l - 1)*4)%sf(i, j, k)
2606 pb_ip(l) = pb_ip(l) + coeff*q_prim_vf(eqn_idx%bub%beg + 2 + (l - 1)*4)%sf(i, j, k)
2607 mv_ip(l) = mv_ip(l) + coeff*q_prim_vf(eqn_idx%bub%beg + 3 + (l - 1)*4)%sf(i, j, k)
2608 end if
2609 end do
2610 end if
2611
2612 if (qbmm) then
2613 do l = 1, nb*nmom
2614 nmom_ip(l) = nmom_ip(l) + coeff*q_prim_vf(eqn_idx%bub%beg - 1 + l)%sf(i, j, k)
2615 end do
2616 if (.not. polytropic) then
2617 do q = 1, nb
2618 do l = 1, nnode
2619 presb_ip((q - 1)*nnode + l) = presb_ip((q - 1)*nnode + l) + coeff*real(pb_in(i, j, k, l, q), &
2620 & kind=wp)
2621 massv_ip((q - 1)*nnode + l) = massv_ip((q - 1)*nnode + l) + coeff*real(mv_in(i, j, k, l, q), &
2622 & kind=wp)
2623 end do
2624 end do
2625 end if
2626 end if
2627 end do
2628 end do
2629 end do
2630
2631 end subroutine s_interpolate_image_point
2632
2633 !> Resets the current indexes of immersed boundaries and replaces them after updating
2634 !> the position of each moving immersed boundary
2635 impure subroutine s_update_mib(num_ibs)
2636
2637 integer, intent(in) :: num_ibs
2638 integer :: i, j, k, z_buff_size
2639
2640 call nvtxstartrange("UPDATE-MIBM")
2641
2642 ! Clears the existing immersed boundary indices over the same halo extent s_get_bounding_indices draws them
2643 z_buff_size = 0; if (p /= 0) z_buff_size = buff_size
2644
2645# 1195 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2646
2647# 1195 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2648#if defined(MFC_OpenACC)
2649# 1195 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2650!$acc parallel loop gang vector default(present) private(i, j, k)
2651# 1195 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2652#elif defined(MFC_OpenMP)
2653# 1195 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2654
2655# 1195 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2656
2657# 1195 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2658
2659# 1195 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2660!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, j, k)
2661# 1195 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2662#endif
2663# 1196 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2664 do i = -buff_size, m + buff_size; do j = -buff_size, n + buff_size; do k = -z_buff_size, p + z_buff_size
2665 ib_markers%sf(i, j, k) = 0._wp
2666 corrected_gps%sf(i, j, k) = 0
2667 end do; end do; end do
2668
2669# 1200 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2670#if defined(MFC_OpenACC)
2671# 1200 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2672!$acc end parallel loop
2673# 1200 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2674#elif defined(MFC_OpenMP)
2675# 1200 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2676
2677# 1200 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2678!$omp end target teams loop
2679# 1200 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2680#endif
2681# 1201 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2682
2683 ! recalulcate the rotation matrix based upon the new angles
2684
2685# 1203 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2686
2687# 1203 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2688#if defined(MFC_OpenACC)
2689# 1203 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2690!$acc parallel loop gang vector default(present) private(i)
2691# 1203 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2692#elif defined(MFC_OpenMP)
2693# 1203 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2694
2695# 1203 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2696
2697# 1203 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2698
2699# 1203 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2700!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i)
2701# 1203 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2702#endif
2703# 1204 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2704 do i = 1, num_ibs
2705 if (patch_ib(i)%moving_ibm /= 0) then
2706 call s_update_ib_rotation_matrix(i)
2707 end if
2708 end do
2709
2710# 1209 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2711#if defined(MFC_OpenACC)
2712# 1209 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2713!$acc end parallel loop
2714# 1209 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2715#elif defined(MFC_OpenMP)
2716# 1209 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2717
2718# 1209 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2719!$omp end target teams loop
2720# 1209 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2721#endif
2722# 1210 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2723
2724 ! recompute the new ib_patch locations and broadcast them.
2725 call nvtxstartrange("APPLY-IB-PATCHES")
2726 call s_apply_ib_patches(ib_markers)
2727 call nvtxendrange
2728
2729 call nvtxstartrange("COMPUTE-GHOST-POINTS")
2730 ! recalculate the ghost point locations and coefficients
2731 call s_find_num_ghost_points(num_gps)
2732 ! overrunning ghost_points is a silent device fault; name the numbers instead
2733 if (int(num_gps, 8) > int(size(ghost_points), 8)) then
2734 print '(A,I0,A,I0,A,I0)', 'ghost points on rank ', proc_rank, ' grew to ', num_gps, ' but the array holds ', &
2735 & size(ghost_points)
2736 call s_mpi_abort('moving immersed boundary: ghost-point count exceeded its allocation')
2737 end if
2738 ! num_gps is a declare-target module variable and bounds every device loop over the ghost points; it is
2739 ! copied to the device once in s_ibm_setup, so without this refresh the kernels keep the setup-time count
2740 ! as a moving body's count changes: stale list entries beyond the current count are written as wall
2741 ! states into cells that are now fluid, new ghost cells are left uncorrected, and which entries those are
2742 ! depends on the atomic fill order -- nondeterministic results run to run (#1886).
2743
2744# 1230 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2745#if defined(MFC_OpenACC)
2746# 1230 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2747!$acc update device(num_gps)
2748# 1230 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2749#elif defined(MFC_OpenMP)
2750# 1230 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2751!$omp target update to(num_gps)
2752# 1230 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2753#endif
2754# 1231 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2755 call s_find_ghost_points()
2756 call nvtxendrange
2757
2758 call nvtxstartrange("COMPUTE-IMAGE-POINTS")
2759 call s_apply_levelset(ghost_points, num_gps)
2760 call s_compute_image_points()
2761 call s_compute_interpolation_coeffs()
2762 call nvtxendrange
2763
2764 call nvtxendrange
2765
2766 end subroutine s_update_mib
2767
2768 !> log(cosh(x)) without overflow for large |x|
2769 pure function f_log_cosh(x) result(y)
2770
2771
2772# 1247 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2773#if MFC_OpenACC
2774# 1247 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2775!$acc routine seq
2776# 1247 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2777#elif MFC_OpenMP
2778# 1247 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2779
2780# 1247 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2781
2782# 1247 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2783!$omp declare target device_type(any)
2784# 1247 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2785#endif
2786# 1248 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2787 real(wp), intent(in) :: x
2788 real(wp) :: y
2789 y = abs(x) + log(0.5_wp*(1._wp + exp(-2._wp*abs(x))))
2790
2791 end function f_log_cosh
2792
2793 !> Prescribed kinematics for patch i at time t: hinged flapping (kin_model = 1) or the Eldredge pitch ramp (kin_model = 2).
2794 !! Flapping rolls phi about the lab x axis through the hinge and pitches theta about the body spanwise (y) axis through it; the
2795 !! ramp holds phi at zero and drives theta alone. Both compose as R = Rx(phi) Ry(theta) (MFC's angle convention) and set the
2796 !! angles, centroid, velocity and lab-frame angular velocity; nothing is integrated, so restarts are exact.
2797 subroutine s_prescribed_kinematics(i, t)
2798
2799
2800# 1260 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2801#if MFC_OpenACC
2802# 1260 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2803!$acc routine seq
2804# 1260 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2805#elif MFC_OpenMP
2806# 1260 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2807
2808# 1260 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2809
2810# 1260 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2811!$omp declare target device_type(any)
2812# 1260 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2813#endif
2814# 1261 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2815
2816 integer, intent(in) :: i
2817 real(wp), intent(in) :: t
2818 real(wp) :: tau, amp, damp, omega, phi, theta, phid, thetad, arg_phi, arg_th, t_p, a_s
2819 real(wp), dimension(3) :: r, c, w, v
2820
2821 tau = t - patch_ib(i)%kin_t0
2822 omega = 2._wp*pi*patch_ib(i)%kin_freq
2823 amp = 0._wp
2824 damp = 0._wp
2825 if (tau >= 0._wp) then
2826 amp = 1._wp
2827 if (patch_ib(i)%kin_ramp > 0._wp .and. tau < patch_ib(i)%kin_ramp) then
2828 amp = 0.5_wp*(1._wp - cos(pi*tau/patch_ib(i)%kin_ramp))
2829 damp = 0.5_wp*pi/patch_ib(i)%kin_ramp*sin(pi*tau/patch_ib(i)%kin_ramp)
2830 end if
2831 end if
2832 arg_phi = omega*tau
2833 arg_th = omega*tau + patch_ib(i)%kin_phase
2834
2835 if (patch_ib(i)%kin_model == 2) then
2836 ! Eldredge smoothed linear pitch ramp and hold (AIAA canonical pitch-up): with a = kin_smooth, nominal rate
2837 ! Omega = kin_pitch_rate and pitch time t_p = theta0/Omega, G = log[cosh(a tau)/cosh(a (tau - t_p))] rises from
2838 ! -a t_p to +a t_p, so theta = mean + theta0/2 [1 + G/(a t_p)] and theta' = Omega/2 [tanh(a tau) - tanh(a (tau - t_p))]
2839 phi = 0._wp
2840 phid = 0._wp
2841 t_p = patch_ib(i)%kin_theta0/patch_ib(i)%kin_pitch_rate
2842 a_s = patch_ib(i)%kin_smooth
2843 theta = patch_ib(i)%kin_theta_mean + 0.5_wp*patch_ib(i)%kin_theta0*(1._wp + (f_log_cosh(a_s*tau) &
2844 & - f_log_cosh(a_s*(tau - t_p)))/(a_s*t_p))
2845 thetad = 0.5_wp*patch_ib(i)%kin_pitch_rate*(tanh(a_s*tau) - tanh(a_s*(tau - t_p)))
2846 else
2847 phi = amp*patch_ib(i)%kin_phi0*sin(arg_phi)
2848 phid = damp*patch_ib(i)%kin_phi0*sin(arg_phi) + amp*patch_ib(i)%kin_phi0*omega*cos(arg_phi)
2849 theta = patch_ib(i)%kin_theta_mean + amp*patch_ib(i)%kin_theta0*sin(arg_th)
2850 thetad = damp*patch_ib(i)%kin_theta0*sin(arg_th) + amp*patch_ib(i)%kin_theta0*omega*cos(arg_th)
2851 end if
2852
2853 ! centroid relative to the hinge: Rx(phi) Ry(theta) offset. kin_offset points to the case-file centroid; a
2854 ! centroid moved to the centre of mass (geometries 4, 5, 11, 12) sits centroid_offset short of it
2855 r = patch_ib(i)%kin_offset - patch_ib(i)%centroid_offset
2856 c(1) = cos(theta)*r(1) + sin(theta)*r(3)
2857 c(2) = r(2)
2858 c(3) = -sin(theta)*r(1) + cos(theta)*r(3)
2859 v(1) = c(1)
2860 v(2) = cos(phi)*c(2) - sin(phi)*c(3)
2861 v(3) = sin(phi)*c(2) + cos(phi)*c(3)
2862 c = v
2863
2864 ! lab-frame angular velocity: phi' e_x + theta' Rx(phi) e_y
2865 w(1) = phid
2866 w(2) = thetad*cos(phi)
2867 w(3) = thetad*sin(phi)
2868 call s_cross_product(w, c, v)
2869
2870 patch_ib(i)%angles(1) = phi
2871 patch_ib(i)%angles(2) = theta
2872 patch_ib(i)%angles(3) = 0._wp
2873 patch_ib(i)%x_centroid = patch_ib(i)%kin_hinge(1) + c(1)
2874 patch_ib(i)%y_centroid = patch_ib(i)%kin_hinge(2) + c(2)
2875 patch_ib(i)%z_centroid = patch_ib(i)%kin_hinge(3) + c(3)
2876 patch_ib(i)%vel = v
2877 patch_ib(i)%angular_vel = w
2878
2879 end subroutine s_prescribed_kinematics
2880
2881 !> Compute pressure and viscous forces and torques on immersed bodies via volume integration
2882 subroutine s_compute_ib_forces(q_prim_vf, fluid_pp)
2883
2884 type(scalar_field), dimension(1:sys_size), intent(in) :: q_prim_vf
2885 type(physical_parameters), dimension(1:num_fluids), intent(in) :: fluid_pp
2886 integer :: i, j, k, l, encoded_ib_idx, xp, yp, zp, ib_idx, ib_idx_temp, fluid_idx
2887 real(wp), dimension(num_ibs, 3) :: forces, torques
2888 ! viscous stress tensor with temp vectors to hold divergence calculations
2889 real(wp), dimension(1:3,1:3) :: viscous_stress
2890 real(wp), dimension(1:3) :: local_force_contribution, radial_vector, local_torque_contribution
2891 real(wp) :: cell_volume, dynamic_viscosity
2892
2893# 1342 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2894 real(wp), dimension(num_fluids) :: dynamic_viscosities
2895# 1344 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2896
2897 call nvtxstartrange("COMPUTE-IB-FORCES")
2898
2899 forces = 0._wp
2900 torques = 0._wp
2901
2902 if (viscous) then
2903 do fluid_idx = 1, num_fluids
2904 if (fluid_pp(fluid_idx)%Re(1) > 0._wp) then
2905 dynamic_viscosities(fluid_idx) = 1._wp/fluid_pp(fluid_idx)%Re(1)
2906 else
2907 dynamic_viscosities(fluid_idx) = 0._wp
2908 end if
2909 end do
2910 end if
2911
2912 ! no kernel over zero-size mapped arrays: a rank can hold no patch after a handoff
2913 if (num_ibs > 0) then
2914
2915# 1362 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2916
2917# 1362 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2918#if defined(MFC_OpenACC)
2919# 1362 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2920!$acc parallel loop collapse(3) gang vector default(present) &
2921# 1362 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2922!$acc& private(i, j, k, l, xp, yp, zp, ib_idx, ib_idx_temp, encoded_ib_idx, fluid_idx, radial_vector, local_force_contribution, cell_volume, local_torque_contribution, dynamic_viscosity, viscous_stress) &
2923# 1362 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2924!$acc& copy(forces, torques) copyin(dynamic_viscosities)
2925# 1362 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2926#elif defined(MFC_OpenMP)
2927# 1362 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2928
2929# 1362 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2930
2931# 1362 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2932
2933# 1362 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2934!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) &
2935# 1362 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2936!$omp& private(i, j, k, l, xp, yp, zp, ib_idx, ib_idx_temp, encoded_ib_idx, fluid_idx, radial_vector, local_force_contribution, cell_volume, local_torque_contribution, dynamic_viscosity, viscous_stress) &
2937# 1362 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2938!$omp& map(tofrom:forces, torques) map(to:dynamic_viscosities)
2939# 1362 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2940#endif
2941# 1365 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
2942 do i = 0, m
2943 do j = 0, n
2944 do k = 0, p
2945 encoded_ib_idx = ib_markers%sf(i, j, k)
2946 if (encoded_ib_idx /= 0) then
2947 call s_decode_patch_periodicity(encoded_ib_idx, ib_idx_temp, xp, yp, zp)
2948 call s_get_neighborhood_idx(ib_idx_temp, ib_idx) ! global patch ID -> local index
2949 if (ib_idx > 0) then
2950 ! get the vector pointing to the grid cell from the IB centroid
2951 radial_vector(1) = x_cc(i) - (patch_ib(ib_idx)%x_centroid + real(xp, &
2952 & wp)*(glb_bounds(1)%end - glb_bounds(1)%beg))
2953 radial_vector(2) = y_cc(j) - (patch_ib(ib_idx)%y_centroid + real(yp, &
2954 & wp)*(glb_bounds(2)%end - glb_bounds(2)%beg))
2955 radial_vector(3) = 0._wp
2956 if (num_dims == 3) radial_vector(3) = z_cc(k) - (patch_ib(ib_idx)%z_centroid + real(zp, &
2957 & wp)*(glb_bounds(3)%end - glb_bounds(3)%beg))
2958
2959 local_force_contribution(:) = 0._wp
2960
2961 ! compute the pressure force component, which is the negative pressure gradient
2962 do l = -fd_number, fd_number
2963 local_force_contribution(1) = local_force_contribution(1) - (fd_coeff_x(l, &
2964 & i)*q_prim_vf(eqn_idx%E)%sf(i + l, j, k))
2965 local_force_contribution(2) = local_force_contribution(2) - (fd_coeff_y(l, &
2966 & j)*q_prim_vf(eqn_idx%E)%sf(i, j + l, k))
2967 if (num_dims == 3) then
2968 local_force_contribution(3) = local_force_contribution(3) - (fd_coeff_z(l, &
2969 & k)*q_prim_vf(eqn_idx%E)%sf(i, j, k + l))
2970 end if
2971 end do
2972
2973 ! get the viscous stress and add its contribution if that is considered
2974 if (viscous) then
2975 ! compute the volume-weighted local dynamic viscosity
2976 dynamic_viscosity = 0._wp
2977 do fluid_idx = 1, num_fluids
2978 ! local dynamic viscosity is the dynamic viscosity of the fluid times alpha of the fluid
2979 dynamic_viscosity = dynamic_viscosity + (q_prim_vf(fluid_idx + eqn_idx%adv%beg - 1)%sf(i, &
2980 & j, k)*dynamic_viscosities(fluid_idx))
2981 end do
2982
2983 do l = -fd_number, fd_number
2984 call s_compute_viscous_stress_tensor(viscous_stress, q_prim_vf, dynamic_viscosity, i + l, &
2985 & j, k)
2986 local_force_contribution(1:3) = local_force_contribution(1:3) + fd_coeff_x(l, &
2987 & i)*viscous_stress(1,1:3)
2988
2989 call s_compute_viscous_stress_tensor(viscous_stress, q_prim_vf, dynamic_viscosity, i, &
2990 & j + l, k)
2991 local_force_contribution(1:3) = local_force_contribution(1:3) + fd_coeff_y(l, &
2992 & j)*viscous_stress(2,1:3)
2993
2994 if (num_dims == 3) then
2995 call s_compute_viscous_stress_tensor(viscous_stress, q_prim_vf, dynamic_viscosity, i, &
2996 & j, k + l)
2997 local_force_contribution(1:3) = local_force_contribution(1:3) + fd_coeff_z(l, &
2998 & k)*viscous_stress(3,1:3)
2999 end if
3000 end do
3001 end if
3002
3003 call s_cross_product(radial_vector, local_force_contribution, local_torque_contribution)
3004
3005 ! Update the force and torque values atomically to prevent race conditions
3006 cell_volume = dx(i)*dy(j)
3007 if (num_dims == 3) cell_volume = cell_volume*dz(k)
3008 do l = 1, num_dims
3009
3010# 1432 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3011#if defined(MFC_OpenACC)
3012# 1432 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3013!$acc atomic update
3014# 1432 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3015#elif defined(MFC_OpenMP)
3016# 1432 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3017!$omp atomic update
3018# 1432 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3019#endif
3020# 1433 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3021 forces(ib_idx, l) = forces(ib_idx, l) + (local_force_contribution(l)*cell_volume)
3022
3023# 1434 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3024#if defined(MFC_OpenACC)
3025# 1434 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3026!$acc atomic update
3027# 1434 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3028#elif defined(MFC_OpenMP)
3029# 1434 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3030!$omp atomic update
3031# 1434 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3032#endif
3033# 1435 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3034 torques(ib_idx, l) = torques(ib_idx, l) + local_torque_contribution(l)*cell_volume
3035 end do
3036 end if ! ib_idx > 0
3037 end if
3038 end do
3039 end do
3040 end do
3041
3042# 1442 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3043#if defined(MFC_OpenACC)
3044# 1442 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3045!$acc end parallel loop
3046# 1442 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3047#elif defined(MFC_OpenMP)
3048# 1442 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3049
3050# 1442 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3051!$omp end target teams loop
3052# 1442 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3053#endif
3054# 1443 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3055 end if
3056
3057 call nvtxstartrange("APPLY-COLLISION-FORCES")
3058 call s_apply_collision_forces(ghost_points, num_gps, ib_markers, forces, torques)
3059 call nvtxendrange
3060
3061 ! reduce the forces across local neighborhood ranks
3062 call nvtxstartrange("COMMUNICATE-IB-FORCES")
3063 call s_communicate_ib_forces(forces, torques)
3064 call nvtxendrange
3065
3066 ! consider body forces after reducing to avoid double counting
3067 do i = 1, num_ibs
3068 if (bf_x) then
3069 forces(i, 1) = forces(i, 1) + accel_bf(1)*patch_ib(i)%mass
3070 end if
3071 if (bf_y) then
3072 forces(i, 2) = forces(i, 2) + accel_bf(2)*patch_ib(i)%mass
3073 end if
3074 if (bf_z) then
3075 forces(i, 3) = forces(i, 3) + accel_bf(3)*patch_ib(i)%mass
3076 end if
3077 end do
3078
3079 ! apply the summed forces
3080 if (num_ibs > 0) then
3081
3082# 1469 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3083
3084# 1469 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3085#if defined(MFC_OpenACC)
3086# 1469 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3087!$acc parallel loop gang vector default(present) private(i, l) copyin(forces, torques)
3088# 1469 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3089#elif defined(MFC_OpenMP)
3090# 1469 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3091
3092# 1469 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3093
3094# 1469 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3095
3096# 1469 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3097!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, l) &
3098# 1469 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3099!$omp& map(to:forces, torques)
3100# 1469 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3101#endif
3102# 1470 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3103 do i = 1, num_ibs
3104 do l = 1, 3
3105 patch_ib(i)%force(l) = forces(i, l)
3106 patch_ib(i)%torque(l) = torques(i, l)
3107 end do
3108 end do
3109
3110# 1476 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3111#if defined(MFC_OpenACC)
3112# 1476 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3113!$acc end parallel loop
3114# 1476 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3115#elif defined(MFC_OpenMP)
3116# 1476 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3117
3118# 1476 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3119!$omp end target teams loop
3120# 1476 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3121#endif
3122# 1477 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3123 end if
3124
3125 call nvtxendrange
3126
3127 end subroutine s_compute_ib_forces
3128
3129 !> Computes the center of mass for IB patch types where we are unable to determine their center of mass analytically.
3130 !> These patches include things like NACA airfoils and STL models
3131 !> Move a moving airfoil/STL patch's centroid to its centre of mass, keeping the difference in centroid_offset. Collective:
3132 !! every rank calls it for every global id, contributing zeros for patches it does not hold.
3133 subroutine s_compute_centroid_offset(gid, ib_marker)
3134
3135 integer, intent(in) :: gid !< global patch id
3136 integer, intent(in) :: ib_marker !< local index on this rank; <= 0 if not held
3137 integer :: i, j, k, num_cells_local, decoded_gbl_id
3138 integer(kind=8) :: num_cells, needs_loc, needs_glb
3139 real(wp), dimension(1:3) :: center_of_mass, center_of_mass_local
3140
3141 needs_loc = 0_8
3142 if (ib_marker > 0) then
3143 if (f_needs_centroid_offset(patch_ib(ib_marker))) needs_loc = 1_8
3144 end if
3145 call s_mpi_allreduce_integer_sum(needs_loc, needs_glb)
3146 if (needs_glb == 0_8) then
3147 if (ib_marker > 0) patch_ib(ib_marker)%centroid_offset(:) = 0._wp
3148 return
3149 end if
3150
3151 center_of_mass_local = [0._wp, 0._wp, 0._wp]
3152 num_cells_local = 0
3153 if (ib_marker > 0) then
3154 ! get the summed mass distribution and number of cells to divide by
3155 do i = 0, m
3156 do j = 0, n
3157 do k = 0, p
3158 if (ib_markers%sf(i, j, k) /= 0) then
3159 call s_decode_patch_periodicity(ib_markers%sf(i, j, k), decoded_gbl_id)
3160 if (decoded_gbl_id == gid) then
3161 num_cells_local = num_cells_local + 1
3162 center_of_mass_local = center_of_mass_local + [x_cc(i), y_cc(j), 0._wp]
3163 if (num_dims == 3) center_of_mass_local(3) = center_of_mass_local(3) + z_cc(k)
3164 end if
3165 end if
3166 end do
3167 end do
3168 end do
3169 end if
3170
3171 ! reduce the mass contribution over all MPI ranks and compute COM
3172 call s_mpi_allreduce_integer_sum(int(num_cells_local, 8), num_cells)
3173 call s_mpi_allreduce_sum(center_of_mass_local(1), center_of_mass(1))
3174 call s_mpi_allreduce_sum(center_of_mass_local(2), center_of_mass(2))
3175 call s_mpi_allreduce_sum(center_of_mass_local(3), center_of_mass(3))
3176 if (ib_marker <= 0) return
3177 if (num_cells == 0) then
3178 patch_ib(ib_marker)%centroid_offset = [0._wp, 0._wp, 0._wp]
3179 return
3180 end if
3181 center_of_mass = center_of_mass/real(num_cells, wp)
3182
3183 ! assign the centroid offset as a vector pointing from the true COM to the "centroid" in the input file and replace the
3184 ! current centroid
3185 patch_ib(ib_marker)%centroid_offset = [patch_ib(ib_marker)%x_centroid, patch_ib(ib_marker)%y_centroid, &
3186 & patch_ib(ib_marker)%z_centroid] - center_of_mass
3187 patch_ib(ib_marker)%x_centroid = center_of_mass(1)
3188 patch_ib(ib_marker)%y_centroid = center_of_mass(2)
3189 patch_ib(ib_marker)%z_centroid = center_of_mass(3)
3190
3191 ! rotate the centroid offset back into the local coords of the IB
3192 patch_ib(ib_marker)%centroid_offset = matmul(patch_ib(ib_marker)%rotation_matrix_inverse, &
3193 & patch_ib(ib_marker)%centroid_offset)
3194
3195 end subroutine s_compute_centroid_offset
3196
3197 !> A moving airfoil or STL (geometries 4, 5, 11, 12) moves its centroid to the centre of mass and keeps a centroid_offset
3198 pure logical function f_needs_centroid_offset(patch)
3199
3200 type(ib_patch_parameters), intent(in) :: patch
3201
3202 f_needs_centroid_offset = patch%moving_ibm /= 0 .and. any(patch%geometry == [4, 5, 11, 12])
3203
3204 end function f_needs_centroid_offset
3205
3206 !> On restart, replace the re-measured centroid offsets with the ones the run was using (restart_data/ib_offset_<step>.dat,
3207 !! written with each checkpoint), and re-place kinematics-driven bodies about them. No file: the measured offsets stand.
3208 impure subroutine s_restore_centroid_offsets(t_init)
3209
3210 real(wp), intent(in) :: t_init
3211 character(len=path_len + 2*name_len) :: file_loc
3212 logical :: file_exist
3213 integer :: gid, i, ios, file_unit, step
3214 real(wp), dimension(3) :: off
3215
3216 step = t_step_start
3217 if (cfl_dt) step = n_start
3218 if (step == 0) return
3219 write (file_loc, '(A,I0,A)') trim(case_dir) // '/restart_data/ib_offset_', step, '.dat'
3220 inquire (file=trim(file_loc), exist=file_exist)
3221 if (.not. file_exist) return
3222 open (newunit=file_unit, file=trim(file_loc), status='old', action='read', iostat=ios)
3223 if (ios /= 0) return
3224 do
3225 read (file_unit, *, iostat=ios) gid, off
3226 if (ios /= 0) exit
3227 call s_get_neighborhood_idx(gid, i)
3228 if (i > 0) then
3229 patch_ib(i)%centroid_offset = off
3230 if (patch_ib(i)%moving_ibm /= 0 .and. patch_ib(i)%kin_model > 0) call s_prescribed_kinematics(i, t_init)
3231 end if
3232 end do
3233 close (file_unit)
3234
3235 end subroutine s_restore_centroid_offsets
3236
3237 !> Computes the moment of inertia for an immersed boundary
3238 subroutine s_compute_moment_of_inertia(patch, axis, moment)
3239
3240
3241# 1594 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3242#if MFC_OpenACC
3243# 1594 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3244!$acc routine seq
3245# 1594 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3246#elif MFC_OpenMP
3247# 1594 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3248
3249# 1594 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3250
3251# 1594 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3252!$omp declare target device_type(any)
3253# 1594 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3254#endif
3255# 1595 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3256
3257 type(ib_patch_parameters), intent(in) :: patch
3258 real(wp), dimension(3), intent(in) :: axis
3259 real(wp), intent(out) :: moment
3260 real(wp) :: distance_to_axis, cell_volume
3261 real(wp), dimension(3) :: position, closest_point_along_axis, vector_to_axis, normal_axis
3262 integer :: i, j, k, count, ib_marker
3263
3264 ! if the IB is in 2D or a 3D sphere, we can compute this exactly
3265 if (patch%geometry == 2) then ! circle
3266 moment = 0.5_wp*patch%mass*(patch%radius)**2
3267 else if (patch%geometry == 3) then ! rectangle
3268 moment = patch%mass*(patch%length_x**2 + patch %length_y**2)/6._wp
3269 else if (patch%geometry == 6) then ! ellipse
3270 moment = 0.0625_wp*patch%mass*(patch%length_x**2 + patch %length_y**2)
3271 else if (patch%geometry == 8) then ! sphere
3272 moment = 0.4*patch%mass*(patch%radius)**2
3273 else ! we do not have an analytic moment of inertia calculation and need to approximate it directly via a sum
3274 count = 0
3275 cell_volume = (x_cc(1) - x_cc(0))*(y_cc(1) - y_cc(0))
3276 ! computed without grid stretching. Update in the loop to perform with stretching
3277 if (p /= 0) then
3278 cell_volume = cell_volume*(z_cc(1) - z_cc(0))
3279 end if
3280
3281 ib_marker = patch%gbl_patch_id
3282
3283 if (p == 0) then
3284 normal_axis = [0, 0, 1]
3285 else if (sqrt(sum(axis**2)) < sgm_eps) then
3286 ! if the object is not actually rotating at this time, return a dummy value and exit
3287 moment = 1._wp
3288 return
3289 else
3290 normal_axis = axis/sqrt(sum(axis**2))
3291 end if
3292
3293 moment = 0._wp
3294
3295 do i = 0, m
3296 do j = 0, n
3297 do k = 0, p
3298 if (ib_markers%sf(i, j, k) == ib_marker) then
3299 count = count + 1 ! increment the count of total cells in the boundary
3300
3301 ! get the position in local coordinates so that the axis passes through 0, 0, 0
3302 if (num_dims < 3) then
3303 position = [x_cc(i), y_cc(j), 0._wp] - [patch%x_centroid, patch%y_centroid, 0._wp]
3304 else
3305 position = [x_cc(i), y_cc(j), z_cc(k)] - [patch%x_centroid, patch%y_centroid, patch%z_centroid]
3306 end if
3307
3308 ! project the position along the axis to find the closest distance to the rotation axis
3309 closest_point_along_axis = normal_axis*dot_product(normal_axis, position)
3310 vector_to_axis = position - closest_point_along_axis
3311 distance_to_axis = dot_product(vector_to_axis, vector_to_axis) ! saves the distance to the axis squared
3312
3313 ! compute the position component of the moment
3314 moment = moment + distance_to_axis
3315 end if
3316 end do
3317 end do
3318 end do
3319
3320 ! write the final moment assuming the points are all uniform density
3321 moment = moment*patch%mass/(count*cell_volume)
3322 end if
3323
3324 end subroutine s_compute_moment_of_inertia
3325
3326 !> Wrap immersed boundary positions across periodic domain boundaries
3327 subroutine s_wrap_periodic_ibs()
3328
3329 integer :: patch_id
3330
3331
3332# 1670 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3333
3334# 1670 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3335#if defined(MFC_OpenACC)
3336# 1670 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3337!$acc parallel loop gang vector default(present) private(patch_id)
3338# 1670 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3339#elif defined(MFC_OpenMP)
3340# 1670 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3341
3342# 1670 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3343
3344# 1670 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3345
3346# 1670 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3347!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(patch_id)
3348# 1670 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3349#endif
3350# 1671 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3351 do patch_id = 1, num_ibs
3352 ! check domain wraps in x, y,
3353# 1674 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3354 if (num_dims >= 1) then
3355# 1675 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3356 ! check for periodicity
3357 if (ib_bc_x%beg == bc_periodic) then
3358# 1677 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3359 ! check if the boundary has left the domain, and then correct
3360 if (patch_ib(patch_id)%x_centroid < glb_bounds(1)%beg) then
3361# 1679 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3362 ! if the boundary exited "left", wrap it back around to the "right"
3363 patch_ib(patch_id)%x_centroid = patch_ib(patch_id)%x_centroid + (glb_bounds(1)%end &
3364# 1681 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3365 & - glb_bounds(1)%beg)
3366# 1682 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3367 else if (patch_ib(patch_id)%x_centroid > glb_bounds(1)%end) then
3368# 1683 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3369 ! if the boundary exited "right", wrap it back around to the "left"
3370 patch_ib(patch_id)%x_centroid = patch_ib(patch_id)%x_centroid - (glb_bounds(1)%end &
3371# 1685 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3372 & - glb_bounds(1)%beg)
3373# 1686 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3374 end if
3375 end if
3376 end if
3377# 1674 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3378 if (num_dims >= 2) then
3379# 1675 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3380 ! check for periodicity
3381 if (ib_bc_y%beg == bc_periodic) then
3382# 1677 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3383 ! check if the boundary has left the domain, and then correct
3384 if (patch_ib(patch_id)%y_centroid < glb_bounds(2)%beg) then
3385# 1679 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3386 ! if the boundary exited "left", wrap it back around to the "right"
3387 patch_ib(patch_id)%y_centroid = patch_ib(patch_id)%y_centroid + (glb_bounds(2)%end &
3388# 1681 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3389 & - glb_bounds(2)%beg)
3390# 1682 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3391 else if (patch_ib(patch_id)%y_centroid > glb_bounds(2)%end) then
3392# 1683 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3393 ! if the boundary exited "right", wrap it back around to the "left"
3394 patch_ib(patch_id)%y_centroid = patch_ib(patch_id)%y_centroid - (glb_bounds(2)%end &
3395# 1685 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3396 & - glb_bounds(2)%beg)
3397# 1686 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3398 end if
3399 end if
3400 end if
3401# 1674 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3402 if (num_dims >= 3) then
3403# 1675 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3404 ! check for periodicity
3405 if (ib_bc_z%beg == bc_periodic) then
3406# 1677 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3407 ! check if the boundary has left the domain, and then correct
3408 if (patch_ib(patch_id)%z_centroid < glb_bounds(3)%beg) then
3409# 1679 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3410 ! if the boundary exited "left", wrap it back around to the "right"
3411 patch_ib(patch_id)%z_centroid = patch_ib(patch_id)%z_centroid + (glb_bounds(3)%end &
3412# 1681 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3413 & - glb_bounds(3)%beg)
3414# 1682 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3415 else if (patch_ib(patch_id)%z_centroid > glb_bounds(3)%end) then
3416# 1683 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3417 ! if the boundary exited "right", wrap it back around to the "left"
3418 patch_ib(patch_id)%z_centroid = patch_ib(patch_id)%z_centroid - (glb_bounds(3)%end &
3419# 1685 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3420 & - glb_bounds(3)%beg)
3421# 1686 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3422 end if
3423 end if
3424 end if
3425# 1690 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3426 end do
3427
3428# 1691 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3429#if defined(MFC_OpenACC)
3430# 1691 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3431!$acc end parallel loop
3432# 1691 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3433#elif defined(MFC_OpenMP)
3434# 1691 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3435
3436# 1691 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3437!$omp end target teams loop
3438# 1691 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3439#endif
3440# 1692 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3441
3442 end subroutine s_wrap_periodic_ibs
3443
3444 !> Builds the ascending list of distinct neighborhood ranks used by s_communicate_ib_forces and allocates its buffers once. The
3445 !! ascending order is the fixed order in which an owner adds its neighbors' contributions.
3446 impure subroutine s_setup_ib_force_comm()
3447
3448#ifdef MFC_MPI
3449 integer :: i, j, r, max_recs
3450 integer, allocatable, dimension(:) :: flat
3451
3452 if (num_procs == 1 .or. allocated(ib_nbrs)) return
3453 if (num_gbl_ibs >= min(2_8**digits(0._wp), int(huge(0), 8))) then
3454# 1704 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3455 call s_prohibit_abort("num_gbl_ibs >= min(2_8**digits(0._wp), int(huge(0), 8))", "Too many IBs to encode their ids exactly in real(wp) force messages")
3456# 1704 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3457 end if
3458# 1706 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3459
3460 ! a rank can fill several table slots (periodicity, few ranks) or be its own neighbor; keep each distinct rank once
3461 flat = reshape(ib_neighbor_ranks, [size(ib_neighbor_ranks)])
3462 allocate (ib_nbrs(size(flat)))
3463 n_ib_nbrs = 0
3464 do i = 1, size(flat)
3465 r = flat(i)
3466 if (r < 0 .or. r == proc_rank) cycle
3467 if (any(ib_nbrs(1:n_ib_nbrs) == r)) cycle
3468 ! insert in ascending order
3469 j = n_ib_nbrs
3470 do while (j > 0)
3471 if (ib_nbrs(j) < r) exit
3472 ib_nbrs(j + 1) = ib_nbrs(j)
3473 j = j - 1
3474 end do
3475 ib_nbrs(j + 1) = r
3476 n_ib_nbrs = n_ib_nbrs + 1
3477 end do
3478
3479 ! an owner never receives more records than it owns, and a broadcast never carries more than its sender owns
3480 max_recs = min(size(patch_ib), num_local_ibs_max)
3481 allocate (ib_send_counts(n_ib_nbrs), ib_reqs(2*n_ib_nbrs), ib_stats(mpi_status_size, 2*n_ib_nbrs))
3482 allocate (ib_send_bufs(ib_rec_len*max_recs, max(n_ib_nbrs, 1)), ib_recv_bufs(ib_rec_len*max_recs, max(n_ib_nbrs, 1)))
3483#endif
3484
3485 end subroutine s_setup_ib_force_comm
3486
3487 !> Position of rank in ib_nbrs
3488 function f_ib_nbr_slot(rank) result(slot)
3489
3490 integer, intent(in) :: rank
3491 integer :: slot
3492
3493 do slot = 1, n_ib_nbrs
3494 if (ib_nbrs(slot) == rank) return
3495 end do
3496 call s_mpi_abort('IB owner is not a neighborhood rank')
3497
3498 end function f_ib_nbr_slot
3499
3500 !> Reduces forces and torques so every rank holding an IB ends with bit-identical values. Each IB has one owner, the rank whose
3501 !! subdomain holds its centroid (patch_ib%owner_rank, stamped at handoff). Phase 1: each rank sends its nonzero partials to the
3502 !! owner, which adds them to its own partial in ascending source-rank order once all have arrived. Phase 2: each owner sends its
3503 !! totals to every neighborhood rank, which overwrite their partials. One nonblocking message per neighbor per phase.
3504 subroutine s_communicate_ib_forces(forces, torques)
3505
3506 real(wp), dimension(num_ibs, 3), intent(inout) :: forces, torques
3507
3508#ifdef MFC_MPI
3509 integer :: i, j, s, r, slot, pos, nvals, ierr
3510
3511 if (num_procs == 1 .or. n_ib_nbrs == 0) return
3512
3513 ! Phase 1: send each partial to the owner of its IB
3514 do s = 1, n_ib_nbrs
3515 call mpi_irecv(ib_recv_bufs(:,s), size(ib_recv_bufs, 1), mpi_p, ib_nbrs(s), 610, mpi_comm_world, ib_reqs(s), ierr)
3516 end do
3517
3518 ib_send_counts = 0
3519 do i = 1, num_ibs
3520 if (patch_ib(i)%owner_rank == proc_rank) cycle
3521 if (all(forces(i,:) == 0._wp) .and. all(torques(i,:) == 0._wp)) cycle
3522 slot = f_ib_nbr_slot(patch_ib(i)%owner_rank)
3523 pos = ib_rec_len*ib_send_counts(slot)
3524 ib_send_bufs(pos + 1, slot) = real(patch_ib(i)%gbl_patch_id, wp)
3525 ib_send_bufs(pos + 2:pos + 4,slot) = forces(i,1:3)
3526 ib_send_bufs(pos + 5:pos + 7,slot) = torques(i,1:3)
3527 ib_send_counts(slot) = ib_send_counts(slot) + 1
3528 end do
3529
3530 ! every neighbor gets a message, possibly empty, so every posted receive completes
3531 do s = 1, n_ib_nbrs
3532 call mpi_isend(ib_send_bufs(:,s), ib_rec_len*ib_send_counts(s), mpi_p, ib_nbrs(s), 610, mpi_comm_world, &
3533 & ib_reqs(n_ib_nbrs + s), ierr)
3534 end do
3535 call mpi_waitall(2*n_ib_nbrs, ib_reqs, ib_stats, ierr)
3536
3537 ! owner sum: own partial first, then each neighbor in ascending rank order, whatever order the messages arrived in
3538 do s = 1, n_ib_nbrs
3539 call mpi_get_count(ib_stats(:,s), mpi_p, nvals, ierr)
3540 do r = 0, nvals/ib_rec_len - 1
3541 pos = ib_rec_len*r
3542 call s_get_neighborhood_idx(nint(ib_recv_bufs(pos + 1, s)), j)
3543 if (.not. (j > 0)) then
3544# 1790 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3545 call s_mpi_abort("m_ibm.fpp:1790: " // "Assertion failed: j > 0. " &
3546# 1790 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3547 & // 'IB force contribution for an IB this rank does not hold')
3548# 1790 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3549 end if
3550# 1791 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3551 if (.not. (patch_ib(j)%owner_rank == proc_rank)) then
3552# 1791 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3553 call s_mpi_abort("m_ibm.fpp:1791: " // "Assertion failed: patch_ib(j)%owner_rank == proc_rank. " &
3554# 1791 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3555 & // 'IB force contribution for an IB this rank does not own')
3556# 1791 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3557 end if
3558# 1792 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3559 forces(j,1:3) = forces(j,1:3) + ib_recv_bufs(pos + 2:pos + 4,s)
3560 torques(j,1:3) = torques(j,1:3) + ib_recv_bufs(pos + 5:pos + 7,s)
3561 end do
3562 end do
3563
3564 ! Phase 2: each owner sends its totals to every neighbor
3565 do s = 1, n_ib_nbrs
3566 call mpi_irecv(ib_recv_bufs(:,s), size(ib_recv_bufs, 1), mpi_p, ib_nbrs(s), 611, mpi_comm_world, ib_reqs(s), ierr)
3567 end do
3568
3569 do i = 1, num_local_ibs
3570 j = local_ib_patch_ids(i)
3571 pos = ib_rec_len*(i - 1)
3572 ib_send_bufs(pos + 1, 1) = real(patch_ib(j)%gbl_patch_id, wp)
3573 ib_send_bufs(pos + 2:pos + 4,1) = forces(j,1:3)
3574 ib_send_bufs(pos + 5:pos + 7,1) = torques(j,1:3)
3575 end do
3576
3577 do s = 1, n_ib_nbrs
3578 call mpi_isend(ib_send_bufs(:,1), ib_rec_len*num_local_ibs, mpi_p, ib_nbrs(s), 611, mpi_comm_world, &
3579 & ib_reqs(n_ib_nbrs + s), ierr)
3580 end do
3581 call mpi_waitall(2*n_ib_nbrs, ib_reqs, ib_stats, ierr)
3582
3583 ! holders take the owner's bits
3584 do s = 1, n_ib_nbrs
3585 call mpi_get_count(ib_stats(:,s), mpi_p, nvals, ierr)
3586 do r = 0, nvals/ib_rec_len - 1
3587 pos = ib_rec_len*r
3588 call s_get_neighborhood_idx(nint(ib_recv_bufs(pos + 1, s)), j)
3589 if (j <= 0) cycle
3590 if (.not. (patch_ib(j)%owner_rank == ib_nbrs(s))) then
3591# 1823 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3592 call s_mpi_abort("m_ibm.fpp:1823: " // "Assertion failed: patch_ib(j)%owner_rank == ib_nbrs(s). " &
3593# 1823 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3594 & // 'IB force total from a rank that does not own the IB')
3595# 1823 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3596 end if
3597# 1824 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3598 forces(j,1:3) = ib_recv_bufs(pos + 2:pos + 4,s)
3599 torques(j,1:3) = ib_recv_bufs(pos + 5:pos + 7,s)
3600 end do
3601 end do
3602#endif
3603
3604 end subroutine s_communicate_ib_forces
3605
3606 !> @brief Swaps ownership of IBs and passes ownership of IBs to neighbor processors
3607 subroutine s_handoff_ib_ownership()
3608
3609 integer :: i, j, k, output_idx, local_output_idx
3610 integer :: old_num_local_ibs, num_ibs_old, num_ibs_pre
3611 integer :: new_count, recv_count
3612 integer :: pack_pos, unpack_pos, buf_size, patch_bytes
3613 integer :: send_neighbor, recv_neighbor, ierr
3614 integer :: dx, dy, dz, tag, nbr_idx, nreqs
3615 real(wp), dimension(3) :: centroid
3616 logical :: is_new
3617 type(ib_patch_parameters) :: tmp_patch
3618 integer, dimension(num_local_ibs_max) :: local_ib_idx_old
3619 integer, dimension(num_ib_patches_max_namelist) :: old_to_new ! old patch_ib slot -> slot after compaction
3620 ! (2R+1)^num_dims - 1 neighbors in the radius-R neighborhood; each gets its own recv buffer
3621 integer :: max_nbrs
3622 character(len=1), allocatable :: send_buf(:), recv_bufs(:,:)
3623 integer, allocatable :: requests(:), recv_neighbor_list(:)
3624
3625#ifdef MFC_MPI
3626 call nvtxstartrange("HANDOFF-IB-OWNERSHIP")
3627
3628 if (num_procs > 1) then
3629 ! save a copy of the local IB's global indices to cross-reference for later.
3630 local_ib_idx_old = 0
3631 old_num_local_ibs = num_local_ibs
3632 do i = 1, num_local_ibs
3633 local_ib_idx_old(i) = patch_ib(local_ib_patch_ids(i))%gbl_patch_id
3634 end do
3635
3636 ! Sync GPU-updated fields (angles, angular_vel, centroids) to host before
3637 ! compaction and MPI packing, which read from host memory.
3638
3639# 1864 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3640#if defined(MFC_OpenACC)
3641# 1864 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3642!$acc update host(patch_ib)
3643# 1864 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3644#elif defined(MFC_OpenMP)
3645# 1864 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3646!$omp target update from(patch_ib)
3647# 1864 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3648#endif
3649# 1865 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3650
3651 ! delete any particles that no longer need to be tracked and coalesce the array
3652 num_ibs_old = num_ibs
3653 output_idx = 0
3654 local_output_idx = 0
3655 do i = 1, num_ibs
3656 old_to_new(i) = -1
3657 centroid = [patch_ib(i)%x_centroid, patch_ib(i)%y_centroid, 0._wp]
3658 if (num_dims == 3) centroid(3) = patch_ib(i)%z_centroid
3659
3660 ! delete if not in neighborhood
3661 if (f_neighborhood_ranks_own_location(centroid)) then
3662 output_idx = output_idx + 1
3663 old_to_new(i) = output_idx
3664 if (i /= output_idx) then
3665 patch_ib(output_idx) = patch_ib(i)
3666 end if
3667
3668 ! check if in local domain
3669 if (f_local_rank_owns_location(centroid, glb_bounds)) then
3670 local_output_idx = local_output_idx + 1
3671 if (local_output_idx > num_local_ibs_max) then
3672# 1886 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3673 call s_prohibit_abort("local_output_idx > num_local_ibs_max", "Too many IBs on a single processor rank. Modify case file or increase limit of num_local_ibs_max to resolve.")
3674# 1886 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3675 end if
3676# 1888 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3677 local_ib_patch_ids(local_output_idx) = output_idx
3678 patch_ib(output_idx)%owner_rank = proc_rank
3679 end if
3680 end if
3681 end do
3682 num_ibs = output_idx
3683 num_local_ibs = local_output_idx
3684 ! num_ibs shrinks here, so refresh it with patch_ib: s_get_neighborhood_idx bounds its search by it on the device
3685
3686# 1896 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3687#if defined(MFC_OpenACC)
3688# 1896 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3689!$acc update device(patch_ib, num_ibs)
3690# 1896 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3691#elif defined(MFC_OpenMP)
3692# 1896 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3693!$omp target update to(patch_ib, num_ibs)
3694# 1896 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3695#endif
3696# 1897 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3697 call s_compact_ib_lookup(old_to_new, num_ibs_old)
3698
3699 ! Broadcast newly-owned patches to all neighborhood neighbors
3700 patch_bytes = storage_size(tmp_patch)/8
3701 ! a rank can hand off at most every global patch, so size by num_gbl_ibs, not num_local_ibs_max
3702 buf_size = storage_size(0)/8 + patch_bytes*max(1, min(num_local_ibs_max, num_gbl_ibs))
3703 max_nbrs = (2*ib_neighborhood_radius + 1)**num_dims - 1
3704 allocate (send_buf(buf_size), recv_bufs(buf_size, max_nbrs), requests(2*max_nbrs), recv_neighbor_list(max_nbrs))
3705
3706 ! Write placeholder count at position 0
3707 pack_pos = 0
3708 call mpi_pack(0, 1, mpi_integer, send_buf, buf_size, pack_pos, mpi_comm_world, ierr)
3709
3710 ! pack new patches and count them
3711 new_count = 0
3712 do i = 1, num_local_ibs
3713 k = local_ib_patch_ids(i)
3714 is_new = .true.
3715 do j = 1, old_num_local_ibs
3716 if (patch_ib(k)%gbl_patch_id == local_ib_idx_old(j)) then
3717 is_new = .false.
3718 exit
3719 end if
3720 end do
3721 if (is_new) then
3722 call mpi_pack(patch_ib(k), patch_bytes, mpi_byte, send_buf, buf_size, pack_pos, mpi_comm_world, ierr)
3723 new_count = new_count + 1
3724 end if
3725 end do
3726
3727 ! Overwrite the placeholder with the real count
3728 pack_pos = 0
3729 call mpi_pack(new_count, 1, mpi_integer, send_buf, buf_size, pack_pos, mpi_comm_world, ierr)
3730 pack_pos = storage_size(0)/8 + new_count*patch_bytes
3731
3732 ! Post all receives first, then sends
3733 nreqs = 0
3734 nbr_idx = 0
3735 do dz = merge(-ib_neighborhood_radius, 0, num_dims == 3), merge(ib_neighborhood_radius, 0, num_dims == 3)
3736 do dy = -ib_neighborhood_radius, ib_neighborhood_radius
3737 do dx = -ib_neighborhood_radius, ib_neighborhood_radius
3738 if (dx == 0 .and. dy == 0 .and. dz == 0) cycle
3739 nbr_idx = nbr_idx + 1
3740 ! one tag per offset in the (2R+1)^3 neighbourhood: a radix-3 encoding collides once R > 1
3741 tag = 200 + ((dx + ib_neighborhood_radius)*(2*ib_neighborhood_radius + 1) + (dy + ib_neighborhood_radius)) &
3742 & *(2*ib_neighborhood_radius + 1) + (dz + ib_neighborhood_radius)
3743 recv_neighbor = ib_neighbor_ranks(-dx, -dy, -dz)
3744 recv_neighbor_list(nbr_idx) = mpi_proc_null
3745 if (recv_neighbor < 0) cycle
3746 recv_neighbor_list(nbr_idx) = recv_neighbor
3747 nreqs = nreqs + 1
3748 call mpi_irecv(recv_bufs(:,nbr_idx), buf_size, mpi_packed, recv_neighbor, tag, mpi_comm_world, &
3749 & requests(nreqs), ierr)
3750 end do
3751 end do
3752 end do
3753
3754 do dz = merge(-ib_neighborhood_radius, 0, num_dims == 3), merge(ib_neighborhood_radius, 0, num_dims == 3)
3755 do dy = -ib_neighborhood_radius, ib_neighborhood_radius
3756 do dx = -ib_neighborhood_radius, ib_neighborhood_radius
3757 if (dx == 0 .and. dy == 0 .and. dz == 0) cycle
3758 ! one tag per offset in the (2R+1)^3 neighbourhood: a radix-3 encoding collides once R > 1
3759 tag = 200 + ((dx + ib_neighborhood_radius)*(2*ib_neighborhood_radius + 1) + (dy + ib_neighborhood_radius)) &
3760 & *(2*ib_neighborhood_radius + 1) + (dz + ib_neighborhood_radius)
3761 send_neighbor = ib_neighbor_ranks(dx, dy, dz)
3762 if (send_neighbor < 0) cycle
3763 nreqs = nreqs + 1
3764 call mpi_isend(send_buf, pack_pos, mpi_packed, send_neighbor, tag, mpi_comm_world, requests(nreqs), ierr)
3765 end do
3766 end do
3767 end do
3768
3769 call mpi_waitall(nreqs, requests, mpi_statuses_ignore, ierr)
3770
3771 ! Unpack all received buffers
3772 num_ibs_pre = num_ibs
3773 do nbr_idx = 1, max_nbrs
3774 if (recv_neighbor_list(nbr_idx) == mpi_proc_null) cycle
3775 unpack_pos = 0
3776 call mpi_unpack(recv_bufs(:,nbr_idx), buf_size, unpack_pos, recv_count, 1, mpi_integer, mpi_comm_world, ierr)
3777 do i = 1, recv_count
3778 call mpi_unpack(recv_bufs(:,nbr_idx), buf_size, unpack_pos, tmp_patch, patch_bytes, mpi_byte, mpi_comm_world, &
3779 & ierr)
3780 call s_get_neighborhood_idx(tmp_patch%gbl_patch_id, j, num_ibs_pre)
3781 if (j < 0) then
3782 num_ibs = num_ibs + 1
3783 if (.not. (num_ibs <= size(patch_ib))) then
3784# 1983 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3785 call s_mpi_abort("m_ibm.fpp:1983: " // "Assertion failed: num_ibs <= size(patch_ib). " &
3786# 1983 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3787 & // 'patch_ib overflow in neighborhood handoff')
3788# 1983 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3789 end if
3790# 1984 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3791 patch_ib(num_ibs) = tmp_patch
3792 else
3793 ! an IB we already hold changed owner: record the new owner, keep our own copy of its state
3794 patch_ib(j)%owner_rank = tmp_patch%owner_rank
3795 end if
3796 end do
3797 end do
3798
3799 deallocate (send_buf, recv_bufs, requests, recv_neighbor_list)
3800
3801# 1993 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3802#if defined(MFC_OpenACC)
3803# 1993 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3804!$acc update device(patch_ib, num_ibs)
3805# 1993 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3806#elif defined(MFC_OpenMP)
3807# 1993 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3808!$omp target update to(patch_ib, num_ibs)
3809# 1993 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3810#endif
3811# 1994 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3812 call s_merge_ib_lookup(num_ibs_pre)
3813 end if
3814
3815 call nvtxendrange()
3816#endif
3817
3818 end subroutine s_handoff_ib_ownership
3819
3820 !> Abort if any namelist immersed boundary marks no cell anywhere. The owning rank is chosen by the patch centroid, so an STL
3821 !! placed with model_translate while its centroid is left elsewhere is handed to ranks that do not hold it, and it vanishes
3822 !! silently. Cloud particles are skipped: a particle's centroid is its body, so it cannot vanish this way.
3823 impure subroutine s_check_every_patch_marked()
3824
3825 integer(kind=8), allocatable :: cnt_loc(:), cnt_glb(:)
3826 integer :: gid, i, j, k, num_checked
3827
3828 num_checked = num_gbl_ibs - sum(particle_cloud(1:num_particle_clouds)%num_particles)
3829 if (num_checked <= 0) return
3830 allocate (cnt_loc(num_checked), cnt_glb(num_checked))
3831 cnt_loc = 0_8
3832 do k = 0, p
3833 do j = 0, n
3834 do i = 0, m
3835 if (ib_markers%sf(i, j, k) /= 0) then
3836 call s_decode_patch_periodicity(ib_markers%sf(i, j, k), gid)
3837 if (gid <= num_checked) cnt_loc(gid) = cnt_loc(gid) + 1_8
3838 end if
3839 end do
3840 end do
3841 end do
3842 call s_mpi_allreduce_integer_sum_vec(cnt_loc, cnt_glb)
3843 if (any(cnt_glb == 0_8)) then
3844# 2025 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3845 call s_prohibit_abort("any(cnt_glb == 0_8)", "An immersed boundary marked no cell anywhere: its centroid decides which rank owns it, so a body placed " // "elsewhere with model_translate is handed to a rank that does not hold it. Set patch_ib%x/y/z_centroid " // "to where the body actually is.")
3846# 2025 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3847 end if
3848# 2029 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3849 deallocate (cnt_loc, cnt_glb)
3850
3851 end subroutine s_check_every_patch_marked
3852
3853 !> Finalize the IBM module
3854 impure subroutine s_finalize_ibm_module()
3855
3856 integer :: i
3857
3858#ifdef MFC_DEBUG
3859# 2038 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3860 block
3861# 2038 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3862 use iso_fortran_env, only: output_unit
3863# 2038 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3864
3865# 2038 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3866 print *, 'm_ibm.fpp:2038: ', '@:DEALLOCATE(ib_markers%sf)'
3867# 2038 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3868
3869# 2038 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3870 call flush (output_unit)
3871# 2038 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3872 end block
3873# 2038 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3874#endif
3875# 2038 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3876
3877# 2038 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3878#if defined(MFC_OpenACC)
3879# 2038 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3880!$acc exit data delete(ib_markers%sf)
3881# 2038 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3882#elif defined(MFC_OpenMP)
3883# 2038 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3884!$omp target exit data map(release:ib_markers%sf)
3885# 2038 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3886#endif
3887# 2038 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3888 deallocate (ib_markers%sf)
3889# 2039 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3890#ifdef MFC_DEBUG
3891# 2039 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3892 block
3893# 2039 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3894 use iso_fortran_env, only: output_unit
3895# 2039 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3896
3897# 2039 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3898 print *, 'm_ibm.fpp:2039: ', '@:DEALLOCATE(corrected_gps%sf)'
3899# 2039 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3900
3901# 2039 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3902 call flush (output_unit)
3903# 2039 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3904 end block
3905# 2039 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3906#endif
3907# 2039 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3908
3909# 2039 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3910#if defined(MFC_OpenACC)
3911# 2039 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3912!$acc exit data delete(corrected_gps%sf)
3913# 2039 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3914#elif defined(MFC_OpenMP)
3915# 2039 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3916!$omp target exit data map(release:corrected_gps%sf)
3917# 2039 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3918#endif
3919# 2039 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3920 deallocate (corrected_gps%sf)
3921# 2040 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3922 do i = 1, num_ib_airfoils_max
3923 if (allocated(ib_airfoil_grids(i)%upper)) then
3924#ifdef MFC_DEBUG
3925# 2042 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3926 block
3927# 2042 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3928 use iso_fortran_env, only: output_unit
3929# 2042 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3930
3931# 2042 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3932 print *, 'm_ibm.fpp:2042: ', '@:DEALLOCATE(ib_airfoil_grids(i)%upper)'
3933# 2042 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3934
3935# 2042 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3936 call flush (output_unit)
3937# 2042 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3938 end block
3939# 2042 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3940#endif
3941# 2042 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3942
3943# 2042 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3944#if defined(MFC_OpenACC)
3945# 2042 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3946!$acc exit data delete(ib_airfoil_grids(i)%upper)
3947# 2042 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3948#elif defined(MFC_OpenMP)
3949# 2042 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3950!$omp target exit data map(release:ib_airfoil_grids(i)%upper)
3951# 2042 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3952#endif
3953# 2042 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3954 deallocate (ib_airfoil_grids(i)%upper)
3955# 2043 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3956#ifdef MFC_DEBUG
3957# 2043 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3958 block
3959# 2043 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3960 use iso_fortran_env, only: output_unit
3961# 2043 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3962
3963# 2043 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3964 print *, 'm_ibm.fpp:2043: ', '@:DEALLOCATE(ib_airfoil_grids(i)%lower)'
3965# 2043 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3966
3967# 2043 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3968 call flush (output_unit)
3969# 2043 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3970 end block
3971# 2043 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3972#endif
3973# 2043 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3974
3975# 2043 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3976#if defined(MFC_OpenACC)
3977# 2043 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3978!$acc exit data delete(ib_airfoil_grids(i)%lower)
3979# 2043 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3980#elif defined(MFC_OpenMP)
3981# 2043 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3982!$omp target exit data map(release:ib_airfoil_grids(i)%lower)
3983# 2043 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3984#endif
3985# 2043 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3986 deallocate (ib_airfoil_grids(i)%lower)
3987# 2044 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3988 end if
3989 end do
3990 if (allocated(models)) then
3991#ifdef MFC_DEBUG
3992# 2047 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3993 block
3994# 2047 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3995 use iso_fortran_env, only: output_unit
3996# 2047 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3997
3998# 2047 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
3999 print *, 'm_ibm.fpp:2047: ', '@:DEALLOCATE(models)'
4000# 2047 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4001
4002# 2047 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4003 call flush (output_unit)
4004# 2047 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4005 end block
4006# 2047 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4007#endif
4008# 2047 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4009
4010# 2047 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4011#if defined(MFC_OpenACC)
4012# 2047 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4013!$acc exit data delete(models)
4014# 2047 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4015#elif defined(MFC_OpenMP)
4016# 2047 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4017!$omp target exit data map(release:models)
4018# 2047 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4019#endif
4020# 2047 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4021 deallocate (models)
4022# 2048 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4023 end if
4024 if (allocated(ghost_points)) then
4025#ifdef MFC_DEBUG
4026# 2050 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4027 block
4028# 2050 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4029 use iso_fortran_env, only: output_unit
4030# 2050 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4031
4032# 2050 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4033 print *, 'm_ibm.fpp:2050: ', '@:DEALLOCATE(ghost_points)'
4034# 2050 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4035
4036# 2050 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4037 call flush (output_unit)
4038# 2050 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4039 end block
4040# 2050 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4041#endif
4042# 2050 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4043
4044# 2050 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4045#if defined(MFC_OpenACC)
4046# 2050 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4047!$acc exit data delete(ghost_points)
4048# 2050 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4049#elif defined(MFC_OpenMP)
4050# 2050 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4051!$omp target exit data map(release:ghost_points)
4052# 2050 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4053#endif
4054# 2050 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4055 deallocate (ghost_points)
4056# 2051 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4057 end if
4058 if (collision_model > 0) call s_finalize_collisions_module()
4059#ifdef MFC_MPI
4060 if (allocated(ib_nbrs)) deallocate (ib_nbrs, ib_send_counts, ib_reqs, ib_stats, ib_send_bufs, ib_recv_bufs)
4061#endif
4062
4063 end subroutine s_finalize_ibm_module
4064
4065 !> Species flux residual at a heterogeneous reacting surface.
4066 subroutine s_surface_species_residual(pres, T_s, d, Ys_IP, Ys_s, W_species, R_species, omega_s, mdot_s)
4067
4068
4069# 2062 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4070#if MFC_OpenACC
4071# 2062 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4072!$acc routine seq
4073# 2062 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4074#elif MFC_OpenMP
4075# 2062 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4076
4077# 2062 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4078
4079# 2062 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4080!$omp declare target device_type(any)
4081# 2062 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4082#endif
4083# 2063 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4084
4085 real(wp), intent(in) :: pres, T_s, d
4086 real(wp), intent(in) :: Ys_IP(num_species), Ys_s(num_species), W_species(num_species)
4087 real(wp), intent(out) :: R_species(num_species), omega_s(num_species), mdot_s
4088 real(wp) :: mw_IP, mw_s, rho_s, sum_BG
4089 real(wp) :: Xs_IP(num_species), Xs_s(num_species)
4090 real(wp) :: D_s(num_species), B_s(num_species), G_s(num_species)
4091 integer :: k
4092
4093 call get_mixture_molecular_weight(ys_ip, mw_ip)
4094 call get_mixture_molecular_weight(ys_s, mw_s)
4095 rho_s = pres*mw_s/(gas_constant*t_s)
4096
4097 do k = 1, num_species
4098 xs_ip(k) = ys_ip(k)*mw_ip/w_species(k)
4099 xs_s(k) = ys_s(k)*mw_s/w_species(k)
4100 end do
4101
4102 call get_species_mass_diffusivities_mixavg(pres, t_s, ys_s, d_s)
4103 call get_surface_net_production_rates(rho_s, t_s, ys_s, omega_s)
4104
4105 mdot_s = 0._wp
4106 do k = 1, num_species
4107 mdot_s = mdot_s + w_species(k)*omega_s(k)
4108 end do
4109
4110 sum_bg = 0._wp
4111 do k = 1, num_species
4112 b_s(k) = rho_s*d_s(k)*w_species(k)/mw_s
4113 g_s(k) = (xs_ip(k) - xs_s(k))/d
4114 sum_bg = sum_bg + b_s(k)*g_s(k)
4115 end do
4116
4117 do k = 1, num_species
4118 r_species(k) = -b_s(k)*g_s(k) + ys_s(k)*(sum_bg + mdot_s) - w_species(k)*omega_s(k)
4119 end do
4120
4121 end subroutine s_surface_species_residual
4122
4123 !> Surface energy residual: gas-side conduction balances heterogeneous reaction heat. Radiation and solid-side conduction are
4124 !! omitted.
4125 subroutine s_surface_energy_residual(pres, T_IP, T_s, d, Ys_s, R_energy)
4126
4127
4128# 2106 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4129#if MFC_OpenACC
4130# 2106 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4131!$acc routine seq
4132# 2106 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4133#elif MFC_OpenMP
4134# 2106 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4135
4136# 2106 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4137
4138# 2106 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4139!$omp declare target device_type(any)
4140# 2106 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4141#endif
4142# 2107 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4143
4144 real(wp), intent(in) :: pres, T_IP, T_s, d
4145 real(wp), intent(in) :: Ys_s(num_species)
4146 real(wp), intent(out) :: R_energy
4147 real(wp) :: mw_s, rho_s, k_s, q_rxn
4148
4149 call get_mixture_molecular_weight(ys_s, mw_s)
4150 rho_s = pres*mw_s/(gas_constant*t_s)
4151
4152 call get_mixture_thermal_conductivity_mixavg(t_s, ys_s, k_s)
4153 call get_surface_reaction_heat_flux(rho_s, t_s, ys_s, q_rxn)
4154
4155 r_energy = k_s*(t_s - t_ip)/d - q_rxn
4156
4157 end subroutine s_surface_energy_residual
4158
4159 !> Assemble the Newton residual for Ns species, with temperature appended only when it is solved. k_bath selects which species
4160 !! balance the sum constraint displaces; see s_solve_surface for why it must be the most abundant one.
4161 subroutine s_surface_residual(pres, T_IP, T_s, d, Ys_IP, Ys_s, W_species, k_bath, solve_temperature, flux_scale, &
4162 & energy_scale, R, R_species, omega_s, mdot_s)
4163
4164
4165# 2128 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4166#if MFC_OpenACC
4167# 2128 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4168!$acc routine seq
4169# 2128 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4170#elif MFC_OpenMP
4171# 2128 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4172
4173# 2128 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4174
4175# 2128 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4176!$omp declare target device_type(any)
4177# 2128 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4178#endif
4179# 2129 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4180
4181 real(wp), intent(in) :: pres, T_IP, T_s, d, flux_scale, energy_scale
4182 real(wp), intent(in) :: Ys_IP(num_species), Ys_s(num_species), W_species(num_species)
4183 integer, intent(in) :: k_bath
4184 logical, intent(in) :: solve_temperature
4185 real(wp), intent(out) :: R(num_species + 1), R_species(num_species), omega_s(num_species), mdot_s
4186 real(wp) :: R_energy
4187 integer :: k
4188
4189 call s_surface_species_residual(pres, t_s, d, ys_ip, ys_s, w_species, r_species, omega_s, mdot_s)
4190
4191 r = 0._wp
4192 do k = 1, num_species
4193 r(k) = r_species(k)/flux_scale
4194 end do
4195 r(k_bath) = sum(ys_s) - 1._wp
4196
4197 if (solve_temperature) then
4198 call s_surface_energy_residual(pres, t_ip, t_s, d, ys_s, r_energy)
4199 r(num_species + 1) = r_energy/energy_scale
4200 end if
4201
4202 end subroutine s_surface_residual
4203
4204 !> Ghost state for an immersed surface, from the image-point state and the surface state.
4205 !!
4206 !! The natural closure is the linear mirror phi_g = 2*phi_s - phi_IP, which is what makes the midpoint of the ghost/image
4207 !! pair reproduce the surface value. Its -1 coefficient on phi_IP is unconditional, though, so it drives strictly positive
4208 !! quantities negative whenever the surface value sits far below the image-point value: a cold wall in hot gas gives a
4209 !! negative ghost temperature and hence a negative ghost density, and a species the surface consumes faster than half the
4210 !! free-stream value gives a negative mass fraction. Rather than clamp the result -- which would break sum(Y) = 1 and hide
4211 !! the excursion -- the ghost state is the largest convex blend of the mirror back toward the surface value that stays
4212 !! physical:
4213 !!
4214 !! phi_g(theta) = phi_s + theta*(phi_s - phi_IP), theta in [0, 1]
4215 !!
4216 !! theta = 1 is the full second-order mirror; theta = 0 is the first-order Dirichlet ghost phi_g = phi_s, which is
4217 !! Gibou et al. (JCP 176:205, 2002) Eq. 17 and is what they likewise fall back to where the linear form is ill-behaved.
4218 !! One theta is shared by every species, so since both endpoints satisfy sum(Y) = 1, so does every blend between them --
4219 !! exactly, with no renormalization. Temperature gets its own theta: nothing couples it to the composition, and sharing
4220 !! would let a trace radical the surface consumes to ~1e-9 (O and OH at a burning carbon wall) drag the thermal mirror
4221 !! down with it, turning a 2100 K ghost into 1205 K over a species whose own excursion was 1e-9. The cost of limiting
4222 !! is small: the mirror already delivers only first-order wall flux (Ezra et al., Int. J. Heat Mass Transfer, 2025),
4223 !! and a boundary closure one order below the interior scheme retains the interior convergence rate (Gustafsson,
4224 !! Math. Comp. 29:396, 1975).
4225 !!
4226 !! This bounds the ghost state; it does not make the boundary strictly conservative. Baskaya et al. (Computers & Fluids
4227 !! 270:106134, 2024) trace the same cold-wall, large-gradient regime in ablation to ghost-cell mass conservation error
4228 !! that surfaces as spurious blowing, and resolve it only by moving to a flux-based cut-cell boundary -- a different
4229 !! discretization from this one, not a tuning of it.
4230 subroutine s_blend_ghost_temperature(T_IP, T_s, T_lo, T_hi, T_g)
4231
4232
4233# 2181 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4234#if MFC_OpenACC
4235# 2181 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4236!$acc routine seq
4237# 2181 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4238#elif MFC_OpenMP
4239# 2181 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4240
4241# 2181 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4242
4243# 2181 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4244!$omp declare target device_type(any)
4245# 2181 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4246#endif
4247# 2182 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4248
4249 real(wp), intent(in) :: T_IP, T_s, T_lo, T_hi
4250 real(wp), intent(out) :: T_g
4251 ! Stop short of the bound so roundoff cannot cross it.
4252 real(wp), parameter :: blend_safety = 0.9_wp
4253 real(wp) :: theta_T
4254
4255 ! Largest theta <= 1 keeping T_g in [T_lo, T_hi]: the NASA fit range with chemistry, positivity without.
4256 ! The bounds are compared before dividing, so T_hi = huge() cannot overflow.
4257 theta_t = 1._wp
4258 if (t_s < t_ip .and. blend_safety*(t_s - t_lo) < t_ip - t_s) then
4259 theta_t = blend_safety*(t_s - t_lo)/(t_ip - t_s)
4260 else if (t_s > t_ip .and. blend_safety*(t_hi - t_s) < t_s - t_ip) then
4261 theta_t = blend_safety*(t_hi - t_s)/(t_s - t_ip)
4262 end if
4263 theta_t = max(theta_t, 0._wp)
4264
4265 t_g = t_s + theta_t*(t_s - t_ip)
4266
4267 end subroutine s_blend_ghost_temperature
4268
4269 !> Species half of the ghost-state blend: one theta shared by every species (see s_blend_ghost_temperature).
4270 subroutine s_blend_ghost_species(Ys_IP, Ys_s, Ys_g)
4271
4272
4273# 2206 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4274#if MFC_OpenACC
4275# 2206 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4276!$acc routine seq
4277# 2206 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4278#elif MFC_OpenMP
4279# 2206 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4280
4281# 2206 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4282
4283# 2206 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4284!$omp declare target device_type(any)
4285# 2206 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4286#endif
4287# 2207 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4288
4289 real(wp), intent(in) :: Ys_IP(num_species), Ys_s(num_species)
4290 real(wp), intent(out) :: Ys_g(num_species)
4291 real(wp), parameter :: blend_safety = 0.9_wp
4292 real(wp) :: theta_Y
4293 integer :: k
4294
4295 theta_y = 1._wp
4296 do k = 1, num_species
4297 if (ys_s(k) < ys_ip(k)) theta_y = min(theta_y, blend_safety*ys_s(k)/(ys_ip(k) - ys_s(k)))
4298 end do
4299
4300 do k = 1, num_species
4301 ys_g(k) = ys_s(k) + theta_y*(ys_s(k) - ys_ip(k))
4302 end do
4303
4304 end subroutine s_blend_ghost_species
4305
4306 !> Index of the most abundant species, i.e. the balance the sum constraint displaces.
4307 !!
4308 !! Closing sum(Y) = 1 by dropping one species' flux balance is standard (Surface CHEMKIN), but the dropped species absorbs
4309 !! the roundoff of every other balance, so it must be the bath gas. Dropping a fixed index instead -- the last species in
4310 !! the mechanism -- lands that error on whatever the mechanism happens to list last, which for the reduced GRI mechanism
4311 !! shipped with the reacting-surface example is H2O2, a trace radical whose own surface balance then goes unenforced.
4312 !! Cantera's solveSP re-scans for the largest species each iteration (evalSurfLarge) for this reason.
4313 subroutine s_pick_bath_species(Ys_s, k_bath)
4314
4315
4316# 2234 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4317#if MFC_OpenACC
4318# 2234 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4319!$acc routine seq
4320# 2234 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4321#elif MFC_OpenMP
4322# 2234 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4323
4324# 2234 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4325
4326# 2234 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4327!$omp declare target device_type(any)
4328# 2234 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4329#endif
4330# 2235 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4331
4332 real(wp), intent(in) :: Ys_s(num_species)
4333 integer, intent(out) :: k_bath
4334 integer :: k
4335
4336 k_bath = 1
4337 do k = 2, num_species
4338 if (ys_s(k) > ys_s(k_bath)) k_bath = k
4339 end do
4340
4341 end subroutine s_pick_bath_species
4342
4343 !> Newton solve for a reacting surface. thermal_bc=0: zero-normal-gradient T; 1: prescribed T; 2: energy balance.
4344 subroutine s_solve_surface(pres, T_IP, T_wall, d, Ys_IP, W_species, thermal_bc, Ys_s, T_s, mdot_s, converged)
4345
4346
4347# 2250 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4348#if MFC_OpenACC
4349# 2250 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4350!$acc routine seq
4351# 2250 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4352#elif MFC_OpenMP
4353# 2250 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4354
4355# 2250 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4356
4357# 2250 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4358!$omp declare target device_type(any)
4359# 2250 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4360#endif
4361# 2251 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4362
4363 real(wp), intent(in) :: pres, T_IP, T_wall, d
4364 real(wp), intent(in) :: Ys_IP(num_species), W_species(num_species)
4365 integer, intent(in) :: thermal_bc
4366 real(wp), intent(out) :: Ys_s(num_species), T_s, mdot_s
4367 logical, intent(out) :: converged
4368 integer, parameter :: max_iter = 30, max_backtrack = 20
4369 real(wp), parameter :: fd_eps_Y = 1.e-7_wp, fd_eps_t = 1.e-6_wp
4370 ! The Jacobian is a forward difference, so the Newton update carries O(sqrt(epsilon)) noise and the iteration turns
4371 ! linear near the root. The residual itself is evaluated exactly, so a tighter tolerance is not unreachable, but it is
4372 ! reached by stalling rather than converging -- which the caller then reports as a failed surface. Cantera's solveSP
4373 ! runs the same 1e-7 forward-difference step against a 1e-4 relative tolerance; 1e-6 here stays well inside that.
4374 real(wp), parameter :: tol = 1.e-6_wp
4375 real(wp), parameter :: Y_tol = 100._wp*epsilon(1._wp)
4376 real(wp) :: A(num_species + 1, num_species + 1), rhs(num_species + 1), delta(num_species + 1)
4377 real(wp) :: R(num_species + 1), R_pert(num_species + 1), R_trial(num_species + 1)
4378 real(wp) :: R_species(num_species), R_species_pert(num_species), R_species_trial(num_species)
4379 real(wp) :: omega_s(num_species), omega_pert(num_species), omega_trial(num_species)
4380 real(wp) :: Ys_pert(num_species), Ys_trial(num_species)
4381 real(wp) :: mdot_pert, mdot_trial, T_pert, T_trial
4382 real(wp) :: flux_scale, energy_scale, R_energy, dx, lambda, norm_R, norm_trial
4383 logical :: solve_temperature, linear_success, accepted
4384 integer :: nsolve, iter, j, iback, k_bath, k_bath_prev
4385
4386 converged = .false.
4387 ys_s = ys_ip
4388
4389 select case (thermal_bc)
4390 case (0)
4391 t_s = t_ip
4392 solve_temperature = .false.
4393 case (1)
4394 t_s = t_wall
4395 solve_temperature = .false.
4396 case (2)
4397 t_s = min(max(t_ip, t_surface_min), t_surface_max)
4398 solve_temperature = .true.
4399 case default
4400 t_s = t_ip
4401 omega_s = 0._wp
4402 mdot_s = 0._wp
4403 return
4404 end select
4405
4406 nsolve = num_species + merge(1, 0, solve_temperature)
4407
4408 call s_pick_bath_species(ys_s, k_bath)
4409
4410 call s_surface_species_residual(pres, t_s, d, ys_ip, ys_s, w_species, r_species, omega_s, mdot_s)
4411 flux_scale = max(maxval(abs(r_species)), 1.e-12_wp)
4412 energy_scale = 1._wp
4413 if (solve_temperature) then
4414 call s_surface_energy_residual(pres, t_ip, t_s, d, ys_s, r_energy)
4415 energy_scale = max(abs(r_energy), 1._wp)
4416 end if
4417
4418 call s_surface_residual(pres, t_ip, t_s, d, ys_ip, ys_s, w_species, k_bath, solve_temperature, flux_scale, energy_scale, &
4419 & r, r_species, omega_s, mdot_s)
4420 norm_r = maxval(abs(r(1:nsolve)))
4421 if (norm_r < tol) then
4422 converged = .true.
4423 return
4424 end if
4425
4426 do iter = 1, max_iter
4427 do j = 1, num_species
4428 ys_pert = ys_s
4429 t_pert = t_s
4430 dx = fd_eps_y*max(abs(ys_s(j)), 1._wp)
4431 if (ys_s(j) + dx > 1._wp) dx = -dx
4432 ys_pert(j) = ys_pert(j) + dx
4433
4434 call s_surface_residual(pres, t_ip, t_pert, d, ys_ip, ys_pert, w_species, k_bath, solve_temperature, flux_scale, &
4435 & energy_scale, r_pert, r_species_pert, omega_pert, mdot_pert)
4436 a(1:nsolve,j) = (r_pert(1:nsolve) - r(1:nsolve))/dx
4437 end do
4438
4439 if (solve_temperature) then
4440 ys_pert = ys_s
4441 dx = fd_eps_t*max(abs(t_s), 1._wp)
4442 t_pert = t_s + dx
4443 if (t_pert > t_surface_max) then
4444 dx = -dx
4445 t_pert = t_s + dx
4446 end if
4447
4448 call s_surface_residual(pres, t_ip, t_pert, d, ys_ip, ys_pert, w_species, k_bath, solve_temperature, flux_scale, &
4449 & energy_scale, r_pert, r_species_pert, omega_pert, mdot_pert)
4450 a(1:nsolve,nsolve) = (r_pert(1:nsolve) - r(1:nsolve))/dx
4451 end if
4452
4453 rhs(1:nsolve) = -r(1:nsolve)
4454 call s_solve_surface_linear_system(a, rhs, delta, nsolve, linear_success)
4455 if (.not. linear_success) return
4456
4457 lambda = 1._wp
4458 accepted = .false.
4459 do iback = 1, max_backtrack
4460 ys_trial = ys_s + lambda*delta(1:num_species)
4461 t_trial = t_s
4462 if (solve_temperature) t_trial = t_s + lambda*delta(nsolve)
4463
4464 if (minval(ys_trial) < -y_tol .or. maxval(ys_trial) > 1._wp + y_tol .or. t_trial < t_surface_min &
4465 & .or. t_trial > t_surface_max) then
4466 lambda = 0.5_wp*lambda
4467 cycle
4468 end if
4469
4470 where (ys_trial < 0._wp) ys_trial = 0._wp
4471 where (ys_trial > 1._wp) ys_trial = 1._wp
4472
4473 call s_surface_residual(pres, t_ip, t_trial, d, ys_ip, ys_trial, w_species, k_bath, solve_temperature, &
4474 & flux_scale, energy_scale, r_trial, r_species_trial, omega_trial, mdot_trial)
4475 norm_trial = maxval(abs(r_trial(1:nsolve)))
4476 if (norm_trial < norm_r) then
4477 accepted = .true.
4478 exit
4479 end if
4480 lambda = 0.5_wp*lambda
4481 end do
4482
4483 if (.not. accepted) return
4484
4485 ys_s = ys_trial
4486 t_s = t_trial
4487 r = r_trial
4488 r_species = r_species_trial
4489 omega_s = omega_trial
4490 mdot_s = mdot_trial
4491 norm_r = norm_trial
4492
4493 ! The composition has moved, so the species the sum constraint displaces may no longer be the most abundant one.
4494 ! Re-closing changes what R means, so the residual is rebuilt on the rare steps where the choice actually changes.
4495 k_bath_prev = k_bath
4496 call s_pick_bath_species(ys_s, k_bath)
4497 if (k_bath /= k_bath_prev) then
4498 call s_surface_residual(pres, t_ip, t_s, d, ys_ip, ys_s, w_species, k_bath, solve_temperature, flux_scale, &
4499 & energy_scale, r, r_species, omega_s, mdot_s)
4500 norm_r = maxval(abs(r(1:nsolve)))
4501 end if
4502
4503 if (norm_r < tol) then
4504 converged = .true.
4505 return
4506 end if
4507 end do
4508
4509 end subroutine s_solve_surface
4510
4511 !> Small dense linear solve with partial pivoting for the local surface Newton system.
4512 subroutine s_solve_surface_linear_system(A, b, x, nsolve, success)
4513
4514
4515# 2403 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4516#if MFC_OpenACC
4517# 2403 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4518!$acc routine seq
4519# 2403 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4520#elif MFC_OpenMP
4521# 2403 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4522
4523# 2403 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4524
4525# 2403 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4526!$omp declare target device_type(any)
4527# 2403 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4528#endif
4529# 2404 "/home/runner/work/MFC/MFC/src/simulation/m_ibm.fpp"
4530
4531 integer, intent(in) :: nsolve
4532 real(wp), intent(inout) :: A(num_species + 1, num_species + 1)
4533 real(wp), intent(inout) :: b(num_species + 1)
4534 real(wp), intent(out) :: x(num_species + 1)
4535 logical, intent(out) :: success
4536 real(wp) :: factor, pivot_value, tmp, row_tmp(num_species + 1)
4537 integer :: i, j, k, pivot
4538
4539 success = .true.
4540 x = 0._wp
4541
4542 do k = 1, nsolve - 1
4543 pivot = k
4544 pivot_value = abs(a(k, k))
4545 do i = k + 1, nsolve
4546 if (abs(a(i, k)) > pivot_value) then
4547 pivot = i
4548 pivot_value = abs(a(i, k))
4549 end if
4550 end do
4551 if (pivot_value <= epsilon(1._wp)) then
4552 success = .false.
4553 return
4554 end if
4555
4556 if (pivot /= k) then
4557 row_tmp(1:nsolve) = a(k,1:nsolve)
4558 a(k,1:nsolve) = a(pivot,1:nsolve)
4559 a(pivot,1:nsolve) = row_tmp(1:nsolve)
4560 tmp = b(k)
4561 b(k) = b(pivot)
4562 b(pivot) = tmp
4563 end if
4564
4565 do i = k + 1, nsolve
4566 factor = a(i, k)/a(k, k)
4567 a(i, k) = 0._wp
4568 do j = k + 1, nsolve
4569 a(i, j) = a(i, j) - factor*a(k, j)
4570 end do
4571 b(i) = b(i) - factor*b(k)
4572 end do
4573 end do
4574
4575 if (abs(a(nsolve, nsolve)) <= epsilon(1._wp)) then
4576 success = .false.
4577 return
4578 end if
4579
4580 x(nsolve) = b(nsolve)/a(nsolve, nsolve)
4581 do i = nsolve - 1, 1, -1
4582 tmp = b(i)
4583 do j = i + 1, nsolve
4584 tmp = tmp - a(i, j)*x(j)
4585 end do
4586 if (abs(a(i, i)) <= epsilon(1._wp)) then
4587 success = .false.
4588 return
4589 end if
4590 x(i) = tmp/a(i, i)
4591 end do
4592
4593 end subroutine s_solve_surface_linear_system
4594
4595end module m_ibm
q_cons_vf
type(scalar_field), dimension(sys_size), intent(inout) q_cons_vf
Definition m_phase_change.fpp.f90:1047
k
integer, intent(in) k
Definition m_phase_change.fpp.f90:1049
j
integer, intent(in) j
Definition m_phase_change.fpp.f90:1049
l
integer, intent(in) l
Definition m_phase_change.fpp.f90:1049
m_collisions
Ghost-node immersed boundary method: locates ghost/image points, computes interpolation coefficients,...
Definition m_collisions.fpp.f90:357
m_compute_levelset
Computes signed-distance level-set fields and surface normals for immersed-boundary patch geometries.
Definition m_compute_levelset.fpp.f90:356
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:356
m_eos
Equations of state in Gamma/Pi form, rho e = Gamma(rho) p + Pi(rho).
Definition m_eos.fpp.f90:402
m_global_parameters
Global parameters for the computational domain, fluid properties, and simulation algorithm configurat...
Definition m_global_parameters.fpp.f90:367
m_global_parameters::buff_size
integer buff_size
Number of ghost cells for boundary condition storage.
Definition m_global_parameters.fpp.f90:853
m_helper_basic
Basic floating-point utilities: approximate equality, default detection, and coordinate bounds.
Definition m_helper_basic.fpp.f90:356
m_helper
Utility routines for bubble model setup, coordinate transforms, array sampling, and special functions...
Definition m_helper.fpp.f90:367
m_ib_patches
Allocate memory and read initial condition data for IC extrusion.
Definition m_ib_patches.fpp.f90:426
m_ibm
Ghost-node immersed boundary method: locates ghost/image points, computes interpolation coefficients,...
Definition m_ibm.fpp.f90:357
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:2636
m_ibm::s_surface_species_residual
subroutine s_surface_species_residual(pres, t_s, d, ys_ip, ys_s, w_species, r_species, omega_s, mdot_s)
Species flux residual at a heterogeneous reacting surface.
Definition m_ibm.fpp.f90:4067
m_ibm::s_compute_centroid_offset
subroutine s_compute_centroid_offset(gid, 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:3134
m_ibm::s_find_ghost_points
subroutine, private s_find_ghost_points()
Locate all ghost points in the domain.
Definition m_ibm.fpp.f90:2089
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, ys_ip)
Interpolate primitive variables to a ghost point's image point using bilinear or trilinear interpolat...
Definition m_ibm.fpp.f90:2419
m_ibm::s_compute_image_points
impure subroutine, private s_compute_image_points()
Compute the image points for each ghost point.
Definition m_ibm.fpp.f90:1845
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:2883
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:763
m_ibm::s_blend_ghost_species
subroutine s_blend_ghost_species(ys_ip, ys_s, ys_g)
Species half of the ghost-state blend: one theta shared by every species (see s_blend_ghost_temperatu...
Definition m_ibm.fpp.f90:4271
m_ibm::n_surface_not_converged
integer n_surface_not_converged
Ghost points at which the reacting-surface Newton solve did not reach its tolerance,...
Definition m_ibm.fpp.f90:387
m_ibm::s_communicate_ib_forces
subroutine s_communicate_ib_forces(forces, torques)
Reduces forces and torques so every rank holding an IB ends with bit-identical values....
Definition m_ibm.fpp.f90:3505
m_ibm::s_pick_bath_species
subroutine s_pick_bath_species(ys_s, k_bath)
Index of the most abundant species, i.e. the balance the sum constraint displaces.
Definition m_ibm.fpp.f90:4314
m_ibm::corrected_gps
type(integer_field) corrected_gps
1 at ghost points whose image point reaches the fluid, 0 elsewhere. Ghost points with a buried image ...
Definition m_ibm.fpp.f90:408
m_ibm::s_handoff_ib_ownership
subroutine s_handoff_ib_ownership()
Swaps ownership of IBs and passes ownership of IBs to neighbor processors.
Definition m_ibm.fpp.f90:3608
m_ibm::s_finalize_ibm_module
impure subroutine, public s_finalize_ibm_module()
Finalize the IBM module.
Definition m_ibm.fpp.f90:3855
m_ibm::s_report_ibm_surface
impure subroutine, public s_report_ibm_surface
One line at the end of a run if the reacting surface ever failed to solve. Silence means every ghost ...
Definition m_ibm.fpp.f90:1820
m_ibm::ghost_points
type(ghost_point), dimension(:), allocatable ghost_points
Definition m_ibm.fpp.f90:422
m_ibm::s_solve_surface
subroutine s_solve_surface(pres, t_ip, t_wall, d, ys_ip, w_species, thermal_bc, ys_s, t_s, mdot_s, converged)
Newton solve for a reacting surface. thermal_bc=0: zero-normal-gradient T; 1: prescribed T; 2: energy...
Definition m_ibm.fpp.f90:4345
m_ibm::num_gps
integer num_gps
Number of ghost points.
Definition m_ibm.fpp.f90:436
m_ibm::ib_send_counts
integer, dimension(:), allocatable, private ib_send_counts
Definition m_ibm.fpp.f90:470
m_ibm::s_compute_moment_of_inertia
subroutine s_compute_moment_of_inertia(patch, axis, moment)
Computes the moment of inertia for an immersed boundary.
Definition m_ibm.fpp.f90:3239
m_ibm::moving_immersed_boundary_flag
logical moving_immersed_boundary_flag
Definition m_ibm.fpp.f90:464
m_ibm::s_compute_ghost_point_velocity
subroutine, private s_compute_ghost_point_velocity(gp, gp_patch_id, radial_vector, vel_ip, pres_ip, vel_gp)
Definition m_ibm.fpp.f90:1050
m_ibm::f_ib_nbr_slot
integer function f_ib_nbr_slot(rank)
Position of rank in ib_nbrs.
Definition m_ibm.fpp.f90:3489
m_ibm::ib_reqs
integer, dimension(:), allocatable, private ib_reqs
Definition m_ibm.fpp.f90:470
m_ibm::s_surface_energy_residual
subroutine s_surface_energy_residual(pres, t_ip, t_s, d, ys_s, r_energy)
Surface energy residual: gas-side conduction balances heterogeneous reaction heat....
Definition m_ibm.fpp.f90:4126
m_ibm::ib_nbrs
integer, dimension(:), allocatable, private ib_nbrs
Definition m_ibm.fpp.f90:470
m_ibm::s_interpolation_weights
subroutine, private s_interpolation_weights(gp, i, j, k, dist, fluid_mask, interp_coeffs)
Fluid-masked image-point weights: bilinear for second-order IBM, else inverse distance squared.
Definition m_ibm.fpp.f90:2356
m_ibm::s_surface_residual
subroutine s_surface_residual(pres, t_ip, t_s, d, ys_ip, ys_s, w_species, k_bath, solve_temperature, flux_scale, energy_scale, r, r_species, omega_s, mdot_s)
Assemble the Newton residual for Ns species, with temperature appended only when it is solved....
Definition m_ibm.fpp.f90:4163
m_ibm::f_log_cosh
pure real(wp) function f_log_cosh(x)
log(cosh(x)) without overflow for large |x|
Definition m_ibm.fpp.f90:2770
m_ibm::ib_send_bufs
real(wp), dimension(:,:), allocatable, private ib_send_bufs
Definition m_ibm.fpp.f90:471
m_ibm::s_prescribed_kinematics
subroutine s_prescribed_kinematics(i, t)
Prescribed kinematics for patch i at time t: hinged flapping (kin_model = 1) or the Eldredge pitch ra...
Definition m_ibm.fpp.f90:2798
m_ibm::ib_markers
type(integer_field), public ib_markers
Definition m_ibm.fpp.f90:392
m_ibm::s_blend_ghost_temperature
subroutine s_blend_ghost_temperature(t_ip, t_s, t_lo, t_hi, t_g)
Ghost state for an immersed surface, from the image-point state and the surface state.
Definition m_ibm.fpp.f90:4231
m_ibm::n_surface_ill_posed
integer n_surface_ill_posed
Ghost points whose levelset distance was not positive, leaving the surface gradient (X_IP - X_s)/d un...
Definition m_ibm.fpp.f90:391
m_ibm::s_check_every_patch_marked
impure subroutine s_check_every_patch_marked()
Abort if any namelist immersed boundary marks no cell anywhere. The owning rank is chosen by the patc...
Definition m_ibm.fpp.f90:3824
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:478
m_ibm::f_needs_centroid_offset
pure logical function f_needs_centroid_offset(patch)
A moving airfoil or STL (geometries 4, 5, 11, 12) moves its centroid to the centre of mass and keeps ...
Definition m_ibm.fpp.f90:3199
m_ibm::s_restore_centroid_offsets
impure subroutine s_restore_centroid_offsets(t_init)
On restart, replace the re-measured centroid offsets with the ones the run was using (restart_data/ib...
Definition m_ibm.fpp.f90:3209
m_ibm::s_solve_surface_linear_system
subroutine s_solve_surface_linear_system(a, b, x, nsolve, success)
Small dense linear solve with partial pivoting for the local surface Newton system.
Definition m_ibm.fpp.f90:4513
m_ibm::s_setup_ib_force_comm
impure subroutine s_setup_ib_force_comm()
Builds the ascending list of distinct neighborhood ranks used by s_communicate_ib_forces and allocate...
Definition m_ibm.fpp.f90:3447
m_ibm::s_wrap_periodic_ibs
subroutine s_wrap_periodic_ibs()
Wrap immersed boundary positions across periodic domain boundaries.
Definition m_ibm.fpp.f90:3328
m_ibm::ib_rec_len
integer, parameter, private ib_rec_len
gbl_patch_id, force(1:3), torque(1:3)
Definition m_ibm.fpp.f90:468
m_ibm::s_compute_ghost_point_pressure
subroutine, private s_compute_ghost_point_pressure(gp, gp_patch_id, rho, pres_ip, pres_gp)
Pressure correction for a moving IB, accounting for the acceleration of the boundary surface....
Definition m_ibm.fpp.f90:1020
m_ibm::ib_recv_bufs
real(wp), dimension(:,:), allocatable, private ib_recv_bufs
Definition m_ibm.fpp.f90:471
m_ibm::centroid_offsets_active
logical centroid_offsets_active
some patch is a moving airfoil or STL, which carries a centroid_offset
Definition m_ibm.fpp.f90:465
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:1129
m_ibm::n_ib_nbrs
integer, private n_ib_nbrs
Definition m_ibm.fpp.f90:469
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:2000
m_ibm::ib_stats
integer, dimension(:,:), allocatable, private ib_stats
Definition m_ibm.fpp.f90:470
m_ibm::s_compute_interpolation_coeffs
subroutine, private s_compute_interpolation_coeffs()
Compute the interpolation coefficients for image points.
Definition m_ibm.fpp.f90:2263
m_model
Binary STL file reader and processor for immersed boundary geometry.
Definition m_model.fpp.f90:357
m_mpi_proxy
MPI halo exchange, domain decomposition, and buffer packing/unpacking for the simulation solver.
Definition m_mpi_proxy.fpp.f90:367
m_patch_geometries
Contains helper functions specific to various patch gemoetries for determining if a grid cell lies in...
Definition m_patch_geometries.fpp.f90:357
m_variables_conversion
Conservative-to-primitive variable conversion, mixture property evaluation, and pressure computation.
Definition m_variables_conversion.fpp.f90:367
m_viscous
Computes viscous stress tensors and diffusive flux contributions for the Navier–Stokes equations.
Definition m_viscous.fpp.f90:368
m_derived_types::ghost_point
Ghost Point for Immersed Boundaries.
Definition m_derived_types.fpp.f90:878
m_derived_types::integer_field
Derived type annexing an integer scalar field (SF).
Definition m_derived_types.fpp.f90:380