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# 15 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
31# 16 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
32# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
33
34# 24 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
35
36# 53 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
37
38# 65 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
39
40# 75 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
41
42# 105 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
43
44# 117 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
45
46# 127 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
47
48# 174 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
49! New line at end of file is required for FYPP
50# 2 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
51# 1 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 1
52# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
53# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
54# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
55# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
56# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
57# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
58
59# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
60# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
61# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
62
63# 15 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
64# 16 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
65# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
66
67# 24 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
68
69# 53 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
70
71# 65 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
72
73# 75 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
74
75# 105 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
76
77# 117 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
78
79# 127 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
80
81# 174 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
82! New line at end of file is required for FYPP
83# 2 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 2
84
85# 4 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
86# 5 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
87# 6 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
88# 7 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
89# 8 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
90
91# 20 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
92
93# 43 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
94
95# 48 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
96
97# 53 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
98
99# 58 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
100
101# 63 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
102
103# 68 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
104
105# 81 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
106
107# 86 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
108
109# 91 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
110
111# 96 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
112
113# 101 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
114
115# 106 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
116
117# 111 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
118
119# 116 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
120
121# 121 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
122
123# 126 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
124
125# 156 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
126
127# 197 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
128
129# 211 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
130
131# 236 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
132
133# 247 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
134
135# 249 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
136# 260 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
137
138# 310 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
139
140# 320 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
141
142# 330 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
143
144# 339 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
145
146# 356 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
147
148# 366 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
149
150# 373 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
151
152# 379 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
153
154# 385 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
155
156# 391 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
157
158# 397 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
159
160# 403 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
161! New line at end of file is required for FYPP
162# 3 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
163# 1 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 1
164# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
165# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
166# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
167# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
168# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
169# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
170
171# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
172# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
173# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
174
175# 15 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
176# 16 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
177# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
178
179# 24 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
180
181# 53 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
182
183# 65 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
184
185# 75 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
186
187# 105 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
188
189# 117 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
190
191# 127 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
192
193# 174 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
194! New line at end of file is required for FYPP
195# 2 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 2
196
197# 7 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
198
199# 17 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
200
201# 22 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
202
203# 27 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
204
205# 32 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
206
207# 37 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
208
209# 42 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
210
211# 47 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
212
213# 52 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
214
215# 57 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
216
217# 62 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
218
219# 73 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
220
221# 78 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
222
223# 83 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
224
225# 88 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
226
227# 103 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
228
229# 131 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
230
231# 160 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
232
233# 175 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
234
235# 193 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
236
237# 215 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
238
239# 244 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
240
241# 259 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
242
243# 269 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
244
245# 278 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
246
247# 294 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
248
249# 304 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
250
251# 311 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
252! New line at end of file is required for FYPP
253# 4 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
254
255! GPU parallel region (scalar reductions, maxval/minval)
256# 23 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
257
258! GPU parallel loop over threads (most common GPU macro)
259# 43 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
260
261! Required closing for GPU_PARALLEL_LOOP
262# 55 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
263
264! Mark routine for device compilation
265# 112 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
266
267! Declare device-resident data
268# 130 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
269
270! Inner loop within a GPU parallel region
271# 145 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
272
273! Scoped GPU data region
274# 164 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
275
276! Host code with device pointers (for MPI with GPU buffers)
277# 193 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
278
279! Allocate device memory (unscoped)
280# 207 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
281
282! Free device memory
283# 219 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
284
285! Atomic operation on device
286# 231 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
287
288! End atomic capture block
289# 242 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
290
291! Copy data between host and device
292# 254 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
293
294! Synchronization barrier
295# 266 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
296
297! Import GPU library module (openacc or omp_lib)
298# 275 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
299
300! Emit code only for AMD compiler
301# 282 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
302
303! Emit code for non-Cray compilers
304# 289 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
305
306! Emit code only for Cray compiler
307# 296 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
308
309! Emit code for non-NVIDIA compilers
310# 303 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
311
312# 305 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
313# 306 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
314! New line at end of file is required for FYPP
315# 2 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 2
316
317# 14 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
318
319! Caution: This macro requires the use of a binding script to set CUDA_VISIBLE_DEVICES, such that we have one GPU device per MPI
320! rank. That's because for both cudaMemAdvise (preferred location) and cudaMemPrefetchAsync we use location = device_id = 0. For an
321! example see misc/nvidia_uvm/bind.sh.
322# 52 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
323
324! Allocate and create GPU device memory
325# 72 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
326
327! Free GPU device memory and deallocate
328# 80 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
329
330! Cray-specific GPU pointer setup for vector fields
331# 104 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
332
333! Cray-specific GPU pointer setup for scalar fields
334# 120 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
335
336! Cray-specific GPU pointer setup for acoustic source spatials
337# 145 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
338
339# 151 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
340
341# 158 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
342! New line at end of file is required for FYPP
343# 7 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp" 2
344
345!> @brief Phase transition relaxation solvers for liquid-vapor flows with cavitation and boiling
346module m_phase_change
347
348 use m_derived_types
349 use m_global_parameters
350 use m_mpi_proxy
351 use m_variables_conversion
352 use m_eos
353 use ieee_arithmetic
354 use m_helper_basic
355 use m_constants, only: model_eqns_6eq
356
357 implicit none
358
359 private
360 public :: s_initialize_phasechange_module, s_relaxation_solver, s_infinite_relaxation_k, s_finalize_relaxation_solver_module
361
362 !> @name Parameters for the first order transition phase change
363 !> @{
364 integer, parameter :: max_iter = 100000 !< max Newton iterations before accepting the last iterate
365 real(wp), parameter :: pcr = 4.94e7_wp !< Critical pressure of water [Pa]
366 real(wp), parameter :: tcr = 385.05_wp + 273.15_wp !< Critical temperature of water [K]
367 integer, parameter :: ptg_ls_max = 30 !< max backtracking-line-search halvings in the pTg solver
368 real(wp), parameter :: mixm = 1.0e-8_wp !< Mixture mass fraction threshold for triggering phase change
369 integer, parameter :: lp = 1 !< index for the liquid phase of the reacting fluid
370 integer, parameter :: vp = 2 !< index for the vapor phase of the reacting fluid
371 !> @}
372
373contains
374
375 !> Dispatch to the correct relaxation solver. Replaces the procedure pointer, which CCE is breaking on.
376 impure subroutine s_relaxation_solver(q_cons_vf)
377
378 type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
379 ! This is empty because in current master the procedure pointer was never assigned
380
381 if (.not. (.false.)) then
382# 44 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
383 call s_mpi_abort("m_phase_change.fpp:44: " // "Assertion failed: .false.. " &
384# 44 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
385 & // "s_relaxation_solver called but it currently does nothing")
386# 44 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
387 end if
388
389 end subroutine s_relaxation_solver
390
391 !> Initialize the phase change module (no module-level state to set up; the pT/pTg relaxation solvers are self-contained)
392 impure subroutine s_initialize_phasechange_module
393
394 end subroutine s_initialize_phasechange_module
395
396 !> Apply pT- or pTg-equilibrium relaxation with mass depletion based on the incoming state conditions.
397 subroutine s_infinite_relaxation_k(q_cons_vf)
398
399 type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
400 real(wp) :: ps !< equilibrium pressure
401 real(wp) :: ts !< equilibrium temperature
402 real(wp) :: rhoe, dyne, rhos !< total internal energy, kinetic energy, and total entropy
403 real(wp) :: rho, rm, m1, m2, mct !< total density, total reacting mass, individual reacting masses
404 real(wp) :: tvf !< total volume fraction
405
406# 66 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
407 real(wp), dimension(num_fluids) :: p_infpt, sk, hk, gk, ek, rhok
408# 68 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
409
410 !> Generic loop iterators
411 integer :: i, j, k, l
412
413#ifdef _CRAYFTN
414#ifdef MFC_OpenACC
415 ! CCE 19 IPA workaround: prevent bring_routine_resident SIGSEGV DIR$ NOINLINE s_infinite_pt_relaxation_k DIR$ NOINLINE
416 ! s_infinite_ptg_relaxation_k DIR$ NOINLINE s_correct_partial_densities
417#endif
418#endif
419
420 ! starting equilibrium solver
421
422
423# 81 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
424
425# 81 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
426#if defined(MFC_OpenACC)
427# 81 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
428!$acc parallel loop collapse(3) gang vector default(present) private(i, j, k, l, p_infpT, sk, hk, gk, ek, rhok, pS, TS, rhoe, dynE, rhos, rho, rM, m1, m2, MCT, TvF)
429# 81 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
430#elif defined(MFC_OpenMP)
431# 81 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
432
433# 81 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
434
435# 81 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
436
437# 81 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
438!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) &
439# 81 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
440!$omp& private(i, j, k, l, p_infpT, sk, hk, gk, ek, rhok, pS, TS, rhoe, dynE, rhos, rho, rM, m1, m2, MCT, TvF)
441# 81 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
442#endif
443# 83 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
444 do j = 0, m
445 do k = 0, n
446 do l = 0, p
447 rho = 0.0_wp; tvf = 0.0_wp
448
449# 87 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
450#if defined(MFC_OpenACC)
451# 87 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
452!$acc loop seq
453# 87 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
454#elif defined(MFC_OpenMP)
455# 87 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
456
457# 87 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
458#endif
459 do i = 1, num_fluids
460 ! Mixture density
461 rho = rho + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)
462
463 ! Total Volume Fraction
464 tvf = tvf + q_cons_vf(i + eqn_idx%adv%beg - 1)%sf(j, k, l)
465 end do
466
467 ! calculating the total reacting mass for the phase change process. By hypothesis, this should not change
468 ! throughout the phase-change process.
469 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)
470
471 ! correcting negative (reacting) mass fraction values in case they happen
472 call s_correct_partial_densities(mct, q_cons_vf, rm, j, k, l)
473
474 ! fixing m1 and m2 AFTER correcting the partial densities. Note that these values must be stored for the phase
475 ! change process that will happen a posteriori
476 m1 = q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l)
477
478 m2 = q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l)
479
480 ! kinetic energy as an auxiliary variable to the calculation of the total internal energy
481 dyne = 0.0_wp
482
483# 111 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
484#if defined(MFC_OpenACC)
485# 111 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
486!$acc loop seq
487# 111 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
488#elif defined(MFC_OpenMP)
489# 111 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
490
491# 111 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
492#endif
493 do i = eqn_idx%mom%beg, eqn_idx%mom%end
494 dyne = dyne + 5.0e-1_wp*q_cons_vf(i)%sf(j, k, l)**2/rho
495 end do
496
497 ! calculating the total energy that MUST be preserved throughout the pT- and pTg-relaxation procedures at each
498 ! of the cells. The internal energy is calculated as the total energy minus the kinetic energy to preserved its
499 ! value at sharp interfaces
500 rhoe = q_cons_vf(eqn_idx%E)%sf(j, k, l) - dyne
501
502 ! Calling pT-equilibrium for either finishing phase-change module, or as an IC for the pTg-equilibrium for this
503 ! case, MFL cannot be either 0 or 1, so I chose it to be 2
504 call s_infinite_pt_relaxation_k(j, k, l, 2, ps, p_infpt, q_cons_vf, rhoe, ts)
505
506 ! Check if pTg-equilibrium needed; only partial densities require updating
507 if ((relax_model == 6) .and. ((q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, &
508 & l) > mixm*rm) .and. (q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, &
509 & l) > mixm*rm)) .and. (ps < pcr) .and. (ts < tcr)) then
510 ! Solve pTg-equilibrium directly on the actual reacting masses. The Newton solver projects
511 ! the liquid mass onto [0, mT], so it recovers the single-phase limits itself (ml -> 0 for
512 ! all-vapor, ml -> mT for all-liquid). The former overheated-vapor / subcooled-liquid pT
513 ! shortcuts were removed: their pT states differ O(1) from the pTg equilibrium, so the
514 ! sub-ULP shortcut/pTg branch decision flipped across backends (CPU vs GPU) near a phase
515 ! boundary and destroyed cross-backend reproducibility.
516 q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) = m1
517 q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) = m2
518
519 call s_infinite_ptg_relaxation_k(j, k, l, ps, rhoe, q_cons_vf, ts)
520 end if
521
522 ! Calculations AFTER equilibrium
523
524
525# 143 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
526#if defined(MFC_OpenACC)
527# 143 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
528!$acc loop seq
529# 143 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
530#elif defined(MFC_OpenMP)
531# 143 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
532
533# 143 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
534#endif
535 do i = 1, num_fluids
536 ! entropy
537 sk(i) = cvs(i)*log((ts**isentrope_n(i))/((ps + isentrope_b(i))**(isentrope_n(i) - 1.0_wp))) + qvps(i)
538
539 ! enthalpy
540 hk(i) = isentrope_n(i)*cvs(i)*ts + qvs(i)
541
542 ! Gibbs-free energy
543 gk(i) = hk(i) - ts*sk(i)
544
545 ! densities
546 rhok(i) = f_sg_thermal(ps, ts, isentrope_n(i), isentrope_b(i), cvs(i))
547
548 ! internal energy
549 ek(i) = (ps + isentrope_n(i)*isentrope_b(i))/(ps + isentrope_b(i))*cvs(i)*ts + qvs(i)
550 end do
551
552 ! calculating volume fractions, internal energies, and total entropy
553 rhos = 0.0_wp
554
555# 163 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
556#if defined(MFC_OpenACC)
557# 163 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
558!$acc loop seq
559# 163 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
560#elif defined(MFC_OpenMP)
561# 163 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
562
563# 163 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
564#endif
565 do i = 1, num_fluids
566 ! volume fractions
567 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)
568
569 ! alpha*rho*e
570 if (model_eqns == model_eqns_6eq) then
571 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, &
572 & l)*ek(i)
573 end if
574
575 ! Total entropy
576 rhos = rhos + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*sk(i)
577 end do
578 end do
579 end do
580 end do
581
582# 180 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
583#if defined(MFC_OpenACC)
584# 180 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
585!$acc end parallel loop
586# 180 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
587#elif defined(MFC_OpenMP)
588# 180 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
589
590# 180 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
591!$omp end target teams loop
592# 180 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
593#endif
594
595 end subroutine s_infinite_relaxation_k
596
597 !> Apply pT-equilibrium relaxation for N fluids
598 !! @param MFL flag: 0=gas, 1=liquid, 2=mixture
599 subroutine s_infinite_pt_relaxation_k(j, k, l, MFL, pS, p_infpT, q_cons_vf, rhoe, TS)
600
601
602# 188 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
603#ifdef _CRAYFTN
604# 188 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
605#if MFC_OpenACC
606# 188 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
607!$acc routine seq
608# 188 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
609#elif MFC_OpenMP
610# 188 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
611
612# 188 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
613
614# 188 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
615!$omp declare target device_type(any)
616# 188 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
617#else
618# 188 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
619!DIR$ NOINLINE s_infinite_pt_relaxation_k
620# 188 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
621#endif
622# 188 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
623#elif MFC_OpenACC
624# 188 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
625!$acc routine seq
626# 188 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
627#elif MFC_OpenMP
628# 188 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
629
630# 188 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
631
632# 188 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
633!$omp declare target device_type(any)
634# 188 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
635#endif
636
637 ! initializing variables
638 integer, intent(in) :: j, k, l, MFL
639 real(wp), intent(out) :: pS
640 real(wp), dimension(1:), intent(out) :: p_infpT
641 type(scalar_field), dimension(sys_size), intent(in) :: q_cons_vf
642 real(wp), intent(in) :: rhoe
643 real(wp), intent(out) :: TS
644 real(wp) :: gp, gpp, hp, pO, mCP, mQ !< variables for the Newton Solver
645 real(wp) :: p_infpT_sum
646 integer :: i, ns !< generic loop iterators
647 ! auxiliary variables for the pT-equilibrium solver
648 mcp = 0.0_wp; mq = 0.0_wp; p_infpt_sum = 0._wp
649
650# 202 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
651#if defined(MFC_OpenACC)
652# 202 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
653!$acc loop seq
654# 202 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
655#elif defined(MFC_OpenMP)
656# 202 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
657
658# 202 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
659#endif
660 do i = 1, num_fluids
661 p_infpt(i) = isentrope_b(i)
662 p_infpt_sum = p_infpt_sum + abs(p_infpt(i))
663 end do
664 ! Performing tests before initializing the pT-equilibrium
665
666# 208 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
667#if defined(MFC_OpenACC)
668# 208 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
669!$acc loop seq
670# 208 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
671#elif defined(MFC_OpenMP)
672# 208 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
673
674# 208 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
675#endif
676 do i = 1, num_fluids
677 ! sum of the total alpha*rho*cp of the system
678 mcp = mcp + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*cvs(i)*isentrope_n(i)
679
680 ! sum of the total alpha*rho*q of the system
681 mq = mq + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*qvs(i)
682 end do
683
684# 225 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
685
686 ! Checking energy constraint
687 if ((rhoe - mq - minval(p_infpt)) < 0.0_wp) then
688 if ((mfl == 0) .or. (mfl == 1)) then
689 ! Assigning zero values for mass depletion cases pressure
690 ps = 0.0_wp
691
692 ! temperature
693 ts = 0.0_wp
694
695 return
696 end if
697 end if
698
699 ! calculating initial estimate for pressure in the pT-relaxation procedure. I will also use this variable to iterate over
700 ! the Newton's solver
701 po = 0.0_wp
702
703 ! Maybe improve this condition afterwards. As long as the initial guess is in between -min(isentrope_B) and infinity, a
704 ! solution
705 ! should be able to be found.
706 ps = 1.0e4_wp
707
708 ! Newton Solver for the pT-equilibrium
709 ns = 0
710 ! change this relative error metric. 1.e4_wp is just arbitrary
711 ! Relative criterion written in multiply form to avoid dividing by pO (pO = 0 on the first pass).
712 do while ((abs(ps - po) > palpha_eps) .and. (abs(ps - po) > (palpha_eps/1.e4_wp)*abs(po)) .or. (ns == 0))
713 ! increasing counter
714 ns = ns + 1
715 ! guard against non-convergence: accept the last iterate rather than looping forever
716 if (ns >= max_iter) exit
717
718 ! updating old pressure
719 po = ps
720
721 ! updating functions used in the Newton's solver
722 gpp = 0.0_wp; gp = 0.0_wp; hp = 0.0_wp
723
724# 263 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
725#if defined(MFC_OpenACC)
726# 263 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
727!$acc loop seq
728# 263 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
729#elif defined(MFC_OpenMP)
730# 263 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
731
732# 263 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
733#endif
734 do i = 1, num_fluids
735 gp = gp + (isentrope_n(i) - 1.0_wp)*q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, &
736 & l)*cvs(i)*(rhoe + ps - mq)/(mcp*(ps + p_infpt(i)))
737
738 gpp = gpp + (isentrope_n(i) - 1.0_wp)*q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, &
739 & l)*cvs(i)*(p_infpt(i) - rhoe + mq)/(mcp*(ps + p_infpt(i))**2)
740 end do
741
742 hp = 1.0_wp/(rhoe + ps - mq) + 1.0_wp/(ps + minval(p_infpt))
743
744 ! updating common pressure for the newton solver
745 ps = po + ((1.0_wp - gp)/gpp)/(1.0_wp - (1.0_wp - gp + abs(1.0_wp - gp))/(2.0_wp*gpp)*hp)
746 end do
747
748 ! common temperature; same closure as f_mixture_temperature in m_eos, written in conservative variables
749 ts = (rhoe + ps - mq)/mcp
750
751 end subroutine s_infinite_pt_relaxation_k
752
753 !> Evaluate the pTg-equilibrium residual R2D and temperature TS at a trial state (ml, pS) WITHOUT mutating q_cons_vf, so the
754 !! Newton driver can line-search. The total reacting mass mT is conserved, so the reacting masses are (ml, mT - ml) and only the
755 !! inert fluids are read from q_cons_vf. Also returns the mixture sums the Jacobian and the final temperature need.
756 subroutine s_compute_ptg_residual(ml, mT, pS, j, k, l, q_cons_vf, rhoe, R2D, TS, mCP, mQ, mCVGP, mCVGP2, mCPD)
757
758
759# 288 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
760#if MFC_OpenACC
761# 288 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
762!$acc routine seq
763# 288 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
764#elif MFC_OpenMP
765# 288 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
766
767# 288 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
768
769# 288 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
770!$omp declare target device_type(any)
771# 288 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
772#endif
773
774 real(wp), intent(in) :: ml, mT, pS, rhoe
775 integer, intent(in) :: j, k, l
776 type(scalar_field), dimension(sys_size), intent(in) :: q_cons_vf
777 real(wp), dimension(2), intent(out) :: R2D
778 real(wp), intent(out) :: TS, mCP, mQ, mCVGP, mCVGP2, mCPD
779 real(wp) :: mQD
780 integer :: i
781
782 ! reacting fluids contribute via (ml, mT - ml); inert fluids are summed from q_cons_vf
783 mcp = ml*cvs(lp)*isentrope_n(lp) + (mt - ml)*cvs(vp)*isentrope_n(vp)
784 mq = ml*qvs(lp) + (mt - ml)*qvs(vp)
785 mcvgp = 0.0_wp; mcvgp2 = 0.0_wp; mcpd = 0.0_wp; mqd = 0.0_wp
786
787# 302 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
788#if defined(MFC_OpenACC)
789# 302 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
790!$acc loop seq
791# 302 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
792#elif defined(MFC_OpenMP)
793# 302 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
794
795# 302 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
796#endif
797 do i = 1, num_fluids
798 if ((i /= lp) .and. (i /= vp)) then
799 mcp = mcp + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*cvs(i)*isentrope_n(i)
800 mq = mq + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*qvs(i)
801 mcvgp = mcvgp + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*cvs(i)*(isentrope_n(i) - 1)/(ps + isentrope_b(i))
802 mcvgp2 = mcvgp2 + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, &
803 & l)*cvs(i)*(isentrope_n(i) - 1)/((ps + isentrope_b(i))**2)
804 mqd = mqd + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*qvs(i)
805 mcpd = mcpd + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*cvs(i)*isentrope_n(i)
806 end if
807 end do
808
809 ts = 1.0_wp/(mt*cvs(vp)*(isentrope_n(vp) - 1)/(ps + isentrope_b(vp)) + ml*(cvs(lp)*(isentrope_n(lp) - 1)/(ps &
810 & + isentrope_b(lp)) - cvs(vp)*(isentrope_n(vp) - 1)/(ps + isentrope_b(vp))) + mcvgp)
811
812 ! (i) Gibbs free-energy equality
813 r2d(1) = ts*((cvs(lp)*isentrope_n(lp) - cvs(vp)*isentrope_n(vp))*(1 - log(ts)) - (qvps(lp) - qvps(vp)) + cvs(lp) &
814 & *(isentrope_n(lp) - 1)*log(ps + isentrope_b(lp)) - cvs(vp)*(isentrope_n(vp) - 1)*log(ps + isentrope_b(vp))) &
815 & + qvs(lp) - qvs(vp)
816
817 ! (ii) constant-energy condition
818 r2d(2) = rhoe + ps + ml*(qvs(vp) - qvs(lp)) - mt*qvs(vp) - mqd + (ml*(isentrope_n(vp)*cvs(vp) - isentrope_n(lp)*cvs(lp)) &
819 & - mt*isentrope_n(vp)*cvs(vp) - mcpd)/(ml*(cvs(lp)*(isentrope_n(lp) - 1)/(ps + isentrope_b(lp)) - cvs(vp) &
820 & *(isentrope_n(vp) - 1)/(ps + isentrope_b(vp))) + mt*cvs(vp)*(isentrope_n(vp) - 1)/(ps + isentrope_b(vp)) + mcvgp)
821
822 end subroutine s_compute_ptg_residual
823
824 !> Apply pTg-equilibrium relaxation: a damped (backtracking line search) Newton solve for the reacting liquid mass ml and
825 !! pressure pS enforcing Gibbs equality and energy conservation, converging on the residual norm (absolute ptgalpha_eps, or the
826 !! rhoe-relative branch). Every step is projected onto the physical bounds 0 <= ml <= mT, pS > pmin. This converges in a handful
827 !! of iterations with a bounded, uniform count (no GPU warp divergence), unlike the former fixed 1e-3 underrelaxation that
828 !! stalled far from the root.
829 subroutine s_infinite_ptg_relaxation_k(j, k, l, pS, rhoe, q_cons_vf, TS)
830
831
832# 337 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
833#ifdef _CRAYFTN
834# 337 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
835#if MFC_OpenACC
836# 337 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
837!$acc routine seq
838# 337 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
839#elif MFC_OpenMP
840# 337 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
841
842# 337 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
843
844# 337 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
845!$omp declare target device_type(any)
846# 337 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
847#else
848# 337 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
849!DIR$ NOINLINE s_infinite_ptg_relaxation_k
850# 337 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
851#endif
852# 337 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
853#elif MFC_OpenACC
854# 337 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
855!$acc routine seq
856# 337 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
857#elif MFC_OpenMP
858# 337 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
859
860# 337 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
861
862# 337 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
863!$omp declare target device_type(any)
864# 337 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
865#endif
866
867 integer, intent(in) :: j, k, l
868 real(wp), intent(inout) :: pS
869 real(wp), intent(in) :: rhoe
870 type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
871 real(wp), intent(inout) :: TS
872 real(wp), dimension(2, 2) :: Jac, InvJac
873 real(wp), dimension(2) :: R2D, R2D_try, DeltamP
874 real(wp) :: mCP, mCPD, mCVGP, mCVGP2, mQ
875 real(wp) :: ml, ml_try, mT, pS_try, pmin, lambda, resnorm, resnorm_try
876 real(wp) :: dFdT, dTdm, dTdp, detJ
877 integer :: ns, ls
878
879 ! total reacting mass is conserved; the liquid mass ml is the primary unknown, vapor mass = mT - ml
880 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)
881 ml = q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l)
882
883 ! recover a physical pressure guess when the incoming pS is non-physical
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 - isentrope_n(lp)*isentrope_b(lp)/(isentrope_n(lp) - 1))/qvs(lp)))) .or. ((ps >= 0.0_wp) &
887 & .and. (ps < 1.0e-1_wp))) then
888 ps = 1.0e4_wp
889 end if
890
891 ! pressure floor (stiffened gas requires pS + isentrope_B > 0 for both phases)
892 pmin = -min(isentrope_b(lp), isentrope_b(vp)) + 1.0_wp
893
894 call s_compute_ptg_residual(ml, mt, ps, j, k, l, q_cons_vf, rhoe, r2d, ts, mcp, mq, mcvgp, mcvgp2, mcpd)
895 resnorm = sqrt(r2d(1)**2 + r2d(2)**2)
896
897 do ns = 1, max_iter
898 ! converged on the absolute residual, or on the rhoe-relative residual (multiply form, rhoe > 0)
899 if ((resnorm <= ptgalpha_eps) .or. (resnorm <= (ptgalpha_eps/1.e6_wp)*rhoe)) exit
900
901 ! 2x2 Jacobian of (Gibbs equality, energy) with respect to (ml, pS) at the current state
902 dfdt = -(cvs(lp)*isentrope_n(lp) - cvs(vp)*isentrope_n(vp))*log(ts) - (qvps(lp) - qvps(vp)) + cvs(lp)*(isentrope_n(lp) &
903 & - 1)*log(ps + isentrope_b(lp)) - cvs(vp)*(isentrope_n(vp) - 1)*log(ps + isentrope_b(vp))
904 dtdm = -(cvs(lp)*(isentrope_n(lp) - 1)/(ps + isentrope_b(lp)) - cvs(vp)*(isentrope_n(vp) - 1)/(ps + isentrope_b(vp))) &
905 & *ts**2
906 dtdp = (mt*cvs(vp)*(isentrope_n(vp) - 1)/(ps + isentrope_b(vp))**2 + ml*(cvs(lp)*(isentrope_n(lp) - 1)/(ps &
907 & + isentrope_b(lp))**2 - cvs(vp)*(isentrope_n(vp) - 1)/(ps + isentrope_b(vp))**2) + mcvgp2)*ts**2
908
909 jac(1, 1) = dfdt*dtdm
910 jac(1, &
911 & 2) = dfdt*dtdp + ts*(cvs(lp)*(isentrope_n(lp) - 1)/(ps + isentrope_b(lp)) - cvs(vp)*(isentrope_n(vp) - 1)/(ps &
912 & + isentrope_b(vp)))
913 jac(2, &
914 & 1) = qvs(vp) - qvs(lp) + (cvs(vp)*isentrope_n(vp) - cvs(lp)*isentrope_n(lp))/(ml*(cvs(lp)*(isentrope_n(lp) - 1) &
915 & /(ps + isentrope_b(lp)) - cvs(vp)*(isentrope_n(vp) - 1)/(ps + isentrope_b(vp))) + mt*cvs(vp)*(isentrope_n(vp) &
916 & - 1)/(ps + isentrope_b(vp)) + mcvgp) - (ml*(cvs(vp)*isentrope_n(vp) - cvs(lp)*isentrope_n(lp)) - mt*cvs(vp) &
917 & *isentrope_n(vp) - mcpd)*(cvs(lp)*(isentrope_n(lp) - 1)/(ps + isentrope_b(lp)) - cvs(vp)*(isentrope_n(vp) - 1) &
918 & /(ps + isentrope_b(vp)))/((ml*(cvs(lp)*(isentrope_n(lp) - 1)/(ps + isentrope_b(lp)) - cvs(vp)*(isentrope_n(vp) &
919 & - 1)/(ps + isentrope_b(vp))) + mt*cvs(vp)*(isentrope_n(vp) - 1)/(ps + isentrope_b(vp)) + mcvgp)**2)
920 jac(2, &
921 & 2) = 1 + (ml*(cvs(vp)*isentrope_n(vp) - cvs(lp)*isentrope_n(lp)) - mt*cvs(vp)*isentrope_n(vp) - mcpd) &
922 & *(ml*(cvs(lp)*(isentrope_n(lp) - 1)/(ps + isentrope_b(lp))**2 - cvs(vp)*(isentrope_n(vp) - 1)/(ps &
923 & + isentrope_b(vp))**2) + mt*cvs(vp)*(isentrope_n(vp) - 1)/(ps + isentrope_b(vp))**2 + mcvgp2)/(ml*(cvs(lp) &
924 & *(isentrope_n(lp) - 1)/(ps + isentrope_b(lp)) - cvs(vp)*(isentrope_n(vp) - 1)/(ps + isentrope_b(vp))) &
925 & + mt*cvs(vp)*(isentrope_n(vp) - 1)/(ps + isentrope_b(vp)) + mcvgp)**2
926
927 detj = jac(1, 1)*jac(2, 2) - jac(1, 2)*jac(2, 1)
928 ! singular Jacobian: no usable Newton direction, accept the current (best) state
929 if (detj == 0.0_wp) exit
930
931 invjac(1, 1) = jac(2, 2)/detj
932 invjac(1, 2) = -jac(1, 2)/detj
933 invjac(2, 1) = -jac(2, 1)/detj
934 invjac(2, 2) = jac(1, 1)/detj
935
936 deltamp(1) = -(invjac(1, 1)*r2d(1) + invjac(1, 2)*r2d(2))
937 deltamp(2) = -(invjac(2, 1)*r2d(1) + invjac(2, 2)*r2d(2))
938
939 ! backtracking line search: halve the step until the residual decreases, keeping the state
940 ! physical (0 <= ml <= mT, pS above the stiffened-gas floor)
941 lambda = 1.0_wp
942 do ls = 1, ptg_ls_max
943 ml_try = min(max(ml + lambda*deltamp(1), 0.0_wp), mt)
944 ps_try = max(ps + lambda*deltamp(2), pmin)
945 call s_compute_ptg_residual(ml_try, mt, ps_try, j, k, l, q_cons_vf, rhoe, r2d_try, ts, mcp, mq, mcvgp, mcvgp2, mcpd)
946 resnorm_try = sqrt(r2d_try(1)**2 + r2d_try(2)**2)
947 if ((resnorm_try < resnorm) .or. (ls == ptg_ls_max)) exit
948 lambda = 0.5_wp*lambda
949 end do
950
951 ! accept the trial state (TS, mCP, mQ, mCVGP, mCVGP2, mCPD already set to it by the last call)
952 ml = ml_try; ps = ps_try; r2d = r2d_try; resnorm = resnorm_try
953 end do
954
955 ! commit the reacting masses (mT conserved) and set the common temperature
956 q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) = ml
957 q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) = mt - ml
958
959 ts = (rhoe + ps - mq)/mcp
960
961 end subroutine s_infinite_ptg_relaxation_k
962
963 !> Correct the partial densities of the reacting fluids in case one of them is negative but their sum is positive. Inert phases
964 !! are not corrected at this moment
965 subroutine s_correct_partial_densities(MCT, q_cons_vf, rM, j, k, l)
966
967
968# 439 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
969#ifdef _CRAYFTN
970# 439 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
971#if MFC_OpenACC
972# 439 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
973!$acc routine seq
974# 439 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
975#elif MFC_OpenMP
976# 439 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
977
978# 439 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
979
980# 439 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
981!$omp declare target device_type(any)
982# 439 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
983#else
984# 439 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
985!DIR$ NOINLINE s_correct_partial_densities
986# 439 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
987#endif
988# 439 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
989#elif MFC_OpenACC
990# 439 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
991!$acc routine seq
992# 439 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
993#elif MFC_OpenMP
994# 439 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
995
996# 439 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
997
998# 439 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
999!$omp declare target device_type(any)
1000# 439 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1001#endif
1002
1003 !> @name variables for the correction of the reacting partial densities
1004 !> @{
1005 real(wp), intent(out) :: mct
1006 type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
1007 real(wp), intent(inout) :: rm
1008 integer, intent(in) :: j, k, l
1009 !> @}
1010 if (rm < 0.0_wp) then
1011 if ((q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, &
1012 & l) >= -1.0_wp*mixm) .and. (q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) >= -1.0_wp*mixm)) then
1013 q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) = 0.0_wp
1014
1015 q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) = 0.0_wp
1016
1017 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)
1018 end if
1019 end if
1020
1021 ! TODO: Consider partitioning partial densities instead of absolute-value correction
1022 mct = 2*mixm
1023
1024 ! correcting the partial densities of the reacting fluids. What to do for the nonreacting ones?
1025 if (q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) < 0.0_wp) then
1026 q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) = mct*rm
1027
1028 q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) = (1.0_wp - mct)*rm
1029 else if (q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) < 0.0_wp) then
1030 q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) = (1.0_wp - mct)*rm
1031
1032 q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) = mct*rm
1033 end if
1034
1035 end subroutine s_correct_partial_densities
1036
1037 !> Finalize the phase change module
1040 end subroutine s_finalize_relaxation_solver_module
1041
1042end module m_phase_change
q_cons_vf
type(scalar_field), dimension(sys_size), intent(inout) q_cons_vf
Definition m_phase_change.fpp.f90:1006
rm
real(wp), intent(inout) rm
Definition m_phase_change.fpp.f90:1007
k
integer, intent(in) k
Definition m_phase_change.fpp.f90:1008
mct
real(wp), intent(out) mct
Definition m_phase_change.fpp.f90:1005
j
integer, intent(in) j
Definition m_phase_change.fpp.f90:1008
l
integer, intent(in) l
Definition m_phase_change.fpp.f90:1008
m_constants
Compile-time constant parameters: default values, tolerances, and physical constants.
Definition m_constants.fpp.f90:7
m_constants::model_eqns_6eq
integer, parameter model_eqns_6eq
Definition m_constants.fpp.f90:149
m_derived_types
Shared derived types for field data, patch geometry, bubble dynamics, and MPI I/O structures.
Definition m_derived_types.fpp.f90:335
m_eos
Equations of state in Gamma/Pi form, rho e = Gamma(rho) p + Pi(rho).
Definition m_eos.fpp.f90:381
m_eos::f_sg_thermal
real(wp) function, public f_sg_thermal(pres, rho_or_t, n, b, cv)
Stiffened-gas thermal law p + B = (n - 1)*cv*rho*T. Pass rho to get T, or T to get rho.
Definition m_eos.fpp.f90:1183
m_global_parameters
Global parameters for the computational domain, fluid properties, and simulation algorithm configurat...
Definition m_global_parameters.fpp.f90:346
m_helper_basic
Basic floating-point utilities: approximate equality, default detection, and coordinate bounds.
Definition m_helper_basic.fpp.f90:335
m_mpi_proxy
MPI halo exchange, domain decomposition, and buffer packing/unpacking for the simulation solver.
Definition m_mpi_proxy.fpp.f90:346
m_phase_change
Phase transition relaxation solvers for liquid-vapor flows with cavitation and boiling.
Definition m_phase_change.fpp.f90:346
m_phase_change::s_finalize_relaxation_solver_module
impure subroutine, public s_finalize_relaxation_solver_module
Finalize the phase change module.
Definition m_phase_change.fpp.f90:1039
m_phase_change::s_infinite_relaxation_k
subroutine, public s_infinite_relaxation_k(q_cons_vf)
Apply pT- or pTg-equilibrium relaxation with mass depletion based on the incoming state conditions.
Definition m_phase_change.fpp.f90:398
m_phase_change::vp
integer, parameter vp
index for the vapor phase of the reacting fluid
Definition m_phase_change.fpp.f90:370
m_phase_change::lp
integer, parameter lp
index for the liquid phase of the reacting fluid
Definition m_phase_change.fpp.f90:369
m_phase_change::s_infinite_pt_relaxation_k
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.
Definition m_phase_change.fpp.f90:600
m_phase_change::s_relaxation_solver
impure subroutine, public s_relaxation_solver(q_cons_vf)
Dispatch to the correct relaxation solver. Replaces the procedure pointer, which CCE is breaking on.
Definition m_phase_change.fpp.f90:377
m_phase_change::tcr
real(wp), parameter tcr
Critical temperature of water [K].
Definition m_phase_change.fpp.f90:366
m_phase_change::s_compute_ptg_residual
subroutine s_compute_ptg_residual(ml, mt, ps, j, k, l, q_cons_vf, rhoe, r2d, ts, mcp, mq, mcvgp, mcvgp2, mcpd)
Evaluate the pTg-equilibrium residual R2D and temperature TS at a trial state (ml,...
Definition m_phase_change.fpp.f90:757
m_phase_change::ptg_ls_max
integer, parameter ptg_ls_max
max backtracking-line-search halvings in the pTg solver
Definition m_phase_change.fpp.f90:367
m_phase_change::max_iter
integer, parameter max_iter
max Newton iterations before accepting the last iterate
Definition m_phase_change.fpp.f90:364
m_phase_change::s_initialize_phasechange_module
impure subroutine, public s_initialize_phasechange_module
Initialize the phase change module (no module-level state to set up; the pT/pTg relaxation solvers ar...
Definition m_phase_change.fpp.f90:393
m_phase_change::mixm
real(wp), parameter mixm
Mixture mass fraction threshold for triggering phase change.
Definition m_phase_change.fpp.f90:368
m_phase_change::s_infinite_ptg_relaxation_k
subroutine s_infinite_ptg_relaxation_k(j, k, l, ps, rhoe, q_cons_vf, ts)
Apply pTg-equilibrium relaxation: a damped (backtracking line search) Newton solve for the reacting l...
Definition m_phase_change.fpp.f90:830
m_phase_change::pcr
real(wp), parameter pcr
Critical pressure of water [Pa].
Definition m_phase_change.fpp.f90:365
m_variables_conversion
Conservative-to-primitive variable conversion, mixture property evaluation, and pressure computation.
Definition m_variables_conversion.fpp.f90:346
m_derived_types::scalar_field
Derived type annexing a scalar field (SF).
Definition m_derived_types.fpp.f90:349