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# 167 "/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# 167 "/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# 167 "/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# 55 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
300
301! Allocate and create GPU device memory
302# 75 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
303
304! Free GPU device memory and deallocate
305# 83 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
306
307! Cray-specific GPU pointer setup for vector fields
308# 107 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
309
310! Cray-specific GPU pointer setup for scalar fields
311# 123 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
312
313! Cray-specific GPU pointer setup for acoustic source spatials
314# 148 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
315
316# 154 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
317
318# 161 "/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) private(i, j, k, di, dj, dk, lb, bub_idx, center, nodecoord, s_coord, cell, cellijk, stddsv, volpart, strength_vel, strength_vol, func, &
621# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
622!$acc& func2, y_kahan, t_kahan, di_beg, di_end, dj_beg, dj_end, dk_beg, dk_end) copyin(smear_x_beg, smear_x_end, smear_y_beg, smear_y_end, smear_z_beg, smear_z_end)
623# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
624#elif defined(MFC_OpenMP)
625# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
626
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!$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, di, dj, &
633# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
634!$omp& 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) &
635# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
636!$omp& map(to:smear_x_beg, smear_x_end, smear_y_beg, smear_y_end, smear_z_beg, smear_z_end)
637# 207 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
638#endif
639# 210 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
640 do k = smear_z_beg, smear_z_end
641 do j = smear_y_beg, smear_y_end
642 do i = smear_x_beg, smear_x_end
643 cellijk(1) = i
644 cellijk(2) = j
645 cellijk(3) = k
646
647 nodecoord(1) = x_cc(i)
648 nodecoord(2) = y_cc(j)
649 nodecoord(3) = 0._wp
650 if (p > 0) nodecoord(3) = z_cc(k)
651
652 ! Neighbor cell range clamped to interior [0:m, 0:n, 0:p]
653 di_beg = max(i - mapcells, 0)
654 di_end = min(i + mapcells, m)
655 dj_beg = max(j - mapcells, 0)
656 dj_end = min(j + mapcells, n)
657 dk_beg = max(k - mapcells, 0)
658 dk_end = min(k + mapcells, p)
659
660
661# 230 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
662#if defined(MFC_OpenACC)
663# 230 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
664!$acc loop seq
665# 230 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
666#elif defined(MFC_OpenMP)
667# 230 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
668
669# 230 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
670#endif
671 do dk = dk_beg, dk_end
672
673# 232 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
674#if defined(MFC_OpenACC)
675# 232 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
676!$acc loop seq
677# 232 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
678#elif defined(MFC_OpenMP)
679# 232 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
680
681# 232 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
682#endif
683 do dj = dj_beg, dj_end
684
685# 234 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
686#if defined(MFC_OpenACC)
687# 234 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
688!$acc loop seq
689# 234 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
690#elif defined(MFC_OpenMP)
691# 234 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
692
693# 234 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
694#endif
695 do di = di_beg, di_end
696
697# 236 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
698#if defined(MFC_OpenACC)
699# 236 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
700!$acc loop seq
701# 236 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
702#elif defined(MFC_OpenMP)
703# 236 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
704
705# 236 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
706#endif
707 do lb = cell_list_start(di, dj, dk), cell_list_start(di, dj, dk) + cell_list_count(di, dj, dk) - 1
708 bub_idx = cell_list_idx(lb)
709
710 ! Bubble properties
711 volpart = 4._wp/3._wp*pi*lbk_rad(bub_idx, 2)**3._wp
712 s_coord(1:3) = lbk_s(bub_idx,1:3,2)
713 call s_get_cell(s_coord, cell)
714 call s_compute_stddsv(cell, volpart, stddsv)
715
716 strength_vol = volpart
717 strength_vel = 4._wp*pi*lbk_rad(bub_idx, 2)**2._wp*lbk_vel(bub_idx, 2)
718
719 center(1:2) = lbk_pos(bub_idx,1:2,2)
720 center(3) = 0._wp
721 if (p > 0) center(3) = lbk_pos(bub_idx, 3, 2)
722
723 call s_applygaussian(center, cellijk, nodecoord, stddsv, 0._wp, func)
724
725 ! Kahan summation for void fraction
726 y_kahan = real(func*strength_vol, kind=wp) - kcomp(1)%sf(i, j, k)
727 t_kahan = updatedvar(1)%sf(i, j, k) + y_kahan
728 kcomp(1)%sf(i, j, k) = (t_kahan - updatedvar(1)%sf(i, j, k)) - y_kahan
729 updatedvar(1)%sf(i, j, k) = t_kahan
730
731 ! Kahan summation for time derivative of void fraction
732 y_kahan = real(func*strength_vel, kind=wp) - kcomp(2)%sf(i, j, k)
733 t_kahan = updatedvar(2)%sf(i, j, k) + y_kahan
734 kcomp(2)%sf(i, j, k) = (t_kahan - updatedvar(2)%sf(i, j, k)) - y_kahan
735 updatedvar(2)%sf(i, j, k) = t_kahan
736
737 if (lag_params%cluster_type >= 4) then
738 call s_applygaussian(center, cellijk, nodecoord, stddsv, 1._wp, func2)
739 y_kahan = real(func2*strength_vol*strength_vel, kind=wp) - kcomp(5)%sf(i, j, k)
740 t_kahan = updatedvar(5)%sf(i, j, k) + y_kahan
741 kcomp(5)%sf(i, j, k) = (t_kahan - updatedvar(5)%sf(i, j, k)) - y_kahan
742 updatedvar(5)%sf(i, j, k) = t_kahan
743 end if
744 end do
745 end do
746 end do
747 end do
748 end do
749 end do
750 end do
751
752# 281 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
753#if defined(MFC_OpenACC)
754# 281 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
755!$acc end parallel loop
756# 281 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
757#elif defined(MFC_OpenMP)
758# 281 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
759
760# 281 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
761!$omp end target teams loop
762# 281 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
763#endif
764
765 end subroutine s_gaussian
766
767 !> Evaluate the Gaussian kernel at a grid node for a given bubble center
768 subroutine s_applygaussian(center, cellaux, nodecoord, stddsv, strength_idx, func)
769
770
771# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
772#ifdef _CRAYFTN
773# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
774#if MFC_OpenACC
775# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
776!$acc routine seq
777# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
778#elif MFC_OpenMP
779# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
780
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!$omp declare target device_type(any)
785# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
786#else
787# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
788!DIR$ INLINEALWAYS s_applygaussian
789# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
790#endif
791# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
792#elif MFC_OpenACC
793# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
794!$acc routine seq
795# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
796#elif MFC_OpenMP
797# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
798
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!$omp declare target device_type(any)
803# 288 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
804#endif
805
806 real(wp), dimension(3), intent(in) :: center
807 integer, dimension(3), intent(in) :: cellaux
808 real(wp), dimension(3), intent(in) :: nodecoord
809 real(wp), intent(in) :: stddsv
810 real(wp), intent(in) :: strength_idx
811 real(wp), intent(out) :: func
812 integer :: i
813 real(wp) :: distance
814 real(wp) :: theta, dtheta, L2, dzp, Lz2, zc
815 real(wp) :: Nr, Nr_count
816
817 distance = sqrt((center(1) - nodecoord(1))**2._wp + (center(2) - nodecoord(2))**2._wp + (center(3) - nodecoord(3))**2._wp)
818
819 if (num_dims == 3) then
820 !> 3D gaussian function
821 func = exp(-0.5_wp*(distance/stddsv)**2._wp)/(sqrt(2._wp*pi)*stddsv)**3._wp
822 else
823 if (cyl_coord) then
824 !> 2D cylindrical function:
825 ! We smear particles in the azimuthal direction for given r
826 theta = 0._wp
827 nr = ceiling(2._wp*pi*nodecoord(2)/(y_cb(cellaux(2)) - y_cb(cellaux(2) - 1)))
828 dtheta = 2._wp*pi/nr
829 l2 = center(2)**2._wp + nodecoord(2)**2._wp - 2._wp*center(2)*nodecoord(2)*cos(theta)
830 distance = sqrt((center(1) - nodecoord(1))**2._wp + l2)
831 ! Factor 2._wp is for symmetry (upper half of the 2D field (+r) is considered)
832 func = dtheta/2._wp/pi*exp(-0.5_wp*(distance/stddsv)**2._wp)/(sqrt(2._wp*pi)*stddsv)**3._wp
833 nr_count = 0._wp
834 do while (nr_count < nr - 1._wp)
835 nr_count = nr_count + 1._wp
836 theta = nr_count*dtheta
837 ! trigonometric relation
838 l2 = center(2)**2._wp + nodecoord(2)**2._wp - 2._wp*center(2)*nodecoord(2)*cos(theta)
839 distance = sqrt((center(1) - nodecoord(1))**2._wp + l2)
840 ! nodecoord(2)*dtheta is the azimuthal width of the cell
841 func = func + dtheta/2._wp/pi*exp(-0.5_wp*(distance/stddsv)**2._wp)/(sqrt(2._wp*pi)*stddsv) &
842 & **(3._wp*(strength_idx + 1._wp))
843 end do
844 else
845 !> 2D cartesian function: Equation (48) from Maeda and Colonius 2018
846 ! We smear particles considering a virtual depth (lag_params%charwidth) with lag_params%charNz cells
847 dzp = (lag_params%charwidth/(lag_params%charNz + 1._wp))
848
849 func = 0._wp
850 do i = 0, lag_params%charNz
851 zc = (-lag_params%charwidth/2._wp + dzp*(0.5_wp + i)) ! Center of virtual cell i in z-direction
852 lz2 = (center(3) - zc)**2._wp
853 distance = sqrt((center(1) - nodecoord(1))**2._wp + (center(2) - nodecoord(2))**2._wp + lz2)
854 func = func + dzp/lag_params%charwidth*exp(-0.5_wp*(distance/stddsv)**2._wp)/(sqrt(2._wp*pi)*stddsv)**3._wp
855 end do
856 end if
857 end if
858
859 end subroutine s_applygaussian
860
861 !> Check if the current cell is outside the computational domain including ghost cells
862 subroutine s_check_celloutside(cellaux, celloutside)
863
864
865# 348 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
866#ifdef _CRAYFTN
867# 348 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
868#if MFC_OpenACC
869# 348 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
870!$acc routine seq
871# 348 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
872#elif MFC_OpenMP
873# 348 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
874
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!$omp declare target device_type(any)
879# 348 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
880#else
881# 348 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
882!DIR$ INLINEALWAYS s_check_celloutside
883# 348 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
884#endif
885# 348 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
886#elif MFC_OpenACC
887# 348 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
888!$acc routine seq
889# 348 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
890#elif MFC_OpenMP
891# 348 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
892
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!$omp declare target device_type(any)
897# 348 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
898#endif
899
900 integer, dimension(3), intent(inout) :: cellaux
901 logical, intent(out) :: celloutside
902
903 celloutside = .false.
904
905 if (num_dims == 2) then
906 if ((cellaux(1) < -buff_size) .or. (cellaux(2) < -buff_size)) then
907 celloutside = .true.
908 end if
909 if (cyl_coord .and. y_cc(cellaux(2)) < 0._wp) then
910 celloutside = .true.
911 end if
912 if ((cellaux(2) > n + buff_size) .or. (cellaux(1) > m + buff_size)) then
913 celloutside = .true.
914 end if
915 else
916 if ((cellaux(3) < -buff_size) .or. (cellaux(1) < -buff_size) .or. (cellaux(2) < -buff_size)) then
917 celloutside = .true.
918 end if
919
920 if ((cellaux(3) > p + buff_size) .or. (cellaux(2) > n + buff_size) .or. (cellaux(1) > m + buff_size)) then
921 celloutside = .true.
922 end if
923 end if
924
925 end subroutine s_check_celloutside
926
927 !> Relocate cells that intersect a symmetric boundary
928 subroutine s_shift_cell_symmetric_bc(cellaux, cell)
929
930
931# 380 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
932#ifdef _CRAYFTN
933# 380 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
934#if MFC_OpenACC
935# 380 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
936!$acc routine seq
937# 380 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
938#elif MFC_OpenMP
939# 380 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
940
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!$omp declare target device_type(any)
945# 380 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
946#else
947# 380 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
948!DIR$ INLINEALWAYS s_shift_cell_symmetric_bc
949# 380 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
950#endif
951# 380 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
952#elif MFC_OpenACC
953# 380 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
954!$acc routine seq
955# 380 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
956#elif MFC_OpenMP
957# 380 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
958
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!$omp declare target device_type(any)
963# 380 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
964#endif
965
966 integer, dimension(3), intent(inout) :: cellaux
967 integer, dimension(3), intent(in) :: cell
968
969 ! x-dir
970 if (bc_x%beg == bc_reflective .and. (cell(1) <= mapcells - 1)) then
971 cellaux(1) = abs(cellaux(1)) - 1
972 end if
973 if (bc_x%end == bc_reflective .and. (cell(1) >= m + 1 - mapcells)) then
974 cellaux(1) = cellaux(1) - (2*(cellaux(1) - m) - 1)
975 end if
976
977 ! y-dir
978 if (bc_y%beg == bc_reflective .and. (cell(2) <= mapcells - 1)) then
979 cellaux(2) = abs(cellaux(2)) - 1
980 end if
981 if (bc_y%end == bc_reflective .and. (cell(2) >= n + 1 - mapcells)) then
982 cellaux(2) = cellaux(2) - (2*(cellaux(2) - n) - 1)
983 end if
984
985 if (p > 0) then
986 ! z-dir
987 if (bc_z%beg == bc_reflective .and. (cell(3) <= mapcells - 1)) then
988 cellaux(3) = abs(cellaux(3)) - 1
989 end if
990 if (bc_z%end == bc_reflective .and. (cell(3) >= p + 1 - mapcells)) then
991 cellaux(3) = cellaux(3) - (2*(cellaux(3) - p) - 1)
992 end if
993 end if
994
995 end subroutine s_shift_cell_symmetric_bc
996
997 !> Calculates the standard deviation of the bubble being smeared in the Eulerian framework.
998 subroutine s_compute_stddsv(cell, volpart, stddsv)
999
1000
1001# 416 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1002#ifdef _CRAYFTN
1003# 416 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1004#if MFC_OpenACC
1005# 416 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1006!$acc routine seq
1007# 416 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1008#elif MFC_OpenMP
1009# 416 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1010
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!$omp declare target device_type(any)
1015# 416 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1016#else
1017# 416 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1018!DIR$ INLINEALWAYS s_compute_stddsv
1019# 416 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1020#endif
1021# 416 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1022#elif MFC_OpenACC
1023# 416 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1024!$acc routine seq
1025# 416 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1026#elif MFC_OpenMP
1027# 416 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1028
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!$omp declare target device_type(any)
1033# 416 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1034#endif
1035
1036 integer, dimension(3), intent(in) :: cell
1037 real(wp), intent(in) :: volpart
1038 real(wp), intent(out) :: stddsv
1039 real(wp) :: chardist, charvol
1040 real(wp) :: rad
1041
1042 !> Compute characteristic distance
1043 chardist = sqrt(dx(cell(1))*dy(cell(2)))
1044 if (p > 0) chardist = (dx(cell(1))*dy(cell(2))*dz(cell(3)))**(1._wp/3._wp)
1045
1046 !> Compute characteristic volume
1047 if (p > 0) then
1048 charvol = dx(cell(1))*dy(cell(2))*dz(cell(3))
1049 else
1050 if (cyl_coord) then
1051 charvol = dx(cell(1))*dy(cell(2))*y_cc(cell(2))*2._wp*pi
1052 else
1053 charvol = dx(cell(1))*dy(cell(2))*lag_params%charwidth
1054 end if
1055 end if
1056
1057 !> Compute Standard deviaton
1058 if ((volpart/charvol) > 0.5_wp*lag_params%valmaxvoid .or. (lag_params%smooth_type == 1)) then
1059 rad = (3._wp*volpart/(4._wp*pi))**(1._wp/3._wp)
1060 stddsv = 1._wp*lag_params%epsilonb*max(chardist, rad)
1061 else
1062 stddsv = 0._wp
1063 end if
1064
1065 end subroutine s_compute_stddsv
1066
1067 !> Compute the characteristic cell volume
1068 subroutine s_get_char_vol(cellx, celly, cellz, Charvol)
1069
1070
1071# 452 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1072#ifdef _CRAYFTN
1073# 452 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1074#if MFC_OpenACC
1075# 452 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1076!$acc routine seq
1077# 452 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1078#elif MFC_OpenMP
1079# 452 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1080
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!$omp declare target device_type(any)
1085# 452 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1086#else
1087# 452 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1088!DIR$ INLINEALWAYS s_get_char_vol
1089# 452 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1090#endif
1091# 452 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1092#elif MFC_OpenACC
1093# 452 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1094!$acc routine seq
1095# 452 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1096#elif MFC_OpenMP
1097# 452 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1098
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!$omp declare target device_type(any)
1103# 452 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1104#endif
1105
1106 integer, intent(in) :: cellx, celly, cellz
1107 real(wp), intent(out) :: Charvol
1108
1109 if (p > 0) then
1110 charvol = dx(cellx)*dy(celly)*dz(cellz)
1111 else
1112 if (cyl_coord) then
1113 charvol = dx(cellx)*dy(celly)*y_cc(celly)*2._wp*pi
1114 else
1115 charvol = dx(cellx)*dy(celly)*lag_params%charwidth
1116 end if
1117 end if
1118
1119 end subroutine s_get_char_vol
1120
1121 !> Convert bubble computational coordinates from real to integer cell indices
1122 subroutine s_get_cell(s_cell, get_cell)
1123
1124
1125# 472 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1126#ifdef _CRAYFTN
1127# 472 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1128#if MFC_OpenACC
1129# 472 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1130!$acc routine seq
1131# 472 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1132#elif MFC_OpenMP
1133# 472 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1134
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!$omp declare target device_type(any)
1139# 472 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1140#else
1141# 472 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1142!DIR$ INLINEALWAYS s_get_cell
1143# 472 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1144#endif
1145# 472 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1146#elif MFC_OpenACC
1147# 472 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1148!$acc routine seq
1149# 472 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1150#elif MFC_OpenMP
1151# 472 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1152
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!$omp declare target device_type(any)
1157# 472 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1158#endif
1159
1160 real(wp), dimension(3), intent(in) :: s_cell
1161 integer, dimension(3), intent(out) :: get_cell
1162 integer :: i
1163
1164 get_cell(:) = int(s_cell(:))
1165 do i = 1, num_dims
1166 if (s_cell(i) < 0._wp) get_cell(i) = get_cell(i) - 1
1167 end do
1168
1169 end subroutine s_get_cell
1170
1171 !> Precompute cell-centered pressure gradients (dp/dx, dp/dy, dp/dz)
1172 subroutine s_compute_pressure_gradients(q_prim_vf)
1173
1174 type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
1175 integer :: i, j, k, r
1176
1177 ! dp/dx at all cell centers
1178
1179
1180# 493 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1181
1182# 493 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1183#if defined(MFC_OpenACC)
1184# 493 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1185!$acc parallel loop collapse(3) gang vector default(present) private(i, j, k, r)
1186# 493 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1187#elif defined(MFC_OpenMP)
1188# 493 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1189
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!$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)
1196# 493 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1197#endif
1198 do k = 0, p
1199 do j = 0, n
1200 do i = 0, m
1201 grad_p_x(i, j, k) = 0._wp
1202
1203# 498 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1204#if defined(MFC_OpenACC)
1205# 498 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1206!$acc loop seq
1207# 498 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1208#elif defined(MFC_OpenMP)
1209# 498 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1210
1211# 498 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1212#endif
1213 do r = -fd_number, fd_number
1214 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)
1215 end do
1216 end do
1217 end do
1218 end do
1219
1220# 505 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1221#if defined(MFC_OpenACC)
1222# 505 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1223!$acc end parallel loop
1224# 505 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1225#elif defined(MFC_OpenMP)
1226# 505 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1227
1228# 505 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1229!$omp end target teams loop
1230# 505 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1231#endif
1232
1233 ! dp/dy at all cell centers
1234 if (n > 0) then
1235
1236# 509 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1237
1238# 509 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1239#if defined(MFC_OpenACC)
1240# 509 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1241!$acc parallel loop collapse(3) gang vector default(present) private(i, j, k, r)
1242# 509 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1243#elif defined(MFC_OpenMP)
1244# 509 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1245
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!$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)
1252# 509 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1253#endif
1254 do k = 0, p
1255 do j = 0, n
1256 do i = 0, m
1257 grad_p_y(i, j, k) = 0._wp
1258
1259# 514 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1260#if defined(MFC_OpenACC)
1261# 514 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1262!$acc loop seq
1263# 514 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1264#elif defined(MFC_OpenMP)
1265# 514 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1266
1267# 514 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1268#endif
1269 do r = -fd_number, fd_number
1270 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)
1271 end do
1272 end do
1273 end do
1274 end do
1275
1276# 521 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1277#if defined(MFC_OpenACC)
1278# 521 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1279!$acc end parallel loop
1280# 521 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1281#elif defined(MFC_OpenMP)
1282# 521 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1283
1284# 521 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1285!$omp end target teams loop
1286# 521 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1287#endif
1288 end if
1289
1290 ! dp/dz at all cell centers
1291 if (p > 0) then
1292
1293# 526 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1294
1295# 526 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1296#if defined(MFC_OpenACC)
1297# 526 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1298!$acc parallel loop collapse(3) gang vector default(present) private(i, j, k, r)
1299# 526 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1300#elif defined(MFC_OpenMP)
1301# 526 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1302
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!$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)
1309# 526 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1310#endif
1311 do k = 0, p
1312 do j = 0, n
1313 do i = 0, m
1314 grad_p_z(i, j, k) = 0._wp
1315
1316# 531 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1317#if defined(MFC_OpenACC)
1318# 531 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1319!$acc loop seq
1320# 531 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1321#elif defined(MFC_OpenMP)
1322# 531 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1323
1324# 531 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1325#endif
1326 do r = -fd_number, fd_number
1327 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)
1328 end do
1329 end do
1330 end do
1331 end do
1332
1333# 538 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1334#if defined(MFC_OpenACC)
1335# 538 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1336!$acc end parallel loop
1337# 538 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1338#elif defined(MFC_OpenMP)
1339# 538 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1340
1341# 538 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1342!$omp end target teams loop
1343# 538 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1344#endif
1345 end if
1346
1347 end subroutine s_compute_pressure_gradients
1348
1349 !! Interpolate the velocity of Eulerian field at the position of the bubble.
1350 function f_interpolate_velocity(pos, cell, i, q_prim_vf) result(v)
1351
1352
1353# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1354#ifdef _CRAYFTN
1355# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1356#if MFC_OpenACC
1357# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1358!$acc routine seq
1359# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1360#elif MFC_OpenMP
1361# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1362
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!$omp declare target device_type(any)
1367# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1368#else
1369# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1370!DIR$ NOINLINE f_interpolate_velocity
1371# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1372#endif
1373# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1374#elif MFC_OpenACC
1375# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1376!$acc routine seq
1377# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1378#elif MFC_OpenMP
1379# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1380
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!$omp declare target device_type(any)
1385# 546 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1386#endif
1387
1388 real(wp), intent(in) :: pos
1389 integer, dimension(3), intent(in) :: cell
1390 integer, intent(in) :: i
1391 type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
1392 real(wp) :: v
1393 real(wp), dimension(5) :: xi, eta, l
1394
1395 if (fd_order == 2) then
1396 if (i == 1) then
1397 xi(1) = x_cc(cell(1) - 1)
1398 eta(1) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1) - 1, cell(2), cell(3))
1399 xi(2) = x_cc(cell(1))
1400 eta(2) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1), cell(2), cell(3))
1401 xi(3) = x_cc(cell(1) + 1)
1402 eta(3) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1) + 1, cell(2), cell(3))
1403 else if (i == 2) then
1404 xi(1) = y_cc(cell(2) - 1)
1405 eta(1) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2) - 1, cell(3))
1406 xi(2) = y_cc(cell(2))
1407 eta(2) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2), cell(3))
1408 xi(3) = y_cc(cell(2) + 1)
1409 eta(3) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2) + 1, cell(3))
1410 else if (i == 3) then
1411 xi(1) = z_cc(cell(3) - 1)
1412 eta(1) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3) - 1)
1413 xi(2) = z_cc(cell(3))
1414 eta(2) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3))
1415 xi(3) = z_cc(cell(3) + 1)
1416 eta(3) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3) + 1)
1417 end if
1418
1419 l(1) = ((pos - xi(2))*(pos - xi(3)))/((xi(1) - xi(2))*(xi(1) - xi(3)))
1420 l(2) = ((pos - xi(1))*(pos - xi(3)))/((xi(2) - xi(1))*(xi(2) - xi(3)))
1421 l(3) = ((pos - xi(1))*(pos - xi(2)))/((xi(3) - xi(1))*(xi(3) - xi(2)))
1422
1423 v = l(1)*eta(1) + l(2)*eta(2) + l(3)*eta(3)
1424 else if (fd_order == 4) then
1425 if (i == 1) then
1426 xi(1) = x_cc(cell(1) - 2)
1427 eta(1) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1) - 2, cell(2), cell(3))
1428 xi(2) = x_cc(cell(1) - 1)
1429 eta(2) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1) - 1, cell(2), cell(3))
1430 xi(3) = x_cc(cell(1))
1431 eta(3) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1), cell(2), cell(3))
1432 xi(4) = x_cc(cell(1) + 1)
1433 eta(4) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1) + 1, cell(2), cell(3))
1434 xi(5) = x_cc(cell(1) + 2)
1435 eta(5) = q_prim_vf(eqn_idx%mom%beg)%sf(cell(1) + 2, cell(2), cell(3))
1436 else if (i == 2) then
1437 xi(1) = y_cc(cell(2) - 2)
1438 eta(1) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2) - 2, cell(3))
1439 xi(2) = y_cc(cell(2) - 1)
1440 eta(2) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2) - 1, cell(3))
1441 xi(3) = y_cc(cell(2))
1442 eta(3) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2), cell(3))
1443 xi(4) = y_cc(cell(2) + 1)
1444 eta(4) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2) + 1, cell(3))
1445 xi(5) = y_cc(cell(2) + 2)
1446 eta(5) = q_prim_vf(eqn_idx%mom%beg + 1)%sf(cell(1), cell(2) + 2, cell(3))
1447 else if (i == 3) then
1448 xi(1) = z_cc(cell(3) - 2)
1449 eta(1) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3) - 2)
1450 xi(2) = z_cc(cell(3) - 1)
1451 eta(2) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3) - 1)
1452 xi(3) = z_cc(cell(3))
1453 eta(3) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3))
1454 xi(4) = z_cc(cell(3) + 1)
1455 eta(4) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3) + 1)
1456 xi(5) = z_cc(cell(3) + 2)
1457 eta(5) = q_prim_vf(eqn_idx%mom%end)%sf(cell(1), cell(2), cell(3) + 2)
1458 end if
1459
1460 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)) &
1461 & *(xi(1) - xi(5)))
1462 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)) &
1463 & *(xi(2) - xi(5)))
1464 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)) &
1465 & *(xi(3) - xi(5)))
1466 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)) &
1467 & *(xi(4) - xi(5)))
1468 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)) &
1469 & *(xi(5) - xi(4)))
1470
1471 v = l(1)*eta(1) + l(2)*eta(2) + l(3)*eta(3) + l(4)*eta(4) + l(5)*eta(5)
1472 end if
1473
1474 end function f_interpolate_velocity
1475
1476 !! Calculate the force on a bubble based on the pressure gradient, velocity, and drag model.
1477 !! @param mg Mass of the gas in the bubble
1478 !! @param mv Mass of the liquid in the bubble
1479 !! @param Re Reynolds number
1480 !! @param i Direction of the velocity (1: x, 2: y, 3: z)
1481 function f_get_bubble_force(pos, rad, rdot, vel, mg, mv, Re, rho, cell, i, q_prim_vf) result(force)
1482
1483
1484# 643 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1485#if MFC_OpenACC
1486# 643 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1487!$acc routine seq
1488# 643 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1489#elif MFC_OpenMP
1490# 643 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1491
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!$omp declare target device_type(any)
1496# 643 "/home/runner/work/MFC/MFC/src/simulation/m_bubbles_EL_kernels.fpp"
1497#endif
1498 real(wp), intent(in) :: pos, rad, rdot, mg, mv, re, rho, vel
1499 integer, dimension(3), intent(in) :: cell
1500 integer, intent(in) :: i
1501 type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
1502 real(wp) :: dp, vol, force
1503 real(wp) :: v_rel
1504
1505 if (fd_order > 1) then
1506 v_rel = vel - f_interpolate_velocity(pos, cell, i, q_prim_vf)
1507 else
1508 v_rel = vel - q_prim_vf(eqn_idx%mom%beg + i - 1)%sf(cell(1), cell(2), cell(3))
1509 end if
1510
1511 force = 0._wp
1512
1513 if (lag_params%drag_model == 1) then ! Free slip Stokes drag
1514 force = force - (4._wp*pi*rad*v_rel)/re
1515 else if (lag_params%drag_model == 2) then ! No slip Stokes drag
1516 force = force - (6._wp*pi*rad*v_rel)/re
1517 else if (lag_params%drag_model == 3) then ! Levich drag
1518 force = force - (12._wp*pi*rad*v_rel)/re
1519 end if
1520
1521 if (lag_pressure_force) then
1522 ! Use precomputed cell-centered pressure gradients
1523 if (i == 1) then
1524 dp = grad_p_x(cell(1), cell(2), cell(3))
1525 else if (i == 2) then
1526 dp = grad_p_y(cell(1), cell(2), cell(3))
1527 else if (i == 3) then
1528 dp = grad_p_z(cell(1), cell(2), cell(3))
1529 end if
1530
1531 vol = (4._wp/3._wp)*pi*(rad**3._wp)
1532 force = force - vol*dp
1533 end if
1534
1535 if (lag_params%gravity_force) then
1536 force = force + (mg + mv)*accel_bf(i)
1537 end if
1538
1539 end function f_get_bubble_force
1540
1541end 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.