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