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