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