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