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