MFC
Exascale flow solver
Loading...
Searching...
No Matches
m_bubbles_EL_kernels.fpp.f90
Go to the documentation of this file.
1# 1 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
2!>
3!! @file
4!! @brief Contains module @ref m_bubbles_el_kernels "m_bubbles_EL_kernels"
5
6# 1 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 1
7# 1 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 1
8# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
9# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
10# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
11# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
12# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
13# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
14
15# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
16# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
17# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
18
19# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
20
21# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
22
23# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
24
25# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
26
27# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
28
29# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
30
31# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
32
33# 145 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
34! New line at end of file is required for FYPP
35# 2 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
36# 1 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 1
37# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
38# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
39# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
40# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
41# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
42# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
43
44# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
45# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
46# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
47
48# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
49
50# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
51
52# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
53
54# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
55
56# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
57
58# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
59
60# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
61
62# 145 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
63! New line at end of file is required for FYPP
64# 2 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 2
65
66# 4 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
67# 5 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
68# 6 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
69# 7 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
70# 8 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
71
72# 20 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
73
74# 43 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
75
76# 48 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
77
78# 53 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
79
80# 58 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
81
82# 63 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
83
84# 68 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
85
86# 76 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
87
88# 81 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
89
90# 86 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
91
92# 91 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
93
94# 96 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
95
96# 101 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
97
98# 106 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
99
100# 111 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
101
102# 116 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
103
104# 121 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
105
106# 151 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
107
108# 192 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
109
110# 206 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
111
112# 231 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
113
114# 242 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
115
116# 244 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
117# 255 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
118
119# 284 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
120
121# 294 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
122
123# 304 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
124
125# 313 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
126
127# 330 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
128
129# 340 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
130
131# 347 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
132
133# 353 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
134
135# 359 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
136
137# 365 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
138
139# 371 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
140
141# 377 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
142! New line at end of file is required for FYPP
143# 3 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
144# 1 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 1
145# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
146# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
147# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
148# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
149# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
150# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
151
152# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
153# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
154# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
155
156# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
157
158# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
159
160# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
161
162# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
163
164# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
165
166# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
167
168# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
169
170# 145 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
171! New line at end of file is required for FYPP
172# 2 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 2
173
174# 7 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
175
176# 17 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
177
178# 22 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
179
180# 27 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
181
182# 32 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
183
184# 37 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
185
186# 42 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
187
188# 47 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
189
190# 52 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
191
192# 57 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
193
194# 62 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
195
196# 73 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
197
198# 78 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
199
200# 83 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
201
202# 88 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
203
204# 103 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
205
206# 131 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
207
208# 160 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
209
210# 175 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
211
212# 193 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
213
214# 215 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
215
216# 244 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
217
218# 259 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
219
220# 269 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
221
222# 278 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
223
224# 294 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
225
226# 304 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
227
228# 311 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
229! New line at end of file is required for FYPP
230# 4 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
231
232! GPU parallel region (scalar reductions, maxval/minval)
233# 23 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
234
235! GPU parallel loop over threads (most common GPU macro)
236# 43 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
237
238! Required closing for GPU_PARALLEL_LOOP
239# 55 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
240
241! Mark routine for device compilation
242# 112 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
243
244! Declare device-resident data
245# 130 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
246
247! Inner loop within a GPU parallel region
248# 145 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
249
250! Scoped GPU data region
251# 164 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
252
253! Host code with device pointers (for MPI with GPU buffers)
254# 193 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
255
256! Allocate device memory (unscoped)
257# 207 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
258
259! Free device memory
260# 219 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
261
262! Atomic operation on device
263# 231 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
264
265! End atomic capture block
266# 242 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
267
268! Copy data between host and device
269# 254 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
270
271! Synchronization barrier
272# 266 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
273
274! Import GPU library module (openacc or omp_lib)
275# 275 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
276
277! Emit code only for AMD compiler
278# 282 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
279
280! Emit code for non-Cray compilers
281# 289 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
282
283! Emit code only for Cray compiler
284# 296 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
285
286! Emit code for non-NVIDIA compilers
287# 303 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
288
289# 305 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
290# 306 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
291! New line at end of file is required for FYPP
292# 2 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 2
293
294# 14 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
295
296! Caution: This macro requires the use of a binding script to set CUDA_VISIBLE_DEVICES, such that we have one GPU device per MPI
297! rank. That's because for both cudaMemAdvise (preferred location) and cudaMemPrefetchAsync we use location = device_id = 0. For an
298! example see misc/nvidia_uvm/bind.sh.
299# 57 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
300
301! Allocate and create GPU device memory
302# 77 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
303
304! Free GPU device memory and deallocate
305# 85 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
306
307! Cray-specific GPU pointer setup for vector fields
308# 109 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
309
310! Cray-specific GPU pointer setup for scalar fields
311# 125 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
312
313! Cray-specific GPU pointer setup for acoustic source spatials
314# 150 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
315
316# 156 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
317
318# 163 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
319! New line at end of file is required for FYPP
320# 6 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp" 2
321
322!> @brief Kernel functions (Gaussian, delta) that smear Lagrangian bubble effects onto the Eulerian grid
324
325 use m_mpi_proxy
326
327 implicit none
328
329 ! Cell-centered pressure gradients (precomputed for translational motion)
330 real(wp), allocatable, dimension(:,:,:) :: grad_p_x, grad_p_y, grad_p_z
331
332# 16 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
333#if defined(MFC_OpenACC)
334# 16 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
335!$acc declare create(grad_p_x, grad_p_y, grad_p_z)
336# 16 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
337#elif defined(MFC_OpenMP)
338# 16 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
339!$omp declare target (grad_p_x, grad_p_y, grad_p_z)
340# 16 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
341#endif
342
343 ! Finite-difference coefficients for pressure gradient computation
344 real(wp), allocatable, dimension(:,:) :: fd_coeff_x_pgrad
345 real(wp), allocatable, dimension(:,:) :: fd_coeff_y_pgrad
346 real(wp), allocatable, dimension(:,:) :: fd_coeff_z_pgrad
347
348# 22 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
349#if defined(MFC_OpenACC)
350# 22 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
351!$acc declare create(fd_coeff_x_pgrad, fd_coeff_y_pgrad, fd_coeff_z_pgrad)
352# 22 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
353#elif defined(MFC_OpenMP)
354# 22 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
355!$omp declare target (fd_coeff_x_pgrad, fd_coeff_y_pgrad, fd_coeff_z_pgrad)
356# 22 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
357#endif
358
359 ! Cell list for bubble-to-cell mapping (rebuilt each RK stage before smearing)
360 integer, allocatable, dimension(:,:,:) :: cell_list_start ! (0:m, 0:n, 0:p)
361 integer, allocatable, dimension(:,:,:) :: cell_list_count ! (0:m, 0:n, 0:p)
362 integer, allocatable, dimension(:) :: cell_list_idx ! (1:nBubs_glb) sorted bubble indices
363
364# 28 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
365#if defined(MFC_OpenACC)
366# 28 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
367!$acc declare create(cell_list_start, cell_list_count, cell_list_idx)
368# 28 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
369#elif defined(MFC_OpenMP)
370# 28 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
371!$omp declare target (cell_list_start, cell_list_count, cell_list_idx)
372# 28 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
373#endif
374
375contains
376
377 !> Smear the Lagrangian bubble effects onto the Eulerian grid using the selected kernel
378 subroutine s_smoothfunction(nBubs, lbk_rad, lbk_vel, lbk_s, lbk_pos, updatedvar, kcomp)
379
380 integer, intent(in) :: nBubs
381 real(wp), dimension(1:lag_params%nBubs_glb,1:3,1:2), intent(in) :: lbk_s, lbk_pos
382 real(wp), dimension(1:lag_params%nBubs_glb,1:2), intent(in) :: lbk_rad, lbk_vel
383 type(scalar_field), dimension(:), intent(inout) :: updatedvar
384 type(scalar_field), dimension(:), intent(inout) :: kcomp
385
386 smoothfunc:select case(lag_params%smooth_type)
387 case (1)
388 call s_gaussian(nbubs, lbk_rad, lbk_vel, lbk_s, lbk_pos, updatedvar, kcomp)
389 case (2)
390 call s_deltafunc(nbubs, lbk_rad, lbk_vel, lbk_s, updatedvar, kcomp)
391 end select smoothfunc
392
393 end subroutine s_smoothfunction
394
395 !> Builds a sorted cell list mapping each interior cell (0:m,0:n,0:p) to its resident bubbles. Uses a counting-sort on the host
396 !! (O(nBubs + N_cells)). Must be called before s_gaussian each RK stage.
397 subroutine s_build_cell_list(nBubs, lbk_s)
398
399 integer, intent(in) :: nBubs
400 real(wp), dimension(1:lag_params%nBubs_glb,1:3,1:2), intent(in) :: lbk_s
401 integer :: l, ci, cj, ck, idx
402 real(wp), dimension(3) :: s_coord
403
404 ! Bring current bubble positions to host
405
406
407# 61 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
408#if defined(MFC_OpenACC)
409# 61 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
410!$acc update host(lbk_s)
411# 61 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
412#elif defined(MFC_OpenMP)
413# 61 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
414!$omp target update from(lbk_s)
415# 61 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
416#endif
417
418 ! Pass 1: zero counts and count bubbles per cell
420 do l = 1, nbubs
421 s_coord(1:3) = lbk_s(l,1:3,2)
422 ci = int(s_coord(1))
423 cj = int(s_coord(2))
424 ck = int(s_coord(3))
425 ! Clamp to interior (bubbles should already be in [0:m,0:n,0:p])
426 ci = max(0, min(ci, m))
427 cj = max(0, min(cj, n))
428 ck = max(0, min(ck, p))
429 cell_list_count(ci, cj, ck) = cell_list_count(ci, cj, ck) + 1
430 end do
431
432 ! Prefix sum to compute start indices (1-based into cell_list_idx)
433 idx = 1
434 do ck = 0, p
435 do cj = 0, n
436 do ci = 0, m
437 cell_list_start(ci, cj, ck) = idx
438 idx = idx + cell_list_count(ci, cj, ck)
439 end do
440 end do
441 end do
442
443 ! Pass 2: place bubble indices into cell_list_idx Temporarily reuse cell_list_count as a running offset
445 do l = 1, nbubs
446 s_coord(1:3) = lbk_s(l,1:3,2)
447 ci = int(s_coord(1))
448 cj = int(s_coord(2))
449 ck = int(s_coord(3))
450 ci = max(0, min(ci, m))
451 cj = max(0, min(cj, n))
452 ck = max(0, min(ck, p))
453 cell_list_idx(cell_list_start(ci, cj, ck) + cell_list_count(ci, cj, ck)) = l
454 cell_list_count(ci, cj, ck) = cell_list_count(ci, cj, ck) + 1
455 end do
456
457 ! Send cell list arrays to GPU
458
459# 103 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
460#if defined(MFC_OpenACC)
461# 103 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
462!$acc update device(cell_list_start, cell_list_count, cell_list_idx)
463# 103 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
464#elif defined(MFC_OpenMP)
465# 103 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
466!$omp target update to(cell_list_start, cell_list_count, cell_list_idx)
467# 103 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
468#endif
469
470 end subroutine s_build_cell_list
471
472 !> Cell-centric delta-function smearing using the cell list (no GPU atomics). Each bubble only affects the cell it resides in.
473 !! The outer GPU loop iterates over interior cells and sums contributions from resident bubbles.
474 subroutine s_deltafunc(nBubs, lbk_rad, lbk_vel, lbk_s, updatedvar, kcomp)
475
476 integer, intent(in) :: nBubs
477 real(wp), dimension(1:lag_params%nBubs_glb,1:3,1:2), intent(in) :: lbk_s
478 real(wp), dimension(1:lag_params%nBubs_glb,1:2), intent(in) :: lbk_rad, lbk_vel
479 type(scalar_field), dimension(:), intent(inout) :: updatedvar
480 type(scalar_field), dimension(:), intent(inout) :: kcomp
481 real(wp) :: strength_vel, strength_vol
482 real(wp) :: volpart, Vol
483 real(wp) :: y_kahan, t_kahan
484 integer :: i, j, k, lb, bub_idx
485
486
487# 121 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
488
489# 121 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
490#if defined(MFC_OpenACC)
491# 121 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
492!$acc parallel loop collapse(3) gang vector default(present) private(i, j, k, lb, bub_idx, volpart, Vol, strength_vel, strength_vol, y_kahan, t_kahan)
493# 121 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
494#elif defined(MFC_OpenMP)
495# 121 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
496
497# 121 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
498
499# 121 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
500
501# 121 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
502!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) &
503# 121 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
504!$omp& private(i, j, k, lb, bub_idx, volpart, Vol, strength_vel, strength_vol, y_kahan, t_kahan)
505# 121 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
506#endif
507# 123 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
508 do k = 0, p
509 do j = 0, n
510 do i = 0, m
511 ! Cell volume
512 if (num_dims == 2) then
513 vol = dx(i)*dy(j)*lag_params%charwidth
514 if (cyl_coord) vol = dx(i)*dy(j)*y_cc(j)*2._wp*pi
515 else
516 vol = dx(i)*dy(j)*dz(k)
517 end if
518
519 ! Loop over bubbles in this cell
520
521# 135 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
522#if defined(MFC_OpenACC)
523# 135 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
524!$acc loop seq
525# 135 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
526#elif defined(MFC_OpenMP)
527# 135 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
528
529# 135 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
530#endif
531 do lb = cell_list_start(i, j, k), cell_list_start(i, j, k) + cell_list_count(i, j, k) - 1
532 bub_idx = cell_list_idx(lb)
533
534 volpart = 4._wp/3._wp*pi*lbk_rad(bub_idx, 2)**3._wp
535 strength_vol = volpart
536 strength_vel = 4._wp*pi*lbk_rad(bub_idx, 2)**2._wp*lbk_vel(bub_idx, 2)
537
538 if (lag_params%kahan_summation) then
539 ! Kahan summation for void fraction
540 y_kahan = real(strength_vol/vol, kind=wp) - kcomp(1)%sf(i, j, k)
541 t_kahan = updatedvar(1)%sf(i, j, k) + y_kahan
542 kcomp(1)%sf(i, j, k) = (t_kahan - updatedvar(1)%sf(i, j, k)) - y_kahan
543 updatedvar(1)%sf(i, j, k) = t_kahan
544
545 ! Kahan summation for time derivative of void fraction
546 y_kahan = real(strength_vel/vol, kind=wp) - kcomp(2)%sf(i, j, k)
547 t_kahan = updatedvar(2)%sf(i, j, k) + y_kahan
548 kcomp(2)%sf(i, j, k) = (t_kahan - updatedvar(2)%sf(i, j, k)) - y_kahan
549 updatedvar(2)%sf(i, j, k) = t_kahan
550 else
551 updatedvar(1)%sf(i, j, k) = updatedvar(1)%sf(i, j, k) + real(strength_vol/vol, kind=wp)
552 updatedvar(2)%sf(i, j, k) = updatedvar(2)%sf(i, j, k) + real(strength_vel/vol, kind=wp)
553 end if
554
555 ! Product of two smeared functions
556 if (lag_params%kahan_summation .and. lag_params%cluster_type >= 4) then
557 y_kahan = real((strength_vol*strength_vel)/vol, kind=wp) - kcomp(5)%sf(i, j, k)
558 t_kahan = updatedvar(5)%sf(i, j, k) + y_kahan
559 kcomp(5)%sf(i, j, k) = (t_kahan - updatedvar(5)%sf(i, j, k)) - y_kahan
560 updatedvar(5)%sf(i, j, k) = t_kahan
561 else if (lag_params%cluster_type >= 4) then
562 updatedvar(5)%sf(i, j, k) = updatedvar(5)%sf(i, j, k) + real((strength_vol*strength_vel)/vol, kind=wp)
563 end if
564 end do
565 end do
566 end do
567 end do
568
569# 173 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
570#if defined(MFC_OpenACC)
571# 173 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
572!$acc end parallel loop
573# 173 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
574#elif defined(MFC_OpenMP)
575# 173 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
576
577# 173 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
578!$omp end target teams loop
579# 173 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
580#endif
581
582 end subroutine s_deltafunc
583
584 !> Cell-centric gaussian smearing using the cell list (no GPU atomics). Each grid cell accumulates contributions from nearby
585 !! bubbles looked up via cell_list_start/count/idx.
586 subroutine s_gaussian(nBubs, lbk_rad, lbk_vel, lbk_s, lbk_pos, updatedvar, kcomp)
587
588 integer, intent(in) :: nBubs
589 real(wp), dimension(1:lag_params%nBubs_glb,1:3,1:2), intent(in) :: lbk_s, lbk_pos
590 real(wp), dimension(1:lag_params%nBubs_glb,1:2), intent(in) :: lbk_rad, lbk_vel
591 type(scalar_field), dimension(:), intent(inout) :: updatedvar
592 type(scalar_field), dimension(:), intent(inout) :: kcomp
593 real(wp), dimension(3) :: center, nodecoord, s_coord
594 integer, dimension(3) :: cell, cellijk
595 real(wp) :: stddsv, volpart
596 real(wp) :: strength_vel, strength_vol
597 real(wp) :: func, func2
598 real(wp) :: y_kahan, t_kahan
599 integer :: i, j, k, di, dj, dk, lb, bub_idx
600 integer :: di_beg, di_end, dj_beg, dj_end, dk_beg, dk_end
601 integer :: smear_x_beg, smear_x_end
602 integer :: smear_y_beg, smear_y_end
603 integer :: smear_z_beg, smear_z_end
604
605 ! Extended grid range for smearing (includes buffer cells for MPI communication)
606
607 smear_x_beg = -mapcells - 1
608 smear_x_end = m + mapcells + 1
609 smear_y_beg = merge(-mapcells - 1, 0, n > 0)
610 smear_y_end = merge(n + mapcells + 1, n, n > 0)
611 smear_z_beg = merge(-mapcells - 1, 0, p > 0)
612 smear_z_end = merge(p + mapcells + 1, p, p > 0)
613
614
615# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
616
617# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
618#if defined(MFC_OpenACC)
619# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
620!$acc parallel loop collapse(3) gang vector default(present) &
621# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
622!$acc& private(i, j, k, di, dj, dk, lb, bub_idx, center, nodecoord, s_coord, cell, cellijk, stddsv, volpart, strength_vel, strength_vol, func, func2, y_kahan, t_kahan, di_beg, di_end, dj_beg, dj_end, dk_beg, dk_end) &
623# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
624!$acc& copyin(smear_x_beg, smear_x_end, smear_y_beg, smear_y_end, smear_z_beg, smear_z_end)
625# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
626#elif defined(MFC_OpenMP)
627# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
628
629# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
630
631# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
632
633# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
634!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) &
635# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
636!$omp& private(i, j, k, di, dj, dk, lb, bub_idx, center, nodecoord, s_coord, cell, cellijk, stddsv, volpart, strength_vel, strength_vol, func, func2, y_kahan, t_kahan, di_beg, di_end, dj_beg, dj_end, dk_beg, dk_end) &
637# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
638!$omp& map(to:smear_x_beg, smear_x_end, smear_y_beg, smear_y_end, smear_z_beg, smear_z_end)
639# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
640#endif
641# 210 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
642 do k = smear_z_beg, smear_z_end
643 do j = smear_y_beg, smear_y_end
644 do i = smear_x_beg, smear_x_end
645 cellijk(1) = i
646 cellijk(2) = j
647 cellijk(3) = k
648
649 nodecoord(1) = x_cc(i)
650 nodecoord(2) = y_cc(j)
651 nodecoord(3) = 0._wp
652 if (p > 0) nodecoord(3) = z_cc(k)
653
654 ! Neighbor cell range clamped to interior [0:m, 0:n, 0:p]
655 di_beg = max(i - mapcells, 0)
656 di_end = min(i + mapcells, m)
657 dj_beg = max(j - mapcells, 0)
658 dj_end = min(j + mapcells, n)
659 dk_beg = max(k - mapcells, 0)
660 dk_end = min(k + mapcells, p)
661
662
663# 230 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
664#if defined(MFC_OpenACC)
665# 230 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
666!$acc loop seq
667# 230 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
668#elif defined(MFC_OpenMP)
669# 230 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
670
671# 230 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
672#endif
673 do dk = dk_beg, dk_end
674
675# 232 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
676#if defined(MFC_OpenACC)
677# 232 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
678!$acc loop seq
679# 232 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
680#elif defined(MFC_OpenMP)
681# 232 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
682
683# 232 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
684#endif
685 do dj = dj_beg, dj_end
686
687# 234 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
688#if defined(MFC_OpenACC)
689# 234 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
690!$acc loop seq
691# 234 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
692#elif defined(MFC_OpenMP)
693# 234 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
694
695# 234 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
696#endif
697 do di = di_beg, di_end
698
699# 236 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
700#if defined(MFC_OpenACC)
701# 236 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
702!$acc loop seq
703# 236 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
704#elif defined(MFC_OpenMP)
705# 236 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
706
707# 236 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
708#endif
709 do lb = cell_list_start(di, dj, dk), cell_list_start(di, dj, dk) + cell_list_count(di, dj, dk) - 1
710 bub_idx = cell_list_idx(lb)
711
712 ! Bubble properties
713 volpart = 4._wp/3._wp*pi*lbk_rad(bub_idx, 2)**3._wp
714 s_coord(1:3) = lbk_s(bub_idx,1:3,2)
715 call s_get_cell(s_coord, cell)
716 call s_compute_stddsv(cell, volpart, stddsv)
717
718 strength_vol = volpart
719 strength_vel = 4._wp*pi*lbk_rad(bub_idx, 2)**2._wp*lbk_vel(bub_idx, 2)
720
721 center(1:2) = lbk_pos(bub_idx,1:2,2)
722 center(3) = 0._wp
723 if (p > 0) center(3) = lbk_pos(bub_idx, 3, 2)
724
725 call s_applygaussian(center, cellijk, nodecoord, stddsv, 0._wp, func)
726
727 ! Kahan summation for void fraction
728 y_kahan = real(func*strength_vol, kind=wp) - kcomp(1)%sf(i, j, k)
729 t_kahan = updatedvar(1)%sf(i, j, k) + y_kahan
730 kcomp(1)%sf(i, j, k) = (t_kahan - updatedvar(1)%sf(i, j, k)) - y_kahan
731 updatedvar(1)%sf(i, j, k) = t_kahan
732
733 ! Kahan summation for time derivative of void fraction
734 y_kahan = real(func*strength_vel, kind=wp) - kcomp(2)%sf(i, j, k)
735 t_kahan = updatedvar(2)%sf(i, j, k) + y_kahan
736 kcomp(2)%sf(i, j, k) = (t_kahan - updatedvar(2)%sf(i, j, k)) - y_kahan
737 updatedvar(2)%sf(i, j, k) = t_kahan
738
739 if (lag_params%cluster_type >= 4) then
740 call s_applygaussian(center, cellijk, nodecoord, stddsv, 1._wp, func2)
741 y_kahan = real(func2*strength_vol*strength_vel, kind=wp) - kcomp(5)%sf(i, j, k)
742 t_kahan = updatedvar(5)%sf(i, j, k) + y_kahan
743 kcomp(5)%sf(i, j, k) = (t_kahan - updatedvar(5)%sf(i, j, k)) - y_kahan
744 updatedvar(5)%sf(i, j, k) = t_kahan
745 end if
746 end do
747 end do
748 end do
749 end do
750 end do
751 end do
752 end do
753
754# 281 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
755#if defined(MFC_OpenACC)
756# 281 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
757!$acc end parallel loop
758# 281 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
759#elif defined(MFC_OpenMP)
760# 281 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
761
762# 281 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
763!$omp end target teams loop
764# 281 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
765#endif
766
767 end subroutine s_gaussian
768
769 !> Evaluate the Gaussian kernel at a grid node for a given bubble center
770 subroutine s_applygaussian(center, cellaux, nodecoord, stddsv, strength_idx, func)
771
772
773# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
774#ifdef _CRAYFTN
775# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
776#if MFC_OpenACC
777# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
778!$acc routine seq
779# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
780#elif MFC_OpenMP
781# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
782
783# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
784
785# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
786!$omp declare target device_type(any)
787# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
788#else
789# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
790!DIR$ INLINEALWAYS s_applygaussian
791# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
792#endif
793# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
794#elif MFC_OpenACC
795# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
796!$acc routine seq
797# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
798#elif MFC_OpenMP
799# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
800
801# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
802
803# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
804!$omp declare target device_type(any)
805# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
806#endif
807
808 real(wp), dimension(3), intent(in) :: center
809 integer, dimension(3), intent(in) :: cellaux
810 real(wp), dimension(3), intent(in) :: nodecoord
811 real(wp), intent(in) :: stddsv
812 real(wp), intent(in) :: strength_idx
813 real(wp), intent(out) :: func
814 integer :: i
815 real(wp) :: distance
816 real(wp) :: theta, dtheta, L2, dzp, Lz2, zc
817 real(wp) :: Nr, Nr_count
818
819 distance = sqrt((center(1) - nodecoord(1))**2._wp + (center(2) - nodecoord(2))**2._wp + (center(3) - nodecoord(3))**2._wp)
820
821 if (num_dims == 3) then
822 !> 3D gaussian function
823 func = exp(-0.5_wp*(distance/stddsv)**2._wp)/(sqrt(2._wp*pi)*stddsv)**3._wp
824 else
825 if (cyl_coord) then
826 !> 2D cylindrical function:
827 ! We smear particles in the azimuthal direction for given r
828 theta = 0._wp
829 nr = ceiling(2._wp*pi*nodecoord(2)/(y_cb(cellaux(2)) - y_cb(cellaux(2) - 1)))
830 dtheta = 2._wp*pi/nr
831 l2 = center(2)**2._wp + nodecoord(2)**2._wp - 2._wp*center(2)*nodecoord(2)*cos(theta)
832 distance = sqrt((center(1) - nodecoord(1))**2._wp + l2)
833 ! Factor 2._wp is for symmetry (upper half of the 2D field (+r) is considered)
834 func = dtheta/2._wp/pi*exp(-0.5_wp*(distance/stddsv)**2._wp)/(sqrt(2._wp*pi)*stddsv)**3._wp
835 nr_count = 0._wp
836 do while (nr_count < nr - 1._wp)
837 nr_count = nr_count + 1._wp
838 theta = nr_count*dtheta
839 ! trigonometric relation
840 l2 = center(2)**2._wp + nodecoord(2)**2._wp - 2._wp*center(2)*nodecoord(2)*cos(theta)
841 distance = sqrt((center(1) - nodecoord(1))**2._wp + l2)
842 ! nodecoord(2)*dtheta is the azimuthal width of the cell
843 func = func + dtheta/2._wp/pi*exp(-0.5_wp*(distance/stddsv)**2._wp)/(sqrt(2._wp*pi)*stddsv) &
844 & **(3._wp*(strength_idx + 1._wp))
845 end do
846 else
847 !> 2D cartesian function: Equation (48) from Maeda and Colonius 2018
848 ! We smear particles considering a virtual depth (lag_params%charwidth) with lag_params%charNz cells
849 dzp = (lag_params%charwidth/(lag_params%charNz + 1._wp))
850
851 func = 0._wp
852 do i = 0, lag_params%charNz
853 zc = (-lag_params%charwidth/2._wp + dzp*(0.5_wp + i)) ! Center of virtual cell i in z-direction
854 lz2 = (center(3) - zc)**2._wp
855 distance = sqrt((center(1) - nodecoord(1))**2._wp + (center(2) - nodecoord(2))**2._wp + lz2)
856 func = func + dzp/lag_params%charwidth*exp(-0.5_wp*(distance/stddsv)**2._wp)/(sqrt(2._wp*pi)*stddsv)**3._wp
857 end do
858 end if
859 end if
860
861 end subroutine s_applygaussian
862
863 !> Check if the current cell is outside the computational domain including ghost cells
864 subroutine s_check_celloutside(cellaux, celloutside)
865
866
867# 348 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
868#ifdef _CRAYFTN
869# 348 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
870#if MFC_OpenACC
871# 348 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
872!$acc routine seq
873# 348 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
874#elif MFC_OpenMP
875# 348 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
876
877# 348 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
878
879# 348 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
880!$omp declare target device_type(any)
881# 348 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
882#else
883# 348 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
884!DIR$ INLINEALWAYS s_check_celloutside
885# 348 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
886#endif
887# 348 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
888#elif MFC_OpenACC
889# 348 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
890!$acc routine seq
891# 348 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
892#elif MFC_OpenMP
893# 348 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
894
895# 348 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
896
897# 348 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
898!$omp declare target device_type(any)
899# 348 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
900#endif
901
902 integer, dimension(3), intent(inout) :: cellaux
903 logical, intent(out) :: celloutside
904
905 celloutside = .false.
906
907 if (num_dims == 2) then
908 if ((cellaux(1) < -buff_size) .or. (cellaux(2) < -buff_size)) then
909 celloutside = .true.
910 end if
911 if (cyl_coord .and. y_cc(cellaux(2)) < 0._wp) then
912 celloutside = .true.
913 end if
914 if ((cellaux(2) > n + buff_size) .or. (cellaux(1) > m + buff_size)) then
915 celloutside = .true.
916 end if
917 else
918 if ((cellaux(3) < -buff_size) .or. (cellaux(1) < -buff_size) .or. (cellaux(2) < -buff_size)) then
919 celloutside = .true.
920 end if
921
922 if ((cellaux(3) > p + buff_size) .or. (cellaux(2) > n + buff_size) .or. (cellaux(1) > m + buff_size)) then
923 celloutside = .true.
924 end if
925 end if
926
927 end subroutine s_check_celloutside
928
929 !> Relocate cells that intersect a symmetric boundary
930 subroutine s_shift_cell_symmetric_bc(cellaux, cell)
931
932
933# 380 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
934#ifdef _CRAYFTN
935# 380 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
936#if MFC_OpenACC
937# 380 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
938!$acc routine seq
939# 380 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
940#elif MFC_OpenMP
941# 380 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
942
943# 380 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
944
945# 380 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
946!$omp declare target device_type(any)
947# 380 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
948#else
949# 380 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
950!DIR$ INLINEALWAYS s_shift_cell_symmetric_bc
951# 380 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
952#endif
953# 380 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
954#elif MFC_OpenACC
955# 380 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
956!$acc routine seq
957# 380 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
958#elif MFC_OpenMP
959# 380 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
960
961# 380 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
962
963# 380 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
964!$omp declare target device_type(any)
965# 380 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
966#endif
967
968 integer, dimension(3), intent(inout) :: cellaux
969 integer, dimension(3), intent(in) :: cell
970
971 ! x-dir
972 if (bc_x%beg == bc_reflective .and. (cell(1) <= mapcells - 1)) then
973 cellaux(1) = abs(cellaux(1)) - 1
974 end if
975 if (bc_x%end == bc_reflective .and. (cell(1) >= m + 1 - mapcells)) then
976 cellaux(1) = cellaux(1) - (2*(cellaux(1) - m) - 1)
977 end if
978
979 ! y-dir
980 if (bc_y%beg == bc_reflective .and. (cell(2) <= mapcells - 1)) then
981 cellaux(2) = abs(cellaux(2)) - 1
982 end if
983 if (bc_y%end == bc_reflective .and. (cell(2) >= n + 1 - mapcells)) then
984 cellaux(2) = cellaux(2) - (2*(cellaux(2) - n) - 1)
985 end if
986
987 if (p > 0) then
988 ! z-dir
989 if (bc_z%beg == bc_reflective .and. (cell(3) <= mapcells - 1)) then
990 cellaux(3) = abs(cellaux(3)) - 1
991 end if
992 if (bc_z%end == bc_reflective .and. (cell(3) >= p + 1 - mapcells)) then
993 cellaux(3) = cellaux(3) - (2*(cellaux(3) - p) - 1)
994 end if
995 end if
996
997 end subroutine s_shift_cell_symmetric_bc
998
999 !> Calculates the standard deviation of the bubble being smeared in the Eulerian framework.
1000 subroutine s_compute_stddsv(cell, volpart, stddsv)
1001
1002
1003# 416 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1004#ifdef _CRAYFTN
1005# 416 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1006#if MFC_OpenACC
1007# 416 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1008!$acc routine seq
1009# 416 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1010#elif MFC_OpenMP
1011# 416 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1012
1013# 416 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1014
1015# 416 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1016!$omp declare target device_type(any)
1017# 416 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1018#else
1019# 416 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1020!DIR$ INLINEALWAYS s_compute_stddsv
1021# 416 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1022#endif
1023# 416 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1024#elif MFC_OpenACC
1025# 416 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1026!$acc routine seq
1027# 416 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1028#elif MFC_OpenMP
1029# 416 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1030
1031# 416 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1032
1033# 416 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1034!$omp declare target device_type(any)
1035# 416 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1036#endif
1037
1038 integer, dimension(3), intent(in) :: cell
1039 real(wp), intent(in) :: volpart
1040 real(wp), intent(out) :: stddsv
1041 real(wp) :: chardist, charvol
1042 real(wp) :: rad
1043
1044 !> Compute characteristic distance
1045 chardist = sqrt(dx(cell(1))*dy(cell(2)))
1046 if (p > 0) chardist = (dx(cell(1))*dy(cell(2))*dz(cell(3)))**(1._wp/3._wp)
1047
1048 !> Compute characteristic volume
1049 if (p > 0) then
1050 charvol = dx(cell(1))*dy(cell(2))*dz(cell(3))
1051 else
1052 if (cyl_coord) then
1053 charvol = dx(cell(1))*dy(cell(2))*y_cc(cell(2))*2._wp*pi
1054 else
1055 charvol = dx(cell(1))*dy(cell(2))*lag_params%charwidth
1056 end if
1057 end if
1058
1059 !> Compute Standard deviaton
1060 if ((volpart/charvol) > 0.5_wp*lag_params%valmaxvoid .or. (lag_params%smooth_type == 1)) then
1061 rad = (3._wp*volpart/(4._wp*pi))**(1._wp/3._wp)
1062 stddsv = 1._wp*lag_params%epsilonb*max(chardist, rad)
1063 else
1064 stddsv = 0._wp
1065 end if
1066
1067 end subroutine s_compute_stddsv
1068
1069 !> Compute the characteristic cell volume
1070 subroutine s_get_char_vol(cellx, celly, cellz, Charvol)
1071
1072
1073# 452 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1074#ifdef _CRAYFTN
1075# 452 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1076#if MFC_OpenACC
1077# 452 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1078!$acc routine seq
1079# 452 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1080#elif MFC_OpenMP
1081# 452 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1082
1083# 452 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1084
1085# 452 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1086!$omp declare target device_type(any)
1087# 452 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1088#else
1089# 452 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1090!DIR$ INLINEALWAYS s_get_char_vol
1091# 452 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1092#endif
1093# 452 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1094#elif MFC_OpenACC
1095# 452 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1096!$acc routine seq
1097# 452 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1098#elif MFC_OpenMP
1099# 452 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1100
1101# 452 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1102
1103# 452 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1104!$omp declare target device_type(any)
1105# 452 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1106#endif
1107
1108 integer, intent(in) :: cellx, celly, cellz
1109 real(wp), intent(out) :: Charvol
1110
1111 if (p > 0) then
1112 charvol = dx(cellx)*dy(celly)*dz(cellz)
1113 else
1114 if (cyl_coord) then
1115 charvol = dx(cellx)*dy(celly)*y_cc(celly)*2._wp*pi
1116 else
1117 charvol = dx(cellx)*dy(celly)*lag_params%charwidth
1118 end if
1119 end if
1120
1121 end subroutine s_get_char_vol
1122
1123 !> Convert bubble computational coordinates from real to integer cell indices
1124 subroutine s_get_cell(s_cell, get_cell)
1125
1126
1127# 472 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1128#ifdef _CRAYFTN
1129# 472 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1130#if MFC_OpenACC
1131# 472 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1132!$acc routine seq
1133# 472 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1134#elif MFC_OpenMP
1135# 472 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1136
1137# 472 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1138
1139# 472 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1140!$omp declare target device_type(any)
1141# 472 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1142#else
1143# 472 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1144!DIR$ INLINEALWAYS s_get_cell
1145# 472 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1146#endif
1147# 472 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1148#elif MFC_OpenACC
1149# 472 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1150!$acc routine seq
1151# 472 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1152#elif MFC_OpenMP
1153# 472 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1154
1155# 472 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1156
1157# 472 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1158!$omp declare target device_type(any)
1159# 472 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1160#endif
1161
1162 real(wp), dimension(3), intent(in) :: s_cell
1163 integer, dimension(3), intent(out) :: get_cell
1164 integer :: i
1165
1166 get_cell(:) = int(s_cell(:))
1167 do i = 1, num_dims
1168 if (s_cell(i) < 0._wp) get_cell(i) = get_cell(i) - 1
1169 end do
1170
1171 end subroutine s_get_cell
1172
1173 !> Precompute cell-centered pressure gradients (dp/dx, dp/dy, dp/dz)
1174 subroutine s_compute_pressure_gradients(q_prim_vf)
1175
1176 type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
1177 integer :: i, j, k, r
1178
1179 ! dp/dx at all cell centers
1180
1181
1182# 493 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1183
1184# 493 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1185#if defined(MFC_OpenACC)
1186# 493 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1187!$acc parallel loop collapse(3) gang vector default(present) private(i, j, k, r)
1188# 493 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1189#elif defined(MFC_OpenMP)
1190# 493 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1191
1192# 493 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1193
1194# 493 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1195
1196# 493 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1197!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, j, k, r)
1198# 493 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1199#endif
1200 do k = 0, p
1201 do j = 0, n
1202 do i = 0, m
1203 grad_p_x(i, j, k) = 0._wp
1204
1205# 498 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1206#if defined(MFC_OpenACC)
1207# 498 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1208!$acc loop seq
1209# 498 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1210#elif defined(MFC_OpenMP)
1211# 498 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1212
1213# 498 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1214#endif
1215 do r = -fd_number, fd_number
1216 grad_p_x(i, j, k) = grad_p_x(i, j, k) + q_prim_vf(eqn_idx%E)%sf(i + r, j, k)*fd_coeff_x_pgrad(r, i)
1217 end do
1218 end do
1219 end do
1220 end do
1221
1222# 505 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1223#if defined(MFC_OpenACC)
1224# 505 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1225!$acc end parallel loop
1226# 505 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1227#elif defined(MFC_OpenMP)
1228# 505 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1229
1230# 505 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1231!$omp end target teams loop
1232# 505 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1233#endif
1234
1235 ! dp/dy at all cell centers
1236 if (n > 0) then
1237
1238# 509 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1239
1240# 509 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1241#if defined(MFC_OpenACC)
1242# 509 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1243!$acc parallel loop collapse(3) gang vector default(present) private(i, j, k, r)
1244# 509 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1245#elif defined(MFC_OpenMP)
1246# 509 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1247
1248# 509 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1249
1250# 509 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1251
1252# 509 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1253!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, j, k, r)
1254# 509 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1255#endif
1256 do k = 0, p
1257 do j = 0, n
1258 do i = 0, m
1259 grad_p_y(i, j, k) = 0._wp
1260
1261# 514 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1262#if defined(MFC_OpenACC)
1263# 514 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1264!$acc loop seq
1265# 514 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1266#elif defined(MFC_OpenMP)
1267# 514 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1268
1269# 514 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1270#endif
1271 do r = -fd_number, fd_number
1272 grad_p_y(i, j, k) = grad_p_y(i, j, k) + q_prim_vf(eqn_idx%E)%sf(i, j + r, k)*fd_coeff_y_pgrad(r, j)
1273 end do
1274 end do
1275 end do
1276 end do
1277
1278# 521 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1279#if defined(MFC_OpenACC)
1280# 521 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1281!$acc end parallel loop
1282# 521 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1283#elif defined(MFC_OpenMP)
1284# 521 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1285
1286# 521 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1287!$omp end target teams loop
1288# 521 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1289#endif
1290 end if
1291
1292 ! dp/dz at all cell centers
1293 if (p > 0) then
1294
1295# 526 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1296
1297# 526 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1298#if defined(MFC_OpenACC)
1299# 526 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1300!$acc parallel loop collapse(3) gang vector default(present) private(i, j, k, r)
1301# 526 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1302#elif defined(MFC_OpenMP)
1303# 526 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1304
1305# 526 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1306
1307# 526 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1308
1309# 526 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1310!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, j, k, r)
1311# 526 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1312#endif
1313 do k = 0, p
1314 do j = 0, n
1315 do i = 0, m
1316 grad_p_z(i, j, k) = 0._wp
1317
1318# 531 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1319#if defined(MFC_OpenACC)
1320# 531 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1321!$acc loop seq
1322# 531 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1323#elif defined(MFC_OpenMP)
1324# 531 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1325
1326# 531 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1327#endif
1328 do r = -fd_number, fd_number
1329 grad_p_z(i, j, k) = grad_p_z(i, j, k) + q_prim_vf(eqn_idx%E)%sf(i, j, k + r)*fd_coeff_z_pgrad(r, k)
1330 end do
1331 end do
1332 end do
1333 end do
1334
1335# 538 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1336#if defined(MFC_OpenACC)
1337# 538 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1338!$acc end parallel loop
1339# 538 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1340#elif defined(MFC_OpenMP)
1341# 538 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1342
1343# 538 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1344!$omp end target teams loop
1345# 538 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1346#endif
1347 end if
1348
1349 end subroutine s_compute_pressure_gradients
1350
1351 !! Interpolate the velocity of Eulerian field at the position of the bubble.
1352 function f_interpolate_velocity(pos, cell, i, q_prim_vf) result(v)
1353
1354
1355# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1356#ifdef _CRAYFTN
1357# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1358#if MFC_OpenACC
1359# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1360!$acc routine seq
1361# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1362#elif MFC_OpenMP
1363# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1364
1365# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1366
1367# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1368!$omp declare target device_type(any)
1369# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1370#else
1371# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1372!DIR$ NOINLINE f_interpolate_velocity
1373# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1374#endif
1375# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1376#elif MFC_OpenACC
1377# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1378!$acc routine seq
1379# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1380#elif MFC_OpenMP
1381# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1382
1383# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1384
1385# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1386!$omp declare target device_type(any)
1387# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1388#endif
1389
1390 real(wp), intent(in) :: pos
1391 integer, dimension(3), intent(in) :: cell
1392 integer, intent(in) :: i
1393 type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
1394 real(wp) :: v
1395 real(wp), dimension(5) :: xi, eta, l
1396
1397 if (fd_order == 2) then
1398 if (i == 1) then
1399 xi(1) = x_cc(cell(1) - 1)
1400 eta(1) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1) - 1, cell(2), cell(3))
1401 xi(2) = x_cc(cell(1))
1402 eta(2) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1), cell(2), cell(3))
1403 xi(3) = x_cc(cell(1) + 1)
1404 eta(3) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1) + 1, cell(2), cell(3))
1405 else if (i == 2) then
1406 xi(1) = y_cc(cell(2) - 1)
1407 eta(1) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2) - 1, cell(3))
1408 xi(2) = y_cc(cell(2))
1409 eta(2) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2), cell(3))
1410 xi(3) = y_cc(cell(2) + 1)
1411 eta(3) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2) + 1, cell(3))
1412 else if (i == 3) then
1413 xi(1) = z_cc(cell(3) - 1)
1414 eta(1) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3) - 1)
1415 xi(2) = z_cc(cell(3))
1416 eta(2) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3))
1417 xi(3) = z_cc(cell(3) + 1)
1418 eta(3) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3) + 1)
1419 end if
1420
1421 l(1) = ((pos - xi(2))*(pos - xi(3)))/((xi(1) - xi(2))*(xi(1) - xi(3)))
1422 l(2) = ((pos - xi(1))*(pos - xi(3)))/((xi(2) - xi(1))*(xi(2) - xi(3)))
1423 l(3) = ((pos - xi(1))*(pos - xi(2)))/((xi(3) - xi(1))*(xi(3) - xi(2)))
1424
1425 v = l(1)*eta(1) + l(2)*eta(2) + l(3)*eta(3)
1426 else if (fd_order == 4) then
1427 if (i == 1) then
1428 xi(1) = x_cc(cell(1) - 2)
1429 eta(1) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1) - 2, cell(2), cell(3))
1430 xi(2) = x_cc(cell(1) - 1)
1431 eta(2) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1) - 1, cell(2), cell(3))
1432 xi(3) = x_cc(cell(1))
1433 eta(3) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1), cell(2), cell(3))
1434 xi(4) = x_cc(cell(1) + 1)
1435 eta(4) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1) + 1, cell(2), cell(3))
1436 xi(5) = x_cc(cell(1) + 2)
1437 eta(5) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1) + 2, cell(2), cell(3))
1438 else if (i == 2) then
1439 xi(1) = y_cc(cell(2) - 2)
1440 eta(1) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2) - 2, cell(3))
1441 xi(2) = y_cc(cell(2) - 1)
1442 eta(2) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2) - 1, cell(3))
1443 xi(3) = y_cc(cell(2))
1444 eta(3) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2), cell(3))
1445 xi(4) = y_cc(cell(2) + 1)
1446 eta(4) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2) + 1, cell(3))
1447 xi(5) = y_cc(cell(2) + 2)
1448 eta(5) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2) + 2, cell(3))
1449 else if (i == 3) then
1450 xi(1) = z_cc(cell(3) - 2)
1451 eta(1) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3) - 2)
1452 xi(2) = z_cc(cell(3) - 1)
1453 eta(2) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3) - 1)
1454 xi(3) = z_cc(cell(3))
1455 eta(3) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3))
1456 xi(4) = z_cc(cell(3) + 1)
1457 eta(4) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3) + 1)
1458 xi(5) = z_cc(cell(3) + 2)
1459 eta(5) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3) + 2)
1460 end if
1461
1462 l(1) = ((pos - xi(2))*(pos - xi(3))*(pos - xi(4))*(pos - xi(5)))/((xi(1) - xi(2))*(xi(1) - xi(3))*(xi(1) - xi(4)) &
1463 & *(xi(1) - xi(5)))
1464 l(2) = ((pos - xi(1))*(pos - xi(3))*(pos - xi(4))*(pos - xi(5)))/((xi(2) - xi(1))*(xi(2) - xi(3))*(xi(2) - xi(4)) &
1465 & *(xi(2) - xi(5)))
1466 l(3) = ((pos - xi(1))*(pos - xi(2))*(pos - xi(4))*(pos - xi(5)))/((xi(3) - xi(1))*(xi(3) - xi(2))*(xi(3) - xi(4)) &
1467 & *(xi(3) - xi(5)))
1468 l(4) = ((pos - xi(1))*(pos - xi(2))*(pos - xi(3))*(pos - xi(5)))/((xi(4) - xi(1))*(xi(4) - xi(2))*(xi(4) - xi(3)) &
1469 & *(xi(4) - xi(5)))
1470 l(5) = ((pos - xi(1))*(pos - xi(2))*(pos - xi(3))*(pos - xi(4)))/((xi(5) - xi(1))*(xi(5) - xi(2))*(xi(5) - xi(3)) &
1471 & *(xi(5) - xi(4)))
1472
1473 v = l(1)*eta(1) + l(2)*eta(2) + l(3)*eta(3) + l(4)*eta(4) + l(5)*eta(5)
1474 end if
1475
1476 end function f_interpolate_velocity
1477
1478 !! Calculate the force on a bubble based on the pressure gradient, velocity, and drag model.
1479 !! @param mg Mass of the gas in the bubble
1480 !! @param mv Mass of the liquid in the bubble
1481 !! @param Re Reynolds number
1482 !! @param i Direction of the velocity (1: x, 2: y, 3: z)
1483 function f_get_bubble_force(pos, rad, rdot, vel, mg, mv, Re, rho, cell, i, q_prim_vf) result(force)
1484
1485
1486# 643 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1487#if MFC_OpenACC
1488# 643 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1489!$acc routine seq
1490# 643 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1491#elif MFC_OpenMP
1492# 643 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1493
1494# 643 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1495
1496# 643 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1497!$omp declare target device_type(any)
1498# 643 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1499#endif
1500 real(wp), intent(in) :: pos, rad, rdot, mg, mv, re, rho, vel
1501 integer, dimension(3), intent(in) :: cell
1502 integer, intent(in) :: i
1503 type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
1504 real(wp) :: dp, vol, force
1505 real(wp) :: v_rel
1506
1507 if (fd_order > 1) then
1508 v_rel = vel - f_interpolate_velocity(pos, cell, i, q_prim_vf)
1509 else
1510 v_rel = vel - q_prim_vf(eqn_idx%mom%beg + i - 1)%sf(cell(1), cell(2), cell(3))
1511 end if
1512
1513 force = 0._wp
1514
1515 if (lag_params%drag_model == 1) then ! Free slip Stokes drag
1516 force = force - (4._wp*pi*rad*v_rel)/re
1517 else if (lag_params%drag_model == 2) then ! No slip Stokes drag
1518 force = force - (6._wp*pi*rad*v_rel)/re
1519 else if (lag_params%drag_model == 3) then ! Levich drag
1520 force = force - (12._wp*pi*rad*v_rel)/re
1521 end if
1522
1523 if (lag_pressure_force) then
1524 ! Use precomputed cell-centered pressure gradients
1525 if (i == 1) then
1526 dp = grad_p_x(cell(1), cell(2), cell(3))
1527 else if (i == 2) then
1528 dp = grad_p_y(cell(1), cell(2), cell(3))
1529 else if (i == 3) then
1530 dp = grad_p_z(cell(1), cell(2), cell(3))
1531 end if
1532
1533 vol = (4._wp/3._wp)*pi*(rad**3._wp)
1534 force = force - vol*dp
1535 end if
1536
1537 if (lag_params%gravity_force) then
1538 force = force + (mg + mv)*accel_bf(i)
1539 end if
1540
1541 end function f_get_bubble_force
1542
1543end module m_bubbles_el_kernels
integer, intent(in) l
Kernel functions (Gaussian, delta) that smear Lagrangian bubble effects onto the Eulerian grid.
subroutine s_build_cell_list(nbubs, lbk_s)
Builds a sorted cell list mapping each interior cell (0:m,0:n,0:p) to its resident bubbles....
real(wp), dimension(:,:), allocatable fd_coeff_z_pgrad
subroutine s_deltafunc(nbubs, lbk_rad, lbk_vel, lbk_s, updatedvar, kcomp)
Cell-centric delta-function smearing using the cell list (no GPU atomics). Each bubble only affects t...
subroutine s_applygaussian(center, cellaux, nodecoord, stddsv, strength_idx, func)
Evaluate the Gaussian kernel at a grid node for a given bubble center.
subroutine s_check_celloutside(cellaux, celloutside)
Check if the current cell is outside the computational domain including ghost cells.
real(wp), dimension(:,:,:), allocatable grad_p_y
subroutine s_compute_pressure_gradients(q_prim_vf)
Precompute cell-centered pressure gradients (dp/dx, dp/dy, dp/dz).
real(wp) function f_get_bubble_force(pos, rad, rdot, vel, mg, mv, re, rho, cell, i, q_prim_vf)
integer, dimension(:), allocatable cell_list_idx
real(wp) function f_interpolate_velocity(pos, cell, i, q_prim_vf)
subroutine s_shift_cell_symmetric_bc(cellaux, cell)
Relocate cells that intersect a symmetric boundary.
integer, dimension(:,:,:), allocatable cell_list_start
subroutine s_compute_stddsv(cell, volpart, stddsv)
Calculates the standard deviation of the bubble being smeared in the Eulerian framework.
real(wp), dimension(:,:,:), allocatable grad_p_x
real(wp), dimension(:,:), allocatable fd_coeff_y_pgrad
real(wp), dimension(:,:,:), allocatable grad_p_z
subroutine s_get_char_vol(cellx, celly, cellz, charvol)
Compute the characteristic cell volume.
subroutine s_smoothfunction(nbubs, lbk_rad, lbk_vel, lbk_s, lbk_pos, updatedvar, kcomp)
Smear the Lagrangian bubble effects onto the Eulerian grid using the selected kernel.
subroutine s_gaussian(nbubs, lbk_rad, lbk_vel, lbk_s, lbk_pos, updatedvar, kcomp)
Cell-centric gaussian smearing using the cell list (no GPU atomics). Each grid cell accumulates contr...
subroutine s_get_cell(s_cell, get_cell)
Convert bubble computational coordinates from real to integer cell indices.
integer, dimension(:,:,:), allocatable cell_list_count
real(wp), dimension(:,:), allocatable fd_coeff_x_pgrad
MPI halo exchange, domain decomposition, and buffer packing/unpacking for the simulation solver.