MFC
Exascale flow solver
Loading...
Searching...
No Matches
m_helper.fpp.f90
Go to the documentation of this file.
1# 1 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
2# 1 "/home/runner/work/MFC/MFC/src/common/include/case.fpp" 1
3! This file exists so that Fypp can be run without generating case.fpp files for
4! each target. This is useful when generating documentation, for example. This
5! should also let MFC be built with CMake directly, without invoking mfc.sh.
6
7! For pre-process.
8# 8 "/home/runner/work/MFC/MFC/src/common/include/case.fpp"
9
10! For moving immersed boundaries in simulation
11# 12 "/home/runner/work/MFC/MFC/src/common/include/case.fpp"
12# 2 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp" 2
13# 1 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 1
14# 1 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 1
15# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
16# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
17# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
18# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
19# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
20# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
21
22# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
23# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
24# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
25
26# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
27
28# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
29
30# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
31
32# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
33
34# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
35
36# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
37
38# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
39
40# 145 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
41! New line at end of file is required for FYPP
42# 2 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
43# 1 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 1
44# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
45# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
46# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
47# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
48# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
49# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
50
51# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
52# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
53# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
54
55# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
56
57# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
58
59# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
60
61# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
62
63# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
64
65# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
66
67# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
68
69# 145 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
70! New line at end of file is required for FYPP
71# 2 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 2
72
73# 4 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
74# 5 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
75# 6 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
76# 7 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
77# 8 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
78
79# 20 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
80
81# 43 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
82
83# 48 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
84
85# 53 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
86
87# 58 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
88
89# 63 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
90
91# 68 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
92
93# 76 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
94
95# 81 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
96
97# 86 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
98
99# 91 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
100
101# 96 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
102
103# 101 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
104
105# 106 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
106
107# 111 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
108
109# 116 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
110
111# 121 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
112
113# 151 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
114
115# 192 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
116
117# 206 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
118
119# 231 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
120
121# 242 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
122
123# 244 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
124# 255 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
125
126# 284 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
127
128# 294 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
129
130# 304 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
131
132# 313 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
133
134# 330 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
135
136# 340 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
137
138# 347 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
139
140# 353 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
141
142# 359 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
143
144# 365 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
145
146# 371 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
147
148# 377 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
149! New line at end of file is required for FYPP
150# 3 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
151# 1 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 1
152# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
153# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
154# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
155# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
156# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
157# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
158
159# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
160# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
161# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
162
163# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
164
165# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
166
167# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
168
169# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
170
171# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
172
173# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
174
175# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
176
177# 145 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
178! New line at end of file is required for FYPP
179# 2 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 2
180
181# 7 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
182
183# 17 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
184
185# 22 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
186
187# 27 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
188
189# 32 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
190
191# 37 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
192
193# 42 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
194
195# 47 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
196
197# 52 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
198
199# 57 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
200
201# 62 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
202
203# 73 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
204
205# 78 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
206
207# 83 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
208
209# 88 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
210
211# 103 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
212
213# 131 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
214
215# 160 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
216
217# 175 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
218
219# 193 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
220
221# 215 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
222
223# 244 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
224
225# 259 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
226
227# 269 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
228
229# 278 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
230
231# 294 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
232
233# 304 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
234
235# 311 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
236! New line at end of file is required for FYPP
237# 4 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
238
239! GPU parallel region (scalar reductions, maxval/minval)
240# 23 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
241
242! GPU parallel loop over threads (most common GPU macro)
243# 43 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
244
245! Required closing for GPU_PARALLEL_LOOP
246# 55 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
247
248! Mark routine for device compilation
249# 112 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
250
251! Declare device-resident data
252# 130 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
253
254! Inner loop within a GPU parallel region
255# 145 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
256
257! Scoped GPU data region
258# 164 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
259
260! Host code with device pointers (for MPI with GPU buffers)
261# 193 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
262
263! Allocate device memory (unscoped)
264# 207 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
265
266! Free device memory
267# 219 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
268
269! Atomic operation on device
270# 231 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
271
272! End atomic capture block
273# 242 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
274
275! Copy data between host and device
276# 254 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
277
278! Synchronization barrier
279# 266 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
280
281! Import GPU library module (openacc or omp_lib)
282# 275 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
283
284! Emit code only for AMD compiler
285# 282 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
286
287! Emit code for non-Cray compilers
288# 289 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
289
290! Emit code only for Cray compiler
291# 296 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
292
293! Emit code for non-NVIDIA compilers
294# 303 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
295
296# 305 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
297# 306 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
298! New line at end of file is required for FYPP
299# 2 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 2
300
301# 14 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
302
303! Caution: This macro requires the use of a binding script to set CUDA_VISIBLE_DEVICES, such that we have one GPU device per MPI
304! rank. That's because for both cudaMemAdvise (preferred location) and cudaMemPrefetchAsync we use location = device_id = 0. For an
305! example see misc/nvidia_uvm/bind.sh.
306# 57 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
307
308! Allocate and create GPU device memory
309# 77 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
310
311! Free GPU device memory and deallocate
312# 85 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
313
314! Cray-specific GPU pointer setup for vector fields
315# 109 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
316
317! Cray-specific GPU pointer setup for scalar fields
318# 125 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
319
320! Cray-specific GPU pointer setup for acoustic source spatials
321# 150 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
322
323# 156 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
324
325# 163 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
326! New line at end of file is required for FYPP
327# 3 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp" 2
328
329!>
330!! @file
331!! @brief Contains module m_helper
332
333!> @brief Utility routines for bubble model setup, coordinate transforms, array sampling, and special functions
334module m_helper
335
336 use m_derived_types
337 use m_global_parameters
338 use ieee_arithmetic !< for checking nan
339
340 implicit none
341
342 private
343 public :: s_comp_n_from_prim, s_comp_n_from_cons, s_initialize_bubbles_model, s_initialize_nonpoly, s_simpson, s_transcoeff, &
344 & s_int_to_str, s_transform_vec, s_transform_triangle, s_transform_model, s_swap, f_cross, f_create_transform_matrix, &
345 & f_create_bbox, s_print_2d_array, f_xor, f_logical_to_int, associated_legendre, real_ylm, double_factorial, factorial, &
346 & f_cut_on, f_cut_off, s_downsample_data, s_upsample_data, s_cross_product, f_unit_vector, s_prng, modmul
347
348contains
349
350 !> Computes the bubble number density n from the primitive variables
351 subroutine s_comp_n_from_prim(vftmp, Rtmp, ntmp, weights)
352
353
354# 28 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
355#if MFC_OpenACC
356# 28 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
357!$acc routine seq
358# 28 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
359#elif MFC_OpenMP
360# 28 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
361
362# 28 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
363
364# 28 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
365!$omp declare target device_type(any)
366# 28 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
367#endif
368 real(wp), intent(in) :: vftmp
369 real(wp), dimension(nb), intent(in) :: rtmp
370 real(wp), intent(out) :: ntmp
371 real(wp), dimension(nb), intent(in) :: weights
372 real(wp) :: r3
373
374 r3 = dot_product(weights, rtmp**3._wp)
375 ntmp = (3._wp/(4._wp*pi))*vftmp/r3
376
377 end subroutine s_comp_n_from_prim
378
379 !> Compute the bubble number density from the conservative void fraction and weighted bubble radii.
380 subroutine s_comp_n_from_cons(vftmp, nRtmp, ntmp, weights)
381
382
383# 43 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
384#if MFC_OpenACC
385# 43 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
386!$acc routine seq
387# 43 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
388#elif MFC_OpenMP
389# 43 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
390
391# 43 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
392
393# 43 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
394!$omp declare target device_type(any)
395# 43 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
396#endif
397 real(wp), intent(in) :: vftmp
398 real(wp), dimension(nb), intent(in) :: nrtmp
399 real(wp), intent(out) :: ntmp
400 real(wp), dimension(nb), intent(in) :: weights
401 real(wp) :: nr3
402
403 nr3 = dot_product(weights, nrtmp**3._wp)
404 ntmp = sqrt((4._wp*pi/3._wp)*nr3/vftmp)
405
406 end subroutine s_comp_n_from_cons
407
408 !> Print a 2D real array to standard output, optionally dividing each element by a given scalar.
409 impure subroutine s_print_2d_array(A, div)
410
411 real(wp), dimension(:,:), intent(in) :: a
412 real(wp), optional, intent(in) :: div
413 integer :: i, j
414 integer :: local_m, local_n
415 real(wp) :: c
416
417 local_m = size(a, 1)
418 local_n = size(a, 2)
419
420 if (present(div)) then
421 c = div
422 else
423 c = 1._wp
424 end if
425
426 print *, local_m, local_n
427
428 do i = 1, local_m
429 do j = 1, local_n
430 write (*, fmt="(F12.4)", advance="no") a(i, j)/c
431 end do
432 write (*, fmt="(A1)") " "
433 end do
434 write (*, fmt="(A1)") " "
435
436 end subroutine s_print_2d_array
437
438 !> Initialize bubble model arrays for Euler or Lagrangian bubbles with polytropic or non-polytropic gas.
439 impure subroutine s_initialize_bubbles_model()
440
441 ! Allocate memory
442 if (bubbles_euler) then
443#ifdef MFC_DEBUG
444# 90 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
445 block
446# 90 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
447 use iso_fortran_env, only: output_unit
448# 90 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
449
450# 90 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
451 print *, 'm_helper.fpp:90: ', '@:ALLOCATE(weight(nb), R0(nb))'
452# 90 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
453
454# 90 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
455 call flush (output_unit)
456# 90 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
457 end block
458# 90 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
459#endif
460# 90 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
461 allocate (weight(nb), r0(nb))
462# 90 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
463
464# 90 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
465
466# 90 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
467
468# 90 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
469#if defined(MFC_OpenACC)
470# 90 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
471!$acc enter data create(weight, R0)
472# 90 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
473#elif defined(MFC_OpenMP)
474# 90 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
475!$omp target enter data map(always,alloc:weight, R0)
476# 90 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
477#endif
478 if (.not. polytropic) then
479#ifdef MFC_DEBUG
480# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
481 block
482# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
483 use iso_fortran_env, only: output_unit
484# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
485
486# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
487 print *, 'm_helper.fpp:92: ', '@:ALLOCATE(pb0(nb), Pe_T(nb), k_g(nb), k_v(nb), mass_g0(nb), mass_v0(nb))'
488# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
489
490# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
491 call flush (output_unit)
492# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
493 end block
494# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
495#endif
496# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
497 allocate (pb0(nb), pe_t(nb), k_g(nb), k_v(nb), mass_g0(nb), mass_v0(nb))
498# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
499
500# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
501
502# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
503
504# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
505
506# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
507
508# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
509
510# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
511
512# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
513#if defined(MFC_OpenACC)
514# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
515!$acc enter data create(pb0, Pe_T, k_g, k_v, mass_g0, mass_v0)
516# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
517#elif defined(MFC_OpenMP)
518# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
519!$omp target enter data map(always,alloc:pb0, Pe_T, k_g, k_v, mass_g0, mass_v0)
520# 92 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
521#endif
522#ifdef MFC_DEBUG
523# 93 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
524 block
525# 93 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
526 use iso_fortran_env, only: output_unit
527# 93 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
528
529# 93 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
530 print *, 'm_helper.fpp:93: ', '@:ALLOCATE(Re_trans_T(nb), Re_trans_c(nb), Im_trans_T(nb), Im_trans_c(nb))'
531# 93 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
532
533# 93 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
534 call flush (output_unit)
535# 93 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
536 end block
537# 93 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
538#endif
539# 93 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
540 allocate (re_trans_t(nb), re_trans_c(nb), im_trans_t(nb), im_trans_c(nb))
541# 93 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
542
543# 93 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
544
545# 93 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
546
547# 93 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
548
549# 93 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
550
551# 93 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
552#if defined(MFC_OpenACC)
553# 93 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
554!$acc enter data create(Re_trans_T, Re_trans_c, Im_trans_T, Im_trans_c)
555# 93 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
556#elif defined(MFC_OpenMP)
557# 93 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
558!$omp target enter data map(always,alloc:Re_trans_T, Re_trans_c, Im_trans_T, Im_trans_c)
559# 93 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
560#endif
561 else if (qbmm) then
562#ifdef MFC_DEBUG
563# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
564 block
565# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
566 use iso_fortran_env, only: output_unit
567# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
568
569# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
570 print *, 'm_helper.fpp:95: ', '@:ALLOCATE(pb0(nb))'
571# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
572
573# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
574 call flush (output_unit)
575# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
576 end block
577# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
578#endif
579# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
580 allocate (pb0(nb))
581# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
582
583# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
584
585# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
586#if defined(MFC_OpenACC)
587# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
588!$acc enter data create(pb0)
589# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
590#elif defined(MFC_OpenMP)
591# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
592!$omp target enter data map(always,alloc:pb0)
593# 95 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
594#endif
595 end if
596
597 ! Compute quadrature weights and nodes for polydisperse simulations
598 if (nb > 1) then
599 call s_simpson(weight, r0)
600 else if (nb == 1) then
601 r0 = 1._wp
602 weight = 1._wp
603 else
604 stop 'Invalid value of nb'
605 end if
606 r0 = r0*bub_pp%R0ref
607 end if
608
609 ! Initialize bubble variables
610 call s_initialize_bubble_vars()
611
612 end subroutine s_initialize_bubbles_model
613
614 !> Set bubble physical parameters and nondimensional numbers from the input configuration.
615 impure subroutine s_initialize_bubble_vars()
616
617 r0ref = bub_pp%R0ref; p0ref = bub_pp%p0ref
618 rho0ref = bub_pp%rho0ref
619 ss = bub_pp%ss; pv = bub_pp%pv; vd = bub_pp%vd
620 mu_l = bub_pp%mu_l; mu_v = bub_pp%mu_v; mu_g = bub_pp%mu_g
621 gam_v = bub_pp%gam_v; gam_g = bub_pp%gam_g
622 if (.not. polytropic) then
623 if (bubbles_euler) then
624 m_v = bub_pp%M_v; m_g = bub_pp%M_g
625 k_v = bub_pp%k_v; k_g = bub_pp%k_g
626 end if
627 r_v = bub_pp%R_v; r_g = bub_pp%R_g
628 tw = bub_pp%T0ref
629 end if
630 if (bubbles_lagrange) then
631 cp_v = bub_pp%cp_v; cp_g = bub_pp%cp_g
632 k_vl = bub_pp%k_v; k_gl = bub_pp%k_g
633 end if
634
635 ! Input quantities
636 if (bubbles_euler .and. (.not. polytropic)) then
637 if (thermal == 2) then
638 gam_m = 1._wp
639 else
640 gam_m = gam_g
641 end if
642 end if
643
644 ! Nondimensional numbers
645 eu = p0ref
646 ca = eu - pv
647 if (.not. f_is_default(bub_pp%ss)) web = 1._wp/ss
648 if (.not. f_is_default(bub_pp%mu_l)) re_inv = mu_l
649 if (.not. polytropic) pe_c = 1._wp/vd
650
651 if (bubbles_euler) then
652 ! Initialize variables for non-polytropic (Preston) model
653 if (.not. polytropic) then
654 call s_initialize_nonpoly()
655 end if
656 ! Initialize pb based on surface tension for qbmm (polytropic)
657 if (qbmm .and. polytropic) then
658 pb0 = eu
659 if (.not. f_is_default(web)) then
660 pb0 = pb0 + 2._wp/web/r0
661 end if
662 end if
663 end if
664
665 end subroutine s_initialize_bubble_vars
666
667 !> Initializes non-polydisperse bubble modeling
668 impure subroutine s_initialize_nonpoly()
669
670 integer :: ir
671 real(wp), dimension(nb) :: chi_vw0, cp_m0, k_m0, rho_m0, x_vw, omegan
672 real(wp), parameter :: k_poly = 1._wp !< polytropic index used to compute isothermal natural frequency
673 ! Chapman-Enskog transport coefficients for vapor-gas mixture, Ando JAS (2010)
674
675 phi_vg = (1._wp + sqrt(mu_v/mu_g)*(m_g/m_v)**(0.25_wp))**2/(sqrt(8._wp)*sqrt(1._wp + m_v/m_g))
676 phi_gv = (1._wp + sqrt(mu_g/mu_v)*(m_v/m_g)**(0.25_wp))**2/(sqrt(8._wp)*sqrt(1._wp + m_g/m_v))
677
678 ! Initial internal bubble pressure (Euler number + Laplace pressure)
679 pb0 = eu + 2._wp/web/r0
680
681 ! Vapor mass fraction at bubble wall, Ando JAS (2010)
682 chi_vw0 = 1._wp/(1._wp + r_v/r_g*(pb0/pv - 1._wp))
683
684 ! Mixture specific heat from mass-weighted vapor/gas contributions
685 cp_m0 = chi_vw0*r_v*gam_v/(gam_v - 1._wp) + (1._wp - chi_vw0)*r_g*gam_g/(gam_g - 1._wp)
686
687 ! mole fraction of vapor (Eq. 2.23 in Ando 2010)
688 x_vw = m_g*chi_vw0/(m_v + (m_g - m_v)*chi_vw0)
689
690 ! thermal conductivity for gas/vapor mixture (Eq. 2.21 in Ando 2010)
691 k_m0 = x_vw*k_v/(x_vw + (1._wp - x_vw)*phi_vg) + (1._wp - x_vw)*k_g/(x_vw*phi_gv + 1._wp - x_vw)
692 k_g(:) = k_g(:)/k_m0(:)
693 k_v(:) = k_v(:)/k_m0(:)
694
695 ! mixture density (Eq. 2.20 in Ando 2010)
696 rho_m0 = pv/(chi_vw0*r_v*tw)
697
698 ! mass of gas/vapor
699 mass_g0(:) = (4._wp*pi/3._wp)*(pb0(:) - pv)/(r_g*tw)*r0(:)**3
700 mass_v0(:) = (4._wp*pi/3._wp)*pv/(r_v*tw)*r0(:)**3
701
702 ! Peclet numbers
703 pe_t(:) = rho_m0*cp_m0(:)/k_m0(:)
704
705 ! Bubble natural frequency, Ando JAS (2010)
706 omegan(:) = sqrt(3._wp*k_poly*ca + 2._wp*(3._wp*k_poly - 1._wp)/(web*r0))/r0/sqrt(rho0ref)
707 do ir = 1, nb
708 call s_transcoeff(omegan(ir)*r0(ir), pe_t(ir)*r0(ir), re_trans_t(ir), im_trans_t(ir))
709 call s_transcoeff(omegan(ir)*r0(ir), pe_c*r0(ir), re_trans_c(ir), im_trans_c(ir))
710 end do
711 im_trans_t = 0._wp
712
713 end subroutine s_initialize_nonpoly
714
715 !> Computes the transfer coefficient for the non-polytropic bubble compression process
716 elemental subroutine s_transcoeff(omega, peclet, Re_trans, Im_trans)
717
718 real(wp), intent(in) :: omega, peclet
719 real(wp), intent(out) :: re_trans, im_trans
720 complex(wp) :: imag, trans, c1, c2, c3
721
722 imag = (0._wp, 1._wp)
723
724 c1 = imag*omega*peclet
725 c2 = sqrt(c1)
726 c3 = (exp(c2) - exp(-c2))/(exp(c2) + exp(-c2)) ! TANH(c2)
727 trans = ((c2/c3 - 1._wp)**(-1) - 3._wp/c1)**(-1) ! transfer function
728
729 re_trans = trans
730 im_trans = aimag(trans)
731
732 end subroutine s_transcoeff
733
734 !> Convert an integer to its trimmed string representation.
735 elemental subroutine s_int_to_str(i, res)
736
737 integer, intent(in) :: i
738 character(len=*), intent(inout) :: res
739
740 write (res, '(I0)') i
741 res = trim(res)
742
743 end subroutine s_int_to_str
744
745 !> Compute the Simpson weights for quadrature
746 subroutine s_simpson(local_weight, local_R0)
747
748 real(wp), dimension(:), intent(inout) :: local_weight
749 real(wp), dimension(:), intent(inout) :: local_r0
750 integer :: ir
751 real(wp) :: r0mn, r0mx, dphi, tmp, sd
752 real(wp), dimension(nb) :: phi
753
754 sd = poly_sigma
755 r0mn = 0.8_wp*exp(-2.8_wp*sd)
756 r0mx = 0.2_wp*exp(9.5_wp*sd) + 1._wp
757
758 ! phi = ln( R0 ) & return R0
759 do ir = 1, nb
760 phi(ir) = log(r0mn) + (ir - 1._wp)*log(r0mx/r0mn)/(nb - 1._wp)
761 local_r0(ir) = exp(phi(ir))
762 end do
763
764# 266 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
765 dphi = phi(2) - phi(1)
766# 268 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
767
768 ! weights for quadrature using Simpson's rule
769 do ir = 2, nb - 1
770 ! Gaussian
771 tmp = exp(-0.5_wp*(phi(ir)/sd)**2)/sqrt(2._wp*pi)/sd
772 if (mod(ir, 2) == 0) then
773 local_weight(ir) = tmp*4._wp*dphi/3._wp
774 else
775 local_weight(ir) = tmp*2._wp*dphi/3._wp
776 end if
777 end do
778 tmp = exp(-0.5_wp*(phi(1)/sd)**2)/sqrt(2._wp*pi)/sd
779 local_weight(1) = tmp*dphi/3._wp
780 tmp = exp(-0.5_wp*(phi(nb)/sd)**2)/sqrt(2._wp*pi)/sd
781 local_weight(nb) = tmp*dphi/3._wp
782
783 end subroutine s_simpson
784
785 !> Compute the cross product of two vectors.
786 pure function f_cross(a, b) result(c)
787
788
789# 289 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
790#if MFC_OpenACC
791# 289 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
792!$acc routine seq
793# 289 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
794#elif MFC_OpenMP
795# 289 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
796
797# 289 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
798
799# 289 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
800!$omp declare target device_type(any)
801# 289 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
802#endif
803
804 real(wp), dimension(3), intent(in) :: a, b
805 real(wp), dimension(3) :: c
806
807 c(1) = a(2)*b(3) - a(3)*b(2)
808 c(2) = a(3)*b(1) - a(1)*b(3)
809 c(3) = a(1)*b(2) - a(2)*b(1)
810
811 end function f_cross
812
813 !> Generate a unit vector uniformly distributed on the sphere from two random parameters.
814 function f_unit_vector(theta, eta) result(vec)
815
816
817# 303 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
818#if MFC_OpenACC
819# 303 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
820!$acc routine seq
821# 303 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
822#elif MFC_OpenMP
823# 303 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
824
825# 303 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
826
827# 303 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
828!$omp declare target device_type(any)
829# 303 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
830#endif
831 real(wp), intent(in) :: theta, eta
832 real(wp) :: zeta, xi
833 real(wp), dimension(3) :: vec
834
835 xi = 2._wp*pi*theta
836 zeta = acos(2._wp*eta - 1._wp)
837 vec(1) = sin(zeta)*cos(xi)
838 vec(2) = sin(zeta)*sin(xi)
839 vec(3) = cos(zeta)
840
841 end function f_unit_vector
842
843 !> Generate a pseudo-random number between 0 and 1 using a linear congruential generator.
844 subroutine s_prng(var, seed)
845
846
847# 319 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
848#if MFC_OpenACC
849# 319 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
850!$acc routine seq
851# 319 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
852#elif MFC_OpenMP
853# 319 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
854
855# 319 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
856
857# 319 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
858!$omp declare target device_type(any)
859# 319 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
860#endif
861 integer, intent(inout) :: seed
862 real(wp), intent(out) :: var
863
864 seed = mod(modmul(seed), modulus)
865 var = seed/real(modulus, wp)
866
867 end subroutine s_prng
868
869 !> Compute a modular multiplication step for the linear congruential pseudo-random number generator.
870 function modmul(a) result(val)
871
872
873# 331 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
874#if MFC_OpenACC
875# 331 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
876!$acc routine seq
877# 331 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
878#elif MFC_OpenMP
879# 331 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
880
881# 331 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
882
883# 331 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
884!$omp declare target device_type(any)
885# 331 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
886#endif
887 integer, intent(in) :: a
888 integer :: val
889 real(wp) :: x, y
890
891 x = (multiplier/real(modulus, wp))*a + (increment/real(modulus, wp))
892 y = nint((x - floor(x))*decimal_trim)/decimal_trim
893 val = nint(y*modulus)
894
895 end function modmul
896
897 !> Compute the cross product c = a x b of two 3D vectors.
898 subroutine s_cross_product(a, b, c)
899
900
901# 345 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
902#if MFC_OpenACC
903# 345 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
904!$acc routine seq
905# 345 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
906#elif MFC_OpenMP
907# 345 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
908
909# 345 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
910
911# 345 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
912!$omp declare target device_type(any)
913# 345 "/home/runner/work/MFC/MFC/src/common/m_helper.fpp"
914#endif
915 real(wp), intent(in) :: a(3), b(3)
916 real(wp), intent(out) :: c(3)
917
918 c(1) = a(2)*b(3) - a(3)*b(2)
919 c(2) = a(3)*b(1) - a(1)*b(3)
920 c(3) = a(1)*b(2) - a(2)*b(1)
921
922 end subroutine s_cross_product
923
924 !> Swap two real numbers.
925 elemental subroutine s_swap(lhs, rhs)
926
927 real(wp), intent(inout) :: lhs, rhs
928 real(wp) :: ltemp
929
930 ltemp = lhs
931 lhs = rhs
932 rhs = ltemp
933
934 end subroutine s_swap
935
936 !> Create a transformation matrix.
937 function f_create_transform_matrix(param, center) result(out_matrix)
938
939 type(ic_model_parameters), intent(in) :: param
940 real(wp), dimension(1:3), optional, intent(in) :: center
941 real(wp), dimension(1:4,1:4) :: sc, rz, rx, ry, tr, t_back, t_to_origin, out_matrix
942
943 sc = transpose(reshape([param%scale(1), 0._wp, 0._wp, 0._wp, 0._wp, param%scale(2), 0._wp, 0._wp, 0._wp, 0._wp, &
944 & param%scale(3), 0._wp, 0._wp, 0._wp, 0._wp, 1._wp], shape(sc)))
945
946 rz = transpose(reshape([cos(param%rotate(3)), -sin(param%rotate(3)), 0._wp, 0._wp, sin(param%rotate(3)), &
947 & cos(param%rotate(3)), 0._wp, 0._wp, 0._wp, 0._wp, 1._wp, 0._wp, 0._wp, 0._wp, 0._wp, 1._wp], shape(rz)))
948
949 rx = transpose(reshape([1._wp, 0._wp, 0._wp, 0._wp, 0._wp, cos(param%rotate(1)), -sin(param%rotate(1)), 0._wp, 0._wp, &
950 & sin(param%rotate(1)), cos(param%rotate(1)), 0._wp, 0._wp, 0._wp, 0._wp, 1._wp], shape(rx)))
951
952 ry = transpose(reshape([cos(param%rotate(2)), 0._wp, sin(param%rotate(2)), 0._wp, 0._wp, 1._wp, 0._wp, 0._wp, &
953 & -sin(param%rotate(2)), 0._wp, cos(param%rotate(2)), 0._wp, 0._wp, 0._wp, 0._wp, 1._wp], shape(ry)))
954
955 tr = transpose(reshape([1._wp, 0._wp, 0._wp, param%translate(1), 0._wp, 1._wp, 0._wp, param%translate(2), 0._wp, 0._wp, &
956 & 1._wp, param%translate(3), 0._wp, 0._wp, 0._wp, 1._wp], shape(tr)))
957
958 if (present(center)) then
959 ! Translation matrix to move center to the origin
960 t_to_origin = transpose(reshape([1._wp, 0._wp, 0._wp, -center(1), 0._wp, 1._wp, 0._wp, -center(2), 0._wp, 0._wp, &
961 & 1._wp, -center(3), 0._wp, 0._wp, 0._wp, 1._wp], shape(tr)))
962
963 ! Translation matrix to move center back to original position
964 t_back = transpose(reshape([1._wp, 0._wp, 0._wp, center(1), 0._wp, 1._wp, 0._wp, center(2), 0._wp, 0._wp, 1._wp, &
965 & center(3), 0._wp, 0._wp, 0._wp, 1._wp], shape(tr)))
966
967 out_matrix = matmul(tr, matmul(t_back, matmul(ry, matmul(rx, matmul(rz, matmul(sc, t_to_origin))))))
968 else
969 out_matrix = matmul(ry, matmul(rx, rz))
970 end if
971
972 end function f_create_transform_matrix
973
974 !> Transform a vector by a matrix.
975 subroutine s_transform_vec(vec, matrix)
976
977 real(wp), dimension(1:3), intent(inout) :: vec
978 real(wp), dimension(1:4,1:4), intent(in) :: matrix
979 real(wp), dimension(1:4) :: tmp
980
981 tmp = matmul(matrix, [vec(1), vec(2), vec(3), 1._wp])
982 vec = tmp(1:3)
983
984 end subroutine s_transform_vec
985
986 !> Transform a triangle by a matrix, one vertex at a time.
987 subroutine s_transform_triangle(triangle, matrix, matrix_n)
988
989 type(t_triangle), intent(inout) :: triangle
990 real(wp), dimension(1:4,1:4), intent(in) :: matrix, matrix_n
991 integer :: i
992
993 do i = 1, 3
994 call s_transform_vec(triangle%v(i,:), matrix)
995 end do
996
997 call s_transform_vec(triangle%n(1:3), matrix_n)
998
999 end subroutine s_transform_triangle
1000
1001 !> Transform a model by a matrix, one triangle at a time.
1002 subroutine s_transform_model(model, matrix, matrix_n)
1003
1004 type(t_model), intent(inout) :: model
1005 real(wp), dimension(1:4,1:4), intent(in) :: matrix, matrix_n
1006 integer :: i
1007
1008 do i = 1, size(model%trs)
1009 call s_transform_triangle(model%trs(i), matrix, matrix_n)
1010 end do
1011
1012 end subroutine s_transform_model
1013
1014 !> Create a bounding box for a model.
1015 function f_create_bbox(model) result(bbox)
1016
1017 type(t_model), intent(in) :: model
1018 type(t_bbox) :: bbox
1019 integer :: i, j
1020
1021 if (size(model%trs) == 0) then
1022 bbox%min = 0._wp
1023 bbox%max = 0._wp
1024 return
1025 end if
1026
1027 bbox%min = model%trs(1)%v(1,:)
1028 bbox%max = model%trs(1)%v(1,:)
1029
1030 do i = 1, size(model%trs)
1031 do j = 1, 3
1032 bbox%min = min(bbox%min, model%trs(i)%v(j,:))
1033 bbox%max = max(bbox%max, model%trs(i)%v(j,:))
1034 end do
1035 end do
1036
1037 end function f_create_bbox
1038
1039 !> Perform XOR on lhs and rhs.
1040 elemental function f_xor(lhs, rhs) result(res)
1041
1042 logical, intent(in) :: lhs, rhs
1043 logical :: res
1044
1045 res = (lhs .and. .not. rhs) .or. (.not. lhs .and. rhs)
1046
1047 end function f_xor
1048
1049 !> Convert a logical to 1 or 0.
1050 elemental function f_logical_to_int(predicate) result(int)
1051
1052 logical, intent(in) :: predicate
1053 integer :: int
1054
1055 if (predicate) then
1056 int = 1
1057 else
1058 int = 0
1059 end if
1060
1061 end function f_logical_to_int
1062
1063 !> Real spherical harmonic Y_lm(theta, phi). theta = polar angle from +z (acos(z/r)), phi = atan2(y,x). Uses associated Legendre
1064 !! P_l^|m|(cos theta). Standard normalisation.
1065 function real_ylm(theta, phi, l, m) result(Y)
1066
1067 integer, intent(in) :: l, m
1068 real(wp), intent(in) :: theta, phi
1069 real(wp) :: y, x, prefac
1070 integer :: m_abs
1071
1072 m_abs = abs(m)
1073 if (m_abs > l) then
1074 y = 0._wp
1075 return
1076 end if
1077 x = cos(theta)
1078 prefac = sqrt((2*l + 1)*real(factorial(l - m_abs), wp)/real(factorial(l + m_abs), wp)/(4._wp*pi))
1079 if (m == 0) then
1080 y = prefac*associated_legendre(x, l, 0)
1081 else if (m > 0) then
1082 y = prefac*sqrt(2._wp)*associated_legendre(x, l, m_abs)*cos(m*phi)
1083 else
1084 y = prefac*sqrt(2._wp)*associated_legendre(x, l, m_abs)*sin(m_abs*phi)
1085 end if
1086
1087 end function real_ylm
1088
1089 !> Associated Legendre polynomial P_l^m(x) (Ferrers function, Condon-Shortley phase). Valid for integer l >= 0, 0 <= m <= l, and
1090 !! x in [-1,1]. Returns 0 for |m| > l or l < 0. Formulas: DLMF 14.10.3 (recurrence in degree), Wikipedia "Associated Legendre
1091 !! polynomials" (P_l^l and P_l^{l-1} identities). Recurrence: (l-m)P_l^m = (2l-1)x P_{l-1}^m - (l+m-1)P_{l-2}^m.
1092 !! @param x argument (typically cos(theta)), should be in [-1,1]
1093 !! @param l degree (>= 0)
1094 !! @param m_order order (0 <= m_order <= l)
1095 recursive function associated_legendre(x, l, m_order) result(result_P)
1096
1097 integer, intent(in) :: l, m_order
1098 real(wp), intent(in) :: x
1099 real(wp) :: result_p
1100 real(wp) :: one_minus_x2
1101
1102 ! Out-of-domain: P_l^m = 0 for |m| > l or l < 0 (standard convention)
1103
1104 if (l < 0 .or. m_order < 0 .or. m_order > l) then
1105 result_p = 0._wp
1106 return
1107 end if
1108
1109 if (m_order <= 0 .and. l <= 0) then
1110 result_p = 1._wp
1111 else if (l == 1 .and. m_order <= 0) then
1112 result_p = x
1113 else if (l == 1 .and. m_order == 1) then
1114 one_minus_x2 = max(0._wp, 1._wp - x**2)
1115 result_p = -sqrt(one_minus_x2)
1116 else if (m_order == l) then
1117 ! P_l^l(x) = (-1)^l (2l-1)!! (1-x^2)^(l/2). Use real exponent for odd l
1118 one_minus_x2 = max(0._wp, 1._wp - x**2)
1119 result_p = (-1)**l*real(double_factorial(2*l - 1), wp)*one_minus_x2**(0.5_wp*real(l, wp))
1120 else if (m_order == l - 1) then
1121 result_p = x*(2*l - 1)*associated_legendre(x, l - 1, l - 1)
1122 else
1123 result_p = ((2*l - 1)*x*associated_legendre(x, l - 1, m_order) - (l + m_order - 1)*associated_legendre(x, l - 2, &
1124 & m_order))/(l - m_order)
1125 end if
1126
1127 end function associated_legendre
1128
1129 !> Calculate the double factorial of an integer
1130 elemental function double_factorial(n_in) result(R_result)
1131
1132 integer, intent(in) :: n_in
1133 integer, parameter :: int64_kind = selected_int_kind(18) !< 18 bytes for 64-bit integer
1134 integer(kind=int64_kind) :: r_result
1135 integer :: i
1136
1137 r_result = product((/(i, i=n_in, 1, -2)/))
1138
1139 end function double_factorial
1140
1141 !> Calculate the factorial of an integer
1142 elemental function factorial(n_in) result(R_result)
1143
1144 integer, intent(in) :: n_in
1145 integer, parameter :: int64_kind = selected_int_kind(18) !< 18 bytes for 64-bit integer
1146 integer(kind=int64_kind) :: r_result
1147 integer :: i
1148
1149 r_result = product((/(i, i=n_in, 1, -1)/))
1150
1151 end function factorial
1152
1153 !> Calculate a smooth cut-on function that is zero for x values smaller than zero and goes to one, for generating smooth initial
1154 !! conditions
1155 function f_cut_on(x, eps) result(fx)
1156
1157 real(wp), intent(in) :: x, eps
1158 real(wp) :: fx
1159
1160 fx = 1 - f_gx(x/eps)/(f_gx(x/eps) + f_gx(1 - x/eps))
1161
1162 end function f_cut_on
1163
1164 !> Calculate a smooth cut-off function that is one for x values smaller than zero and goes to zero, for generating smooth
1165 !! initial conditions
1166 function f_cut_off(x, eps) result(fx)
1167
1168 real(wp), intent(in) :: x, eps
1169 real(wp) :: fx
1170
1171 fx = f_gx(x/eps)/(f_gx(x/eps) + f_gx(1 - x/eps))
1172
1173 end function f_cut_off
1174
1175 !> Helper function for f_cut_on and f_cut_off
1176 function f_gx(x) result(gx)
1177
1178 real(wp), intent(in) :: x
1179 real(wp) :: gx
1180
1181 if (x > 0) then
1182 gx = exp(-1._wp/x)
1183 else
1184 gx = 0._wp
1185 end if
1186
1187 end function f_gx
1188
1189 !> Downsample conservative variable fields by a factor of 3 in each direction using volume averaging.
1190 subroutine s_downsample_data(q_cons_vf, q_cons_temp, m_ds, n_ds, p_ds, m_glb_ds, n_glb_ds, p_glb_ds)
1191
1192 type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf, q_cons_temp
1193
1194 ! Down sampling variables
1195 integer :: i, j, k, l
1196 integer :: ix, iy, iz, x_id, y_id, z_id
1197 integer, intent(inout) :: m_ds, n_ds, p_ds, m_glb_ds, n_glb_ds, p_glb_ds
1198
1199 m_ds = int((m + 1)/3) - 1
1200 n_ds = int((n + 1)/3) - 1
1201 p_ds = int((p + 1)/3) - 1
1202
1203 m_glb_ds = int((m_glb + 1)/3) - 1
1204 n_glb_ds = int((n_glb + 1)/3) - 1
1205 p_glb_ds = int((p_glb + 1)/3) - 1
1206
1207 do l = -1, p_ds + 1
1208 do k = -1, n_ds + 1
1209 do j = -1, m_ds + 1
1210 x_id = 3*j + 1
1211 y_id = 3*k + 1
1212 z_id = 3*l + 1
1213 do i = 1, sys_size
1214 q_cons_temp(i)%sf(j, k, l) = 0
1215
1216 do iz = -1, 1
1217 do iy = -1, 1
1218 do ix = -1, 1
1219 q_cons_temp(i)%sf(j, k, l) = q_cons_temp(i)%sf(j, k, &
1220 & l) + (1._wp/27._wp)*q_cons_vf(i)%sf(x_id + ix, y_id + iy, z_id + iz)
1221 end do
1222 end do
1223 end do
1224 end do
1225 end do
1226 end do
1227 end do
1228
1229 end subroutine s_downsample_data
1230
1231 !> Upsample conservative variable fields from a coarsened grid back to the original resolution using interpolation.
1232 subroutine s_upsample_data(q_cons_vf, q_cons_temp)
1233
1234 type(scalar_field), intent(inout), dimension(sys_size) :: q_cons_vf, q_cons_temp
1235 integer :: i, j, k, l
1236 integer :: ix, iy, iz
1237 integer :: x_id, y_id, z_id
1238 real(wp), dimension(4) :: temp
1239
1240 do l = 0, p
1241 do k = 0, n
1242 do j = 0, m
1243 do i = 1, sys_size
1244 ix = int(j/3._wp)
1245 iy = int(k/3._wp)
1246 iz = int(l/3._wp)
1247
1248 x_id = j - int(3*ix) - 1
1249 y_id = k - int(3*iy) - 1
1250 z_id = l - int(3*iz) - 1
1251
1252 temp(1) = (2._wp/3._wp)*q_cons_temp(i)%sf(ix, iy, iz) + (1._wp/3._wp)*q_cons_temp(i)%sf(ix + x_id, iy, iz)
1253 temp(2) = (2._wp/3._wp)*q_cons_temp(i)%sf(ix, iy + y_id, iz) + (1._wp/3._wp)*q_cons_temp(i)%sf(ix + x_id, &
1254 & iy + y_id, iz)
1255 temp(3) = (2._wp/3._wp)*temp(1) + (1._wp/3._wp)*temp(2)
1256
1257 temp(1) = (2._wp/3._wp)*q_cons_temp(i)%sf(ix, iy, iz + z_id) + (1._wp/3._wp)*q_cons_temp(i)%sf(ix + x_id, &
1258 & iy, iz + z_id)
1259 temp(2) = (2._wp/3._wp)*q_cons_temp(i)%sf(ix, iy + y_id, &
1260 & iz + z_id) + (1._wp/3._wp)*q_cons_temp(i)%sf(ix + x_id, iy + y_id, iz + z_id)
1261 temp(4) = (2._wp/3._wp)*temp(1) + (1._wp/3._wp)*temp(2)
1262
1263 q_cons_vf(i)%sf(j, k, l) = (2._wp/3._wp)*temp(3) + (1._wp/3._wp)*temp(4)
1264 end do
1265 end do
1266 end do
1267 end do
1268
1269 end subroutine s_upsample_data
1270
1271end module m_helper
q_cons_vf
type(scalar_field), dimension(sys_size), intent(inout) q_cons_vf
Definition m_phase_change.fpp.f90:988
k
integer, intent(in) k
Definition m_phase_change.fpp.f90:990
j
integer, intent(in) j
Definition m_phase_change.fpp.f90:990
l
integer, intent(in) l
Definition m_phase_change.fpp.f90:990
m_derived_types
Shared derived types for field data, patch geometry, bubble dynamics, and MPI I/O structures.
Definition m_derived_types.fpp.f90:323
m_global_parameters
Global parameters for the post-process: domain geometry, equation of state, and output database setti...
Definition m_global_parameters.fpp.f90:19
m_global_parameters::weight
real(wp), dimension(:), allocatable weight
Definition m_global_parameters.fpp.f90:146
m_global_parameters::r0
real(wp), dimension(:), allocatable r0
Definition m_global_parameters.fpp.f90:146
m_helper
Utility routines for bubble model setup, coordinate transforms, array sampling, and special functions...
Definition m_helper.fpp.f90:334
m_helper::s_comp_n_from_prim
subroutine, public s_comp_n_from_prim(vftmp, rtmp, ntmp, weights)
Computes the bubble number density n from the primitive variables.
Definition m_helper.fpp.f90:352
m_helper::s_cross_product
subroutine, public s_cross_product(a, b, c)
Compute the cross product c = a x b of two 3D vectors.
Definition m_helper.fpp.f90:899
m_helper::modmul
integer function, public modmul(a)
Compute a modular multiplication step for the linear congruential pseudo-random number generator.
Definition m_helper.fpp.f90:871
m_helper::s_initialize_nonpoly
impure subroutine, public s_initialize_nonpoly()
Initializes non-polydisperse bubble modeling.
Definition m_helper.fpp.f90:669
m_helper::associated_legendre
recursive real(wp) function, public associated_legendre(x, l, m_order)
Associated Legendre polynomial P_l^m(x) (Ferrers function, Condon-Shortley phase)....
Definition m_helper.fpp.f90:1096
m_helper::s_transform_triangle
subroutine, public s_transform_triangle(triangle, matrix, matrix_n)
Transform a triangle by a matrix, one vertex at a time.
Definition m_helper.fpp.f90:988
m_helper::f_unit_vector
real(wp) function, dimension(3), public f_unit_vector(theta, eta)
Generate a unit vector uniformly distributed on the sphere from two random parameters.
Definition m_helper.fpp.f90:815
m_helper::f_cut_on
real(wp) function, public f_cut_on(x, eps)
Calculate a smooth cut-on function that is zero for x values smaller than zero and goes to one,...
Definition m_helper.fpp.f90:1156
m_helper::s_transform_model
subroutine, public s_transform_model(model, matrix, matrix_n)
Transform a model by a matrix, one triangle at a time.
Definition m_helper.fpp.f90:1003
m_helper::f_create_bbox
type(t_bbox) function, public f_create_bbox(model)
Create a bounding box for a model.
Definition m_helper.fpp.f90:1016
m_helper::f_gx
real(wp) function f_gx(x)
Helper function for f_cut_on and f_cut_off.
Definition m_helper.fpp.f90:1177
m_helper::s_simpson
subroutine, public s_simpson(local_weight, local_r0)
Compute the Simpson weights for quadrature.
Definition m_helper.fpp.f90:747
m_helper::double_factorial
elemental integer(kind=int64_kind) function, public double_factorial(n_in)
Calculate the double factorial of an integer.
Definition m_helper.fpp.f90:1131
m_helper::s_initialize_bubble_vars
impure subroutine s_initialize_bubble_vars()
Set bubble physical parameters and nondimensional numbers from the input configuration.
Definition m_helper.fpp.f90:616
m_helper::factorial
elemental integer(kind=int64_kind) function, public factorial(n_in)
Calculate the factorial of an integer.
Definition m_helper.fpp.f90:1143
m_helper::s_upsample_data
subroutine, public s_upsample_data(q_cons_vf, q_cons_temp)
Upsample conservative variable fields from a coarsened grid back to the original resolution using int...
Definition m_helper.fpp.f90:1233
m_helper::f_create_transform_matrix
real(wp) function, dimension(1:4, 1:4), public f_create_transform_matrix(param, center)
Create a transformation matrix.
Definition m_helper.fpp.f90:938
m_helper::f_cut_off
real(wp) function, public f_cut_off(x, eps)
Calculate a smooth cut-off function that is one for x values smaller than zero and goes to zero,...
Definition m_helper.fpp.f90:1167
m_helper::s_transform_vec
subroutine, public s_transform_vec(vec, matrix)
Transform a vector by a matrix.
Definition m_helper.fpp.f90:976
m_helper::s_initialize_bubbles_model
impure subroutine, public s_initialize_bubbles_model()
Initialize bubble model arrays for Euler or Lagrangian bubbles with polytropic or non-polytropic gas.
Definition m_helper.fpp.f90:440
m_helper::s_print_2d_array
impure subroutine, public s_print_2d_array(a, div)
Print a 2D real array to standard output, optionally dividing each element by a given scalar.
Definition m_helper.fpp.f90:410
m_helper::s_transcoeff
elemental subroutine, public s_transcoeff(omega, peclet, re_trans, im_trans)
Computes the transfer coefficient for the non-polytropic bubble compression process.
Definition m_helper.fpp.f90:717
m_helper::s_prng
subroutine, public s_prng(var, seed)
Generate a pseudo-random number between 0 and 1 using a linear congruential generator.
Definition m_helper.fpp.f90:845
m_helper::f_cross
pure real(wp) function, dimension(3), public f_cross(a, b)
Compute the cross product of two vectors.
Definition m_helper.fpp.f90:787
m_helper::s_downsample_data
subroutine, public s_downsample_data(q_cons_vf, q_cons_temp, m_ds, n_ds, p_ds, m_glb_ds, n_glb_ds, p_glb_ds)
Downsample conservative variable fields by a factor of 3 in each direction using volume averaging.
Definition m_helper.fpp.f90:1191
m_helper::f_xor
elemental logical function, public f_xor(lhs, rhs)
Perform XOR on lhs and rhs.
Definition m_helper.fpp.f90:1041
m_helper::s_comp_n_from_cons
subroutine, public s_comp_n_from_cons(vftmp, nrtmp, ntmp, weights)
Compute the bubble number density from the conservative void fraction and weighted bubble radii.
Definition m_helper.fpp.f90:381
m_helper::real_ylm
real(wp) function, public real_ylm(theta, phi, l, m)
Real spherical harmonic Y_lm(theta, phi). theta = polar angle from +z (acos(z/r)),...
Definition m_helper.fpp.f90:1066
m_helper::f_logical_to_int
elemental integer function, public f_logical_to_int(predicate)
Convert a logical to 1 or 0.
Definition m_helper.fpp.f90:1051
m_helper::s_int_to_str
elemental subroutine, public s_int_to_str(i, res)
Convert an integer to its trimmed string representation.
Definition m_helper.fpp.f90:736
m_helper::s_swap
elemental subroutine, public s_swap(lhs, rhs)
Swap two real numbers.
Definition m_helper.fpp.f90:926