MFC
Exascale flow solver
Loading...
Searching...
No Matches
m_eos.fpp.f90
Go to the documentation of this file.
1# 1 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2!>
3!! @file
4!! @brief Contains module m_eos
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# 12 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
20# 1 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp" 1
21! Default for Fypp runs outside the toolchain (e.g. documentation). Searched last, so the
22! build/include/<target>/thermochem.fpp the toolchain writes from m_thermochem takes precedence.
23# 4 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp"
24# 5 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp"
25# 13 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 2
26
27# 16 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
28# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
29
30# 24 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
31
32# 53 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
33
34# 65 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
35
36# 75 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
37
38# 105 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
39
40# 117 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
41
42# 127 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
43
44# 174 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
45! New line at end of file is required for FYPP
46# 2 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
47# 1 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 1
48# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
49# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
50# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
51# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
52# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
53# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
54
55# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
56# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
57# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
58
59# 12 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
60# 1 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp" 1
61! Default for Fypp runs outside the toolchain (e.g. documentation). Searched last, so the
62! build/include/<target>/thermochem.fpp the toolchain writes from m_thermochem takes precedence.
63# 4 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp"
64# 5 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp"
65# 13 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 2
66
67# 16 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
68# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
69
70# 24 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
71
72# 53 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
73
74# 65 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
75
76# 75 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
77
78# 105 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
79
80# 117 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
81
82# 127 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
83
84# 174 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
85! New line at end of file is required for FYPP
86# 2 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 2
87
88# 4 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
89# 5 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
90# 6 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
91# 7 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
92# 8 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
93
94# 20 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
95
96# 43 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
97
98# 48 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
99
100# 53 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
101
102# 58 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
103
104# 63 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
105
106# 68 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
107
108# 81 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
109
110# 95 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
111
112# 100 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
113
114# 105 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
115
116# 110 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
117
118# 115 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
119
120# 120 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
121
122# 125 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
123
124# 130 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
125
126# 135 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
127
128# 165 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
129
130# 219 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
131
132# 233 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
133
134# 258 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
135
136# 269 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
137
138# 271 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
139# 282 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
140
141# 332 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
142
143# 342 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
144
145# 352 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
146
147# 361 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
148
149# 378 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
150
151# 388 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
152
153# 395 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
154
155# 401 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
156
157# 407 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
158
159# 413 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
160
161# 419 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
162
163# 425 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
164! New line at end of file is required for FYPP
165# 3 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
166# 1 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 1
167# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
168# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
169# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
170# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
171# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
172# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
173
174# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
175# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
176# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
177
178# 12 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
179# 1 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp" 1
180! Default for Fypp runs outside the toolchain (e.g. documentation). Searched last, so the
181! build/include/<target>/thermochem.fpp the toolchain writes from m_thermochem takes precedence.
182# 4 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp"
183# 5 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp"
184# 13 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 2
185
186# 16 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
187# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
188
189# 24 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
190
191# 53 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
192
193# 65 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
194
195# 75 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
196
197# 105 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
198
199# 117 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
200
201# 127 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
202
203# 174 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
204! New line at end of file is required for FYPP
205# 2 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 2
206
207# 7 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
208
209# 17 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
210
211# 22 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
212
213# 27 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
214
215# 32 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
216
217# 37 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
218
219# 42 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
220
221# 47 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
222
223# 52 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
224
225# 57 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
226
227# 62 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
228
229# 73 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
230
231# 78 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
232
233# 83 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
234
235# 88 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
236
237# 103 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
238
239# 131 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
240
241# 160 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
242
243# 175 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
244
245# 193 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
246
247# 215 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
248
249# 244 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
250
251# 259 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
252
253# 269 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
254
255# 278 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
256
257# 294 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
258
259# 304 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
260
261# 311 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
262! New line at end of file is required for FYPP
263# 4 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
264
265! GPU parallel region (scalar reductions, maxval/minval)
266# 23 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
267
268! GPU parallel loop over threads (most common GPU macro)
269# 43 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
270
271! Required closing for GPU_PARALLEL_LOOP
272# 55 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
273
274! Mark routine for device compilation
275# 112 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
276
277! Declare device-resident data
278# 130 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
279
280! Inner loop within a GPU parallel region
281# 145 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
282
283! Scoped GPU data region
284# 164 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
285
286! Host code with device pointers (for MPI with GPU buffers)
287# 193 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
288
289! Allocate device memory (unscoped)
290# 207 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
291
292! Free device memory
293# 219 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
294
295! Atomic operation on device
296# 231 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
297
298! End atomic capture block
299# 242 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
300
301! Copy data between host and device
302# 254 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
303
304! Synchronization barrier
305# 266 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
306
307! Import GPU library module (openacc or omp_lib)
308# 275 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
309
310! Emit code only for AMD compiler
311# 282 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
312
313! Emit code for non-Cray compilers
314# 289 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
315
316! Emit code only for Cray compiler
317# 296 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
318
319! Emit code for non-NVIDIA compilers
320# 303 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
321
322# 305 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
323# 306 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
324! New line at end of file is required for FYPP
325# 2 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 2
326
327# 14 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
328
329! Caution: This macro requires the use of a binding script to set CUDA_VISIBLE_DEVICES, such that we have one GPU device per MPI
330! rank. That's because for both cudaMemAdvise (preferred location) and cudaMemPrefetchAsync we use location = device_id = 0. For an
331! example see misc/nvidia_uvm/bind.sh.
332# 52 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
333
334! Allocate and create GPU device memory
335# 72 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
336
337! Free GPU device memory and deallocate
338# 80 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
339
340! Cray-specific GPU pointer setup for vector fields
341# 104 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
342
343! Cray-specific GPU pointer setup for scalar fields
344# 120 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
345
346! Cray-specific GPU pointer setup for acoustic source spatials
347# 145 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
348
349# 151 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
350
351# 158 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
352! New line at end of file is required for FYPP
353# 6 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp" 2
354# 1 "/home/runner/work/MFC/MFC/src/common/include/case.fpp" 1
355! This file exists so that Fypp can be run without generating case.fpp files for
356! each target. This is useful when generating documentation, for example. This
357! should also let MFC be built with CMake directly, without invoking mfc.sh.
358
359! For pre-process.
360# 8 "/home/runner/work/MFC/MFC/src/common/include/case.fpp"
361
362! For moving immersed boundaries in simulation
363# 12 "/home/runner/work/MFC/MFC/src/common/include/case.fpp"
364# 7 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp" 2
365# 1 "/home/runner/work/MFC/MFC/build/include/pre_process/generated_eos.fpp" 1
366! AUTO-GENERATED - do not edit directly. Regenerate: cmake reconfigure
367!
368# 4 "/home/runner/work/MFC/MFC/build/include/pre_process/generated_eos.fpp"
369
370# 6 "/home/runner/work/MFC/MFC/build/include/pre_process/generated_eos.fpp"
371# 9 "/home/runner/work/MFC/MFC/build/include/pre_process/generated_eos.fpp"
372
373# 12 "/home/runner/work/MFC/MFC/build/include/pre_process/generated_eos.fpp"
374# 15 "/home/runner/work/MFC/MFC/build/include/pre_process/generated_eos.fpp"
375
376# 17 "/home/runner/work/MFC/MFC/build/include/pre_process/generated_eos.fpp"
377# 31 "/home/runner/work/MFC/MFC/build/include/pre_process/generated_eos.fpp"
378
379# 34 "/home/runner/work/MFC/MFC/build/include/pre_process/generated_eos.fpp"
380# 58 "/home/runner/work/MFC/MFC/build/include/pre_process/generated_eos.fpp"
381# 8 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp" 2
382
383!> @brief Equations of state in Gamma/Pi form, rho e = Gamma(rho) p + Pi(rho).
384!!
385!! Stiffened and ideal gas keep constant coefficients, resolved once at start-up. The
386!! state-dependent families (Mie-Gruneisen, JWL, Vinet) evaluate theirs per cell from a
387!! reference curve.
388!!
389!! The mixture closure rules that combine the phases -- Wood's law, the six-equation mean, the
390!! bubbly branch -- are not themselves equations of state, but s_compute_mixture_coefficients,
391!! s_compute_speed_of_sound and their _dt/_avg variants live here rather than in
392!! m_variables_conversion, and must stay here. They are the hot-path callers of the phase
393!! chain (s_phase_coefficients -> s_eos_coefficients -> s_reference_curve), and on NVHPC that
394!! inlining only happens within a single file: the cross-file inliner refuses any device routine
395!! with a subroutine call in its call tree. Splitting them from the chain costs ~25% of grind
396!! time on NVHPC and nothing on the other backends, so it fails quietly. The solver kernels never
397!! inlined these four in the first place, which is why the module boundary is drawn above them and
398!! not below. See docs/documentation/gpuParallelization.md, "Module boundaries and NVHPC inlining".
399!!
400!! This module is a leaf: it directly uses only m_derived_types, m_constants, and
401!! m_global_parameters_common. Adding an EOS family means one case in s_reference_curve.
402module m_eos
403
409
410 implicit none
411
412 private
413
419
420contains
421
422 !> Resolve every fluid's EOS coefficients once, before any conversion runs.
423 impure subroutine s_initialize_eos_module()
424
425 integer :: i
426 logical :: state_dependent !< Whether this case's fluids need a density-dependent EOS
427
428#ifdef MFC_DEBUG
429# 54 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
430 block
431# 54 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
432 use iso_fortran_env, only: output_unit
433# 54 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
434
435# 54 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
436 print *, 'm_eos.fpp:54: ', '@:ALLOCATE(gammas (1:num_fluids))'
437# 54 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
438
439# 54 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
440 call flush (output_unit)
441# 54 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
442 end block
443# 54 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
444#endif
445# 54 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
446 allocate (gammas(1:num_fluids))
447# 54 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
448
449# 54 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
450
451# 54 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
452#if defined(MFC_OpenACC)
453# 54 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
454!$acc enter data create(gammas)
455# 54 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
456#elif defined(MFC_OpenMP)
457# 54 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
458!$omp target enter data map(always,alloc:gammas)
459# 54 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
460#endif
461# 55 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
462#ifdef MFC_DEBUG
463# 55 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
464 block
465# 55 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
466 use iso_fortran_env, only: output_unit
467# 55 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
468
469# 55 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
470 print *, 'm_eos.fpp:55: ', '@:ALLOCATE(eoss (1:num_fluids))'
471# 55 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
472
473# 55 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
474 call flush (output_unit)
475# 55 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
476 end block
477# 55 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
478#endif
479# 55 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
480 allocate (eoss(1:num_fluids))
481# 55 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
482
483# 55 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
484
485# 55 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
486#if defined(MFC_OpenACC)
487# 55 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
488!$acc enter data create(eoss)
489# 55 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
490#elif defined(MFC_OpenMP)
491# 55 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
492!$omp target enter data map(always,alloc:eoss)
493# 55 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
494#endif
495# 56 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
496#ifdef MFC_DEBUG
497# 56 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
498 block
499# 56 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
500 use iso_fortran_env, only: output_unit
501# 56 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
502
503# 56 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
504 print *, 'm_eos.fpp:56: ', '@:ALLOCATE(isentrope_n (1:num_fluids))'
505# 56 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
506
507# 56 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
508 call flush (output_unit)
509# 56 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
510 end block
511# 56 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
512#endif
513# 56 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
514 allocate (isentrope_n(1:num_fluids))
515# 56 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
516
517# 56 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
518
519# 56 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
520#if defined(MFC_OpenACC)
521# 56 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
522!$acc enter data create(isentrope_n)
523# 56 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
524#elif defined(MFC_OpenMP)
525# 56 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
526!$omp target enter data map(always,alloc:isentrope_n)
527# 56 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
528#endif
529# 57 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
530#ifdef MFC_DEBUG
531# 57 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
532 block
533# 57 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
534 use iso_fortran_env, only: output_unit
535# 57 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
536
537# 57 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
538 print *, 'm_eos.fpp:57: ', '@:ALLOCATE(pi_infs(1:num_fluids))'
539# 57 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
540
541# 57 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
542 call flush (output_unit)
543# 57 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
544 end block
545# 57 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
546#endif
547# 57 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
548 allocate (pi_infs(1:num_fluids))
549# 57 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
550
551# 57 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
552
553# 57 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
554#if defined(MFC_OpenACC)
555# 57 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
556!$acc enter data create(pi_infs)
557# 57 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
558#elif defined(MFC_OpenMP)
559# 57 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
560!$omp target enter data map(always,alloc:pi_infs)
561# 57 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
562#endif
563# 58 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
564#ifdef MFC_DEBUG
565# 58 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
566 block
567# 58 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
568 use iso_fortran_env, only: output_unit
569# 58 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
570
571# 58 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
572 print *, 'm_eos.fpp:58: ', '@:ALLOCATE(isentrope_B(1:num_fluids))'
573# 58 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
574
575# 58 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
576 call flush (output_unit)
577# 58 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
578 end block
579# 58 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
580#endif
581# 58 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
582 allocate (isentrope_b(1:num_fluids))
583# 58 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
584
585# 58 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
586
587# 58 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
588#if defined(MFC_OpenACC)
589# 58 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
590!$acc enter data create(isentrope_B)
591# 58 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
592#elif defined(MFC_OpenMP)
593# 58 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
594!$omp target enter data map(always,alloc:isentrope_B)
595# 58 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
596#endif
597# 59 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
598#ifdef MFC_DEBUG
599# 59 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
600 block
601# 59 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
602 use iso_fortran_env, only: output_unit
603# 59 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
604
605# 59 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
606 print *, 'm_eos.fpp:59: ', '@:ALLOCATE(cvs (1:num_fluids))'
607# 59 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
608
609# 59 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
610 call flush (output_unit)
611# 59 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
612 end block
613# 59 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
614#endif
615# 59 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
616 allocate (cvs(1:num_fluids))
617# 59 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
618
619# 59 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
620
621# 59 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
622#if defined(MFC_OpenACC)
623# 59 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
624!$acc enter data create(cvs)
625# 59 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
626#elif defined(MFC_OpenMP)
627# 59 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
628!$omp target enter data map(always,alloc:cvs)
629# 59 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
630#endif
631# 60 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
632#ifdef MFC_DEBUG
633# 60 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
634 block
635# 60 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
636 use iso_fortran_env, only: output_unit
637# 60 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
638
639# 60 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
640 print *, 'm_eos.fpp:60: ', '@:ALLOCATE(qvs (1:num_fluids))'
641# 60 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
642
643# 60 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
644 call flush (output_unit)
645# 60 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
646 end block
647# 60 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
648#endif
649# 60 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
650 allocate (qvs(1:num_fluids))
651# 60 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
652
653# 60 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
654
655# 60 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
656#if defined(MFC_OpenACC)
657# 60 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
658!$acc enter data create(qvs)
659# 60 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
660#elif defined(MFC_OpenMP)
661# 60 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
662!$omp target enter data map(always,alloc:qvs)
663# 60 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
664#endif
665# 61 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
666#ifdef MFC_DEBUG
667# 61 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
668 block
669# 61 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
670 use iso_fortran_env, only: output_unit
671# 61 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
672
673# 61 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
674 print *, 'm_eos.fpp:61: ', '@:ALLOCATE(qvps (1:num_fluids))'
675# 61 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
676
677# 61 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
678 call flush (output_unit)
679# 61 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
680 end block
681# 61 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
682#endif
683# 61 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
684 allocate (qvps(1:num_fluids))
685# 61 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
686
687# 61 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
688
689# 61 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
690#if defined(MFC_OpenACC)
691# 61 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
692!$acc enter data create(qvps)
693# 61 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
694#elif defined(MFC_OpenMP)
695# 61 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
696!$omp target enter data map(always,alloc:qvps)
697# 61 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
698#endif
699# 62 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
700
701 state_dependent = .false.
702 do i = 1, num_fluids
703 gammas(i) = fluid_pp(i)%gamma
704 isentrope_n(i) = f_isentrope_exponent(gammas(i))
705
706 ! Each EOS supplies its own coefficients. Resolved once here, not per cell: a branch in the mixture loop costs
707 ! registers in the Riemann kernels. An EOS whose coefficients depend on state must move to per-cell evaluation.
708 select case (fluid_pp(i)%eos)
709 case (eos_ideal_gas)
710 pi_infs(i) = 0._wp
711 case default
712 pi_infs(i) = fluid_pp(i)%pi_inf
713 end select
714 isentrope_b(i) = f_isentrope_pressure(pi_infs(i), gammas(i))
715 cvs(i) = fluid_pp(i)%cv
716 qvs(i) = fluid_pp(i)%qv
717 qvps(i) = fluid_pp(i)%qvp
718 eoss(i) = fluid_pp(i)%eos
719 ! Every fluid's single-source coefficients, mu_max among them: where a cubic Hugoniot fit turns over.
720 ! mu(u_p) peaks where c0 = s2 u_p^2 + 2 s3 u_p^3, and past it no shock state exists, so the Newton below
721 ! would wander. Solved once here, on the host.
722 eos_coeffs(i)%c0 = fluid_pp(i)%mg_c0
723# 84 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
724 eos_coeffs(i)%s = fluid_pp(i)%mg_s
725# 84 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
726 eos_coeffs(i)%s2 = fluid_pp(i)%mg_s2
727# 84 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
728 eos_coeffs(i)%s3 = fluid_pp(i)%mg_s3
729# 84 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
730 eos_coeffs(i)%mu_max = f_hugoniot_compression_limit(fluid_pp(i)%mg_c0, fluid_pp(i)%mg_s, fluid_pp(i)%mg_s2, &
731# 84 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
732 & fluid_pp(i)%mg_s3)
733# 84 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
734 eos_coeffs(i)%a = fluid_pp(i)%jwl_a
735# 84 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
736 eos_coeffs(i)%b = fluid_pp(i)%jwl_b
737# 84 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
738 eos_coeffs(i)%r1 = fluid_pp(i)%jwl_r1
739# 84 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
740 eos_coeffs(i)%r2 = fluid_pp(i)%jwl_r2
741# 84 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
742 eos_coeffs(i)%k0 = fluid_pp(i)%vinet_k0
743# 84 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
744 eos_coeffs(i)%k0p = fluid_pp(i)%vinet_k0p
745# 85 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
746 ! One reference state and Gruneisen closure for every family; the user-facing names keep their prefix.
747 select case (fluid_pp(i)%eos)
748# 86 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
749 case (eos_mie_gruneisen)
750# 86 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
751 eos_coeffs(i)%rho0 = fluid_pp(i)%mg_rho0
752# 86 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
753 eos_coeffs(i)%t0 = fluid_pp(i)%mg_t0
754# 86 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
755 eos_coeffs(i)%gruneisen0 = fluid_pp(i)%mg_gruneisen
756# 86 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
757 eos_coeffs(i)%gruneisen_a = fluid_pp(i)%mg_gruneisen_a
758# 86 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
759 case (eos_jwl)
760# 86 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
761 eos_coeffs(i)%rho0 = fluid_pp(i)%jwl_rho0
762# 86 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
763 eos_coeffs(i)%t0 = fluid_pp(i)%jwl_t0
764# 86 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
765 eos_coeffs(i)%gruneisen0 = fluid_pp(i)%jwl_omega
766# 86 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
767 eos_coeffs(i)%gruneisen_a = 0._wp
768# 86 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
769 case (eos_vinet)
770# 86 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
771 eos_coeffs(i)%rho0 = fluid_pp(i)%vinet_rho0
772# 86 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
773 eos_coeffs(i)%t0 = fluid_pp(i)%vinet_t0
774# 86 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
775 eos_coeffs(i)%gruneisen0 = fluid_pp(i)%vinet_gruneisen
776# 86 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
777 eos_coeffs(i)%gruneisen_a = fluid_pp(i)%vinet_gruneisen_a
778# 86 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
779 case default
780# 86 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
781 eos_coeffs(i)%rho0 = dflt_real
782# 86 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
783 eos_coeffs(i)%t0 = dflt_real
784# 86 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
785 eos_coeffs(i)%gruneisen0 = dflt_real
786# 86 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
787 eos_coeffs(i)%gruneisen_a = 0._wp
788# 86 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
789 end select
790# 87 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
791 if (f_is_state_dependent(i)) state_dependent = .true.
792 end do
793# 95 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
794 any_state_dependent_eos = state_dependent
795# 97 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
796
797# 97 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
798#if defined(MFC_OpenACC)
799# 97 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
800!$acc update device(gammas, isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps, eoss, eos_coeffs)
801# 97 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
802#elif defined(MFC_OpenMP)
803# 97 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
804!$omp target update to(gammas, isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps, eoss, eos_coeffs)
805# 97 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
806#endif
807# 98 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
808# 99 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
809
810# 99 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
811#if defined(MFC_OpenACC)
812# 99 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
813!$acc update device(any_state_dependent_eos)
814# 99 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
815#elif defined(MFC_OpenMP)
816# 99 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
817!$omp target update to(any_state_dependent_eos)
818# 99 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
819#endif
820# 100 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
821# 101 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
822
823 end subroutine s_initialize_eos_module
824
825 !> Deallocate the fluid property arrays allocated in s_initialize_eos_module.
826 impure subroutine s_finalize_eos_module()
827
828#ifdef MFC_DEBUG
829# 107 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
830 block
831# 107 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
832 use iso_fortran_env, only: output_unit
833# 107 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
834
835# 107 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
836 print *, 'm_eos.fpp:107: ', '@:DEALLOCATE(gammas, isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps, eoss)'
837# 107 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
838
839# 107 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
840 call flush (output_unit)
841# 107 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
842 end block
843# 107 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
844#endif
845# 107 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
846
847# 107 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
848#if defined(MFC_OpenACC)
849# 107 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
850!$acc exit data delete(gammas, isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps, eoss)
851# 107 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
852#elif defined(MFC_OpenMP)
853# 107 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
854!$omp target exit data map(release:gammas, isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps, eoss)
855# 107 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
856#endif
857# 107 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
858 deallocate (gammas, isentrope_n, pi_infs, isentrope_b, cvs, qvs, qvps, eoss)
859# 108 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
860
861 end subroutine s_finalize_eos_module
862
863 !> 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
864 !! adds one case here and nothing else.
865 subroutine s_reference_curve(rho, i, p_ref, e_ref, dp_drho, de_drho, G0, dG0)
866
867
868# 115 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
869#if MFC_OpenACC
870# 115 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
871!$acc routine seq
872# 115 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
873#elif MFC_OpenMP
874# 115 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
875
876# 115 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
877
878# 115 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
879!$omp declare target device_type(any)
880# 115 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
881#endif
882# 116 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
883
884 real(wp), intent(in) :: rho
885 integer, intent(in) :: i
886 real(wp), intent(out) :: p_ref, e_ref, dp_drho, de_drho, G0, dG0
887 real(wp) :: mu, d, V, ea, eb, up, us, dus, dup_dmu, x, ex, dp_dmu, de_dmu
888 integer :: iter
889
890 mu = rho/eos_coeffs(i)%rho0 - 1._wp
891 ! Past the fit's turnover there is no shock state to find; clamp rather than let the Newton below wander
892 ! off and return a silently wrong pressure. mu_max is huge for the linear fit, so this is a no-op there.
893 ! Bounded here so the step finishes and the host-side check in s_write_run_time_information can report it;
894 ! past the turnover there is no shock state and the Newton below would wander.
895 if (eoss(i) == eos_mie_gruneisen .and. mu > eos_coeffs(i)%mu_max) mu = eos_coeffs(i)%mu_max
896 select case (eoss(i))
897 case (eos_mie_gruneisen)
898 ! 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
899 ! energy e_H = p_H mu/(2 rho0 (1 + mu)); linear on release. Pole at mu = 1/(s - 1) for the linear fit;
900 ! the validator refuses initial states outside the EOS.
901 if (mu < 0._wp) then
902 p_ref = eos_coeffs(i)%rho0*eos_coeffs(i)%c0**2*mu
903 dp_dmu = eos_coeffs(i)%rho0*eos_coeffs(i)%c0**2
904 else if (eos_coeffs(i)%s2 == 0._wp .and. eos_coeffs(i)%s3 == 0._wp) then
905 d = 1._wp - (eos_coeffs(i)%s - 1._wp)*mu
906 p_ref = eos_coeffs(i)%rho0*eos_coeffs(i)%c0**2*mu*(1._wp + mu)/(d*d)
907 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 &
908 & + mu))/(d*d*d)
909 else
910 ! u_p solves u_s(u_p) mu = u_p (1 + mu): Newton from the linear fit, then implicit differentiation
911 up = eos_coeffs(i)%c0*mu/(1._wp - (eos_coeffs(i)%s - 1._wp)*mu)
912
913# 145 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
914#if defined(MFC_OpenACC)
915# 145 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
916!$acc loop seq
917# 145 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
918#elif defined(MFC_OpenMP)
919# 145 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
920
921# 145 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
922#endif
923# 146 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
924 do iter = 1, 8
925 us = eos_coeffs(i)%c0 + up*(eos_coeffs(i)%s + up*(eos_coeffs(i)%s2 + up*eos_coeffs(i)%s3))
926 dus = eos_coeffs(i)%s + up*(2._wp*eos_coeffs(i)%s2 + 3._wp*eos_coeffs(i)%s3*up)
927 up = up - (us*mu - up*(1._wp + mu))/(dus*mu - (1._wp + mu))
928 end do
929 us = eos_coeffs(i)%c0 + up*(eos_coeffs(i)%s + up*(eos_coeffs(i)%s2 + up*eos_coeffs(i)%s3))
930 dus = eos_coeffs(i)%s + up*(2._wp*eos_coeffs(i)%s2 + 3._wp*eos_coeffs(i)%s3*up)
931 dup_dmu = (up - us)/(dus*mu - (1._wp + mu))
932 p_ref = eos_coeffs(i)%rho0*us*up
933 dp_dmu = eos_coeffs(i)%rho0*(dus*up + us)*dup_dmu
934 end if
935 e_ref = p_ref*mu/(2._wp*eos_coeffs(i)%rho0*(1._wp + mu))
936 de_dmu = (dp_dmu*mu*(1._wp + mu) + p_ref)/(2._wp*eos_coeffs(i)%rho0*(1._wp + mu)**2)
937 dp_drho = dp_dmu/eos_coeffs(i)%rho0
938 de_drho = de_dmu/eos_coeffs(i)%rho0
939 case (eos_jwl)
940 ! 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).
941 v = eos_coeffs(i)%rho0/rho
942 ea = eos_coeffs(i)%a*exp(-eos_coeffs(i)%r1*v)
943 eb = eos_coeffs(i)%b*exp(-eos_coeffs(i)%r2*v)
944 p_ref = ea + eb
945 e_ref = (ea/eos_coeffs(i)%r1 + eb/eos_coeffs(i)%r2)/eos_coeffs(i)%rho0
946 dp_drho = (eos_coeffs(i)%rho0/rho**2)*(eos_coeffs(i)%r1*ea + eos_coeffs(i)%r2*eb)
947 de_drho = p_ref/rho**2
948 case (eos_vinet)
949 ! 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,
950 ! an isentrope like JWL (its energy integrates in closed form).
951 d = 1.5_wp*(eos_coeffs(i)%k0p - 1._wp)
952 x = (eos_coeffs(i)%rho0/rho)**(1._wp/3._wp)
953 ex = exp(d*(1._wp - x))
954 p_ref = 3._wp*eos_coeffs(i)%k0*(1._wp - x)/x**2*ex
955 e_ref = 9._wp*eos_coeffs(i)%k0/(eos_coeffs(i)%rho0*d**2)*(1._wp - (1._wp - d*(1._wp - x))*ex)
956 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))
957 de_drho = p_ref/rho**2
958 end select
959 g0 = eos_coeffs(i)%gruneisen0 + eos_coeffs(i)%gruneisen_a*mu
960 dg0 = eos_coeffs(i)%gruneisen_a/eos_coeffs(i)%rho0
961
962 end subroutine s_reference_curve
963
964 !> Whether the EOS of fluid i is a family whose coefficients vary with density.
965 function f_is_state_dependent(i) result(yes)
966
967
968# 189 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
969#ifdef _CRAYFTN
970# 189 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
971#if MFC_OpenACC
972# 189 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
973!$acc routine seq
974# 189 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
975#elif MFC_OpenMP
976# 189 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
977
978# 189 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
979
980# 189 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
981!$omp declare target device_type(any)
982# 189 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
983#else
984# 189 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
985!DIR$ INLINEALWAYS f_is_state_dependent
986# 189 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
987#endif
988# 189 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
989#elif MFC_OpenACC
990# 189 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
991!$acc routine seq
992# 189 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
993#elif MFC_OpenMP
994# 189 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
995
996# 189 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
997
998# 189 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
999!$omp declare target device_type(any)
1000# 189 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1001#endif
1002# 190 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1003
1004 integer, intent(in) :: i
1005 logical :: yes
1006
1007 yes = eoss(i) == eos_mie_gruneisen .or. eoss(i) == eos_jwl .or. eoss(i) == eos_vinet
1008# 195 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1009
1010 end function f_is_state_dependent
1011
1012 !> 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
1013 !! for the Mie-Gruneisen Hugoniot) and its Gruneisen coefficient is constant. Those two together make the isentrope through any
1014 !! state closed-form, so it never has to be integrated.
1015 function f_has_isentropic_reference(i) result(yes)
1016
1017
1018# 203 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1019#ifdef _CRAYFTN
1020# 203 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1021#if MFC_OpenACC
1022# 203 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1023!$acc routine seq
1024# 203 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1025#elif MFC_OpenMP
1026# 203 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1027
1028# 203 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1029
1030# 203 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1031!$omp declare target device_type(any)
1032# 203 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1033#else
1034# 203 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1035!DIR$ INLINEALWAYS f_has_isentropic_reference
1036# 203 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1037#endif
1038# 203 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1039#elif MFC_OpenACC
1040# 203 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1041!$acc routine seq
1042# 203 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1043#elif MFC_OpenMP
1044# 203 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1045
1046# 203 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1047
1048# 203 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1049!$omp declare target device_type(any)
1050# 203 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1051#endif
1052# 204 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1053
1054 integer, intent(in) :: i
1055 logical :: yes
1056
1057 yes = (eoss(i) == eos_jwl .or. eoss(i) == eos_vinet) .and. eos_coeffs(i)%gruneisen_a == 0._wp
1058# 209 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1059
1060 end function f_has_isentropic_reference
1061
1062 !> 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
1063 !! u_p^3, and falls after; only the rising branch is a physical shock. Returns a huge value for the linear fit, which never
1064 !! turns over. Host-side: called once per fluid at initialization.
1065 impure function f_hugoniot_compression_limit(c0, s, s2, s3) result(mu_max)
1066
1067 real(wp), intent(in) :: c0, s, s2, s3
1068 real(wp) :: mu_max, up, f, df, us
1069 integer :: iter
1070
1071 if (s2 == 0._wp .and. s3 == 0._wp) then
1072 mu_max = huge(1._wp)
1073 return
1074 end if
1075
1076 ! Newton on c0 - s2 u^2 - 2 s3 u^3 = 0, from a guess that brackets the physical range
1077 up = c0
1078 do iter = 1, 100
1079 f = c0 - s2*up**2 - 2._wp*s3*up**3
1080 df = -2._wp*s2*up - 6._wp*s3*up**2
1081 if (abs(df) < verysmall) exit
1082 up = max(up - f/df, verysmall)
1083 end do
1084 us = c0 + up*(s + up*(s2 + up*s3))
1085 mu_max = up/max(us - up, verysmall)
1086
1087 end function f_hugoniot_compression_limit
1088
1089 !> Gamma, Pi, dPi/drho and dGamma/drho of fluid i at density rho, the coefficients of rho e = Gamma p + Pi(rho). Stiffened and
1090 !! ideal gas keep the constants resolved at init, bit for bit.
1091 subroutine s_eos_coefficients(rho, i, gamma, pi_inf, dpi, dgamma)
1092
1093
1094# 243 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1095#if MFC_OpenACC
1096# 243 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1097!$acc routine seq
1098# 243 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1099#elif MFC_OpenMP
1100# 243 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1101
1102# 243 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1103
1104# 243 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1105!$omp declare target device_type(any)
1106# 243 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1107#endif
1108# 244 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1109
1110 real(wp), intent(in) :: rho
1111 integer, intent(in) :: i
1112 real(wp), intent(out) :: gamma, pi_inf, dpi, dgamma
1113 real(wp) :: p_ref, e_ref, dp_drho, de_drho, G0, dG0
1114
1115 if (.not. f_is_state_dependent(i)) then
1116 gamma = gammas(i)
1117 pi_inf = pi_infs(i)
1118 dpi = 0._wp
1119 dgamma = 0._wp
1120 return
1121 end if
1122 call s_reference_curve(rho, i, p_ref, e_ref, dp_drho, de_drho, g0, dg0)
1123 gamma = 1._wp/g0
1124 pi_inf = rho*e_ref - p_ref/g0
1125 dpi = e_ref + rho*de_drho - dp_drho/g0 + p_ref*dg0/g0**2
1126 dgamma = -dg0/g0**2
1127
1128 end subroutine s_eos_coefficients
1129
1130 !> Exponent of the stiffened-gas isentrope p + B = const rho**n. Precomputed per fluid as isentrope_n.
1131 function f_isentrope_exponent(gamma) result(n)
1132
1133
1134# 268 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1135#ifdef _CRAYFTN
1136# 268 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1137#if MFC_OpenACC
1138# 268 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1139!$acc routine seq
1140# 268 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1141#elif MFC_OpenMP
1142# 268 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1143
1144# 268 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1145
1146# 268 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1147!$omp declare target device_type(any)
1148# 268 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1149#else
1150# 268 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1151!DIR$ INLINEALWAYS f_isentrope_exponent
1152# 268 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1153#endif
1154# 268 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1155#elif MFC_OpenACC
1156# 268 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1157!$acc routine seq
1158# 268 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1159#elif MFC_OpenMP
1160# 268 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1161
1162# 268 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1163
1164# 268 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1165!$omp declare target device_type(any)
1166# 268 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1167#endif
1168# 269 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1169
1170 real(wp), intent(in) :: gamma
1171 real(wp) :: n
1172
1173 n = 1._wp/gamma + 1._wp
1174
1175 end function f_isentrope_exponent
1176
1177 !> Reference pressure of that isentrope. Precomputed per fluid as isentrope_B.
1178 function f_isentrope_pressure(pi_inf, gamma) result(B)
1179
1180
1181# 280 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1182#ifdef _CRAYFTN
1183# 280 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1184#if MFC_OpenACC
1185# 280 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1186!$acc routine seq
1187# 280 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1188#elif MFC_OpenMP
1189# 280 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1190
1191# 280 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1192
1193# 280 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1194!$omp declare target device_type(any)
1195# 280 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1196#else
1197# 280 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1198!DIR$ INLINEALWAYS f_isentrope_pressure
1199# 280 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1200#endif
1201# 280 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1202#elif MFC_OpenACC
1203# 280 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1204!$acc routine seq
1205# 280 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1206#elif MFC_OpenMP
1207# 280 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1208
1209# 280 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1210
1211# 280 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1212!$omp declare target device_type(any)
1213# 280 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1214#endif
1215# 281 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1216
1217 real(wp), intent(in) :: pi_inf, gamma
1218 real(wp) :: b
1219
1220 b = pi_inf/(1._wp + gamma)
1221
1222 end function f_isentrope_pressure
1223
1224 !> Stiffened-gas thermal law p + B = (n - 1)*cv*rho*T. Pass rho to get T, or T to get rho.
1225 function f_sg_thermal(pres, rho_or_T, n, B, cv) result(T_or_rho)
1226
1227
1228# 292 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1229#ifdef _CRAYFTN
1230# 292 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1231#if MFC_OpenACC
1232# 292 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1233!$acc routine seq
1234# 292 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1235#elif MFC_OpenMP
1236# 292 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1237
1238# 292 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1239
1240# 292 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1241!$omp declare target device_type(any)
1242# 292 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1243#else
1244# 292 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1245!DIR$ INLINEALWAYS f_sg_thermal
1246# 292 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1247#endif
1248# 292 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1249#elif MFC_OpenACC
1250# 292 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1251!$acc routine seq
1252# 292 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1253#elif MFC_OpenMP
1254# 292 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1255
1256# 292 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1257
1258# 292 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1259!$omp declare target device_type(any)
1260# 292 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1261#endif
1262# 293 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1263
1264 real(wp), intent(in) :: pres, rho_or_t, n, b, cv
1265 real(wp) :: t_or_rho
1266
1267 t_or_rho = (pres + b)/((n - 1._wp)*cv*rho_or_t)
1268
1269 end function f_sg_thermal
1270
1271 !> Thermal-equilibrium mixture temperature for stiffened gas, from primitives. Algebraically identical to the conservative form
1272 !! 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
1273 !! rho*e = gamma_mix*p + pi_inf_mix + sum(alpha_rho_i*qv_i) in MFC's stored variables.
1274 function f_mixture_temperature(alpha_rho_K, pres, gamma_K, pi_inf_K) result(T)
1275
1276
1277# 306 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1278#ifdef _CRAYFTN
1279# 306 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1280#if MFC_OpenACC
1281# 306 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1282!$acc routine seq
1283# 306 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1284#elif MFC_OpenMP
1285# 306 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1286
1287# 306 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1288
1289# 306 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1290!$omp declare target device_type(any)
1291# 306 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1292#else
1293# 306 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1294!DIR$ INLINEALWAYS f_mixture_temperature
1295# 306 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1296#endif
1297# 306 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1298#elif MFC_OpenACC
1299# 306 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1300!$acc routine seq
1301# 306 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1302#elif MFC_OpenMP
1303# 306 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1304
1305# 306 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1306
1307# 306 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1308!$omp declare target device_type(any)
1309# 306 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1310#endif
1311# 307 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1312
1313# 311 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1314 real(wp), dimension(num_fluids), intent(in) :: alpha_rho_k
1315# 313 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1316 real(wp), intent(in) :: pres, gamma_k, pi_inf_k
1317 real(wp) :: t
1318 real(wp) :: mcp !< sum of alpha_rho_i*cp_i; cp_i = n_i*cv_i
1319 integer :: i
1320
1321 mcp = 0._wp
1322
1323# 319 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1324#if defined(MFC_OpenACC)
1325# 319 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1326!$acc loop seq
1327# 319 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1328#elif defined(MFC_OpenMP)
1329# 319 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1330
1331# 319 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1332#endif
1333# 320 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1334 do i = 1, num_fluids
1335 mcp = mcp + alpha_rho_k(i)*cvs(i)*isentrope_n(i)
1336 end do
1337
1338 t = ((gamma_k + 1._wp)*pres + pi_inf_k)/max(mcp, sgm_eps)
1339
1340 end function f_mixture_temperature
1341
1342 !> Coefficients of phase i at its own density alpha_rho/alpha: the per-cell dispatch when some fluid's EOS is state dependent,
1343 !! the constants resolved at init otherwise (bit for bit).
1344 subroutine s_phase_coefficients(alpha_rho, alpha, i, rho, gamma, pi_inf, dpi, dgamma)
1345
1346
1347# 332 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1348#ifdef _CRAYFTN
1349# 332 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1350#if MFC_OpenACC
1351# 332 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1352!$acc routine seq
1353# 332 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1354#elif MFC_OpenMP
1355# 332 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1356
1357# 332 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1358
1359# 332 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1360!$omp declare target device_type(any)
1361# 332 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1362#else
1363# 332 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1364!DIR$ INLINEALWAYS s_phase_coefficients
1365# 332 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1366#endif
1367# 332 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1368#elif MFC_OpenACC
1369# 332 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1370!$acc routine seq
1371# 332 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1372#elif MFC_OpenMP
1373# 332 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1374
1375# 332 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1376
1377# 332 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1378!$omp declare target device_type(any)
1379# 332 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1380#endif
1381# 333 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1382
1383 real(wp), intent(in) :: alpha_rho, alpha
1384 integer, intent(in) :: i
1385 real(wp), intent(out) :: rho, gamma, pi_inf, dpi, dgamma
1386
1387 rho = max(alpha_rho, sgm_eps)/max(alpha, sgm_eps)
1388 if (any_state_dependent_eos) then
1389 call s_eos_coefficients(rho, i, gamma, pi_inf, dpi, dgamma)
1390 else
1391 gamma = gammas(i)
1392 pi_inf = pi_infs(i)
1393 dpi = 0._wp
1394 dgamma = 0._wp
1395 end if
1396
1397 end subroutine s_phase_coefficients
1398
1399 !> c^2 = [((Gamma + 1) p + Pi)/rho - dPi/drho - p dGamma/drho]/Gamma, the frozen speed of one phase.
1400 function f_c2_from_coefficients(rho, pres, gamma, pi_inf, dpi, dgamma) result(c2)
1401
1402
1403# 353 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1404#ifdef _CRAYFTN
1405# 353 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1406#if MFC_OpenACC
1407# 353 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1408!$acc routine seq
1409# 353 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1410#elif MFC_OpenMP
1411# 353 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1412
1413# 353 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1414
1415# 353 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1416!$omp declare target device_type(any)
1417# 353 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1418#else
1419# 353 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1420!DIR$ INLINEALWAYS f_c2_from_coefficients
1421# 353 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1422#endif
1423# 353 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1424#elif MFC_OpenACC
1425# 353 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1426!$acc routine seq
1427# 353 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1428#elif MFC_OpenMP
1429# 353 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1430
1431# 353 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1432
1433# 353 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1434!$omp declare target device_type(any)
1435# 353 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1436#endif
1437# 354 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1438
1439 real(wp), intent(in) :: rho, pres, gamma, pi_inf, dpi, dgamma
1440 real(wp) :: c2
1441
1442 c2 = (((gamma + 1._wp)*pres + pi_inf)/rho - dpi - pres*dgamma)/gamma
1443
1444 end function f_c2_from_coefficients
1445
1446 !> Frozen sound speed squared of one phase at (rho, p) from its own coefficients. These helpers are subroutines, not functions:
1447 !! a device function that calls a device subroutine is a pattern no other backend-tested code in MFC uses.
1448 subroutine s_phase_c2(rho, pres, i, c2)
1449
1450
1451# 366 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1452#if MFC_OpenACC
1453# 366 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1454!$acc routine seq
1455# 366 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1456#elif MFC_OpenMP
1457# 366 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1458
1459# 366 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1460
1461# 366 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1462!$omp declare target device_type(any)
1463# 366 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1464#endif
1465# 367 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1466
1467 real(wp), intent(in) :: rho, pres
1468 integer, intent(in) :: i
1469 real(wp), intent(out) :: c2
1470 real(wp) :: gamma, pi_inf, dpi, dgamma
1471
1472 call s_eos_coefficients(rho, i, gamma, pi_inf, dpi, dgamma)
1473 c2 = f_c2_from_coefficients(rho, pres, gamma, pi_inf, dpi, dgamma)
1474
1475 end subroutine s_phase_c2
1476
1477 !> 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 =
1478 !! (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).
1479 subroutine s_ode_slope(kind, i, x, y, dydx)
1480
1481
1482# 382 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1483#if MFC_OpenACC
1484# 382 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1485!$acc routine seq
1486# 382 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1487#elif MFC_OpenMP
1488# 382 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1489
1490# 382 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1491
1492# 382 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1493!$omp declare target device_type(any)
1494# 382 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1495#endif
1496# 383 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1497
1498 integer, intent(in) :: kind, i
1499 real(wp), intent(in) :: x, y
1500 real(wp), intent(out) :: dydx
1501 real(wp) :: p_ref, e_ref, dp_drho, de_drho, G0, dG0
1502
1503 if (kind == ode_isentrope) then
1504 call s_phase_c2(x, y, i, dydx)
1505 else
1506 call s_reference_curve(1._wp/x, i, p_ref, e_ref, dp_drho, de_drho, g0, dg0)
1507 dydx = (p_ref - de_drho/x**2)/cvs(i) - g0*y/x
1508 end if
1509
1510 end subroutine s_ode_slope
1511
1512 !> Fixed-step classical RK4 for the ODE `kind` from (x0, y0) to x1.
1513 subroutine s_rk4(kind, i, x0, y0, x1, y)
1514
1515
1516# 401 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1517#if MFC_OpenACC
1518# 401 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1519!$acc routine seq
1520# 401 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1521#elif MFC_OpenMP
1522# 401 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1523
1524# 401 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1525
1526# 401 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1527!$omp declare target device_type(any)
1528# 401 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1529#endif
1530# 402 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1531
1532 integer, intent(in) :: kind, i
1533 real(wp), intent(in) :: x0, y0, x1
1534 real(wp), intent(out) :: y
1535 real(wp) :: x, h, k1, k2, k3, k4
1536 integer :: step
1537
1538 x = x0
1539 y = y0
1540 h = (x1 - x0)/eos_rk4_steps
1541
1542# 412 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1543#if defined(MFC_OpenACC)
1544# 412 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1545!$acc loop seq
1546# 412 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1547#elif defined(MFC_OpenMP)
1548# 412 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1549
1550# 412 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1551#endif
1552# 413 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1553 do step = 1, eos_rk4_steps
1554 call s_ode_slope(kind, i, x, y, k1)
1555 call s_ode_slope(kind, i, x + 0.5_wp*h, y + 0.5_wp*h*k1, k2)
1556 call s_ode_slope(kind, i, x + 0.5_wp*h, y + 0.5_wp*h*k2, k3)
1557 call s_ode_slope(kind, i, x + h, y + h*k3, k4)
1558 y = y + h*(k1 + 2._wp*(k2 + k3) + k4)/6._wp
1559 x = x + h
1560 end do
1561
1562 end subroutine s_rk4
1563
1564 !> Pressure of phase i after the isentropic density change rho -> xi rho: closed form for the constant-coefficient families,
1565 !! integrated for a state-dependent EOS (the star states it serves are close to rho).
1566 subroutine s_phase_pressure_on_isentrope(pres, rho, xi, i, p_isen)
1567
1568
1569# 428 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1570#ifdef _CRAYFTN
1571# 428 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1572#if MFC_OpenACC
1573# 428 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1574!$acc routine seq
1575# 428 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1576#elif MFC_OpenMP
1577# 428 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1578
1579# 428 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1580
1581# 428 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1582!$omp declare target device_type(any)
1583# 428 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1584#else
1585# 428 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1586!DIR$ INLINEALWAYS s_phase_pressure_on_isentrope
1587# 428 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1588#endif
1589# 428 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1590#elif MFC_OpenACC
1591# 428 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1592!$acc routine seq
1593# 428 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1594#elif MFC_OpenMP
1595# 428 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1596
1597# 428 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1598
1599# 428 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1600!$omp declare target device_type(any)
1601# 428 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1602#endif
1603# 429 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1604
1605 real(wp), intent(in) :: pres, rho, xi
1606 integer, intent(in) :: i
1607 real(wp), intent(out) :: p_isen
1608 real(wp) :: p_ref_from, p_ref_to, e_ref, dp_drho, de_drho, g0, dg0
1609
1610 if (.not. f_is_state_dependent(i)) then
1611 p_isen = (pres + isentrope_b(i))*xi**isentrope_n(i) - isentrope_b(i)
1612 else if (f_has_isentropic_reference(i)) then
1613 ! Exact: the offset from an isentropic reference obeys dDelta/Delta = Gamma drho/rho, so
1614 ! p - p_ref scales as (rho'/rho)**(1 + Gamma). Integrating it instead costs a decimal per
1615 ! doubling of the expansion and turns the pressure negative past roughly twentyfold.
1616 call s_reference_curve(rho, i, p_ref_from, e_ref, dp_drho, de_drho, g0, dg0)
1617 call s_reference_curve(xi*rho, i, p_ref_to, e_ref, dp_drho, de_drho, g0, dg0)
1618 p_isen = p_ref_to + (pres - p_ref_from)*xi**(1._wp + eos_coeffs(i)%gruneisen0)
1619 else
1620 call s_rk4(ode_isentrope, i, rho, pres, xi*rho, p_isen)
1621 end if
1622
1623 end subroutine s_phase_pressure_on_isentrope
1624
1625 !> Temperature of phase i at (rho, p): the stiffened-gas relation, or T_ref(rho) + (e - e_ref)/c_v.
1626 subroutine s_phase_temperature(rho, pres, i, T)
1627
1628
1629# 453 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1630#ifdef _CRAYFTN
1631# 453 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1632#if MFC_OpenACC
1633# 453 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1634!$acc routine seq
1635# 453 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1636#elif MFC_OpenMP
1637# 453 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1638
1639# 453 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1640
1641# 453 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1642!$omp declare target device_type(any)
1643# 453 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1644#else
1645# 453 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1646!DIR$ INLINEALWAYS s_phase_temperature
1647# 453 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1648#endif
1649# 453 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1650#elif MFC_OpenACC
1651# 453 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1652!$acc routine seq
1653# 453 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1654#elif MFC_OpenMP
1655# 453 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1656
1657# 453 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1658
1659# 453 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1660!$omp declare target device_type(any)
1661# 453 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1662#endif
1663# 454 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1664
1665 real(wp), intent(in) :: rho, pres
1666 integer, intent(in) :: i
1667 real(wp), intent(out) :: t
1668 real(wp) :: p_ref, e_ref, dp_drho, de_drho, g0, dg0, t0, t_ref
1669
1670 if (f_is_state_dependent(i)) then
1671 call s_reference_curve(rho, i, p_ref, e_ref, dp_drho, de_drho, g0, dg0)
1672 t0 = eos_coeffs(i)%t0
1673 call s_rk4(ode_reference_temperature, i, 1._wp/eos_coeffs(i)%rho0, t0, 1._wp/rho, t_ref)
1674 t = t_ref + (pres - p_ref)/(rho*g0*cvs(i))
1675 else
1676 t = (pres + isentrope_b(i))/((isentrope_n(i) - 1._wp)*cvs(i)*rho)
1677 end if
1678
1679 end subroutine s_phase_temperature
1680
1681 !> Density of phase i at (p, T), the inverse of s_phase_temperature (constant-coefficient EOS only).
1682 subroutine s_phase_density(pres, T, i, rho)
1683
1684
1685# 474 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1686#ifdef _CRAYFTN
1687# 474 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1688#if MFC_OpenACC
1689# 474 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1690!$acc routine seq
1691# 474 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1692#elif MFC_OpenMP
1693# 474 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1694
1695# 474 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1696
1697# 474 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1698!$omp declare target device_type(any)
1699# 474 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1700#else
1701# 474 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1702!DIR$ INLINEALWAYS s_phase_density
1703# 474 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1704#endif
1705# 474 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1706#elif MFC_OpenACC
1707# 474 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1708!$acc routine seq
1709# 474 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1710#elif MFC_OpenMP
1711# 474 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1712
1713# 474 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1714
1715# 474 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1716!$omp declare target device_type(any)
1717# 474 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1718#endif
1719# 475 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1720
1721 real(wp), intent(in) :: pres, t
1722 integer, intent(in) :: i
1723 real(wp), intent(out) :: rho
1724
1725 rho = (pres + isentrope_b(i))/((isentrope_n(i) - 1._wp)*cvs(i)*t)
1726
1727 end subroutine s_phase_density
1728
1729 !> Density of phase i on the isentrope through (rho_from, p_from) at p_to, and c^2 there: Newton on the pressure integrator,
1730 !! whose slope is c^2. The relaxation's own Newton wraps this, so a few steps suffice.
1731 subroutine s_phase_density_on_isentrope(i, rho_from, p_from, p_to, rho_to, c2_to)
1732
1733
1734# 488 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1735#if MFC_OpenACC
1736# 488 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1737!$acc routine seq
1738# 488 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1739#elif MFC_OpenMP
1740# 488 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1741
1742# 488 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1743
1744# 488 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1745!$omp declare target device_type(any)
1746# 488 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1747#endif
1748# 489 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1749
1750 integer, intent(in) :: i
1751 real(wp), intent(in) :: rho_from, p_from, p_to
1752 real(wp), intent(out) :: rho_to, c2_to
1753 real(wp) :: p_at, c2_at
1754 integer :: iter
1755
1756 rho_to = rho_from
1757
1758# 497 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1759#if defined(MFC_OpenACC)
1760# 497 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1761!$acc loop seq
1762# 497 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1763#elif defined(MFC_OpenMP)
1764# 497 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1765
1766# 497 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1767#endif
1768# 498 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1769 do iter = 1, 4
1770 call s_phase_pressure_on_isentrope(p_from, rho_from, rho_to/rho_from, i, p_at)
1771 call s_phase_c2(rho_to, p_at, i, c2_at)
1772 rho_to = rho_to - (p_at - p_to)/c2_at
1773 end do
1774 call s_phase_c2(rho_to, p_to, i, c2_to)
1775
1776 end subroutine s_phase_density_on_isentrope
1777
1778 !> Internal energy per unit volume of phase i at pressure pres: alpha (Gamma p + Pi) + alpha_rho qv, with the coefficients at
1779 !! the phase's own density.
1780 subroutine s_phase_internal_energy(pres, alpha, alpha_rho, i, e_phase)
1781
1782
1783# 511 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1784#ifdef _CRAYFTN
1785# 511 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1786#if MFC_OpenACC
1787# 511 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1788!$acc routine seq
1789# 511 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1790#elif MFC_OpenMP
1791# 511 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1792
1793# 511 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1794
1795# 511 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1796!$omp declare target device_type(any)
1797# 511 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1798#else
1799# 511 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1800!DIR$ INLINEALWAYS s_phase_internal_energy
1801# 511 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1802#endif
1803# 511 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1804#elif MFC_OpenACC
1805# 511 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1806!$acc routine seq
1807# 511 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1808#elif MFC_OpenMP
1809# 511 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1810
1811# 511 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1812
1813# 511 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1814!$omp declare target device_type(any)
1815# 511 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1816#endif
1817# 512 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1818
1819 real(wp), intent(in) :: pres, alpha, alpha_rho
1820 integer, intent(in) :: i
1821 real(wp), intent(out) :: e_phase
1822 real(wp) :: rho, gamma, pi_inf, dpi, dgamma
1823
1824 call s_phase_coefficients(alpha_rho, alpha, i, rho, gamma, pi_inf, dpi, dgamma)
1825 e_phase = alpha*(gamma*pres + pi_inf) + alpha_rho*qvs(i)
1826
1827 end subroutine s_phase_internal_energy
1828
1829 !> Bulk modulus rho c^2 of phase i at pressure pres: f_bulk_modulus for a constant-coefficient fluid, bit for bit, minus the
1830 !! reference-curve terms rho (dPi/drho + p dGamma/drho)/Gamma otherwise.
1831 subroutine s_phase_bulk_modulus(pres, alpha, alpha_rho, i, blkmod)
1832
1833
1834# 527 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1835#ifdef _CRAYFTN
1836# 527 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1837#if MFC_OpenACC
1838# 527 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1839!$acc routine seq
1840# 527 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1841#elif MFC_OpenMP
1842# 527 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1843
1844# 527 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1845
1846# 527 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1847!$omp declare target device_type(any)
1848# 527 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1849#else
1850# 527 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1851!DIR$ INLINEALWAYS s_phase_bulk_modulus
1852# 527 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1853#endif
1854# 527 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1855#elif MFC_OpenACC
1856# 527 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1857!$acc routine seq
1858# 527 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1859#elif MFC_OpenMP
1860# 527 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1861
1862# 527 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1863
1864# 527 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1865!$omp declare target device_type(any)
1866# 527 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1867#endif
1868# 528 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1869
1870 real(wp), intent(in) :: alpha_rho, alpha, pres
1871 integer, intent(in) :: i
1872 real(wp), intent(out) :: blkmod
1873 real(wp) :: rho, gamma, pi_inf, dpi, dgamma
1874
1875 call s_phase_coefficients(alpha_rho, alpha, i, rho, gamma, pi_inf, dpi, dgamma)
1876 blkmod = f_bulk_modulus(pres, gamma, pi_inf) - rho*(dpi + pres*dgamma)/gamma
1877
1878 end subroutine s_phase_bulk_modulus
1879
1880 !> Pressure of a stiffened gas from its internal energy density - the inverse of s_compute_energy. Callers subtract the kinetic,
1881 !! magnetic and elastic energy first; none of those are equation-of-state terms.
1882 function f_pressure(e_int, gamma, pi_inf, qv) result(pres)
1883
1884
1885# 543 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1886#ifdef _CRAYFTN
1887# 543 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1888#if MFC_OpenACC
1889# 543 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1890!$acc routine seq
1891# 543 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1892#elif MFC_OpenMP
1893# 543 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1894
1895# 543 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1896
1897# 543 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1898!$omp declare target device_type(any)
1899# 543 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1900#else
1901# 543 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1902!DIR$ INLINEALWAYS f_pressure
1903# 543 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1904#endif
1905# 543 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1906#elif MFC_OpenACC
1907# 543 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1908!$acc routine seq
1909# 543 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1910#elif MFC_OpenMP
1911# 543 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1912
1913# 543 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1914
1915# 543 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1916!$omp declare target device_type(any)
1917# 543 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1918#endif
1919# 544 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1920
1921 real(wp), intent(in) :: e_int, gamma, pi_inf, qv
1922 real(wp) :: pres
1923
1924 pres = (e_int - pi_inf - qv)/gamma
1925
1926 end function f_pressure
1927
1928 !> Isentropic bulk modulus. Takes coefficients rather than a fluid index, so a mixture - whose effective gamma and pi_inf come
1929 !! from s_compute_mixture_coefficients - is the same call as a single fluid. Elastic callers add their own shear term.
1930 function f_bulk_modulus(pres, gamma, pi_inf) result(blkmod)
1931
1932
1933# 556 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1934#ifdef _CRAYFTN
1935# 556 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1936#if MFC_OpenACC
1937# 556 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1938!$acc routine seq
1939# 556 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1940#elif MFC_OpenMP
1941# 556 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1942
1943# 556 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1944
1945# 556 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1946!$omp declare target device_type(any)
1947# 556 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1948#else
1949# 556 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1950!DIR$ INLINEALWAYS f_bulk_modulus
1951# 556 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1952#endif
1953# 556 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1954#elif MFC_OpenACC
1955# 556 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1956!$acc routine seq
1957# 556 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1958#elif MFC_OpenMP
1959# 556 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1960
1961# 556 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1962
1963# 556 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1964!$omp declare target device_type(any)
1965# 556 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1966#endif
1967# 557 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1968
1969 real(wp), intent(in) :: pres, gamma, pi_inf
1970 real(wp) :: blkmod
1971
1972 blkmod = ((gamma + 1._wp)*pres + pi_inf)/gamma
1973
1974 end function f_bulk_modulus
1975
1976 !> Relativistic specific enthalpy, h = 1 + (Gamma + 1)p/rho. Ideal gas only: the stiffness does not appear, so a fluid with a
1977 !! nonzero pi_inf is not represented here (the validator refuses that combination).
1978 function f_relativistic_enthalpy(pres, rho, gamma) result(H)
1979
1980
1981# 569 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1982#ifdef _CRAYFTN
1983# 569 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1984#if MFC_OpenACC
1985# 569 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1986!$acc routine seq
1987# 569 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1988#elif MFC_OpenMP
1989# 569 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1990
1991# 569 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1992
1993# 569 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1994!$omp declare target device_type(any)
1995# 569 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1996#else
1997# 569 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1998!DIR$ INLINEALWAYS f_relativistic_enthalpy
1999# 569 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2000#endif
2001# 569 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2002#elif MFC_OpenACC
2003# 569 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2004!$acc routine seq
2005# 569 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2006#elif MFC_OpenMP
2007# 569 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2008
2009# 569 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2010
2011# 569 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2012!$omp declare target device_type(any)
2013# 569 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2014#endif
2015# 570 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2016
2017 real(wp), intent(in) :: pres, rho, gamma
2018 real(wp) :: h
2019
2020 h = 1._wp + (gamma + 1._wp)*pres/rho
2021
2022 end function f_relativistic_enthalpy
2023
2024 !> Mixture coefficients of one state. Under bubbles_euler with num_fluids == 1 the sole advection slot aliases the void fraction
2025 !! (eqn_idx%alf == eqn_idx%adv%end), so alpha is not a composition there and the coefficients are the liquid's. Clipping stays
2026 !! with callers; it differs between solvers and cannot coincide with that case, as mpp_lim requires num_fluids > 1.
2027 subroutine s_compute_mixture_coefficients(alpha_rho_K, alpha_K, rho_K, gamma_K, pi_inf_K, qv_K)
2028
2029
2030# 583 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2031#ifdef _CRAYFTN
2032# 583 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2033#if MFC_OpenACC
2034# 583 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2035!$acc routine seq
2036# 583 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2037#elif MFC_OpenMP
2038# 583 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2039
2040# 583 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2041
2042# 583 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2043!$omp declare target device_type(any)
2044# 583 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2045#else
2046# 583 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2047!DIR$ INLINEALWAYS s_compute_mixture_coefficients
2048# 583 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2049#endif
2050# 583 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2051#elif MFC_OpenACC
2052# 583 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2053!$acc routine seq
2054# 583 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2055#elif MFC_OpenMP
2056# 583 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2057
2058# 583 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2059
2060# 583 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2061!$omp declare target device_type(any)
2062# 583 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2063#endif
2064# 584 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2065
2066# 588 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2067 real(wp), dimension(num_fluids), intent(in) :: alpha_rho_k, alpha_k
2068# 590 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2069 real(wp), intent(out) :: rho_k, gamma_k, pi_inf_k, qv_k
2070 real(wp) :: gamma_i, pi_inf_i, dpi_i, dgamma_i
2071 real(wp) :: rho_i, alpha_i, alpha_rho_i
2072 integer :: i !< Loop iterator over fluids
2073
2074 ! The bubbly closure is written for one carrier liquid, which keeps its own coefficients
2075 ! undiluted: Gamma_l*p_l = (E - rho|u|^2/2)/(1 - alf) - Pi_inf_l, the void entering only through
2076 ! the (1 - alf) that s_compute_pressure applies. There is nothing to sum - the last advection
2077 ! slot is the void, not a material - and the checker holds num_fluids <= 2 here.
2078 if (bubbles_euler) then
2079 rho_k = alpha_rho_k(1)
2080 gamma_k = gammas(1)
2081 pi_inf_k = pi_infs(1)
2082 ! Energy per unit volume, as below: alpha_rho_K(1) is the liquid partial density
2083 qv_k = alpha_rho_k(1)*qvs(1)
2084 else
2085 rho_k = 0._wp
2086 gamma_k = 0._wp
2087 pi_inf_k = 0._wp
2088 qv_k = 0._wp
2089
2090
2091# 611 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2092#if defined(MFC_OpenACC)
2093# 611 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2094!$acc loop seq
2095# 611 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2096#elif defined(MFC_OpenMP)
2097# 611 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2098
2099# 611 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2100#endif
2101# 612 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2102 do i = 1, num_fluids
2103 rho_k = rho_k + alpha_rho_k(i)
2104 alpha_rho_i = alpha_rho_k(i)
2105 alpha_i = alpha_k(i)
2106 call s_phase_coefficients(alpha_rho_i, alpha_i, i, rho_i, gamma_i, pi_inf_i, dpi_i, dgamma_i)
2107 gamma_k = gamma_k + alpha_k(i)*gamma_i
2108 pi_inf_k = pi_inf_k + alpha_k(i)*pi_inf_i
2109 qv_k = qv_k + alpha_rho_k(i)*qvs(i)
2110 end do
2111 end if
2112
2113 end subroutine s_compute_mixture_coefficients
2114
2115 !> Time derivative of the mixture coefficients, mirroring s_compute_mixture_coefficients.
2116 subroutine s_compute_mixture_coefficients_dt(dalpha_rho_dt, dadv_dt, alpha_rho, adv, drho_dt, dgamma_dt, dpi_inf_dt, dqv_dt)
2117
2118
2119# 628 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2120#ifdef _CRAYFTN
2121# 628 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2122#if MFC_OpenACC
2123# 628 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2124!$acc routine seq
2125# 628 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2126#elif MFC_OpenMP
2127# 628 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2128
2129# 628 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2130
2131# 628 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2132!$omp declare target device_type(any)
2133# 628 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2134#else
2135# 628 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2136!DIR$ INLINEALWAYS s_compute_mixture_coefficients_dt
2137# 628 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2138#endif
2139# 628 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2140#elif MFC_OpenACC
2141# 628 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2142!$acc routine seq
2143# 628 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2144#elif MFC_OpenMP
2145# 628 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2146
2147# 628 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2148
2149# 628 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2150!$omp declare target device_type(any)
2151# 628 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2152#endif
2153# 629 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2154
2155# 633 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2156 real(wp), dimension(num_fluids), intent(in) :: dalpha_rho_dt, dadv_dt, alpha_rho, adv
2157# 635 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2158 real(wp), intent(out) :: drho_dt, dgamma_dt, dpi_inf_dt, dqv_dt
2159 real(wp) :: rho_i, gamma_i, pi_inf_i, dpi_i, dgamma_i, alpha_i, alpha_rho_i
2160 integer :: i !< Loop iterator over fluids
2161
2162 dgamma_dt = 0._wp
2163 dpi_inf_dt = 0._wp
2164 dqv_dt = 0._wp
2165
2166 if (num_fluids == 1 .and. bubbles_euler) then
2167 ! Fluid 1's coefficients are constants here, so only rho varies.
2168 drho_dt = dalpha_rho_dt(1)
2169 else
2170 drho_dt = 0._wp
2171
2172
2173# 649 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2174#if defined(MFC_OpenACC)
2175# 649 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2176!$acc loop seq
2177# 649 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2178#elif defined(MFC_OpenMP)
2179# 649 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2180
2181# 649 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2182#endif
2183# 650 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2184 do i = 1, num_fluids
2185 drho_dt = drho_dt + dalpha_rho_dt(i)
2186 alpha_rho_i = alpha_rho(i)
2187 alpha_i = adv(i)
2188 call s_phase_coefficients(alpha_rho_i, alpha_i, i, rho_i, gamma_i, pi_inf_i, dpi_i, dgamma_i)
2189 ! d(alpha X(rho_i))/dt with rho_i = alpha_rho/alpha; the alpha in dX/dt cancels
2190 dgamma_dt = dgamma_dt + dadv_dt(i)*gamma_i + dgamma_i*(dalpha_rho_dt(i) - rho_i*dadv_dt(i))
2191 dpi_inf_dt = dpi_inf_dt + dadv_dt(i)*pi_inf_i + dpi_i*(dalpha_rho_dt(i) - rho_i*dadv_dt(i))
2192 dqv_dt = dqv_dt + dalpha_rho_dt(i)*qvs(i)
2193 end do
2194 end if
2195
2197
2198 !> 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
2199 !! = ((Gamma + 1)p + Pi)/(Gamma rho). Averaged states, whose enthalpy is a free input, use the _avg variant.
2200 subroutine s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c, alpha_rho)
2201
2202
2203# 668 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2204#if MFC_OpenACC
2205# 668 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2206!$acc routine seq
2207# 668 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2208#elif MFC_OpenMP
2209# 668 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2210
2211# 668 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2212
2213# 668 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2214!$omp declare target device_type(any)
2215# 668 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2216#endif
2217# 669 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2218
2219 real(wp), intent(in) :: pres, rho, gamma, pi_inf
2220# 674 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2221 real(wp), dimension(num_fluids), intent(in) :: adv
2222# 676 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2223 real(wp), intent(out) :: c
2224# 680 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2225 real(wp), dimension(num_fluids), intent(in), optional :: alpha_rho
2226# 682 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2227 real(wp) :: alf !< Subgrid void fraction; dilute by construction
2228 real(wp) :: blkmod_q, alpha_q, alpha_rho_q, gamma_q, pi_inf_q
2229 integer :: q
2230
2231 if (chemistry) then ! Reacting mixture sound speed
2232 c = sqrt((1.0_wp + 1.0_wp/gamma)*pres/rho)
2233 else if (relativity) then ! Relativistic sound speed, whose enthalpy is 1 + (Gamma + 1)p/rho
2234 c = sqrt((1._wp + 1._wp/gamma)*pres/rho/f_relativistic_enthalpy(pres, rho, gamma))
2235 else
2236 ! Every case below is a bulk modulus over a density. The equation of state enters
2237 ! only through f_bulk_modulus; the cases differ in how the phases are mixed.
2238 if (any_state_dependent_eos .and. present(alpha_rho)) then ! frozen mixing: each phase's modulus at its own density
2239 c = 0._wp
2240
2241# 695 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2242#if defined(MFC_OpenACC)
2243# 695 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2244!$acc loop seq
2245# 695 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2246#elif defined(MFC_OpenMP)
2247# 695 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2248
2249# 695 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2250#endif
2251# 696 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2252 do q = 1, num_fluids
2253 alpha_q = adv(q)
2254 alpha_rho_q = alpha_rho(q)
2255 call s_phase_bulk_modulus(pres, alpha_q, alpha_rho_q, q, blkmod_q)
2256 if (alt_soundspeed) then
2257 c = c + adv(q)/blkmod_q
2258 else
2259 c = c + adv(q)*blkmod_q
2260 end if
2261 end do
2262 if (alt_soundspeed) then
2263 c = 1._wp/(rho*c)
2264 else
2265 c = c/rho
2266 end if
2267 else if (alt_soundspeed) then ! Wood's law: volume-weighted harmonic mean
2268 c = 0._wp
2269
2270# 713 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2271#if defined(MFC_OpenACC)
2272# 713 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2273!$acc loop seq
2274# 713 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2275#elif defined(MFC_OpenMP)
2276# 713 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2277
2278# 713 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2279#endif
2280# 714 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2281 do q = 1, num_fluids
2282 gamma_q = gammas(q)
2283 pi_inf_q = pi_infs(q)
2284 c = c + adv(q)/f_bulk_modulus(pres, gamma_q, pi_inf_q)
2285 end do
2286 c = 1._wp/(rho*c)
2287 else if (model_eqns == model_eqns_6eq) then ! volume-weighted arithmetic mean
2288 c = 0._wp
2289
2290# 722 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2291#if defined(MFC_OpenACC)
2292# 722 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2293!$acc loop seq
2294# 722 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2295#elif defined(MFC_OpenMP)
2296# 722 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2297
2298# 722 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2299#endif
2300# 723 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2301 do q = 1, num_fluids
2302 gamma_q = gammas(q)
2303 pi_inf_q = pi_infs(q)
2304 c = c + adv(q)*f_bulk_modulus(pres, gamma_q, pi_inf_q)
2305 end do
2306 c = c/rho
2307 else ! the mixture coefficients already carry the mixing
2308 c = f_bulk_modulus(pres, gamma, pi_inf)/rho
2309
2310 ! Subgrid bubbles: c = c_l/(1 - alf), the carrier-liquid speed with an O(alf) void
2311 ! correction. alf is dilute by construction; near one means a wrong index or an
2312 ! out-of-regime case, which the toolchain warns about at case load (#1793).
2313 if (model_eqns == model_eqns_5eq .and. bubbles_euler .and. .not. (mpp_lim .and. num_fluids > 1)) then
2314 alf = adv(num_fluids)
2315 c = c/(1._wp - alf)
2316 end if
2317 end if
2318
2319 if (mixture_err .and. c < 0._wp) then
2320 c = 100._wp*sgm_eps
2321 else
2322 c = sqrt(c)
2323 end if
2324 end if
2325
2326 end subroutine s_compute_speed_of_sound
2327
2328 !> Speed of sound of an interface-averaged state. An average of two states is not a state - its enthalpy is not the one its
2329 !! pressure and density imply - so the caller supplies H, |u|^2 and qv. Only the enthalpy-reading branches differ from
2330 !! s_compute_speed_of_sound; keep the condition below in step with the branch list there.
2331 subroutine s_compute_speed_of_sound_avg(pres, rho, gamma, pi_inf, qv, vel_sum, H, c_c, adv, c, alpha_rho)
2332
2333
2334# 755 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2335#if MFC_OpenACC
2336# 755 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2337!$acc routine seq
2338# 755 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2339#elif MFC_OpenMP
2340# 755 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2341
2342# 755 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2343
2344# 755 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2345!$omp declare target device_type(any)
2346# 755 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2347#endif
2348# 756 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2349
2350 real(wp), intent(in) :: pres, rho, gamma, pi_inf, qv, vel_sum, h, c_c
2351# 761 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2352 real(wp), dimension(num_fluids), intent(in) :: adv
2353# 763 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2354 real(wp), intent(out) :: c
2355# 767 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2356 real(wp), dimension(num_fluids), intent(in), optional :: alpha_rho
2357# 769 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2358
2359 if (chemistry) then ! Reacting mixture sound speed
2360 if (avg_state == avg_state_roe .and. abs(c_c) > verysmall) then
2361 c = sqrt(c_c - (gamma - 1.0_wp)*(vel_sum - h))
2362 else
2363 call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c, alpha_rho)
2364 end if
2365 else if (relativity) then ! Relativistic sound speed
2366 c = sqrt((1._wp + 1._wp/gamma)*pres/rho/h)
2367 else if (alt_soundspeed .or. model_eqns == model_eqns_6eq .or. (model_eqns == model_eqns_5eq .and. bubbles_euler) &
2368 & .or. any_state_dependent_eos) then
2369 call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c, alpha_rho)
2370 else ! Stiffened-gas mixture, the one branch where the averaged enthalpy survives
2371 c = (h - 5.e-1*vel_sum - qv/rho)/gamma
2372
2373 if (mixture_err .and. c < 0._wp) then
2374 c = 100._wp*sgm_eps
2375 else
2376 c = sqrt(c)
2377 end if
2378 end if
2379
2380 end subroutine s_compute_speed_of_sound_avg
2381
2382end module m_eos
Compile-time constant parameters: default values, tolerances, and physical constants.
integer, parameter eos_mie_gruneisen
integer, parameter eos_rk4_steps
Equation of state per fluid (eos_stiffened_gas, eos_ideal_gas, eos_mie_gruneisen, eos_jwl,...
real(wp), parameter sgm_eps
Segmentation tolerance.
real(wp), parameter dflt_real
Default real value.
integer, parameter eos_vinet
integer, parameter eos_ideal_gas
integer, parameter ode_isentrope
integer, parameter ode_reference_temperature
the two ODEs s_rk4 integrates
real(wp), parameter verysmall
Very small number.
integer, parameter eos_jwl
Shared derived types for field data, patch geometry, bubble dynamics, and MPI I/O structures.
Equations of state in Gamma/Pi form, rho e = Gamma(rho) p + Pi(rho).
impure subroutine, public s_finalize_eos_module()
Deallocate the fluid property arrays allocated in s_initialize_eos_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,...
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 ...
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_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_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...
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...
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...
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_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, 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...
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_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.
logical function, public f_is_state_dependent(i)
Whether the EOS of fluid i is a family whose coefficients vary with density.
subroutine s_rk4(kind, i, x0, y0, x1, y)
Fixed-step classical RK4 for the ODE kind from (x0, y0) to x1.
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 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_phase_density(pres, t, i, rho)
Density of phase i at (p, T), the inverse of s_phase_temperature (constant-coefficient EOS only).
impure subroutine, public s_initialize_eos_module()
Resolve every fluid's EOS coefficients once, before any conversion runs.
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.
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 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 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),...
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) function, public f_isentrope_pressure(pi_inf, gamma)
Reference pressure of that isentrope. Precomputed per fluid as isentrope_B.
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_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_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,...
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.
Shared global parameters and equation-index setup for all three executables. Each per-target m_global...
real(wp), dimension(:), allocatable gammas
MPI communication layer: domain decomposition, halo exchange, reductions, and parallel I/O setup.
impure subroutine s_prohibit_abort(condition, message)
Print a case file error with the prohibited condition and message, then abort execution.