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 + contxb - 1)%sf(j, k, l)
469
470 ! Total Volume Fraction
471 tvf = tvf + q_cons_vf(i + advxb - 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 + contxb - 1)%sf(j, k, l) + q_cons_vf(vp + contxb - 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 + contxb - 1)%sf(j, k, l)
484
485 m2 = q_cons_vf(vp + contxb - 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 = momxb, momxe
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(e_idx)%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 + contxb - 1)%sf(j, k, &
515 & l) > mixm*rm) .and. (q_cons_vf(vp + contxb - 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 + contxb - 1)%sf(j, k, l) = mixm*rm
522
523 ! transferring the total mass to vapor
524 q_cons_vf(vp + contxb - 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 + contxb - 1)%sf(j, k, l) = (1.0_wp - mixm)*rm
534
535 ! depleting the mass of vapor
536 q_cons_vf(vp + contxb - 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 + contxb - 1)%sf(j, k, l) = mixm*rm
554
555 ! correcting the vapor partial density
556 q_cons_vf(vp + contxb - 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 + contxb - 1)%sf(j, k, l) = (1.0_wp - mixm)*rm
566
567 ! correcting the vapor partial density
568 q_cons_vf(vp + contxb - 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 + contxb - 1)%sf(j, k, l) = m1
572
573 ! vapor
574 q_cons_vf(vp + contxb - 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 + advxb - 1)%sf(j, k, l) = q_cons_vf(i + contxb - 1)%sf(j, k, l)/rhok(i)
627
628 ! alpha*rho*e
629 if (model_eqns == 3) then
630 q_cons_vf(i + intxb - 1)%sf(j, k, l) = q_cons_vf(i + contxb - 1)%sf(j, k, l)*ek(i)
631 end if
632
633 ! Total entropy
634 rhos = rhos + q_cons_vf(i + contxb - 1)%sf(j, k, l)*sk(i)
635 end do
636 end do
637 end do
638 end do
639
640# 248 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
641#if defined(MFC_OpenACC)
642# 248 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
643!$acc end parallel loop
644# 248 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
645#elif defined(MFC_OpenMP)
646# 248 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
647
648# 248 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
649!$omp end target teams loop
650# 248 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
651#endif
652
653 end subroutine s_infinite_relaxation_k
654
655 !> Apply pT-equilibrium relaxation for N fluids
656 !! @param MFL flag: 0=gas, 1=liquid, 2=mixture
657 subroutine s_infinite_pt_relaxation_k(j, k, l, MFL, pS, p_infpT, q_cons_vf, rhoe, TS)
658
659
660# 256 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
661#ifdef _CRAYFTN
662# 256 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
663#if MFC_OpenACC
664# 256 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
665!$acc routine seq
666# 256 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
667#elif MFC_OpenMP
668# 256 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
669
670# 256 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
671
672# 256 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
673!$omp declare target device_type(any)
674# 256 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
675#else
676# 256 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
677!DIR$ NOINLINE s_infinite_pt_relaxation_k
678# 256 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
679#endif
680# 256 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
681#elif MFC_OpenACC
682# 256 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
683!$acc routine seq
684# 256 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
685#elif MFC_OpenMP
686# 256 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
687
688# 256 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
689
690# 256 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
691!$omp declare target device_type(any)
692# 256 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
693#endif
694
695 ! initializing variables
696 integer, intent(in) :: j, k, l, MFL
697 real(wp), intent(out) :: pS
698 real(wp), dimension(1:), intent(out) :: p_infpT
699 type(scalar_field), dimension(sys_size), intent(in) :: q_cons_vf
700 real(wp), intent(in) :: rhoe
701 real(wp), intent(out) :: TS
702 real(wp) :: gp, gpp, hp, pO, mCP, mQ !< variables for the Newton Solver
703 real(wp) :: p_infpT_sum
704 integer :: i, ns !< generic loop iterators
705 ! auxiliary variables for the pT-equilibrium solver
706 mcp = 0.0_wp; mq = 0.0_wp; p_infpt_sum = 0._wp
707
708# 270 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
709#if defined(MFC_OpenACC)
710# 270 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
711!$acc loop seq
712# 270 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
713#elif defined(MFC_OpenMP)
714# 270 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
715
716# 270 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
717#endif
718 do i = 1, num_fluids
719 p_infpt(i) = ps_inf(i)
720 p_infpt_sum = p_infpt_sum + abs(p_infpt(i))
721 end do
722 ! Performing tests before initializing the pT-equilibrium
723
724# 276 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
725#if defined(MFC_OpenACC)
726# 276 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
727!$acc loop seq
728# 276 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
729#elif defined(MFC_OpenMP)
730# 276 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
731
732# 276 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
733#endif
734 do i = 1, num_fluids
735 ! sum of the total alpha*rho*cp of the system
736 mcp = mcp + q_cons_vf(i + contxb - 1)%sf(j, k, l)*cvs(i)*gs_min(i)
737
738 ! sum of the total alpha*rho*q of the system
739 mq = mq + q_cons_vf(i + contxb - 1)%sf(j, k, l)*qvs(i)
740 end do
741
742# 293 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
743
744 ! Checking energy constraint
745 if ((rhoe - mq - minval(p_infpt)) < 0.0_wp) then
746 if ((mfl == 0) .or. (mfl == 1)) then
747 ! Assigning zero values for mass depletion cases pressure
748 ps = 0.0_wp
749
750 ! temperature
751 ts = 0.0_wp
752
753 return
754 end if
755 end if
756
757 ! calculating initial estimate for pressure in the pT-relaxation procedure. I will also use this variable to iterate over
758 ! the Newton's solver
759 po = 0.0_wp
760
761 ! Maybe improve this condition afterwards. As long as the initial guess is in between -min(ps_inf) and infinity, a solution
762 ! should be able to be found.
763 ps = 1.0e4_wp
764
765 ! Newton Solver for the pT-equilibrium
766 ns = 0
767 ! change this relative error metric. 1.e4_wp is just arbitrary
768 do while ((abs(ps - po) > palpha_eps) .and. (abs((ps - po)/po) > palpha_eps/1.e4_wp) .or. (ns == 0))
769 ! increasing counter
770 ns = ns + 1
771
772 ! updating old pressure
773 po = ps
774
775 ! updating functions used in the Newton's solver
776 gpp = 0.0_wp; gp = 0.0_wp; hp = 0.0_wp
777
778# 327 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
779#if defined(MFC_OpenACC)
780# 327 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
781!$acc loop seq
782# 327 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
783#elif defined(MFC_OpenMP)
784# 327 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
785
786# 327 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
787#endif
788 do i = 1, num_fluids
789 gp = gp + (gs_min(i) - 1.0_wp)*q_cons_vf(i + contxb - 1)%sf(j, k, l)*cvs(i)*(rhoe + ps - mq)/(mcp*(ps + p_infpt(i)))
790
791 gpp = gpp + (gs_min(i) - 1.0_wp)*q_cons_vf(i + contxb - 1)%sf(j, k, &
792 & l)*cvs(i)*(p_infpt(i) - rhoe + mq)/(mcp*(ps + p_infpt(i))**2)
793 end do
794
795 hp = 1.0_wp/(rhoe + ps - mq) + 1.0_wp/(ps + minval(p_infpt))
796
797 ! updating common pressure for the newton solver
798 ps = po + ((1.0_wp - gp)/gpp)/(1.0_wp - (1.0_wp - gp + abs(1.0_wp - gp))/(2.0_wp*gpp)*hp)
799 end do
800
801 ! common temperature
802 ts = (rhoe + ps - mq)/mcp
803
804 end subroutine s_infinite_pt_relaxation_k
805
806 !> Apply pTg-equilibrium relaxation for N fluids under pT and 2 fluids under pTg-equilibrium. There is a final common p and T
807 !! during relaxation
808 subroutine s_infinite_ptg_relaxation_k(j, k, l, pS, p_infpT, rhoe, q_cons_vf, TS)
809
810
811# 350 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
812#ifdef _CRAYFTN
813# 350 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
814#if MFC_OpenACC
815# 350 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
816!$acc routine seq
817# 350 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
818#elif MFC_OpenMP
819# 350 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
820
821# 350 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
822
823# 350 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
824!$omp declare target device_type(any)
825# 350 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
826#else
827# 350 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
828!DIR$ NOINLINE s_infinite_ptg_relaxation_k
829# 350 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
830#endif
831# 350 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
832#elif MFC_OpenACC
833# 350 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
834!$acc routine seq
835# 350 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
836#elif MFC_OpenMP
837# 350 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
838
839# 350 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
840
841# 350 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
842!$omp declare target device_type(any)
843# 350 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
844#endif
845
846 integer, intent(in) :: j, k, l
847 real(wp), intent(inout) :: pS
848 real(wp), dimension(1:), intent(in) :: p_infpT
849 real(wp), intent(in) :: rhoe
850 type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
851 real(wp), intent(inout) :: TS
852# 361 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
853 real(wp), dimension(num_fluids) :: p_infpTg !< stiffness for the participating fluids for pTg-equilibrium
854# 363 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
855 real(wp), dimension(2, 2) :: Jac, InvJac, TJac !< matrices for the Newton Solver
856 real(wp), dimension(2) :: R2D, DeltamP !< residual and correction array
857 real(wp) :: Om !< underrelaxation factor
858 real(wp) :: mCP, mCPD, mCVGP, mCVGP2, mQ, mQD !< auxiliary variables for the pTg-solver
859 real(wp) :: ml, mT, dFdT, dTdm, dTdp
860
861 !> Generic loop iterators
862 integer :: i, ns
863 ! pTg-equilibrium solution procedure Newton Solver parameters counter
864 ns = 0
865
866 ! Relaxation factor
867 om = 1.0e-3_wp
868
869 p_infptg = p_infpt
870
871 if (((ps < 0.0_wp) .and. ((q_cons_vf(lp + contxb - 1)%sf(j, k, l) + q_cons_vf(vp + contxb - 1)%sf(j, k, &
872 & 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
873 ! improve this initial condition
874 ps = 1.0e4_wp
875 end if
876
877 ! Loop until the solution for F(X) is satisfied Check whether I need to use both absolute and relative values for the
878 ! residual, and how to do it adequately. Dummy guess to start the pTg-equilibrium problem. improve this initial condition
879 r2d(1) = 0.0_wp; r2d(2) = 0.0_wp
880 deltamp(1) = 0.0_wp; deltamp(2) = 0.0_wp
881 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) &
882 & )) .or. (ns == 0))
883
884 ! Updating counter for the iterative procedure
885 ns = ns + 1
886
887 ! Auxiliary variables to help in the calculation of the residue
888 mcp = 0.0_wp; mcpd = 0.0_wp; mcvgp = 0.0_wp; mcvgp2 = 0.0_wp; mq = 0.0_wp; mqd = 0.0_wp
889 ! Those must be updated through the iterations, as they either depend on the partial masses for all fluids, or on the
890 ! equilibrium pressure
891
892# 399 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
893#if defined(MFC_OpenACC)
894# 399 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
895!$acc loop seq
896# 399 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
897#elif defined(MFC_OpenMP)
898# 399 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
899
900# 399 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
901#endif
902 do i = 1, num_fluids
903 ! sum of the total alpha*rho*cp of the system
904 mcp = mcp + q_cons_vf(i + contxb - 1)%sf(j, k, l)*cvs(i)*gs_min(i)
905
906 ! sum of the total alpha*rho*q of the system
907 mq = mq + q_cons_vf(i + contxb - 1)%sf(j, k, l)*qvs(i)
908
909 ! These auxiliary variables now need to be updated, as the partial densities now vary at every iteration
910 if ((i /= lp) .and. (i /= vp)) then
911 mcvgp = mcvgp + q_cons_vf(i + contxb - 1)%sf(j, k, l)*cvs(i)*(gs_min(i) - 1)/(ps + ps_inf(i))
912
913 mcvgp2 = mcvgp2 + q_cons_vf(i + contxb - 1)%sf(j, k, l)*cvs(i)*(gs_min(i) - 1)/((ps + ps_inf(i))**2)
914
915 mqd = mqd + q_cons_vf(i + contxb - 1)%sf(j, k, l)*qvs(i)
916
917 ! sum of the total alpha*rho*cp of the system
918 mcpd = mcpd + q_cons_vf(i + contxb - 1)%sf(j, k, l)*cvs(i)*gs_min(i)
919 end if
920 end do
921
922 ! calculating the (2D) Jacobian Matrix used in the solution of the pTg-quilibrium model
923
924 ! mass of the reacting liquid
925 ml = q_cons_vf(lp + contxb - 1)%sf(j, k, l)
926
927 ! mass of the two participating fluids
928 mt = q_cons_vf(lp + contxb - 1)%sf(j, k, l) + q_cons_vf(vp + contxb - 1)%sf(j, k, l)
929
930 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) &
931 & *(gs_min(vp) - 1)/(ps + ps_inf(vp))) + mcvgp)
932
933 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 &
934 & + ps_inf(lp)) - cvs(vp)*(gs_min(vp) - 1)*log(ps + ps_inf(vp))
935
936 dtdm = -(cvs(lp)*(gs_min(lp) - 1)/(ps + ps_inf(lp)) - cvs(vp)*(gs_min(vp) - 1)/(ps + ps_inf(vp)))*ts**2
937
938 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) &
939 & *(gs_min(vp) - 1)/(ps + ps_inf(vp))**2) + mcvgp2)*ts**2
940
941 ! F = (F1,F2) is the function whose roots we are looking for x = (m1, p) are the independent variables. m1 = mass of the
942 ! first participant fluid, p = pressure F1 = 0 is the Gibbs free energy quality F2 = 0 is the enforcement of the
943 ! thermodynamic (total - kinectic) energy dF1dm
944 jac(1, 1) = dfdt*dtdm
945
946 ! dF1dp
947 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)))
948
949 ! dF2dm
950 jac(2, &
951 & 1) = (qvs(vp) - qvs(lp) + (cvs(vp)*gs_min(vp) - cvs(lp)*gs_min(lp))/(ml*(cvs(lp)*(gs_min(lp) - 1)/(ps &
952 & + ps_inf(lp)) - cvs(vp)*(gs_min(vp) - 1)/(ps + ps_inf(vp))) + mt*cvs(vp)*(gs_min(vp) - 1)/(ps + ps_inf(vp)) &
953 & + mcvgp) - (ml*(cvs(vp)*gs_min(vp) - cvs(lp)*gs_min(lp)) - mt*cvs(vp)*gs_min(vp) - mcpd)*(cvs(lp)*(gs_min(lp) &
954 & - 1)/(ps + ps_inf(lp)) - cvs(vp)*(gs_min(vp) - 1)/(ps + ps_inf(vp)))/((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)**2))/1
957 ! dF2dp
958 jac(2, &
959 & 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) &
960 & - 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 &
961 & + ps_inf(vp))**2 + mcvgp2)/(ml*(cvs(lp)*(gs_min(lp) - 1)/(ps + ps_inf(lp)) - cvs(vp)*(gs_min(vp) - 1)/(ps &
962 & + ps_inf(vp))) + mt*cvs(vp)*(gs_min(vp) - 1)/(ps + ps_inf(vp)) + mcvgp)**2)/1
963
964 ! intermediate elements of J^{-1}
965 invjac(1, 1) = jac(2, 2)
966 invjac(1, 2) = -1.0_wp*jac(1, 2)
967 invjac(2, 1) = -1.0_wp*jac(2, 1)
968 invjac(2, 2) = jac(1, 1)
969
970 ! elements of J^{T}
971 tjac(1, 1) = jac(1, 1)
972 tjac(1, 2) = jac(2, 1)
973 tjac(2, 1) = jac(1, 2)
974 tjac(2, 2) = jac(2, 2)
975
976 ! dividing by det(J)
977 invjac = invjac/(jac(1, 1)*jac(2, 2) - jac(1, 2)*jac(2, 1))
978
979 ! calculating correction array for Newton's method
980 deltamp(1) = -1.0_wp*(invjac(1, 1)*r2d(1) + invjac(1, 2)*r2d(2))
981 deltamp(2) = -1.0_wp*(invjac(2, 1)*r2d(1) + invjac(2, 2)*r2d(2))
982
983 ! updating two reacting 'masses'. Recall that inert 'masses' do not change during the phase change liquid
984 q_cons_vf(lp + contxb - 1)%sf(j, k, l) = q_cons_vf(lp + contxb - 1)%sf(j, k, l) + om*deltamp(1)
985
986 ! gas
987 q_cons_vf(vp + contxb - 1)%sf(j, k, l) = q_cons_vf(vp + contxb - 1)%sf(j, k, l) - om*deltamp(1)
988
989 ! updating pressure
990 ps = ps + om*deltamp(2)
991
992 ! calculating residuals, which are (i) the difference between the Gibbs Free energy of the gas and the liquid and (ii)
993 ! the energy before and after the phase-change process.
994
995 ! mass of the reacting liquid
996 ml = q_cons_vf(lp + contxb - 1)%sf(j, k, l)
997
998 ! mass of the two participating fluids
999 mt = q_cons_vf(lp + contxb - 1)%sf(j, k, l) + q_cons_vf(vp + contxb - 1)%sf(j, k, l)
1000
1001 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) &
1002 & *(gs_min(vp) - 1)/(ps + ps_inf(vp))) + mcvgp)
1003
1004 ! Gibbs Free Energy Equality condition (DG)
1005 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) &
1006 & - 1)*log(ps + ps_inf(lp)) - cvs(vp)*(gs_min(vp) - 1)*log(ps + ps_inf(vp))) + qvs(lp) - qvs(vp)
1007
1008 ! Constant Energy Process condition (DE)
1009 r2d(2) = (rhoe + ps + ml*(qvs(vp) - qvs(lp)) - mt*qvs(vp) - mqd + (ml*(gs_min(vp)*cvs(vp) - gs_min(lp)*cvs(lp)) &
1010 & - 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 &
1011 & + ps_inf(vp))) + mt*cvs(vp)*(gs_min(vp) - 1)/(ps + ps_inf(vp)) + mcvgp))/1
1012 end do
1013
1014 ! common temperature
1015 ts = (rhoe + ps - mq)/mcp
1016
1017 end subroutine s_infinite_ptg_relaxation_k
1018
1019 !> Correct the partial densities of the reacting fluids in case one of them is negative but their sum is positive. Inert phases
1020 !! are not corrected at this moment
1021 subroutine s_correct_partial_densities(MCT, q_cons_vf, rM, j, k, l)
1022
1023
1024# 521 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1025#ifdef _CRAYFTN
1026# 521 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1027#if MFC_OpenACC
1028# 521 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1029!$acc routine seq
1030# 521 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1031#elif MFC_OpenMP
1032# 521 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1033
1034# 521 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1035
1036# 521 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1037!$omp declare target device_type(any)
1038# 521 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1039#else
1040# 521 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1041!DIR$ NOINLINE s_correct_partial_densities
1042# 521 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1043#endif
1044# 521 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1045#elif MFC_OpenACC
1046# 521 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1047!$acc routine seq
1048# 521 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1049#elif MFC_OpenMP
1050# 521 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1051
1052# 521 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1053
1054# 521 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1055!$omp declare target device_type(any)
1056# 521 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1057#endif
1058
1059 !> @name variables for the correction of the reacting partial densities
1060 !> @{
1061 real(wp), intent(out) :: mct
1062 type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
1063 real(wp), intent(inout) :: rm
1064 integer, intent(in) :: j, k, l
1065 !> @}
1066 if (rm < 0.0_wp) then
1067 if ((q_cons_vf(lp + contxb - 1)%sf(j, k, l) >= -1.0_wp*mixm) .and. (q_cons_vf(vp + contxb - 1)%sf(j, k, &
1068 & l) >= -1.0_wp*mixm)) then
1069 q_cons_vf(lp + contxb - 1)%sf(j, k, l) = 0.0_wp
1070
1071 q_cons_vf(vp + contxb - 1)%sf(j, k, l) = 0.0_wp
1072
1073 rm = q_cons_vf(lp + contxb - 1)%sf(j, k, l) + q_cons_vf(vp + contxb - 1)%sf(j, k, l)
1074 end if
1075 end if
1076
1077 ! TODO: Consider partitioning partial densities instead of absolute-value correction
1078 mct = 2*mixm
1079
1080 ! correcting the partial densities of the reacting fluids. What to do for the nonreacting ones?
1081 if (q_cons_vf(lp + contxb - 1)%sf(j, k, l) < 0.0_wp) then
1082 q_cons_vf(lp + contxb - 1)%sf(j, k, l) = mct*rm
1083
1084 q_cons_vf(vp + contxb - 1)%sf(j, k, l) = (1.0_wp - mct)*rm
1085 else if (q_cons_vf(vp + contxb - 1)%sf(j, k, l) < 0.0_wp) then
1086 q_cons_vf(lp + contxb - 1)%sf(j, k, l) = (1.0_wp - mct)*rm
1087
1088 q_cons_vf(vp + contxb - 1)%sf(j, k, l) = mct*rm
1089 end if
1090
1091 end subroutine s_correct_partial_densities
1092
1093 !> Find the saturation temperature for a given saturation pressure using a Newton solver
1094 elemental subroutine s_tsat(pSat, TSat, TSIn)
1095
1096
1097# 560 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1098#ifdef _CRAYFTN
1099# 560 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1100#if MFC_OpenACC
1101# 560 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1102!$acc routine seq
1103# 560 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1104#elif MFC_OpenMP
1105# 560 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1106
1107# 560 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1108
1109# 560 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1110!$omp declare target device_type(any)
1111# 560 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1112#else
1113# 560 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1114!DIR$ NOINLINE s_TSat
1115# 560 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1116#endif
1117# 560 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1118#elif MFC_OpenACC
1119# 560 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1120!$acc routine seq
1121# 560 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1122#elif MFC_OpenMP
1123# 560 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1124
1125# 560 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1126
1127# 560 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1128!$omp declare target device_type(any)
1129# 560 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1130#endif
1131
1132 real(wp), intent(in) :: psat
1133 real(wp), intent(out) :: tsat
1134 real(wp), intent(in) :: tsin
1135 real(wp) :: dfdt, ft, om !< auxiliary variables
1136 ! Generic loop iterators
1137 integer :: ns
1138
1139 if ((f_approx_equal(psat, 0.0_wp)) .and. (f_approx_equal(tsin, 0.0_wp))) then
1140 ! assigning Saturation temperature
1141 tsat = 0.0_wp
1142 else
1143 ! calculating initial estimate for temperature in the TSat procedure. I will also use this variable to iterate over the
1144 ! Newton's solver
1145 tsat = tsin
1146
1147 ! iteration counter
1148 ns = 0
1149
1150 ! underrelaxation factor
1151 om = 1.0e-3_wp
1152
1153 ! FT must be initialized before the do while condition is evaluated. Fortran .or. is not short-circuit: abs(FT) is
1154 ! always evaluated even when ns == 0, so FT must have a defined value here.
1155 ft = huge(1.0_wp)
1156
1157 do while ((abs(ft) > ptgalpha_eps) .or. (ns == 0))
1158 ! increasing counter
1159 ns = ns + 1
1160
1161 ! calculating residual
1162 ft = tsat*((cvs(lp)*gs_min(lp) - cvs(vp)*gs_min(vp))*(1 - log(tsat)) - (qvps(lp) - qvps(vp)) + cvs(lp)*(gs_min(lp) &
1163 & - 1)*log(psat + ps_inf(lp)) - cvs(vp)*(gs_min(vp) - 1)*log(psat + ps_inf(vp))) + qvs(lp) - qvs(vp)
1164
1165 ! calculating the jacobian
1166 dfdt = -(cvs(lp)*gs_min(lp) - cvs(vp)*gs_min(vp))*log(tsat) - (qvps(lp) - qvps(vp)) + cvs(lp)*(gs_min(lp) - 1) &
1167 & *log(psat + ps_inf(lp)) - cvs(vp)*(gs_min(vp) - 1)*log(psat + ps_inf(vp))
1168
1169 ! updating saturation temperature
1170 tsat = tsat - om*ft/dfdt
1171
1172 if (abs(ft) <= ptgalpha_eps) exit
1173 end do
1174 end if
1175
1176 end subroutine s_tsat
1177
1178 !> Finalize the phase change module
1182#endif
1183end 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.
Defines global parameters for the computational domain, simulation algorithm, and initial conditions.
integer num_fluids
Number of different fluids present in the flow.
real(wp) ptgalpha_eps
trigger parameter for the pTg relaxation procedure, phase change model
integer relax_model
Relax Model.
integer model_eqns
Multicomponent flow model.
integer e_idx
Index of total energy equation.
real(wp) palpha_eps
trigger parameter for the p relaxation procedure, phase change model
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.
Broadcasts user inputs and decomposes the domain across MPI ranks for pre-processing.
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.
real(wp), dimension(:), allocatable, public ps_inf
real(wp), dimension(:), allocatable, public gs_min
real(wp), dimension(:), allocatable, public qvs
real(wp), dimension(:), allocatable, public cvs
real(wp), dimension(:), allocatable, public qvps
Derived type annexing a scalar field (SF).