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