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