MFC
Exascale flow solver
Loading...
Searching...
No Matches
m_riemann_solver_hll.fpp.f90
Go to the documentation of this file.
1# 1 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2!>
3!! @file
4!! @brief Contains module m_riemann_solver_hll
5
6!> @brief HLL approximate Riemann solver, Harten et al. SIAM Review (1983)
7# 1 "/home/runner/work/MFC/MFC/src/common/include/case.fpp" 1
8! This file exists so that Fypp can be run without generating case.fpp files for
9! each target. This is useful when generating documentation, for example. This
10! should also let MFC be built with CMake directly, without invoking mfc.sh.
11
12! For pre-process.
13# 8 "/home/runner/work/MFC/MFC/src/common/include/case.fpp"
14
15! For moving immersed boundaries in simulation
16# 12 "/home/runner/work/MFC/MFC/src/common/include/case.fpp"
17# 7 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp" 2
18# 1 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 1
19# 1 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 1
20# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
21# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
22# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
23# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
24# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
25# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
26
27# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
28# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
29# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
30
31# 12 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
32# 1 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp" 1
33! Default for Fypp runs outside the toolchain (e.g. documentation). Searched last, so the
34! build/include/<target>/thermochem.fpp the toolchain writes from m_thermochem takes precedence.
35# 4 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp"
36# 5 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp"
37# 13 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 2
38
39# 16 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
40# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
41
42# 24 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
43
44# 53 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
45
46# 65 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
47
48# 75 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
49
50# 105 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
51
52# 117 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
53
54# 127 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
55
56# 174 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
57! New line at end of file is required for FYPP
58# 2 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
59# 1 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 1
60# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
61# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
62# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
63# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
64# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
65# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
66
67# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
68# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
69# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
70
71# 12 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
72# 1 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp" 1
73! Default for Fypp runs outside the toolchain (e.g. documentation). Searched last, so the
74! build/include/<target>/thermochem.fpp the toolchain writes from m_thermochem takes precedence.
75# 4 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp"
76# 5 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp"
77# 13 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 2
78
79# 16 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
80# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
81
82# 24 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
83
84# 53 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
85
86# 65 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
87
88# 75 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
89
90# 105 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
91
92# 117 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
93
94# 127 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
95
96# 174 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
97! New line at end of file is required for FYPP
98# 2 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 2
99
100# 4 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
101# 5 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
102# 6 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
103# 7 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
104# 8 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
105
106# 20 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
107
108# 43 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
109
110# 48 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
111
112# 53 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
113
114# 58 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
115
116# 63 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
117
118# 68 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
119
120# 81 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
121
122# 95 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
123
124# 100 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
125
126# 105 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
127
128# 110 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
129
130# 115 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
131
132# 120 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
133
134# 125 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
135
136# 130 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
137
138# 135 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
139
140# 165 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
141
142# 219 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
143
144# 233 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
145
146# 258 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
147
148# 269 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
149
150# 271 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
151# 282 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
152
153# 332 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
154
155# 342 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
156
157# 352 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
158
159# 361 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
160
161# 378 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
162
163# 388 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
164
165# 395 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
166
167# 401 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
168
169# 407 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
170
171# 413 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
172
173# 419 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
174
175# 425 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
176! New line at end of file is required for FYPP
177# 3 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
178# 1 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 1
179# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
180# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
181# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
182# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
183# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
184# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
185
186# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
187# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
188# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
189
190# 12 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
191# 1 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp" 1
192! Default for Fypp runs outside the toolchain (e.g. documentation). Searched last, so the
193! build/include/<target>/thermochem.fpp the toolchain writes from m_thermochem takes precedence.
194# 4 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp"
195# 5 "/home/runner/work/MFC/MFC/src/common/include/defaults/thermochem.fpp"
196# 13 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 2
197
198# 16 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
199# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
200
201# 24 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
202
203# 53 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
204
205# 65 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
206
207# 75 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
208
209# 105 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
210
211# 117 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
212
213# 127 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
214
215# 174 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
216! New line at end of file is required for FYPP
217# 2 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 2
218
219# 7 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
220
221# 17 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
222
223# 22 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
224
225# 27 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
226
227# 32 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
228
229# 37 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
230
231# 42 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
232
233# 47 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
234
235# 52 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
236
237# 57 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
238
239# 62 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
240
241# 73 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
242
243# 78 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
244
245# 83 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
246
247# 88 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
248
249# 103 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
250
251# 131 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
252
253# 160 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
254
255# 175 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
256
257# 193 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
258
259# 215 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
260
261# 244 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
262
263# 259 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
264
265# 269 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
266
267# 278 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
268
269# 294 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
270
271# 304 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
272
273# 311 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
274! New line at end of file is required for FYPP
275# 4 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
276
277! GPU parallel region (scalar reductions, maxval/minval)
278# 23 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
279
280! GPU parallel loop over threads (most common GPU macro)
281# 43 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
282
283! Required closing for GPU_PARALLEL_LOOP
284# 55 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
285
286! Mark routine for device compilation
287# 112 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
288
289! Declare device-resident data
290# 130 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
291
292! Inner loop within a GPU parallel region
293# 145 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
294
295! Scoped GPU data region
296# 164 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
297
298! Host code with device pointers (for MPI with GPU buffers)
299# 193 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
300
301! Allocate device memory (unscoped)
302# 207 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
303
304! Free device memory
305# 219 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
306
307! Atomic operation on device
308# 231 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
309
310! End atomic capture block
311# 242 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
312
313! Copy data between host and device
314# 254 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
315
316! Synchronization barrier
317# 266 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
318
319! Import GPU library module (openacc or omp_lib)
320# 275 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
321
322! Emit code only for AMD compiler
323# 282 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
324
325! Emit code for non-Cray compilers
326# 289 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
327
328! Emit code only for Cray compiler
329# 296 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
330
331! Emit code for non-NVIDIA compilers
332# 303 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
333
334# 305 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
335# 306 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
336! New line at end of file is required for FYPP
337# 2 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 2
338
339# 14 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
340
341! Caution: This macro requires the use of a binding script to set CUDA_VISIBLE_DEVICES, such that we have one GPU device per MPI
342! rank. That's because for both cudaMemAdvise (preferred location) and cudaMemPrefetchAsync we use location = device_id = 0. For an
343! example see misc/nvidia_uvm/bind.sh.
344# 52 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
345
346! Allocate and create GPU device memory
347# 72 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
348
349! Free GPU device memory and deallocate
350# 80 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
351
352! Cray-specific GPU pointer setup for vector fields
353# 104 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
354
355! Cray-specific GPU pointer setup for scalar fields
356# 120 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
357
358! Cray-specific GPU pointer setup for acoustic source spatials
359# 145 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
360
361# 151 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
362
363# 158 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
364! New line at end of file is required for FYPP
365# 8 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp" 2
366
367module m_riemann_solver_hll
368
369 use m_derived_types
370 use m_global_parameters
371 use m_variables_conversion
372 use m_eos
373 use m_constants, only: riemann_solver_hll, riemann_solver_hllc, riemann_solver_lax_friedrichs, avg_state_roe, &
374 & avg_state_arithmetic, wave_speeds_direct, wave_speeds_pressure
375 use m_chemistry
376 use m_thermochem, only: gas_constant, get_mixture_molecular_weight, get_species_enthalpies_rt, molecular_weights
377 use m_thermochem_state, only: get_mixture_caloric_state
378 use m_riemann_state
379
380 implicit none
381
382contains
383
384 !> HLL approximate Riemann solver, Harten et al. SIAM Review (1983)
385 subroutine s_hll_riemann_solver(qL_prim_rsx_vf, dqL_prim_dx_vf, dqL_prim_dy_vf, dqL_prim_dz_vf, qL_prim_vf, qR_prim_rsx_vf, &
386 & dqR_prim_dx_vf, dqR_prim_dy_vf, dqR_prim_dz_vf, qR_prim_vf, q_prim_vf, flux_vf, &
387 & flux_src_vf, flux_gsrc_vf, norm_dir, ix, iy, iz)
388
389 real(wp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:), intent(inout) :: qL_prim_rsx_vf, qR_prim_rsx_vf
390 type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
391 type(scalar_field), allocatable, dimension(:), intent(inout) :: qL_prim_vf, qR_prim_vf
392 type(scalar_field), allocatable, dimension(:), intent(inout) :: dqL_prim_dx_vf, dqR_prim_dx_vf, dqL_prim_dy_vf, &
393 & dqR_prim_dy_vf, dqL_prim_dz_vf, dqR_prim_dz_vf
394
395 ! Intercell fluxes
396 type(scalar_field), dimension(sys_size), intent(inout) :: flux_vf, flux_src_vf, flux_gsrc_vf
397 real(wp) :: flux_tau_L, flux_tau_R
398 integer, intent(in) :: norm_dir
399 type(int_bounds_info), intent(in) :: ix, iy, iz
400
401# 48 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
402 real(wp), dimension(num_fluids) :: alpha_rho_L, alpha_rho_R
403 real(wp), dimension(num_vels) :: vel_L, vel_R
404 real(wp), dimension(num_fluids) :: alpha_L, alpha_R
405# 52 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
406 real(wp), dimension(1) :: Ys_L, Ys_R, R_species, h_iL, h_iR
407# 53 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
408 real(wp), dimension(1) :: Cp_iL, Cp_iR, Xs_L, Xs_R, Gamma_iL, Gamma_iR
409# 54 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
410 real(wp) :: rho_L, rho_R
411 real(wp) :: pres_L, pres_R
412 real(wp) :: E_L, E_R
413 real(wp) :: H_L, H_R
414 real(wp) :: c_sum_Yi_Phi
415 real(wp) :: T_L, T_R
416 real(wp) :: Y_L, Y_R
417 real(wp) :: MW_L, MW_R
418 real(wp) :: R_gas_L, R_gas_R
419 real(wp) :: Cp_L, Cp_R
420 real(wp) :: Cv_L, Cv_R
421 real(wp) :: Gamm_L, Gamm_R
422 real(wp) :: gamma_L, gamma_R
423 real(wp) :: pi_inf_L, pi_inf_R
424 real(wp) :: qv_L, qv_R
425 real(wp) :: c_L, c_R
426 real(wp), dimension(6) :: tau_e_L, tau_e_R
427 real(wp) :: G_L, G_R
428 real(wp) :: damage_L, damage_R
429 real(wp) :: solid_partial_density_L, solid_partial_density_R
430 real(wp), dimension(2) :: Re_L, Re_R
431 real(wp) :: rho_avg
432 real(wp) :: H_avg
433 real(wp) :: qv_avg
434 real(wp) :: gamma_avg
435 real(wp) :: c_avg
436 real(wp) :: s_L, s_R, s_M, s_P, s_S
437 real(wp) :: xi_M, xi_P
438 real(wp) :: ptilde_L, ptilde_R
439 real(wp) :: vel_L_rms, vel_R_rms, vel_avg_rms
440 real(wp) :: Ms_L, Ms_R, pres_SL, pres_SR
441 real(wp) :: alpha_L_sum, alpha_R_sum
442 real(wp) :: pcorr !< low Mach number correction
443 type(riemann_states) :: c_fast, pres_mag
444 type(riemann_states_vec3) :: B
445 type(riemann_states) :: Ga !< Gamma (Lorentz factor)
446 type(riemann_states) :: vdotB, B2
447 type(riemann_states_vec3) :: b4 !< 4-magnetic field components (spatial: b4x, b4y, b4z)
448 type(riemann_states_vec3) :: cm !< Conservative momentum variables
449 integer :: i, j, k, l !< Generic loop iterators
450 integer :: Re_size_loc1, Re_size_loc2 !< host copies of Re_size; amdflang reads the declare-target original stale cross-TU
451 ! Populating the buffers of the left and right Riemann problem states variables, based on the choice of boundary conditions
452
453 call s_populate_riemann_states_variables_buffers(ql_prim_rsx_vf, dql_prim_dx_vf, dql_prim_dy_vf, dql_prim_dz_vf, &
454 & qr_prim_rsx_vf, dqr_prim_dx_vf, dqr_prim_dy_vf, dqr_prim_dz_vf, norm_dir, ix, iy, iz)
455
456 ! Reshaping inputted data based on dimensional splitting direction
457 call s_initialize_riemann_solver(flux_src_vf, norm_dir)
458 re_size_loc1 = re_size(1); re_size_loc2 = re_size(2)
459# 107 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
460# 108 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
461# 109 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
462 if (norm_dir == 1) then
463# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
464
465# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
466
467# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
468#if defined(MFC_OpenACC)
469# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
470!$acc parallel loop collapse(3) gang vector default(present) private(i, j, k, l, alpha_rho_L, alpha_rho_R, vel_L, vel_R, alpha_L, alpha_R, tau_e_L, tau_e_R, Re_L, Re_R, s_L, s_R, s_M, s_P, s_S, xi_M, &
471# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
472!$acc& xi_P, Ys_L, Ys_R, R_species, h_iL, h_iR, Cp_iL, Cp_iR, Xs_L, Xs_R, Gamma_iL, Gamma_iR, c_fast, pres_mag, B, Ga, vdotB, B2, b4, cm, pcorr, rho_L, rho_R, pres_L, pres_R, E_L, E_R, H_L, H_R, &
473# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
474!$acc& c_sum_Yi_Phi, T_L, T_R, Y_L, Y_R, MW_L, MW_R, R_gas_L, R_gas_R, Cp_L, Cp_R, Cv_L, Cv_R, Gamm_L, Gamm_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, qv_L, qv_R, qv_avg, c_L, c_R, G_L, G_R, damage_L, &
475# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
476!$acc& damage_R, solid_partial_density_L, solid_partial_density_R, rho_avg, H_avg, c_avg, gamma_avg, ptilde_L, ptilde_R, vel_L_rms, vel_R_rms, vel_avg_rms, Ms_L, Ms_R, pres_SL, pres_SR, alpha_L_sum, &
477# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
478!$acc& alpha_R_sum, flux_tau_L, flux_tau_R) firstprivate(Re_size_loc1, Re_size_loc2) copyin(norm_dir)
479# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
480#elif defined(MFC_OpenMP)
481# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
482
483# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
484
485# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
486
487# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
488!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, j, k, l, &
489# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
490!$omp& alpha_rho_L, alpha_rho_R, vel_L, vel_R, alpha_L, alpha_R, tau_e_L, tau_e_R, Re_L, Re_R, s_L, s_R, s_M, s_P, s_S, xi_M, xi_P, Ys_L, Ys_R, R_species, h_iL, h_iR, Cp_iL, Cp_iR, Xs_L, Xs_R, &
491# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
492!$omp& Gamma_iL, Gamma_iR, c_fast, pres_mag, B, Ga, vdotB, B2, b4, cm, pcorr, rho_L, rho_R, pres_L, pres_R, E_L, E_R, H_L, H_R, c_sum_Yi_Phi, T_L, T_R, Y_L, Y_R, MW_L, MW_R, R_gas_L, R_gas_R, Cp_L, &
493# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
494!$omp& Cp_R, Cv_L, Cv_R, Gamm_L, Gamm_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, qv_L, qv_R, qv_avg, c_L, c_R, G_L, G_R, damage_L, damage_R, solid_partial_density_L, solid_partial_density_R, rho_avg, &
495# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
496!$omp& H_avg, c_avg, gamma_avg, ptilde_L, ptilde_R, vel_L_rms, vel_R_rms, vel_avg_rms, Ms_L, Ms_R, pres_SL, pres_SR, alpha_L_sum, alpha_R_sum, flux_tau_L, flux_tau_R) &
497# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
498!$omp& firstprivate(Re_size_loc1, Re_size_loc2) map(to:norm_dir)
499# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
500#endif
501# 120 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
502 do l = is3%beg, is3%end
503# 121 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
504 do k = is2%beg, is2%end
505# 122 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
506 do j = is1%beg, is1%end
507# 123 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
508
509# 123 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
510#if defined(MFC_OpenACC)
511# 123 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
512!$acc loop seq
513# 123 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
514#elif defined(MFC_OpenMP)
515# 123 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
516
517# 123 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
518#endif
519# 124 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
520 do i = 1, eqn_idx%cont%end
521 alpha_rho_l(i) = ql_prim_rsx_vf(j, k, l, i)
522# 126 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
523 alpha_rho_r(i) = qr_prim_rsx_vf(j + 1, k, l, i)
524# 127 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
525 end do
526
527 vel_l_rms = 0._wp; vel_r_rms = 0._wp
528
529
530# 131 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
531#if defined(MFC_OpenACC)
532# 131 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
533!$acc loop seq
534# 131 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
535#elif defined(MFC_OpenMP)
536# 131 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
537
538# 131 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
539#endif
540# 132 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
541 do i = 1, num_vels
542 vel_l(i) = ql_prim_rsx_vf(j, k, l, eqn_idx%cont%end + i)
543# 134 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
544 vel_r(i) = qr_prim_rsx_vf(j + 1, k, l, eqn_idx%cont%end + i)
545# 135 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
546 vel_l_rms = vel_l_rms + vel_l(i)**2._wp
547 vel_r_rms = vel_r_rms + vel_r(i)**2._wp
548 end do
549
550
551# 139 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
552#if defined(MFC_OpenACC)
553# 139 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
554!$acc loop seq
555# 139 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
556#elif defined(MFC_OpenMP)
557# 139 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
558
559# 139 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
560#endif
561# 140 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
562 do i = 1, num_fluids
563 alpha_l(i) = ql_prim_rsx_vf(j, k, l, eqn_idx%E + i)
564# 142 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
565 alpha_r(i) = qr_prim_rsx_vf(j + 1, k, l, eqn_idx%E + i)
566# 143 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
567 end do
568
569 pres_l = ql_prim_rsx_vf(j, k, l, eqn_idx%E)
570# 146 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
571 pres_r = qr_prim_rsx_vf(j + 1, k, l, eqn_idx%E)
572# 147 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
573
574 if (mhd) then
575 if (n == 0) then ! 1D: constant Bx; By, Bz as variables
576 b%L(1) = bx0
577 b%R(1) = bx0
578 b%L(2) = ql_prim_rsx_vf(j, k, l, eqn_idx%B%beg)
579# 153 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
580 b%R(2) = qr_prim_rsx_vf(j + 1, k, l, eqn_idx%B%beg)
581# 154 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
582 b%L(3) = ql_prim_rsx_vf(j, k, l, eqn_idx%B%beg + 1)
583# 155 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
584 b%R(3) = qr_prim_rsx_vf(j + 1, k, l, eqn_idx%B%beg + 1)
585# 156 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
586 else ! 2D/3D: Bx, By, Bz as variables
587 b%L(1) = ql_prim_rsx_vf(j, k, l, eqn_idx%B%beg)
588# 158 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
589 b%R(1) = qr_prim_rsx_vf(j + 1, k, l, eqn_idx%B%beg)
590# 159 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
591 b%L(2) = ql_prim_rsx_vf(j, k, l, eqn_idx%B%beg + 1)
592# 160 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
593 b%R(2) = qr_prim_rsx_vf(j + 1, k, l, eqn_idx%B%beg + 1)
594# 161 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
595 b%L(3) = ql_prim_rsx_vf(j, k, l, eqn_idx%B%beg + 2)
596# 162 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
597 b%R(3) = qr_prim_rsx_vf(j + 1, k, l, eqn_idx%B%beg + 2)
598# 163 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
599 end if
600 end if
601
602 rho_l = 0._wp
603 gamma_l = 0._wp
604 pi_inf_l = 0._wp
605 qv_l = 0._wp
606
607 rho_r = 0._wp
608 gamma_r = 0._wp
609 pi_inf_r = 0._wp
610 qv_r = 0._wp
611
612 alpha_l_sum = 0._wp
613 alpha_r_sum = 0._wp
614
615 pres_mag%L = 0._wp
616 pres_mag%R = 0._wp
617
618 if (mpp_lim) then
619
620# 183 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
621#if defined(MFC_OpenACC)
622# 183 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
623!$acc loop seq
624# 183 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
625#elif defined(MFC_OpenMP)
626# 183 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
627
628# 183 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
629#endif
630# 184 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
631 do i = 1, num_fluids
632 alpha_rho_l(i) = max(0._wp, alpha_rho_l(i))
633 alpha_l(i) = min(max(0._wp, alpha_l(i)), 1._wp)
634 alpha_l_sum = alpha_l_sum + alpha_l(i)
635 alpha_rho_r(i) = max(0._wp, alpha_rho_r(i))
636 alpha_r(i) = min(max(0._wp, alpha_r(i)), 1._wp)
637 alpha_r_sum = alpha_r_sum + alpha_r(i)
638 end do
639
640 alpha_l = alpha_l/max(alpha_l_sum, sgm_eps)
641 alpha_r = alpha_r/max(alpha_r_sum, sgm_eps)
642 end if
643
644 call s_compute_mixture_coefficients(alpha_rho_l, alpha_l, rho_l, gamma_l, pi_inf_l, qv_l)
645 call s_compute_mixture_coefficients(alpha_rho_r, alpha_r, rho_r, gamma_r, pi_inf_r, qv_r)
646
647 if (viscous) then
648 call s_compute_interface_reynolds(alpha_l, re_l, re_size_loc1, re_size_loc2)
649 call s_compute_interface_reynolds(alpha_r, re_r, re_size_loc1, re_size_loc2)
650 end if
651
652 if (chemistry) then
653 ! Only the Roe path writes this; zero it so the arithmetic path reads nothing undefined.
654 c_sum_yi_phi = 0._wp
655
656# 208 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
657#if defined(MFC_OpenACC)
658# 208 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
659!$acc loop seq
660# 208 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
661#elif defined(MFC_OpenMP)
662# 208 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
663
664# 208 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
665#endif
666# 209 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
667 do i = eqn_idx%species%beg, eqn_idx%species%end
668 ys_l(i - eqn_idx%species%beg + 1) = ql_prim_rsx_vf(j, k, l, i)
669# 211 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
670 ys_r(i - eqn_idx%species%beg + 1) = qr_prim_rsx_vf(j + 1, k, l, i)
671# 212 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
672 end do
673
674 call get_mixture_molecular_weight(ys_l, mw_l)
675 call get_mixture_molecular_weight(ys_r, mw_r)
676 xs_l(1:num_species) = ys_l(1:num_species)*mw_l/molecular_weights(:)
677 xs_r(1:num_species) = ys_r(1:num_species)*mw_r/molecular_weights(:)
678
679 r_gas_l = gas_constant/mw_l
680 r_gas_r = gas_constant/mw_r
681 t_l = pres_l/rho_l/r_gas_l
682 t_r = pres_r/rho_r/r_gas_r
683
684 call get_mixture_caloric_state(t_l, ys_l, cp_il, cp_l, cv_l, e_l)
685 call get_mixture_caloric_state(t_r, ys_r, cp_ir, cp_r, cv_r, e_r)
686
687 if (chem_params%gamma_method == 1) then
688 ! gamma_method = 1: Ref. Section 2.3.1 Formulation of doi:10.7907/ZKW8-ES97.
689 gamma_il(1:num_species) = cp_il(1:num_species)/(cp_il(1:num_species) - 1.0_wp)
690 gamma_ir(1:num_species) = cp_ir(1:num_species)/(cp_ir(1:num_species) - 1.0_wp)
691
692 gamma_l = sum(xs_l(1:num_species)/(gamma_il(1:num_species) - 1.0_wp))
693 gamma_r = sum(xs_r(1:num_species)/(gamma_ir(1:num_species) - 1.0_wp))
694 else if (chem_params%gamma_method == 2) then
695 ! gamma_method = 2: c_p / c_v where c_p, c_v are specific heats.
696 gamm_l = cp_l/cv_l
697 gamma_l = 1.0_wp/(gamm_l - 1.0_wp)
698 gamm_r = cp_r/cv_r
699 gamma_r = 1.0_wp/(gamm_r - 1.0_wp)
700 end if
701
702 e_l = rho_l*e_l + 5.e-1*rho_l*vel_l_rms
703 e_r = rho_r*e_r + 5.e-1*rho_r*vel_r_rms
704 h_l = (e_l + pres_l)/rho_l
705 h_r = (e_r + pres_r)/rho_r
706 else if (mhd .and. relativity) then
707 ga%L = 1._wp/sqrt(1._wp - vel_l_rms)
708 ga%R = 1._wp/sqrt(1._wp - vel_r_rms)
709# 250 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
710 vdotb%L = vel_l(1)*b%L(1) + vel_l(2)*b%L(2) + vel_l(3)*b%L(3)
711 vdotb%R = vel_r(1)*b%R(1) + vel_r(2)*b%R(2) + vel_r(3)*b%R(3)
712
713 b4%L(1:3) = b%L(1:3)/ga%L + ga%L*vel_l(1:3)*vdotb%L
714 b4%R(1:3) = b%R(1:3)/ga%R + ga%R*vel_r(1:3)*vdotb%R
715 b2%L = b%L(1)**2._wp + b%L(2)**2._wp + b%L(3)**2._wp
716 b2%R = b%R(1)**2._wp + b%R(2)**2._wp + b%R(3)**2._wp
717# 258 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
718
719 pres_mag%L = 0.5_wp*(b2%L/ga%L**2._wp + vdotb%L**2._wp)
720 pres_mag%R = 0.5_wp*(b2%R/ga%R**2._wp + vdotb%R**2._wp)
721
722 h_l = f_relativistic_enthalpy(pres_l, rho_l, gamma_l)
723 h_r = f_relativistic_enthalpy(pres_r, rho_r, gamma_r)
724# 265 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
725 cm%L(1:3) = (rho_l*h_l*ga%L**2 + b2%L)*vel_l(1:3) - vdotb%L*b%L(1:3)
726 cm%R(1:3) = (rho_r*h_r*ga%R**2 + b2%R)*vel_r(1:3) - vdotb%R*b%R(1:3)
727# 268 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
728
729 e_l = rho_l*h_l*ga%L**2 - pres_l + 0.5_wp*(b2%L + vel_l_rms*b2%L - vdotb%L**2._wp) - rho_l*ga%L
730 e_r = rho_r*h_r*ga%R**2 - pres_r + 0.5_wp*(b2%R + vel_r_rms*b2%R - vdotb%R**2._wp) - rho_r*ga%R
731 else if (mhd .and. .not. relativity) then
732# 273 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
733 pres_mag%L = 0.5_wp*(b%L(1)**2._wp + b%L(2)**2._wp + b%L(3)**2._wp)
734 pres_mag%R = 0.5_wp*(b%R(1)**2._wp + b%R(2)**2._wp + b%R(3)**2._wp)
735# 276 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
736 call s_compute_energy(pres_l, alpha_rho_l, alpha_l, vel_l_rms, e_l)
737 e_l = e_l + pres_mag%L
738 ! includes magnetic energy
739 call s_compute_energy(pres_r, alpha_rho_r, alpha_r, vel_r_rms, e_r)
740 e_r = e_r + pres_mag%R
741 h_l = (e_l + pres_l - pres_mag%L)/rho_l
742 ! stagnation enthalpy here excludes magnetic energy (only used to find speed of sound)
743 h_r = (e_r + pres_r - pres_mag%R)/rho_r
744 else
745 call s_compute_energy(pres_l, alpha_rho_l, alpha_l, vel_l_rms, e_l)
746 call s_compute_energy(pres_r, alpha_rho_r, alpha_r, vel_r_rms, e_r)
747 h_l = (e_l + pres_l)/rho_l
748 h_r = (e_r + pres_r)/rho_r
749 end if
750
751 ! elastic energy update
752 ! HLL handles hypoelasticity as inline branches like this one - one of three code
753 ! shapes (cf. HLLC inline, and the separate HLLD dual-pass module). See s_riemann_solver.
754 if (hypoelasticity) then
755
756# 295 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
757#if defined(MFC_OpenACC)
758# 295 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
759!$acc loop seq
760# 295 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
761#elif defined(MFC_OpenMP)
762# 295 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
763
764# 295 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
765#endif
766# 296 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
767 do i = 1, eqn_idx%stress%end - eqn_idx%stress%beg + 1
768 tau_e_l(i) = ql_prim_rsx_vf(j, k, l, eqn_idx%stress%beg - 1 + i)
769# 298 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
770 tau_e_r(i) = qr_prim_rsx_vf(j + 1, k, l, eqn_idx%stress%beg - 1 + i)
771# 299 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
772 end do
773
774 damage_l = 0._wp; damage_r = 0._wp
775 if (cont_damage) then
776 damage_l = ql_prim_rsx_vf(j, k, l, eqn_idx%damage)
777# 304 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
778 damage_r = qr_prim_rsx_vf(j + 1, k, l, eqn_idx%damage)
779# 305 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
780 end if
781
782 call s_compute_hypoelastic_interface_energy(num_fluids, alpha_l, alpha_r, damage_l, damage_r, &
783 & tau_e_l, tau_e_r, g_l, g_r, e_l, e_r)
784 end if
785
786 ! Only the pressure-based wave-speed estimate reads the averaged state, and the Roe
787 ! average costs eight square roots per face.
788 if (wave_speeds == wave_speeds_pressure) then
789 call s_compute_average_state(rho_l, rho_r, vel_l, vel_r, h_l, h_r, gamma_l, gamma_r, qv_l, qv_r, &
790 & rho_avg, vel_avg_rms, h_avg, gamma_avg, qv_avg)
791 if (chemistry .and. avg_state == avg_state_roe) then
792 r_species(1:num_species) = gas_constant/molecular_weights
793 call get_species_enthalpies_rt(t_l, h_il)
794 call get_species_enthalpies_rt(t_r, h_ir)
795 h_il(1:num_species) = h_il(1:num_species)*r_species(1:num_species)*t_l
796 h_ir(1:num_species) = h_ir(1:num_species)*r_species(1:num_species)*t_r
797 call s_compute_chemistry_average_state(rho_l, rho_r, t_l, t_r, ys_l, ys_r, r_species, h_il, &
798 & h_ir, cp_il, cp_ir, vel_avg_rms, gamma_avg, &
799 & c_sum_yi_phi)
800 end if
801 end if
802
803 call s_compute_speed_of_sound(pres_l, rho_l, gamma_l, pi_inf_l, alpha_l, c_l, alpha_rho_l)
804
805 call s_compute_speed_of_sound(pres_r, rho_r, gamma_r, pi_inf_r, alpha_r, c_r, alpha_rho_r)
806
807 if (wave_speeds == wave_speeds_pressure) then
808 call s_compute_speed_of_sound_avg(pres_r, rho_avg, gamma_avg, pi_inf_r, qv_avg, vel_avg_rms, &
809 & h_avg, c_sum_yi_phi, alpha_r, c_avg, alpha_rho_r)
810 end if
811
812 if (mhd) then
813 call s_compute_fast_magnetosonic_speed(rho_l, c_l, b%L, norm_dir, c_fast%L, h_l)
814 call s_compute_fast_magnetosonic_speed(rho_r, c_r, b%R, norm_dir, c_fast%R, h_r)
815 end if
816
817 if (viscous) then
818 if (chemistry) then
819 call compute_viscosity_and_inversion(t_l, ys_l, t_r, ys_r, re_l(1), re_r(1))
820 end if
821
822# 346 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
823#if defined(MFC_OpenACC)
824# 346 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
825!$acc loop seq
826# 346 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
827#elif defined(MFC_OpenMP)
828# 346 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
829
830# 346 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
831#endif
832# 347 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
833 do i = 1, 2
834 re_avg_rsx_vf(j, k, l, i) = 2._wp/(1._wp/re_l(i) + 1._wp/re_r(i))
835# 349 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
836 end do
837 end if
838
839 ! Wave speed estimates (wave_speeds=1: direct, wave_speeds=2: pressure-based)
840 if (wave_speeds == wave_speeds_direct) then
841 if (mhd) then
842 ! MHD: use fast magnetosonic speed
843 s_l = min(vel_l(dir_idx(1)) - c_fast%L, vel_r(dir_idx(1)) - c_fast%R)
844 s_r = max(vel_r(dir_idx(1)) + c_fast%R, vel_l(dir_idx(1)) + c_fast%L)
845 else if (hypoelasticity) then
846 ! Elastic wave speed, Rodriguez et al. JCP (2019)
847 s_l = min(vel_l(dir_idx(1)) - f_elastic_signal_speed(c_l, g_l, tau_e_l(dir_idx_tau(1)), &
848 & rho_l), vel_r(dir_idx(1)) - f_elastic_signal_speed(c_r, g_r, &
849 & tau_e_r(dir_idx_tau(1)), rho_r))
850 s_r = max(vel_r(dir_idx(1)) + f_elastic_signal_speed(c_r, g_r, tau_e_r(dir_idx_tau(1)), &
851 & rho_r), vel_l(dir_idx(1)) + f_elastic_signal_speed(c_l, g_l, &
852 & tau_e_l(dir_idx_tau(1)), rho_l))
853 else
854 s_l = min(vel_l(dir_idx(1)) - c_l, vel_r(dir_idx(1)) - c_r)
855 s_r = max(vel_r(dir_idx(1)) + c_r, vel_l(dir_idx(1)) + c_l)
856 end if
857
858 if (hyper_cleaning) then
859 ! Dedner GLM divergence cleaning, Dedner et al. JCP (2002)
860 s_l = min(s_l, -hyper_cleaning_speed)
861 s_r = max(s_r, hyper_cleaning_speed)
862 end if
863
864 s_s = (pres_r - pres_l + rho_l*vel_l(dir_idx(1))*(s_l - vel_l(dir_idx(1))) &
865 & - rho_r*vel_r(dir_idx(1))*(s_r - vel_r(dir_idx(1))))/(rho_l*(s_l - vel_l(dir_idx(1))) &
866 & - rho_r*(s_r - vel_r(dir_idx(1))))
867 else if (wave_speeds == wave_speeds_pressure) then
868 pres_sl = 5.e-1_wp*(pres_l + pres_r + rho_avg*c_avg*(vel_l(dir_idx(1)) - vel_r(dir_idx(1))))
869
870 pres_sr = pres_sl
871
872 ! Low Mach correction: Thornber et al. JCP (2008)
873 ms_l = max(1._wp, &
874 & sqrt(1._wp + 5.e-1_wp*(f_isentrope_exponent(gamma_l) + 1._wp) &
875 & /f_isentrope_exponent(gamma_l)*(pres_sl - pres_l)/(pres_l &
876 & + f_isentrope_pressure(pi_inf_l, gamma_l))))
877 ms_r = max(1._wp, &
878 & sqrt(1._wp + 5.e-1_wp*(f_isentrope_exponent(gamma_r) + 1._wp) &
879 & /f_isentrope_exponent(gamma_r)*(pres_sr - pres_r)/(pres_r &
880 & + f_isentrope_pressure(pi_inf_r, gamma_r))))
881
882 s_l = vel_l(dir_idx(1)) - c_l*ms_l
883 s_r = vel_r(dir_idx(1)) + c_r*ms_r
884
885 s_s = 5.e-1_wp*((vel_l(dir_idx(1)) + vel_r(dir_idx(1))) + (pres_l - pres_r)/(rho_avg*c_avg))
886 end if
887
888 s_m = min(0._wp, s_l); s_p = max(0._wp, s_r)
889
890 xi_m = (5.e-1_wp + sign(5.e-1_wp, s_l)) + (5.e-1_wp - sign(5.e-1_wp, s_l))*(5.e-1_wp + sign(5.e-1_wp, &
891 & s_r))
892 xi_p = (5.e-1_wp - sign(5.e-1_wp, s_r)) + (5.e-1_wp - sign(5.e-1_wp, s_l))*(5.e-1_wp + sign(5.e-1_wp, &
893 & s_r))
894
895 ! HLL intercell flux: F* = (s_R*F_L - s_L*F_R + s_L*s_R*(U_R - U_L)) / (s_R - s_L) Low Mach correction
896 pcorr = f_low_mach_pcorr_hll(vel_l_rms, vel_r_rms, c_l, c_r, rho_l, rho_r, s_m, s_p)
897
898 ! Mass
899 if (.not. relativity) then
900
901# 413 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
902#if defined(MFC_OpenACC)
903# 413 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
904!$acc loop seq
905# 413 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
906#elif defined(MFC_OpenMP)
907# 413 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
908
909# 413 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
910#endif
911# 414 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
912 do i = 1, eqn_idx%cont%end
913 flux_rsx_vf(j, k, l, &
914# 416 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
915 & i) = (s_m*alpha_rho_r(i)*vel_r(norm_dir) - s_p*alpha_rho_l(i)*vel_l(norm_dir) &
916 & + s_m*s_p*(alpha_rho_l(i) - alpha_rho_r(i)))/(s_m - s_p)
917 end do
918 else if (relativity) then
919
920# 420 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
921#if defined(MFC_OpenACC)
922# 420 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
923!$acc loop seq
924# 420 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
925#elif defined(MFC_OpenMP)
926# 420 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
927
928# 420 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
929#endif
930# 421 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
931 do i = 1, eqn_idx%cont%end
932 flux_rsx_vf(j, k, l, &
933# 423 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
934 & i) = (s_m*ga%R*alpha_rho_r(i)*vel_r(norm_dir) - s_p*ga%L*alpha_rho_l(i) &
935 & *vel_l(norm_dir) + s_m*s_p*(ga%L*alpha_rho_l(i) - ga%R*alpha_rho_r(i)))/(s_m &
936 & - s_p)
937 end do
938 end if
939
940 ! Momentum
941 if (mhd .and. (.not. relativity)) then
942
943# 431 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
944#if defined(MFC_OpenACC)
945# 431 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
946!$acc loop seq
947# 431 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
948#elif defined(MFC_OpenMP)
949# 431 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
950
951# 431 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
952#endif
953# 432 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
954 do i = 1, 3
955 ! Flux of rho*v_i in the x direction = rho * v_i * v_x - B_i * B_x +
956# 434 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
957 ! delta_(x,i) * p_tot
958# 435 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
959 flux_rsx_vf(j, k, l, &
960# 436 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
961 & eqn_idx%cont%end + i) = (s_m*(rho_r*vel_r(i)*vel_r(norm_dir) - b%R(i) &
962 & *b%R(norm_dir) + dir_flg(i)*(pres_r + pres_mag%R)) - s_p*(rho_l*vel_l(i) &
963 & *vel_l(norm_dir) - b%L(i)*b%L(norm_dir) + dir_flg(i)*(pres_l + pres_mag%L)) &
964 & + s_m*s_p*(rho_l*vel_l(i) - rho_r*vel_r(i)))/(s_m - s_p)
965 end do
966 else if (mhd .and. relativity) then
967
968# 442 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
969#if defined(MFC_OpenACC)
970# 442 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
971!$acc loop seq
972# 442 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
973#elif defined(MFC_OpenMP)
974# 442 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
975
976# 442 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
977#endif
978# 443 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
979 do i = 1, 3
980 ! Flux of m_i in the x direction = m_i * v_x - b_i/Gamma * B_x +
981# 445 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
982 ! delta_(x,i) * p_tot
983# 446 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
984 flux_rsx_vf(j, k, l, &
985# 447 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
986 & eqn_idx%cont%end + i) = (s_m*(cm%R(i)*vel_r(norm_dir) - b4%R(i) &
987 & /ga%R*b%R(norm_dir) + dir_flg(i)*(pres_r + pres_mag%R)) - s_p*(cm%L(i) &
988 & *vel_l(norm_dir) - b4%L(i)/ga%L*b%L(norm_dir) + dir_flg(i)*(pres_l + pres_mag%L) &
989 & ) + s_m*s_p*(cm%L(i) - cm%R(i)))/(s_m - s_p)
990 end do
991 else if (bubbles_euler) then
992
993# 453 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
994#if defined(MFC_OpenACC)
995# 453 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
996!$acc loop seq
997# 453 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
998#elif defined(MFC_OpenMP)
999# 453 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1000
1001# 453 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1002#endif
1003# 454 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1004 do i = 1, num_vels
1005 flux_rsx_vf(j, k, l, &
1006# 456 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1007 & eqn_idx%cont%end + dir_idx(i)) = (s_m*(rho_r*vel_r(dir_idx(1))*vel_r(dir_idx(i)) &
1008 & + dir_flg(dir_idx(i))*(pres_r - ptilde_r)) - s_p*(rho_l*vel_l(dir_idx(1)) &
1009 & *vel_l(dir_idx(i)) + dir_flg(dir_idx(i))*(pres_l - ptilde_l)) &
1010 & + s_m*s_p*(rho_l*vel_l(dir_idx(i)) - rho_r*vel_r(dir_idx(i))))/(s_m - s_p) &
1011 & + (s_m/s_l)*(s_p/s_r)*pcorr*(vel_r(dir_idx(i)) - vel_l(dir_idx(i)))
1012 end do
1013 else if (hypoelasticity) then
1014
1015# 463 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1016#if defined(MFC_OpenACC)
1017# 463 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1018!$acc loop seq
1019# 463 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1020#elif defined(MFC_OpenMP)
1021# 463 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1022
1023# 463 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1024#endif
1025# 464 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1026 do i = 1, num_vels
1027 flux_rsx_vf(j, k, l, &
1028# 466 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1029 & eqn_idx%cont%end + dir_idx(i)) = (s_m*(rho_r*vel_r(dir_idx(1))*vel_r(dir_idx(i)) &
1030 & + dir_flg(dir_idx(i))*pres_r - tau_e_r(dir_idx_tau(i))) &
1031 & - s_p*(rho_l*vel_l(dir_idx(1))*vel_l(dir_idx(i)) + dir_flg(dir_idx(i))*pres_l &
1032 & - tau_e_l(dir_idx_tau(i))) + s_m*s_p*(rho_l*vel_l(dir_idx(i)) &
1033 & - rho_r*vel_r(dir_idx(i))))/(s_m - s_p)
1034 end do
1035 else
1036
1037# 473 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1038#if defined(MFC_OpenACC)
1039# 473 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1040!$acc loop seq
1041# 473 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1042#elif defined(MFC_OpenMP)
1043# 473 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1044
1045# 473 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1046#endif
1047# 474 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1048 do i = 1, num_vels
1049 flux_rsx_vf(j, k, l, &
1050# 476 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1051 & eqn_idx%cont%end + dir_idx(i)) = (s_m*(rho_r*vel_r(dir_idx(1))*vel_r(dir_idx(i)) &
1052 & + dir_flg(dir_idx(i))*pres_r) - s_p*(rho_l*vel_l(dir_idx(1))*vel_l(dir_idx(i)) &
1053 & + dir_flg(dir_idx(i))*pres_l) + s_m*s_p*(rho_l*vel_l(dir_idx(i)) &
1054 & - rho_r*vel_r(dir_idx(i))))/(s_m - s_p) + (s_m/s_l)*(s_p/s_r) &
1055 & *pcorr*(vel_r(dir_idx(i)) - vel_l(dir_idx(i)))
1056 end do
1057 end if
1058
1059 ! Energy
1060 if (mhd .and. (.not. relativity)) then
1061 ! energy flux = (E + p + p_mag) * v_x - B_x * (v_x*B_x + v_y*B_y + v_z*B_z)
1062# 487 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1063# 488 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1064 flux_rsx_vf(j, k, l, &
1065# 489 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1066 & eqn_idx%E) = (s_m*(vel_r(norm_dir)*(e_r + pres_r + pres_mag%R) - b%R(norm_dir) &
1067 & *(vel_r(1)*b%R(1) + vel_r(2)*b%R(2) + vel_r(3)*b%R(3))) - s_p*(vel_l(norm_dir) &
1068 & *(e_l + pres_l + pres_mag%L) - b%L(norm_dir)*(vel_l(1)*b%L(1) + vel_l(2)*b%L(2) &
1069 & + vel_l(3)*b%L(3))) + s_m*s_p*(e_l - e_r))/(s_m - s_p)
1070# 494 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1071 else if (mhd .and. relativity) then
1072 ! energy flux = m_x - mass flux Hard-coded for single-component for now
1073# 496 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1074 flux_rsx_vf(j, k, l, &
1075# 497 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1076 & eqn_idx%E) = (s_m*(cm%R(norm_dir) - ga%R*alpha_rho_r(1)*vel_r(norm_dir)) &
1077 & - s_p*(cm%L(norm_dir) - ga%L*alpha_rho_l(1)*vel_l(norm_dir)) + s_m*s_p*(e_l - e_r)) &
1078 & /(s_m - s_p)
1079 else if (bubbles_euler) then
1080 flux_rsx_vf(j, k, l, &
1081# 502 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1082 & eqn_idx%E) = (s_m*vel_r(dir_idx(1))*(e_r + pres_r - ptilde_r) - s_p*vel_l(dir_idx(1) &
1083 & )*(e_l + pres_l - ptilde_l) + s_m*s_p*(e_l - e_r))/(s_m - s_p) + (s_m/s_l)*(s_p/s_r) &
1084 & *pcorr*(vel_r_rms - vel_l_rms)/2._wp
1085 else if (hypoelasticity) then
1086 flux_tau_l = 0._wp; flux_tau_r = 0._wp
1087
1088# 507 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1089#if defined(MFC_OpenACC)
1090# 507 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1091!$acc loop seq
1092# 507 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1093#elif defined(MFC_OpenMP)
1094# 507 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1095
1096# 507 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1097#endif
1098# 508 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1099 do i = 1, num_dims
1100 flux_tau_l = flux_tau_l + tau_e_l(dir_idx_tau(i))*vel_l(dir_idx(i))
1101 flux_tau_r = flux_tau_r + tau_e_r(dir_idx_tau(i))*vel_r(dir_idx(i))
1102 end do
1103 flux_rsx_vf(j, k, l, &
1104# 513 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1105 & eqn_idx%E) = (s_m*(vel_r(dir_idx(1))*(e_r + pres_r) - flux_tau_r) &
1106 & - s_p*(vel_l(dir_idx(1))*(e_l + pres_l) - flux_tau_l) + s_m*s_p*(e_l - e_r))/(s_m &
1107 & - s_p)
1108 else
1109 flux_rsx_vf(j, k, l, &
1110# 518 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1111 & eqn_idx%E) = (s_m*vel_r(dir_idx(1))*(e_r + pres_r) - s_p*vel_l(dir_idx(1))*(e_l &
1112 & + pres_l) + s_m*s_p*(e_l - e_r))/(s_m - s_p) + (s_m/s_l)*(s_p/s_r)*pcorr*(vel_r_rms &
1113 & - vel_l_rms)/2._wp
1114 end if
1115
1116 ! Elastic Stresses
1117 if (hypoelasticity) then
1118 do i = 1, eqn_idx%stress%end - eqn_idx%stress%beg + 1 ! TODO: this indexing may be slow
1119 flux_rsx_vf(j, k, l, &
1120# 527 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1121 & eqn_idx%stress%beg - 1 + i) = (s_m*(rho_r*vel_r(dir_idx(1))*tau_e_r(i)) &
1122 & - s_p*(rho_l*vel_l(dir_idx(1))*tau_e_l(i)) + s_m*s_p*(rho_l*tau_e_l(i) &
1123 & - rho_r*tau_e_r(i)))/(s_m - s_p)
1124 end do
1125 end if
1126
1127 ! Damage flux: U_D = m_s*D (damageable-solid partial mass)
1128 if (cont_damage) then
1129 solid_partial_density_l = 0._wp; solid_partial_density_r = 0._wp
1130
1131# 536 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1132#if defined(MFC_OpenACC)
1133# 536 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1134!$acc loop seq
1135# 536 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1136#elif defined(MFC_OpenMP)
1137# 536 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1138
1139# 536 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1140#endif
1141# 537 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1142 do i = 1, num_fluids
1143 if (gs_rs(i) > verysmall) then
1144 solid_partial_density_l = solid_partial_density_l + alpha_rho_l(i)
1145 solid_partial_density_r = solid_partial_density_r + alpha_rho_r(i)
1146 end if
1147 end do
1148 flux_rsx_vf(j, k, l, &
1149# 544 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1150 & eqn_idx%damage) = (s_m*(solid_partial_density_r*vel_r(dir_idx(1))*damage_r) &
1151 & - s_p*(solid_partial_density_l*vel_l(dir_idx(1))*damage_l) &
1152 & + s_m*s_p*(solid_partial_density_l*damage_l - solid_partial_density_r*damage_r)) &
1153 & /(s_m - s_p)
1154 end if
1155
1156 ! Export interface velocity for NC RHS
1157 if (hypo_nc_mode == hypo_nc_mode_interface .or. (alt_soundspeed .and. .not. hll_u_interface)) then
1158
1159# 552 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1160#if defined(MFC_OpenACC)
1161# 552 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1162!$acc loop seq
1163# 552 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1164#elif defined(MFC_OpenMP)
1165# 552 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1166
1167# 552 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1168#endif
1169# 553 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1170 do i = 1, num_dims
1171 if (0._wp <= s_l) then
1172 nc_iface_vel_rsx_vf(j, k, l, dir_idx(i)) = vel_l(dir_idx(i))
1173# 556 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1174 else if (s_r <= 0._wp) then
1175 nc_iface_vel_rsx_vf(j, k, l, dir_idx(i)) = vel_r(dir_idx(i))
1176# 558 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1177 else
1178 nc_iface_vel_rsx_vf(j, k, l, &
1179# 560 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1180 & dir_idx(i)) = (s_r*vel_l(dir_idx(i)) - s_l*vel_r(dir_idx(i)))/(s_r &
1181 & - s_l)
1182 end if
1183 end do
1184 end if
1185
1186 if (.not. hll_u_interface) then ! HLL Method 1: per-fluid alpha interface flux
1187 ! Branchless: s_M/s_P fold the upwinding, so signs of s_L/s_R need no
1188 ! per-thread branches (which diverge across GPU warps)
1189
1190# 569 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1191#if defined(MFC_OpenACC)
1192# 569 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1193!$acc loop seq
1194# 569 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1195#elif defined(MFC_OpenMP)
1196# 569 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1197
1198# 569 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1199#endif
1200# 570 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1201 do i = eqn_idx%adv%beg, eqn_idx%adv%end
1202 flux_rsx_vf(j, k, l, i) = (ql_prim_rsx_vf(j, k, l, i) - qr_prim_rsx_vf(j + 1, k, l, &
1203# 572 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1204 & i))*s_m*s_p/(s_m - s_p)
1205 flux_src_rsx_vf(j, k, l, i) = (s_m*qr_prim_rsx_vf(j + 1, k, l, &
1206# 574 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1207 & i) - s_p*ql_prim_rsx_vf(j, k, l, i))/(s_m - s_p)
1208# 575 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1209 end do
1210 else ! HLL Method 2: shared velocity interface flux
1211
1212# 577 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1213#if defined(MFC_OpenACC)
1214# 577 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1215!$acc loop seq
1216# 577 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1217#elif defined(MFC_OpenMP)
1218# 577 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1219
1220# 577 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1221#endif
1222# 578 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1223 do i = eqn_idx%adv%beg, eqn_idx%adv%end
1224 if (0._wp <= s_l) then
1225 flux_rsx_vf(j, k, l, i) = ql_prim_rsx_vf(j, k, l, i)*vel_l(dir_idx(1))
1226# 581 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1227 else if (s_r <= 0._wp) then
1228 flux_rsx_vf(j, k, l, i) = qr_prim_rsx_vf(j + 1, k, l, i)*vel_r(dir_idx(1))
1229# 583 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1230 else
1231 flux_rsx_vf(j, k, l, i) = (s_r*ql_prim_rsx_vf(j, k, l, &
1232# 585 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1233 & i)*vel_l(dir_idx(1)) - s_l*qr_prim_rsx_vf(j + 1, k, l, &
1234# 586 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1235 & i)*vel_r(dir_idx(1)) + s_l*s_r*(qr_prim_rsx_vf(j + 1, k, l, &
1236# 587 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1237 & i) - ql_prim_rsx_vf(j, k, l, i)))/(s_r - s_l)
1238# 588 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1239 end if
1240 end do
1241 if (0._wp <= s_l) then
1242 flux_src_rsx_vf(j, k, l, eqn_idx%adv%beg) = vel_l(dir_idx(1))
1243# 592 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1244 else if (s_r <= 0._wp) then
1245 flux_src_rsx_vf(j, k, l, eqn_idx%adv%beg) = vel_r(dir_idx(1))
1246# 594 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1247 else
1248 flux_src_rsx_vf(j, k, l, &
1249# 596 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1250 & eqn_idx%adv%beg) = (s_r*vel_l(dir_idx(1)) - s_l*vel_r(dir_idx(1)))/(s_r - s_l)
1251 end if
1252 end if
1253
1254 if (bubbles_euler) then
1255 ! From HLLC: Kills mass transport @ bubble gas density
1256 if (num_fluids > 1) then
1257 flux_rsx_vf(j, k, l, eqn_idx%cont%end) = 0._wp
1258# 604 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1259 end if
1260 end if
1261
1262 if (chemistry) then
1263
1264# 608 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1265#if defined(MFC_OpenACC)
1266# 608 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1267!$acc loop seq
1268# 608 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1269#elif defined(MFC_OpenMP)
1270# 608 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1271
1272# 608 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1273#endif
1274# 609 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1275 do i = eqn_idx%species%beg, eqn_idx%species%end
1276 y_l = ql_prim_rsx_vf(j, k, l, i)
1277# 611 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1278 y_r = qr_prim_rsx_vf(j + 1, k, l, i)
1279# 612 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1280
1281 flux_rsx_vf(j, k, l, &
1282# 614 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1283 & i) = (s_m*y_r*rho_r*vel_r(dir_idx(1)) - s_p*y_l*rho_l*vel_l(dir_idx(1)) &
1284 & + s_m*s_p*(y_l*rho_l - y_r*rho_r))/(s_m - s_p)
1285 flux_src_rsx_vf(j, k, l, i) = 0._wp
1286# 617 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1287 end do
1288 end if
1289
1290 ! MHD: magnetic flux and Maxwell stress contributions
1291 if (mhd) then
1292 if (n == 0) then ! 1D: d/dx flux only & Bx = Bx0 = const.
1293 ! B_y flux = v_x * B_y - v_y * Bx0 B_z flux = v_x * B_z - v_z * Bx0
1294
1295# 624 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1296#if defined(MFC_OpenACC)
1297# 624 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1298!$acc loop seq
1299# 624 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1300#elif defined(MFC_OpenMP)
1301# 624 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1302
1303# 624 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1304#endif
1305# 625 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1306 do i = 0, 1
1307 flux_rsx_vf(j, k, l, &
1308 & eqn_idx%B%beg + i) = (s_m*(vel_r(1)*b%R(2 + i) - vel_r(2 + i)*bx0) &
1309 & - s_p*(vel_l(1)*b%L(2 + i) - vel_l(2 + i)*bx0) + s_m*s_p*(b%L(2 + i) &
1310 & - b%R(2 + i)))/(s_m - s_p)
1311 end do
1312 else ! 2D/3D: Bx, By, Bz /= const. but zero flux component in the same direction
1313 ! B_x d/dx flux = (1 - delta(x,x)) * (v_x * B_x - v_x * B_x) B_y
1314# 633 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1315 ! d/dx flux = (1 - delta(y,x)) * (v_x * B_y - v_y * B_x) B_z d/dx
1316# 634 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1317 ! flux = (1 - delta(z,x)) * (v_x * B_z - v_z * B_x)
1318# 635 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1319
1320# 635 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1321#if defined(MFC_OpenACC)
1322# 635 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1323!$acc loop seq
1324# 635 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1325#elif defined(MFC_OpenMP)
1326# 635 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1327
1328# 635 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1329#endif
1330# 636 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1331 do i = 0, 2
1332 flux_rsx_vf(j, k, l, &
1333# 638 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1334 & eqn_idx%B%beg + i) = (s_m*(vel_r(dir_idx(1))*b%R(i + 1) - vel_r(i + 1) &
1335 & *b%R(norm_dir)) - s_p*(vel_l(dir_idx(1))*b%L(i + 1) - vel_l(i + 1) &
1336 & *b%L(norm_dir)) + s_m*s_p*(b%L(i + 1) - b%R(i + 1)))/(s_m - s_p)
1337 end do
1338
1339 if (hyper_cleaning) then
1340 ! propagate magnetic field divergence as a wave
1341 flux_rsx_vf(j, k, l, eqn_idx%B%beg + norm_dir - 1) = flux_rsx_vf(j, k, l, &
1342# 646 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1343 & eqn_idx%B%beg + norm_dir - 1) + (s_m*qr_prim_rsx_vf(j + 1, k, l, &
1344# 647 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1345 & eqn_idx%psi) - s_p*ql_prim_rsx_vf(j, k, l, eqn_idx%psi))/(s_m - s_p)
1346# 648 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1347
1348 flux_rsx_vf(j, k, l, &
1349# 650 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1350 & eqn_idx%psi) = (hyper_cleaning_speed**2*(s_m*b%R(norm_dir) &
1351 & - s_p*b%L(norm_dir)) + s_m*s_p*(ql_prim_rsx_vf(j, k, l, &
1352# 652 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1353 & eqn_idx%psi) - qr_prim_rsx_vf(j + 1, k, l, eqn_idx%psi)))/(s_m - s_p)
1354# 653 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1355 else
1356 ! Without hyperbolic cleaning, make sure flux of B_normal is identically zero
1357 flux_rsx_vf(j, k, l, eqn_idx%B%beg + norm_dir - 1) = 0._wp
1358# 656 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1359 end if
1360 end if
1361 flux_src_rsx_vf(j, k, l, eqn_idx%adv%beg) = 0._wp
1362# 659 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1363 end if
1364
1365# 698 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1366 end do
1367 end do
1368 end do
1369
1370# 701 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1371#if defined(MFC_OpenACC)
1372# 701 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1373!$acc end parallel loop
1374# 701 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1375#elif defined(MFC_OpenMP)
1376# 701 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1377
1378# 701 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1379!$omp end target teams loop
1380# 701 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1381#endif
1382# 702 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1383 end if
1384# 107 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1385# 108 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1386# 109 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1387 if (norm_dir == 2) then
1388# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1389
1390# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1391
1392# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1393#if defined(MFC_OpenACC)
1394# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1395!$acc parallel loop collapse(3) gang vector default(present) private(i, j, k, l, alpha_rho_L, alpha_rho_R, vel_L, vel_R, alpha_L, alpha_R, tau_e_L, tau_e_R, Re_L, Re_R, s_L, s_R, s_M, s_P, s_S, xi_M, &
1396# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1397!$acc& xi_P, Ys_L, Ys_R, R_species, h_iL, h_iR, Cp_iL, Cp_iR, Xs_L, Xs_R, Gamma_iL, Gamma_iR, c_fast, pres_mag, B, Ga, vdotB, B2, b4, cm, pcorr, rho_L, rho_R, pres_L, pres_R, E_L, E_R, H_L, H_R, &
1398# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1399!$acc& c_sum_Yi_Phi, T_L, T_R, Y_L, Y_R, MW_L, MW_R, R_gas_L, R_gas_R, Cp_L, Cp_R, Cv_L, Cv_R, Gamm_L, Gamm_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, qv_L, qv_R, qv_avg, c_L, c_R, G_L, G_R, damage_L, &
1400# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1401!$acc& damage_R, solid_partial_density_L, solid_partial_density_R, rho_avg, H_avg, c_avg, gamma_avg, ptilde_L, ptilde_R, vel_L_rms, vel_R_rms, vel_avg_rms, Ms_L, Ms_R, pres_SL, pres_SR, alpha_L_sum, &
1402# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1403!$acc& alpha_R_sum, flux_tau_L, flux_tau_R) firstprivate(Re_size_loc1, Re_size_loc2) copyin(norm_dir)
1404# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1405#elif defined(MFC_OpenMP)
1406# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1407
1408# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1409
1410# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1411
1412# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1413!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, j, k, l, &
1414# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1415!$omp& alpha_rho_L, alpha_rho_R, vel_L, vel_R, alpha_L, alpha_R, tau_e_L, tau_e_R, Re_L, Re_R, s_L, s_R, s_M, s_P, s_S, xi_M, xi_P, Ys_L, Ys_R, R_species, h_iL, h_iR, Cp_iL, Cp_iR, Xs_L, Xs_R, &
1416# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1417!$omp& Gamma_iL, Gamma_iR, c_fast, pres_mag, B, Ga, vdotB, B2, b4, cm, pcorr, rho_L, rho_R, pres_L, pres_R, E_L, E_R, H_L, H_R, c_sum_Yi_Phi, T_L, T_R, Y_L, Y_R, MW_L, MW_R, R_gas_L, R_gas_R, Cp_L, &
1418# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1419!$omp& Cp_R, Cv_L, Cv_R, Gamm_L, Gamm_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, qv_L, qv_R, qv_avg, c_L, c_R, G_L, G_R, damage_L, damage_R, solid_partial_density_L, solid_partial_density_R, rho_avg, &
1420# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1421!$omp& H_avg, c_avg, gamma_avg, ptilde_L, ptilde_R, vel_L_rms, vel_R_rms, vel_avg_rms, Ms_L, Ms_R, pres_SL, pres_SR, alpha_L_sum, alpha_R_sum, flux_tau_L, flux_tau_R) &
1422# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1423!$omp& firstprivate(Re_size_loc1, Re_size_loc2) map(to:norm_dir)
1424# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1425#endif
1426# 120 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1427 do l = is3%beg, is3%end
1428# 121 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1429 do k = is1%beg, is1%end
1430# 122 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1431 do j = is2%beg, is2%end
1432# 123 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1433
1434# 123 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1435#if defined(MFC_OpenACC)
1436# 123 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1437!$acc loop seq
1438# 123 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1439#elif defined(MFC_OpenMP)
1440# 123 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1441
1442# 123 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1443#endif
1444# 124 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1445 do i = 1, eqn_idx%cont%end
1446 alpha_rho_l(i) = ql_prim_rsx_vf(j, k, l, i)
1447# 126 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1448 alpha_rho_r(i) = qr_prim_rsx_vf(j, k + 1, l, i)
1449# 127 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1450 end do
1451
1452 vel_l_rms = 0._wp; vel_r_rms = 0._wp
1453
1454
1455# 131 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1456#if defined(MFC_OpenACC)
1457# 131 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1458!$acc loop seq
1459# 131 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1460#elif defined(MFC_OpenMP)
1461# 131 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1462
1463# 131 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1464#endif
1465# 132 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1466 do i = 1, num_vels
1467 vel_l(i) = ql_prim_rsx_vf(j, k, l, eqn_idx%cont%end + i)
1468# 134 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1469 vel_r(i) = qr_prim_rsx_vf(j, k + 1, l, eqn_idx%cont%end + i)
1470# 135 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1471 vel_l_rms = vel_l_rms + vel_l(i)**2._wp
1472 vel_r_rms = vel_r_rms + vel_r(i)**2._wp
1473 end do
1474
1475
1476# 139 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1477#if defined(MFC_OpenACC)
1478# 139 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1479!$acc loop seq
1480# 139 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1481#elif defined(MFC_OpenMP)
1482# 139 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1483
1484# 139 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1485#endif
1486# 140 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1487 do i = 1, num_fluids
1488 alpha_l(i) = ql_prim_rsx_vf(j, k, l, eqn_idx%E + i)
1489# 142 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1490 alpha_r(i) = qr_prim_rsx_vf(j, k + 1, l, eqn_idx%E + i)
1491# 143 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1492 end do
1493
1494 pres_l = ql_prim_rsx_vf(j, k, l, eqn_idx%E)
1495# 146 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1496 pres_r = qr_prim_rsx_vf(j, k + 1, l, eqn_idx%E)
1497# 147 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1498
1499 if (mhd) then
1500 if (n == 0) then ! 1D: constant Bx; By, Bz as variables
1501 b%L(1) = bx0
1502 b%R(1) = bx0
1503 b%L(2) = ql_prim_rsx_vf(j, k, l, eqn_idx%B%beg)
1504# 153 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1505 b%R(2) = qr_prim_rsx_vf(j, k + 1, l, eqn_idx%B%beg)
1506# 154 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1507 b%L(3) = ql_prim_rsx_vf(j, k, l, eqn_idx%B%beg + 1)
1508# 155 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1509 b%R(3) = qr_prim_rsx_vf(j, k + 1, l, eqn_idx%B%beg + 1)
1510# 156 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1511 else ! 2D/3D: Bx, By, Bz as variables
1512 b%L(1) = ql_prim_rsx_vf(j, k, l, eqn_idx%B%beg)
1513# 158 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1514 b%R(1) = qr_prim_rsx_vf(j, k + 1, l, eqn_idx%B%beg)
1515# 159 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1516 b%L(2) = ql_prim_rsx_vf(j, k, l, eqn_idx%B%beg + 1)
1517# 160 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1518 b%R(2) = qr_prim_rsx_vf(j, k + 1, l, eqn_idx%B%beg + 1)
1519# 161 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1520 b%L(3) = ql_prim_rsx_vf(j, k, l, eqn_idx%B%beg + 2)
1521# 162 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1522 b%R(3) = qr_prim_rsx_vf(j, k + 1, l, eqn_idx%B%beg + 2)
1523# 163 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1524 end if
1525 end if
1526
1527 rho_l = 0._wp
1528 gamma_l = 0._wp
1529 pi_inf_l = 0._wp
1530 qv_l = 0._wp
1531
1532 rho_r = 0._wp
1533 gamma_r = 0._wp
1534 pi_inf_r = 0._wp
1535 qv_r = 0._wp
1536
1537 alpha_l_sum = 0._wp
1538 alpha_r_sum = 0._wp
1539
1540 pres_mag%L = 0._wp
1541 pres_mag%R = 0._wp
1542
1543 if (mpp_lim) then
1544
1545# 183 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1546#if defined(MFC_OpenACC)
1547# 183 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1548!$acc loop seq
1549# 183 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1550#elif defined(MFC_OpenMP)
1551# 183 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1552
1553# 183 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1554#endif
1555# 184 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1556 do i = 1, num_fluids
1557 alpha_rho_l(i) = max(0._wp, alpha_rho_l(i))
1558 alpha_l(i) = min(max(0._wp, alpha_l(i)), 1._wp)
1559 alpha_l_sum = alpha_l_sum + alpha_l(i)
1560 alpha_rho_r(i) = max(0._wp, alpha_rho_r(i))
1561 alpha_r(i) = min(max(0._wp, alpha_r(i)), 1._wp)
1562 alpha_r_sum = alpha_r_sum + alpha_r(i)
1563 end do
1564
1565 alpha_l = alpha_l/max(alpha_l_sum, sgm_eps)
1566 alpha_r = alpha_r/max(alpha_r_sum, sgm_eps)
1567 end if
1568
1569 call s_compute_mixture_coefficients(alpha_rho_l, alpha_l, rho_l, gamma_l, pi_inf_l, qv_l)
1570 call s_compute_mixture_coefficients(alpha_rho_r, alpha_r, rho_r, gamma_r, pi_inf_r, qv_r)
1571
1572 if (viscous) then
1573 call s_compute_interface_reynolds(alpha_l, re_l, re_size_loc1, re_size_loc2)
1574 call s_compute_interface_reynolds(alpha_r, re_r, re_size_loc1, re_size_loc2)
1575 end if
1576
1577 if (chemistry) then
1578 ! Only the Roe path writes this; zero it so the arithmetic path reads nothing undefined.
1579 c_sum_yi_phi = 0._wp
1580
1581# 208 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1582#if defined(MFC_OpenACC)
1583# 208 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1584!$acc loop seq
1585# 208 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1586#elif defined(MFC_OpenMP)
1587# 208 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1588
1589# 208 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1590#endif
1591# 209 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1592 do i = eqn_idx%species%beg, eqn_idx%species%end
1593 ys_l(i - eqn_idx%species%beg + 1) = ql_prim_rsx_vf(j, k, l, i)
1594# 211 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1595 ys_r(i - eqn_idx%species%beg + 1) = qr_prim_rsx_vf(j, k + 1, l, i)
1596# 212 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1597 end do
1598
1599 call get_mixture_molecular_weight(ys_l, mw_l)
1600 call get_mixture_molecular_weight(ys_r, mw_r)
1601 xs_l(1:num_species) = ys_l(1:num_species)*mw_l/molecular_weights(:)
1602 xs_r(1:num_species) = ys_r(1:num_species)*mw_r/molecular_weights(:)
1603
1604 r_gas_l = gas_constant/mw_l
1605 r_gas_r = gas_constant/mw_r
1606 t_l = pres_l/rho_l/r_gas_l
1607 t_r = pres_r/rho_r/r_gas_r
1608
1609 call get_mixture_caloric_state(t_l, ys_l, cp_il, cp_l, cv_l, e_l)
1610 call get_mixture_caloric_state(t_r, ys_r, cp_ir, cp_r, cv_r, e_r)
1611
1612 if (chem_params%gamma_method == 1) then
1613 ! gamma_method = 1: Ref. Section 2.3.1 Formulation of doi:10.7907/ZKW8-ES97.
1614 gamma_il(1:num_species) = cp_il(1:num_species)/(cp_il(1:num_species) - 1.0_wp)
1615 gamma_ir(1:num_species) = cp_ir(1:num_species)/(cp_ir(1:num_species) - 1.0_wp)
1616
1617 gamma_l = sum(xs_l(1:num_species)/(gamma_il(1:num_species) - 1.0_wp))
1618 gamma_r = sum(xs_r(1:num_species)/(gamma_ir(1:num_species) - 1.0_wp))
1619 else if (chem_params%gamma_method == 2) then
1620 ! gamma_method = 2: c_p / c_v where c_p, c_v are specific heats.
1621 gamm_l = cp_l/cv_l
1622 gamma_l = 1.0_wp/(gamm_l - 1.0_wp)
1623 gamm_r = cp_r/cv_r
1624 gamma_r = 1.0_wp/(gamm_r - 1.0_wp)
1625 end if
1626
1627 e_l = rho_l*e_l + 5.e-1*rho_l*vel_l_rms
1628 e_r = rho_r*e_r + 5.e-1*rho_r*vel_r_rms
1629 h_l = (e_l + pres_l)/rho_l
1630 h_r = (e_r + pres_r)/rho_r
1631 else if (mhd .and. relativity) then
1632 ga%L = 1._wp/sqrt(1._wp - vel_l_rms)
1633 ga%R = 1._wp/sqrt(1._wp - vel_r_rms)
1634# 250 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1635 vdotb%L = vel_l(1)*b%L(1) + vel_l(2)*b%L(2) + vel_l(3)*b%L(3)
1636 vdotb%R = vel_r(1)*b%R(1) + vel_r(2)*b%R(2) + vel_r(3)*b%R(3)
1637
1638 b4%L(1:3) = b%L(1:3)/ga%L + ga%L*vel_l(1:3)*vdotb%L
1639 b4%R(1:3) = b%R(1:3)/ga%R + ga%R*vel_r(1:3)*vdotb%R
1640 b2%L = b%L(1)**2._wp + b%L(2)**2._wp + b%L(3)**2._wp
1641 b2%R = b%R(1)**2._wp + b%R(2)**2._wp + b%R(3)**2._wp
1642# 258 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1643
1644 pres_mag%L = 0.5_wp*(b2%L/ga%L**2._wp + vdotb%L**2._wp)
1645 pres_mag%R = 0.5_wp*(b2%R/ga%R**2._wp + vdotb%R**2._wp)
1646
1647 h_l = f_relativistic_enthalpy(pres_l, rho_l, gamma_l)
1648 h_r = f_relativistic_enthalpy(pres_r, rho_r, gamma_r)
1649# 265 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1650 cm%L(1:3) = (rho_l*h_l*ga%L**2 + b2%L)*vel_l(1:3) - vdotb%L*b%L(1:3)
1651 cm%R(1:3) = (rho_r*h_r*ga%R**2 + b2%R)*vel_r(1:3) - vdotb%R*b%R(1:3)
1652# 268 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1653
1654 e_l = rho_l*h_l*ga%L**2 - pres_l + 0.5_wp*(b2%L + vel_l_rms*b2%L - vdotb%L**2._wp) - rho_l*ga%L
1655 e_r = rho_r*h_r*ga%R**2 - pres_r + 0.5_wp*(b2%R + vel_r_rms*b2%R - vdotb%R**2._wp) - rho_r*ga%R
1656 else if (mhd .and. .not. relativity) then
1657# 273 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1658 pres_mag%L = 0.5_wp*(b%L(1)**2._wp + b%L(2)**2._wp + b%L(3)**2._wp)
1659 pres_mag%R = 0.5_wp*(b%R(1)**2._wp + b%R(2)**2._wp + b%R(3)**2._wp)
1660# 276 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1661 call s_compute_energy(pres_l, alpha_rho_l, alpha_l, vel_l_rms, e_l)
1662 e_l = e_l + pres_mag%L
1663 ! includes magnetic energy
1664 call s_compute_energy(pres_r, alpha_rho_r, alpha_r, vel_r_rms, e_r)
1665 e_r = e_r + pres_mag%R
1666 h_l = (e_l + pres_l - pres_mag%L)/rho_l
1667 ! stagnation enthalpy here excludes magnetic energy (only used to find speed of sound)
1668 h_r = (e_r + pres_r - pres_mag%R)/rho_r
1669 else
1670 call s_compute_energy(pres_l, alpha_rho_l, alpha_l, vel_l_rms, e_l)
1671 call s_compute_energy(pres_r, alpha_rho_r, alpha_r, vel_r_rms, e_r)
1672 h_l = (e_l + pres_l)/rho_l
1673 h_r = (e_r + pres_r)/rho_r
1674 end if
1675
1676 ! elastic energy update
1677 ! HLL handles hypoelasticity as inline branches like this one - one of three code
1678 ! shapes (cf. HLLC inline, and the separate HLLD dual-pass module). See s_riemann_solver.
1679 if (hypoelasticity) then
1680
1681# 295 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1682#if defined(MFC_OpenACC)
1683# 295 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1684!$acc loop seq
1685# 295 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1686#elif defined(MFC_OpenMP)
1687# 295 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1688
1689# 295 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1690#endif
1691# 296 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1692 do i = 1, eqn_idx%stress%end - eqn_idx%stress%beg + 1
1693 tau_e_l(i) = ql_prim_rsx_vf(j, k, l, eqn_idx%stress%beg - 1 + i)
1694# 298 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1695 tau_e_r(i) = qr_prim_rsx_vf(j, k + 1, l, eqn_idx%stress%beg - 1 + i)
1696# 299 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1697 end do
1698
1699 damage_l = 0._wp; damage_r = 0._wp
1700 if (cont_damage) then
1701 damage_l = ql_prim_rsx_vf(j, k, l, eqn_idx%damage)
1702# 304 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1703 damage_r = qr_prim_rsx_vf(j, k + 1, l, eqn_idx%damage)
1704# 305 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1705 end if
1706
1707 call s_compute_hypoelastic_interface_energy(num_fluids, alpha_l, alpha_r, damage_l, damage_r, &
1708 & tau_e_l, tau_e_r, g_l, g_r, e_l, e_r)
1709 end if
1710
1711 ! Only the pressure-based wave-speed estimate reads the averaged state, and the Roe
1712 ! average costs eight square roots per face.
1713 if (wave_speeds == wave_speeds_pressure) then
1714 call s_compute_average_state(rho_l, rho_r, vel_l, vel_r, h_l, h_r, gamma_l, gamma_r, qv_l, qv_r, &
1715 & rho_avg, vel_avg_rms, h_avg, gamma_avg, qv_avg)
1716 if (chemistry .and. avg_state == avg_state_roe) then
1717 r_species(1:num_species) = gas_constant/molecular_weights
1718 call get_species_enthalpies_rt(t_l, h_il)
1719 call get_species_enthalpies_rt(t_r, h_ir)
1720 h_il(1:num_species) = h_il(1:num_species)*r_species(1:num_species)*t_l
1721 h_ir(1:num_species) = h_ir(1:num_species)*r_species(1:num_species)*t_r
1722 call s_compute_chemistry_average_state(rho_l, rho_r, t_l, t_r, ys_l, ys_r, r_species, h_il, &
1723 & h_ir, cp_il, cp_ir, vel_avg_rms, gamma_avg, &
1724 & c_sum_yi_phi)
1725 end if
1726 end if
1727
1728 call s_compute_speed_of_sound(pres_l, rho_l, gamma_l, pi_inf_l, alpha_l, c_l, alpha_rho_l)
1729
1730 call s_compute_speed_of_sound(pres_r, rho_r, gamma_r, pi_inf_r, alpha_r, c_r, alpha_rho_r)
1731
1732 if (wave_speeds == wave_speeds_pressure) then
1733 call s_compute_speed_of_sound_avg(pres_r, rho_avg, gamma_avg, pi_inf_r, qv_avg, vel_avg_rms, &
1734 & h_avg, c_sum_yi_phi, alpha_r, c_avg, alpha_rho_r)
1735 end if
1736
1737 if (mhd) then
1738 call s_compute_fast_magnetosonic_speed(rho_l, c_l, b%L, norm_dir, c_fast%L, h_l)
1739 call s_compute_fast_magnetosonic_speed(rho_r, c_r, b%R, norm_dir, c_fast%R, h_r)
1740 end if
1741
1742 if (viscous) then
1743 if (chemistry) then
1744 call compute_viscosity_and_inversion(t_l, ys_l, t_r, ys_r, re_l(1), re_r(1))
1745 end if
1746
1747# 346 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1748#if defined(MFC_OpenACC)
1749# 346 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1750!$acc loop seq
1751# 346 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1752#elif defined(MFC_OpenMP)
1753# 346 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1754
1755# 346 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1756#endif
1757# 347 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1758 do i = 1, 2
1759 re_avg_rsx_vf(j, k, l, i) = 2._wp/(1._wp/re_l(i) + 1._wp/re_r(i))
1760# 349 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1761 end do
1762 end if
1763
1764 ! Wave speed estimates (wave_speeds=1: direct, wave_speeds=2: pressure-based)
1765 if (wave_speeds == wave_speeds_direct) then
1766 if (mhd) then
1767 ! MHD: use fast magnetosonic speed
1768 s_l = min(vel_l(dir_idx(1)) - c_fast%L, vel_r(dir_idx(1)) - c_fast%R)
1769 s_r = max(vel_r(dir_idx(1)) + c_fast%R, vel_l(dir_idx(1)) + c_fast%L)
1770 else if (hypoelasticity) then
1771 ! Elastic wave speed, Rodriguez et al. JCP (2019)
1772 s_l = min(vel_l(dir_idx(1)) - f_elastic_signal_speed(c_l, g_l, tau_e_l(dir_idx_tau(1)), &
1773 & rho_l), vel_r(dir_idx(1)) - f_elastic_signal_speed(c_r, g_r, &
1774 & tau_e_r(dir_idx_tau(1)), rho_r))
1775 s_r = max(vel_r(dir_idx(1)) + f_elastic_signal_speed(c_r, g_r, tau_e_r(dir_idx_tau(1)), &
1776 & rho_r), vel_l(dir_idx(1)) + f_elastic_signal_speed(c_l, g_l, &
1777 & tau_e_l(dir_idx_tau(1)), rho_l))
1778 else
1779 s_l = min(vel_l(dir_idx(1)) - c_l, vel_r(dir_idx(1)) - c_r)
1780 s_r = max(vel_r(dir_idx(1)) + c_r, vel_l(dir_idx(1)) + c_l)
1781 end if
1782
1783 if (hyper_cleaning) then
1784 ! Dedner GLM divergence cleaning, Dedner et al. JCP (2002)
1785 s_l = min(s_l, -hyper_cleaning_speed)
1786 s_r = max(s_r, hyper_cleaning_speed)
1787 end if
1788
1789 s_s = (pres_r - pres_l + rho_l*vel_l(dir_idx(1))*(s_l - vel_l(dir_idx(1))) &
1790 & - rho_r*vel_r(dir_idx(1))*(s_r - vel_r(dir_idx(1))))/(rho_l*(s_l - vel_l(dir_idx(1))) &
1791 & - rho_r*(s_r - vel_r(dir_idx(1))))
1792 else if (wave_speeds == wave_speeds_pressure) then
1793 pres_sl = 5.e-1_wp*(pres_l + pres_r + rho_avg*c_avg*(vel_l(dir_idx(1)) - vel_r(dir_idx(1))))
1794
1795 pres_sr = pres_sl
1796
1797 ! Low Mach correction: Thornber et al. JCP (2008)
1798 ms_l = max(1._wp, &
1799 & sqrt(1._wp + 5.e-1_wp*(f_isentrope_exponent(gamma_l) + 1._wp) &
1800 & /f_isentrope_exponent(gamma_l)*(pres_sl - pres_l)/(pres_l &
1801 & + f_isentrope_pressure(pi_inf_l, gamma_l))))
1802 ms_r = max(1._wp, &
1803 & sqrt(1._wp + 5.e-1_wp*(f_isentrope_exponent(gamma_r) + 1._wp) &
1804 & /f_isentrope_exponent(gamma_r)*(pres_sr - pres_r)/(pres_r &
1805 & + f_isentrope_pressure(pi_inf_r, gamma_r))))
1806
1807 s_l = vel_l(dir_idx(1)) - c_l*ms_l
1808 s_r = vel_r(dir_idx(1)) + c_r*ms_r
1809
1810 s_s = 5.e-1_wp*((vel_l(dir_idx(1)) + vel_r(dir_idx(1))) + (pres_l - pres_r)/(rho_avg*c_avg))
1811 end if
1812
1813 s_m = min(0._wp, s_l); s_p = max(0._wp, s_r)
1814
1815 xi_m = (5.e-1_wp + sign(5.e-1_wp, s_l)) + (5.e-1_wp - sign(5.e-1_wp, s_l))*(5.e-1_wp + sign(5.e-1_wp, &
1816 & s_r))
1817 xi_p = (5.e-1_wp - sign(5.e-1_wp, s_r)) + (5.e-1_wp - sign(5.e-1_wp, s_l))*(5.e-1_wp + sign(5.e-1_wp, &
1818 & s_r))
1819
1820 ! HLL intercell flux: F* = (s_R*F_L - s_L*F_R + s_L*s_R*(U_R - U_L)) / (s_R - s_L) Low Mach correction
1821 pcorr = f_low_mach_pcorr_hll(vel_l_rms, vel_r_rms, c_l, c_r, rho_l, rho_r, s_m, s_p)
1822
1823 ! Mass
1824 if (.not. relativity) then
1825
1826# 413 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1827#if defined(MFC_OpenACC)
1828# 413 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1829!$acc loop seq
1830# 413 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1831#elif defined(MFC_OpenMP)
1832# 413 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1833
1834# 413 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1835#endif
1836# 414 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1837 do i = 1, eqn_idx%cont%end
1838 flux_rsx_vf(j, k, l, &
1839# 416 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1840 & i) = (s_m*alpha_rho_r(i)*vel_r(norm_dir) - s_p*alpha_rho_l(i)*vel_l(norm_dir) &
1841 & + s_m*s_p*(alpha_rho_l(i) - alpha_rho_r(i)))/(s_m - s_p)
1842 end do
1843 else if (relativity) then
1844
1845# 420 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1846#if defined(MFC_OpenACC)
1847# 420 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1848!$acc loop seq
1849# 420 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1850#elif defined(MFC_OpenMP)
1851# 420 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1852
1853# 420 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1854#endif
1855# 421 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1856 do i = 1, eqn_idx%cont%end
1857 flux_rsx_vf(j, k, l, &
1858# 423 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1859 & i) = (s_m*ga%R*alpha_rho_r(i)*vel_r(norm_dir) - s_p*ga%L*alpha_rho_l(i) &
1860 & *vel_l(norm_dir) + s_m*s_p*(ga%L*alpha_rho_l(i) - ga%R*alpha_rho_r(i)))/(s_m &
1861 & - s_p)
1862 end do
1863 end if
1864
1865 ! Momentum
1866 if (mhd .and. (.not. relativity)) then
1867
1868# 431 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1869#if defined(MFC_OpenACC)
1870# 431 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1871!$acc loop seq
1872# 431 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1873#elif defined(MFC_OpenMP)
1874# 431 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1875
1876# 431 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1877#endif
1878# 432 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1879 do i = 1, 3
1880 ! Flux of rho*v_i in the y direction = rho * v_i * v_y - B_i * B_y +
1881# 434 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1882 ! delta_(y,i) * p_tot
1883# 435 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1884 flux_rsx_vf(j, k, l, &
1885# 436 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1886 & eqn_idx%cont%end + i) = (s_m*(rho_r*vel_r(i)*vel_r(norm_dir) - b%R(i) &
1887 & *b%R(norm_dir) + dir_flg(i)*(pres_r + pres_mag%R)) - s_p*(rho_l*vel_l(i) &
1888 & *vel_l(norm_dir) - b%L(i)*b%L(norm_dir) + dir_flg(i)*(pres_l + pres_mag%L)) &
1889 & + s_m*s_p*(rho_l*vel_l(i) - rho_r*vel_r(i)))/(s_m - s_p)
1890 end do
1891 else if (mhd .and. relativity) then
1892
1893# 442 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1894#if defined(MFC_OpenACC)
1895# 442 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1896!$acc loop seq
1897# 442 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1898#elif defined(MFC_OpenMP)
1899# 442 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1900
1901# 442 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1902#endif
1903# 443 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1904 do i = 1, 3
1905 ! Flux of m_i in the y direction = m_i * v_y - b_i/Gamma * B_y +
1906# 445 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1907 ! delta_(y,i) * p_tot
1908# 446 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1909 flux_rsx_vf(j, k, l, &
1910# 447 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1911 & eqn_idx%cont%end + i) = (s_m*(cm%R(i)*vel_r(norm_dir) - b4%R(i) &
1912 & /ga%R*b%R(norm_dir) + dir_flg(i)*(pres_r + pres_mag%R)) - s_p*(cm%L(i) &
1913 & *vel_l(norm_dir) - b4%L(i)/ga%L*b%L(norm_dir) + dir_flg(i)*(pres_l + pres_mag%L) &
1914 & ) + s_m*s_p*(cm%L(i) - cm%R(i)))/(s_m - s_p)
1915 end do
1916 else if (bubbles_euler) then
1917
1918# 453 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1919#if defined(MFC_OpenACC)
1920# 453 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1921!$acc loop seq
1922# 453 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1923#elif defined(MFC_OpenMP)
1924# 453 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1925
1926# 453 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1927#endif
1928# 454 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1929 do i = 1, num_vels
1930 flux_rsx_vf(j, k, l, &
1931# 456 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1932 & eqn_idx%cont%end + dir_idx(i)) = (s_m*(rho_r*vel_r(dir_idx(1))*vel_r(dir_idx(i)) &
1933 & + dir_flg(dir_idx(i))*(pres_r - ptilde_r)) - s_p*(rho_l*vel_l(dir_idx(1)) &
1934 & *vel_l(dir_idx(i)) + dir_flg(dir_idx(i))*(pres_l - ptilde_l)) &
1935 & + s_m*s_p*(rho_l*vel_l(dir_idx(i)) - rho_r*vel_r(dir_idx(i))))/(s_m - s_p) &
1936 & + (s_m/s_l)*(s_p/s_r)*pcorr*(vel_r(dir_idx(i)) - vel_l(dir_idx(i)))
1937 end do
1938 else if (hypoelasticity) then
1939
1940# 463 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1941#if defined(MFC_OpenACC)
1942# 463 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1943!$acc loop seq
1944# 463 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1945#elif defined(MFC_OpenMP)
1946# 463 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1947
1948# 463 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1949#endif
1950# 464 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1951 do i = 1, num_vels
1952 flux_rsx_vf(j, k, l, &
1953# 466 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1954 & eqn_idx%cont%end + dir_idx(i)) = (s_m*(rho_r*vel_r(dir_idx(1))*vel_r(dir_idx(i)) &
1955 & + dir_flg(dir_idx(i))*pres_r - tau_e_r(dir_idx_tau(i))) &
1956 & - s_p*(rho_l*vel_l(dir_idx(1))*vel_l(dir_idx(i)) + dir_flg(dir_idx(i))*pres_l &
1957 & - tau_e_l(dir_idx_tau(i))) + s_m*s_p*(rho_l*vel_l(dir_idx(i)) &
1958 & - rho_r*vel_r(dir_idx(i))))/(s_m - s_p)
1959 end do
1960 else
1961
1962# 473 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1963#if defined(MFC_OpenACC)
1964# 473 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1965!$acc loop seq
1966# 473 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1967#elif defined(MFC_OpenMP)
1968# 473 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1969
1970# 473 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1971#endif
1972# 474 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1973 do i = 1, num_vels
1974 flux_rsx_vf(j, k, l, &
1975# 476 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1976 & eqn_idx%cont%end + dir_idx(i)) = (s_m*(rho_r*vel_r(dir_idx(1))*vel_r(dir_idx(i)) &
1977 & + dir_flg(dir_idx(i))*pres_r) - s_p*(rho_l*vel_l(dir_idx(1))*vel_l(dir_idx(i)) &
1978 & + dir_flg(dir_idx(i))*pres_l) + s_m*s_p*(rho_l*vel_l(dir_idx(i)) &
1979 & - rho_r*vel_r(dir_idx(i))))/(s_m - s_p) + (s_m/s_l)*(s_p/s_r) &
1980 & *pcorr*(vel_r(dir_idx(i)) - vel_l(dir_idx(i)))
1981 end do
1982 end if
1983
1984 ! Energy
1985 if (mhd .and. (.not. relativity)) then
1986 ! energy flux = (E + p + p_mag) * v_y - B_y * (v_x*B_x + v_y*B_y + v_z*B_z)
1987# 487 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1988# 488 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1989 flux_rsx_vf(j, k, l, &
1990# 489 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1991 & eqn_idx%E) = (s_m*(vel_r(norm_dir)*(e_r + pres_r + pres_mag%R) - b%R(norm_dir) &
1992 & *(vel_r(1)*b%R(1) + vel_r(2)*b%R(2) + vel_r(3)*b%R(3))) - s_p*(vel_l(norm_dir) &
1993 & *(e_l + pres_l + pres_mag%L) - b%L(norm_dir)*(vel_l(1)*b%L(1) + vel_l(2)*b%L(2) &
1994 & + vel_l(3)*b%L(3))) + s_m*s_p*(e_l - e_r))/(s_m - s_p)
1995# 494 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1996 else if (mhd .and. relativity) then
1997 ! energy flux = m_y - mass flux Hard-coded for single-component for now
1998# 496 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
1999 flux_rsx_vf(j, k, l, &
2000# 497 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2001 & eqn_idx%E) = (s_m*(cm%R(norm_dir) - ga%R*alpha_rho_r(1)*vel_r(norm_dir)) &
2002 & - s_p*(cm%L(norm_dir) - ga%L*alpha_rho_l(1)*vel_l(norm_dir)) + s_m*s_p*(e_l - e_r)) &
2003 & /(s_m - s_p)
2004 else if (bubbles_euler) then
2005 flux_rsx_vf(j, k, l, &
2006# 502 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2007 & eqn_idx%E) = (s_m*vel_r(dir_idx(1))*(e_r + pres_r - ptilde_r) - s_p*vel_l(dir_idx(1) &
2008 & )*(e_l + pres_l - ptilde_l) + s_m*s_p*(e_l - e_r))/(s_m - s_p) + (s_m/s_l)*(s_p/s_r) &
2009 & *pcorr*(vel_r_rms - vel_l_rms)/2._wp
2010 else if (hypoelasticity) then
2011 flux_tau_l = 0._wp; flux_tau_r = 0._wp
2012
2013# 507 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2014#if defined(MFC_OpenACC)
2015# 507 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2016!$acc loop seq
2017# 507 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2018#elif defined(MFC_OpenMP)
2019# 507 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2020
2021# 507 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2022#endif
2023# 508 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2024 do i = 1, num_dims
2025 flux_tau_l = flux_tau_l + tau_e_l(dir_idx_tau(i))*vel_l(dir_idx(i))
2026 flux_tau_r = flux_tau_r + tau_e_r(dir_idx_tau(i))*vel_r(dir_idx(i))
2027 end do
2028 flux_rsx_vf(j, k, l, &
2029# 513 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2030 & eqn_idx%E) = (s_m*(vel_r(dir_idx(1))*(e_r + pres_r) - flux_tau_r) &
2031 & - s_p*(vel_l(dir_idx(1))*(e_l + pres_l) - flux_tau_l) + s_m*s_p*(e_l - e_r))/(s_m &
2032 & - s_p)
2033 else
2034 flux_rsx_vf(j, k, l, &
2035# 518 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2036 & eqn_idx%E) = (s_m*vel_r(dir_idx(1))*(e_r + pres_r) - s_p*vel_l(dir_idx(1))*(e_l &
2037 & + pres_l) + s_m*s_p*(e_l - e_r))/(s_m - s_p) + (s_m/s_l)*(s_p/s_r)*pcorr*(vel_r_rms &
2038 & - vel_l_rms)/2._wp
2039 end if
2040
2041 ! Elastic Stresses
2042 if (hypoelasticity) then
2043 do i = 1, eqn_idx%stress%end - eqn_idx%stress%beg + 1 ! TODO: this indexing may be slow
2044 flux_rsx_vf(j, k, l, &
2045# 527 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2046 & eqn_idx%stress%beg - 1 + i) = (s_m*(rho_r*vel_r(dir_idx(1))*tau_e_r(i)) &
2047 & - s_p*(rho_l*vel_l(dir_idx(1))*tau_e_l(i)) + s_m*s_p*(rho_l*tau_e_l(i) &
2048 & - rho_r*tau_e_r(i)))/(s_m - s_p)
2049 end do
2050 end if
2051
2052 ! Damage flux: U_D = m_s*D (damageable-solid partial mass)
2053 if (cont_damage) then
2054 solid_partial_density_l = 0._wp; solid_partial_density_r = 0._wp
2055
2056# 536 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2057#if defined(MFC_OpenACC)
2058# 536 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2059!$acc loop seq
2060# 536 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2061#elif defined(MFC_OpenMP)
2062# 536 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2063
2064# 536 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2065#endif
2066# 537 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2067 do i = 1, num_fluids
2068 if (gs_rs(i) > verysmall) then
2069 solid_partial_density_l = solid_partial_density_l + alpha_rho_l(i)
2070 solid_partial_density_r = solid_partial_density_r + alpha_rho_r(i)
2071 end if
2072 end do
2073 flux_rsx_vf(j, k, l, &
2074# 544 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2075 & eqn_idx%damage) = (s_m*(solid_partial_density_r*vel_r(dir_idx(1))*damage_r) &
2076 & - s_p*(solid_partial_density_l*vel_l(dir_idx(1))*damage_l) &
2077 & + s_m*s_p*(solid_partial_density_l*damage_l - solid_partial_density_r*damage_r)) &
2078 & /(s_m - s_p)
2079 end if
2080
2081 ! Export interface velocity for NC RHS
2082 if (hypo_nc_mode == hypo_nc_mode_interface .or. (alt_soundspeed .and. .not. hll_u_interface)) then
2083
2084# 552 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2085#if defined(MFC_OpenACC)
2086# 552 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2087!$acc loop seq
2088# 552 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2089#elif defined(MFC_OpenMP)
2090# 552 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2091
2092# 552 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2093#endif
2094# 553 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2095 do i = 1, num_dims
2096 if (0._wp <= s_l) then
2097 nc_iface_vel_rsx_vf(j, k, l, dir_idx(i)) = vel_l(dir_idx(i))
2098# 556 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2099 else if (s_r <= 0._wp) then
2100 nc_iface_vel_rsx_vf(j, k, l, dir_idx(i)) = vel_r(dir_idx(i))
2101# 558 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2102 else
2103 nc_iface_vel_rsx_vf(j, k, l, &
2104# 560 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2105 & dir_idx(i)) = (s_r*vel_l(dir_idx(i)) - s_l*vel_r(dir_idx(i)))/(s_r &
2106 & - s_l)
2107 end if
2108 end do
2109 end if
2110
2111 if (.not. hll_u_interface) then ! HLL Method 1: per-fluid alpha interface flux
2112 ! Branchless: s_M/s_P fold the upwinding, so signs of s_L/s_R need no
2113 ! per-thread branches (which diverge across GPU warps)
2114
2115# 569 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2116#if defined(MFC_OpenACC)
2117# 569 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2118!$acc loop seq
2119# 569 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2120#elif defined(MFC_OpenMP)
2121# 569 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2122
2123# 569 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2124#endif
2125# 570 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2126 do i = eqn_idx%adv%beg, eqn_idx%adv%end
2127 flux_rsx_vf(j, k, l, i) = (ql_prim_rsx_vf(j, k, l, i) - qr_prim_rsx_vf(j, k + 1, l, &
2128# 572 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2129 & i))*s_m*s_p/(s_m - s_p)
2130 flux_src_rsx_vf(j, k, l, i) = (s_m*qr_prim_rsx_vf(j, k + 1, l, &
2131# 574 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2132 & i) - s_p*ql_prim_rsx_vf(j, k, l, i))/(s_m - s_p)
2133# 575 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2134 end do
2135 else ! HLL Method 2: shared velocity interface flux
2136
2137# 577 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2138#if defined(MFC_OpenACC)
2139# 577 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2140!$acc loop seq
2141# 577 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2142#elif defined(MFC_OpenMP)
2143# 577 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2144
2145# 577 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2146#endif
2147# 578 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2148 do i = eqn_idx%adv%beg, eqn_idx%adv%end
2149 if (0._wp <= s_l) then
2150 flux_rsx_vf(j, k, l, i) = ql_prim_rsx_vf(j, k, l, i)*vel_l(dir_idx(1))
2151# 581 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2152 else if (s_r <= 0._wp) then
2153 flux_rsx_vf(j, k, l, i) = qr_prim_rsx_vf(j, k + 1, l, i)*vel_r(dir_idx(1))
2154# 583 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2155 else
2156 flux_rsx_vf(j, k, l, i) = (s_r*ql_prim_rsx_vf(j, k, l, &
2157# 585 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2158 & i)*vel_l(dir_idx(1)) - s_l*qr_prim_rsx_vf(j, k + 1, l, &
2159# 586 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2160 & i)*vel_r(dir_idx(1)) + s_l*s_r*(qr_prim_rsx_vf(j, k + 1, l, &
2161# 587 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2162 & i) - ql_prim_rsx_vf(j, k, l, i)))/(s_r - s_l)
2163# 588 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2164 end if
2165 end do
2166 if (0._wp <= s_l) then
2167 flux_src_rsx_vf(j, k, l, eqn_idx%adv%beg) = vel_l(dir_idx(1))
2168# 592 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2169 else if (s_r <= 0._wp) then
2170 flux_src_rsx_vf(j, k, l, eqn_idx%adv%beg) = vel_r(dir_idx(1))
2171# 594 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2172 else
2173 flux_src_rsx_vf(j, k, l, &
2174# 596 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2175 & eqn_idx%adv%beg) = (s_r*vel_l(dir_idx(1)) - s_l*vel_r(dir_idx(1)))/(s_r - s_l)
2176 end if
2177 end if
2178
2179 if (bubbles_euler) then
2180 ! From HLLC: Kills mass transport @ bubble gas density
2181 if (num_fluids > 1) then
2182 flux_rsx_vf(j, k, l, eqn_idx%cont%end) = 0._wp
2183# 604 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2184 end if
2185 end if
2186
2187 if (chemistry) then
2188
2189# 608 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2190#if defined(MFC_OpenACC)
2191# 608 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2192!$acc loop seq
2193# 608 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2194#elif defined(MFC_OpenMP)
2195# 608 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2196
2197# 608 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2198#endif
2199# 609 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2200 do i = eqn_idx%species%beg, eqn_idx%species%end
2201 y_l = ql_prim_rsx_vf(j, k, l, i)
2202# 611 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2203 y_r = qr_prim_rsx_vf(j, k + 1, l, i)
2204# 612 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2205
2206 flux_rsx_vf(j, k, l, &
2207# 614 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2208 & i) = (s_m*y_r*rho_r*vel_r(dir_idx(1)) - s_p*y_l*rho_l*vel_l(dir_idx(1)) &
2209 & + s_m*s_p*(y_l*rho_l - y_r*rho_r))/(s_m - s_p)
2210 flux_src_rsx_vf(j, k, l, i) = 0._wp
2211# 617 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2212 end do
2213 end if
2214
2215 ! MHD: magnetic flux and Maxwell stress contributions
2216 if (mhd) then
2217 if (n == 0) then ! 1D: d/dx flux only & Bx = Bx0 = const.
2218 ! B_y flux = v_x * B_y - v_y * Bx0 B_z flux = v_x * B_z - v_z * Bx0
2219
2220# 624 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2221#if defined(MFC_OpenACC)
2222# 624 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2223!$acc loop seq
2224# 624 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2225#elif defined(MFC_OpenMP)
2226# 624 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2227
2228# 624 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2229#endif
2230# 625 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2231 do i = 0, 1
2232 flux_rsx_vf(j, k, l, &
2233 & eqn_idx%B%beg + i) = (s_m*(vel_r(1)*b%R(2 + i) - vel_r(2 + i)*bx0) &
2234 & - s_p*(vel_l(1)*b%L(2 + i) - vel_l(2 + i)*bx0) + s_m*s_p*(b%L(2 + i) &
2235 & - b%R(2 + i)))/(s_m - s_p)
2236 end do
2237 else ! 2D/3D: Bx, By, Bz /= const. but zero flux component in the same direction
2238 ! B_x d/dy flux = (1 - delta(x,y)) * (v_y * B_x - v_x * B_y) B_y
2239# 633 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2240 ! d/dy flux = (1 - delta(y,y)) * (v_y * B_y - v_y * B_y) B_z d/dy
2241# 634 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2242 ! flux = (1 - delta(z,y)) * (v_y * B_z - v_z * B_y)
2243# 635 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2244
2245# 635 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2246#if defined(MFC_OpenACC)
2247# 635 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2248!$acc loop seq
2249# 635 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2250#elif defined(MFC_OpenMP)
2251# 635 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2252
2253# 635 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2254#endif
2255# 636 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2256 do i = 0, 2
2257 flux_rsx_vf(j, k, l, &
2258# 638 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2259 & eqn_idx%B%beg + i) = (s_m*(vel_r(dir_idx(1))*b%R(i + 1) - vel_r(i + 1) &
2260 & *b%R(norm_dir)) - s_p*(vel_l(dir_idx(1))*b%L(i + 1) - vel_l(i + 1) &
2261 & *b%L(norm_dir)) + s_m*s_p*(b%L(i + 1) - b%R(i + 1)))/(s_m - s_p)
2262 end do
2263
2264 if (hyper_cleaning) then
2265 ! propagate magnetic field divergence as a wave
2266 flux_rsx_vf(j, k, l, eqn_idx%B%beg + norm_dir - 1) = flux_rsx_vf(j, k, l, &
2267# 646 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2268 & eqn_idx%B%beg + norm_dir - 1) + (s_m*qr_prim_rsx_vf(j, k + 1, l, &
2269# 647 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2270 & eqn_idx%psi) - s_p*ql_prim_rsx_vf(j, k, l, eqn_idx%psi))/(s_m - s_p)
2271# 648 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2272
2273 flux_rsx_vf(j, k, l, &
2274# 650 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2275 & eqn_idx%psi) = (hyper_cleaning_speed**2*(s_m*b%R(norm_dir) &
2276 & - s_p*b%L(norm_dir)) + s_m*s_p*(ql_prim_rsx_vf(j, k, l, &
2277# 652 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2278 & eqn_idx%psi) - qr_prim_rsx_vf(j, k + 1, l, eqn_idx%psi)))/(s_m - s_p)
2279# 653 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2280 else
2281 ! Without hyperbolic cleaning, make sure flux of B_normal is identically zero
2282 flux_rsx_vf(j, k, l, eqn_idx%B%beg + norm_dir - 1) = 0._wp
2283# 656 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2284 end if
2285 end if
2286 flux_src_rsx_vf(j, k, l, eqn_idx%adv%beg) = 0._wp
2287# 659 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2288 end if
2289
2290# 662 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2291 if (cyl_coord) then
2292 ! Substituting the advective flux into the inviscid geometrical source flux
2293
2294# 664 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2295#if defined(MFC_OpenACC)
2296# 664 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2297!$acc loop seq
2298# 664 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2299#elif defined(MFC_OpenMP)
2300# 664 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2301
2302# 664 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2303#endif
2304# 665 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2305 do i = 1, eqn_idx%E
2306 flux_gsrc_rsx_vf(j, k, l, i) = flux_rsx_vf(j, k, l, i)
2307# 667 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2308 end do
2309 ! Recalculating the radial momentum geometric source flux
2310 flux_gsrc_rsx_vf(j, k, l, eqn_idx%cont%end + 2) = flux_rsx_vf(j, k, l, &
2311# 670 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2312 & eqn_idx%cont%end + 2) - (s_m*pres_r - s_p*pres_l)/(s_m - s_p)
2313 ! Geometrical source of the void fraction(s) is zero
2314
2315# 672 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2316#if defined(MFC_OpenACC)
2317# 672 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2318!$acc loop seq
2319# 672 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2320#elif defined(MFC_OpenMP)
2321# 672 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2322
2323# 672 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2324#endif
2325# 673 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2326 do i = eqn_idx%adv%beg, eqn_idx%adv%end
2327 if (.not. hll_u_interface) then
2328 flux_gsrc_rsx_vf(j, k, l, i) = flux_rsx_vf(j, k, l, i)
2329# 676 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2330 else
2331 flux_gsrc_rsx_vf(j, k, l, i) = 0._wp
2332# 678 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2333 end if
2334 end do
2335 end if
2336
2337 if (cyl_coord .and. hypoelasticity) then
2338 ! += tau_sigmasigma using HLL
2339 flux_gsrc_rsx_vf(j, k, l, eqn_idx%cont%end + 2) = flux_gsrc_rsx_vf(j, k, l, &
2340# 685 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2341 & eqn_idx%cont%end + 2) + (s_m*tau_e_r(4) - s_p*tau_e_l(4))/(s_m - s_p)
2342
2343
2344# 687 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2345#if defined(MFC_OpenACC)
2346# 687 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2347!$acc loop seq
2348# 687 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2349#elif defined(MFC_OpenMP)
2350# 687 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2351
2352# 687 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2353#endif
2354# 688 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2355 do i = eqn_idx%stress%beg, eqn_idx%stress%end
2356 flux_gsrc_rsx_vf(j, k, l, i) = flux_rsx_vf(j, k, l, i)
2357# 690 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2358 end do
2359
2360 ! Damage advects conservatively, so its geometric source is the advective flux
2361 if (cont_damage) then
2362 flux_gsrc_rsx_vf(j, k, l, eqn_idx%damage) = flux_rsx_vf(j, k, l, eqn_idx%damage)
2363# 695 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2364 end if
2365 end if
2366# 698 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2367 end do
2368 end do
2369 end do
2370
2371# 701 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2372#if defined(MFC_OpenACC)
2373# 701 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2374!$acc end parallel loop
2375# 701 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2376#elif defined(MFC_OpenMP)
2377# 701 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2378
2379# 701 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2380!$omp end target teams loop
2381# 701 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2382#endif
2383# 702 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2384 end if
2385# 107 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2386# 108 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2387# 109 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2388 if (norm_dir == 3) then
2389# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2390
2391# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2392
2393# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2394#if defined(MFC_OpenACC)
2395# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2396!$acc parallel loop collapse(3) gang vector default(present) private(i, j, k, l, alpha_rho_L, alpha_rho_R, vel_L, vel_R, alpha_L, alpha_R, tau_e_L, tau_e_R, Re_L, Re_R, s_L, s_R, s_M, s_P, s_S, xi_M, &
2397# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2398!$acc& xi_P, Ys_L, Ys_R, R_species, h_iL, h_iR, Cp_iL, Cp_iR, Xs_L, Xs_R, Gamma_iL, Gamma_iR, c_fast, pres_mag, B, Ga, vdotB, B2, b4, cm, pcorr, rho_L, rho_R, pres_L, pres_R, E_L, E_R, H_L, H_R, &
2399# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2400!$acc& c_sum_Yi_Phi, T_L, T_R, Y_L, Y_R, MW_L, MW_R, R_gas_L, R_gas_R, Cp_L, Cp_R, Cv_L, Cv_R, Gamm_L, Gamm_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, qv_L, qv_R, qv_avg, c_L, c_R, G_L, G_R, damage_L, &
2401# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2402!$acc& damage_R, solid_partial_density_L, solid_partial_density_R, rho_avg, H_avg, c_avg, gamma_avg, ptilde_L, ptilde_R, vel_L_rms, vel_R_rms, vel_avg_rms, Ms_L, Ms_R, pres_SL, pres_SR, alpha_L_sum, &
2403# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2404!$acc& alpha_R_sum, flux_tau_L, flux_tau_R) firstprivate(Re_size_loc1, Re_size_loc2) copyin(norm_dir)
2405# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2406#elif defined(MFC_OpenMP)
2407# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2408
2409# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2410
2411# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2412
2413# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2414!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(i, j, k, l, &
2415# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2416!$omp& alpha_rho_L, alpha_rho_R, vel_L, vel_R, alpha_L, alpha_R, tau_e_L, tau_e_R, Re_L, Re_R, s_L, s_R, s_M, s_P, s_S, xi_M, xi_P, Ys_L, Ys_R, R_species, h_iL, h_iR, Cp_iL, Cp_iR, Xs_L, Xs_R, &
2417# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2418!$omp& Gamma_iL, Gamma_iR, c_fast, pres_mag, B, Ga, vdotB, B2, b4, cm, pcorr, rho_L, rho_R, pres_L, pres_R, E_L, E_R, H_L, H_R, c_sum_Yi_Phi, T_L, T_R, Y_L, Y_R, MW_L, MW_R, R_gas_L, R_gas_R, Cp_L, &
2419# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2420!$omp& Cp_R, Cv_L, Cv_R, Gamm_L, Gamm_R, gamma_L, gamma_R, pi_inf_L, pi_inf_R, qv_L, qv_R, qv_avg, c_L, c_R, G_L, G_R, damage_L, damage_R, solid_partial_density_L, solid_partial_density_R, rho_avg, &
2421# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2422!$omp& H_avg, c_avg, gamma_avg, ptilde_L, ptilde_R, vel_L_rms, vel_R_rms, vel_avg_rms, Ms_L, Ms_R, pres_SL, pres_SR, alpha_L_sum, alpha_R_sum, flux_tau_L, flux_tau_R) &
2423# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2424!$omp& firstprivate(Re_size_loc1, Re_size_loc2) map(to:norm_dir)
2425# 110 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2426#endif
2427# 120 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2428 do l = is1%beg, is1%end
2429# 121 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2430 do k = is2%beg, is2%end
2431# 122 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2432 do j = is3%beg, is3%end
2433# 123 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2434
2435# 123 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2436#if defined(MFC_OpenACC)
2437# 123 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2438!$acc loop seq
2439# 123 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2440#elif defined(MFC_OpenMP)
2441# 123 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2442
2443# 123 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2444#endif
2445# 124 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2446 do i = 1, eqn_idx%cont%end
2447 alpha_rho_l(i) = ql_prim_rsx_vf(j, k, l, i)
2448# 126 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2449 alpha_rho_r(i) = qr_prim_rsx_vf(j, k, l + 1, i)
2450# 127 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2451 end do
2452
2453 vel_l_rms = 0._wp; vel_r_rms = 0._wp
2454
2455
2456# 131 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2457#if defined(MFC_OpenACC)
2458# 131 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2459!$acc loop seq
2460# 131 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2461#elif defined(MFC_OpenMP)
2462# 131 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2463
2464# 131 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2465#endif
2466# 132 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2467 do i = 1, num_vels
2468 vel_l(i) = ql_prim_rsx_vf(j, k, l, eqn_idx%cont%end + i)
2469# 134 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2470 vel_r(i) = qr_prim_rsx_vf(j, k, l + 1, eqn_idx%cont%end + i)
2471# 135 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2472 vel_l_rms = vel_l_rms + vel_l(i)**2._wp
2473 vel_r_rms = vel_r_rms + vel_r(i)**2._wp
2474 end do
2475
2476
2477# 139 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2478#if defined(MFC_OpenACC)
2479# 139 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2480!$acc loop seq
2481# 139 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2482#elif defined(MFC_OpenMP)
2483# 139 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2484
2485# 139 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2486#endif
2487# 140 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2488 do i = 1, num_fluids
2489 alpha_l(i) = ql_prim_rsx_vf(j, k, l, eqn_idx%E + i)
2490# 142 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2491 alpha_r(i) = qr_prim_rsx_vf(j, k, l + 1, eqn_idx%E + i)
2492# 143 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2493 end do
2494
2495 pres_l = ql_prim_rsx_vf(j, k, l, eqn_idx%E)
2496# 146 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2497 pres_r = qr_prim_rsx_vf(j, k, l + 1, eqn_idx%E)
2498# 147 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2499
2500 if (mhd) then
2501 if (n == 0) then ! 1D: constant Bx; By, Bz as variables
2502 b%L(1) = bx0
2503 b%R(1) = bx0
2504 b%L(2) = ql_prim_rsx_vf(j, k, l, eqn_idx%B%beg)
2505# 153 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2506 b%R(2) = qr_prim_rsx_vf(j, k, l + 1, eqn_idx%B%beg)
2507# 154 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2508 b%L(3) = ql_prim_rsx_vf(j, k, l, eqn_idx%B%beg + 1)
2509# 155 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2510 b%R(3) = qr_prim_rsx_vf(j, k, l + 1, eqn_idx%B%beg + 1)
2511# 156 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2512 else ! 2D/3D: Bx, By, Bz as variables
2513 b%L(1) = ql_prim_rsx_vf(j, k, l, eqn_idx%B%beg)
2514# 158 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2515 b%R(1) = qr_prim_rsx_vf(j, k, l + 1, eqn_idx%B%beg)
2516# 159 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2517 b%L(2) = ql_prim_rsx_vf(j, k, l, eqn_idx%B%beg + 1)
2518# 160 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2519 b%R(2) = qr_prim_rsx_vf(j, k, l + 1, eqn_idx%B%beg + 1)
2520# 161 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2521 b%L(3) = ql_prim_rsx_vf(j, k, l, eqn_idx%B%beg + 2)
2522# 162 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2523 b%R(3) = qr_prim_rsx_vf(j, k, l + 1, eqn_idx%B%beg + 2)
2524# 163 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2525 end if
2526 end if
2527
2528 rho_l = 0._wp
2529 gamma_l = 0._wp
2530 pi_inf_l = 0._wp
2531 qv_l = 0._wp
2532
2533 rho_r = 0._wp
2534 gamma_r = 0._wp
2535 pi_inf_r = 0._wp
2536 qv_r = 0._wp
2537
2538 alpha_l_sum = 0._wp
2539 alpha_r_sum = 0._wp
2540
2541 pres_mag%L = 0._wp
2542 pres_mag%R = 0._wp
2543
2544 if (mpp_lim) then
2545
2546# 183 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2547#if defined(MFC_OpenACC)
2548# 183 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2549!$acc loop seq
2550# 183 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2551#elif defined(MFC_OpenMP)
2552# 183 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2553
2554# 183 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2555#endif
2556# 184 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2557 do i = 1, num_fluids
2558 alpha_rho_l(i) = max(0._wp, alpha_rho_l(i))
2559 alpha_l(i) = min(max(0._wp, alpha_l(i)), 1._wp)
2560 alpha_l_sum = alpha_l_sum + alpha_l(i)
2561 alpha_rho_r(i) = max(0._wp, alpha_rho_r(i))
2562 alpha_r(i) = min(max(0._wp, alpha_r(i)), 1._wp)
2563 alpha_r_sum = alpha_r_sum + alpha_r(i)
2564 end do
2565
2566 alpha_l = alpha_l/max(alpha_l_sum, sgm_eps)
2567 alpha_r = alpha_r/max(alpha_r_sum, sgm_eps)
2568 end if
2569
2570 call s_compute_mixture_coefficients(alpha_rho_l, alpha_l, rho_l, gamma_l, pi_inf_l, qv_l)
2571 call s_compute_mixture_coefficients(alpha_rho_r, alpha_r, rho_r, gamma_r, pi_inf_r, qv_r)
2572
2573 if (viscous) then
2574 call s_compute_interface_reynolds(alpha_l, re_l, re_size_loc1, re_size_loc2)
2575 call s_compute_interface_reynolds(alpha_r, re_r, re_size_loc1, re_size_loc2)
2576 end if
2577
2578 if (chemistry) then
2579 ! Only the Roe path writes this; zero it so the arithmetic path reads nothing undefined.
2580 c_sum_yi_phi = 0._wp
2581
2582# 208 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2583#if defined(MFC_OpenACC)
2584# 208 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2585!$acc loop seq
2586# 208 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2587#elif defined(MFC_OpenMP)
2588# 208 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2589
2590# 208 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2591#endif
2592# 209 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2593 do i = eqn_idx%species%beg, eqn_idx%species%end
2594 ys_l(i - eqn_idx%species%beg + 1) = ql_prim_rsx_vf(j, k, l, i)
2595# 211 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2596 ys_r(i - eqn_idx%species%beg + 1) = qr_prim_rsx_vf(j, k, l + 1, i)
2597# 212 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2598 end do
2599
2600 call get_mixture_molecular_weight(ys_l, mw_l)
2601 call get_mixture_molecular_weight(ys_r, mw_r)
2602 xs_l(1:num_species) = ys_l(1:num_species)*mw_l/molecular_weights(:)
2603 xs_r(1:num_species) = ys_r(1:num_species)*mw_r/molecular_weights(:)
2604
2605 r_gas_l = gas_constant/mw_l
2606 r_gas_r = gas_constant/mw_r
2607 t_l = pres_l/rho_l/r_gas_l
2608 t_r = pres_r/rho_r/r_gas_r
2609
2610 call get_mixture_caloric_state(t_l, ys_l, cp_il, cp_l, cv_l, e_l)
2611 call get_mixture_caloric_state(t_r, ys_r, cp_ir, cp_r, cv_r, e_r)
2612
2613 if (chem_params%gamma_method == 1) then
2614 ! gamma_method = 1: Ref. Section 2.3.1 Formulation of doi:10.7907/ZKW8-ES97.
2615 gamma_il(1:num_species) = cp_il(1:num_species)/(cp_il(1:num_species) - 1.0_wp)
2616 gamma_ir(1:num_species) = cp_ir(1:num_species)/(cp_ir(1:num_species) - 1.0_wp)
2617
2618 gamma_l = sum(xs_l(1:num_species)/(gamma_il(1:num_species) - 1.0_wp))
2619 gamma_r = sum(xs_r(1:num_species)/(gamma_ir(1:num_species) - 1.0_wp))
2620 else if (chem_params%gamma_method == 2) then
2621 ! gamma_method = 2: c_p / c_v where c_p, c_v are specific heats.
2622 gamm_l = cp_l/cv_l
2623 gamma_l = 1.0_wp/(gamm_l - 1.0_wp)
2624 gamm_r = cp_r/cv_r
2625 gamma_r = 1.0_wp/(gamm_r - 1.0_wp)
2626 end if
2627
2628 e_l = rho_l*e_l + 5.e-1*rho_l*vel_l_rms
2629 e_r = rho_r*e_r + 5.e-1*rho_r*vel_r_rms
2630 h_l = (e_l + pres_l)/rho_l
2631 h_r = (e_r + pres_r)/rho_r
2632 else if (mhd .and. relativity) then
2633 ga%L = 1._wp/sqrt(1._wp - vel_l_rms)
2634 ga%R = 1._wp/sqrt(1._wp - vel_r_rms)
2635# 250 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2636 vdotb%L = vel_l(1)*b%L(1) + vel_l(2)*b%L(2) + vel_l(3)*b%L(3)
2637 vdotb%R = vel_r(1)*b%R(1) + vel_r(2)*b%R(2) + vel_r(3)*b%R(3)
2638
2639 b4%L(1:3) = b%L(1:3)/ga%L + ga%L*vel_l(1:3)*vdotb%L
2640 b4%R(1:3) = b%R(1:3)/ga%R + ga%R*vel_r(1:3)*vdotb%R
2641 b2%L = b%L(1)**2._wp + b%L(2)**2._wp + b%L(3)**2._wp
2642 b2%R = b%R(1)**2._wp + b%R(2)**2._wp + b%R(3)**2._wp
2643# 258 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2644
2645 pres_mag%L = 0.5_wp*(b2%L/ga%L**2._wp + vdotb%L**2._wp)
2646 pres_mag%R = 0.5_wp*(b2%R/ga%R**2._wp + vdotb%R**2._wp)
2647
2648 h_l = f_relativistic_enthalpy(pres_l, rho_l, gamma_l)
2649 h_r = f_relativistic_enthalpy(pres_r, rho_r, gamma_r)
2650# 265 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2651 cm%L(1:3) = (rho_l*h_l*ga%L**2 + b2%L)*vel_l(1:3) - vdotb%L*b%L(1:3)
2652 cm%R(1:3) = (rho_r*h_r*ga%R**2 + b2%R)*vel_r(1:3) - vdotb%R*b%R(1:3)
2653# 268 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2654
2655 e_l = rho_l*h_l*ga%L**2 - pres_l + 0.5_wp*(b2%L + vel_l_rms*b2%L - vdotb%L**2._wp) - rho_l*ga%L
2656 e_r = rho_r*h_r*ga%R**2 - pres_r + 0.5_wp*(b2%R + vel_r_rms*b2%R - vdotb%R**2._wp) - rho_r*ga%R
2657 else if (mhd .and. .not. relativity) then
2658# 273 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2659 pres_mag%L = 0.5_wp*(b%L(1)**2._wp + b%L(2)**2._wp + b%L(3)**2._wp)
2660 pres_mag%R = 0.5_wp*(b%R(1)**2._wp + b%R(2)**2._wp + b%R(3)**2._wp)
2661# 276 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2662 call s_compute_energy(pres_l, alpha_rho_l, alpha_l, vel_l_rms, e_l)
2663 e_l = e_l + pres_mag%L
2664 ! includes magnetic energy
2665 call s_compute_energy(pres_r, alpha_rho_r, alpha_r, vel_r_rms, e_r)
2666 e_r = e_r + pres_mag%R
2667 h_l = (e_l + pres_l - pres_mag%L)/rho_l
2668 ! stagnation enthalpy here excludes magnetic energy (only used to find speed of sound)
2669 h_r = (e_r + pres_r - pres_mag%R)/rho_r
2670 else
2671 call s_compute_energy(pres_l, alpha_rho_l, alpha_l, vel_l_rms, e_l)
2672 call s_compute_energy(pres_r, alpha_rho_r, alpha_r, vel_r_rms, e_r)
2673 h_l = (e_l + pres_l)/rho_l
2674 h_r = (e_r + pres_r)/rho_r
2675 end if
2676
2677 ! elastic energy update
2678 ! HLL handles hypoelasticity as inline branches like this one - one of three code
2679 ! shapes (cf. HLLC inline, and the separate HLLD dual-pass module). See s_riemann_solver.
2680 if (hypoelasticity) then
2681
2682# 295 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2683#if defined(MFC_OpenACC)
2684# 295 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2685!$acc loop seq
2686# 295 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2687#elif defined(MFC_OpenMP)
2688# 295 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2689
2690# 295 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2691#endif
2692# 296 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2693 do i = 1, eqn_idx%stress%end - eqn_idx%stress%beg + 1
2694 tau_e_l(i) = ql_prim_rsx_vf(j, k, l, eqn_idx%stress%beg - 1 + i)
2695# 298 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2696 tau_e_r(i) = qr_prim_rsx_vf(j, k, l + 1, eqn_idx%stress%beg - 1 + i)
2697# 299 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2698 end do
2699
2700 damage_l = 0._wp; damage_r = 0._wp
2701 if (cont_damage) then
2702 damage_l = ql_prim_rsx_vf(j, k, l, eqn_idx%damage)
2703# 304 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2704 damage_r = qr_prim_rsx_vf(j, k, l + 1, eqn_idx%damage)
2705# 305 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2706 end if
2707
2708 call s_compute_hypoelastic_interface_energy(num_fluids, alpha_l, alpha_r, damage_l, damage_r, &
2709 & tau_e_l, tau_e_r, g_l, g_r, e_l, e_r)
2710 end if
2711
2712 ! Only the pressure-based wave-speed estimate reads the averaged state, and the Roe
2713 ! average costs eight square roots per face.
2714 if (wave_speeds == wave_speeds_pressure) then
2715 call s_compute_average_state(rho_l, rho_r, vel_l, vel_r, h_l, h_r, gamma_l, gamma_r, qv_l, qv_r, &
2716 & rho_avg, vel_avg_rms, h_avg, gamma_avg, qv_avg)
2717 if (chemistry .and. avg_state == avg_state_roe) then
2718 r_species(1:num_species) = gas_constant/molecular_weights
2719 call get_species_enthalpies_rt(t_l, h_il)
2720 call get_species_enthalpies_rt(t_r, h_ir)
2721 h_il(1:num_species) = h_il(1:num_species)*r_species(1:num_species)*t_l
2722 h_ir(1:num_species) = h_ir(1:num_species)*r_species(1:num_species)*t_r
2723 call s_compute_chemistry_average_state(rho_l, rho_r, t_l, t_r, ys_l, ys_r, r_species, h_il, &
2724 & h_ir, cp_il, cp_ir, vel_avg_rms, gamma_avg, &
2725 & c_sum_yi_phi)
2726 end if
2727 end if
2728
2729 call s_compute_speed_of_sound(pres_l, rho_l, gamma_l, pi_inf_l, alpha_l, c_l, alpha_rho_l)
2730
2731 call s_compute_speed_of_sound(pres_r, rho_r, gamma_r, pi_inf_r, alpha_r, c_r, alpha_rho_r)
2732
2733 if (wave_speeds == wave_speeds_pressure) then
2734 call s_compute_speed_of_sound_avg(pres_r, rho_avg, gamma_avg, pi_inf_r, qv_avg, vel_avg_rms, &
2735 & h_avg, c_sum_yi_phi, alpha_r, c_avg, alpha_rho_r)
2736 end if
2737
2738 if (mhd) then
2739 call s_compute_fast_magnetosonic_speed(rho_l, c_l, b%L, norm_dir, c_fast%L, h_l)
2740 call s_compute_fast_magnetosonic_speed(rho_r, c_r, b%R, norm_dir, c_fast%R, h_r)
2741 end if
2742
2743 if (viscous) then
2744 if (chemistry) then
2745 call compute_viscosity_and_inversion(t_l, ys_l, t_r, ys_r, re_l(1), re_r(1))
2746 end if
2747
2748# 346 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2749#if defined(MFC_OpenACC)
2750# 346 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2751!$acc loop seq
2752# 346 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2753#elif defined(MFC_OpenMP)
2754# 346 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2755
2756# 346 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2757#endif
2758# 347 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2759 do i = 1, 2
2760 re_avg_rsx_vf(j, k, l, i) = 2._wp/(1._wp/re_l(i) + 1._wp/re_r(i))
2761# 349 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2762 end do
2763 end if
2764
2765 ! Wave speed estimates (wave_speeds=1: direct, wave_speeds=2: pressure-based)
2766 if (wave_speeds == wave_speeds_direct) then
2767 if (mhd) then
2768 ! MHD: use fast magnetosonic speed
2769 s_l = min(vel_l(dir_idx(1)) - c_fast%L, vel_r(dir_idx(1)) - c_fast%R)
2770 s_r = max(vel_r(dir_idx(1)) + c_fast%R, vel_l(dir_idx(1)) + c_fast%L)
2771 else if (hypoelasticity) then
2772 ! Elastic wave speed, Rodriguez et al. JCP (2019)
2773 s_l = min(vel_l(dir_idx(1)) - f_elastic_signal_speed(c_l, g_l, tau_e_l(dir_idx_tau(1)), &
2774 & rho_l), vel_r(dir_idx(1)) - f_elastic_signal_speed(c_r, g_r, &
2775 & tau_e_r(dir_idx_tau(1)), rho_r))
2776 s_r = max(vel_r(dir_idx(1)) + f_elastic_signal_speed(c_r, g_r, tau_e_r(dir_idx_tau(1)), &
2777 & rho_r), vel_l(dir_idx(1)) + f_elastic_signal_speed(c_l, g_l, &
2778 & tau_e_l(dir_idx_tau(1)), rho_l))
2779 else
2780 s_l = min(vel_l(dir_idx(1)) - c_l, vel_r(dir_idx(1)) - c_r)
2781 s_r = max(vel_r(dir_idx(1)) + c_r, vel_l(dir_idx(1)) + c_l)
2782 end if
2783
2784 if (hyper_cleaning) then
2785 ! Dedner GLM divergence cleaning, Dedner et al. JCP (2002)
2786 s_l = min(s_l, -hyper_cleaning_speed)
2787 s_r = max(s_r, hyper_cleaning_speed)
2788 end if
2789
2790 s_s = (pres_r - pres_l + rho_l*vel_l(dir_idx(1))*(s_l - vel_l(dir_idx(1))) &
2791 & - rho_r*vel_r(dir_idx(1))*(s_r - vel_r(dir_idx(1))))/(rho_l*(s_l - vel_l(dir_idx(1))) &
2792 & - rho_r*(s_r - vel_r(dir_idx(1))))
2793 else if (wave_speeds == wave_speeds_pressure) then
2794 pres_sl = 5.e-1_wp*(pres_l + pres_r + rho_avg*c_avg*(vel_l(dir_idx(1)) - vel_r(dir_idx(1))))
2795
2796 pres_sr = pres_sl
2797
2798 ! Low Mach correction: Thornber et al. JCP (2008)
2799 ms_l = max(1._wp, &
2800 & sqrt(1._wp + 5.e-1_wp*(f_isentrope_exponent(gamma_l) + 1._wp) &
2801 & /f_isentrope_exponent(gamma_l)*(pres_sl - pres_l)/(pres_l &
2802 & + f_isentrope_pressure(pi_inf_l, gamma_l))))
2803 ms_r = max(1._wp, &
2804 & sqrt(1._wp + 5.e-1_wp*(f_isentrope_exponent(gamma_r) + 1._wp) &
2805 & /f_isentrope_exponent(gamma_r)*(pres_sr - pres_r)/(pres_r &
2806 & + f_isentrope_pressure(pi_inf_r, gamma_r))))
2807
2808 s_l = vel_l(dir_idx(1)) - c_l*ms_l
2809 s_r = vel_r(dir_idx(1)) + c_r*ms_r
2810
2811 s_s = 5.e-1_wp*((vel_l(dir_idx(1)) + vel_r(dir_idx(1))) + (pres_l - pres_r)/(rho_avg*c_avg))
2812 end if
2813
2814 s_m = min(0._wp, s_l); s_p = max(0._wp, s_r)
2815
2816 xi_m = (5.e-1_wp + sign(5.e-1_wp, s_l)) + (5.e-1_wp - sign(5.e-1_wp, s_l))*(5.e-1_wp + sign(5.e-1_wp, &
2817 & s_r))
2818 xi_p = (5.e-1_wp - sign(5.e-1_wp, s_r)) + (5.e-1_wp - sign(5.e-1_wp, s_l))*(5.e-1_wp + sign(5.e-1_wp, &
2819 & s_r))
2820
2821 ! HLL intercell flux: F* = (s_R*F_L - s_L*F_R + s_L*s_R*(U_R - U_L)) / (s_R - s_L) Low Mach correction
2822 pcorr = f_low_mach_pcorr_hll(vel_l_rms, vel_r_rms, c_l, c_r, rho_l, rho_r, s_m, s_p)
2823
2824 ! Mass
2825 if (.not. relativity) then
2826
2827# 413 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2828#if defined(MFC_OpenACC)
2829# 413 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2830!$acc loop seq
2831# 413 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2832#elif defined(MFC_OpenMP)
2833# 413 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2834
2835# 413 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2836#endif
2837# 414 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2838 do i = 1, eqn_idx%cont%end
2839 flux_rsx_vf(j, k, l, &
2840# 416 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2841 & i) = (s_m*alpha_rho_r(i)*vel_r(norm_dir) - s_p*alpha_rho_l(i)*vel_l(norm_dir) &
2842 & + s_m*s_p*(alpha_rho_l(i) - alpha_rho_r(i)))/(s_m - s_p)
2843 end do
2844 else if (relativity) then
2845
2846# 420 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2847#if defined(MFC_OpenACC)
2848# 420 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2849!$acc loop seq
2850# 420 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2851#elif defined(MFC_OpenMP)
2852# 420 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2853
2854# 420 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2855#endif
2856# 421 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2857 do i = 1, eqn_idx%cont%end
2858 flux_rsx_vf(j, k, l, &
2859# 423 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2860 & i) = (s_m*ga%R*alpha_rho_r(i)*vel_r(norm_dir) - s_p*ga%L*alpha_rho_l(i) &
2861 & *vel_l(norm_dir) + s_m*s_p*(ga%L*alpha_rho_l(i) - ga%R*alpha_rho_r(i)))/(s_m &
2862 & - s_p)
2863 end do
2864 end if
2865
2866 ! Momentum
2867 if (mhd .and. (.not. relativity)) then
2868
2869# 431 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2870#if defined(MFC_OpenACC)
2871# 431 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2872!$acc loop seq
2873# 431 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2874#elif defined(MFC_OpenMP)
2875# 431 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2876
2877# 431 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2878#endif
2879# 432 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2880 do i = 1, 3
2881 ! Flux of rho*v_i in the z direction = rho * v_i * v_z - B_i * B_z +
2882# 434 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2883 ! delta_(z,i) * p_tot
2884# 435 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2885 flux_rsx_vf(j, k, l, &
2886# 436 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2887 & eqn_idx%cont%end + i) = (s_m*(rho_r*vel_r(i)*vel_r(norm_dir) - b%R(i) &
2888 & *b%R(norm_dir) + dir_flg(i)*(pres_r + pres_mag%R)) - s_p*(rho_l*vel_l(i) &
2889 & *vel_l(norm_dir) - b%L(i)*b%L(norm_dir) + dir_flg(i)*(pres_l + pres_mag%L)) &
2890 & + s_m*s_p*(rho_l*vel_l(i) - rho_r*vel_r(i)))/(s_m - s_p)
2891 end do
2892 else if (mhd .and. relativity) then
2893
2894# 442 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2895#if defined(MFC_OpenACC)
2896# 442 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2897!$acc loop seq
2898# 442 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2899#elif defined(MFC_OpenMP)
2900# 442 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2901
2902# 442 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2903#endif
2904# 443 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2905 do i = 1, 3
2906 ! Flux of m_i in the z direction = m_i * v_z - b_i/Gamma * B_z +
2907# 445 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2908 ! delta_(z,i) * p_tot
2909# 446 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2910 flux_rsx_vf(j, k, l, &
2911# 447 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2912 & eqn_idx%cont%end + i) = (s_m*(cm%R(i)*vel_r(norm_dir) - b4%R(i) &
2913 & /ga%R*b%R(norm_dir) + dir_flg(i)*(pres_r + pres_mag%R)) - s_p*(cm%L(i) &
2914 & *vel_l(norm_dir) - b4%L(i)/ga%L*b%L(norm_dir) + dir_flg(i)*(pres_l + pres_mag%L) &
2915 & ) + s_m*s_p*(cm%L(i) - cm%R(i)))/(s_m - s_p)
2916 end do
2917 else if (bubbles_euler) then
2918
2919# 453 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2920#if defined(MFC_OpenACC)
2921# 453 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2922!$acc loop seq
2923# 453 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2924#elif defined(MFC_OpenMP)
2925# 453 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2926
2927# 453 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2928#endif
2929# 454 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2930 do i = 1, num_vels
2931 flux_rsx_vf(j, k, l, &
2932# 456 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2933 & eqn_idx%cont%end + dir_idx(i)) = (s_m*(rho_r*vel_r(dir_idx(1))*vel_r(dir_idx(i)) &
2934 & + dir_flg(dir_idx(i))*(pres_r - ptilde_r)) - s_p*(rho_l*vel_l(dir_idx(1)) &
2935 & *vel_l(dir_idx(i)) + dir_flg(dir_idx(i))*(pres_l - ptilde_l)) &
2936 & + s_m*s_p*(rho_l*vel_l(dir_idx(i)) - rho_r*vel_r(dir_idx(i))))/(s_m - s_p) &
2937 & + (s_m/s_l)*(s_p/s_r)*pcorr*(vel_r(dir_idx(i)) - vel_l(dir_idx(i)))
2938 end do
2939 else if (hypoelasticity) then
2940
2941# 463 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2942#if defined(MFC_OpenACC)
2943# 463 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2944!$acc loop seq
2945# 463 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2946#elif defined(MFC_OpenMP)
2947# 463 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2948
2949# 463 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2950#endif
2951# 464 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2952 do i = 1, num_vels
2953 flux_rsx_vf(j, k, l, &
2954# 466 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2955 & eqn_idx%cont%end + dir_idx(i)) = (s_m*(rho_r*vel_r(dir_idx(1))*vel_r(dir_idx(i)) &
2956 & + dir_flg(dir_idx(i))*pres_r - tau_e_r(dir_idx_tau(i))) &
2957 & - s_p*(rho_l*vel_l(dir_idx(1))*vel_l(dir_idx(i)) + dir_flg(dir_idx(i))*pres_l &
2958 & - tau_e_l(dir_idx_tau(i))) + s_m*s_p*(rho_l*vel_l(dir_idx(i)) &
2959 & - rho_r*vel_r(dir_idx(i))))/(s_m - s_p)
2960 end do
2961 else
2962
2963# 473 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2964#if defined(MFC_OpenACC)
2965# 473 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2966!$acc loop seq
2967# 473 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2968#elif defined(MFC_OpenMP)
2969# 473 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2970
2971# 473 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2972#endif
2973# 474 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2974 do i = 1, num_vels
2975 flux_rsx_vf(j, k, l, &
2976# 476 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2977 & eqn_idx%cont%end + dir_idx(i)) = (s_m*(rho_r*vel_r(dir_idx(1))*vel_r(dir_idx(i)) &
2978 & + dir_flg(dir_idx(i))*pres_r) - s_p*(rho_l*vel_l(dir_idx(1))*vel_l(dir_idx(i)) &
2979 & + dir_flg(dir_idx(i))*pres_l) + s_m*s_p*(rho_l*vel_l(dir_idx(i)) &
2980 & - rho_r*vel_r(dir_idx(i))))/(s_m - s_p) + (s_m/s_l)*(s_p/s_r) &
2981 & *pcorr*(vel_r(dir_idx(i)) - vel_l(dir_idx(i)))
2982 end do
2983 end if
2984
2985 ! Energy
2986 if (mhd .and. (.not. relativity)) then
2987 ! energy flux = (E + p + p_mag) * v_z - B_z * (v_x*B_x + v_y*B_y + v_z*B_z)
2988# 487 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2989# 488 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2990 flux_rsx_vf(j, k, l, &
2991# 489 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2992 & eqn_idx%E) = (s_m*(vel_r(norm_dir)*(e_r + pres_r + pres_mag%R) - b%R(norm_dir) &
2993 & *(vel_r(1)*b%R(1) + vel_r(2)*b%R(2) + vel_r(3)*b%R(3))) - s_p*(vel_l(norm_dir) &
2994 & *(e_l + pres_l + pres_mag%L) - b%L(norm_dir)*(vel_l(1)*b%L(1) + vel_l(2)*b%L(2) &
2995 & + vel_l(3)*b%L(3))) + s_m*s_p*(e_l - e_r))/(s_m - s_p)
2996# 494 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
2997 else if (mhd .and. relativity) then
2998 ! energy flux = m_z - mass flux Hard-coded for single-component for now
2999# 496 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3000 flux_rsx_vf(j, k, l, &
3001# 497 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3002 & eqn_idx%E) = (s_m*(cm%R(norm_dir) - ga%R*alpha_rho_r(1)*vel_r(norm_dir)) &
3003 & - s_p*(cm%L(norm_dir) - ga%L*alpha_rho_l(1)*vel_l(norm_dir)) + s_m*s_p*(e_l - e_r)) &
3004 & /(s_m - s_p)
3005 else if (bubbles_euler) then
3006 flux_rsx_vf(j, k, l, &
3007# 502 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3008 & eqn_idx%E) = (s_m*vel_r(dir_idx(1))*(e_r + pres_r - ptilde_r) - s_p*vel_l(dir_idx(1) &
3009 & )*(e_l + pres_l - ptilde_l) + s_m*s_p*(e_l - e_r))/(s_m - s_p) + (s_m/s_l)*(s_p/s_r) &
3010 & *pcorr*(vel_r_rms - vel_l_rms)/2._wp
3011 else if (hypoelasticity) then
3012 flux_tau_l = 0._wp; flux_tau_r = 0._wp
3013
3014# 507 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3015#if defined(MFC_OpenACC)
3016# 507 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3017!$acc loop seq
3018# 507 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3019#elif defined(MFC_OpenMP)
3020# 507 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3021
3022# 507 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3023#endif
3024# 508 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3025 do i = 1, num_dims
3026 flux_tau_l = flux_tau_l + tau_e_l(dir_idx_tau(i))*vel_l(dir_idx(i))
3027 flux_tau_r = flux_tau_r + tau_e_r(dir_idx_tau(i))*vel_r(dir_idx(i))
3028 end do
3029 flux_rsx_vf(j, k, l, &
3030# 513 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3031 & eqn_idx%E) = (s_m*(vel_r(dir_idx(1))*(e_r + pres_r) - flux_tau_r) &
3032 & - s_p*(vel_l(dir_idx(1))*(e_l + pres_l) - flux_tau_l) + s_m*s_p*(e_l - e_r))/(s_m &
3033 & - s_p)
3034 else
3035 flux_rsx_vf(j, k, l, &
3036# 518 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3037 & eqn_idx%E) = (s_m*vel_r(dir_idx(1))*(e_r + pres_r) - s_p*vel_l(dir_idx(1))*(e_l &
3038 & + pres_l) + s_m*s_p*(e_l - e_r))/(s_m - s_p) + (s_m/s_l)*(s_p/s_r)*pcorr*(vel_r_rms &
3039 & - vel_l_rms)/2._wp
3040 end if
3041
3042 ! Elastic Stresses
3043 if (hypoelasticity) then
3044 do i = 1, eqn_idx%stress%end - eqn_idx%stress%beg + 1 ! TODO: this indexing may be slow
3045 flux_rsx_vf(j, k, l, &
3046# 527 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3047 & eqn_idx%stress%beg - 1 + i) = (s_m*(rho_r*vel_r(dir_idx(1))*tau_e_r(i)) &
3048 & - s_p*(rho_l*vel_l(dir_idx(1))*tau_e_l(i)) + s_m*s_p*(rho_l*tau_e_l(i) &
3049 & - rho_r*tau_e_r(i)))/(s_m - s_p)
3050 end do
3051 end if
3052
3053 ! Damage flux: U_D = m_s*D (damageable-solid partial mass)
3054 if (cont_damage) then
3055 solid_partial_density_l = 0._wp; solid_partial_density_r = 0._wp
3056
3057# 536 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3058#if defined(MFC_OpenACC)
3059# 536 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3060!$acc loop seq
3061# 536 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3062#elif defined(MFC_OpenMP)
3063# 536 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3064
3065# 536 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3066#endif
3067# 537 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3068 do i = 1, num_fluids
3069 if (gs_rs(i) > verysmall) then
3070 solid_partial_density_l = solid_partial_density_l + alpha_rho_l(i)
3071 solid_partial_density_r = solid_partial_density_r + alpha_rho_r(i)
3072 end if
3073 end do
3074 flux_rsx_vf(j, k, l, &
3075# 544 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3076 & eqn_idx%damage) = (s_m*(solid_partial_density_r*vel_r(dir_idx(1))*damage_r) &
3077 & - s_p*(solid_partial_density_l*vel_l(dir_idx(1))*damage_l) &
3078 & + s_m*s_p*(solid_partial_density_l*damage_l - solid_partial_density_r*damage_r)) &
3079 & /(s_m - s_p)
3080 end if
3081
3082 ! Export interface velocity for NC RHS
3083 if (hypo_nc_mode == hypo_nc_mode_interface .or. (alt_soundspeed .and. .not. hll_u_interface)) then
3084
3085# 552 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3086#if defined(MFC_OpenACC)
3087# 552 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3088!$acc loop seq
3089# 552 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3090#elif defined(MFC_OpenMP)
3091# 552 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3092
3093# 552 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3094#endif
3095# 553 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3096 do i = 1, num_dims
3097 if (0._wp <= s_l) then
3098 nc_iface_vel_rsx_vf(j, k, l, dir_idx(i)) = vel_l(dir_idx(i))
3099# 556 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3100 else if (s_r <= 0._wp) then
3101 nc_iface_vel_rsx_vf(j, k, l, dir_idx(i)) = vel_r(dir_idx(i))
3102# 558 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3103 else
3104 nc_iface_vel_rsx_vf(j, k, l, &
3105# 560 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3106 & dir_idx(i)) = (s_r*vel_l(dir_idx(i)) - s_l*vel_r(dir_idx(i)))/(s_r &
3107 & - s_l)
3108 end if
3109 end do
3110 end if
3111
3112 if (.not. hll_u_interface) then ! HLL Method 1: per-fluid alpha interface flux
3113 ! Branchless: s_M/s_P fold the upwinding, so signs of s_L/s_R need no
3114 ! per-thread branches (which diverge across GPU warps)
3115
3116# 569 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3117#if defined(MFC_OpenACC)
3118# 569 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3119!$acc loop seq
3120# 569 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3121#elif defined(MFC_OpenMP)
3122# 569 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3123
3124# 569 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3125#endif
3126# 570 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3127 do i = eqn_idx%adv%beg, eqn_idx%adv%end
3128 flux_rsx_vf(j, k, l, i) = (ql_prim_rsx_vf(j, k, l, i) - qr_prim_rsx_vf(j, k, l + 1, &
3129# 572 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3130 & i))*s_m*s_p/(s_m - s_p)
3131 flux_src_rsx_vf(j, k, l, i) = (s_m*qr_prim_rsx_vf(j, k, l + 1, &
3132# 574 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3133 & i) - s_p*ql_prim_rsx_vf(j, k, l, i))/(s_m - s_p)
3134# 575 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3135 end do
3136 else ! HLL Method 2: shared velocity interface flux
3137
3138# 577 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3139#if defined(MFC_OpenACC)
3140# 577 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3141!$acc loop seq
3142# 577 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3143#elif defined(MFC_OpenMP)
3144# 577 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3145
3146# 577 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3147#endif
3148# 578 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3149 do i = eqn_idx%adv%beg, eqn_idx%adv%end
3150 if (0._wp <= s_l) then
3151 flux_rsx_vf(j, k, l, i) = ql_prim_rsx_vf(j, k, l, i)*vel_l(dir_idx(1))
3152# 581 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3153 else if (s_r <= 0._wp) then
3154 flux_rsx_vf(j, k, l, i) = qr_prim_rsx_vf(j, k, l + 1, i)*vel_r(dir_idx(1))
3155# 583 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3156 else
3157 flux_rsx_vf(j, k, l, i) = (s_r*ql_prim_rsx_vf(j, k, l, &
3158# 585 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3159 & i)*vel_l(dir_idx(1)) - s_l*qr_prim_rsx_vf(j, k, l + 1, &
3160# 586 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3161 & i)*vel_r(dir_idx(1)) + s_l*s_r*(qr_prim_rsx_vf(j, k, l + 1, &
3162# 587 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3163 & i) - ql_prim_rsx_vf(j, k, l, i)))/(s_r - s_l)
3164# 588 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3165 end if
3166 end do
3167 if (0._wp <= s_l) then
3168 flux_src_rsx_vf(j, k, l, eqn_idx%adv%beg) = vel_l(dir_idx(1))
3169# 592 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3170 else if (s_r <= 0._wp) then
3171 flux_src_rsx_vf(j, k, l, eqn_idx%adv%beg) = vel_r(dir_idx(1))
3172# 594 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3173 else
3174 flux_src_rsx_vf(j, k, l, &
3175# 596 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3176 & eqn_idx%adv%beg) = (s_r*vel_l(dir_idx(1)) - s_l*vel_r(dir_idx(1)))/(s_r - s_l)
3177 end if
3178 end if
3179
3180 if (bubbles_euler) then
3181 ! From HLLC: Kills mass transport @ bubble gas density
3182 if (num_fluids > 1) then
3183 flux_rsx_vf(j, k, l, eqn_idx%cont%end) = 0._wp
3184# 604 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3185 end if
3186 end if
3187
3188 if (chemistry) then
3189
3190# 608 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3191#if defined(MFC_OpenACC)
3192# 608 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3193!$acc loop seq
3194# 608 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3195#elif defined(MFC_OpenMP)
3196# 608 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3197
3198# 608 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3199#endif
3200# 609 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3201 do i = eqn_idx%species%beg, eqn_idx%species%end
3202 y_l = ql_prim_rsx_vf(j, k, l, i)
3203# 611 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3204 y_r = qr_prim_rsx_vf(j, k, l + 1, i)
3205# 612 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3206
3207 flux_rsx_vf(j, k, l, &
3208# 614 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3209 & i) = (s_m*y_r*rho_r*vel_r(dir_idx(1)) - s_p*y_l*rho_l*vel_l(dir_idx(1)) &
3210 & + s_m*s_p*(y_l*rho_l - y_r*rho_r))/(s_m - s_p)
3211 flux_src_rsx_vf(j, k, l, i) = 0._wp
3212# 617 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3213 end do
3214 end if
3215
3216 ! MHD: magnetic flux and Maxwell stress contributions
3217 if (mhd) then
3218 if (n == 0) then ! 1D: d/dx flux only & Bx = Bx0 = const.
3219 ! B_y flux = v_x * B_y - v_y * Bx0 B_z flux = v_x * B_z - v_z * Bx0
3220
3221# 624 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3222#if defined(MFC_OpenACC)
3223# 624 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3224!$acc loop seq
3225# 624 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3226#elif defined(MFC_OpenMP)
3227# 624 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3228
3229# 624 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3230#endif
3231# 625 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3232 do i = 0, 1
3233 flux_rsx_vf(j, k, l, &
3234 & eqn_idx%B%beg + i) = (s_m*(vel_r(1)*b%R(2 + i) - vel_r(2 + i)*bx0) &
3235 & - s_p*(vel_l(1)*b%L(2 + i) - vel_l(2 + i)*bx0) + s_m*s_p*(b%L(2 + i) &
3236 & - b%R(2 + i)))/(s_m - s_p)
3237 end do
3238 else ! 2D/3D: Bx, By, Bz /= const. but zero flux component in the same direction
3239 ! B_x d/dz flux = (1 - delta(x,z)) * (v_z * B_x - v_x * B_z) B_y
3240# 633 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3241 ! d/dz flux = (1 - delta(y,z)) * (v_z * B_y - v_y * B_z) B_z d/dz
3242# 634 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3243 ! flux = (1 - delta(z,z)) * (v_z * B_z - v_z * B_z)
3244# 635 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3245
3246# 635 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3247#if defined(MFC_OpenACC)
3248# 635 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3249!$acc loop seq
3250# 635 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3251#elif defined(MFC_OpenMP)
3252# 635 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3253
3254# 635 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3255#endif
3256# 636 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3257 do i = 0, 2
3258 flux_rsx_vf(j, k, l, &
3259# 638 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3260 & eqn_idx%B%beg + i) = (s_m*(vel_r(dir_idx(1))*b%R(i + 1) - vel_r(i + 1) &
3261 & *b%R(norm_dir)) - s_p*(vel_l(dir_idx(1))*b%L(i + 1) - vel_l(i + 1) &
3262 & *b%L(norm_dir)) + s_m*s_p*(b%L(i + 1) - b%R(i + 1)))/(s_m - s_p)
3263 end do
3264
3265 if (hyper_cleaning) then
3266 ! propagate magnetic field divergence as a wave
3267 flux_rsx_vf(j, k, l, eqn_idx%B%beg + norm_dir - 1) = flux_rsx_vf(j, k, l, &
3268# 646 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3269 & eqn_idx%B%beg + norm_dir - 1) + (s_m*qr_prim_rsx_vf(j, k, l + 1, &
3270# 647 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3271 & eqn_idx%psi) - s_p*ql_prim_rsx_vf(j, k, l, eqn_idx%psi))/(s_m - s_p)
3272# 648 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3273
3274 flux_rsx_vf(j, k, l, &
3275# 650 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3276 & eqn_idx%psi) = (hyper_cleaning_speed**2*(s_m*b%R(norm_dir) &
3277 & - s_p*b%L(norm_dir)) + s_m*s_p*(ql_prim_rsx_vf(j, k, l, &
3278# 652 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3279 & eqn_idx%psi) - qr_prim_rsx_vf(j, k, l + 1, eqn_idx%psi)))/(s_m - s_p)
3280# 653 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3281 else
3282 ! Without hyperbolic cleaning, make sure flux of B_normal is identically zero
3283 flux_rsx_vf(j, k, l, eqn_idx%B%beg + norm_dir - 1) = 0._wp
3284# 656 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3285 end if
3286 end if
3287 flux_src_rsx_vf(j, k, l, eqn_idx%adv%beg) = 0._wp
3288# 659 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3289 end if
3290
3291# 698 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3292 end do
3293 end do
3294 end do
3295
3296# 701 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3297#if defined(MFC_OpenACC)
3298# 701 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3299!$acc end parallel loop
3300# 701 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3301#elif defined(MFC_OpenMP)
3302# 701 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3303
3304# 701 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3305!$omp end target teams loop
3306# 701 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3307#endif
3308# 702 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3309 end if
3310# 704 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hll.fpp"
3311
3312 if (viscous) then
3313 if (weno_re_flux) then
3314 call s_compute_viscous_source_flux(ql_prim_vf(eqn_idx%mom%beg:eqn_idx%mom%end), &
3315 & dql_prim_dx_vf(eqn_idx%mom%beg:eqn_idx%mom%end), &
3316 & dql_prim_dy_vf(eqn_idx%mom%beg:eqn_idx%mom%end), &
3317 & dql_prim_dz_vf(eqn_idx%mom%beg:eqn_idx%mom%end), &
3318 & qr_prim_vf(eqn_idx%mom%beg:eqn_idx%mom%end), &
3319 & dqr_prim_dx_vf(eqn_idx%mom%beg:eqn_idx%mom%end), &
3320 & dqr_prim_dy_vf(eqn_idx%mom%beg:eqn_idx%mom%end), &
3321 & dqr_prim_dz_vf(eqn_idx%mom%beg:eqn_idx%mom%end), flux_src_vf, q_prim_vf, &
3322 & norm_dir, ix, iy, iz)
3323 else
3324 call s_compute_viscous_source_flux(q_prim_vf(eqn_idx%mom%beg:eqn_idx%mom%end), &
3325 & dql_prim_dx_vf(eqn_idx%mom%beg:eqn_idx%mom%end), &
3326 & dql_prim_dy_vf(eqn_idx%mom%beg:eqn_idx%mom%end), &
3327 & dql_prim_dz_vf(eqn_idx%mom%beg:eqn_idx%mom%end), &
3328 & q_prim_vf(eqn_idx%mom%beg:eqn_idx%mom%end), &
3329 & dqr_prim_dx_vf(eqn_idx%mom%beg:eqn_idx%mom%end), &
3330 & dqr_prim_dy_vf(eqn_idx%mom%beg:eqn_idx%mom%end), &
3331 & dqr_prim_dz_vf(eqn_idx%mom%beg:eqn_idx%mom%end), flux_src_vf, q_prim_vf, &
3332 & norm_dir, ix, iy, iz)
3333 end if
3334 end if
3335
3336 call s_finalize_riemann_solver(flux_vf, flux_src_vf, flux_gsrc_vf, norm_dir)
3337
3338 end subroutine s_hll_riemann_solver
3339
3340end module m_riemann_solver_hll
m_chemistry
Multi-species chemistry interface for thermodynamic properties, reaction rates, and transport coeffic...
Definition m_chemistry.fpp.f90:368
m_chemistry::compute_viscosity_and_inversion
subroutine compute_viscosity_and_inversion(t_l, ys_l, t_r, ys_r, re_l, re_r)
Compute mixture viscosities for left and right states and invert them for use as reciprocal Reynolds ...
Definition m_chemistry.fpp.f90:411
m_constants
Compile-time constant parameters: default values, tolerances, and physical constants.
Definition m_constants.fpp.f90:7
m_constants::avg_state_roe
integer, parameter avg_state_roe
Definition m_constants.fpp.f90:138
m_constants::wave_speeds_direct
integer, parameter wave_speeds_direct
Definition m_constants.fpp.f90:176
m_constants::riemann_solver_hll
integer, parameter riemann_solver_hll
Definition m_constants.fpp.f90:169
m_constants::sgm_eps
real(wp), parameter sgm_eps
Segmentation tolerance.
Definition m_constants.fpp.f90:13
m_constants::riemann_solver_hllc
integer, parameter riemann_solver_hllc
Definition m_constants.fpp.f90:170
m_constants::wave_speeds_pressure
integer, parameter wave_speeds_pressure
Definition m_constants.fpp.f90:177
m_constants::riemann_solver_lax_friedrichs
integer, parameter riemann_solver_lax_friedrichs
Definition m_constants.fpp.f90:172
m_constants::verysmall
real(wp), parameter verysmall
Very small number.
Definition m_constants.fpp.f90:17
m_constants::avg_state_arithmetic
integer, parameter avg_state_arithmetic
Definition m_constants.fpp.f90:139
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::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_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_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::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
Global parameters for the computational domain, fluid properties, and simulation algorithm configurat...
Definition m_global_parameters.fpp.f90:367
m_global_parameters::hypo_nc_mode_interface
integer, parameter hypo_nc_mode_interface
interface-velocity export for the velocity-gradient tensor (HLL Method 2, HLLC)
Definition m_global_parameters.fpp.f90:474
m_global_parameters::hypo_nc_mode
integer hypo_nc_mode
Definition m_global_parameters.fpp.f90:476
m_global_parameters::dir_idx
integer, dimension(3) dir_idx
Definition m_global_parameters.fpp.f90:834
m_global_parameters::dir_idx_tau
integer, dimension(3) dir_idx_tau
(nn, nt, nt2) stress indices for wave speeds and momentum flux
Definition m_global_parameters.fpp.f90:836
m_global_parameters::dir_flg
real(wp), dimension(3) dir_flg
Definition m_global_parameters.fpp.f90:835
m_riemann_solver_hll
HLL approximate Riemann solver, Harten et al. SIAM Review (1983).
Definition m_riemann_solver_hll.fpp.f90:367
m_riemann_solver_hll::s_hll_riemann_solver
subroutine s_hll_riemann_solver(ql_prim_rsx_vf, dql_prim_dx_vf, dql_prim_dy_vf, dql_prim_dz_vf, ql_prim_vf, qr_prim_rsx_vf, dqr_prim_dx_vf, dqr_prim_dy_vf, dqr_prim_dz_vf, qr_prim_vf, q_prim_vf, flux_vf, flux_src_vf, flux_gsrc_vf, norm_dir, ix, iy, iz)
HLL approximate Riemann solver, Harten et al. SIAM Review (1983).
Definition m_riemann_solver_hll.fpp.f90:388
m_riemann_state
Shared Riemann-solver module state and the per-sweep setup, state-buffer population,...
Definition m_riemann_state.fpp.f90:368
m_riemann_state::flux_src_rsx_vf
real(wp), dimension(:,:,:,:), allocatable flux_src_rsx_vf
Definition m_riemann_state.fpp.f90:381
m_riemann_state::s_compute_average_state
subroutine s_compute_average_state(rho_l, rho_r, vel_l, vel_r, h_l, h_r, gamma_l, gamma_r, qv_l, qv_r, rho_avg, vel_avg_rms, h_avg, gamma_avg, qv_avg)
Interface-averaged state that the pressure-based wave-speed estimate reads. avg_state selects between...
Definition m_riemann_state.fpp.f90:836
m_riemann_state::f_elastic_signal_speed
real(wp) function f_elastic_signal_speed(c, g, tau, rho)
Elastic signal speed of Rodriguez et al. JCP (2019): the acoustic speed stiffened by the shear modulu...
Definition m_riemann_state.fpp.f90:569
m_riemann_state::is3
type(int_bounds_info) is3
Definition m_riemann_state.fpp.f90:519
m_riemann_state::flux_rsx_vf
real(wp), dimension(:,:,:,:), allocatable flux_rsx_vf
The cell-boundary values of the fluxes (src - source) that are computed through the chosen Riemann pr...
Definition m_riemann_state.fpp.f90:381
m_riemann_state::s_compute_chemistry_average_state
subroutine s_compute_chemistry_average_state(rho_l, rho_r, t_l, t_r, ys_l, ys_r, r_species, h_il, h_ir, cp_il, cp_ir, vel_avg_rms, gamma_avg, c_sum_yi_phi)
Roe-averaged reacting-mixture quantities: replaces gamma_avg with the mixture Cp/Cv and builds the c_...
Definition m_riemann_state.fpp.f90:940
m_riemann_state::s_initialize_riemann_solver
subroutine s_initialize_riemann_solver(flux_src_vf, norm_dir)
Set up the chosen Riemann solver algorithm for the current direction.
Definition m_riemann_state.fpp.f90:2178
m_riemann_state::s_compute_viscous_source_flux
subroutine s_compute_viscous_source_flux(vell_vf, dvell_dx_vf, dvell_dy_vf, dvell_dz_vf, velr_vf, dvelr_dx_vf, dvelr_dy_vf, dvelr_dz_vf, flux_src_vf, q_prim_vf, norm_dir, ix, iy, iz)
Dispatch to the subroutines that are utilized to compute the viscous source fluxes for either Cartesi...
Definition m_riemann_state.fpp.f90:1023
m_riemann_state::nc_iface_vel_rsx_vf
real(wp), dimension(:,:,:,:), allocatable nc_iface_vel_rsx_vf
Definition m_riemann_state.fpp.f90:413
m_riemann_state::gs_rs
real(wp), dimension(:), allocatable gs_rs
Definition m_riemann_state.fpp.f90:536
m_riemann_state::s_populate_riemann_states_variables_buffers
subroutine s_populate_riemann_states_variables_buffers(ql_prim_rsx_vf, dql_prim_dx_vf, dql_prim_dy_vf, dql_prim_dz_vf, qr_prim_rsx_vf, dqr_prim_dx_vf, dqr_prim_dy_vf, dqr_prim_dz_vf, norm_dir, ix, iy, iz)
Populate the left and right Riemann state variable buffers based on boundary conditions.
Definition m_riemann_state.fpp.f90:1045
m_riemann_state::f_low_mach_pcorr_hll
real(wp) function f_low_mach_pcorr_hll(vel_l_rms, vel_r_rms, c_l, c_r, rho_l, rho_r, s_m, s_p)
Low-Mach pressure correction added to the HLL and Lax-Friedrichs fluxes, which restores the pressure ...
Definition m_riemann_state.fpp.f90:666
m_riemann_state::s_compute_hypoelastic_interface_energy
subroutine s_compute_hypoelastic_interface_energy(nf, alpha_l, alpha_r, damage_l, damage_r, tau_e_l, tau_e_r, g_l, g_r, e_l, e_r)
Accumulate the hypoelastic stress contribution to the energies of the left and right Riemann states: ...
Definition m_riemann_state.fpp.f90:3263
m_riemann_state::is2
type(int_bounds_info) is2
Definition m_riemann_state.fpp.f90:519
m_riemann_state::s_finalize_riemann_solver
subroutine s_finalize_riemann_solver(flux_vf, flux_src_vf, flux_gsrc_vf, norm_dir)
Reshape and copy the Riemann-solver flux buffers back to the physical-space output arrays for the sel...
Definition m_riemann_state.fpp.f90:3432
m_riemann_state::s_compute_interface_reynolds
subroutine s_compute_interface_reynolds(alpha_k, re_k, re_size_loc1, re_size_loc2)
Compute the shear and volume Reynolds numbers of one Riemann state by inverse-weighting the fluid Rey...
Definition m_riemann_state.fpp.f90:3170
m_riemann_state::flux_gsrc_rsx_vf
real(wp), dimension(:,:,:,:), allocatable flux_gsrc_rsx_vf
The cell-boundary values of the geometrical source flux that are computed through the chosen Riemann ...
Definition m_riemann_state.fpp.f90:399
m_riemann_state::re_avg_rsx_vf
real(wp), dimension(:,:,:,:), allocatable re_avg_rsx_vf
Definition m_riemann_state.fpp.f90:503
m_riemann_state::is1
type(int_bounds_info) is1
Definition m_riemann_state.fpp.f90:519
m_variables_conversion
Conservative-to-primitive variable conversion, mixture property evaluation, and pressure computation.
Definition m_variables_conversion.fpp.f90:367
m_variables_conversion::s_compute_energy
subroutine, public s_compute_energy(pres, alpha_rho_k, alpha_k, vel_sum, e)
Total energy per unit volume, thermodynamic terms only. Callers add magnetic and elastic energy,...
Definition m_variables_conversion.fpp.f90:2424
m_variables_conversion::s_compute_fast_magnetosonic_speed
subroutine, public s_compute_fast_magnetosonic_speed(rho, c, b, norm, c_fast, h)
Compute the fast magnetosonic wave speed from the sound speed, density, and magnetic field components...
Definition m_variables_conversion.fpp.f90:2550
m_derived_types::int_bounds_info
Integer bounds for variables.
Definition m_derived_types.fpp.f90:450
m_derived_types::riemann_states_vec3
Left and right Riemann states for 3-component vectors.
Definition m_derived_types.fpp.f90:438
m_derived_types::riemann_states
Left and right Riemann states.
Definition m_derived_types.fpp.f90:432
m_derived_types::scalar_field
Derived type annexing a scalar field (SF).
Definition m_derived_types.fpp.f90:370