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/simulation/generated_eos.fpp" 1
366! AUTO-GENERATED - do not edit directly. Regenerate: cmake reconfigure
367!
368# 4 "/home/runner/work/MFC/MFC/build/include/simulation/generated_eos.fpp"
369
370# 6 "/home/runner/work/MFC/MFC/build/include/simulation/generated_eos.fpp"
371# 9 "/home/runner/work/MFC/MFC/build/include/simulation/generated_eos.fpp"
372
373# 12 "/home/runner/work/MFC/MFC/build/include/simulation/generated_eos.fpp"
374# 15 "/home/runner/work/MFC/MFC/build/include/simulation/generated_eos.fpp"
375
376# 17 "/home/runner/work/MFC/MFC/build/include/simulation/generated_eos.fpp"
377# 31 "/home/runner/work/MFC/MFC/build/include/simulation/generated_eos.fpp"
378
379# 34 "/home/runner/work/MFC/MFC/build/include/simulation/generated_eos.fpp"
380# 58 "/home/runner/work/MFC/MFC/build/include/simulation/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
404 use m_derived_types
405 use m_global_parameters_common
406 use m_mpi_common, only: s_prohibit_abort
407 use m_constants, only: eos_ideal_gas, eos_mie_gruneisen, eos_jwl, eos_vinet, eos_rk4_steps, ode_isentrope, &
408 & ode_reference_temperature, sgm_eps, verysmall, dflt_real
409
410 implicit none
411
412 private
413
414 public :: s_compute_mixture_coefficients, s_compute_mixture_coefficients_dt, s_compute_speed_of_sound, &
415 & s_compute_speed_of_sound_avg, s_initialize_eos_module, s_finalize_eos_module, f_pressure, f_bulk_modulus, &
416 & f_relativistic_enthalpy, f_isentrope_exponent, f_isentrope_pressure, f_sg_thermal, f_mixture_temperature, &
417 & f_is_state_dependent, s_phase_coefficients, s_phase_pressure_on_isentrope, s_phase_temperature, &
418 & s_phase_density_on_isentrope, s_phase_internal_energy, s_phase_bulk_modulus
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 on the isentrope through (rho_from, p_from) at p_to, and c^2 there: Newton on the pressure integrator,
1682 !! whose slope is c^2. The relaxation's own Newton wraps this, so a few steps suffice.
1683 subroutine s_phase_density_on_isentrope(i, rho_from, p_from, p_to, rho_to, c2_to)
1684
1685
1686# 475 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1687#if MFC_OpenACC
1688# 475 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1689!$acc routine seq
1690# 475 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1691#elif MFC_OpenMP
1692# 475 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1693
1694# 475 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1695
1696# 475 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1697!$omp declare target device_type(any)
1698# 475 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1699#endif
1700# 476 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1701
1702 integer, intent(in) :: i
1703 real(wp), intent(in) :: rho_from, p_from, p_to
1704 real(wp), intent(out) :: rho_to, c2_to
1705 real(wp) :: p_at, c2_at
1706 integer :: iter
1707
1708 rho_to = rho_from
1709
1710# 484 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1711#if defined(MFC_OpenACC)
1712# 484 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1713!$acc loop seq
1714# 484 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1715#elif defined(MFC_OpenMP)
1716# 484 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1717
1718# 484 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1719#endif
1720# 485 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1721 do iter = 1, 4
1722 call s_phase_pressure_on_isentrope(p_from, rho_from, rho_to/rho_from, i, p_at)
1723 call s_phase_c2(rho_to, p_at, i, c2_at)
1724 rho_to = rho_to - (p_at - p_to)/c2_at
1725 end do
1726 call s_phase_c2(rho_to, p_to, i, c2_to)
1727
1728 end subroutine s_phase_density_on_isentrope
1729
1730 !> Internal energy per unit volume of phase i at pressure pres: alpha (Gamma p + Pi) + alpha_rho qv, with the coefficients at
1731 !! the phase's own density.
1732 subroutine s_phase_internal_energy(pres, alpha, alpha_rho, i, e_phase)
1733
1734
1735# 498 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1736#ifdef _CRAYFTN
1737# 498 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1738#if MFC_OpenACC
1739# 498 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1740!$acc routine seq
1741# 498 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1742#elif MFC_OpenMP
1743# 498 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1744
1745# 498 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1746
1747# 498 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1748!$omp declare target device_type(any)
1749# 498 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1750#else
1751# 498 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1752!DIR$ INLINEALWAYS s_phase_internal_energy
1753# 498 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1754#endif
1755# 498 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1756#elif MFC_OpenACC
1757# 498 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1758!$acc routine seq
1759# 498 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1760#elif MFC_OpenMP
1761# 498 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1762
1763# 498 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1764
1765# 498 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1766!$omp declare target device_type(any)
1767# 498 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1768#endif
1769# 499 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1770
1771 real(wp), intent(in) :: pres, alpha, alpha_rho
1772 integer, intent(in) :: i
1773 real(wp), intent(out) :: e_phase
1774 real(wp) :: rho, gamma, pi_inf, dpi, dgamma
1775
1776 call s_phase_coefficients(alpha_rho, alpha, i, rho, gamma, pi_inf, dpi, dgamma)
1777 e_phase = alpha*(gamma*pres + pi_inf) + alpha_rho*qvs(i)
1778
1779 end subroutine s_phase_internal_energy
1780
1781 !> Bulk modulus rho c^2 of phase i at pressure pres: f_bulk_modulus for a constant-coefficient fluid, bit for bit, minus the
1782 !! reference-curve terms rho (dPi/drho + p dGamma/drho)/Gamma otherwise.
1783 subroutine s_phase_bulk_modulus(pres, alpha, alpha_rho, i, blkmod)
1784
1785
1786# 514 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1787#ifdef _CRAYFTN
1788# 514 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1789#if MFC_OpenACC
1790# 514 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1791!$acc routine seq
1792# 514 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1793#elif MFC_OpenMP
1794# 514 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1795
1796# 514 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1797
1798# 514 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1799!$omp declare target device_type(any)
1800# 514 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1801#else
1802# 514 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1803!DIR$ INLINEALWAYS s_phase_bulk_modulus
1804# 514 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1805#endif
1806# 514 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1807#elif MFC_OpenACC
1808# 514 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1809!$acc routine seq
1810# 514 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1811#elif MFC_OpenMP
1812# 514 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1813
1814# 514 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1815
1816# 514 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1817!$omp declare target device_type(any)
1818# 514 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1819#endif
1820# 515 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1821
1822 real(wp), intent(in) :: alpha_rho, alpha, pres
1823 integer, intent(in) :: i
1824 real(wp), intent(out) :: blkmod
1825 real(wp) :: rho, gamma, pi_inf, dpi, dgamma
1826
1827 call s_phase_coefficients(alpha_rho, alpha, i, rho, gamma, pi_inf, dpi, dgamma)
1828 blkmod = f_bulk_modulus(pres, gamma, pi_inf) - rho*(dpi + pres*dgamma)/gamma
1829
1830 end subroutine s_phase_bulk_modulus
1831
1832 !> Pressure of a stiffened gas from its internal energy density - the inverse of s_compute_energy. Callers subtract the kinetic,
1833 !! magnetic and elastic energy first; none of those are equation-of-state terms.
1834 function f_pressure(e_int, gamma, pi_inf, qv) result(pres)
1835
1836
1837# 530 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1838#ifdef _CRAYFTN
1839# 530 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1840#if MFC_OpenACC
1841# 530 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1842!$acc routine seq
1843# 530 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1844#elif MFC_OpenMP
1845# 530 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1846
1847# 530 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1848
1849# 530 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1850!$omp declare target device_type(any)
1851# 530 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1852#else
1853# 530 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1854!DIR$ INLINEALWAYS f_pressure
1855# 530 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1856#endif
1857# 530 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1858#elif MFC_OpenACC
1859# 530 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1860!$acc routine seq
1861# 530 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1862#elif MFC_OpenMP
1863# 530 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1864
1865# 530 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1866
1867# 530 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1868!$omp declare target device_type(any)
1869# 530 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1870#endif
1871# 531 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1872
1873 real(wp), intent(in) :: e_int, gamma, pi_inf, qv
1874 real(wp) :: pres
1875
1876 pres = (e_int - pi_inf - qv)/gamma
1877
1878 end function f_pressure
1879
1880 !> Isentropic bulk modulus. Takes coefficients rather than a fluid index, so a mixture - whose effective gamma and pi_inf come
1881 !! from s_compute_mixture_coefficients - is the same call as a single fluid. Elastic callers add their own shear term.
1882 function f_bulk_modulus(pres, gamma, pi_inf) result(blkmod)
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_bulk_modulus
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) :: pres, gamma, pi_inf
1922 real(wp) :: blkmod
1923
1924 blkmod = ((gamma + 1._wp)*pres + pi_inf)/gamma
1925
1926 end function f_bulk_modulus
1927
1928 !> Relativistic specific enthalpy, h = 1 + (Gamma + 1)p/rho. Ideal gas only: the stiffness does not appear, so a fluid with a
1929 !! nonzero pi_inf is not represented here (the validator refuses that combination).
1930 function f_relativistic_enthalpy(pres, rho, gamma) result(H)
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_relativistic_enthalpy
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, rho, gamma
1970 real(wp) :: h
1971
1972 h = 1._wp + (gamma + 1._wp)*pres/rho
1973
1974 end function f_relativistic_enthalpy
1975
1976 !> Mixture coefficients of one state. Under bubbles_euler with num_fluids == 1 the sole advection slot aliases the void fraction
1977 !! (eqn_idx%alf == eqn_idx%adv%end), so alpha is not a composition there and the coefficients are the liquid's. Clipping stays
1978 !! with callers; it differs between solvers and cannot coincide with that case, as mpp_lim requires num_fluids > 1.
1979 subroutine s_compute_mixture_coefficients(alpha_rho_K, alpha_K, rho_K, gamma_K, pi_inf_K, qv_K)
1980
1981
1982# 570 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1983#ifdef _CRAYFTN
1984# 570 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1985#if MFC_OpenACC
1986# 570 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1987!$acc routine seq
1988# 570 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1989#elif MFC_OpenMP
1990# 570 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1991
1992# 570 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1993
1994# 570 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1995!$omp declare target device_type(any)
1996# 570 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1997#else
1998# 570 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
1999!DIR$ INLINEALWAYS s_compute_mixture_coefficients
2000# 570 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2001#endif
2002# 570 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2003#elif MFC_OpenACC
2004# 570 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2005!$acc routine seq
2006# 570 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2007#elif MFC_OpenMP
2008# 570 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2009
2010# 570 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2011
2012# 570 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2013!$omp declare target device_type(any)
2014# 570 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2015#endif
2016# 571 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2017
2018# 575 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2019 real(wp), dimension(num_fluids), intent(in) :: alpha_rho_k, alpha_k
2020# 577 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2021 real(wp), intent(out) :: rho_k, gamma_k, pi_inf_k, qv_k
2022 real(wp) :: gamma_i, pi_inf_i, dpi_i, dgamma_i
2023 real(wp) :: rho_i, alpha_i, alpha_rho_i
2024 integer :: i !< Loop iterator over fluids
2025
2026 ! The bubbly closure is written for one carrier liquid, which keeps its own coefficients
2027 ! undiluted: Gamma_l*p_l = (E - rho|u|^2/2)/(1 - alf) - Pi_inf_l, the void entering only through
2028 ! the (1 - alf) that s_compute_pressure applies. There is nothing to sum - the last advection
2029 ! slot is the void, not a material - and the checker holds num_fluids <= 2 here.
2030 if (bubbles_euler) then
2031 rho_k = alpha_rho_k(1)
2032 gamma_k = gammas(1)
2033 pi_inf_k = pi_infs(1)
2034 ! Energy per unit volume, as below: alpha_rho_K(1) is the liquid partial density
2035 qv_k = alpha_rho_k(1)*qvs(1)
2036 else
2037 rho_k = 0._wp
2038 gamma_k = 0._wp
2039 pi_inf_k = 0._wp
2040 qv_k = 0._wp
2041
2042
2043# 598 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2044#if defined(MFC_OpenACC)
2045# 598 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2046!$acc loop seq
2047# 598 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2048#elif defined(MFC_OpenMP)
2049# 598 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2050
2051# 598 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2052#endif
2053# 599 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2054 do i = 1, num_fluids
2055 rho_k = rho_k + alpha_rho_k(i)
2056 alpha_rho_i = alpha_rho_k(i)
2057 alpha_i = alpha_k(i)
2058 call s_phase_coefficients(alpha_rho_i, alpha_i, i, rho_i, gamma_i, pi_inf_i, dpi_i, dgamma_i)
2059 gamma_k = gamma_k + alpha_k(i)*gamma_i
2060 pi_inf_k = pi_inf_k + alpha_k(i)*pi_inf_i
2061 qv_k = qv_k + alpha_rho_k(i)*qvs(i)
2062 end do
2063 end if
2064
2065 end subroutine s_compute_mixture_coefficients
2066
2067 !> Time derivative of the mixture coefficients, mirroring s_compute_mixture_coefficients.
2068 subroutine s_compute_mixture_coefficients_dt(dalpha_rho_dt, dadv_dt, alpha_rho, adv, drho_dt, dgamma_dt, dpi_inf_dt, dqv_dt)
2069
2070
2071# 615 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2072#ifdef _CRAYFTN
2073# 615 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2074#if MFC_OpenACC
2075# 615 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2076!$acc routine seq
2077# 615 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2078#elif MFC_OpenMP
2079# 615 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2080
2081# 615 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2082
2083# 615 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2084!$omp declare target device_type(any)
2085# 615 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2086#else
2087# 615 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2088!DIR$ INLINEALWAYS s_compute_mixture_coefficients_dt
2089# 615 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2090#endif
2091# 615 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2092#elif MFC_OpenACC
2093# 615 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2094!$acc routine seq
2095# 615 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2096#elif MFC_OpenMP
2097# 615 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2098
2099# 615 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2100
2101# 615 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2102!$omp declare target device_type(any)
2103# 615 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2104#endif
2105# 616 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2106
2107# 620 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2108 real(wp), dimension(num_fluids), intent(in) :: dalpha_rho_dt, dadv_dt, alpha_rho, adv
2109# 622 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2110 real(wp), intent(out) :: drho_dt, dgamma_dt, dpi_inf_dt, dqv_dt
2111 real(wp) :: rho_i, gamma_i, pi_inf_i, dpi_i, dgamma_i, alpha_i, alpha_rho_i
2112 integer :: i !< Loop iterator over fluids
2113
2114 dgamma_dt = 0._wp
2115 dpi_inf_dt = 0._wp
2116 dqv_dt = 0._wp
2117
2118 if (num_fluids == 1 .and. bubbles_euler) then
2119 ! Fluid 1's coefficients are constants here, so only rho varies.
2120 drho_dt = dalpha_rho_dt(1)
2121 else
2122 drho_dt = 0._wp
2123
2124
2125# 636 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2126#if defined(MFC_OpenACC)
2127# 636 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2128!$acc loop seq
2129# 636 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2130#elif defined(MFC_OpenMP)
2131# 636 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2132
2133# 636 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2134#endif
2135# 637 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2136 do i = 1, num_fluids
2137 drho_dt = drho_dt + dalpha_rho_dt(i)
2138 alpha_rho_i = alpha_rho(i)
2139 alpha_i = adv(i)
2140 call s_phase_coefficients(alpha_rho_i, alpha_i, i, rho_i, gamma_i, pi_inf_i, dpi_i, dgamma_i)
2141 ! d(alpha X(rho_i))/dt with rho_i = alpha_rho/alpha; the alpha in dX/dt cancels
2142 dgamma_dt = dgamma_dt + dadv_dt(i)*gamma_i + dgamma_i*(dalpha_rho_dt(i) - rho_i*dadv_dt(i))
2143 dpi_inf_dt = dpi_inf_dt + dadv_dt(i)*pi_inf_i + dpi_i*(dalpha_rho_dt(i) - rho_i*dadv_dt(i))
2144 dqv_dt = dqv_dt + dalpha_rho_dt(i)*qvs(i)
2145 end do
2146 end if
2147
2148 end subroutine s_compute_mixture_coefficients_dt
2149
2150 !> 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
2151 !! = ((Gamma + 1)p + Pi)/(Gamma rho). Averaged states, whose enthalpy is a free input, use the _avg variant.
2152 subroutine s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c, alpha_rho)
2153
2154
2155# 655 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2156#if MFC_OpenACC
2157# 655 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2158!$acc routine seq
2159# 655 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2160#elif MFC_OpenMP
2161# 655 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2162
2163# 655 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2164
2165# 655 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2166!$omp declare target device_type(any)
2167# 655 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2168#endif
2169# 656 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2170
2171 real(wp), intent(in) :: pres, rho, gamma, pi_inf
2172# 661 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2173 real(wp), dimension(num_fluids), intent(in) :: adv
2174# 663 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2175 real(wp), intent(out) :: c
2176# 667 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2177 real(wp), dimension(num_fluids), intent(in), optional :: alpha_rho
2178# 669 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2179 real(wp) :: alf !< Subgrid void fraction; dilute by construction
2180 real(wp) :: blkmod_q, alpha_q, alpha_rho_q, gamma_q, pi_inf_q
2181 integer :: q
2182
2183 if (chemistry) then ! Reacting mixture sound speed
2184 c = sqrt((1.0_wp + 1.0_wp/gamma)*pres/rho)
2185 else if (relativity) then ! Relativistic sound speed, whose enthalpy is 1 + (Gamma + 1)p/rho
2186 c = sqrt((1._wp + 1._wp/gamma)*pres/rho/f_relativistic_enthalpy(pres, rho, gamma))
2187 else
2188 ! Every case below is a bulk modulus over a density. The equation of state enters
2189 ! only through f_bulk_modulus; the cases differ in how the phases are mixed.
2190 if (any_state_dependent_eos .and. present(alpha_rho)) then ! frozen mixing: each phase's modulus at its own density
2191 c = 0._wp
2192
2193# 682 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2194#if defined(MFC_OpenACC)
2195# 682 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2196!$acc loop seq
2197# 682 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2198#elif defined(MFC_OpenMP)
2199# 682 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2200
2201# 682 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2202#endif
2203# 683 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2204 do q = 1, num_fluids
2205 alpha_q = adv(q)
2206 alpha_rho_q = alpha_rho(q)
2207 call s_phase_bulk_modulus(pres, alpha_q, alpha_rho_q, q, blkmod_q)
2208 if (alt_soundspeed) then
2209 c = c + adv(q)/blkmod_q
2210 else
2211 c = c + adv(q)*blkmod_q
2212 end if
2213 end do
2214 if (alt_soundspeed) then
2215 c = 1._wp/(rho*c)
2216 else
2217 c = c/rho
2218 end if
2219 else if (alt_soundspeed) then ! Wood's law: volume-weighted harmonic mean
2220 c = 0._wp
2221
2222# 700 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2223#if defined(MFC_OpenACC)
2224# 700 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2225!$acc loop seq
2226# 700 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2227#elif defined(MFC_OpenMP)
2228# 700 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2229
2230# 700 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2231#endif
2232# 701 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2233 do q = 1, num_fluids
2234 gamma_q = gammas(q)
2235 pi_inf_q = pi_infs(q)
2236 c = c + adv(q)/f_bulk_modulus(pres, gamma_q, pi_inf_q)
2237 end do
2238 c = 1._wp/(rho*c)
2239 else if (model_eqns == model_eqns_6eq) then ! volume-weighted arithmetic mean
2240 c = 0._wp
2241
2242# 709 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2243#if defined(MFC_OpenACC)
2244# 709 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2245!$acc loop seq
2246# 709 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2247#elif defined(MFC_OpenMP)
2248# 709 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2249
2250# 709 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2251#endif
2252# 710 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2253 do q = 1, num_fluids
2254 gamma_q = gammas(q)
2255 pi_inf_q = pi_infs(q)
2256 c = c + adv(q)*f_bulk_modulus(pres, gamma_q, pi_inf_q)
2257 end do
2258 c = c/rho
2259 else ! the mixture coefficients already carry the mixing
2260 c = f_bulk_modulus(pres, gamma, pi_inf)/rho
2261
2262 ! Subgrid bubbles: c = c_l/(1 - alf), the carrier-liquid speed with an O(alf) void
2263 ! correction. alf is dilute by construction; near one means a wrong index or an
2264 ! out-of-regime case, which the toolchain warns about at case load (#1793).
2265 if (model_eqns == model_eqns_5eq .and. bubbles_euler .and. .not. (mpp_lim .and. num_fluids > 1)) then
2266 alf = adv(num_fluids)
2267 c = c/(1._wp - alf)
2268 end if
2269 end if
2270
2271 if (mixture_err .and. c < 0._wp) then
2272 c = 100._wp*sgm_eps
2273 else
2274 c = sqrt(c)
2275 end if
2276 end if
2277
2278 end subroutine s_compute_speed_of_sound
2279
2280 !> Speed of sound of an interface-averaged state. An average of two states is not a state - its enthalpy is not the one its
2281 !! pressure and density imply - so the caller supplies H, |u|^2 and qv. Only the enthalpy-reading branches differ from
2282 !! s_compute_speed_of_sound; keep the condition below in step with the branch list there.
2283 subroutine s_compute_speed_of_sound_avg(pres, rho, gamma, pi_inf, qv, vel_sum, H, c_c, adv, c, alpha_rho)
2284
2285
2286# 742 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2287#if MFC_OpenACC
2288# 742 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2289!$acc routine seq
2290# 742 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2291#elif MFC_OpenMP
2292# 742 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2293
2294# 742 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2295
2296# 742 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2297!$omp declare target device_type(any)
2298# 742 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2299#endif
2300# 743 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2301
2302 real(wp), intent(in) :: pres, rho, gamma, pi_inf, qv, vel_sum, h, c_c
2303# 748 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2304 real(wp), dimension(num_fluids), intent(in) :: adv
2305# 750 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2306 real(wp), intent(out) :: c
2307# 754 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2308 real(wp), dimension(num_fluids), intent(in), optional :: alpha_rho
2309# 756 "/home/runner/work/MFC/MFC/src/common/m_eos.fpp"
2310
2311 if (chemistry) then ! Reacting mixture sound speed
2312 if (avg_state == avg_state_roe .and. abs(c_c) > verysmall) then
2313 c = sqrt(c_c - (gamma - 1.0_wp)*(vel_sum - h))
2314 else
2315 call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c, alpha_rho)
2316 end if
2317 else if (relativity) then ! Relativistic sound speed
2318 c = sqrt((1._wp + 1._wp/gamma)*pres/rho/h)
2319 else if (alt_soundspeed .or. model_eqns == model_eqns_6eq .or. (model_eqns == model_eqns_5eq .and. bubbles_euler) &
2320 & .or. any_state_dependent_eos) then
2321 call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c, alpha_rho)
2322 else ! Stiffened-gas mixture, the one branch where the averaged enthalpy survives
2323 c = (h - 5.e-1*vel_sum - qv/rho)/gamma
2324
2325 if (mixture_err .and. c < 0._wp) then
2326 c = 100._wp*sgm_eps
2327 else
2328 c = sqrt(c)
2329 end if
2330 end if
2331
2332 end subroutine s_compute_speed_of_sound_avg
2333
2334end module m_eos
m_constants
Compile-time constant parameters: default values, tolerances, and physical constants.
Definition m_constants.fpp.f90:7
m_constants::eos_mie_gruneisen
integer, parameter eos_mie_gruneisen
Definition m_constants.fpp.f90:146
m_constants::eos_rk4_steps
integer, parameter eos_rk4_steps
Equation of state per fluid (eos_stiffened_gas, eos_ideal_gas, eos_mie_gruneisen, eos_jwl,...
Definition m_constants.fpp.f90:126
m_constants::sgm_eps
real(wp), parameter sgm_eps
Segmentation tolerance.
Definition m_constants.fpp.f90:13
m_constants::dflt_real
real(wp), parameter dflt_real
Default real value.
Definition m_constants.fpp.f90:12
m_constants::eos_vinet
integer, parameter eos_vinet
Definition m_constants.fpp.f90:148
m_constants::eos_ideal_gas
integer, parameter eos_ideal_gas
Definition m_constants.fpp.f90:145
m_constants::ode_isentrope
integer, parameter ode_isentrope
Definition m_constants.fpp.f90:127
m_constants::ode_reference_temperature
integer, parameter ode_reference_temperature
the two ODEs s_rk4 integrates
Definition m_constants.fpp.f90:127
m_constants::verysmall
real(wp), parameter verysmall
Very small number.
Definition m_constants.fpp.f90:17
m_constants::eos_jwl
integer, parameter eos_jwl
Definition m_constants.fpp.f90:147
m_derived_types
Shared derived types for field data, patch geometry, bubble dynamics, and MPI I/O structures.
Definition m_derived_types.fpp.f90:356
m_eos
Equations of state in Gamma/Pi form, rho e = Gamma(rho) p + Pi(rho).
Definition m_eos.fpp.f90:402
m_eos::s_finalize_eos_module
impure subroutine, public s_finalize_eos_module()
Deallocate the fluid property arrays allocated in s_initialize_eos_module.
Definition m_eos.fpp.f90:827
m_eos::f_hugoniot_compression_limit
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,...
Definition m_eos.fpp.f90:1066
m_eos::s_phase_density_on_isentrope
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 ...
Definition m_eos.fpp.f90:1684
m_eos::f_relativistic_enthalpy
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...
Definition m_eos.fpp.f90:1931
m_eos::s_phase_bulk_modulus
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,...
Definition m_eos.fpp.f90:1784
m_eos::s_phase_pressure_on_isentrope
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...
Definition m_eos.fpp.f90:1567
m_eos::s_reference_curve
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...
Definition m_eos.fpp.f90:866
m_eos::f_mixture_temperature
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...
Definition m_eos.fpp.f90:1275
m_eos::s_compute_speed_of_sound_avg
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...
Definition m_eos.fpp.f90:2284
m_eos::f_sg_thermal
real(wp) function, public f_sg_thermal(pres, rho_or_t, n, b, cv)
Stiffened-gas thermal law p + B = (n - 1)*cv*rho*T. Pass rho to get T, or T to get rho.
Definition m_eos.fpp.f90:1226
m_eos::s_phase_coefficients
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...
Definition m_eos.fpp.f90:1345
m_eos::s_phase_temperature
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.
Definition m_eos.fpp.f90:1627
m_eos::s_compute_mixture_coefficients_dt
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.
Definition m_eos.fpp.f90:2069
m_eos::f_is_state_dependent
logical function, public f_is_state_dependent(i)
Whether the EOS of fluid i is a family whose coefficients vary with density.
Definition m_eos.fpp.f90:966
m_eos::s_rk4
subroutine s_rk4(kind, i, x0, y0, x1, y)
Fixed-step classical RK4 for the ODE kind from (x0, y0) to x1.
Definition m_eos.fpp.f90:1514
m_eos::f_has_isentropic_reference
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),...
Definition m_eos.fpp.f90:1016
m_eos::s_phase_c2
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...
Definition m_eos.fpp.f90:1449
m_eos::s_initialize_eos_module
impure subroutine, public s_initialize_eos_module()
Resolve every fluid's EOS coefficients once, before any conversion runs.
Definition m_eos.fpp.f90:424
m_eos::f_c2_from_coefficients
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.
Definition m_eos.fpp.f90:1401
m_eos::f_bulk_modulus
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...
Definition m_eos.fpp.f90:1883
m_eos::s_eos_coefficients
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 + ...
Definition m_eos.fpp.f90:1092
m_eos::s_ode_slope
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),...
Definition m_eos.fpp.f90:1480
m_eos::s_phase_internal_energy
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,...
Definition m_eos.fpp.f90:1733
m_eos::f_isentrope_pressure
real(wp) function, public f_isentrope_pressure(pi_inf, gamma)
Reference pressure of that isentrope. Precomputed per fluid as isentrope_B.
Definition m_eos.fpp.f90:1179
m_eos::s_compute_mixture_coefficients
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...
Definition m_eos.fpp.f90:1980
m_eos::f_pressure
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....
Definition m_eos.fpp.f90:1835
m_eos::s_compute_speed_of_sound
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,...
Definition m_eos.fpp.f90:2153
m_eos::f_isentrope_exponent
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.
Definition m_eos.fpp.f90:1132
m_global_parameters_common
Shared global parameters and equation-index setup for all three executables. Each per-target m_global...
Definition m_global_parameters_common.fpp.f90:368
m_global_parameters_common::gammas
real(wp), dimension(:), allocatable gammas
Definition m_global_parameters_common.fpp.f90:1726
m_global_parameters_common::num_fluids
integer num_fluids
Definition m_global_parameters_common.fpp.f90:1426
m_mpi_common
MPI communication layer: domain decomposition, halo exchange, reductions, and parallel I/O setup.
Definition m_mpi_common.fpp.f90:367
m_mpi_common::s_prohibit_abort
impure subroutine s_prohibit_abort(condition, message)
Print a case file error with the prohibited condition and message, then abort execution.
Definition m_mpi_common.fpp.f90:780