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