MFC
Exascale flow solver
Loading...
Searching...
No Matches
m_variables_conversion.fpp.f90
Go to the documentation of this file.
1# 1 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2!>
3!! @file
4!! @brief Contains module m_variables_conversion
5
6# 1 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 1
7# 1 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 1
8# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
9# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
10# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
11# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
12# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
13# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
14
15# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
16# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
17# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
18
19# 15 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
20# 16 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
21# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
22
23# 24 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
24
25# 53 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
26
27# 65 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
28
29# 75 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
30
31# 105 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
32
33# 117 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
34
35# 127 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
36
37# 174 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
38! New line at end of file is required for FYPP
39# 2 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
40# 1 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 1
41# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
42# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
43# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
44# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
45# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
46# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
47
48# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
49# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
50# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
51
52# 15 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
53# 16 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
54# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
55
56# 24 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
57
58# 53 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
59
60# 65 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
61
62# 75 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
63
64# 105 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
65
66# 117 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
67
68# 127 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
69
70# 174 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
71! New line at end of file is required for FYPP
72# 2 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 2
73
74# 4 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
75# 5 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
76# 6 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
77# 7 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
78# 8 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
79
80# 20 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
81
82# 43 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
83
84# 48 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
85
86# 53 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
87
88# 58 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
89
90# 63 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
91
92# 68 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
93
94# 81 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
95
96# 86 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
97
98# 91 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
99
100# 96 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
101
102# 101 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
103
104# 106 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
105
106# 111 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
107
108# 116 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
109
110# 121 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
111
112# 126 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
113
114# 156 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
115
116# 197 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
117
118# 211 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
119
120# 236 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
121
122# 247 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
123
124# 249 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
125# 260 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
126
127# 310 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
128
129# 320 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
130
131# 330 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
132
133# 339 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
134
135# 356 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
136
137# 366 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
138
139# 373 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
140
141# 379 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
142
143# 385 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
144
145# 391 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
146
147# 397 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
148
149# 403 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
150! New line at end of file is required for FYPP
151# 3 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
152# 1 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 1
153# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
154# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
155# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
156# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
157# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
158# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
159
160# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
161# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
162# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
163
164# 15 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
165# 16 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
166# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
167
168# 24 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
169
170# 53 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
171
172# 65 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
173
174# 75 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
175
176# 105 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
177
178# 117 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
179
180# 127 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
181
182# 174 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
183! New line at end of file is required for FYPP
184# 2 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 2
185
186# 7 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
187
188# 17 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
189
190# 22 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
191
192# 27 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
193
194# 32 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
195
196# 37 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
197
198# 42 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
199
200# 47 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
201
202# 52 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
203
204# 57 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
205
206# 62 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
207
208# 73 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
209
210# 78 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
211
212# 83 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
213
214# 88 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
215
216# 103 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
217
218# 131 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
219
220# 160 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
221
222# 175 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
223
224# 193 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
225
226# 215 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
227
228# 244 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
229
230# 259 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
231
232# 269 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
233
234# 278 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
235
236# 294 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
237
238# 304 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
239
240# 311 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
241! New line at end of file is required for FYPP
242# 4 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
243
244! GPU parallel region (scalar reductions, maxval/minval)
245# 23 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
246
247! GPU parallel loop over threads (most common GPU macro)
248# 43 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
249
250! Required closing for GPU_PARALLEL_LOOP
251# 55 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
252
253! Mark routine for device compilation
254# 112 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
255
256! Declare device-resident data
257# 130 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
258
259! Inner loop within a GPU parallel region
260# 145 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
261
262! Scoped GPU data region
263# 164 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
264
265! Host code with device pointers (for MPI with GPU buffers)
266# 193 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
267
268! Allocate device memory (unscoped)
269# 207 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
270
271! Free device memory
272# 219 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
273
274! Atomic operation on device
275# 231 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
276
277! End atomic capture block
278# 242 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
279
280! Copy data between host and device
281# 254 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
282
283! Synchronization barrier
284# 266 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
285
286! Import GPU library module (openacc or omp_lib)
287# 275 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
288
289! Emit code only for AMD compiler
290# 282 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
291
292! Emit code for non-Cray compilers
293# 289 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
294
295! Emit code only for Cray compiler
296# 296 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
297
298! Emit code for non-NVIDIA compilers
299# 303 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
300
301# 305 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
302# 306 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
303! New line at end of file is required for FYPP
304# 2 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 2
305
306# 14 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
307
308! Caution: This macro requires the use of a binding script to set CUDA_VISIBLE_DEVICES, such that we have one GPU device per MPI
309! rank. That's because for both cudaMemAdvise (preferred location) and cudaMemPrefetchAsync we use location = device_id = 0. For an
310! example see misc/nvidia_uvm/bind.sh.
311# 52 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
312
313! Allocate and create GPU device memory
314# 72 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
315
316! Free GPU device memory and deallocate
317# 80 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
318
319! Cray-specific GPU pointer setup for vector fields
320# 104 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
321
322! Cray-specific GPU pointer setup for scalar fields
323# 120 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
324
325! Cray-specific GPU pointer setup for acoustic source spatials
326# 145 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
327
328# 151 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
329
330# 158 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
331! New line at end of file is required for FYPP
332# 6 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp" 2
333# 1 "/home/runner/work/MFC/MFC/src/common/include/case.fpp" 1
334! This file exists so that Fypp can be run without generating case.fpp files for
335! each target. This is useful when generating documentation, for example. This
336! should also let MFC be built with CMake directly, without invoking mfc.sh.
337
338! For pre-process.
339# 8 "/home/runner/work/MFC/MFC/src/common/include/case.fpp"
340
341! For moving immersed boundaries in simulation
342# 12 "/home/runner/work/MFC/MFC/src/common/include/case.fpp"
343# 7 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp" 2
344
345!> @brief Conservative-to-primitive variable conversion, mixture property evaluation, and pressure computation
347
350 use m_mpi_proxy
352 use m_helper
355 use m_thermochem, only: num_species, get_temperature, get_pressure, gas_constant, get_mixture_molecular_weight, &
356 & get_mixture_energy_mass
357
358 implicit none
359
360 private
370 & s_finalize_variables_conversion_module, gammas, isentrope_n, pi_infs, isentrope_b, cvs, qvs, qvps
371
372 real(wp), allocatable, dimension(:) :: gs_vc
373 integer, allocatable, dimension(:) :: bubrs_vc
374 real(wp), allocatable, dimension(:,:) :: res_vc
375
376# 38 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
377#if defined(MFC_OpenACC)
378# 38 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
379!$acc declare create(bubrs_vc, Gs_vc, Res_vc)
380# 38 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
381#elif defined(MFC_OpenMP)
382# 38 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
383!$omp declare target (bubrs_vc, Gs_vc, Res_vc)
384# 38 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
385#endif
386
387 integer :: is1b, is2b, is3b, is1e, is2e, is3e
388
389# 41 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
390#if defined(MFC_OpenACC)
391# 41 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
392!$acc declare create(is1b, is2b, is3b, is1e, is2e, is3e)
393# 41 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
394#elif defined(MFC_OpenMP)
395# 41 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
396!$omp declare target (is1b, is2b, is3b, is1e, is2e, is3e)
397# 41 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
398#endif
399
400 logical :: enforce_density_floor_vc = .false.
401 logical :: preserve_qbmm_number_vc = .false.
403
404# 46 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
405#if defined(MFC_OpenACC)
406# 46 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
407!$acc declare create(enforce_density_floor_vc, preserve_qbmm_number_vc, lagrange_beta_index_vc)
408# 46 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
409#elif defined(MFC_OpenMP)
410# 46 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
411!$omp declare target (enforce_density_floor_vc, preserve_qbmm_number_vc, lagrange_beta_index_vc)
412# 46 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
413#endif
414
415 real(wp), allocatable, dimension(:,:,:), public :: rho_sf !< Scalar density function
416 real(wp), allocatable, dimension(:,:,:), public :: gamma_sf !< Scalar sp. heat ratio function
417 real(wp), allocatable, dimension(:,:,:), public :: pi_inf_sf !< Scalar liquid stiffness function
418
419contains
420
421 !> Dispatch to the s_convert_mixture_to_mixture_variables and s_convert_species_to_mixture_variables subroutines. Replaces a
422 !! procedure pointer.
423 subroutine s_convert_to_mixture_variables(q_vf, i, j, k, rho, gamma, pi_inf, qv, Re_K, G_K, G)
424
425 type(scalar_field), dimension(sys_size), intent(in) :: q_vf
426 integer, intent(in) :: i, j, k
427 real(wp), intent(out), target :: rho, gamma, pi_inf, qv
428 real(wp), optional, dimension(2), intent(out) :: re_k
429 real(wp), optional, intent(out) :: g_k
430 real(wp), optional, dimension(num_fluids), intent(in) :: g
431
432 if (model_eqns == model_eqns_gamma_law) then ! Gamma/pi_inf model
433 call s_convert_mixture_to_mixture_variables(q_vf, i, j, k, rho, gamma, pi_inf, qv)
434 else ! Volume fraction model
435 call s_convert_species_to_mixture_variables(q_vf, i, j, k, rho, gamma, pi_inf, qv, re_k, g_k, g)
436 end if
437
438 end subroutine s_convert_to_mixture_variables
439
440 !> Compute the pressure from the appropriate equation of state
441 subroutine s_compute_pressure(energy, alf, dyn_p, pi_inf, gamma, rho, qv, rhoYks, pres, T, E_e_in, pres_mag)
442
443
444# 76 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
445#ifdef _CRAYFTN
446# 76 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
447#if MFC_OpenACC
448# 76 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
449!$acc routine seq
450# 76 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
451#elif MFC_OpenMP
452# 76 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
453
454# 76 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
455
456# 76 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
457!$omp declare target device_type(any)
458# 76 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
459#else
460# 76 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
461!DIR$ NOINLINE s_compute_pressure
462# 76 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
463#endif
464# 76 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
465#elif MFC_OpenACC
466# 76 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
467!$acc routine seq
468# 76 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
469#elif MFC_OpenMP
470# 76 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
471
472# 76 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
473
474# 76 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
475!$omp declare target device_type(any)
476# 76 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
477#endif
478
479 real(stp), intent(in) :: energy, alf
480 real(wp), intent(in) :: dyn_p
481 real(wp), intent(in) :: pi_inf, gamma, rho, qv
482 real(wp), intent(out) :: pres
483 real(wp), intent(inout) :: t
484 real(wp), intent(in), optional :: e_e_in, pres_mag
485
486 ! Chemistry
487 real(wp), dimension(1:num_species), intent(in) :: rhoyks
488 real(wp), dimension(1:num_species) :: y_rs
489 real(wp) :: e_int
490 real(wp) :: e_per_kg, pdyn_per_kg
491 real(wp) :: t_guess
492# 92 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
493 ! What is the internal energy? The magnetic, elastic and kinetic parts are model
494 ! bookkeeping rather than equation of state, so they come off here and the inversion runs once.
495 if (mhd) then
496 ! MHD: the magnetic energy is not an equation-of-state term
497 e_int = energy - dyn_p - pres_mag
498 else if (bubbles_euler .neqv. .true.) then
499 ! Gamma/pi_inf model or five-equation model (Allaire et al. JCP 2002)
500 e_int = energy - dyn_p
501 else
502 ! Bubble-augmented; qv comes off before the division rather than being scaled by it
503 e_int = (energy - dyn_p - qv)/(1._wp - alf) + qv
504 end if
505
506 if (hypoelasticity .and. present(e_e_in)) then
507 ! Subtract elastic strain energy before computing pressure (hypoelastic model)
508 e_int = energy - dyn_p - e_e_in
509 end if
510
511 pres = f_pressure(e_int, gamma, pi_inf, qv)
512# 122 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
513
514 end subroutine s_compute_pressure
515
516 !> Convert mixture variables to density, gamma, pi_inf, and qv for the gamma/pi_inf model. Given conservative or primitive
517 !! variables, transfers the density, specific heat ratio function and the liquid stiffness function from q_vf to rho, gamma and
518 !! pi_inf.
519 subroutine s_convert_mixture_to_mixture_variables(q_vf, i, j, k, rho, gamma, pi_inf, qv)
520
521 type(scalar_field), dimension(sys_size), intent(in) :: q_vf
522 integer, intent(in) :: i, j, k
523 real(wp), intent(out), target :: rho
524 real(wp), intent(out), target :: gamma
525 real(wp), intent(out), target :: pi_inf
526 real(wp), intent(out), target :: qv
527
528 ! Transferring the density, the specific heat ratio function and the liquid stiffness function, respectively
529
530 rho = q_vf(1)%sf(i, j, k)
531 gamma = q_vf(eqn_idx%gamma)%sf(i, j, k)
532 pi_inf = q_vf(eqn_idx%pi_inf)%sf(i, j, k)
533 qv = 0._wp ! keep this value nil for now. For future adjustment
534
535 ! Store derived mixture fields when requested during module initialization.
536 if (allocated(rho_sf)) then
537 rho_sf(i, j, k) = rho
538 gamma_sf(i, j, k) = gamma
539 pi_inf_sf(i, j, k) = pi_inf
540 end if
541
543
544 !> Convert species volume fractions and partial densities to mixture density, gamma, pi_inf, and qv. Given conservative or
545 !! primitive variables, computes the density, the specific heat ratio function and the liquid stiffness function from q_vf and
546 !! stores the results into rho, gamma and pi_inf.
547 subroutine s_convert_species_to_mixture_variables(q_vf, k, l, r, rho, gamma, pi_inf, qv, Re_K, G_K, G)
548
549 type(scalar_field), dimension(sys_size), intent(in) :: q_vf
550 integer, intent(in) :: k, l, r
551 real(wp), intent(out), target :: rho
552 real(wp), intent(out), target :: gamma
553 real(wp), intent(out), target :: pi_inf
554 real(wp), intent(out), target :: qv
555 real(wp), optional, dimension(2), intent(out) :: re_k
556 real(wp), optional, intent(out) :: g_k
557 real(wp), dimension(num_fluids) :: alpha_rho_k, alpha_k
558 real(wp), optional, dimension(num_fluids), intent(in) :: g
559 integer :: i, j !< Generic loop iterator
560 ! Computing the density, the specific heat ratio function and the liquid stiffness function, respectively
561
562 call s_compute_species_fraction(q_vf, k, l, r, alpha_rho_k, alpha_k)
563
564 ! Use the same scalar kernel on host and device so mixture semantics do not depend on the executable or accelerator backend.
565 ! Absent optional dummies forward as absent, so the optional arguments need no dispatch here.
566 call s_convert_species_to_mixture_variables_kernel(rho, gamma, pi_inf, qv, alpha_k, alpha_rho_k, re_k, g_k, g)
567
568 ! Store derived mixture fields when requested during module initialization.
569 if (allocated(rho_sf)) then
570 rho_sf(k, l, r) = rho
571 gamma_sf(k, l, r) = gamma
572 pi_inf_sf(k, l, r) = pi_inf
573 end if
574
576
577 !> Host- and device-callable conversion kernel for species and mixture variables.
578 subroutine s_convert_species_to_mixture_variables_kernel(rho_K, gamma_K, pi_inf_K, qv_K, alpha_K, alpha_rho_K, Re_K, G_K, G)
579
580
581# 189 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
582#ifdef _CRAYFTN
583# 189 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
584#if MFC_OpenACC
585# 189 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
586!$acc routine seq
587# 189 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
588#elif MFC_OpenMP
589# 189 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
590
591# 189 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
592
593# 189 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
594!$omp declare target device_type(any)
595# 189 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
596#else
597# 189 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
598!DIR$ NOINLINE s_convert_species_to_mixture_variables_kernel
599# 189 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
600#endif
601# 189 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
602#elif MFC_OpenACC
603# 189 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
604!$acc routine seq
605# 189 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
606#elif MFC_OpenMP
607# 189 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
608
609# 189 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
610
611# 189 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
612!$omp declare target device_type(any)
613# 189 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
614#endif
615
616 real(wp), intent(out) :: rho_k, gamma_k, pi_inf_k, qv_k
617# 196 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
618 real(wp), dimension(num_fluids), intent(inout) :: alpha_rho_k, alpha_k
619 real(wp), optional, dimension(num_fluids), intent(in) :: g
620# 199 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
621 real(wp), optional, dimension(2), intent(out) :: re_k
622 real(wp), optional, intent(out) :: g_k
623 real(wp) :: alpha_k_sum
624 integer :: i, j !< Generic loop iterators
625
626 rho_k = 0._wp
627 gamma_k = 0._wp
628 pi_inf_k = 0._wp
629 qv_k = 0._wp
630 if (present(re_k)) re_k = dflt_real
631 if (present(g_k)) g_k = 0._wp
632
633 ! Constrain partial densities and volume fractions within physical bounds
634 if (mpp_lim) then
635 alpha_k_sum = 0._wp
636 do i = 1, num_fluids
637 alpha_rho_k(i) = max(0._wp, alpha_rho_k(i))
638 alpha_k(i) = min(max(0._wp, alpha_k(i)), 1._wp)
639 alpha_k_sum = alpha_k_sum + alpha_k(i)
640 end do
641 alpha_k = alpha_k/max(alpha_k_sum, sgm_eps)
642 end if
643 call s_compute_mixture_coefficients(alpha_rho_k, alpha_k, rho_k, gamma_k, pi_inf_k, qv_k)
644
645 if (present(g_k)) then
646 g_k = 0._wp
647 do i = 1, num_fluids
648 ! TODO: change to use Gs_vc directly here? TODO: Make this change as well for GPUs
649 g_k = g_k + alpha_k(i)*g(i)
650 end do
651 g_k = max(0._wp, g_k)
652 end if
653
654 if (viscous .and. present(re_k)) then
655 do i = 1, 2
656 re_k(i) = dflt_real
657
658 if (re_size(i) > 0) re_k(i) = 0._wp
659
660 do j = 1, re_size(i)
661 re_k(i) = alpha_k(re_idx(i, j))/res_vc(i, j) + re_k(i)
662 end do
663
664 re_k(i) = 1._wp/max(re_k(i), sgm_eps)
665 end do
666 end if
667
669
670 !> Initialize the variables conversion module.
671 impure subroutine s_initialize_variables_conversion_module(store_mixture_fields, enforce_density_floor, preserve_qbmm_number, &
672 & lagrange_beta_index)
673
674 integer :: i, j
675 logical, optional, intent(in) :: store_mixture_fields
676 logical, optional, intent(in) :: enforce_density_floor, preserve_qbmm_number
677 integer, optional, intent(in) :: lagrange_beta_index
678 logical :: allocate_mixture_fields
679 logical :: state_dependent !< Whether this case's fluids need a density-dependent EOS
680
681 allocate_mixture_fields = .false.
682 if (present(store_mixture_fields)) allocate_mixture_fields = store_mixture_fields
684 if (present(enforce_density_floor)) enforce_density_floor_vc = enforce_density_floor
686 if (present(preserve_qbmm_number)) preserve_qbmm_number_vc = preserve_qbmm_number
688 if (present(lagrange_beta_index)) lagrange_beta_index_vc = lagrange_beta_index
689
690
691# 268 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
692#if defined(MFC_OpenACC)
693# 268 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
694!$acc enter data copyin(is1b, is1e, is2b, is2e, is3b, is3e)
695# 268 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
696#elif defined(MFC_OpenMP)
697# 268 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
698!$omp target enter data map(to:is1b, is1e, is2b, is2e, is3b, is3e)
699# 268 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
700#endif
701
702# 269 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
703#if defined(MFC_OpenACC)
704# 269 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
705!$acc update device(enforce_density_floor_vc, preserve_qbmm_number_vc, lagrange_beta_index_vc)
706# 269 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
707#elif defined(MFC_OpenMP)
708# 269 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
709!$omp target update to(enforce_density_floor_vc, preserve_qbmm_number_vc, lagrange_beta_index_vc)
710# 269 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
711#endif
712
713#ifdef MFC_DEBUG
714# 271 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
715 block
716# 271 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
717 use iso_fortran_env, only: output_unit
718# 271 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
719
720# 271 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
721 print *, 'm_variables_conversion.fpp:271: ', '@:ALLOCATE(gammas (1:num_fluids))'
722# 271 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
723
724# 271 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
725 call flush (output_unit)
726# 271 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
727 end block
728# 271 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
729#endif
730# 271 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
731 allocate (gammas(1:num_fluids))
732# 271 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
733
734# 271 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
735
736# 271 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
737#if defined(MFC_OpenACC)
738# 271 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
739!$acc enter data create(gammas)
740# 271 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
741#elif defined(MFC_OpenMP)
742# 271 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
743!$omp target enter data map(always,alloc:gammas)
744# 271 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
745#endif
746#ifdef MFC_DEBUG
747# 272 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
748 block
749# 272 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
750 use iso_fortran_env, only: output_unit
751# 272 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
752
753# 272 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
754 print *, 'm_variables_conversion.fpp:272: ', '@:ALLOCATE(eoss (1:num_fluids))'
755# 272 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
756
757# 272 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
758 call flush (output_unit)
759# 272 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
760 end block
761# 272 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
762#endif
763# 272 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
764 allocate (eoss(1:num_fluids))
765# 272 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
766
767# 272 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
768
769# 272 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
770#if defined(MFC_OpenACC)
771# 272 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
772!$acc enter data create(eoss)
773# 272 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
774#elif defined(MFC_OpenMP)
775# 272 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
776!$omp target enter data map(always,alloc:eoss)
777# 272 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
778#endif
779#ifdef MFC_DEBUG
780# 273 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
781 block
782# 273 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
783 use iso_fortran_env, only: output_unit
784# 273 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
785
786# 273 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
787 print *, 'm_variables_conversion.fpp:273: ', '@:ALLOCATE(isentrope_n (1:num_fluids))'
788# 273 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
789
790# 273 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
791 call flush (output_unit)
792# 273 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
793 end block
794# 273 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
795#endif
796# 273 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
797 allocate (isentrope_n(1:num_fluids))
798# 273 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
799
800# 273 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
801
802# 273 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
803#if defined(MFC_OpenACC)
804# 273 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
805!$acc enter data create(isentrope_n)
806# 273 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
807#elif defined(MFC_OpenMP)
808# 273 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
809!$omp target enter data map(always,alloc:isentrope_n)
810# 273 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
811#endif
812#ifdef MFC_DEBUG
813# 274 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
814 block
815# 274 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
816 use iso_fortran_env, only: output_unit
817# 274 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
818
819# 274 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
820 print *, 'm_variables_conversion.fpp:274: ', '@:ALLOCATE(pi_infs(1:num_fluids))'
821# 274 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
822
823# 274 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
824 call flush (output_unit)
825# 274 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
826 end block
827# 274 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
828#endif
829# 274 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
830 allocate (pi_infs(1:num_fluids))
831# 274 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
832
833# 274 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
834
835# 274 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
836#if defined(MFC_OpenACC)
837# 274 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
838!$acc enter data create(pi_infs)
839# 274 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
840#elif defined(MFC_OpenMP)
841# 274 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
842!$omp target enter data map(always,alloc:pi_infs)
843# 274 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
844#endif
845#ifdef MFC_DEBUG
846# 275 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
847 block
848# 275 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
849 use iso_fortran_env, only: output_unit
850# 275 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
851
852# 275 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
853 print *, 'm_variables_conversion.fpp:275: ', '@:ALLOCATE(isentrope_B(1:num_fluids))'
854# 275 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
855
856# 275 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
857 call flush (output_unit)
858# 275 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
859 end block
860# 275 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
861#endif
862# 275 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
863 allocate (isentrope_b(1:num_fluids))
864# 275 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
865
866# 275 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
867
868# 275 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
869#if defined(MFC_OpenACC)
870# 275 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
871!$acc enter data create(isentrope_B)
872# 275 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
873#elif defined(MFC_OpenMP)
874# 275 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
875!$omp target enter data map(always,alloc:isentrope_B)
876# 275 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
877#endif
878#ifdef MFC_DEBUG
879# 276 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
880 block
881# 276 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
882 use iso_fortran_env, only: output_unit
883# 276 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
884
885# 276 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
886 print *, 'm_variables_conversion.fpp:276: ', '@:ALLOCATE(cvs (1:num_fluids))'
887# 276 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
888
889# 276 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
890 call flush (output_unit)
891# 276 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
892 end block
893# 276 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
894#endif
895# 276 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
896 allocate (cvs(1:num_fluids))
897# 276 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
898
899# 276 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
900
901# 276 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
902#if defined(MFC_OpenACC)
903# 276 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
904!$acc enter data create(cvs)
905# 276 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
906#elif defined(MFC_OpenMP)
907# 276 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
908!$omp target enter data map(always,alloc:cvs)
909# 276 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
910#endif
911#ifdef MFC_DEBUG
912# 277 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
913 block
914# 277 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
915 use iso_fortran_env, only: output_unit
916# 277 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
917
918# 277 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
919 print *, 'm_variables_conversion.fpp:277: ', '@:ALLOCATE(qvs (1:num_fluids))'
920# 277 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
921
922# 277 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
923 call flush (output_unit)
924# 277 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
925 end block
926# 277 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
927#endif
928# 277 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
929 allocate (qvs(1:num_fluids))
930# 277 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
931
932# 277 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
933
934# 277 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
935#if defined(MFC_OpenACC)
936# 277 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
937!$acc enter data create(qvs)
938# 277 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
939#elif defined(MFC_OpenMP)
940# 277 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
941!$omp target enter data map(always,alloc:qvs)
942# 277 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
943#endif
944#ifdef MFC_DEBUG
945# 278 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
946 block
947# 278 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
948 use iso_fortran_env, only: output_unit
949# 278 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
950
951# 278 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
952 print *, 'm_variables_conversion.fpp:278: ', '@:ALLOCATE(qvps (1:num_fluids))'
953# 278 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
954
955# 278 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
956 call flush (output_unit)
957# 278 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
958 end block
959# 278 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
960#endif
961# 278 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
962 allocate (qvps(1:num_fluids))
963# 278 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
964
965# 278 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
966
967# 278 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
968#if defined(MFC_OpenACC)
969# 278 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
970!$acc enter data create(qvps)
971# 278 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
972#elif defined(MFC_OpenMP)
973# 278 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
974!$omp target enter data map(always,alloc:qvps)
975# 278 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
976#endif
977#ifdef MFC_DEBUG
978# 279 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
979 block
980# 279 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
981 use iso_fortran_env, only: output_unit
982# 279 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
983
984# 279 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
985 print *, 'm_variables_conversion.fpp:279: ', '@:ALLOCATE(Gs_vc (1:num_fluids))'
986# 279 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
987
988# 279 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
989 call flush (output_unit)
990# 279 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
991 end block
992# 279 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
993#endif
994# 279 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
995 allocate (gs_vc(1:num_fluids))
996# 279 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
997
998# 279 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
999
1000# 279 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1001#if defined(MFC_OpenACC)
1002# 279 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1003!$acc enter data create(Gs_vc)
1004# 279 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1005#elif defined(MFC_OpenMP)
1006# 279 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1007!$omp target enter data map(always,alloc:Gs_vc)
1008# 279 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1009#endif
1010#ifdef MFC_DEBUG
1011# 280 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1012 block
1013# 280 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1014 use iso_fortran_env, only: output_unit
1015# 280 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1016
1017# 280 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1018 print *, 'm_variables_conversion.fpp:280: ', '@:ALLOCATE(fluid_k_therm(1:num_fluids))'
1019# 280 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1020
1021# 280 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1022 call flush (output_unit)
1023# 280 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1024 end block
1025# 280 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1026#endif
1027# 280 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1028 allocate (fluid_k_therm(1:num_fluids))
1029# 280 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1030
1031# 280 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1032
1033# 280 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1034#if defined(MFC_OpenACC)
1035# 280 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1036!$acc enter data create(fluid_k_therm)
1037# 280 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1038#elif defined(MFC_OpenMP)
1039# 280 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1040!$omp target enter data map(always,alloc:fluid_k_therm)
1041# 280 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1042#endif
1043
1044 state_dependent = .false.
1045 do i = 1, num_fluids
1046 gammas(i) = fluid_pp(i)%gamma
1047 isentrope_n(i) = f_isentrope_exponent(gammas(i))
1048
1049 ! Each EOS supplies its own coefficients. Resolved once here, not per cell: a branch in the mixture loop costs
1050 ! registers in the Riemann kernels. An EOS whose coefficients depend on state must move to per-cell evaluation.
1051 select case (fluid_pp(i)%eos)
1052 case (eos_ideal_gas)
1053 pi_infs(i) = 0._wp
1054 case default
1055 pi_infs(i) = fluid_pp(i)%pi_inf
1056 end select
1057 gs_vc(i) = fluid_pp(i)%G
1058 isentrope_b(i) = f_isentrope_pressure(pi_infs(i), gammas(i))
1059 cvs(i) = fluid_pp(i)%cv
1060 qvs(i) = fluid_pp(i)%qv
1061 qvps(i) = fluid_pp(i)%qvp
1062 fluid_k_therm(i) = fluid_pp(i)%k_therm
1063 eoss(i) = fluid_pp(i)%eos
1064 eos_coeffs(i)%c0 = fluid_pp(i)%mg_c0
1065 eos_coeffs(i)%s = fluid_pp(i)%mg_s
1066 eos_coeffs(i)%s2 = fluid_pp(i)%mg_s2
1067 eos_coeffs(i)%s3 = fluid_pp(i)%mg_s3
1068 ! Where a cubic Hugoniot fit turns over: mu(u_p) peaks where c0 = s2 u_p^2 + 2 s3 u_p^3, and past it
1069 ! no shock state exists, so the Newton below would wander. Solved once here, on the host.
1070 eos_coeffs(i)%mu_max = f_hugoniot_compression_limit(fluid_pp(i)%mg_c0, fluid_pp(i)%mg_s, fluid_pp(i)%mg_s2, &
1071 & fluid_pp(i)%mg_s3)
1072 eos_coeffs(i)%a = fluid_pp(i)%jwl_a
1073 eos_coeffs(i)%b = fluid_pp(i)%jwl_b
1074 eos_coeffs(i)%r1 = fluid_pp(i)%jwl_r1
1075 eos_coeffs(i)%r2 = fluid_pp(i)%jwl_r2
1076 eos_coeffs(i)%k0 = fluid_pp(i)%vinet_k0
1077 eos_coeffs(i)%k0p = fluid_pp(i)%vinet_k0p
1078 ! One reference state and Gruneisen closure for every family; the user-facing names keep their prefix.
1079 select case (fluid_pp(i)%eos)
1080 case (eos_mie_gruneisen)
1081 eos_coeffs(i)%rho0 = fluid_pp(i)%mg_rho0
1082 eos_coeffs(i)%t0 = fluid_pp(i)%mg_t0
1083 eos_coeffs(i)%gruneisen0 = fluid_pp(i)%mg_gruneisen
1084 eos_coeffs(i)%gruneisen_a = fluid_pp(i)%mg_gruneisen_a
1085 case (eos_jwl)
1086 eos_coeffs(i)%rho0 = fluid_pp(i)%jwl_rho0
1087 eos_coeffs(i)%t0 = fluid_pp(i)%jwl_t0
1088 eos_coeffs(i)%gruneisen0 = fluid_pp(i)%jwl_omega
1089 eos_coeffs(i)%gruneisen_a = 0._wp
1090 case default
1091 eos_coeffs(i)%rho0 = dflt_real
1092 eos_coeffs(i)%t0 = dflt_real
1093 eos_coeffs(i)%gruneisen0 = dflt_real
1094 eos_coeffs(i)%gruneisen_a = 0._wp
1095 case (eos_vinet)
1096 eos_coeffs(i)%rho0 = fluid_pp(i)%vinet_rho0
1097 eos_coeffs(i)%t0 = fluid_pp(i)%vinet_t0
1098 eos_coeffs(i)%gruneisen0 = fluid_pp(i)%vinet_gruneisen
1099 eos_coeffs(i)%gruneisen_a = fluid_pp(i)%vinet_gruneisen_a
1100 end select
1101 if (f_is_state_dependent(i)) state_dependent = .true.
1102 end do
1103# 347 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1104 any_state_dependent_eos = state_dependent
1105# 349 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1106
1107# 349 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1108#if defined(MFC_OpenACC)
1109# 349 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1110!$acc update device(gammas, isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps, Gs_vc, eoss, eos_coeffs, fluid_k_therm, heat_conduction)
1111# 349 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1112#elif defined(MFC_OpenMP)
1113# 349 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1114!$omp target update to(gammas, isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps, Gs_vc, eoss, eos_coeffs, fluid_k_therm, heat_conduction)
1115# 349 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1116#endif
1117# 351 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1118
1119# 351 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1120#if defined(MFC_OpenACC)
1121# 351 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1122!$acc update device(any_state_dependent_eos)
1123# 351 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1124#elif defined(MFC_OpenMP)
1125# 351 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1126!$omp target update to(any_state_dependent_eos)
1127# 351 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1128#endif
1129# 353 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1130
1131#ifdef MFC_DEBUG
1132# 354 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1133 block
1134# 354 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1135 use iso_fortran_env, only: output_unit
1136# 354 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1137
1138# 354 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1139 print *, 'm_variables_conversion.fpp:354: ', '@:ALLOCATE(Res_vc(1:2, 1:max(1, Re_size_max)))'
1140# 354 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1141
1142# 354 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1143 call flush (output_unit)
1144# 354 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1145 end block
1146# 354 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1147#endif
1148# 354 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1149 allocate (res_vc(1:2, 1:max(1, re_size_max)))
1150# 354 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1151
1152# 354 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1153
1154# 354 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1155#if defined(MFC_OpenACC)
1156# 354 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1157!$acc enter data create(Res_vc)
1158# 354 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1159#elif defined(MFC_OpenMP)
1160# 354 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1161!$omp target enter data map(always,alloc:Res_vc)
1162# 354 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1163#endif
1164 res_vc = dflt_real
1165 if (allocated(re_idx)) then
1166 do i = 1, 2
1167 do j = 1, re_size(i)
1168 res_vc(i, j) = fluid_pp(re_idx(i, j))%Re(i)
1169 end do
1170 end do
1171
1172# 362 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1173#if defined(MFC_OpenACC)
1174# 362 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1175!$acc update device(Re_idx)
1176# 362 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1177#elif defined(MFC_OpenMP)
1178# 362 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1179!$omp target update to(Re_idx)
1180# 362 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1181#endif
1182 end if
1183
1184# 364 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1185#if defined(MFC_OpenACC)
1186# 364 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1187!$acc update device(Res_vc, Re_size)
1188# 364 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1189#elif defined(MFC_OpenMP)
1190# 364 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1191!$omp target update to(Res_vc, Re_size)
1192# 364 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1193#endif
1194
1195 if (bubbles_euler) then
1196#ifdef MFC_DEBUG
1197# 367 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1198 block
1199# 367 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1200 use iso_fortran_env, only: output_unit
1201# 367 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1202
1203# 367 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1204 print *, 'm_variables_conversion.fpp:367: ', '@:ALLOCATE(bubrs_vc(1:nb))'
1205# 367 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1206
1207# 367 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1208 call flush (output_unit)
1209# 367 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1210 end block
1211# 367 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1212#endif
1213# 367 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1214 allocate (bubrs_vc(1:nb))
1215# 367 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1216
1217# 367 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1218
1219# 367 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1220#if defined(MFC_OpenACC)
1221# 367 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1222!$acc enter data create(bubrs_vc)
1223# 367 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1224#elif defined(MFC_OpenMP)
1225# 367 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1226!$omp target enter data map(always,alloc:bubrs_vc)
1227# 367 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1228#endif
1229 do i = 1, nb
1230 bubrs_vc(i) = qbmm_idx%rs(i)
1231 end do
1232
1233# 371 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1234#if defined(MFC_OpenACC)
1235# 371 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1236!$acc update device(bubrs_vc)
1237# 371 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1238#elif defined(MFC_OpenMP)
1239# 371 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1240!$omp target update to(bubrs_vc)
1241# 371 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1242#endif
1243 end if
1244
1245 if (allocate_mixture_fields) then
1246 ! Allocate derived mixture fields over the available grid storage.
1247 if (n > 0) then
1248 if (p > 0) then
1249 allocate (rho_sf(-buff_size:m + buff_size,-buff_size:n + buff_size,-buff_size:p + buff_size))
1250 allocate (gamma_sf(-buff_size:m + buff_size,-buff_size:n + buff_size,-buff_size:p + buff_size))
1251 allocate (pi_inf_sf(-buff_size:m + buff_size,-buff_size:n + buff_size,-buff_size:p + buff_size))
1252 else
1253 allocate (rho_sf(-buff_size:m + buff_size,-buff_size:n + buff_size,0:0))
1254 allocate (gamma_sf(-buff_size:m + buff_size,-buff_size:n + buff_size,0:0))
1255 allocate (pi_inf_sf(-buff_size:m + buff_size,-buff_size:n + buff_size,0:0))
1256 end if
1257 else
1258 allocate (rho_sf(-buff_size:m + buff_size,0:0,0:0))
1259 allocate (gamma_sf(-buff_size:m + buff_size,0:0,0:0))
1260 allocate (pi_inf_sf(-buff_size:m + buff_size,0:0,0:0))
1261 end if
1262 end if
1263
1265
1266 !> Initialize bubble mass-vapor values at quadrature nodes from the conserved moment statistics.
1267 subroutine s_initialize_mv(qK_cons_vf, mv)
1268
1269 type(scalar_field), dimension(sys_size), intent(in) :: qk_cons_vf
1270 real(stp), dimension(idwint(1)%beg:,idwint(2)%beg:,idwint(3)%beg:,1:,1:), intent(inout) :: mv
1271 integer :: i, j, k, l
1272 real(wp) :: mu, sig, nbub_sc
1273
1274 do l = idwint(3)%beg, idwint(3)%end
1275 do k = idwint(2)%beg, idwint(2)%end
1276 do j = idwint(1)%beg, idwint(1)%end
1277 nbub_sc = qk_cons_vf(eqn_idx%bub%beg)%sf(j, k, l)
1278
1279
1280# 408 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1281#if defined(MFC_OpenACC)
1282# 408 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1283!$acc loop seq
1284# 408 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1285#elif defined(MFC_OpenMP)
1286# 408 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1287
1288# 408 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1289#endif
1290 do i = 1, nb
1291 mu = qk_cons_vf(eqn_idx%bub%beg + 1 + (i - 1)*nmom)%sf(j, k, l)/nbub_sc
1292 sig = (qk_cons_vf(eqn_idx%bub%beg + 3 + (i - 1)*nmom)%sf(j, k, l)/nbub_sc - mu**2)**0.5_wp
1293
1294 mv(j, k, l, 1, i) = (mass_v0(i))*(mu - sig)**(3._wp)/(r0(i)**(3._wp))
1295 mv(j, k, l, 2, i) = (mass_v0(i))*(mu - sig)**(3._wp)/(r0(i)**(3._wp))
1296 mv(j, k, l, 3, i) = (mass_v0(i))*(mu + sig)**(3._wp)/(r0(i)**(3._wp))
1297 mv(j, k, l, 4, i) = (mass_v0(i))*(mu + sig)**(3._wp)/(r0(i)**(3._wp))
1298 end do
1299 end do
1300 end do
1301 end do
1302
1303 end subroutine s_initialize_mv
1304
1305 !> Initialize bubble internal pressures at quadrature nodes using isothermal relations from the Preston model.
1306 subroutine s_initialize_pb(qK_cons_vf, mv, pb)
1307
1308 type(scalar_field), dimension(sys_size), intent(in) :: qk_cons_vf
1309 real(stp), dimension(idwint(1)%beg:,idwint(2)%beg:,idwint(3)%beg:,1:,1:), intent(in) :: mv
1310 real(stp), dimension(idwint(1)%beg:,idwint(2)%beg:,idwint(3)%beg:,1:,1:), intent(inout) :: pb
1311 integer :: i, j, k, l
1312 real(wp) :: mu, sig, nbub_sc
1313
1314 do l = idwint(3)%beg, idwint(3)%end
1315 do k = idwint(2)%beg, idwint(2)%end
1316 do j = idwint(1)%beg, idwint(1)%end
1317 nbub_sc = qk_cons_vf(eqn_idx%bub%beg)%sf(j, k, l)
1318
1319
1320# 438 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1321#if defined(MFC_OpenACC)
1322# 438 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1323!$acc loop seq
1324# 438 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1325#elif defined(MFC_OpenMP)
1326# 438 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1327
1328# 438 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1329#endif
1330 do i = 1, nb
1331 mu = qk_cons_vf(eqn_idx%bub%beg + 1 + (i - 1)*nmom)%sf(j, k, l)/nbub_sc
1332 sig = (qk_cons_vf(eqn_idx%bub%beg + 3 + (i - 1)*nmom)%sf(j, k, l)/nbub_sc - mu**2)**0.5_wp
1333
1334 ! PRESTON (ISOTHERMAL)
1335 pb(j, k, l, 1, i) = (pb0(i))*(r0(i)**(3._wp))*(mass_g0(i) + mv(j, k, l, 1, &
1336 & i))/(mu - sig)**(3._wp)/(mass_g0(i) + mass_v0(i))
1337 pb(j, k, l, 2, i) = (pb0(i))*(r0(i)**(3._wp))*(mass_g0(i) + mv(j, k, l, 2, &
1338 & i))/(mu - sig)**(3._wp)/(mass_g0(i) + mass_v0(i))
1339 pb(j, k, l, 3, i) = (pb0(i))*(r0(i)**(3._wp))*(mass_g0(i) + mv(j, k, l, 3, &
1340 & i))/(mu + sig)**(3._wp)/(mass_g0(i) + mass_v0(i))
1341 pb(j, k, l, 4, i) = (pb0(i))*(r0(i)**(3._wp))*(mass_g0(i) + mv(j, k, l, 4, &
1342 & i))/(mu + sig)**(3._wp)/(mass_g0(i) + mass_v0(i))
1343 end do
1344 end do
1345 end do
1346 end do
1347
1348 end subroutine s_initialize_pb
1349
1350 !> Convert conserved variables (rho*alpha, rho*u, E, alpha) to primitives (rho, u, p, alpha). Conversion depends on model_eqns:
1351 !! each model has different variable sets and EOS.
1352 subroutine s_convert_conservative_to_primitive_variables(qK_cons_vf, q_T_sf, qK_prim_vf, ibounds)
1353
1354 use m_global_parameters_common, only: shear_indices ! Performance fix with AMDFlang
1355
1356 type(scalar_field), dimension(sys_size), intent(in) :: qk_cons_vf
1357 type(scalar_field), intent(inout) :: q_t_sf
1358 type(scalar_field), dimension(sys_size), intent(inout) :: qk_prim_vf
1359 type(int_bounds_info), dimension(1:3), intent(in) :: ibounds
1360
1361# 475 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1362 real(wp), dimension(num_fluids) :: alpha_k, alpha_rho_k
1363 real(wp), dimension(nb) :: nrtmp
1364 real(wp) :: rhoyks(1:num_species)
1365# 479 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1366 real(wp), dimension(2) :: re_k
1367 real(wp) :: rho_k, gamma_k, pi_inf_k, qv_k, dyn_pres_k
1368 real(wp) :: vftmp, nbub_sc
1369 real(wp) :: g_k
1370 real(wp) :: solid_partial_density
1371 real(wp) :: pres
1372 integer :: i, j, k, l !< Generic loop iterators
1373 real(wp) :: t
1374 real(wp) :: pres_mag
1375 real(wp) :: ga !< Lorentz factor (gamma in relativity)
1376 real(wp) :: b2 !< Magnetic field magnitude squared
1377 real(wp) :: b(3) !< Magnetic field components
1378 real(wp) :: m2 !< Relativistic momentum magnitude squared
1379 real(wp) :: s !< Dot product of the magnetic field and the relativistic momentum
1380 real(wp) :: w, dw !< W := rho*v*Ga**2; f = f(W) in Newton-Raphson
1381 real(wp) :: e, d !< Prim/Cons variables within Newton-Raphson iteration
1382 real(wp) :: f, dga_dw, dp_dw, df_dw !< Functions within Newton-Raphson iteration
1383 integer :: iter !< Newton-Raphson iteration counter
1384
1385
1386# 498 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1387
1388# 498 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1389#if defined(MFC_OpenACC)
1390# 498 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1391!$acc parallel loop collapse(3) gang vector default(present) private(alpha_K, alpha_rho_K, Re_K, nRtmp, rho_K, gamma_K, pi_inf_K, qv_K, dyn_pres_K, rhoYks, B, pres, vftmp, nbub_sc, G_K, &
1392# 498 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1393!$acc& solid_partial_density, T, pres_mag, Ga, B2, m2, S, W, dW, E, D, f, dGa_dW, dp_dW, df_dW, iter)
1394# 498 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1395#elif defined(MFC_OpenMP)
1396# 498 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1397
1398# 498 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1399
1400# 498 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1401
1402# 498 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1403!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(alpha_K, &
1404# 498 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1405!$omp& alpha_rho_K, Re_K, nRtmp, rho_K, gamma_K, pi_inf_K, qv_K, dyn_pres_K, rhoYks, B, pres, vftmp, nbub_sc, G_K, solid_partial_density, T, pres_mag, Ga, B2, m2, S, W, dW, E, D, f, dGa_dW, dp_dW, &
1406# 498 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1407!$omp& df_dW, iter)
1408# 498 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1409#endif
1410# 501 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1411 do l = ibounds(3)%beg, ibounds(3)%end
1412 do k = ibounds(2)%beg, ibounds(2)%end
1413 do j = ibounds(1)%beg, ibounds(1)%end
1414 dyn_pres_k = 0._wp
1415
1416 call s_compute_species_fraction(qk_cons_vf, j, k, l, alpha_rho_k, alpha_k)
1417
1418#ifdef MFC_GPU
1419 ! Device regions call the device-compiled scalar kernel directly.
1420 if (hypoelasticity) then
1421 call s_convert_species_to_mixture_variables_kernel(rho_k, gamma_k, pi_inf_k, qv_k, alpha_k, alpha_rho_k, &
1422 & re_k, g_k, gs_vc)
1423 else
1424 call s_convert_species_to_mixture_variables_kernel(rho_k, gamma_k, pi_inf_k, qv_k, alpha_k, alpha_rho_k, &
1425 & re_k)
1426 end if
1427#else
1428 ! Host execution uses the wrapper, which also stores requested diagnostics.
1429 if (hypoelasticity) then
1430 call s_convert_to_mixture_variables(qk_cons_vf, j, k, l, rho_k, gamma_k, pi_inf_k, qv_k, re_k, g_k, &
1431 & fluid_pp(:)%G)
1432 else
1433 call s_convert_to_mixture_variables(qk_cons_vf, j, k, l, rho_k, gamma_k, pi_inf_k, qv_k)
1434 end if
1435#endif
1436
1437 ! Relativistic MHD primitive variable recovery, Mignone & Bodo A&A (2006)
1438 if (relativity) then
1439 if (n == 0) then
1440 b(1) = bx0
1441 b(2) = qk_cons_vf(eqn_idx%B%beg)%sf(j, k, l)
1442 b(3) = qk_cons_vf(eqn_idx%B%beg + 1)%sf(j, k, l)
1443 else
1444 b(1) = qk_cons_vf(eqn_idx%B%beg)%sf(j, k, l)
1445 b(2) = qk_cons_vf(eqn_idx%B%beg + 1)%sf(j, k, l)
1446 b(3) = qk_cons_vf(eqn_idx%B%beg + 2)%sf(j, k, l)
1447 end if
1448 b2 = b(1)**2 + b(2)**2 + b(3)**2
1449
1450 m2 = 0._wp
1451
1452# 541 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1453#if defined(MFC_OpenACC)
1454# 541 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1455!$acc loop seq
1456# 541 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1457#elif defined(MFC_OpenMP)
1458# 541 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1459
1460# 541 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1461#endif
1462 do i = eqn_idx%mom%beg, eqn_idx%mom%end
1463 m2 = m2 + qk_cons_vf(i)%sf(j, k, l)**2
1464 end do
1465
1466 s = 0._wp
1467
1468# 547 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1469#if defined(MFC_OpenACC)
1470# 547 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1471!$acc loop seq
1472# 547 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1473#elif defined(MFC_OpenMP)
1474# 547 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1475
1476# 547 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1477#endif
1478 do i = 1, 3
1479 s = s + qk_cons_vf(eqn_idx%mom%beg + i - 1)%sf(j, k, l)*b(i)
1480 end do
1481
1482 e = qk_cons_vf(eqn_idx%E)%sf(j, k, l)
1483
1484 d = 0._wp
1485
1486# 555 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1487#if defined(MFC_OpenACC)
1488# 555 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1489!$acc loop seq
1490# 555 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1491#elif defined(MFC_OpenMP)
1492# 555 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1493
1494# 555 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1495#endif
1496 do i = 1, eqn_idx%cont%end
1497 d = d + qk_cons_vf(i)%sf(j, k, l)
1498 end do
1499
1500 ! Newton-Raphson
1501 w = e + d
1502
1503# 562 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1504#if defined(MFC_OpenACC)
1505# 562 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1506!$acc loop seq
1507# 562 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1508#elif defined(MFC_OpenMP)
1509# 562 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1510
1511# 562 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1512#endif
1513 do iter = 1, relativity_cons_to_prim_max_iter
1514 ! Lorentz factor from total enthalpy and magnetic field
1515 ga = (w + b2)*w/sqrt((w + b2)**2*w**2 - (m2*w**2 + s**2*(2*w + b2)))
1516 ! Thermal pressure from EOS
1517 pres = (w - d*ga)/((gamma_k + 1)*ga**2)
1518 f = w - pres + (1 - 1/(2*ga**2))*b2 - s**2/(2*w**2) - e - d
1519
1520 ! The first equation below corrects a typo in (Mignone & Bodo, 2006) m2*W**2 -> 2*m2*W**2, which would
1521 ! cancel with the 2* in other terms This corrected version is not used as the second equation
1522 ! empirically converges faster. First equation is kept for further investigation. dGa_dW = -Ga**3 * (
1523 ! S**2*(3*W**2+3*W*B2+B2**2) + m2*W**2 ) / (W**3 * (W+B2)**3) ! first (corrected)
1524 dga_dw = -ga**3*(2*s**2*(3*w**2 + 3*w*b2 + b2**2) + m2*w**2)/(2*w**3*(w + b2)**3) ! second (in paper)
1525
1526 dp_dw = (ga*(1 + d*dga_dw) - 2*w*dga_dw)/((gamma_k + 1)*ga**3)
1527 df_dw = 1 - dp_dw + (b2/ga**3)*dga_dw + s**2/w**3
1528
1529 dw = -f/df_dw
1530 w = w + dw
1531 if (abs(dw) < 1.e-12_wp*w) exit ! Relative convergence criterion
1532 end do
1533
1534 ! Recalculate pressure using converged W
1535 ga = (w + b2)*w/sqrt((w + b2)**2*w**2 - (m2*w**2 + s**2*(2*w + b2)))
1536 qk_prim_vf(eqn_idx%E)%sf(j, k, l) = (w - d*ga)/((gamma_k + 1)*ga**2)
1537
1538 ! Recover the other primitive variables
1539
1540# 589 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1541#if defined(MFC_OpenACC)
1542# 589 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1543!$acc loop seq
1544# 589 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1545#elif defined(MFC_OpenMP)
1546# 589 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1547
1548# 589 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1549#endif
1550 do i = 1, 3
1551 qk_prim_vf(eqn_idx%mom%beg + i - 1)%sf(j, k, l) = (qk_cons_vf(eqn_idx%mom%beg + i - 1)%sf(j, k, &
1552 & l) + (s/w)*b(i))/(w + b2)
1553 end do
1554 qk_prim_vf(1)%sf(j, k, l) = d/ga ! Hard-coded for single-component for now
1555
1556
1557# 596 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1558#if defined(MFC_OpenACC)
1559# 596 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1560!$acc loop seq
1561# 596 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1562#elif defined(MFC_OpenMP)
1563# 596 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1564
1565# 596 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1566#endif
1567 do i = eqn_idx%B%beg, eqn_idx%B%end
1568 qk_prim_vf(i)%sf(j, k, l) = qk_cons_vf(i)%sf(j, k, l)
1569 end do
1570
1571 cycle ! skip all the non-relativistic conversions below
1572 end if
1573
1574 if (chemistry) then
1575 ! Reacting flow: recover density from species partial densities, compute mass fractions Y_k = rhoY_k / rho
1576 rho_k = 0._wp
1577
1578# 607 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1579#if defined(MFC_OpenACC)
1580# 607 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1581!$acc loop seq
1582# 607 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1583#elif defined(MFC_OpenMP)
1584# 607 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1585
1586# 607 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1587#endif
1588 do i = eqn_idx%species%beg, eqn_idx%species%end
1589 rho_k = rho_k + max(0._wp, qk_cons_vf(i)%sf(j, k, l))
1590 end do
1591
1592
1593# 612 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1594#if defined(MFC_OpenACC)
1595# 612 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1596!$acc loop seq
1597# 612 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1598#elif defined(MFC_OpenMP)
1599# 612 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1600
1601# 612 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1602#endif
1603 do i = 1, eqn_idx%cont%end
1604 qk_prim_vf(i)%sf(j, k, l) = rho_k
1605 end do
1606
1607
1608# 617 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1609#if defined(MFC_OpenACC)
1610# 617 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1611!$acc loop seq
1612# 617 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1613#elif defined(MFC_OpenMP)
1614# 617 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1615
1616# 617 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1617#endif
1618 do i = eqn_idx%species%beg, eqn_idx%species%end
1619 qk_prim_vf(i)%sf(j, k, l) = max(0._wp, qk_cons_vf(i)%sf(j, k, l)/rho_k)
1620 end do
1621 else
1622 ! Non-reacting: partial densities are directly primitive (alpha_i * rho_i)
1623
1624# 623 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1625#if defined(MFC_OpenACC)
1626# 623 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1627!$acc loop seq
1628# 623 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1629#elif defined(MFC_OpenMP)
1630# 623 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1631
1632# 623 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1633#endif
1634 do i = 1, eqn_idx%cont%end
1635 qk_prim_vf(i)%sf(j, k, l) = qk_cons_vf(i)%sf(j, k, l)
1636 end do
1637 end if
1638
1639 if (enforce_density_floor_vc) rho_k = max(rho_k, sgm_eps)
1640
1641 ! Recover velocity from momentum: u = rho*u / rho, and accumulate dynamic pressure 0.5*rho*|u|^2
1642
1643# 632 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1644#if defined(MFC_OpenACC)
1645# 632 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1646!$acc loop seq
1647# 632 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1648#elif defined(MFC_OpenMP)
1649# 632 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1650
1651# 632 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1652#endif
1653 do i = eqn_idx%mom%beg, eqn_idx%mom%end
1654 qk_prim_vf(i)%sf(j, k, l) = qk_cons_vf(i)%sf(j, k, l)/rho_k
1655 dyn_pres_k = dyn_pres_k + 5.e-1_wp*qk_cons_vf(i)%sf(j, k, l)*qk_prim_vf(i)%sf(j, k, l)
1656 end do
1657
1658 if (chemistry) then
1659
1660# 639 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1661#if defined(MFC_OpenACC)
1662# 639 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1663!$acc loop seq
1664# 639 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1665#elif defined(MFC_OpenMP)
1666# 639 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1667
1668# 639 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1669#endif
1670 do i = 1, num_species
1671 rhoyks(i) = qk_cons_vf(eqn_idx%species%beg + i - 1)%sf(j, k, l)
1672 end do
1673
1674 t = q_t_sf%sf(j, k, l)
1675 end if
1676
1677 if (mhd) then
1678 if (n == 0) then
1679 pres_mag = 0.5_wp*(bx0**2 + qk_cons_vf(eqn_idx%B%beg)%sf(j, k, &
1680 & l)**2 + qk_cons_vf(eqn_idx%B%beg + 1)%sf(j, k, l)**2)
1681 else
1682 pres_mag = 0.5_wp*(qk_cons_vf(eqn_idx%B%beg)%sf(j, k, l)**2 + qk_cons_vf(eqn_idx%B%beg + 1)%sf(j, k, &
1683 & l)**2 + qk_cons_vf(eqn_idx%B%beg + 2)%sf(j, k, l)**2)
1684 end if
1685 else
1686 pres_mag = 0._wp
1687 end if
1688
1689 call s_compute_pressure(qk_cons_vf(eqn_idx%E)%sf(j, k, l), qk_cons_vf(eqn_idx%alf)%sf(j, k, l), dyn_pres_k, &
1690 & pi_inf_k, gamma_k, rho_k, qv_k, rhoyks, pres, t, pres_mag=pres_mag)
1691
1692 qk_prim_vf(eqn_idx%E)%sf(j, k, l) = pres
1693
1694 if (chemistry) then
1695 q_t_sf%sf(j, k, l) = t
1696 else if (heat_conduction) then
1697 q_t_sf%sf(j, k, l) = f_mixture_temperature(alpha_rho_k, pres, gamma_k, pi_inf_k)
1698 end if
1699
1700 if (bubbles_euler) then
1701 ! Recover bubble primitive variables: divide conserved moments by bubble number density
1702
1703# 672 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1704#if defined(MFC_OpenACC)
1705# 672 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1706!$acc loop seq
1707# 672 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1708#elif defined(MFC_OpenMP)
1709# 672 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1710
1711# 672 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1712#endif
1713 do i = 1, nb
1714 nrtmp(i) = qk_cons_vf(bubrs_vc(i))%sf(j, k, l)
1715 end do
1716
1717 vftmp = qk_cons_vf(eqn_idx%alf)%sf(j, k, l)
1718
1719 if (qbmm) then
1720 ! Get nb (constant across all R0 bins)
1721 nbub_sc = qk_cons_vf(eqn_idx%bub%beg)%sf(j, k, l)
1722
1723 ! Convert cons to prim
1724
1725# 684 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1726#if defined(MFC_OpenACC)
1727# 684 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1728!$acc loop seq
1729# 684 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1730#elif defined(MFC_OpenMP)
1731# 684 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1732
1733# 684 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1734#endif
1735 do i = eqn_idx%bub%beg, eqn_idx%bub%end
1736 qk_prim_vf(i)%sf(j, k, l) = qk_cons_vf(i)%sf(j, k, l)/nbub_sc
1737 end do
1738 ! Need to keep track of nb in the primitive variable list (converted back to true value before output)
1739 if (preserve_qbmm_number_vc) then
1740 qk_prim_vf(eqn_idx%bub%beg)%sf(j, k, l) = qk_cons_vf(eqn_idx%bub%beg)%sf(j, k, l)
1741 end if
1742 else
1743 if (adv_n) then
1744 qk_prim_vf(eqn_idx%n)%sf(j, k, l) = qk_cons_vf(eqn_idx%n)%sf(j, k, l)
1745 nbub_sc = qk_prim_vf(eqn_idx%n)%sf(j, k, l)
1746 else
1747 call s_comp_n_from_cons(vftmp, nrtmp, nbub_sc, weight)
1748 end if
1749
1750
1751# 700 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1752#if defined(MFC_OpenACC)
1753# 700 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1754!$acc loop seq
1755# 700 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1756#elif defined(MFC_OpenMP)
1757# 700 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1758
1759# 700 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1760#endif
1761 do i = eqn_idx%bub%beg, eqn_idx%bub%end
1762 qk_prim_vf(i)%sf(j, k, l) = qk_cons_vf(i)%sf(j, k, l)/nbub_sc
1763 end do
1764 end if
1765 end if
1766
1767 if (mhd) then
1768
1769# 708 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1770#if defined(MFC_OpenACC)
1771# 708 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1772!$acc loop seq
1773# 708 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1774#elif defined(MFC_OpenMP)
1775# 708 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1776
1777# 708 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1778#endif
1779 do i = eqn_idx%B%beg, eqn_idx%B%end
1780 qk_prim_vf(i)%sf(j, k, l) = qk_cons_vf(i)%sf(j, k, l)
1781 end do
1782 end if
1783
1784 if (hypoelasticity) then
1785
1786# 715 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1787#if defined(MFC_OpenACC)
1788# 715 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1789!$acc loop seq
1790# 715 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1791#elif defined(MFC_OpenMP)
1792# 715 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1793
1794# 715 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1795#endif
1796 do i = eqn_idx%stress%beg, eqn_idx%stress%end
1797 qk_prim_vf(i)%sf(j, k, l) = qk_cons_vf(i)%sf(j, k, l)/rho_k
1798 end do
1799 end if
1800
1801 if (cont_damage) then
1802 ! Recover D = U_D/m_s (damageable-solid partial mass), clamped to [0, 1]
1803 solid_partial_density = 0._wp
1804
1805# 724 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1806#if defined(MFC_OpenACC)
1807# 724 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1808!$acc loop seq
1809# 724 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1810#elif defined(MFC_OpenMP)
1811# 724 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1812
1813# 724 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1814#endif
1815 do i = 1, num_fluids
1816 if (gs_vc(i) > verysmall) then
1817 solid_partial_density = solid_partial_density + qk_cons_vf(eqn_idx%cont%beg + i - 1)%sf(j, k, l)
1818 end if
1819 end do
1820 qk_prim_vf(eqn_idx%damage)%sf(j, k, l) = min(max(qk_cons_vf(eqn_idx%damage)%sf(j, k, &
1821 & l)/max(solid_partial_density, verysmall), 0._wp), 1._wp)
1822 end if
1823
1824 if (hypoelasticity) then
1825 ! Elastic energy uses the undamaged modulus; tau^2/(4 G0 (1-D)) diverges as D -> 1
1826
1827# 736 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1828#if defined(MFC_OpenACC)
1829# 736 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1830!$acc loop seq
1831# 736 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1832#elif defined(MFC_OpenMP)
1833# 736 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1834
1835# 736 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1836#endif
1837 do i = eqn_idx%stress%beg, eqn_idx%stress%end
1838 qk_prim_vf(eqn_idx%E)%sf(j, k, l) = qk_prim_vf(eqn_idx%E)%sf(j, k, &
1839 & l) - f_elastic_energy(real(qk_prim_vf(i)%sf(j, k, l), wp), g_k, &
1840 & any(i == shear_indices))/gamma_k
1841 end do
1842 end if
1843
1844 if (.not. igr .or. num_fluids > 1) then
1845
1846# 745 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1847#if defined(MFC_OpenACC)
1848# 745 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1849!$acc loop seq
1850# 745 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1851#elif defined(MFC_OpenMP)
1852# 745 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1853
1854# 745 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1855#endif
1856 do i = eqn_idx%adv%beg, eqn_idx%adv%end
1857 qk_prim_vf(i)%sf(j, k, l) = qk_cons_vf(i)%sf(j, k, l)
1858 end do
1859 end if
1860
1861 if (surface_tension) then
1862 qk_prim_vf(eqn_idx%c)%sf(j, k, l) = qk_cons_vf(eqn_idx%c)%sf(j, k, l)
1863 end if
1864
1865 if (hyper_cleaning) qk_prim_vf(eqn_idx%psi)%sf(j, k, l) = qk_cons_vf(eqn_idx%psi)%sf(j, k, l)
1866 if (bubbles_lagrange .and. lagrange_beta_index_vc > 0) then
1867 qk_prim_vf(lagrange_beta_index_vc)%sf(j, k, l) = qk_cons_vf(lagrange_beta_index_vc)%sf(j, k, l)
1868 end if
1869 end do
1870 end do
1871 end do
1872
1873# 762 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1874#if defined(MFC_OpenACC)
1875# 762 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1876!$acc end parallel loop
1877# 762 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1878#elif defined(MFC_OpenMP)
1879# 762 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1880
1881# 762 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1882!$omp end target teams loop
1883# 762 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
1884#endif
1885
1887
1888 !> Convert primitives (rho, u, p, alpha) to conserved variables (rho*alpha, rho*u, E, alpha).
1889 impure subroutine s_convert_primitive_to_conservative_variables(q_prim_vf, q_cons_vf)
1890
1891 type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
1892 type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf
1893
1894 ! Density, specific heat ratio function, liquid stiffness function and dynamic pressure, as defined in the incompressible
1895 ! flow sense, respectively
1896 real(wp) :: rho
1897 real(wp) :: gamma
1898 real(wp) :: pi_inf
1899 real(wp) :: pres_i, alpha_i, alpha_rho_i, e_i
1900 real(wp) :: qv
1901 real(wp) :: dyn_pres
1902 real(wp) :: nbub, r3tmp
1903 real(wp), dimension(nb) :: rtmp
1904 real(wp) :: g
1905 real(wp) :: solid_partial_density
1906 real(wp), dimension(2) :: re_k
1907 integer :: i, j, k, l !< Generic loop iterators
1908 real(wp), dimension(num_species) :: ys
1909 real(wp) :: e_mix, mix_mol_weight, t
1910 real(wp) :: pres_mag
1911 real(wp) :: ga !< Lorentz factor (gamma in relativity)
1912 real(wp) :: h !< relativistic enthalpy
1913 real(wp) :: v2 !< Square of the velocity magnitude
1914 real(wp) :: b2 !< Square of the magnetic field magnitude
1915 real(wp) :: vdotb !< Dot product of the velocity and magnetic field vectors
1916 real(wp) :: b(3) !< Magnetic field components
1917
1918 pres_mag = 0._wp
1919
1920 g = 0._wp
1921
1922 ! Converting the primitive variables to the conservative variables
1923 do l = 0, p
1924 do k = 0, n
1925 do j = 0, m
1926 ! Obtaining the density, specific heat ratio function and the liquid stiffness function, respectively
1927 call s_convert_to_mixture_variables(q_prim_vf, j, k, l, rho, gamma, pi_inf, qv, re_k, g, fluid_pp(:)%G)
1928
1929 if (.not. igr .or. num_fluids > 1) then
1930 ! Transferring the advection equation(s) variable(s)
1931 do i = eqn_idx%adv%beg, eqn_idx%adv%end
1932 q_cons_vf(i)%sf(j, k, l) = q_prim_vf(i)%sf(j, k, l)
1933 end do
1934 end if
1935
1936 if (relativity) then
1937 if (n == 0) then
1938 b(1) = bx0
1939 b(2) = q_prim_vf(eqn_idx%B%beg)%sf(j, k, l)
1940 b(3) = q_prim_vf(eqn_idx%B%beg + 1)%sf(j, k, l)
1941 else
1942 b(1) = q_prim_vf(eqn_idx%B%beg)%sf(j, k, l)
1943 b(2) = q_prim_vf(eqn_idx%B%beg + 1)%sf(j, k, l)
1944 b(3) = q_prim_vf(eqn_idx%B%beg + 2)%sf(j, k, l)
1945 end if
1946
1947 v2 = 0._wp
1948 do i = eqn_idx%mom%beg, eqn_idx%mom%end
1949 v2 = v2 + q_prim_vf(i)%sf(j, k, l)**2
1950 end do
1951 if (v2 >= 1._wp) call s_mpi_abort('Error: v squared > 1 in s_convert_primitive_to_conservative_variables')
1952
1953 ga = 1._wp/sqrt(1._wp - v2)
1954
1955 h = 1._wp + (gamma + 1)*q_prim_vf(eqn_idx%E)%sf(j, k, l)/rho ! Assume perfect gas for now
1956
1957 b2 = 0._wp
1958 do i = eqn_idx%B%beg, eqn_idx%B%end
1959 b2 = b2 + q_prim_vf(i)%sf(j, k, l)**2
1960 end do
1961 if (n == 0) b2 = b2 + bx0**2
1962
1963 vdotb = 0._wp
1964 do i = 1, 3
1965 vdotb = vdotb + q_prim_vf(eqn_idx%mom%beg + i - 1)%sf(j, k, l)*b(i)
1966 end do
1967
1968 do i = 1, eqn_idx%cont%end
1969 q_cons_vf(i)%sf(j, k, l) = ga*q_prim_vf(i)%sf(j, k, l)
1970 end do
1971
1972 do i = eqn_idx%mom%beg, eqn_idx%mom%end
1973 q_cons_vf(i)%sf(j, k, l) = (rho*h*ga**2 + b2)*q_prim_vf(i)%sf(j, k, &
1974 & l) - vdotb*b(i - eqn_idx%mom%beg + 1)
1975 end do
1976
1977 q_cons_vf(eqn_idx%E)%sf(j, k, l) = rho*h*ga**2 - q_prim_vf(eqn_idx%E)%sf(j, k, &
1978 & l) + 0.5_wp*(b2 + v2*b2 - vdotb**2)
1979 ! Remove rest energy
1980 do i = 1, eqn_idx%cont%end
1981 q_cons_vf(eqn_idx%E)%sf(j, k, l) = q_cons_vf(eqn_idx%E)%sf(j, k, l) - q_cons_vf(i)%sf(j, k, l)
1982 end do
1983
1984 do i = eqn_idx%B%beg, eqn_idx%B%end
1985 q_cons_vf(i)%sf(j, k, l) = q_prim_vf(i)%sf(j, k, l)
1986 end do
1987
1988 cycle ! skip all the non-relativistic conversions below
1989 end if
1990
1991 ! Transferring the continuity equation(s) variable(s)
1992 do i = 1, eqn_idx%cont%end
1993 q_cons_vf(i)%sf(j, k, l) = q_prim_vf(i)%sf(j, k, l)
1994 end do
1995
1996 ! Zeroing out the dynamic pressure since it is computed iteratively by cycling through the velocity equations
1997 dyn_pres = 0._wp
1998
1999 ! Computing momenta and dynamic pressure from velocity
2000 do i = eqn_idx%mom%beg, eqn_idx%mom%end
2001 q_cons_vf(i)%sf(j, k, l) = rho*q_prim_vf(i)%sf(j, k, l)
2002 dyn_pres = dyn_pres + q_cons_vf(i)%sf(j, k, l)*q_prim_vf(i)%sf(j, k, l)/2._wp
2003 end do
2004
2005 if (chemistry) then
2006 ! Reacting mixture: compute conserved energy from species mass fractions and temperature
2007 do i = eqn_idx%species%beg, eqn_idx%species%end
2008 ys(i - eqn_idx%species%beg + 1) = q_prim_vf(i)%sf(j, k, l)
2009 q_cons_vf(i)%sf(j, k, l) = rho*q_prim_vf(i)%sf(j, k, l)
2010 end do
2011
2012 call get_mixture_molecular_weight(ys, mix_mol_weight)
2013 t = q_prim_vf(eqn_idx%E)%sf(j, k, l)*mix_mol_weight/(gas_constant*rho)
2014 call get_mixture_energy_mass(t, ys, e_mix)
2015
2016 q_cons_vf(eqn_idx%E)%sf(j, k, l) = dyn_pres + rho*e_mix
2017 else
2018 ! Computing the energy from the pressure
2019 if (mhd) then
2020 if (n == 0) then
2021 pres_mag = 0.5_wp*(bx0**2 + q_prim_vf(eqn_idx%B%beg)%sf(j, k, &
2022 & l)**2 + q_prim_vf(eqn_idx%B%beg + 1)%sf(j, k, l)**2)
2023 else
2024 pres_mag = 0.5_wp*(q_prim_vf(eqn_idx%B%beg)%sf(j, k, l)**2 + q_prim_vf(eqn_idx%B%beg + 1)%sf(j, &
2025 & k, l)**2 + q_prim_vf(eqn_idx%B%beg + 2)%sf(j, k, l)**2)
2026 end if
2027 ! MHD energy includes magnetic pressure contribution
2028 q_cons_vf(eqn_idx%E)%sf(j, k, l) = gamma*q_prim_vf(eqn_idx%E)%sf(j, k, &
2029 & l) + dyn_pres + pres_mag + pi_inf + qv
2030 else if (bubbles_euler .neqv. .true.) then
2031 ! Five-equation model (Allaire et al. JCP 2002): E = Gamma*p + 0.5*rho*|u|^2 + pi_inf + qv
2032 q_cons_vf(eqn_idx%E)%sf(j, k, l) = gamma*q_prim_vf(eqn_idx%E)%sf(j, k, l) + dyn_pres + pi_inf + qv
2033 else
2034 ! Bubble-augmented energy; qv is an energy density and is not diluted
2035 q_cons_vf(eqn_idx%E)%sf(j, k, l) = dyn_pres + (1._wp - q_prim_vf(eqn_idx%alf)%sf(j, k, &
2036 & l))*(gamma*q_prim_vf(eqn_idx%E)%sf(j, k, l) + pi_inf) + qv
2037 end if
2038 end if
2039
2040 ! Six-equation model (Saurel et al. JCP 2009): compute per-phase internal energies
2041 if (model_eqns == model_eqns_6eq) then
2042 do i = 1, num_fluids
2043 pres_i = q_prim_vf(eqn_idx%E)%sf(j, k, l)
2044 alpha_i = q_cons_vf(i + eqn_idx%adv%beg - 1)%sf(j, k, l)
2045 alpha_rho_i = q_cons_vf(i + eqn_idx%cont%beg - 1)%sf(j, k, l)
2046 call s_phase_internal_energy(pres_i, alpha_i, alpha_rho_i, i, e_i)
2047 q_cons_vf(i + eqn_idx%int_en%beg - 1)%sf(j, k, l) = e_i
2048 end do
2049 end if
2050
2051 if (bubbles_euler) then
2052 ! From prim: Compute nbub = (3/4pi) * \alpha / \bar{R^3}
2053 do i = 1, nb
2054 rtmp(i) = q_prim_vf(qbmm_idx%rs(i))%sf(j, k, l)
2055 end do
2056
2057 if (.not. qbmm) then
2058 if (adv_n) then
2059 q_cons_vf(eqn_idx%n)%sf(j, k, l) = q_prim_vf(eqn_idx%n)%sf(j, k, l)
2060 nbub = q_prim_vf(eqn_idx%n)%sf(j, k, l)
2061 else
2062 call s_comp_n_from_prim(real(q_prim_vf(eqn_idx%alf)%sf(j, k, l), kind=wp), rtmp, nbub, weight)
2063 end if
2064 else
2065 ! Initialize R3 averaging over R0 and R directions
2066 r3tmp = 0._wp
2067 do i = 1, nb
2068 r3tmp = r3tmp + weight(i)*0.5_wp*(rtmp(i) + sigr)**3._wp
2069 r3tmp = r3tmp + weight(i)*0.5_wp*(rtmp(i) - sigr)**3._wp
2070 end do
2071 ! Initialize nb
2072 nbub = 3._wp*q_prim_vf(eqn_idx%alf)%sf(j, k, l)/(4._wp*pi*r3tmp)
2073 end if
2074
2075 do i = eqn_idx%bub%beg, eqn_idx%bub%end
2076 q_cons_vf(i)%sf(j, k, l) = q_prim_vf(i)%sf(j, k, l)*nbub
2077 end do
2078 end if
2079
2080 if (mhd) then
2081 do i = eqn_idx%B%beg, eqn_idx%B%end
2082 q_cons_vf(i)%sf(j, k, l) = q_prim_vf(i)%sf(j, k, l)
2083 end do
2084 end if
2085
2086 if (hypoelasticity) then
2087 ! adding the elastic contribution Multiply \tau to \rho \tau
2088 do i = eqn_idx%stress%beg, eqn_idx%stress%end
2089 q_cons_vf(i)%sf(j, k, l) = rho*q_prim_vf(i)%sf(j, k, l)
2090 end do
2091 end if
2092
2093 if (hypoelasticity) then
2094 ! Elastic energy uses the undamaged modulus
2095 do i = eqn_idx%stress%beg, eqn_idx%stress%end
2096 ! Elastic energy addition (guard skips when G near zero from alpha undershoot)
2097 if (g > verysmall) then
2098 q_cons_vf(eqn_idx%E)%sf(j, k, l) = q_cons_vf(eqn_idx%E)%sf(j, k, l) + (q_prim_vf(i)%sf(j, k, &
2099 & l)**2._wp)/max(4._wp*g, verysmall)
2100 ! Double for shear stresses
2101 if (any(i == shear_indices)) then
2102 q_cons_vf(eqn_idx%E)%sf(j, k, l) = q_cons_vf(eqn_idx%E)%sf(j, k, l) + (q_prim_vf(i)%sf(j, k, &
2103 & l)**2._wp)/max(4._wp*g, verysmall)
2104 end if
2105 end if
2106 end do
2107 end if
2108
2109 if (surface_tension) then
2110 q_cons_vf(eqn_idx%c)%sf(j, k, l) = q_prim_vf(eqn_idx%c)%sf(j, k, l)
2111 end if
2112
2113 if (cont_damage) then
2114 ! U_D = m_s*D (damageable-solid partial mass)
2115 solid_partial_density = 0._wp
2116 do i = 1, num_fluids
2117 if (fluid_pp(i)%G > verysmall) then
2118 solid_partial_density = solid_partial_density + q_prim_vf(eqn_idx%cont%beg + i - 1)%sf(j, k, l)
2119 end if
2120 end do
2121 q_cons_vf(eqn_idx%damage)%sf(j, k, l) = solid_partial_density*q_prim_vf(eqn_idx%damage)%sf(j, k, l)
2122 end if
2123
2124 if (hyper_cleaning) q_cons_vf(eqn_idx%psi)%sf(j, k, l) = q_prim_vf(eqn_idx%psi)%sf(j, k, l)
2125 end do
2126 end do
2127 end do
2128
2130
2131 !> Convert primitive variables to Eulerian flux variables.
2132 subroutine s_convert_primitive_to_flux_variables(qK_prim_vf, FK_vf, FK_src_vf, is1, is2, is3, s2b, s3b, dir_idx_in, &
2133 & dir_flg_in, hll_u_interface_in)
2134
2135 integer, intent(in) :: s2b, s3b
2136 !> Working-direction mapping, passed explicitly: it is simulation state (m_global_parameters), and use-associating it into
2137 !! this common kernel spills registers on AMD OpenMP offload.
2138 integer, dimension(3), intent(in) :: dir_idx_in
2139 real(wp), dimension(3), intent(in) :: dir_flg_in
2140 logical, intent(in) :: hll_u_interface_in
2141 real(wp), dimension(0:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:), intent(in) :: qk_prim_vf
2142 real(wp), dimension(0:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:), intent(inout) :: fk_vf
2143 real(wp), dimension(0:,idwbuff(2)%beg:,idwbuff(3)%beg:,eqn_idx%adv%beg:), intent(inout) :: fk_src_vf
2144 type(int_bounds_info), intent(in) :: is1, is2, is3
2145
2146 ! Partial densities, density, velocity, pressure, energy, advection variables, the specific heat ratio and liquid stiffness
2147 ! functions, the shear and volume Reynolds numbers and the Weber numbers
2148
2149# 1033 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2150 real(wp), dimension(num_fluids) :: alpha_rho_k
2151 real(wp), dimension(num_fluids) :: alpha_k
2152 real(wp), dimension(num_vels) :: vel_k
2153 real(wp), dimension(num_species) :: y_k
2154# 1038 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2155 real(wp) :: rho_k
2156 real(wp) :: vel_k_sum
2157 real(wp) :: pres_k
2158 real(wp) :: e_k
2159 real(wp) :: gamma_k
2160 real(wp) :: pi_inf_k
2161 real(wp) :: qv_k
2162 real(wp), dimension(2) :: re_k
2163 real(wp) :: g_k
2164 real(wp) :: blkmod1_k, blkmod2_k, k_k
2165 real(wp) :: t_k, mix_mol_weight, r_gas
2166 integer :: i, j, k, l !< Generic loop iterators
2167
2168 is1b = is1%beg; is1e = is1%end
2169 is2b = is2%beg; is2e = is2%end
2170 is3b = is3%beg; is3e = is3%end
2171
2172
2173# 1055 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2174#if defined(MFC_OpenACC)
2175# 1055 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2176!$acc update device(is1b, is2b, is3b, is1e, is2e, is3e)
2177# 1055 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2178#elif defined(MFC_OpenMP)
2179# 1055 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2180!$omp target update to(is1b, is2b, is3b, is1e, is2e, is3e)
2181# 1055 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2182#endif
2183
2184 ! Computing the flux variables from the primitive variables, without accounting for the contribution of either viscosity or
2185 ! capillarity
2186
2187# 1059 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2188
2189# 1059 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2190#if defined(MFC_OpenACC)
2191# 1059 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2192!$acc parallel loop collapse(3) gang vector default(present) &
2193# 1059 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2194!$acc& private(alpha_rho_K, vel_K, alpha_K, Re_K, Y_K, rho_K, vel_K_sum, pres_K, E_K, gamma_K, pi_inf_K, qv_K, G_K, blkmod1_K, blkmod2_K, K_K, T_K, mix_mol_weight, R_gas) &
2195# 1059 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2196!$acc& copyin(dir_idx_in, dir_flg_in, hll_u_interface_in)
2197# 1059 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2198#elif defined(MFC_OpenMP)
2199# 1059 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2200
2201# 1059 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2202
2203# 1059 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2204
2205# 1059 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2206!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) &
2207# 1059 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2208!$omp& private(alpha_rho_K, vel_K, alpha_K, Re_K, Y_K, rho_K, vel_K_sum, pres_K, E_K, gamma_K, pi_inf_K, qv_K, G_K, blkmod1_K, blkmod2_K, K_K, T_K, mix_mol_weight, R_gas) &
2209# 1059 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2210!$omp& map(to:dir_idx_in, dir_flg_in, hll_u_interface_in)
2211# 1059 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2212#endif
2213# 1062 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2214 do l = is3b, is3e
2215 do k = is2b, is2e
2216 do j = is1b, is1e
2217
2218# 1065 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2219#if defined(MFC_OpenACC)
2220# 1065 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2221!$acc loop seq
2222# 1065 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2223#elif defined(MFC_OpenMP)
2224# 1065 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2225
2226# 1065 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2227#endif
2228 do i = 1, eqn_idx%cont%end
2229 alpha_rho_k(i) = qk_prim_vf(j, k, l, i)
2230 end do
2231
2232
2233# 1070 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2234#if defined(MFC_OpenACC)
2235# 1070 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2236!$acc loop seq
2237# 1070 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2238#elif defined(MFC_OpenMP)
2239# 1070 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2240
2241# 1070 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2242#endif
2243 do i = eqn_idx%adv%beg, eqn_idx%adv%end
2244 alpha_k(i - eqn_idx%E) = qk_prim_vf(j, k, l, i)
2245 end do
2246
2247
2248# 1075 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2249#if defined(MFC_OpenACC)
2250# 1075 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2251!$acc loop seq
2252# 1075 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2253#elif defined(MFC_OpenMP)
2254# 1075 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2255
2256# 1075 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2257#endif
2258 do i = 1, num_vels
2259 vel_k(i) = qk_prim_vf(j, k, l, eqn_idx%cont%end + i)
2260 end do
2261
2262 vel_k_sum = 0._wp
2263
2264# 1081 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2265#if defined(MFC_OpenACC)
2266# 1081 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2267!$acc loop seq
2268# 1081 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2269#elif defined(MFC_OpenMP)
2270# 1081 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2271
2272# 1081 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2273#endif
2274 do i = 1, num_vels
2275 vel_k_sum = vel_k_sum + vel_k(i)**2._wp
2276 end do
2277
2278 pres_k = qk_prim_vf(j, k, l, eqn_idx%E)
2279 if (hypoelasticity) then
2280 call s_convert_species_to_mixture_variables_kernel(rho_k, gamma_k, pi_inf_k, qv_k, alpha_k, alpha_rho_k, &
2281 & re_k, g_k, gs_vc)
2282 else
2283 call s_convert_species_to_mixture_variables_kernel(rho_k, gamma_k, pi_inf_k, qv_k, alpha_k, alpha_rho_k, &
2284 & re_k)
2285 end if
2286
2287 ! Computing the energy from the pressure
2288
2289 if (chemistry) then
2290
2291# 1098 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2292#if defined(MFC_OpenACC)
2293# 1098 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2294!$acc loop seq
2295# 1098 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2296#elif defined(MFC_OpenMP)
2297# 1098 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2298
2299# 1098 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2300#endif
2301 do i = eqn_idx%species%beg, eqn_idx%species%end
2302 y_k(i - eqn_idx%species%beg + 1) = qk_prim_vf(j, k, l, i)
2303 end do
2304 ! Computing the energy from the internal energy of the mixture
2305 call get_mixture_molecular_weight(y_k, mix_mol_weight)
2306 r_gas = gas_constant/mix_mol_weight
2307 t_k = pres_k/rho_k/r_gas
2308 call get_mixture_energy_mass(t_k, y_k, e_k)
2309 e_k = rho_k*e_k + 5.e-1_wp*rho_k*vel_k_sum
2310 else
2311 ! Computing the energy from the pressure
2312 e_k = gamma_k*pres_k + pi_inf_k + 5.e-1_wp*rho_k*vel_k_sum + qv_k
2313 end if
2314
2315 ! mass flux, this should be \alpha_i \rho_i u_i
2316
2317# 1114 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2318#if defined(MFC_OpenACC)
2319# 1114 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2320!$acc loop seq
2321# 1114 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2322#elif defined(MFC_OpenMP)
2323# 1114 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2324
2325# 1114 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2326#endif
2327 do i = 1, eqn_idx%cont%end
2328 fk_vf(j, k, l, i) = alpha_rho_k(i)*vel_k(dir_idx_in(1))
2329 end do
2330
2331
2332# 1119 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2333#if defined(MFC_OpenACC)
2334# 1119 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2335!$acc loop seq
2336# 1119 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2337#elif defined(MFC_OpenMP)
2338# 1119 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2339
2340# 1119 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2341#endif
2342 do i = 1, num_vels
2343 fk_vf(j, k, l, &
2344 & eqn_idx%cont%end + dir_idx_in(i)) = rho_k*vel_k(dir_idx_in(1))*vel_k(dir_idx_in(i)) &
2345 & + pres_k*dir_flg_in(dir_idx_in(i))
2346 end do
2347
2348 ! energy flux, u(E+p)
2349 fk_vf(j, k, l, eqn_idx%E) = vel_k(dir_idx_in(1))*(e_k + pres_k)
2350
2351 ! Species advection Flux, \rho*u*Y
2352 if (chemistry) then
2353
2354# 1131 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2355#if defined(MFC_OpenACC)
2356# 1131 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2357!$acc loop seq
2358# 1131 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2359#elif defined(MFC_OpenMP)
2360# 1131 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2361
2362# 1131 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2363#endif
2364 do i = 1, num_species
2365 fk_vf(j, k, l, i - 1 + eqn_idx%species%beg) = vel_k(dir_idx_in(1))*(rho_k*y_k(i))
2366 end do
2367 end if
2368
2369 ! Match the volume-fraction flux representation exported by the Riemann solver. HLL Method 1: zero alpha
2370 ! flux plus per-fluid interface-alpha source traces. Hypoelastic HLLD folds every non-conservative term
2371 ! into its augmented flux (adv_src_mode_none), so its source trace is zero; for this cell-local conversion
2372 ! the fold collapses exactly to -/+ K*u_n on the two volume-fraction rows (K = 0 without alt_soundspeed),
2373 ! with the same two-fluid longitudinal-modulus K as the HLLD kernel (num_fluids = 2 is checker-enforced).
2374 ! MHD HLLD keeps the per-fluid-trace representation it has always used. HLL Method 2, HLLC, and LF use the
2375 ! shared-velocity representation below.
2376 if (riemann_solver == riemann_solver_hlld) then
2377 if (hypoelasticity) then
2378 k_k = 0._wp
2379 ! The fluid-2 subscripts must not be compiled when case optimization
2380 ! bakes num_fluids = 1 (amdflang rejects them at compile time); the
2381 ! checker prohibits hypoelastic HLLD there, so the block is dead code.
2382# 1151 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2383 if (alt_soundspeed) then
2384 blkmod1_k = f_bulk_modulus(pres_k, gammas(1), pi_infs(1)) + (4._wp/3._wp)*gs_vc(1)
2385 blkmod2_k = f_bulk_modulus(pres_k, gammas(2), pi_infs(2)) + (4._wp/3._wp)*gs_vc(2)
2386 k_k = alpha_k(1)*alpha_k(2)*(blkmod2_k - blkmod1_k)/(alpha_k(1)*blkmod2_k + alpha_k(2) &
2387 & *blkmod1_k + verysmall)
2388 end if
2389# 1158 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2390
2391# 1158 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2392#if defined(MFC_OpenACC)
2393# 1158 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2394!$acc loop seq
2395# 1158 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2396#elif defined(MFC_OpenMP)
2397# 1158 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2398
2399# 1158 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2400#endif
2401 do i = eqn_idx%adv%beg, eqn_idx%adv%end
2402 fk_vf(j, k, l, i) = 0._wp
2403 fk_src_vf(j, k, l, i) = 0._wp
2404 end do
2405 fk_vf(j, k, l, eqn_idx%adv%beg) = -k_k*vel_k(dir_idx_in(1))
2406 fk_vf(j, k, l, eqn_idx%adv%end) = k_k*vel_k(dir_idx_in(1))
2407 else
2408
2409# 1166 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2410#if defined(MFC_OpenACC)
2411# 1166 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2412!$acc loop seq
2413# 1166 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2414#elif defined(MFC_OpenMP)
2415# 1166 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2416
2417# 1166 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2418#endif
2419 do i = eqn_idx%adv%beg, eqn_idx%adv%end
2420 fk_vf(j, k, l, i) = 0._wp
2421 fk_src_vf(j, k, l, i) = alpha_k(i - eqn_idx%E)
2422 end do
2423 end if
2424 else if (riemann_solver == riemann_solver_hll .and. .not. hll_u_interface_in) then
2425
2426# 1173 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2427#if defined(MFC_OpenACC)
2428# 1173 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2429!$acc loop seq
2430# 1173 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2431#elif defined(MFC_OpenMP)
2432# 1173 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2433
2434# 1173 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2435#endif
2436 do i = eqn_idx%adv%beg, eqn_idx%adv%end
2437 fk_vf(j, k, l, i) = 0._wp
2438 fk_src_vf(j, k, l, i) = alpha_k(i - eqn_idx%E)
2439 end do
2440 else
2441 ! Could be bubbles_euler!
2442
2443# 1180 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2444#if defined(MFC_OpenACC)
2445# 1180 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2446!$acc loop seq
2447# 1180 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2448#elif defined(MFC_OpenMP)
2449# 1180 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2450
2451# 1180 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2452#endif
2453 do i = eqn_idx%adv%beg, eqn_idx%adv%end
2454 fk_vf(j, k, l, i) = vel_k(dir_idx_in(1))*alpha_k(i - eqn_idx%E)
2455 end do
2456
2457
2458# 1185 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2459#if defined(MFC_OpenACC)
2460# 1185 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2461!$acc loop seq
2462# 1185 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2463#elif defined(MFC_OpenMP)
2464# 1185 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2465
2466# 1185 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2467#endif
2468 do i = eqn_idx%adv%beg, eqn_idx%adv%end
2469 fk_src_vf(j, k, l, i) = vel_k(dir_idx_in(1))
2470 end do
2471 end if
2472 end do
2473 end do
2474 end do
2475
2476# 1193 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2477#if defined(MFC_OpenACC)
2478# 1193 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2479!$acc end parallel loop
2480# 1193 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2481#elif defined(MFC_OpenMP)
2482# 1193 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2483
2484# 1193 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2485!$omp end target teams loop
2486# 1193 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2487#endif
2488
2490
2491 !> Compute partial densities and volume fractions
2492 subroutine s_compute_species_fraction(q_vf, k, l, r, alpha_rho_K, alpha_K)
2493
2494
2495# 1200 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2496#ifdef _CRAYFTN
2497# 1200 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2498#if MFC_OpenACC
2499# 1200 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2500!$acc routine seq
2501# 1200 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2502#elif MFC_OpenMP
2503# 1200 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2504
2505# 1200 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2506
2507# 1200 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2508!$omp declare target device_type(any)
2509# 1200 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2510#else
2511# 1200 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2512!DIR$ NOINLINE s_compute_species_fraction
2513# 1200 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2514#endif
2515# 1200 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2516#elif MFC_OpenACC
2517# 1200 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2518!$acc routine seq
2519# 1200 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2520#elif MFC_OpenMP
2521# 1200 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2522
2523# 1200 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2524
2525# 1200 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2526!$omp declare target device_type(any)
2527# 1200 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2528#endif
2529 type(scalar_field), dimension(sys_size), intent(in) :: q_vf
2530 integer, intent(in) :: k, l, r
2531# 1206 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2532 real(wp), dimension(num_fluids), intent(out) :: alpha_rho_k, alpha_k
2533# 1208 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2534 integer :: i
2535 real(wp) :: alpha_k_sum
2536
2537 if (num_fluids == 1) then
2538 alpha_rho_k(1) = q_vf(eqn_idx%cont%beg)%sf(k, l, r)
2539 if (igr .or. bubbles_euler) then
2540 alpha_k(1) = 1._wp
2541 else
2542 alpha_k(1) = q_vf(eqn_idx%adv%beg)%sf(k, l, r)
2543 end if
2544 else
2545 if (igr) then
2546 do i = 1, num_fluids - 1
2547 alpha_rho_k(i) = q_vf(i)%sf(k, l, r)
2548 alpha_k(i) = q_vf(eqn_idx%adv%beg + i - 1)%sf(k, l, r)
2549 end do
2550 alpha_rho_k(num_fluids) = q_vf(num_fluids)%sf(k, l, r)
2551 alpha_k(num_fluids) = 1._wp - sum(alpha_k(1:num_fluids - 1))
2552 else
2553 do i = 1, num_fluids
2554 alpha_rho_k(i) = q_vf(i)%sf(k, l, r)
2555 alpha_k(i) = q_vf(eqn_idx%adv%beg + i - 1)%sf(k, l, r)
2556 end do
2557 end if
2558 end if
2559
2560 if (mpp_lim) then
2561 alpha_k_sum = 0._wp
2562 do i = 1, num_fluids
2563 alpha_rho_k(i) = max(0._wp, alpha_rho_k(i))
2564 alpha_k(i) = min(max(0._wp, alpha_k(i)), 1._wp)
2565 alpha_k_sum = alpha_k_sum + alpha_k(i)
2566 end do
2567 alpha_k = alpha_k/max(alpha_k_sum, 1.e-16_wp)
2568 end if
2569
2570 if (num_fluids == 1 .and. bubbles_euler) alpha_k(1) = q_vf(eqn_idx%adv%beg)%sf(k, l, r)
2571
2572 end subroutine s_compute_species_fraction
2573
2574 !> Deallocate fluid property arrays and post-processing fields allocated during module initialization.
2576
2577 if (allocated(rho_sf)) deallocate (rho_sf, gamma_sf, pi_inf_sf)
2578
2579#ifdef MFC_DEBUG
2580# 1253 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2581 block
2582# 1253 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2583 use iso_fortran_env, only: output_unit
2584# 1253 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2585
2586# 1253 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2587 print *, 'm_variables_conversion.fpp:1253: ', '@:DEALLOCATE(gammas, isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps, Gs_vc, eoss, fluid_k_therm)'
2588# 1253 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2589
2590# 1253 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2591 call flush (output_unit)
2592# 1253 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2593 end block
2594# 1253 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2595#endif
2596# 1253 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2597
2598# 1253 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2599#if defined(MFC_OpenACC)
2600# 1253 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2601!$acc exit data delete(gammas, isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps, Gs_vc, eoss, fluid_k_therm)
2602# 1253 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2603#elif defined(MFC_OpenMP)
2604# 1253 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2605!$omp target exit data map(release:gammas, isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps, Gs_vc, eoss, fluid_k_therm)
2606# 1253 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2607#endif
2608# 1253 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2609 deallocate (gammas, isentrope_n, pi_infs, isentrope_b, cvs, qvs, qvps, gs_vc, eoss, fluid_k_therm)
2610 if (allocated(bubrs_vc)) then
2611#ifdef MFC_DEBUG
2612# 1255 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2613 block
2614# 1255 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2615 use iso_fortran_env, only: output_unit
2616# 1255 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2617
2618# 1255 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2619 print *, 'm_variables_conversion.fpp:1255: ', '@:DEALLOCATE(bubrs_vc)'
2620# 1255 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2621
2622# 1255 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2623 call flush (output_unit)
2624# 1255 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2625 end block
2626# 1255 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2627#endif
2628# 1255 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2629
2630# 1255 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2631#if defined(MFC_OpenACC)
2632# 1255 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2633!$acc exit data delete(bubrs_vc)
2634# 1255 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2635#elif defined(MFC_OpenMP)
2636# 1255 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2637!$omp target exit data map(release:bubrs_vc)
2638# 1255 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2639#endif
2640# 1255 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2641 deallocate (bubrs_vc)
2642 end if
2643 if (allocated(res_vc)) then
2644#ifdef MFC_DEBUG
2645# 1258 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2646 block
2647# 1258 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2648 use iso_fortran_env, only: output_unit
2649# 1258 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2650
2651# 1258 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2652 print *, 'm_variables_conversion.fpp:1258: ', '@:DEALLOCATE(Res_vc)'
2653# 1258 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2654
2655# 1258 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2656 call flush (output_unit)
2657# 1258 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2658 end block
2659# 1258 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2660#endif
2661# 1258 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2662
2663# 1258 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2664#if defined(MFC_OpenACC)
2665# 1258 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2666!$acc exit data delete(Res_vc)
2667# 1258 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2668#elif defined(MFC_OpenMP)
2669# 1258 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2670!$omp target exit data map(release:Res_vc)
2671# 1258 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2672#endif
2673# 1258 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2674 deallocate (res_vc)
2675 end if
2676
2678
2679 !> Mixture coefficients of one state. Under bubbles_euler with num_fluids == 1 the sole advection slot aliases the void fraction
2680 !! (eqn_idx%alf == eqn_idx%adv%end), so alpha is not a composition there and the coefficients are the liquid's. Clipping stays
2681 !! with callers; it differs between solvers and cannot coincide with that case, as mpp_lim requires num_fluids > 1.
2682 subroutine s_compute_mixture_coefficients(alpha_rho_K, alpha_K, rho_K, gamma_K, pi_inf_K, qv_K)
2683
2684
2685# 1268 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2686#ifdef _CRAYFTN
2687# 1268 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2688#if MFC_OpenACC
2689# 1268 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2690!$acc routine seq
2691# 1268 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2692#elif MFC_OpenMP
2693# 1268 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2694
2695# 1268 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2696
2697# 1268 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2698!$omp declare target device_type(any)
2699# 1268 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2700#else
2701# 1268 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2702!DIR$ INLINEALWAYS s_compute_mixture_coefficients
2703# 1268 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2704#endif
2705# 1268 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2706#elif MFC_OpenACC
2707# 1268 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2708!$acc routine seq
2709# 1268 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2710#elif MFC_OpenMP
2711# 1268 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2712
2713# 1268 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2714
2715# 1268 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2716!$omp declare target device_type(any)
2717# 1268 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2718#endif
2719
2720# 1273 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2721 real(wp), dimension(num_fluids), intent(in) :: alpha_rho_k, alpha_k
2722# 1275 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2723 real(wp), intent(out) :: rho_k, gamma_k, pi_inf_k, qv_k
2724 real(wp) :: gamma_i, pi_inf_i, dpi_i, dgamma_i
2725 real(wp) :: rho_i, alpha_i, alpha_rho_i
2726 integer :: i !< Loop iterator over fluids
2727
2728 ! The bubbly closure is written for one carrier liquid, which keeps its own coefficients
2729 ! undiluted: Gamma_l*p_l = (E - rho|u|^2/2)/(1 - alf) - Pi_inf_l, the void entering only through
2730 ! the (1 - alf) that s_compute_pressure applies. There is nothing to sum - the last advection
2731 ! slot is the void, not a material - and the checker holds num_fluids <= 2 here.
2732 if (bubbles_euler) then
2733 rho_k = alpha_rho_k(1)
2734 gamma_k = gammas(1)
2735 pi_inf_k = pi_infs(1)
2736 ! Energy per unit volume, as below: alpha_rho_K(1) is the liquid partial density
2737 qv_k = alpha_rho_k(1)*qvs(1)
2738 else
2739 rho_k = 0._wp
2740 gamma_k = 0._wp
2741 pi_inf_k = 0._wp
2742 qv_k = 0._wp
2743
2744
2745# 1296 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2746#if defined(MFC_OpenACC)
2747# 1296 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2748!$acc loop seq
2749# 1296 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2750#elif defined(MFC_OpenMP)
2751# 1296 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2752
2753# 1296 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2754#endif
2755 do i = 1, num_fluids
2756 rho_k = rho_k + alpha_rho_k(i)
2757 alpha_rho_i = alpha_rho_k(i)
2758 alpha_i = alpha_k(i)
2759 call s_phase_coefficients(alpha_rho_i, alpha_i, i, rho_i, gamma_i, pi_inf_i, dpi_i, dgamma_i)
2760 gamma_k = gamma_k + alpha_k(i)*gamma_i
2761 pi_inf_k = pi_inf_k + alpha_k(i)*pi_inf_i
2762 qv_k = qv_k + alpha_rho_k(i)*qvs(i)
2763 end do
2764 end if
2765
2766 end subroutine s_compute_mixture_coefficients
2767
2768 !> Time derivative of the mixture coefficients, mirroring s_compute_mixture_coefficients.
2769 subroutine s_compute_mixture_coefficients_dt(dalpha_rho_dt, dadv_dt, alpha_rho, adv, drho_dt, dgamma_dt, dpi_inf_dt, dqv_dt)
2770
2771
2772# 1313 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2773#ifdef _CRAYFTN
2774# 1313 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2775#if MFC_OpenACC
2776# 1313 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2777!$acc routine seq
2778# 1313 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2779#elif MFC_OpenMP
2780# 1313 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2781
2782# 1313 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2783
2784# 1313 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2785!$omp declare target device_type(any)
2786# 1313 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2787#else
2788# 1313 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2789!DIR$ INLINEALWAYS s_compute_mixture_coefficients_dt
2790# 1313 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2791#endif
2792# 1313 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2793#elif MFC_OpenACC
2794# 1313 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2795!$acc routine seq
2796# 1313 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2797#elif MFC_OpenMP
2798# 1313 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2799
2800# 1313 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2801
2802# 1313 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2803!$omp declare target device_type(any)
2804# 1313 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2805#endif
2806
2807# 1318 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2808 real(wp), dimension(num_fluids), intent(in) :: dalpha_rho_dt, dadv_dt, alpha_rho, adv
2809# 1320 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2810 real(wp), intent(out) :: drho_dt, dgamma_dt, dpi_inf_dt, dqv_dt
2811 real(wp) :: rho_i, gamma_i, pi_inf_i, dpi_i, dgamma_i, alpha_i, alpha_rho_i
2812 integer :: i !< Loop iterator over fluids
2813
2814 dgamma_dt = 0._wp
2815 dpi_inf_dt = 0._wp
2816 dqv_dt = 0._wp
2817
2818 if (num_fluids == 1 .and. bubbles_euler) then
2819 ! Fluid 1's coefficients are constants here, so only rho varies.
2820 drho_dt = dalpha_rho_dt(1)
2821 else
2822 drho_dt = 0._wp
2823
2824
2825# 1334 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2826#if defined(MFC_OpenACC)
2827# 1334 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2828!$acc loop seq
2829# 1334 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2830#elif defined(MFC_OpenMP)
2831# 1334 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2832
2833# 1334 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2834#endif
2835 do i = 1, num_fluids
2836 drho_dt = drho_dt + dalpha_rho_dt(i)
2837 alpha_rho_i = alpha_rho(i)
2838 alpha_i = adv(i)
2839 call s_phase_coefficients(alpha_rho_i, alpha_i, i, rho_i, gamma_i, pi_inf_i, dpi_i, dgamma_i)
2840 ! d(alpha X(rho_i))/dt with rho_i = alpha_rho/alpha; the alpha in dX/dt cancels
2841 dgamma_dt = dgamma_dt + dadv_dt(i)*gamma_i + dgamma_i*(dalpha_rho_dt(i) - rho_i*dadv_dt(i))
2842 dpi_inf_dt = dpi_inf_dt + dadv_dt(i)*pi_inf_i + dpi_i*(dalpha_rho_dt(i) - rho_i*dadv_dt(i))
2843 dqv_dt = dqv_dt + dalpha_rho_dt(i)*qvs(i)
2844 end do
2845 end if
2846
2848
2849 !> Total energy per unit volume, thermodynamic terms only. Callers add magnetic and elastic energy, which are not
2850 !! equation-of-state terms. The chemistry and relativistic branches use a different relation and stay open-coded.
2851 subroutine s_compute_energy(pres, alpha_rho_K, alpha_K, vel_sum, E)
2852
2853
2854# 1353 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2855#ifdef _CRAYFTN
2856# 1353 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2857#if MFC_OpenACC
2858# 1353 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2859!$acc routine seq
2860# 1353 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2861#elif MFC_OpenMP
2862# 1353 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2863
2864# 1353 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2865
2866# 1353 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2867!$omp declare target device_type(any)
2868# 1353 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2869#else
2870# 1353 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2871!DIR$ INLINEALWAYS s_compute_energy
2872# 1353 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2873#endif
2874# 1353 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2875#elif MFC_OpenACC
2876# 1353 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2877!$acc routine seq
2878# 1353 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2879#elif MFC_OpenMP
2880# 1353 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2881
2882# 1353 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2883
2884# 1353 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2885!$omp declare target device_type(any)
2886# 1353 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2887#endif
2888
2889# 1358 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2890 real(wp), dimension(num_fluids), intent(in) :: alpha_rho_k, alpha_k
2891# 1360 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2892 real(wp), intent(in) :: pres, vel_sum
2893 real(wp), intent(out) :: e
2894 real(wp) :: rho, gamma, pi_inf, qv
2895
2896 call s_compute_mixture_coefficients(alpha_rho_k, alpha_k, rho, gamma, pi_inf, qv)
2897
2898 ! E = dyn_p + (1 - alf)(gamma p + pi_inf) + qv. Only the liquid's internal energy is diluted;
2899 ! qv is already an energy density.
2900 e = gamma*pres + pi_inf
2901 if (bubbles_euler) e = e*(1._wp - alpha_k(num_fluids))
2902 e = e + qv + 5.e-1_wp*rho*vel_sum
2903
2904 end subroutine s_compute_energy
2905
2906 !> The reference curve of a state-dependent EOS at rho: p_ref, e_ref, their d/drho, and Gamma_G with its d/drho. A new family
2907 !! adds one case here and nothing else.
2908 subroutine s_reference_curve(rho, i, p_ref, e_ref, dp_drho, de_drho, G0, dG0)
2909
2910
2911# 1378 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2912#if MFC_OpenACC
2913# 1378 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2914!$acc routine seq
2915# 1378 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2916#elif MFC_OpenMP
2917# 1378 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2918
2919# 1378 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2920
2921# 1378 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2922!$omp declare target device_type(any)
2923# 1378 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2924#endif
2925
2926 real(wp), intent(in) :: rho
2927 integer, intent(in) :: i
2928 real(wp), intent(out) :: p_ref, e_ref, dp_drho, de_drho, G0, dG0
2929 real(wp) :: mu, d, V, ea, eb, up, us, dus, dup_dmu, x, ex, dp_dmu, de_dmu
2930 integer :: iter
2931
2932 mu = rho/eos_coeffs(i)%rho0 - 1._wp
2933 ! Past the fit's turnover there is no shock state to find; clamp rather than let the Newton below wander
2934 ! off and return a silently wrong pressure. mu_max is huge for the linear fit, so this is a no-op there.
2935 ! Bounded here so the step finishes and the host-side check in s_write_run_time_information can report it;
2936 ! past the turnover there is no shock state and the Newton below would wander.
2937 if (eoss(i) == eos_mie_gruneisen .and. mu > eos_coeffs(i)%mu_max) mu = eos_coeffs(i)%mu_max
2938 select case (eoss(i))
2939 case (eos_mie_gruneisen)
2940 ! Hugoniot reference u_s = c0 + s u_p + s2 u_p^2 + s3 u_p^3, with p_H = rho0 u_s u_p and the Hugoniot
2941 ! energy e_H = p_H mu/(2 rho0 (1 + mu)); linear on release. Pole at mu = 1/(s - 1) for the linear fit;
2942 ! the validator refuses initial states outside the EOS.
2943 if (mu < 0._wp) then
2944 p_ref = eos_coeffs(i)%rho0*eos_coeffs(i)%c0**2*mu
2945 dp_dmu = eos_coeffs(i)%rho0*eos_coeffs(i)%c0**2
2946 else if (eos_coeffs(i)%s2 == 0._wp .and. eos_coeffs(i)%s3 == 0._wp) then
2947 d = 1._wp - (eos_coeffs(i)%s - 1._wp)*mu
2948 p_ref = eos_coeffs(i)%rho0*eos_coeffs(i)%c0**2*mu*(1._wp + mu)/(d*d)
2949 dp_dmu = eos_coeffs(i)%rho0*eos_coeffs(i)%c0**2*((1._wp + 2._wp*mu)*d + 2._wp*(eos_coeffs(i)%s - 1._wp)*mu*(1._wp &
2950 & + mu))/(d*d*d)
2951 else
2952 ! u_p solves u_s(u_p) mu = u_p (1 + mu): Newton from the linear fit, then implicit differentiation
2953 up = eos_coeffs(i)%c0*mu/(1._wp - (eos_coeffs(i)%s - 1._wp)*mu)
2954
2955# 1408 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2956#if defined(MFC_OpenACC)
2957# 1408 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2958!$acc loop seq
2959# 1408 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2960#elif defined(MFC_OpenMP)
2961# 1408 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2962
2963# 1408 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
2964#endif
2965 do iter = 1, 8
2966 us = eos_coeffs(i)%c0 + up*(eos_coeffs(i)%s + up*(eos_coeffs(i)%s2 + up*eos_coeffs(i)%s3))
2967 dus = eos_coeffs(i)%s + up*(2._wp*eos_coeffs(i)%s2 + 3._wp*eos_coeffs(i)%s3*up)
2968 up = up - (us*mu - up*(1._wp + mu))/(dus*mu - (1._wp + mu))
2969 end do
2970 us = eos_coeffs(i)%c0 + up*(eos_coeffs(i)%s + up*(eos_coeffs(i)%s2 + up*eos_coeffs(i)%s3))
2971 dus = eos_coeffs(i)%s + up*(2._wp*eos_coeffs(i)%s2 + 3._wp*eos_coeffs(i)%s3*up)
2972 dup_dmu = (up - us)/(dus*mu - (1._wp + mu))
2973 p_ref = eos_coeffs(i)%rho0*us*up
2974 dp_dmu = eos_coeffs(i)%rho0*(dus*up + us)*dup_dmu
2975 end if
2976 e_ref = p_ref*mu/(2._wp*eos_coeffs(i)%rho0*(1._wp + mu))
2977 de_dmu = (dp_dmu*mu*(1._wp + mu) + p_ref)/(2._wp*eos_coeffs(i)%rho0*(1._wp + mu)**2)
2978 dp_drho = dp_dmu/eos_coeffs(i)%rho0
2979 de_drho = de_dmu/eos_coeffs(i)%rho0
2980 case (eos_jwl)
2981 ! JWL: p_ref = A exp(-R1 V) + B exp(-R2 V), V = rho0/rho. The curve is itself an isentrope, so de_ref = -p_ref d(1/rho).
2982 v = eos_coeffs(i)%rho0/rho
2983 ea = eos_coeffs(i)%a*exp(-eos_coeffs(i)%r1*v)
2984 eb = eos_coeffs(i)%b*exp(-eos_coeffs(i)%r2*v)
2985 p_ref = ea + eb
2986 e_ref = (ea/eos_coeffs(i)%r1 + eb/eos_coeffs(i)%r2)/eos_coeffs(i)%rho0
2987 dp_drho = (eos_coeffs(i)%rho0/rho**2)*(eos_coeffs(i)%r1*ea + eos_coeffs(i)%r2*eb)
2988 de_drho = p_ref/rho**2
2989 case (eos_vinet)
2990 ! Vinet cold curve: p_c = 3 K0 (1 - x)/x^2 exp(eta (1 - x)), x = (rho0/rho)^(1/3), eta = 3 (K0' - 1)/2,
2991 ! an isentrope like JWL (its energy integrates in closed form).
2992 d = 1.5_wp*(eos_coeffs(i)%k0p - 1._wp)
2993 x = (eos_coeffs(i)%rho0/rho)**(1._wp/3._wp)
2994 ex = exp(d*(1._wp - x))
2995 p_ref = 3._wp*eos_coeffs(i)%k0*(1._wp - x)/x**2*ex
2996 e_ref = 9._wp*eos_coeffs(i)%k0/(eos_coeffs(i)%rho0*d**2)*(1._wp - (1._wp - d*(1._wp - x))*ex)
2997 dp_drho = 3._wp*eos_coeffs(i)%k0*ex*(-1._wp/x**2 - 2._wp*(1._wp - x)/x**3 - d*(1._wp - x)/x**2)*(-x/(3._wp*rho))
2998 de_drho = p_ref/rho**2
2999 end select
3000 g0 = eos_coeffs(i)%gruneisen0 + eos_coeffs(i)%gruneisen_a*mu
3001 dg0 = eos_coeffs(i)%gruneisen_a/eos_coeffs(i)%rho0
3002
3003 end subroutine s_reference_curve
3004
3005 !> Whether the EOS of fluid i is a family whose coefficients vary with density.
3006 function f_is_state_dependent(i) result(yes)
3007
3008
3009# 1452 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3010#ifdef _CRAYFTN
3011# 1452 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3012#if MFC_OpenACC
3013# 1452 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3014!$acc routine seq
3015# 1452 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3016#elif MFC_OpenMP
3017# 1452 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3018
3019# 1452 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3020
3021# 1452 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3022!$omp declare target device_type(any)
3023# 1452 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3024#else
3025# 1452 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3026!DIR$ INLINEALWAYS f_is_state_dependent
3027# 1452 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3028#endif
3029# 1452 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3030#elif MFC_OpenACC
3031# 1452 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3032!$acc routine seq
3033# 1452 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3034#elif MFC_OpenMP
3035# 1452 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3036
3037# 1452 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3038
3039# 1452 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3040!$omp declare target device_type(any)
3041# 1452 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3042#endif
3043
3044 integer, intent(in) :: i
3045 logical :: yes
3046
3047 yes = eoss(i) == eos_mie_gruneisen .or. eoss(i) == eos_jwl .or. eoss(i) == eos_vinet
3048
3049 end function f_is_state_dependent
3050
3051 !> True when fluid i's reference curve is itself an isentrope (de_ref = -p_ref d(1/rho), which holds for JWL and Vinet but not
3052 !! for the Mie-Gruneisen Hugoniot) and its Gruneisen coefficient is constant. Those two together make the isentrope through any
3053 !! state closed-form, so it never has to be integrated.
3054 function f_has_isentropic_reference(i) result(yes)
3055
3056
3057# 1466 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3058#ifdef _CRAYFTN
3059# 1466 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3060#if MFC_OpenACC
3061# 1466 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3062!$acc routine seq
3063# 1466 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3064#elif MFC_OpenMP
3065# 1466 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3066
3067# 1466 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3068
3069# 1466 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3070!$omp declare target device_type(any)
3071# 1466 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3072#else
3073# 1466 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3074!DIR$ INLINEALWAYS f_has_isentropic_reference
3075# 1466 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3076#endif
3077# 1466 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3078#elif MFC_OpenACC
3079# 1466 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3080!$acc routine seq
3081# 1466 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3082#elif MFC_OpenMP
3083# 1466 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3084
3085# 1466 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3086
3087# 1466 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3088!$omp declare target device_type(any)
3089# 1466 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3090#endif
3091
3092 integer, intent(in) :: i
3093 logical :: yes
3094
3095 yes = (eoss(i) == eos_jwl .or. eoss(i) == eos_vinet) .and. eos_coeffs(i)%gruneisen_a == 0._wp
3096
3097 end function f_has_isentropic_reference
3098
3099 !> The largest compression a cubic Hugoniot fit can represent. mu(u_p) = u_p/(u_s - u_p) rises, peaks where c0 = s2 u_p^2 + 2 s3
3100 !! u_p^3, and falls after; only the rising branch is a physical shock. Returns a huge value for the linear fit, which never
3101 !! turns over. Host-side: called once per fluid at initialization.
3102 impure function f_hugoniot_compression_limit(c0, s, s2, s3) result(mu_max)
3103
3104 real(wp), intent(in) :: c0, s, s2, s3
3105 real(wp) :: mu_max, up, f, df, us
3106 integer :: iter
3107
3108 if (s2 == 0._wp .and. s3 == 0._wp) then
3109 mu_max = huge(1._wp)
3110 return
3111 end if
3112
3113 ! Newton on c0 - s2 u^2 - 2 s3 u^3 = 0, from a guess that brackets the physical range
3114 up = c0
3115 do iter = 1, 100
3116 f = c0 - s2*up**2 - 2._wp*s3*up**3
3117 df = -2._wp*s2*up - 6._wp*s3*up**2
3118 if (abs(df) < verysmall) exit
3119 up = max(up - f/df, verysmall)
3120 end do
3121 us = c0 + up*(s + up*(s2 + up*s3))
3122 mu_max = up/max(us - up, verysmall)
3123
3124 end function f_hugoniot_compression_limit
3125
3126 !> Gamma, Pi, dPi/drho and dGamma/drho of fluid i at density rho, the coefficients of rho e = Gamma p + Pi(rho). Stiffened and
3127 !! ideal gas keep the constants resolved at init, bit for bit.
3128 subroutine s_eos_coefficients(rho, i, gamma, pi_inf, dpi, dgamma)
3129
3130
3131# 1506 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3132#if MFC_OpenACC
3133# 1506 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3134!$acc routine seq
3135# 1506 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3136#elif MFC_OpenMP
3137# 1506 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3138
3139# 1506 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3140
3141# 1506 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3142!$omp declare target device_type(any)
3143# 1506 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3144#endif
3145
3146 real(wp), intent(in) :: rho
3147 integer, intent(in) :: i
3148 real(wp), intent(out) :: gamma, pi_inf, dpi, dgamma
3149 real(wp) :: p_ref, e_ref, dp_drho, de_drho, g0, dg0
3150
3151 if (.not. f_is_state_dependent(i)) then
3152 gamma = gammas(i)
3153 pi_inf = pi_infs(i)
3154 dpi = 0._wp
3155 dgamma = 0._wp
3156 return
3157 end if
3158 call s_reference_curve(rho, i, p_ref, e_ref, dp_drho, de_drho, g0, dg0)
3159 gamma = 1._wp/g0
3160 pi_inf = rho*e_ref - p_ref/g0
3161 dpi = e_ref + rho*de_drho - dp_drho/g0 + p_ref*dg0/g0**2
3162 dgamma = -dg0/g0**2
3163
3164 end subroutine s_eos_coefficients
3165
3166 !> Exponent of the stiffened-gas isentrope p + B = const rho**n. Precomputed per fluid as isentrope_n.
3167 function f_isentrope_exponent(gamma) result(n)
3168
3169
3170# 1531 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3171#ifdef _CRAYFTN
3172# 1531 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3173#if MFC_OpenACC
3174# 1531 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3175!$acc routine seq
3176# 1531 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3177#elif MFC_OpenMP
3178# 1531 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3179
3180# 1531 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3181
3182# 1531 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3183!$omp declare target device_type(any)
3184# 1531 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3185#else
3186# 1531 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3187!DIR$ INLINEALWAYS f_isentrope_exponent
3188# 1531 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3189#endif
3190# 1531 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3191#elif MFC_OpenACC
3192# 1531 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3193!$acc routine seq
3194# 1531 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3195#elif MFC_OpenMP
3196# 1531 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3197
3198# 1531 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3199
3200# 1531 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3201!$omp declare target device_type(any)
3202# 1531 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3203#endif
3204
3205 real(wp), intent(in) :: gamma
3206 real(wp) :: n
3207
3208 n = 1._wp/gamma + 1._wp
3209
3210 end function f_isentrope_exponent
3211
3212 !> Reference pressure of that isentrope. Precomputed per fluid as isentrope_B.
3213 function f_isentrope_pressure(pi_inf, gamma) result(B)
3214
3215
3216# 1543 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3217#ifdef _CRAYFTN
3218# 1543 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3219#if MFC_OpenACC
3220# 1543 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3221!$acc routine seq
3222# 1543 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3223#elif MFC_OpenMP
3224# 1543 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3225
3226# 1543 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3227
3228# 1543 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3229!$omp declare target device_type(any)
3230# 1543 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3231#else
3232# 1543 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3233!DIR$ INLINEALWAYS f_isentrope_pressure
3234# 1543 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3235#endif
3236# 1543 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3237#elif MFC_OpenACC
3238# 1543 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3239!$acc routine seq
3240# 1543 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3241#elif MFC_OpenMP
3242# 1543 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3243
3244# 1543 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3245
3246# 1543 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3247!$omp declare target device_type(any)
3248# 1543 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3249#endif
3250
3251 real(wp), intent(in) :: pi_inf, gamma
3252 real(wp) :: b
3253
3254 b = pi_inf/(1._wp + gamma)
3255
3256 end function f_isentrope_pressure
3257
3258 !> Stiffened-gas thermal law p + B = (n - 1)*cv*rho*T. Pass rho to get T, or T to get rho.
3259 function f_sg_thermal(pres, rho_or_T, n, B, cv) result(T_or_rho)
3260
3261
3262# 1555 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3263#ifdef _CRAYFTN
3264# 1555 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3265#if MFC_OpenACC
3266# 1555 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3267!$acc routine seq
3268# 1555 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3269#elif MFC_OpenMP
3270# 1555 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3271
3272# 1555 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3273
3274# 1555 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3275!$omp declare target device_type(any)
3276# 1555 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3277#else
3278# 1555 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3279!DIR$ INLINEALWAYS f_sg_thermal
3280# 1555 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3281#endif
3282# 1555 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3283#elif MFC_OpenACC
3284# 1555 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3285!$acc routine seq
3286# 1555 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3287#elif MFC_OpenMP
3288# 1555 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3289
3290# 1555 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3291
3292# 1555 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3293!$omp declare target device_type(any)
3294# 1555 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3295#endif
3296
3297 real(wp), intent(in) :: pres, rho_or_t, n, b, cv
3298 real(wp) :: t_or_rho
3299
3300 t_or_rho = (pres + b)/((n - 1._wp)*cv*rho_or_t)
3301
3302 end function f_sg_thermal
3303
3304 !> Thermal-equilibrium mixture temperature for stiffened gas, from primitives. Algebraically identical to the conservative form
3305 !! in m_phase_change's s_infinite_pt_relaxation_k, T = (rho*e + p - sum(alpha_rho_i*qv_i)) / sum(alpha_rho_i*cv_i*n_i), because
3306 !! rho*e = gamma_mix*p + pi_inf_mix + sum(alpha_rho_i*qv_i) in MFC's stored variables.
3307 function f_mixture_temperature(alpha_rho_K, pres, gamma_K, pi_inf_K) result(T)
3308
3309
3310# 1569 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3311#ifdef _CRAYFTN
3312# 1569 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3313#if MFC_OpenACC
3314# 1569 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3315!$acc routine seq
3316# 1569 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3317#elif MFC_OpenMP
3318# 1569 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3319
3320# 1569 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3321
3322# 1569 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3323!$omp declare target device_type(any)
3324# 1569 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3325#else
3326# 1569 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3327!DIR$ INLINEALWAYS f_mixture_temperature
3328# 1569 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3329#endif
3330# 1569 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3331#elif MFC_OpenACC
3332# 1569 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3333!$acc routine seq
3334# 1569 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3335#elif MFC_OpenMP
3336# 1569 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3337
3338# 1569 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3339
3340# 1569 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3341!$omp declare target device_type(any)
3342# 1569 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3343#endif
3344
3345# 1574 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3346 real(wp), dimension(num_fluids), intent(in) :: alpha_rho_k
3347# 1576 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3348 real(wp), intent(in) :: pres, gamma_k, pi_inf_k
3349 real(wp) :: t
3350 real(wp) :: mcp !< sum of alpha_rho_i*cp_i; cp_i = n_i*cv_i
3351 integer :: i
3352
3353 mcp = 0._wp
3354
3355# 1582 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3356#if defined(MFC_OpenACC)
3357# 1582 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3358!$acc loop seq
3359# 1582 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3360#elif defined(MFC_OpenMP)
3361# 1582 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3362
3363# 1582 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3364#endif
3365 do i = 1, num_fluids
3366 mcp = mcp + alpha_rho_k(i)*cvs(i)*isentrope_n(i)
3367 end do
3368
3369 t = ((gamma_k + 1._wp)*pres + pi_inf_k)/max(mcp, sgm_eps)
3370
3371 end function f_mixture_temperature
3372
3373 !> Coefficients of phase i at its own density alpha_rho/alpha: the per-cell dispatch when some fluid's EOS is state dependent,
3374 !! the constants resolved at init otherwise (bit for bit).
3375 subroutine s_phase_coefficients(alpha_rho, alpha, i, rho, gamma, pi_inf, dpi, dgamma)
3376
3377
3378# 1595 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3379#ifdef _CRAYFTN
3380# 1595 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3381#if MFC_OpenACC
3382# 1595 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3383!$acc routine seq
3384# 1595 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3385#elif MFC_OpenMP
3386# 1595 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3387
3388# 1595 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3389
3390# 1595 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3391!$omp declare target device_type(any)
3392# 1595 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3393#else
3394# 1595 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3395!DIR$ INLINEALWAYS s_phase_coefficients
3396# 1595 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3397#endif
3398# 1595 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3399#elif MFC_OpenACC
3400# 1595 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3401!$acc routine seq
3402# 1595 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3403#elif MFC_OpenMP
3404# 1595 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3405
3406# 1595 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3407
3408# 1595 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3409!$omp declare target device_type(any)
3410# 1595 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3411#endif
3412
3413 real(wp), intent(in) :: alpha_rho, alpha
3414 integer, intent(in) :: i
3415 real(wp), intent(out) :: rho, gamma, pi_inf, dpi, dgamma
3416
3417 rho = max(alpha_rho, sgm_eps)/max(alpha, sgm_eps)
3418 if (any_state_dependent_eos) then
3419 call s_eos_coefficients(rho, i, gamma, pi_inf, dpi, dgamma)
3420 else
3421 gamma = gammas(i)
3422 pi_inf = pi_infs(i)
3423 dpi = 0._wp
3424 dgamma = 0._wp
3425 end if
3426
3427 end subroutine s_phase_coefficients
3428
3429 !> c^2 = [((Gamma + 1) p + Pi)/rho - dPi/drho - p dGamma/drho]/Gamma, the frozen speed of one phase.
3430 function f_c2_from_coefficients(rho, pres, gamma, pi_inf, dpi, dgamma) result(c2)
3431
3432
3433# 1616 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3434#ifdef _CRAYFTN
3435# 1616 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3436#if MFC_OpenACC
3437# 1616 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3438!$acc routine seq
3439# 1616 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3440#elif MFC_OpenMP
3441# 1616 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3442
3443# 1616 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3444
3445# 1616 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3446!$omp declare target device_type(any)
3447# 1616 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3448#else
3449# 1616 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3450!DIR$ INLINEALWAYS f_c2_from_coefficients
3451# 1616 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3452#endif
3453# 1616 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3454#elif MFC_OpenACC
3455# 1616 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3456!$acc routine seq
3457# 1616 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3458#elif MFC_OpenMP
3459# 1616 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3460
3461# 1616 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3462
3463# 1616 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3464!$omp declare target device_type(any)
3465# 1616 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3466#endif
3467
3468 real(wp), intent(in) :: rho, pres, gamma, pi_inf, dpi, dgamma
3469 real(wp) :: c2
3470
3471 c2 = (((gamma + 1._wp)*pres + pi_inf)/rho - dpi - pres*dgamma)/gamma
3472
3473 end function f_c2_from_coefficients
3474
3475 !> Frozen sound speed squared of one phase at (rho, p) from its own coefficients. These helpers are subroutines, not functions:
3476 !! a device function that calls a device subroutine is a pattern no other backend-tested code in MFC uses.
3477 subroutine s_phase_c2(rho, pres, i, c2)
3478
3479
3480# 1629 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3481#if MFC_OpenACC
3482# 1629 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3483!$acc routine seq
3484# 1629 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3485#elif MFC_OpenMP
3486# 1629 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3487
3488# 1629 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3489
3490# 1629 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3491!$omp declare target device_type(any)
3492# 1629 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3493#endif
3494
3495 real(wp), intent(in) :: rho, pres
3496 integer, intent(in) :: i
3497 real(wp), intent(out) :: c2
3498 real(wp) :: gamma, pi_inf, dpi, dgamma
3499
3500 call s_eos_coefficients(rho, i, gamma, pi_inf, dpi, dgamma)
3501 c2 = f_c2_from_coefficients(rho, pres, gamma, pi_inf, dpi, dgamma)
3502
3503 end subroutine s_phase_c2
3504
3505 !> Slope of the ODE `kind` for fluid i: dp/drho = c^2 along an isentrope (x = rho, y = p), or the reference temperature dT/dV =
3506 !! (de_ref/dV + p_ref)/c_v - Gamma_G T/V (x = V, y = T), the Maxwell relation applied to e = e_ref + c_v (T - T_ref).
3507 subroutine s_ode_slope(kind, i, x, y, dydx)
3508
3509
3510# 1645 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3511#if MFC_OpenACC
3512# 1645 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3513!$acc routine seq
3514# 1645 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3515#elif MFC_OpenMP
3516# 1645 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3517
3518# 1645 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3519
3520# 1645 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3521!$omp declare target device_type(any)
3522# 1645 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3523#endif
3524
3525 integer, intent(in) :: kind, i
3526 real(wp), intent(in) :: x, y
3527 real(wp), intent(out) :: dydx
3528 real(wp) :: p_ref, e_ref, dp_drho, de_drho, G0, dG0
3529
3530 if (kind == ode_isentrope) then
3531 call s_phase_c2(x, y, i, dydx)
3532 else
3533 call s_reference_curve(1._wp/x, i, p_ref, e_ref, dp_drho, de_drho, g0, dg0)
3534 dydx = (p_ref - de_drho/x**2)/cvs(i) - g0*y/x
3535 end if
3536
3537 end subroutine s_ode_slope
3538
3539 !> Fixed-step classical RK4 for the ODE `kind` from (x0, y0) to x1.
3540 subroutine s_rk4(kind, i, x0, y0, x1, y)
3541
3542
3543# 1664 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3544#if MFC_OpenACC
3545# 1664 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3546!$acc routine seq
3547# 1664 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3548#elif MFC_OpenMP
3549# 1664 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3550
3551# 1664 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3552
3553# 1664 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3554!$omp declare target device_type(any)
3555# 1664 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3556#endif
3557
3558 integer, intent(in) :: kind, i
3559 real(wp), intent(in) :: x0, y0, x1
3560 real(wp), intent(out) :: y
3561 real(wp) :: x, h, k1, k2, k3, k4
3562 integer :: step
3563
3564 x = x0
3565 y = y0
3566 h = (x1 - x0)/eos_rk4_steps
3567
3568# 1675 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3569#if defined(MFC_OpenACC)
3570# 1675 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3571!$acc loop seq
3572# 1675 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3573#elif defined(MFC_OpenMP)
3574# 1675 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3575
3576# 1675 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3577#endif
3578 do step = 1, eos_rk4_steps
3579 call s_ode_slope(kind, i, x, y, k1)
3580 call s_ode_slope(kind, i, x + 0.5_wp*h, y + 0.5_wp*h*k1, k2)
3581 call s_ode_slope(kind, i, x + 0.5_wp*h, y + 0.5_wp*h*k2, k3)
3582 call s_ode_slope(kind, i, x + h, y + h*k3, k4)
3583 y = y + h*(k1 + 2._wp*(k2 + k3) + k4)/6._wp
3584 x = x + h
3585 end do
3586
3587 end subroutine s_rk4
3588
3589 !> Pressure of phase i after the isentropic density change rho -> xi rho: closed form for the constant-coefficient families,
3590 !! integrated for a state-dependent EOS (the star states it serves are close to rho).
3591 subroutine s_phase_pressure_on_isentrope(pres, rho, xi, i, p_isen)
3592
3593
3594# 1691 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3595#if MFC_OpenACC
3596# 1691 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3597!$acc routine seq
3598# 1691 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3599#elif MFC_OpenMP
3600# 1691 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3601
3602# 1691 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3603
3604# 1691 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3605!$omp declare target device_type(any)
3606# 1691 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3607#endif
3608
3609 real(wp), intent(in) :: pres, rho, xi
3610 integer, intent(in) :: i
3611 real(wp), intent(out) :: p_isen
3612 real(wp) :: p_ref_from, p_ref_to, e_ref, dp_drho, de_drho, g0, dg0
3613
3614 if (.not. f_is_state_dependent(i)) then
3615 p_isen = (pres + isentrope_b(i))*xi**isentrope_n(i) - isentrope_b(i)
3616 else if (f_has_isentropic_reference(i)) then
3617 ! Exact: the offset from an isentropic reference obeys dDelta/Delta = Gamma drho/rho, so
3618 ! p - p_ref scales as (rho'/rho)**(1 + Gamma). Integrating it instead costs a decimal per
3619 ! doubling of the expansion and turns the pressure negative past roughly twentyfold.
3620 call s_reference_curve(rho, i, p_ref_from, e_ref, dp_drho, de_drho, g0, dg0)
3621 call s_reference_curve(xi*rho, i, p_ref_to, e_ref, dp_drho, de_drho, g0, dg0)
3622 p_isen = p_ref_to + (pres - p_ref_from)*xi**(1._wp + eos_coeffs(i)%gruneisen0)
3623 else
3624 call s_rk4(ode_isentrope, i, rho, pres, xi*rho, p_isen)
3625 end if
3626
3627 end subroutine s_phase_pressure_on_isentrope
3628
3629 !> Temperature of phase i at (rho, p): the stiffened-gas relation, or T_ref(rho) + (e - e_ref)/c_v.
3630 subroutine s_phase_temperature(rho, pres, i, T)
3631
3632
3633# 1716 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3634#if MFC_OpenACC
3635# 1716 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3636!$acc routine seq
3637# 1716 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3638#elif MFC_OpenMP
3639# 1716 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3640
3641# 1716 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3642
3643# 1716 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3644!$omp declare target device_type(any)
3645# 1716 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3646#endif
3647
3648 real(wp), intent(in) :: rho, pres
3649 integer, intent(in) :: i
3650 real(wp), intent(out) :: t
3651 real(wp) :: p_ref, e_ref, dp_drho, de_drho, g0, dg0, t0, t_ref
3652
3653 if (f_is_state_dependent(i)) then
3654 call s_reference_curve(rho, i, p_ref, e_ref, dp_drho, de_drho, g0, dg0)
3655 t0 = eos_coeffs(i)%t0
3656 call s_rk4(ode_reference_temperature, i, 1._wp/eos_coeffs(i)%rho0, t0, 1._wp/rho, t_ref)
3657 t = t_ref + (pres - p_ref)/(rho*g0*cvs(i))
3658 else
3659 t = (pres + isentrope_b(i))/((isentrope_n(i) - 1._wp)*cvs(i)*rho)
3660 end if
3661
3662 end subroutine s_phase_temperature
3663
3664 !> Density of phase i on the isentrope through (rho_from, p_from) at p_to, and c^2 there: Newton on the pressure integrator,
3665 !! whose slope is c^2. The relaxation's own Newton wraps this, so a few steps suffice.
3666 subroutine s_phase_density_on_isentrope(i, rho_from, p_from, p_to, rho_to, c2_to)
3667
3668
3669# 1738 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3670#if MFC_OpenACC
3671# 1738 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3672!$acc routine seq
3673# 1738 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3674#elif MFC_OpenMP
3675# 1738 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3676
3677# 1738 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3678
3679# 1738 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3680!$omp declare target device_type(any)
3681# 1738 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3682#endif
3683
3684 integer, intent(in) :: i
3685 real(wp), intent(in) :: rho_from, p_from, p_to
3686 real(wp), intent(out) :: rho_to, c2_to
3687 real(wp) :: p_at, c2_at
3688 integer :: iter
3689
3690 rho_to = rho_from
3691
3692# 1747 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3693#if defined(MFC_OpenACC)
3694# 1747 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3695!$acc loop seq
3696# 1747 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3697#elif defined(MFC_OpenMP)
3698# 1747 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3699
3700# 1747 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3701#endif
3702 do iter = 1, 4
3703 call s_phase_pressure_on_isentrope(p_from, rho_from, rho_to/rho_from, i, p_at)
3704 call s_phase_c2(rho_to, p_at, i, c2_at)
3705 rho_to = rho_to - (p_at - p_to)/c2_at
3706 end do
3707 call s_phase_c2(rho_to, p_to, i, c2_to)
3708
3709 end subroutine s_phase_density_on_isentrope
3710
3711 !> Internal energy per unit volume of phase i at pressure pres: alpha (Gamma p + Pi) + alpha_rho qv, with the coefficients at
3712 !! the phase's own density.
3713 subroutine s_phase_internal_energy(pres, alpha, alpha_rho, i, e_phase)
3714
3715
3716# 1761 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3717#ifdef _CRAYFTN
3718# 1761 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3719#if MFC_OpenACC
3720# 1761 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3721!$acc routine seq
3722# 1761 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3723#elif MFC_OpenMP
3724# 1761 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3725
3726# 1761 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3727
3728# 1761 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3729!$omp declare target device_type(any)
3730# 1761 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3731#else
3732# 1761 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3733!DIR$ INLINEALWAYS s_phase_internal_energy
3734# 1761 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3735#endif
3736# 1761 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3737#elif MFC_OpenACC
3738# 1761 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3739!$acc routine seq
3740# 1761 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3741#elif MFC_OpenMP
3742# 1761 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3743
3744# 1761 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3745
3746# 1761 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3747!$omp declare target device_type(any)
3748# 1761 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3749#endif
3750
3751 real(wp), intent(in) :: pres, alpha, alpha_rho
3752 integer, intent(in) :: i
3753 real(wp), intent(out) :: e_phase
3754 real(wp) :: rho, gamma, pi_inf, dpi, dgamma
3755
3756 call s_phase_coefficients(alpha_rho, alpha, i, rho, gamma, pi_inf, dpi, dgamma)
3757 e_phase = alpha*(gamma*pres + pi_inf) + alpha_rho*qvs(i)
3758
3759 end subroutine s_phase_internal_energy
3760
3761 !> Bulk modulus rho c^2 of phase i at pressure pres: f_bulk_modulus for a constant-coefficient fluid, bit for bit, minus the
3762 !! reference-curve terms rho (dPi/drho + p dGamma/drho)/Gamma otherwise.
3763 subroutine s_phase_bulk_modulus(pres, alpha, alpha_rho, i, blkmod)
3764
3765
3766# 1777 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3767#ifdef _CRAYFTN
3768# 1777 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3769#if MFC_OpenACC
3770# 1777 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3771!$acc routine seq
3772# 1777 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3773#elif MFC_OpenMP
3774# 1777 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3775
3776# 1777 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3777
3778# 1777 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3779!$omp declare target device_type(any)
3780# 1777 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3781#else
3782# 1777 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3783!DIR$ INLINEALWAYS s_phase_bulk_modulus
3784# 1777 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3785#endif
3786# 1777 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3787#elif MFC_OpenACC
3788# 1777 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3789!$acc routine seq
3790# 1777 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3791#elif MFC_OpenMP
3792# 1777 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3793
3794# 1777 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3795
3796# 1777 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3797!$omp declare target device_type(any)
3798# 1777 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3799#endif
3800
3801 real(wp), intent(in) :: alpha_rho, alpha, pres
3802 integer, intent(in) :: i
3803 real(wp), intent(out) :: blkmod
3804 real(wp) :: rho, gamma, pi_inf, dpi, dgamma
3805
3806 call s_phase_coefficients(alpha_rho, alpha, i, rho, gamma, pi_inf, dpi, dgamma)
3807 blkmod = f_bulk_modulus(pres, gamma, pi_inf) - rho*(dpi + pres*dgamma)/gamma
3808
3809 end subroutine s_phase_bulk_modulus
3810
3811 !> Elastic strain energy of one stress component, doubled for a shear component: the tensor stores it once, the energy counts
3812 !! both off-diagonal entries. Zero without a shear modulus.
3813 function f_elastic_energy(tau, G, is_shear) result(dE)
3814
3815
3816# 1793 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3817#ifdef _CRAYFTN
3818# 1793 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3819#if MFC_OpenACC
3820# 1793 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3821!$acc routine seq
3822# 1793 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3823#elif MFC_OpenMP
3824# 1793 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3825
3826# 1793 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3827
3828# 1793 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3829!$omp declare target device_type(any)
3830# 1793 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3831#else
3832# 1793 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3833!DIR$ INLINEALWAYS f_elastic_energy
3834# 1793 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3835#endif
3836# 1793 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3837#elif MFC_OpenACC
3838# 1793 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3839!$acc routine seq
3840# 1793 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3841#elif MFC_OpenMP
3842# 1793 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3843
3844# 1793 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3845
3846# 1793 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3847!$omp declare target device_type(any)
3848# 1793 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3849#endif
3850
3851 real(wp), intent(in) :: tau, g
3852 logical, intent(in) :: is_shear
3853 real(wp) :: de
3854
3855 de = 0._wp
3856 if (g > verysmall) then
3857 de = (tau*tau)/max(4._wp*g, verysmall)
3858 if (is_shear) de = de + (tau*tau)/max(4._wp*g, verysmall)
3859 end if
3860
3861 end function f_elastic_energy
3862
3863 !> Hypoelastic strain energy at one cell, summed over the stress components.
3864 function f_hypoelastic_energy(q_cons_vf, j, k, l, rho, G) result(E_e)
3865
3866 type(scalar_field), dimension(sys_size), intent(in) :: q_cons_vf
3867 integer, intent(in) :: j, k, l
3868 real(wp), intent(in) :: rho, g
3869 real(wp) :: e_e
3870 integer :: s
3871
3872 e_e = 0._wp
3873 do s = eqn_idx%stress%beg, eqn_idx%stress%end
3874 e_e = e_e + f_elastic_energy(real(q_cons_vf(s)%sf(j, k, l), wp)/rho, g, any(s == shear_indices))
3875 end do
3876
3877 end function f_hypoelastic_energy
3878
3879 !> Pressure of a stiffened gas from its internal energy density - the inverse of s_compute_energy. Callers subtract the kinetic,
3880 !! magnetic and elastic energy first; none of those are equation-of-state terms.
3881 function f_pressure(e_int, gamma, pi_inf, qv) result(pres)
3882
3883
3884# 1827 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3885#ifdef _CRAYFTN
3886# 1827 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3887#if MFC_OpenACC
3888# 1827 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3889!$acc routine seq
3890# 1827 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3891#elif MFC_OpenMP
3892# 1827 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3893
3894# 1827 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3895
3896# 1827 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3897!$omp declare target device_type(any)
3898# 1827 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3899#else
3900# 1827 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3901!DIR$ INLINEALWAYS f_pressure
3902# 1827 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3903#endif
3904# 1827 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3905#elif MFC_OpenACC
3906# 1827 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3907!$acc routine seq
3908# 1827 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3909#elif MFC_OpenMP
3910# 1827 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3911
3912# 1827 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3913
3914# 1827 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3915!$omp declare target device_type(any)
3916# 1827 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3917#endif
3918
3919 real(wp), intent(in) :: e_int, gamma, pi_inf, qv
3920 real(wp) :: pres
3921
3922 pres = (e_int - pi_inf - qv)/gamma
3923
3924 end function f_pressure
3925
3926 !> Isentropic bulk modulus. Takes coefficients rather than a fluid index, so a mixture - whose effective gamma and pi_inf come
3927 !! from s_compute_mixture_coefficients - is the same call as a single fluid. Elastic callers add their own shear term.
3928 function f_bulk_modulus(pres, gamma, pi_inf) result(blkmod)
3929
3930
3931# 1840 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3932#ifdef _CRAYFTN
3933# 1840 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3934#if MFC_OpenACC
3935# 1840 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3936!$acc routine seq
3937# 1840 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3938#elif MFC_OpenMP
3939# 1840 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3940
3941# 1840 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3942
3943# 1840 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3944!$omp declare target device_type(any)
3945# 1840 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3946#else
3947# 1840 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3948!DIR$ INLINEALWAYS f_bulk_modulus
3949# 1840 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3950#endif
3951# 1840 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3952#elif MFC_OpenACC
3953# 1840 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3954!$acc routine seq
3955# 1840 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3956#elif MFC_OpenMP
3957# 1840 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3958
3959# 1840 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3960
3961# 1840 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3962!$omp declare target device_type(any)
3963# 1840 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3964#endif
3965
3966 real(wp), intent(in) :: pres, gamma, pi_inf
3967 real(wp) :: blkmod
3968
3969 blkmod = ((gamma + 1._wp)*pres + pi_inf)/gamma
3970
3971 end function f_bulk_modulus
3972
3973 !> Relativistic specific enthalpy, h = 1 + (Gamma + 1)p/rho. Ideal gas only: the stiffness does not appear, so a fluid with a
3974 !! nonzero pi_inf is not represented here (the validator refuses that combination).
3975 function f_relativistic_enthalpy(pres, rho, gamma) result(H)
3976
3977
3978# 1853 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3979#ifdef _CRAYFTN
3980# 1853 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3981#if MFC_OpenACC
3982# 1853 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3983!$acc routine seq
3984# 1853 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3985#elif MFC_OpenMP
3986# 1853 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3987
3988# 1853 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3989
3990# 1853 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3991!$omp declare target device_type(any)
3992# 1853 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3993#else
3994# 1853 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3995!DIR$ INLINEALWAYS f_relativistic_enthalpy
3996# 1853 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3997#endif
3998# 1853 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
3999#elif MFC_OpenACC
4000# 1853 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4001!$acc routine seq
4002# 1853 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4003#elif MFC_OpenMP
4004# 1853 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4005
4006# 1853 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4007
4008# 1853 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4009!$omp declare target device_type(any)
4010# 1853 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4011#endif
4012
4013 real(wp), intent(in) :: pres, rho, gamma
4014 real(wp) :: h
4015
4016 h = 1._wp + (gamma + 1._wp)*pres/rho
4017
4018 end function f_relativistic_enthalpy
4019
4020 !> Speed of sound of a thermodynamic state. Enthalpy is not an argument: for a real state H, |u|^2 and qv all cancel out of c^2
4021 !! = ((Gamma + 1)p + Pi)/(Gamma rho). Averaged states, whose enthalpy is a free input, use the _avg variant.
4022 subroutine s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c, alpha_rho)
4023
4024
4025# 1866 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4026#if MFC_OpenACC
4027# 1866 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4028!$acc routine seq
4029# 1866 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4030#elif MFC_OpenMP
4031# 1866 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4032
4033# 1866 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4034
4035# 1866 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4036!$omp declare target device_type(any)
4037# 1866 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4038#endif
4039
4040 real(wp), intent(in) :: pres, rho, gamma, pi_inf
4041# 1872 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4042 real(wp), dimension(num_fluids), intent(in) :: adv
4043# 1874 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4044 real(wp), intent(out) :: c
4045# 1878 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4046 real(wp), dimension(num_fluids), intent(in), optional :: alpha_rho
4047# 1880 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4048 real(wp) :: alf !< Subgrid void fraction; dilute by construction
4049 real(wp) :: blkmod_q, alpha_q, alpha_rho_q, gamma_q, pi_inf_q
4050 integer :: q
4051
4052 if (chemistry) then ! Reacting mixture sound speed
4053 c = sqrt((1.0_wp + 1.0_wp/gamma)*pres/rho)
4054 else if (relativity) then ! Relativistic sound speed, whose enthalpy is 1 + (Gamma + 1)p/rho
4055 c = sqrt((1._wp + 1._wp/gamma)*pres/rho/f_relativistic_enthalpy(pres, rho, gamma))
4056 else
4057 ! Every case below is a bulk modulus over a density. The equation of state enters
4058 ! only through f_bulk_modulus; the cases differ in how the phases are mixed.
4059 if (any_state_dependent_eos .and. present(alpha_rho)) then ! frozen mixing: each phase's modulus at its own density
4060 c = 0._wp
4061
4062# 1893 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4063#if defined(MFC_OpenACC)
4064# 1893 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4065!$acc loop seq
4066# 1893 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4067#elif defined(MFC_OpenMP)
4068# 1893 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4069
4070# 1893 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4071#endif
4072 do q = 1, num_fluids
4073 alpha_q = adv(q)
4074 alpha_rho_q = alpha_rho(q)
4075 call s_phase_bulk_modulus(pres, alpha_q, alpha_rho_q, q, blkmod_q)
4076 if (alt_soundspeed) then
4077 c = c + adv(q)/blkmod_q
4078 else
4079 c = c + adv(q)*blkmod_q
4080 end if
4081 end do
4082 if (alt_soundspeed) then
4083 c = 1._wp/(rho*c)
4084 else
4085 c = c/rho
4086 end if
4087 else if (alt_soundspeed) then ! Wood's law: volume-weighted harmonic mean
4088 c = 0._wp
4089
4090# 1911 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4091#if defined(MFC_OpenACC)
4092# 1911 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4093!$acc loop seq
4094# 1911 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4095#elif defined(MFC_OpenMP)
4096# 1911 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4097
4098# 1911 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4099#endif
4100 do q = 1, num_fluids
4101 gamma_q = gammas(q)
4102 pi_inf_q = pi_infs(q)
4103 c = c + adv(q)/f_bulk_modulus(pres, gamma_q, pi_inf_q)
4104 end do
4105 c = 1._wp/(rho*c)
4106 else if (model_eqns == model_eqns_6eq) then ! volume-weighted arithmetic mean
4107 c = 0._wp
4108
4109# 1920 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4110#if defined(MFC_OpenACC)
4111# 1920 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4112!$acc loop seq
4113# 1920 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4114#elif defined(MFC_OpenMP)
4115# 1920 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4116
4117# 1920 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4118#endif
4119 do q = 1, num_fluids
4120 gamma_q = gammas(q)
4121 pi_inf_q = pi_infs(q)
4122 c = c + adv(q)*f_bulk_modulus(pres, gamma_q, pi_inf_q)
4123 end do
4124 c = c/rho
4125 else ! the mixture coefficients already carry the mixing
4126 c = f_bulk_modulus(pres, gamma, pi_inf)/rho
4127
4128 ! Subgrid bubbles: c = c_l/(1 - alf), the carrier-liquid speed with an O(alf) void
4129 ! correction. alf is dilute by construction; near one means a wrong index or an
4130 ! out-of-regime case, which the toolchain warns about at case load (#1793).
4131 if (model_eqns == model_eqns_5eq .and. bubbles_euler .and. .not. (mpp_lim .and. num_fluids > 1)) then
4132 alf = adv(num_fluids)
4133 c = c/(1._wp - alf)
4134 end if
4135 end if
4136
4137 if (mixture_err .and. c < 0._wp) then
4138 c = 100._wp*sgm_eps
4139 else
4140 c = sqrt(c)
4141 end if
4142 end if
4143
4144 end subroutine s_compute_speed_of_sound
4145
4146 !> Speed of sound of an interface-averaged state. An average of two states is not a state - its enthalpy is not the one its
4147 !! pressure and density imply - so the caller supplies H, |u|^2 and qv. Only the enthalpy-reading branches differ from
4148 !! s_compute_speed_of_sound; keep the condition below in step with the branch list there.
4149 subroutine s_compute_speed_of_sound_avg(pres, rho, gamma, pi_inf, qv, vel_sum, H, c_c, adv, c, alpha_rho)
4150
4151
4152# 1953 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4153#if MFC_OpenACC
4154# 1953 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4155!$acc routine seq
4156# 1953 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4157#elif MFC_OpenMP
4158# 1953 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4159
4160# 1953 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4161
4162# 1953 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4163!$omp declare target device_type(any)
4164# 1953 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4165#endif
4166
4167 real(wp), intent(in) :: pres, rho, gamma, pi_inf, qv, vel_sum, h, c_c
4168# 1959 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4169 real(wp), dimension(num_fluids), intent(in) :: adv
4170# 1961 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4171 real(wp), intent(out) :: c
4172# 1965 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4173 real(wp), dimension(num_fluids), intent(in), optional :: alpha_rho
4174# 1967 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4175
4176 if (chemistry) then ! Reacting mixture sound speed
4177 if (avg_state == avg_state_roe .and. abs(c_c) > verysmall) then
4178 c = sqrt(c_c - (gamma - 1.0_wp)*(vel_sum - h))
4179 else
4180 call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c, alpha_rho)
4181 end if
4182 else if (relativity) then ! Relativistic sound speed
4183 c = sqrt((1._wp + 1._wp/gamma)*pres/rho/h)
4184 else if (alt_soundspeed .or. model_eqns == model_eqns_6eq .or. (model_eqns == model_eqns_5eq .and. bubbles_euler) &
4185 & .or. any_state_dependent_eos) then
4186 call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c, alpha_rho)
4187 else ! Stiffened-gas mixture, the one branch where the averaged enthalpy survives
4188 c = (h - 5.e-1*vel_sum - qv/rho)/gamma
4189
4190 if (mixture_err .and. c < 0._wp) then
4191 c = 100._wp*sgm_eps
4192 else
4193 c = sqrt(c)
4194 end if
4195 end if
4196
4197 end subroutine s_compute_speed_of_sound_avg
4198
4199 !> Compute the fast magnetosonic wave speed from the sound speed, density, and magnetic field components.
4200 subroutine s_compute_fast_magnetosonic_speed(rho, c, B, norm, c_fast, h)
4201
4202
4203# 1994 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4204#ifdef _CRAYFTN
4205# 1994 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4206#if MFC_OpenACC
4207# 1994 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4208!$acc routine seq
4209# 1994 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4210#elif MFC_OpenMP
4211# 1994 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4212
4213# 1994 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4214
4215# 1994 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4216!$omp declare target device_type(any)
4217# 1994 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4218#else
4219# 1994 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4220!DIR$ NOINLINE s_compute_fast_magnetosonic_speed
4221# 1994 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4222#endif
4223# 1994 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4224#elif MFC_OpenACC
4225# 1994 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4226!$acc routine seq
4227# 1994 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4228#elif MFC_OpenMP
4229# 1994 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4230
4231# 1994 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4232
4233# 1994 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4234!$omp declare target device_type(any)
4235# 1994 "/home/runner/work/MFC/MFC/src/common/m_variables_conversion.fpp"
4236#endif
4237
4238 real(wp), intent(in) :: b(3), rho, c
4239 real(wp), intent(in) :: h !< only used for relativity
4240 real(wp), intent(out) :: c_fast
4241 integer, intent(in) :: norm
4242 real(wp) :: b2, term, disc
4243
4244 b2 = sum(b**2)
4245
4246 if (.not. relativity) then
4247 term = c**2 + b2/rho
4248 disc = term**2 - 4*c**2*(b(norm)**2/rho)
4249 else
4250 ! Note: this is approximation for the non-relatisitic limit; accurate solution requires solving a quartic equation
4251 term = (c**2*(b(norm)**2 + rho*h) + b2)/(rho*h + b2)
4252 disc = term**2 - 4*c**2*b(norm)**2/(rho*h + b2)
4253 end if
4254
4255#ifdef MFC_DEBUG
4256 if (disc < 0._wp) then
4257 print *, 'rho, c, Bx, By, Bz, h, term, disc:', rho, c, b(1), b(2), b(3), h, term, disc
4258 ! s_mpi_abort is a host routine and cannot be called from device code
4259 ! (this is a GPU routine); on GPU builds, emit the diagnostic print only.
4260#ifndef MFC_GPU
4261 call s_mpi_abort('Error: negative discriminant in s_compute_fast_magnetosonic_speed')
4262#endif
4263 end if
4264#endif
4265
4266 c_fast = sqrt(0.5_wp*(term + sqrt(disc)))
4267
4269
4270end module m_variables_conversion
type(scalar_field), dimension(sys_size), intent(inout) q_cons_vf
integer, intent(in) k
integer, intent(in) j
integer, intent(in) l
Compile-time constant parameters: default values, tolerances, and physical constants.
integer, parameter model_eqns_5eq
integer, parameter avg_state_roe
integer, parameter riemann_solver_hll
integer, parameter riemann_solver_hlld
real(wp), parameter sgm_eps
Segmentation tolerance.
real(wp), parameter dflt_real
Default real value.
integer, parameter model_eqns_6eq
integer, parameter model_eqns_gamma_law
Shared derived types for field data, patch geometry, bubble dynamics, and MPI I/O structures.
Shared global parameters and equation-index setup for all three executables. Each per-target m_global...
type(physical_parameters), dimension(num_fluids_max) fluid_pp
Per-fluid stiffened-gas EOS parameters, Reynolds numbers, and shear modulus.
logical heat_conduction
any_state_dependent_eos is declared with the case-optimization block above: a parameter when the case...
type(eqn_idx_info) eqn_idx
All conserved-variable equation index ranges and scalars.
integer, dimension(3) shear_indices
Indices of the stress components that represent shear stress.
Defines global parameters for the computational domain, simulation algorithm, and initial conditions.
Basic floating-point utilities: approximate equality, default detection, and coordinate bounds.
Utility routines for bubble model setup, coordinate transforms, array sampling, and special functions...
Broadcasts user inputs and decomposes the domain across MPI ranks for pre-processing.
Conservative-to-primitive variable conversion, mixture property evaluation, and pressure computation.
real(wp) function, public f_bulk_modulus(pres, gamma, pi_inf)
Isentropic bulk modulus. Takes coefficients rather than a fluid index, so a mixture - whose effective...
subroutine, public s_compute_speed_of_sound_avg(pres, rho, gamma, pi_inf, qv, vel_sum, h, c_c, adv, c, alpha_rho)
Speed of sound of an interface-averaged state. An average of two states is not a state - its enthalpy...
real(wp) function, public f_elastic_energy(tau, g, is_shear)
Elastic strain energy of one stress component, doubled for a shear component: the tensor stores it on...
subroutine, public s_phase_temperature(rho, pres, i, t)
Temperature of phase i at (rho, p): the stiffened-gas relation, or T_ref(rho) + (e - e_ref)/c_v.
subroutine, public s_compute_species_fraction(q_vf, k, l, r, alpha_rho_k, alpha_k)
Compute partial densities and volume fractions.
subroutine s_reference_curve(rho, i, p_ref, e_ref, dp_drho, de_drho, g0, dg0)
The reference curve of a state-dependent EOS at rho: p_ref, e_ref, their d/drho, and Gamma_G with its...
subroutine, public s_phase_internal_energy(pres, alpha, alpha_rho, i, e_phase)
Internal energy per unit volume of phase i at pressure pres: alpha (Gamma p + Pi) + alpha_rho qv,...
real(wp), dimension(:,:), allocatable res_vc
impure subroutine, public s_finalize_variables_conversion_module()
Deallocate fluid property arrays and post-processing fields allocated during module initialization.
subroutine, public s_compute_mixture_coefficients(alpha_rho_k, alpha_k, rho_k, gamma_k, pi_inf_k, qv_k)
Mixture coefficients of one state. Under bubbles_euler with num_fluids == 1 the sole advection slot a...
real(wp) function, public f_hypoelastic_energy(q_cons_vf, j, k, l, rho, g)
Hypoelastic strain energy at one cell, summed over the stress components.
subroutine, public s_initialize_mv(qk_cons_vf, mv)
Initialize bubble mass-vapor values at quadrature nodes from the conserved moment statistics.
subroutine, public s_initialize_pb(qk_cons_vf, mv, pb)
Initialize bubble internal pressures at quadrature nodes using isothermal relations from the Preston ...
impure subroutine, public s_convert_primitive_to_conservative_variables(q_prim_vf, q_cons_vf)
Convert primitives (rho, u, p, alpha) to conserved variables (rho*alpha, rho*u, E,...
subroutine s_ode_slope(kind, i, x, y, dydx)
Slope of the ODE kind for fluid i: dp/drho = c^2 along an isentrope (x = rho, y = p),...
real(wp) function, public f_sg_thermal(pres, rho_or_t, n, b, cv)
Stiffened-gas thermal law p + B = (n - 1)*cv*rho*T. Pass rho to get T, or T to get rho.
subroutine s_phase_c2(rho, pres, i, c2)
Frozen sound speed squared of one phase at (rho, p) from its own coefficients. These helpers are subr...
subroutine, public s_convert_conservative_to_primitive_variables(qk_cons_vf, q_t_sf, qk_prim_vf, ibounds)
Convert conserved variables (rho*alpha, rho*u, E, alpha) to primitives (rho, u, p,...
real(wp) function, public f_relativistic_enthalpy(pres, rho, gamma)
Relativistic specific enthalpy, h = 1 + (Gamma + 1)p/rho. Ideal gas only: the stiffness does not appe...
subroutine, public s_convert_species_to_mixture_variables(q_vf, k, l, r, rho, gamma, pi_inf, qv, re_k, g_k, g)
Convert species volume fractions and partial densities to mixture density, gamma, pi_inf,...
subroutine, public s_compute_pressure(energy, alf, dyn_p, pi_inf, gamma, rho, qv, rhoyks, pres, t, e_e_in, pres_mag)
Compute the pressure from the appropriate equation of state.
subroutine, public s_convert_mixture_to_mixture_variables(q_vf, i, j, k, rho, gamma, pi_inf, qv)
Convert mixture variables to density, gamma, pi_inf, and qv for the gamma/pi_inf model....
real(wp), dimension(:,:,:), allocatable, public pi_inf_sf
Scalar liquid stiffness function.
logical function f_has_isentropic_reference(i)
True when fluid i's reference curve is itself an isentrope (de_ref = -p_ref d(1/rho),...
subroutine, public s_compute_energy(pres, alpha_rho_k, alpha_k, vel_sum, e)
Total energy per unit volume, thermodynamic terms only. Callers add magnetic and elastic energy,...
subroutine, public s_phase_bulk_modulus(pres, alpha, alpha_rho, i, blkmod)
Bulk modulus rho c^2 of phase i at pressure pres: f_bulk_modulus for a constant-coefficient fluid,...
subroutine, public s_compute_mixture_coefficients_dt(dalpha_rho_dt, dadv_dt, alpha_rho, adv, drho_dt, dgamma_dt, dpi_inf_dt, dqv_dt)
Time derivative of the mixture coefficients, mirroring s_compute_mixture_coefficients.
subroutine, public s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c, alpha_rho)
Speed of sound of a thermodynamic state. Enthalpy is not an argument: for a real state H,...
subroutine, public s_convert_primitive_to_flux_variables(qk_prim_vf, fk_vf, fk_src_vf, is1, is2, is3, s2b, s3b, dir_idx_in, dir_flg_in, hll_u_interface_in)
Convert primitive variables to Eulerian flux variables.
subroutine, public s_compute_fast_magnetosonic_speed(rho, c, b, norm, c_fast, h)
Compute the fast magnetosonic wave speed from the sound speed, density, and magnetic field components...
real(wp), dimension(:,:,:), allocatable, public gamma_sf
Scalar sp. heat ratio function.
real(wp) function, public f_isentrope_exponent(gamma)
Exponent of the stiffened-gas isentrope p + B = const rho**n. Precomputed per fluid as isentrope_n.
subroutine, public s_phase_density_on_isentrope(i, rho_from, p_from, p_to, rho_to, c2_to)
Density of phase i on the isentrope through (rho_from, p_from) at p_to, and c^2 there: Newton on the ...
impure subroutine, public s_initialize_variables_conversion_module(store_mixture_fields, enforce_density_floor, preserve_qbmm_number, lagrange_beta_index)
Initialize the variables conversion module.
impure real(wp) function f_hugoniot_compression_limit(c0, s, s2, s3)
The largest compression a cubic Hugoniot fit can represent. mu(u_p) = u_p/(u_s - u_p) rises,...
real(wp), dimension(:,:,:), allocatable, public rho_sf
Scalar density function.
real(wp) function, public f_isentrope_pressure(pi_inf, gamma)
Reference pressure of that isentrope. Precomputed per fluid as isentrope_B.
integer, dimension(:), allocatable bubrs_vc
real(wp) function, public f_mixture_temperature(alpha_rho_k, pres, gamma_k, pi_inf_k)
Thermal-equilibrium mixture temperature for stiffened gas, from primitives. Algebraically identical t...
real(wp), dimension(:), allocatable gs_vc
subroutine s_rk4(kind, i, x0, y0, x1, y)
Fixed-step classical RK4 for the ODE kind from (x0, y0) to x1.
subroutine, public s_convert_species_to_mixture_variables_kernel(rho_k, gamma_k, pi_inf_k, qv_k, alpha_k, alpha_rho_k, re_k, g_k, g)
Host- and device-callable conversion kernel for species and mixture variables.
logical function, public f_is_state_dependent(i)
Whether the EOS of fluid i is a family whose coefficients vary with density.
real(wp) function f_c2_from_coefficients(rho, pres, gamma, pi_inf, dpi, dgamma)
c^2 = [((Gamma + 1) p + Pi)/rho - dPi/drho - p dGamma/drho]/Gamma, the frozen speed of one phase.
subroutine, public s_convert_to_mixture_variables(q_vf, i, j, k, rho, gamma, pi_inf, qv, re_k, g_k, g)
Dispatch to the s_convert_mixture_to_mixture_variables and s_convert_species_to_mixture_variables sub...
subroutine, public s_eos_coefficients(rho, i, gamma, pi_inf, dpi, dgamma)
Gamma, Pi, dPi/drho and dGamma/drho of fluid i at density rho, the coefficients of rho e = Gamma p + ...
subroutine, public s_phase_pressure_on_isentrope(pres, rho, xi, i, p_isen)
Pressure of phase i after the isentropic density change rho -> xi rho: closed form for the constant-c...
real(wp) function, public f_pressure(e_int, gamma, pi_inf, qv)
Pressure of a stiffened gas from its internal energy density - the inverse of s_compute_energy....
subroutine, public s_phase_coefficients(alpha_rho, alpha, i, rho, gamma, pi_inf, dpi, dgamma)
Coefficients of phase i at its own density alpha_rho/alpha: the per-cell dispatch when some fluid's E...
Derived type annexing a scalar field (SF).