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 ieee_arithmetic
353 use m_helper_basic
354 use m_constants, only: model_eqns_6eq
355
356 implicit none
357
358 private
359 public :: s_initialize_phasechange_module, s_relaxation_solver, s_infinite_relaxation_k, s_finalize_relaxation_solver_module
360
361 !> @name Parameters for the first order transition phase change
362 !> @{
363 integer, parameter :: max_iter = 100000 !< max Newton iterations before accepting the last iterate
364 real(wp), parameter :: pcr = 4.94e7_wp !< Critical pressure of water [Pa]
365 real(wp), parameter :: tcr = 385.05_wp + 273.15_wp !< Critical temperature of water [K]
366 integer, parameter :: ptg_ls_max = 30 !< max backtracking-line-search halvings in the pTg solver
367 real(wp), parameter :: mixm = 1.0e-8_wp !< Mixture mass fraction threshold for triggering phase change
368 integer, parameter :: lp = 1 !< index for the liquid phase of the reacting fluid
369 integer, parameter :: vp = 2 !< index for the vapor phase of the reacting fluid
370 !> @}
371
372contains
373
374 !> Dispatch to the correct relaxation solver. Replaces the procedure pointer, which CCE is breaking on.
375 impure subroutine s_relaxation_solver(q_cons_vf)
376
377 type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
378 ! This is empty because in current master the procedure pointer was never assigned
379
380 if (.not. (.false.)) then
381# 43 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
382 call s_mpi_abort("m_phase_change.fpp:43: " // "Assertion failed: .false.. " &
383# 43 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
384 & // "s_relaxation_solver called but it currently does nothing")
385# 43 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
386 end if
387
388 end subroutine s_relaxation_solver
389
390 !> Initialize the phase change module (no module-level state to set up; the pT/pTg relaxation solvers are self-contained)
391 impure subroutine s_initialize_phasechange_module
392
393 end subroutine s_initialize_phasechange_module
394
395 !> Apply pT- or pTg-equilibrium relaxation with mass depletion based on the incoming state conditions.
396 subroutine s_infinite_relaxation_k(q_cons_vf)
397
398 type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
399 real(wp) :: ps !< equilibrium pressure
400 real(wp) :: ts !< equilibrium temperature
401 real(wp) :: rhoe, dyne, rhos !< total internal energy, kinetic energy, and total entropy
402 real(wp) :: rho, rm, m1, m2, mct !< total density, total reacting mass, individual reacting masses
403 real(wp) :: tvf !< total volume fraction
404
405# 65 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
406 real(wp), dimension(num_fluids) :: p_infpt, sk, hk, gk, ek, rhok
407# 67 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
408
409 !> Generic loop iterators
410 integer :: i, j, k, l
411
412#ifdef _CRAYFTN
413#ifdef MFC_OpenACC
414 ! CCE 19 IPA workaround: prevent bring_routine_resident SIGSEGV DIR$ NOINLINE s_infinite_pt_relaxation_k DIR$ NOINLINE
415 ! s_infinite_ptg_relaxation_k DIR$ NOINLINE s_correct_partial_densities
416#endif
417#endif
418
419 ! starting equilibrium solver
420
421
422# 80 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
423
424# 80 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
425#if defined(MFC_OpenACC)
426# 80 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
427!$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)
428# 80 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
429#elif defined(MFC_OpenMP)
430# 80 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
431
432# 80 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
433
434# 80 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
435
436# 80 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
437!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) &
438# 80 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
439!$omp& private(i, j, k, l, p_infpT, sk, hk, gk, ek, rhok, pS, TS, rhoe, dynE, rhos, rho, rM, m1, m2, MCT, TvF)
440# 80 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
441#endif
442# 82 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
443 do j = 0, m
444 do k = 0, n
445 do l = 0, p
446 rho = 0.0_wp; tvf = 0.0_wp
447
448# 86 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
449#if defined(MFC_OpenACC)
450# 86 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
451!$acc loop seq
452# 86 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
453#elif defined(MFC_OpenMP)
454# 86 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
455
456# 86 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
457#endif
458 do i = 1, num_fluids
459 ! Mixture density
460 rho = rho + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)
461
462 ! Total Volume Fraction
463 tvf = tvf + q_cons_vf(i + eqn_idx%adv%beg - 1)%sf(j, k, l)
464 end do
465
466 ! calculating the total reacting mass for the phase change process. By hypothesis, this should not change
467 ! throughout the phase-change process.
468 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)
469
470 ! correcting negative (reacting) mass fraction values in case they happen
471 call s_correct_partial_densities(mct, q_cons_vf, rm, j, k, l)
472
473 ! fixing m1 and m2 AFTER correcting the partial densities. Note that these values must be stored for the phase
474 ! change process that will happen a posteriori
475 m1 = q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l)
476
477 m2 = q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l)
478
479 ! kinetic energy as an auxiliary variable to the calculation of the total internal energy
480 dyne = 0.0_wp
481
482# 110 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
483#if defined(MFC_OpenACC)
484# 110 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
485!$acc loop seq
486# 110 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
487#elif defined(MFC_OpenMP)
488# 110 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
489
490# 110 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
491#endif
492 do i = eqn_idx%mom%beg, eqn_idx%mom%end
493 dyne = dyne + 5.0e-1_wp*q_cons_vf(i)%sf(j, k, l)**2/rho
494 end do
495
496 ! calculating the total energy that MUST be preserved throughout the pT- and pTg-relaxation procedures at each
497 ! of the cells. The internal energy is calculated as the total energy minus the kinetic energy to preserved its
498 ! value at sharp interfaces
499 rhoe = q_cons_vf(eqn_idx%E)%sf(j, k, l) - dyne
500
501 ! Calling pT-equilibrium for either finishing phase-change module, or as an IC for the pTg-equilibrium for this
502 ! case, MFL cannot be either 0 or 1, so I chose it to be 2
503 call s_infinite_pt_relaxation_k(j, k, l, 2, ps, p_infpt, q_cons_vf, rhoe, ts)
504
505 ! Check if pTg-equilibrium needed; only partial densities require updating
506 if ((relax_model == 6) .and. ((q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, &
507 & l) > mixm*rm) .and. (q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, &
508 & l) > mixm*rm)) .and. (ps < pcr) .and. (ts < tcr)) then
509 ! Solve pTg-equilibrium directly on the actual reacting masses. The Newton solver projects
510 ! the liquid mass onto [0, mT], so it recovers the single-phase limits itself (ml -> 0 for
511 ! all-vapor, ml -> mT for all-liquid). The former overheated-vapor / subcooled-liquid pT
512 ! shortcuts were removed: their pT states differ O(1) from the pTg equilibrium, so the
513 ! sub-ULP shortcut/pTg branch decision flipped across backends (CPU vs GPU) near a phase
514 ! boundary and destroyed cross-backend reproducibility.
515 q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) = m1
516 q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) = m2
517
518 call s_infinite_ptg_relaxation_k(j, k, l, ps, rhoe, q_cons_vf, ts)
519 end if
520
521 ! Calculations AFTER equilibrium
522
523
524# 142 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
525#if defined(MFC_OpenACC)
526# 142 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
527!$acc loop seq
528# 142 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
529#elif defined(MFC_OpenMP)
530# 142 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
531
532# 142 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
533#endif
534 do i = 1, num_fluids
535 ! entropy
536 sk(i) = cvs(i)*log((ts**isentrope_n(i))/((ps + isentrope_b(i))**(isentrope_n(i) - 1.0_wp))) + qvps(i)
537
538 ! enthalpy
539 hk(i) = isentrope_n(i)*cvs(i)*ts + qvs(i)
540
541 ! Gibbs-free energy
542 gk(i) = hk(i) - ts*sk(i)
543
544 ! densities
545 rhok(i) = f_sg_thermal(ps, ts, isentrope_n(i), isentrope_b(i), cvs(i))
546
547 ! internal energy
548 ek(i) = (ps + isentrope_n(i)*isentrope_b(i))/(ps + isentrope_b(i))*cvs(i)*ts + qvs(i)
549 end do
550
551 ! calculating volume fractions, internal energies, and total entropy
552 rhos = 0.0_wp
553
554# 162 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
555#if defined(MFC_OpenACC)
556# 162 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
557!$acc loop seq
558# 162 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
559#elif defined(MFC_OpenMP)
560# 162 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
561
562# 162 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
563#endif
564 do i = 1, num_fluids
565 ! volume fractions
566 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)
567
568 ! alpha*rho*e
569 if (model_eqns == model_eqns_6eq) then
570 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, &
571 & l)*ek(i)
572 end if
573
574 ! Total entropy
575 rhos = rhos + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*sk(i)
576 end do
577 end do
578 end do
579 end do
580
581# 179 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
582#if defined(MFC_OpenACC)
583# 179 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
584!$acc end parallel loop
585# 179 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
586#elif defined(MFC_OpenMP)
587# 179 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
588
589# 179 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
590!$omp end target teams loop
591# 179 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
592#endif
593
594 end subroutine s_infinite_relaxation_k
595
596 !> Apply pT-equilibrium relaxation for N fluids
597 !! @param MFL flag: 0=gas, 1=liquid, 2=mixture
598 subroutine s_infinite_pt_relaxation_k(j, k, l, MFL, pS, p_infpT, q_cons_vf, rhoe, TS)
599
600
601# 187 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
602#ifdef _CRAYFTN
603# 187 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
604#if MFC_OpenACC
605# 187 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
606!$acc routine seq
607# 187 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
608#elif MFC_OpenMP
609# 187 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
610
611# 187 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
612
613# 187 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
614!$omp declare target device_type(any)
615# 187 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
616#else
617# 187 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
618!DIR$ NOINLINE s_infinite_pt_relaxation_k
619# 187 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
620#endif
621# 187 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
622#elif MFC_OpenACC
623# 187 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
624!$acc routine seq
625# 187 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
626#elif MFC_OpenMP
627# 187 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
628
629# 187 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
630
631# 187 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
632!$omp declare target device_type(any)
633# 187 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
634#endif
635
636 ! initializing variables
637 integer, intent(in) :: j, k, l, MFL
638 real(wp), intent(out) :: pS
639 real(wp), dimension(1:), intent(out) :: p_infpT
640 type(scalar_field), dimension(sys_size), intent(in) :: q_cons_vf
641 real(wp), intent(in) :: rhoe
642 real(wp), intent(out) :: TS
643 real(wp) :: gp, gpp, hp, pO, mCP, mQ !< variables for the Newton Solver
644 real(wp) :: p_infpT_sum
645 integer :: i, ns !< generic loop iterators
646 ! auxiliary variables for the pT-equilibrium solver
647 mcp = 0.0_wp; mq = 0.0_wp; p_infpt_sum = 0._wp
648
649# 201 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
650#if defined(MFC_OpenACC)
651# 201 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
652!$acc loop seq
653# 201 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
654#elif defined(MFC_OpenMP)
655# 201 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
656
657# 201 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
658#endif
659 do i = 1, num_fluids
660 p_infpt(i) = isentrope_b(i)
661 p_infpt_sum = p_infpt_sum + abs(p_infpt(i))
662 end do
663 ! Performing tests before initializing the pT-equilibrium
664
665# 207 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
666#if defined(MFC_OpenACC)
667# 207 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
668!$acc loop seq
669# 207 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
670#elif defined(MFC_OpenMP)
671# 207 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
672
673# 207 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
674#endif
675 do i = 1, num_fluids
676 ! sum of the total alpha*rho*cp of the system
677 mcp = mcp + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*cvs(i)*isentrope_n(i)
678
679 ! sum of the total alpha*rho*q of the system
680 mq = mq + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*qvs(i)
681 end do
682
683# 224 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
684
685 ! Checking energy constraint
686 if ((rhoe - mq - minval(p_infpt)) < 0.0_wp) then
687 if ((mfl == 0) .or. (mfl == 1)) then
688 ! Assigning zero values for mass depletion cases pressure
689 ps = 0.0_wp
690
691 ! temperature
692 ts = 0.0_wp
693
694 return
695 end if
696 end if
697
698 ! calculating initial estimate for pressure in the pT-relaxation procedure. I will also use this variable to iterate over
699 ! the Newton's solver
700 po = 0.0_wp
701
702 ! Maybe improve this condition afterwards. As long as the initial guess is in between -min(isentrope_B) and infinity, a
703 ! solution
704 ! should be able to be found.
705 ps = 1.0e4_wp
706
707 ! Newton Solver for the pT-equilibrium
708 ns = 0
709 ! change this relative error metric. 1.e4_wp is just arbitrary
710 ! Relative criterion written in multiply form to avoid dividing by pO (pO = 0 on the first pass).
711 do while ((abs(ps - po) > palpha_eps) .and. (abs(ps - po) > (palpha_eps/1.e4_wp)*abs(po)) .or. (ns == 0))
712 ! increasing counter
713 ns = ns + 1
714 ! guard against non-convergence: accept the last iterate rather than looping forever
715 if (ns >= max_iter) exit
716
717 ! updating old pressure
718 po = ps
719
720 ! updating functions used in the Newton's solver
721 gpp = 0.0_wp; gp = 0.0_wp; hp = 0.0_wp
722
723# 262 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
724#if defined(MFC_OpenACC)
725# 262 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
726!$acc loop seq
727# 262 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
728#elif defined(MFC_OpenMP)
729# 262 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
730
731# 262 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
732#endif
733 do i = 1, num_fluids
734 gp = gp + (isentrope_n(i) - 1.0_wp)*q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, &
735 & l)*cvs(i)*(rhoe + ps - mq)/(mcp*(ps + p_infpt(i)))
736
737 gpp = gpp + (isentrope_n(i) - 1.0_wp)*q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, &
738 & l)*cvs(i)*(p_infpt(i) - rhoe + mq)/(mcp*(ps + p_infpt(i))**2)
739 end do
740
741 hp = 1.0_wp/(rhoe + ps - mq) + 1.0_wp/(ps + minval(p_infpt))
742
743 ! updating common pressure for the newton solver
744 ps = po + ((1.0_wp - gp)/gpp)/(1.0_wp - (1.0_wp - gp + abs(1.0_wp - gp))/(2.0_wp*gpp)*hp)
745 end do
746
747 ! common temperature; same closure as f_mixture_temperature in m_variables_conversion, written in conservative variables
748 ts = (rhoe + ps - mq)/mcp
749
750 end subroutine s_infinite_pt_relaxation_k
751
752 !> Evaluate the pTg-equilibrium residual R2D and temperature TS at a trial state (ml, pS) WITHOUT mutating q_cons_vf, so the
753 !! Newton driver can line-search. The total reacting mass mT is conserved, so the reacting masses are (ml, mT - ml) and only the
754 !! inert fluids are read from q_cons_vf. Also returns the mixture sums the Jacobian and the final temperature need.
755 subroutine s_compute_ptg_residual(ml, mT, pS, j, k, l, q_cons_vf, rhoe, R2D, TS, mCP, mQ, mCVGP, mCVGP2, mCPD)
756
757
758# 287 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
759#if MFC_OpenACC
760# 287 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
761!$acc routine seq
762# 287 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
763#elif MFC_OpenMP
764# 287 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
765
766# 287 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
767
768# 287 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
769!$omp declare target device_type(any)
770# 287 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
771#endif
772
773 real(wp), intent(in) :: ml, mT, pS, rhoe
774 integer, intent(in) :: j, k, l
775 type(scalar_field), dimension(sys_size), intent(in) :: q_cons_vf
776 real(wp), dimension(2), intent(out) :: R2D
777 real(wp), intent(out) :: TS, mCP, mQ, mCVGP, mCVGP2, mCPD
778 real(wp) :: mQD
779 integer :: i
780
781 ! reacting fluids contribute via (ml, mT - ml); inert fluids are summed from q_cons_vf
782 mcp = ml*cvs(lp)*isentrope_n(lp) + (mt - ml)*cvs(vp)*isentrope_n(vp)
783 mq = ml*qvs(lp) + (mt - ml)*qvs(vp)
784 mcvgp = 0.0_wp; mcvgp2 = 0.0_wp; mcpd = 0.0_wp; mqd = 0.0_wp
785
786# 301 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
787#if defined(MFC_OpenACC)
788# 301 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
789!$acc loop seq
790# 301 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
791#elif defined(MFC_OpenMP)
792# 301 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
793
794# 301 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
795#endif
796 do i = 1, num_fluids
797 if ((i /= lp) .and. (i /= vp)) then
798 mcp = mcp + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*cvs(i)*isentrope_n(i)
799 mq = mq + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*qvs(i)
800 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))
801 mcvgp2 = mcvgp2 + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, &
802 & l)*cvs(i)*(isentrope_n(i) - 1)/((ps + isentrope_b(i))**2)
803 mqd = mqd + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*qvs(i)
804 mcpd = mcpd + q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)*cvs(i)*isentrope_n(i)
805 end if
806 end do
807
808 ts = 1.0_wp/(mt*cvs(vp)*(isentrope_n(vp) - 1)/(ps + isentrope_b(vp)) + ml*(cvs(lp)*(isentrope_n(lp) - 1)/(ps &
809 & + isentrope_b(lp)) - cvs(vp)*(isentrope_n(vp) - 1)/(ps + isentrope_b(vp))) + mcvgp)
810
811 ! (i) Gibbs free-energy equality
812 r2d(1) = ts*((cvs(lp)*isentrope_n(lp) - cvs(vp)*isentrope_n(vp))*(1 - log(ts)) - (qvps(lp) - qvps(vp)) + cvs(lp) &
813 & *(isentrope_n(lp) - 1)*log(ps + isentrope_b(lp)) - cvs(vp)*(isentrope_n(vp) - 1)*log(ps + isentrope_b(vp))) &
814 & + qvs(lp) - qvs(vp)
815
816 ! (ii) constant-energy condition
817 r2d(2) = rhoe + ps + ml*(qvs(vp) - qvs(lp)) - mt*qvs(vp) - mqd + (ml*(isentrope_n(vp)*cvs(vp) - isentrope_n(lp)*cvs(lp)) &
818 & - mt*isentrope_n(vp)*cvs(vp) - mcpd)/(ml*(cvs(lp)*(isentrope_n(lp) - 1)/(ps + isentrope_b(lp)) - cvs(vp) &
819 & *(isentrope_n(vp) - 1)/(ps + isentrope_b(vp))) + mt*cvs(vp)*(isentrope_n(vp) - 1)/(ps + isentrope_b(vp)) + mcvgp)
820
821 end subroutine s_compute_ptg_residual
822
823 !> Apply pTg-equilibrium relaxation: a damped (backtracking line search) Newton solve for the reacting liquid mass ml and
824 !! pressure pS enforcing Gibbs equality and energy conservation, converging on the residual norm (absolute ptgalpha_eps, or the
825 !! rhoe-relative branch). Every step is projected onto the physical bounds 0 <= ml <= mT, pS > pmin. This converges in a handful
826 !! of iterations with a bounded, uniform count (no GPU warp divergence), unlike the former fixed 1e-3 underrelaxation that
827 !! stalled far from the root.
828 subroutine s_infinite_ptg_relaxation_k(j, k, l, pS, rhoe, q_cons_vf, TS)
829
830
831# 336 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
832#ifdef _CRAYFTN
833# 336 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
834#if MFC_OpenACC
835# 336 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
836!$acc routine seq
837# 336 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
838#elif MFC_OpenMP
839# 336 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
840
841# 336 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
842
843# 336 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
844!$omp declare target device_type(any)
845# 336 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
846#else
847# 336 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
848!DIR$ NOINLINE s_infinite_ptg_relaxation_k
849# 336 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
850#endif
851# 336 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
852#elif MFC_OpenACC
853# 336 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
854!$acc routine seq
855# 336 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
856#elif MFC_OpenMP
857# 336 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
858
859# 336 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
860
861# 336 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
862!$omp declare target device_type(any)
863# 336 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
864#endif
865
866 integer, intent(in) :: j, k, l
867 real(wp), intent(inout) :: pS
868 real(wp), intent(in) :: rhoe
869 type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
870 real(wp), intent(inout) :: TS
871 real(wp), dimension(2, 2) :: Jac, InvJac
872 real(wp), dimension(2) :: R2D, R2D_try, DeltamP
873 real(wp) :: mCP, mCPD, mCVGP, mCVGP2, mQ
874 real(wp) :: ml, ml_try, mT, pS_try, pmin, lambda, resnorm, resnorm_try
875 real(wp) :: dFdT, dTdm, dTdp, detJ
876 integer :: ns, ls
877
878 ! total reacting mass is conserved; the liquid mass ml is the primary unknown, vapor mass = mT - ml
879 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)
880 ml = q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l)
881
882 ! recover a physical pressure guess when the incoming pS is non-physical
883 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, &
884 & k, &
885 & l)) > ((rhoe - isentrope_n(lp)*isentrope_b(lp)/(isentrope_n(lp) - 1))/qvs(lp)))) .or. ((ps >= 0.0_wp) &
886 & .and. (ps < 1.0e-1_wp))) then
887 ps = 1.0e4_wp
888 end if
889
890 ! pressure floor (stiffened gas requires pS + isentrope_B > 0 for both phases)
891 pmin = -min(isentrope_b(lp), isentrope_b(vp)) + 1.0_wp
892
893 call s_compute_ptg_residual(ml, mt, ps, j, k, l, q_cons_vf, rhoe, r2d, ts, mcp, mq, mcvgp, mcvgp2, mcpd)
894 resnorm = sqrt(r2d(1)**2 + r2d(2)**2)
895
896 do ns = 1, max_iter
897 ! converged on the absolute residual, or on the rhoe-relative residual (multiply form, rhoe > 0)
898 if ((resnorm <= ptgalpha_eps) .or. (resnorm <= (ptgalpha_eps/1.e6_wp)*rhoe)) exit
899
900 ! 2x2 Jacobian of (Gibbs equality, energy) with respect to (ml, pS) at the current state
901 dfdt = -(cvs(lp)*isentrope_n(lp) - cvs(vp)*isentrope_n(vp))*log(ts) - (qvps(lp) - qvps(vp)) + cvs(lp)*(isentrope_n(lp) &
902 & - 1)*log(ps + isentrope_b(lp)) - cvs(vp)*(isentrope_n(vp) - 1)*log(ps + isentrope_b(vp))
903 dtdm = -(cvs(lp)*(isentrope_n(lp) - 1)/(ps + isentrope_b(lp)) - cvs(vp)*(isentrope_n(vp) - 1)/(ps + isentrope_b(vp))) &
904 & *ts**2
905 dtdp = (mt*cvs(vp)*(isentrope_n(vp) - 1)/(ps + isentrope_b(vp))**2 + ml*(cvs(lp)*(isentrope_n(lp) - 1)/(ps &
906 & + isentrope_b(lp))**2 - cvs(vp)*(isentrope_n(vp) - 1)/(ps + isentrope_b(vp))**2) + mcvgp2)*ts**2
907
908 jac(1, 1) = dfdt*dtdm
909 jac(1, &
910 & 2) = dfdt*dtdp + ts*(cvs(lp)*(isentrope_n(lp) - 1)/(ps + isentrope_b(lp)) - cvs(vp)*(isentrope_n(vp) - 1)/(ps &
911 & + isentrope_b(vp)))
912 jac(2, &
913 & 1) = qvs(vp) - qvs(lp) + (cvs(vp)*isentrope_n(vp) - cvs(lp)*isentrope_n(lp))/(ml*(cvs(lp)*(isentrope_n(lp) - 1) &
914 & /(ps + isentrope_b(lp)) - cvs(vp)*(isentrope_n(vp) - 1)/(ps + isentrope_b(vp))) + mt*cvs(vp)*(isentrope_n(vp) &
915 & - 1)/(ps + isentrope_b(vp)) + mcvgp) - (ml*(cvs(vp)*isentrope_n(vp) - cvs(lp)*isentrope_n(lp)) - mt*cvs(vp) &
916 & *isentrope_n(vp) - mcpd)*(cvs(lp)*(isentrope_n(lp) - 1)/(ps + isentrope_b(lp)) - cvs(vp)*(isentrope_n(vp) - 1) &
917 & /(ps + isentrope_b(vp)))/((ml*(cvs(lp)*(isentrope_n(lp) - 1)/(ps + isentrope_b(lp)) - cvs(vp)*(isentrope_n(vp) &
918 & - 1)/(ps + isentrope_b(vp))) + mt*cvs(vp)*(isentrope_n(vp) - 1)/(ps + isentrope_b(vp)) + mcvgp)**2)
919 jac(2, &
920 & 2) = 1 + (ml*(cvs(vp)*isentrope_n(vp) - cvs(lp)*isentrope_n(lp)) - mt*cvs(vp)*isentrope_n(vp) - mcpd) &
921 & *(ml*(cvs(lp)*(isentrope_n(lp) - 1)/(ps + isentrope_b(lp))**2 - cvs(vp)*(isentrope_n(vp) - 1)/(ps &
922 & + isentrope_b(vp))**2) + mt*cvs(vp)*(isentrope_n(vp) - 1)/(ps + isentrope_b(vp))**2 + mcvgp2)/(ml*(cvs(lp) &
923 & *(isentrope_n(lp) - 1)/(ps + isentrope_b(lp)) - cvs(vp)*(isentrope_n(vp) - 1)/(ps + isentrope_b(vp))) &
924 & + mt*cvs(vp)*(isentrope_n(vp) - 1)/(ps + isentrope_b(vp)) + mcvgp)**2
925
926 detj = jac(1, 1)*jac(2, 2) - jac(1, 2)*jac(2, 1)
927 ! singular Jacobian: no usable Newton direction, accept the current (best) state
928 if (detj == 0.0_wp) exit
929
930 invjac(1, 1) = jac(2, 2)/detj
931 invjac(1, 2) = -jac(1, 2)/detj
932 invjac(2, 1) = -jac(2, 1)/detj
933 invjac(2, 2) = jac(1, 1)/detj
934
935 deltamp(1) = -(invjac(1, 1)*r2d(1) + invjac(1, 2)*r2d(2))
936 deltamp(2) = -(invjac(2, 1)*r2d(1) + invjac(2, 2)*r2d(2))
937
938 ! backtracking line search: halve the step until the residual decreases, keeping the state
939 ! physical (0 <= ml <= mT, pS above the stiffened-gas floor)
940 lambda = 1.0_wp
941 do ls = 1, ptg_ls_max
942 ml_try = min(max(ml + lambda*deltamp(1), 0.0_wp), mt)
943 ps_try = max(ps + lambda*deltamp(2), pmin)
944 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)
945 resnorm_try = sqrt(r2d_try(1)**2 + r2d_try(2)**2)
946 if ((resnorm_try < resnorm) .or. (ls == ptg_ls_max)) exit
947 lambda = 0.5_wp*lambda
948 end do
949
950 ! accept the trial state (TS, mCP, mQ, mCVGP, mCVGP2, mCPD already set to it by the last call)
951 ml = ml_try; ps = ps_try; r2d = r2d_try; resnorm = resnorm_try
952 end do
953
954 ! commit the reacting masses (mT conserved) and set the common temperature
955 q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) = ml
956 q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) = mt - ml
957
958 ts = (rhoe + ps - mq)/mcp
959
960 end subroutine s_infinite_ptg_relaxation_k
961
962 !> Correct the partial densities of the reacting fluids in case one of them is negative but their sum is positive. Inert phases
963 !! are not corrected at this moment
964 subroutine s_correct_partial_densities(MCT, q_cons_vf, rM, j, k, l)
965
966
967# 438 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
968#ifdef _CRAYFTN
969# 438 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
970#if MFC_OpenACC
971# 438 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
972!$acc routine seq
973# 438 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
974#elif MFC_OpenMP
975# 438 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
976
977# 438 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
978
979# 438 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
980!$omp declare target device_type(any)
981# 438 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
982#else
983# 438 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
984!DIR$ NOINLINE s_correct_partial_densities
985# 438 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
986#endif
987# 438 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
988#elif MFC_OpenACC
989# 438 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
990!$acc routine seq
991# 438 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
992#elif MFC_OpenMP
993# 438 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
994
995# 438 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
996
997# 438 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
998!$omp declare target device_type(any)
999# 438 "/home/runner/work/MFC/MFC/src/common/m_phase_change.fpp"
1000#endif
1001
1002 !> @name variables for the correction of the reacting partial densities
1003 !> @{
1004 real(wp), intent(out) :: mct
1005 type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
1006 real(wp), intent(inout) :: rm
1007 integer, intent(in) :: j, k, l
1008 !> @}
1009 if (rm < 0.0_wp) then
1010 if ((q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, &
1011 & l) >= -1.0_wp*mixm) .and. (q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) >= -1.0_wp*mixm)) then
1012 q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) = 0.0_wp
1013
1014 q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) = 0.0_wp
1015
1016 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)
1017 end if
1018 end if
1019
1020 ! TODO: Consider partitioning partial densities instead of absolute-value correction
1021 mct = 2*mixm
1022
1023 ! correcting the partial densities of the reacting fluids. What to do for the nonreacting ones?
1024 if (q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) < 0.0_wp) then
1025 q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) = mct*rm
1026
1027 q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) = (1.0_wp - mct)*rm
1028 else if (q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) < 0.0_wp) then
1029 q_cons_vf(lp + eqn_idx%cont%beg - 1)%sf(j, k, l) = (1.0_wp - mct)*rm
1030
1031 q_cons_vf(vp + eqn_idx%cont%beg - 1)%sf(j, k, l) = mct*rm
1032 end if
1033
1034 end subroutine s_correct_partial_densities
1035
1036 !> Finalize the phase change module
1039 end subroutine s_finalize_relaxation_solver_module
1040
1041end module m_phase_change
q_cons_vf
type(scalar_field), dimension(sys_size), intent(inout) q_cons_vf
Definition m_phase_change.fpp.f90:1005
rm
real(wp), intent(inout) rm
Definition m_phase_change.fpp.f90:1006
k
integer, intent(in) k
Definition m_phase_change.fpp.f90:1007
mct
real(wp), intent(out) mct
Definition m_phase_change.fpp.f90:1004
j
integer, intent(in) j
Definition m_phase_change.fpp.f90:1007
l
integer, intent(in) l
Definition m_phase_change.fpp.f90:1007
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:148
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_global_parameters
Global parameters for the post-process: domain geometry, equation of state, and output database setti...
Definition m_global_parameters.fpp.f90:19
m_helper_basic
Basic floating-point utilities: approximate equality, default detection, and coordinate bounds.
Definition m_helper_basic.fpp.f90:335
m_mpi_proxy
MPI gather and scatter operations for distributing post-process grid and flow-variable data.
Definition m_mpi_proxy.fpp.f90:7
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:1038
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:397
m_phase_change::vp
integer, parameter vp
index for the vapor phase of the reacting fluid
Definition m_phase_change.fpp.f90:369
m_phase_change::lp
integer, parameter lp
index for the liquid phase of the reacting fluid
Definition m_phase_change.fpp.f90:368
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:599
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:376
m_phase_change::tcr
real(wp), parameter tcr
Critical temperature of water [K].
Definition m_phase_change.fpp.f90:365
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:756
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:366
m_phase_change::max_iter
integer, parameter max_iter
max Newton iterations before accepting the last iterate
Definition m_phase_change.fpp.f90:363
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:392
m_phase_change::mixm
real(wp), parameter mixm
Mixture mass fraction threshold for triggering phase change.
Definition m_phase_change.fpp.f90:367
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:829
m_phase_change::pcr
real(wp), parameter pcr
Critical pressure of water [Pa].
Definition m_phase_change.fpp.f90:364
m_variables_conversion
Conservative-to-primitive variable conversion, mixture property evaluation, and pressure computation.
Definition m_variables_conversion.fpp.f90:346
m_variables_conversion::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_variables_conversion.fpp.f90:3260
m_derived_types::scalar_field
Derived type annexing a scalar field (SF).
Definition m_derived_types.fpp.f90:349