MFC
Exascale flow solver
Loading...
Searching...
No Matches
m_phase_change.fpp.f90
Go to the documentation of this file.
1# 1 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
2!>
3!! @file
4!! @brief Contains module m_phase_change
5
6# 1 "/home/runner/work/MFC/MFC/src/common/include/case.fpp" 1
7! This file exists so that Fypp can be run without generating case.fpp files for
8! each target. This is useful when generating documentation, for example. This
9! should also let MFC be built with CMake directly, without invoking mfc.sh.
10
11! For pre-process.
12# 9 "/home/runner/work/MFC/MFC/src/common/include/case.fpp"
13
14! For moving immersed boundaries in simulation
15# 14 "/home/runner/work/MFC/MFC/src/common/include/case.fpp"
16# 6 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp" 2
17# 1 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 1
18# 1 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 1
19# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
20# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
21# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
22# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
23# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
24# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
25
26# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
27# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
28# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
29
30# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
31
32# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
33
34# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
35
36# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
37
38# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
39
40# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
41
42# 120 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
43! New line at end of file is required for FYPP
44# 2 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
45# 1 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 1
46# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
47# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
48# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
49# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
50# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
51# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
52
53# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
54# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
55# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
56
57# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
58
59# 46 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
60
61# 58 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
62
63# 68 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
64
65# 98 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
66
67# 110 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
68
69# 120 "/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! New line at end of file is required for FYPP
177# 2 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 2
178
179# 7 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
180
181# 17 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
182
183# 22 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
184
185# 27 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
186
187# 32 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
188
189# 37 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
190
191# 42 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
192
193# 47 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
194
195# 52 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
196
197# 57 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
198
199# 62 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
200
201# 73 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
202
203# 78 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
204
205# 83 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
206
207# 88 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
208
209# 103 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
210
211# 131 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
212
213# 160 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
214
215# 175 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
216
217# 193 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
218
219# 215 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
220
221# 244 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
222
223# 259 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
224
225# 269 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
226
227# 278 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
228
229# 294 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
230
231# 304 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
232
233# 311 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
234! New line at end of file is required for FYPP
235# 4 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
236
237! GPU parallel region (scalar reductions, maxval/minval)
238# 23 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
239
240! GPU parallel loop over threads (most common GPU macro)
241# 43 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
242
243! Required closing for GPU_PARALLEL_LOOP
244# 55 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
245
246! Mark routine for device compilation
247# 112 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
248
249! Declare device-resident data
250# 130 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
251
252! Inner loop within a GPU parallel region
253# 145 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
254
255! Scoped GPU data region
256# 164 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
257
258! Host code with device pointers (for MPI with GPU buffers)
259# 193 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
260
261! Allocate device memory (unscoped)
262# 207 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
263
264! Free device memory
265# 219 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
266
267! Atomic operation on device
268# 231 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
269
270! End atomic capture block
271# 242 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
272
273! Copy data between host and device
274# 254 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
275
276! Synchronization barrier
277# 266 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
278
279! Import GPU library module (openacc or omp_lib)
280# 275 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
281
282! Emit code only for AMD compiler
283# 282 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
284
285! Emit code for non-Cray compilers
286# 289 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
287
288! Emit code only for Cray compiler
289# 296 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
290
291! Emit code for non-NVIDIA compilers
292# 303 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
293
294# 305 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
295# 306 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
296! New line at end of file is required for FYPP
297# 2 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 2
298
299# 14 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
300
301! Caution: This macro requires the use of a binding script to set CUDA_VISIBLE_DEVICES, such that we have one GPU device per MPI
302! rank. That's because for both cudaMemAdvise (preferred location) and cudaMemPrefetchAsync we use location = device_id = 0. For an
303! example see misc/nvidia_uvm/bind.sh. NVIDIA unified memory page placement hint
304# 57 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
305
306! Allocate and create GPU device memory
307# 77 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
308
309! Free GPU device memory and deallocate
310# 85 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
311
312! Cray-specific GPU pointer setup for vector fields
313# 109 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
314
315! Cray-specific GPU pointer setup for scalar fields
316# 125 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
317
318! Cray-specific GPU pointer setup for acoustic source spatials
319# 150 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
320
321# 156 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
322
323# 163 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
324! New line at end of file is required for FYPP
325# 7 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp" 2
326
327!> @brief Phase transition relaxation solvers for liquid-vapor flows with cavitation and boiling
329
330#ifndef MFC_POST_PROCESS
333 use m_mpi_proxy
335 use ieee_arithmetic
337
338 implicit none
339
340 private
342
343 !> @name Parameters for the first order transition phase change
344 !> @{
345 integer, parameter :: max_iter = 1e8_wp !< max # of iterations
346 real(wp), parameter :: pcr = 4.94e7_wp !< Critical pressure of water [Pa]
347 real(wp), parameter :: tcr = 385.05_wp + 273.15_wp !< Critical temperature of water [K]
348 real(wp), parameter :: mixm = 1.0e-8_wp !< Mixture mass fraction threshold for triggering phase change
349 integer, parameter :: lp = 1 !< index for the liquid phase of the reacting fluid
350 integer, parameter :: vp = 2 !< index for the vapor phase of the reacting fluid
351 !> @}
352
353 !> @name Gibbs free energy phase change parameters
354 !> @{
355 real(wp) :: a, b, c, d
356 !> @}
357
358
359# 39 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
360#if defined(MFC_OpenACC)
361# 39 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
362!$acc declare create(A, B, C, D)
363# 39 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
364#elif defined(MFC_OpenMP)
365# 39 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
366!$omp declare target (A, B, C, D)
367# 39 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
368#endif
369
370contains
371
372 !> Dispatch to the correct relaxation solver. Replaces the procedure pointer, which CCE is breaking on.
373 impure subroutine s_relaxation_solver(q_cons_vf)
374
375 type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
376 ! This is empty because in current master the procedure pointer was never assigned
377
378 if (.not. (.false.)) then
379# 49 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
380 call s_mpi_abort("m_phase_change.fpp:49: " // "Assertion failed: .false.. " &
381# 49 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
382 & // "s_relaxation_solver called but it currently does nothing")
383# 49 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
384 end if
385
386 end subroutine s_relaxation_solver
387
388 !> Initialize the phase change module by setting saturation curve coefficients for pT- or pTg-equilibrium
390
391 ! Saturation curve coefficients via stiffened gas EOS. Saurel et al. JCP (2008), Le Metayer et al. JFE (2004)
392 a = (gs_min(lp)*cvs(lp) - gs_min(vp)*cvs(vp) + qvps(vp) - qvps(lp))/((gs_min(vp) - 1.0_wp)*cvs(vp))
393
394 b = (qvs(lp) - qvs(vp))/((gs_min(vp) - 1.0_wp)*cvs(vp))
395
396 c = (gs_min(vp)*cvs(vp) - gs_min(lp)*cvs(lp))/((gs_min(vp) - 1.0_wp)*cvs(vp))
397
398 d = ((gs_min(lp) - 1.0_wp)*cvs(lp))/((gs_min(vp) - 1.0_wp)*cvs(vp))
399
401
402 !> Apply pT- or pTg-equilibrium relaxation with mass depletion based on the incoming state conditions.
403 subroutine s_infinite_relaxation_k(q_cons_vf)
404
405 type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
406 real(wp) :: ps, psov, pssl !< equilibrium pressure for mixture, overheated vapor, and subcooled liquid
407 real(wp) :: ts, tsov, tssl, tsatov, tsatsl !< Equilibrium and saturation temperatures
408 real(wp) :: rhoe, dyne, rhos !< total internal energy, kinetic energy, and total entropy
409 real(wp) :: rho, rm, m1, m2, mct !< total density, total reacting mass, individual reacting masses
410 real(wp) :: tvf !< total volume fraction
411 ! $:GPU_DECLARE(create='[pS,pSOV,pSSL,TS,TSOV,TSSL,TSatOV,TSatSL]')
412 ! $:GPU_DECLARE(create='[rhoe,dynE,rhos,rho,rM,m1,m2,MCT,TvF]')
413
414# 82 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
415 real(wp), dimension(num_fluids) :: p_infov, p_infpt, p_infsl, sk, hk, gk, ek, rhok
416# 84 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
417 ! $:GPU_DECLARE(create='[p_infOV,p_infpT,p_infSL,sk,hk,gk,ek,rhok]')
418
419 !> Generic loop iterators
420 integer :: i, j, k, l
421
422#ifdef _CRAYFTN
423#ifdef MFC_OpenACC
424 ! CCE 19 IPA workaround: prevent bring_routine_resident SIGSEGV DIR$ NOINLINE s_infinite_pt_relaxation_k DIR$ NOINLINE
425 ! s_infinite_ptg_relaxation_k DIR$ NOINLINE s_correct_partial_densities DIR$ NOINLINE s_TSat
426#endif
427#endif
428
429 ! starting equilibrium solver
430
431
432# 98 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
433
434# 98 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
435#if defined(MFC_OpenACC)
436# 98 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
437!$acc parallel loop collapse(3) gang vector default(present) private(i, j, k, l, p_infOV, p_infpT, p_infSL, sk, hk, gk, ek, rhok, pS, pSOV, pSSL, TS, TSOV, TSatOV, TSatSL, TSSL, rhoe, dynE, rhos, rho, rM, m1, m2, MCT, TvF)
438# 98 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
439#elif defined(MFC_OpenMP)
440# 98 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
441
442# 98 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
443
444# 98 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
445
446# 98 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
447!$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, l, p_infOV, p_infpT, p_infSL, sk, hk, gk, ek, rhok, pS, pSOV, pSSL, TS, TSOV, TSatOV, TSatSL, TSSL, rhoe, dynE, rhos, rho, rM, m1, m2, MCT, TvF)
448# 98 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
449#endif
450# 100 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
451 do j = 0, m
452 do k = 0, n
453 do l = 0, p
454 rho = 0.0_wp; tvf = 0.0_wp
455
456# 104 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
457#if defined(MFC_OpenACC)
458# 104 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
459!$acc loop seq
460# 104 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
461#elif defined(MFC_OpenMP)
462# 104 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
463
464# 104 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
465#endif
466 do i = 1, num_fluids
467 ! Mixture density
468 rho = rho + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)
469
470 ! Total Volume Fraction
471 tvf = tvf + q_cons_vf(i + eqn_idx%adv%beg - 1)%sf(j, k, l)
472 end do
473
474 ! calculating the total reacting mass for the phase change process. By hypothesis, this should not change
475 ! throughout the phase-change process.
476 rm = q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) + q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l)
477
478 ! correcting negative (reacting) mass fraction values in case they happen
479 call s_correct_partial_densities(mct, q_cons_vf, rm, j, k, l)
480
481 ! fixing m1 and m2 AFTER correcting the partial densities. Note that these values must be stored for the phase
482 ! change process that will happen a posteriori
483 m1 = q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l)
484
485 m2 = q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l)
486
487 ! kinetic energy as an auxiliary variable to the calculation of the total internal energy
488 dyne = 0.0_wp
489
490# 128 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
491#if defined(MFC_OpenACC)
492# 128 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
493!$acc loop seq
494# 128 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
495#elif defined(MFC_OpenMP)
496# 128 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
497
498# 128 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
499#endif
500 do i = eqn_idx%mom%beg, eqn_idx%mom%end
501 dyne = dyne + 5.0e-1_wp*q_cons_vf(i)%sf(j, k, l)**2/rho
502 end do
503
504 ! calculating the total energy that MUST be preserved throughout the pT- and pTg-relaxation procedures at each
505 ! of the cells. The internal energy is calculated as the total energy minus the kinetic energy to preserved its
506 ! value at sharp interfaces
507 rhoe = q_cons_vf(eqn_idx%E)%sf(j, k, l) - dyne
508
509 ! Calling pT-equilibrium for either finishing phase-change module, or as an IC for the pTg-equilibrium for this
510 ! case, MFL cannot be either 0 or 1, so I chose it to be 2
511 call s_infinite_pt_relaxation_k(j, k, l, 2, ps, p_infpt, q_cons_vf, rhoe, ts)
512
513 ! Check if pTg-equilibrium needed; only partial densities require updating
514 if ((relax_model == 6) .and. ((q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, &
515 & l) > mixm*rm) .and. (q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, &
516 & l) > mixm*rm)) .and. (ps < pcr) .and. (ts < tcr)) then
517 ! Checking if phase change is needed, by checking whether the final solution is either subcoooled liquid or
518 ! overheated vapor.
519
520 ! overheated vapor case depleting the mass of liquid
521 q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) = mixm*rm
522
523 ! transferring the total mass to vapor
524 q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) = (1.0_wp - mixm)*rm
525
526 ! calling pT-equilibrium for overheated vapor, which is MFL = 0
527 call s_infinite_pt_relaxation_k(j, k, l, 0, psov, p_infov, q_cons_vf, rhoe, tsov)
528
529 ! calculating Saturation temperature
530 call s_tsat(psov, tsatov, tsov)
531
532 ! subcooled liquid case transferring the total mass to liquid
533 q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) = (1.0_wp - mixm)*rm
534
535 ! depleting the mass of vapor
536 q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) = mixm*rm
537
538 ! calling pT-equilibrium for subcooled liquid, which is MFL = 1
539 call s_infinite_pt_relaxation_k(j, k, l, 1, pssl, p_infsl, q_cons_vf, rhoe, tssl)
540
541 ! calculating Saturation temperature
542 call s_tsat(pssl, tsatsl, tssl)
543
544 ! checking the conditions for overheated vapor and subcooled liquide
545 if (tsov > tsatov) then
546 ! Assigning pressure
547 ps = psov
548
549 ! Assigning Temperature
550 ts = tsov
551
552 ! correcting the liquid partial density
553 q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) = mixm*rm
554
555 ! correcting the vapor partial density
556 q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) = (1.0_wp - mixm)*rm
557 else if (tssl < tsatsl) then
558 ! Assigning pressure
559 ps = pssl
560
561 ! Assigning Temperature
562 ts = tssl
563
564 ! correcting the liquid partial density
565 q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) = (1.0_wp - mixm)*rm
566
567 ! correcting the vapor partial density
568 q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) = mixm*rm
569 else
570 ! returning partial pressures to what they were from the homogeneous solver liquid
571 q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) = m1
572
573 ! vapor
574 q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) = m2
575
576 ! calling the pTg-equilibrium solver
577 call s_infinite_ptg_relaxation_k(j, k, l, ps, p_infpt, rhoe, q_cons_vf, ts)
578 end if
579 end if
580
581 ! Calculations AFTER equilibrium
582
583
584# 212 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
585#if defined(MFC_OpenACC)
586# 212 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
587!$acc loop seq
588# 212 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
589#elif defined(MFC_OpenMP)
590# 212 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
591
592# 212 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
593#endif
594 do i = 1, num_fluids
595 ! entropy
596 sk(i) = cvs(i)*log((ts**gs_min(i))/((ps + ps_inf(i))**(gs_min(i) - 1.0_wp))) + qvps(i)
597
598 ! enthalpy
599 hk(i) = gs_min(i)*cvs(i)*ts + qvs(i)
600
601 ! Gibbs-free energy
602 gk(i) = hk(i) - ts*sk(i)
603
604 ! densities
605 rhok(i) = (ps + ps_inf(i))/((gs_min(i) - 1)*cvs(i)*ts)
606
607 ! internal energy
608 ek(i) = (ps + gs_min(i)*ps_inf(i))/(ps + ps_inf(i))*cvs(i)*ts + qvs(i)
609 end do
610
611 ! calculating volume fractions, internal energies, and total entropy
612 rhos = 0.0_wp
613
614# 232 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
615#if defined(MFC_OpenACC)
616# 232 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
617!$acc loop seq
618# 232 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
619#elif defined(MFC_OpenMP)
620# 232 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
621
622# 232 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
623#endif
624 do i = 1, num_fluids
625 ! volume fractions
626 q_cons_vf(i + eqn_idx%adv%beg - 1)%sf(j, k, l) = q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)/rhok(i)
627
628 ! alpha*rho*e
629 if (model_eqns == 3) then
630 q_cons_vf(i + eqn_idx%int_en%beg - 1)%sf(j, k, l) = q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, &
631 & l)*ek(i)
632 end if
633
634 ! Total entropy
635 rhos = rhos + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*sk(i)
636 end do
637 end do
638 end do
639 end do
640
641# 249 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
642#if defined(MFC_OpenACC)
643# 249 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
644!$acc end parallel loop
645# 249 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
646#elif defined(MFC_OpenMP)
647# 249 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
648
649# 249 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
650!$omp end target teams loop
651# 249 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
652#endif
653
654 end subroutine s_infinite_relaxation_k
655
656 !> Apply pT-equilibrium relaxation for N fluids
657 !! @param MFL flag: 0=gas, 1=liquid, 2=mixture
658 subroutine s_infinite_pt_relaxation_k(j, k, l, MFL, pS, p_infpT, q_cons_vf, rhoe, TS)
659
660
661# 257 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
662#ifdef _CRAYFTN
663# 257 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
664#if MFC_OpenACC
665# 257 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
666!$acc routine seq
667# 257 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
668#elif MFC_OpenMP
669# 257 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
670
671# 257 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
672
673# 257 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
674!$omp declare target device_type(any)
675# 257 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
676#else
677# 257 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
678!DIR$ NOINLINE s_infinite_pt_relaxation_k
679# 257 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
680#endif
681# 257 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
682#elif MFC_OpenACC
683# 257 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
684!$acc routine seq
685# 257 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
686#elif MFC_OpenMP
687# 257 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
688
689# 257 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
690
691# 257 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
692!$omp declare target device_type(any)
693# 257 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
694#endif
695
696 ! initializing variables
697 integer, intent(in) :: j, k, l, MFL
698 real(wp), intent(out) :: pS
699 real(wp), dimension(1:), intent(out) :: p_infpT
700 type(scalar_field), dimension(sys_size), intent(in) :: q_cons_vf
701 real(wp), intent(in) :: rhoe
702 real(wp), intent(out) :: TS
703 real(wp) :: gp, gpp, hp, pO, mCP, mQ !< variables for the Newton Solver
704 real(wp) :: p_infpT_sum
705 integer :: i, ns !< generic loop iterators
706 ! auxiliary variables for the pT-equilibrium solver
707 mcp = 0.0_wp; mq = 0.0_wp; p_infpt_sum = 0._wp
708
709# 271 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
710#if defined(MFC_OpenACC)
711# 271 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
712!$acc loop seq
713# 271 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
714#elif defined(MFC_OpenMP)
715# 271 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
716
717# 271 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
718#endif
719 do i = 1, num_fluids
720 p_infpt(i) = ps_inf(i)
721 p_infpt_sum = p_infpt_sum + abs(p_infpt(i))
722 end do
723 ! Performing tests before initializing the pT-equilibrium
724
725# 277 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
726#if defined(MFC_OpenACC)
727# 277 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
728!$acc loop seq
729# 277 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
730#elif defined(MFC_OpenMP)
731# 277 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
732
733# 277 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
734#endif
735 do i = 1, num_fluids
736 ! sum of the total alpha*rho*cp of the system
737 mcp = mcp + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*cvs(i)*gs_min(i)
738
739 ! sum of the total alpha*rho*q of the system
740 mq = mq + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*qvs(i)
741 end do
742
743# 294 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
744
745 ! Checking energy constraint
746 if ((rhoe - mq - minval(p_infpt)) < 0.0_wp) then
747 if ((mfl == 0) .or. (mfl == 1)) then
748 ! Assigning zero values for mass depletion cases pressure
749 ps = 0.0_wp
750
751 ! temperature
752 ts = 0.0_wp
753
754 return
755 end if
756 end if
757
758 ! calculating initial estimate for pressure in the pT-relaxation procedure. I will also use this variable to iterate over
759 ! the Newton's solver
760 po = 0.0_wp
761
762 ! Maybe improve this condition afterwards. As long as the initial guess is in between -min(ps_inf) and infinity, a solution
763 ! should be able to be found.
764 ps = 1.0e4_wp
765
766 ! Newton Solver for the pT-equilibrium
767 ns = 0
768 ! change this relative error metric. 1.e4_wp is just arbitrary
769 do while ((abs(ps - po) > palpha_eps) .and. (abs((ps - po)/po) > palpha_eps/1.e4_wp) .or. (ns == 0))
770 ! increasing counter
771 ns = ns + 1
772
773 ! updating old pressure
774 po = ps
775
776 ! updating functions used in the Newton's solver
777 gpp = 0.0_wp; gp = 0.0_wp; hp = 0.0_wp
778
779# 328 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
780#if defined(MFC_OpenACC)
781# 328 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
782!$acc loop seq
783# 328 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
784#elif defined(MFC_OpenMP)
785# 328 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
786
787# 328 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
788#endif
789 do i = 1, num_fluids
790 gp = gp + (gs_min(i) - 1.0_wp)*q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, &
791 & l)*cvs(i)*(rhoe + ps - mq)/(mcp*(ps + p_infpt(i)))
792
793 gpp = gpp + (gs_min(i) - 1.0_wp)*q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, &
794 & l)*cvs(i)*(p_infpt(i) - rhoe + mq)/(mcp*(ps + p_infpt(i))**2)
795 end do
796
797 hp = 1.0_wp/(rhoe + ps - mq) + 1.0_wp/(ps + minval(p_infpt))
798
799 ! updating common pressure for the newton solver
800 ps = po + ((1.0_wp - gp)/gpp)/(1.0_wp - (1.0_wp - gp + abs(1.0_wp - gp))/(2.0_wp*gpp)*hp)
801 end do
802
803 ! common temperature
804 ts = (rhoe + ps - mq)/mcp
805
806 end subroutine s_infinite_pt_relaxation_k
807
808 !> Apply pTg-equilibrium relaxation for N fluids under pT and 2 fluids under pTg-equilibrium. There is a final common p and T
809 !! during relaxation
810 subroutine s_infinite_ptg_relaxation_k(j, k, l, pS, p_infpT, rhoe, q_cons_vf, TS)
811
812
813# 352 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
814#ifdef _CRAYFTN
815# 352 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
816#if MFC_OpenACC
817# 352 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
818!$acc routine seq
819# 352 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
820#elif MFC_OpenMP
821# 352 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
822
823# 352 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
824
825# 352 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
826!$omp declare target device_type(any)
827# 352 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
828#else
829# 352 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
830!DIR$ NOINLINE s_infinite_ptg_relaxation_k
831# 352 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
832#endif
833# 352 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
834#elif MFC_OpenACC
835# 352 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
836!$acc routine seq
837# 352 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
838#elif MFC_OpenMP
839# 352 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
840
841# 352 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
842
843# 352 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
844!$omp declare target device_type(any)
845# 352 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
846#endif
847
848 integer, intent(in) :: j, k, l
849 real(wp), intent(inout) :: pS
850 real(wp), dimension(1:), intent(in) :: p_infpT
851 real(wp), intent(in) :: rhoe
852 type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
853 real(wp), intent(inout) :: TS
854# 363 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
855 real(wp), dimension(num_fluids) :: p_infpTg !< stiffness for the participating fluids for pTg-equilibrium
856# 365 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
857 real(wp), dimension(2, 2) :: Jac, InvJac, TJac !< matrices for the Newton Solver
858 real(wp), dimension(2) :: R2D, DeltamP !< residual and correction array
859 real(wp) :: Om !< underrelaxation factor
860 real(wp) :: mCP, mCPD, mCVGP, mCVGP2, mQ, mQD !< auxiliary variables for the pTg-solver
861 real(wp) :: ml, mT, dFdT, dTdm, dTdp
862
863 !> Generic loop iterators
864 integer :: i, ns
865 ! pTg-equilibrium solution procedure Newton Solver parameters counter
866 ns = 0
867
868 ! Relaxation factor
869 om = 1.0e-3_wp
870
871 p_infptg = p_infpt
872
873 if (((ps < 0.0_wp) .and. ((q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) + q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, &
874 & k, &
875 & l)) > ((rhoe - gs_min(lp)*ps_inf(lp)/(gs_min(lp) - 1))/qvs(lp)))) .or. ((ps >= 0.0_wp) .and. (ps < 1.0e-1_wp))) then
876 ! improve this initial condition
877 ps = 1.0e4_wp
878 end if
879
880 ! Loop until the solution for F(X) is satisfied Check whether I need to use both absolute and relative values for the
881 ! residual, and how to do it adequately. Dummy guess to start the pTg-equilibrium problem. improve this initial condition
882 r2d(1) = 0.0_wp; r2d(2) = 0.0_wp
883 deltamp(1) = 0.0_wp; deltamp(2) = 0.0_wp
884 do while (((sqrt(r2d(1)**2 + r2d(2)**2) > ptgalpha_eps) .and. ((sqrt(r2d(1)**2 + r2d(2)**2)/rhoe) > (ptgalpha_eps/1.e6_wp) &
885 & )) .or. (ns == 0))
886
887 ! Updating counter for the iterative procedure
888 ns = ns + 1
889
890 ! Auxiliary variables to help in the calculation of the residue
891 mcp = 0.0_wp; mcpd = 0.0_wp; mcvgp = 0.0_wp; mcvgp2 = 0.0_wp; mq = 0.0_wp; mqd = 0.0_wp
892 ! Those must be updated through the iterations, as they either depend on the partial masses for all fluids, or on the
893 ! equilibrium pressure
894
895# 402 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
896#if defined(MFC_OpenACC)
897# 402 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
898!$acc loop seq
899# 402 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
900#elif defined(MFC_OpenMP)
901# 402 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
902
903# 402 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
904#endif
905 do i = 1, num_fluids
906 ! sum of the total alpha*rho*cp of the system
907 mcp = mcp + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*cvs(i)*gs_min(i)
908
909 ! sum of the total alpha*rho*q of the system
910 mq = mq + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*qvs(i)
911
912 ! These auxiliary variables now need to be updated, as the partial densities now vary at every iteration
913 if ((i /= lp) .and. (i /= vp)) then
914 mcvgp = mcvgp + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*cvs(i)*(gs_min(i) - 1)/(ps + ps_inf(i))
915
916 mcvgp2 = mcvgp2 + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*cvs(i)*(gs_min(i) - 1)/((ps + ps_inf(i))**2)
917
918 mqd = mqd + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*qvs(i)
919
920 ! sum of the total alpha*rho*cp of the system
921 mcpd = mcpd + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*cvs(i)*gs_min(i)
922 end if
923 end do
924
925 ! calculating the (2D) Jacobian Matrix used in the solution of the pTg-quilibrium model
926
927 ! mass of the reacting liquid
928 ml = q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l)
929
930 ! mass of the two participating fluids
931 mt = q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) + q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l)
932
933 ts = 1/(mt*cvs(vp)*(gs_min(vp) - 1)/(ps + ps_inf(vp)) + ml*(cvs(lp)*(gs_min(lp) - 1)/(ps + ps_inf(lp)) - cvs(vp) &
934 & *(gs_min(vp) - 1)/(ps + ps_inf(vp))) + mcvgp)
935
936 dfdt = -(cvs(lp)*gs_min(lp) - cvs(vp)*gs_min(vp))*log(ts) - (qvps(lp) - qvps(vp)) + cvs(lp)*(gs_min(lp) - 1)*log(ps &
937 & + ps_inf(lp)) - cvs(vp)*(gs_min(vp) - 1)*log(ps + ps_inf(vp))
938
939 dtdm = -(cvs(lp)*(gs_min(lp) - 1)/(ps + ps_inf(lp)) - cvs(vp)*(gs_min(vp) - 1)/(ps + ps_inf(vp)))*ts**2
940
941 dtdp = (mt*cvs(vp)*(gs_min(vp) - 1)/(ps + ps_inf(vp))**2 + ml*(cvs(lp)*(gs_min(lp) - 1)/(ps + ps_inf(lp))**2 - cvs(vp) &
942 & *(gs_min(vp) - 1)/(ps + ps_inf(vp))**2) + mcvgp2)*ts**2
943
944 ! F = (F1,F2) is the function whose roots we are looking for x = (m1, p) are the independent variables. m1 = mass of the
945 ! first participant fluid, p = pressure F1 = 0 is the Gibbs free energy quality F2 = 0 is the enforcement of the
946 ! thermodynamic (total - kinectic) energy dF1dm
947 jac(1, 1) = dfdt*dtdm
948
949 ! dF1dp
950 jac(1, 2) = dfdt*dtdp + ts*(cvs(lp)*(gs_min(lp) - 1)/(ps + ps_inf(lp)) - cvs(vp)*(gs_min(vp) - 1)/(ps + ps_inf(vp)))
951
952 ! dF2dm
953 jac(2, &
954 & 1) = (qvs(vp) - qvs(lp) + (cvs(vp)*gs_min(vp) - cvs(lp)*gs_min(lp))/(ml*(cvs(lp)*(gs_min(lp) - 1)/(ps &
955 & + ps_inf(lp)) - cvs(vp)*(gs_min(vp) - 1)/(ps + ps_inf(vp))) + mt*cvs(vp)*(gs_min(vp) - 1)/(ps + ps_inf(vp)) &
956 & + mcvgp) - (ml*(cvs(vp)*gs_min(vp) - cvs(lp)*gs_min(lp)) - mt*cvs(vp)*gs_min(vp) - mcpd)*(cvs(lp)*(gs_min(lp) &
957 & - 1)/(ps + ps_inf(lp)) - cvs(vp)*(gs_min(vp) - 1)/(ps + ps_inf(vp)))/((ml*(cvs(lp)*(gs_min(lp) - 1)/(ps &
958 & + ps_inf(lp)) - cvs(vp)*(gs_min(vp) - 1)/(ps + ps_inf(vp))) + mt*cvs(vp)*(gs_min(vp) - 1)/(ps + ps_inf(vp)) &
959 & + mcvgp)**2))/1
960 ! dF2dp
961 jac(2, &
962 & 2) = (1 + (ml*(cvs(vp)*gs_min(vp) - cvs(lp)*gs_min(lp)) - mt*cvs(vp)*gs_min(vp) - mcpd)*(ml*(cvs(lp)*(gs_min(lp) &
963 & - 1)/(ps + ps_inf(lp))**2 - cvs(vp)*(gs_min(vp) - 1)/(ps + ps_inf(vp))**2) + mt*cvs(vp)*(gs_min(vp) - 1)/(ps &
964 & + ps_inf(vp))**2 + mcvgp2)/(ml*(cvs(lp)*(gs_min(lp) - 1)/(ps + ps_inf(lp)) - cvs(vp)*(gs_min(vp) - 1)/(ps &
965 & + ps_inf(vp))) + mt*cvs(vp)*(gs_min(vp) - 1)/(ps + ps_inf(vp)) + mcvgp)**2)/1
966
967 ! intermediate elements of J^{-1}
968 invjac(1, 1) = jac(2, 2)
969 invjac(1, 2) = -1.0_wp*jac(1, 2)
970 invjac(2, 1) = -1.0_wp*jac(2, 1)
971 invjac(2, 2) = jac(1, 1)
972
973 ! elements of J^{T}
974 tjac(1, 1) = jac(1, 1)
975 tjac(1, 2) = jac(2, 1)
976 tjac(2, 1) = jac(1, 2)
977 tjac(2, 2) = jac(2, 2)
978
979 ! dividing by det(J)
980 invjac = invjac/(jac(1, 1)*jac(2, 2) - jac(1, 2)*jac(2, 1))
981
982 ! calculating correction array for Newton's method
983 deltamp(1) = -1.0_wp*(invjac(1, 1)*r2d(1) + invjac(1, 2)*r2d(2))
984 deltamp(2) = -1.0_wp*(invjac(2, 1)*r2d(1) + invjac(2, 2)*r2d(2))
985
986 ! updating two reacting 'masses'. Recall that inert 'masses' do not change during the phase change liquid
987 q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) = q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) + om*deltamp(1)
988
989 ! gas
990 q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) = q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) - om*deltamp(1)
991
992 ! updating pressure
993 ps = ps + om*deltamp(2)
994
995 ! calculating residuals, which are (i) the difference between the Gibbs Free energy of the gas and the liquid and (ii)
996 ! the energy before and after the phase-change process.
997
998 ! mass of the reacting liquid
999 ml = q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l)
1000
1001 ! mass of the two participating fluids
1002 mt = q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) + q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l)
1003
1004 ts = 1/(mt*cvs(vp)*(gs_min(vp) - 1)/(ps + ps_inf(vp)) + ml*(cvs(lp)*(gs_min(lp) - 1)/(ps + ps_inf(lp)) - cvs(vp) &
1005 & *(gs_min(vp) - 1)/(ps + ps_inf(vp))) + mcvgp)
1006
1007 ! Gibbs Free Energy Equality condition (DG)
1008 r2d(1) = ts*((cvs(lp)*gs_min(lp) - cvs(vp)*gs_min(vp))*(1 - log(ts)) - (qvps(lp) - qvps(vp)) + cvs(lp)*(gs_min(lp) &
1009 & - 1)*log(ps + ps_inf(lp)) - cvs(vp)*(gs_min(vp) - 1)*log(ps + ps_inf(vp))) + qvs(lp) - qvs(vp)
1010
1011 ! Constant Energy Process condition (DE)
1012 r2d(2) = (rhoe + ps + ml*(qvs(vp) - qvs(lp)) - mt*qvs(vp) - mqd + (ml*(gs_min(vp)*cvs(vp) - gs_min(lp)*cvs(lp)) &
1013 & - mt*gs_min(vp)*cvs(vp) - mcpd)/(ml*(cvs(lp)*(gs_min(lp) - 1)/(ps + ps_inf(lp)) - cvs(vp)*(gs_min(vp) - 1)/(ps &
1014 & + ps_inf(vp))) + mt*cvs(vp)*(gs_min(vp) - 1)/(ps + ps_inf(vp)) + mcvgp))/1
1015 end do
1016
1017 ! common temperature
1018 ts = (rhoe + ps - mq)/mcp
1019
1020 end subroutine s_infinite_ptg_relaxation_k
1021
1022 !> Correct the partial densities of the reacting fluids in case one of them is negative but their sum is positive. Inert phases
1023 !! are not corrected at this moment
1024 subroutine s_correct_partial_densities(MCT, q_cons_vf, rM, j, k, l)
1025
1026
1027# 524 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1028#ifdef _CRAYFTN
1029# 524 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1030#if MFC_OpenACC
1031# 524 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1032!$acc routine seq
1033# 524 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1034#elif MFC_OpenMP
1035# 524 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1036
1037# 524 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1038
1039# 524 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1040!$omp declare target device_type(any)
1041# 524 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1042#else
1043# 524 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1044!DIR$ NOINLINE s_correct_partial_densities
1045# 524 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1046#endif
1047# 524 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1048#elif MFC_OpenACC
1049# 524 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1050!$acc routine seq
1051# 524 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1052#elif MFC_OpenMP
1053# 524 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1054
1055# 524 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1056
1057# 524 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1058!$omp declare target device_type(any)
1059# 524 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1060#endif
1061
1062 !> @name variables for the correction of the reacting partial densities
1063 !> @{
1064 real(wp), intent(out) :: mct
1065 type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
1066 real(wp), intent(inout) :: rm
1067 integer, intent(in) :: j, k, l
1068 !> @}
1069 if (rm < 0.0_wp) then
1070 if ((q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, &
1071 & l) >= -1.0_wp*mixm) .and. (q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) >= -1.0_wp*mixm)) then
1072 q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) = 0.0_wp
1073
1074 q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) = 0.0_wp
1075
1076 rm = q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) + q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l)
1077 end if
1078 end if
1079
1080 ! TODO: Consider partitioning partial densities instead of absolute-value correction
1081 mct = 2*mixm
1082
1083 ! correcting the partial densities of the reacting fluids. What to do for the nonreacting ones?
1084 if (q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) < 0.0_wp) then
1085 q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) = mct*rm
1086
1087 q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) = (1.0_wp - mct)*rm
1088 else if (q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) < 0.0_wp) then
1089 q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) = (1.0_wp - mct)*rm
1090
1091 q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) = mct*rm
1092 end if
1093
1094 end subroutine s_correct_partial_densities
1095
1096 !> Find the saturation temperature for a given saturation pressure using a Newton solver
1097 elemental subroutine s_tsat(pSat, TSat, TSIn)
1098
1099
1100# 563 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1101#ifdef _CRAYFTN
1102# 563 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1103#if MFC_OpenACC
1104# 563 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1105!$acc routine seq
1106# 563 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1107#elif MFC_OpenMP
1108# 563 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1109
1110# 563 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1111
1112# 563 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1113!$omp declare target device_type(any)
1114# 563 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1115#else
1116# 563 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1117!DIR$ NOINLINE s_TSat
1118# 563 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1119#endif
1120# 563 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1121#elif MFC_OpenACC
1122# 563 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1123!$acc routine seq
1124# 563 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1125#elif MFC_OpenMP
1126# 563 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1127
1128# 563 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1129
1130# 563 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1131!$omp declare target device_type(any)
1132# 563 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1133#endif
1134
1135 real(wp), intent(in) :: psat
1136 real(wp), intent(out) :: tsat
1137 real(wp), intent(in) :: tsin
1138 real(wp) :: dfdt, ft, om !< auxiliary variables
1139 ! Generic loop iterators
1140 integer :: ns
1141
1142 if ((f_approx_equal(psat, 0.0_wp)) .and. (f_approx_equal(tsin, 0.0_wp))) then
1143 ! assigning Saturation temperature
1144 tsat = 0.0_wp
1145 else
1146 ! calculating initial estimate for temperature in the TSat procedure. I will also use this variable to iterate over the
1147 ! Newton's solver
1148 tsat = tsin
1149
1150 ! iteration counter
1151 ns = 0
1152
1153 ! underrelaxation factor
1154 om = 1.0e-3_wp
1155
1156 ! FT must be initialized before the do while condition is evaluated. Fortran .or. is not short-circuit: abs(FT) is
1157 ! always evaluated even when ns == 0, so FT must have a defined value here.
1158 ft = huge(1.0_wp)
1159
1160 do while ((abs(ft) > ptgalpha_eps) .or. (ns == 0))
1161 ! increasing counter
1162 ns = ns + 1
1163
1164 ! calculating residual
1165 ft = tsat*((cvs(lp)*gs_min(lp) - cvs(vp)*gs_min(vp))*(1 - log(tsat)) - (qvps(lp) - qvps(vp)) + cvs(lp)*(gs_min(lp) &
1166 & - 1)*log(psat + ps_inf(lp)) - cvs(vp)*(gs_min(vp) - 1)*log(psat + ps_inf(vp))) + qvs(lp) - qvs(vp)
1167
1168 ! calculating the jacobian
1169 dfdt = -(cvs(lp)*gs_min(lp) - cvs(vp)*gs_min(vp))*log(tsat) - (qvps(lp) - qvps(vp)) + cvs(lp)*(gs_min(lp) - 1) &
1170 & *log(psat + ps_inf(lp)) - cvs(vp)*(gs_min(vp) - 1)*log(psat + ps_inf(vp))
1171
1172 ! updating saturation temperature
1173 tsat = tsat - om*ft/dfdt
1174
1175 if (abs(ft) <= ptgalpha_eps) exit
1176 end do
1177 end if
1178
1179 end subroutine s_tsat
1180
1181 !> Finalize the phase change module
1185#endif
1186end module m_phase_change
type(scalar_field), dimension(sys_size), intent(inout) q_cons_vf
real(wp), intent(inout) rm
integer, intent(in) k
real(wp), intent(out) mct
integer, intent(in) j
integer, intent(in) l
Shared derived types for field data, patch geometry, bubble dynamics, and MPI I/O structures.
Global parameters for the computational domain, fluid properties, and simulation algorithm configurat...
integer num_fluids
number of fluids in the simulation
real(wp) ptgalpha_eps
trigger parameter for the pTg relaxation procedure, phase change model
integer relax_model
Relaxation model.
integer model_eqns
Multicomponent flow model.
real(wp), dimension(:), allocatable ps_inf
real(wp), dimension(:), allocatable cvs
real(wp), dimension(:), allocatable qvps
real(wp), dimension(:), allocatable qvs
real(wp) palpha_eps
trigger parameter for the p relaxation procedure, phase change model
real(wp), dimension(:), allocatable gs_min
type(eqn_idx_info) eqn_idx
All conserved-variable equation index ranges and scalars.
Basic floating-point utilities: approximate equality, default detection, and coordinate bounds.
logical elemental function, public f_approx_equal(a, b, tol_input)
Check if two floating point numbers of wp are within tolerance.
MPI halo exchange, domain decomposition, and buffer packing/unpacking for the simulation solver.
Phase transition relaxation solvers for liquid-vapor flows with cavitation and boiling.
impure subroutine, public s_finalize_relaxation_solver_module
Finalize the phase change module.
subroutine, public s_infinite_relaxation_k(q_cons_vf)
Apply pT- or pTg-equilibrium relaxation with mass depletion based on the incoming state conditions.
elemental subroutine s_tsat(psat, tsat, tsin)
Find the saturation temperature for a given saturation pressure using a Newton solver.
integer, parameter vp
index for the vapor phase of the reacting fluid
integer, parameter lp
index for the liquid phase of the reacting fluid
subroutine s_infinite_pt_relaxation_k(j, k, l, mfl, ps, p_infpt, q_cons_vf, rhoe, ts)
Apply pT-equilibrium relaxation for N fluids.
impure subroutine, public s_relaxation_solver(q_cons_vf)
Dispatch to the correct relaxation solver. Replaces the procedure pointer, which CCE is breaking on.
real(wp), parameter tcr
Critical temperature of water [K].
integer, parameter max_iter
max # of iterations
impure subroutine, public s_initialize_phasechange_module
Initialize the phase change module by setting saturation curve coefficients for pT- or pTg-equilibriu...
real(wp), parameter mixm
Mixture mass fraction threshold for triggering phase change.
subroutine s_infinite_ptg_relaxation_k(j, k, l, ps, p_infpt, rhoe, q_cons_vf, ts)
Apply pTg-equilibrium relaxation for N fluids under pT and 2 fluids under pTg-equilibrium....
real(wp), parameter pcr
Critical pressure of water [Pa].
Conservative-to-primitive variable conversion, mixture property evaluation, and pressure computation.
Derived type annexing a scalar field (SF).