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