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