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