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