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