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