MFC
Exascale flow solver
Loading...
Searching...
No Matches
m_riemann_solver_hlld.fpp.f90
Go to the documentation of this file.
1# 1 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
2!>
3!! @file
4!! @brief Contains module m_riemann_solver_hlld
5
6!> @brief HLLD approximate Riemann solver for MHD, Miyoshi & Kusano JCP (2005)
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_hlld.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# 15 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
32# 16 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
33# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
34
35# 24 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
36
37# 53 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
38
39# 65 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
40
41# 75 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
42
43# 105 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
44
45# 117 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
46
47# 127 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
48
49# 174 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
50! New line at end of file is required for FYPP
51# 2 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
52# 1 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 1
53# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
54# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
55# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
56# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
57# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
58# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
59
60# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
61# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
62# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
63
64# 15 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
65# 16 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
66# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
67
68# 24 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
69
70# 53 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
71
72# 65 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
73
74# 75 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
75
76# 105 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
77
78# 117 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
79
80# 127 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
81
82# 174 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
83! New line at end of file is required for FYPP
84# 2 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp" 2
85
86# 4 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
87# 5 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
88# 6 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
89# 7 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
90# 8 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
91
92# 20 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
93
94# 43 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
95
96# 48 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
97
98# 53 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
99
100# 58 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
101
102# 63 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
103
104# 68 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
105
106# 81 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
107
108# 86 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
109
110# 91 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
111
112# 96 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
113
114# 101 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
115
116# 106 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
117
118# 111 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
119
120# 116 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
121
122# 121 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
123
124# 126 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
125
126# 156 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
127
128# 197 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
129
130# 211 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
131
132# 236 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
133
134# 247 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
135
136# 249 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
137# 260 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
138
139# 310 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
140
141# 320 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
142
143# 330 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
144
145# 339 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
146
147# 356 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
148
149# 366 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
150
151# 373 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
152
153# 379 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
154
155# 385 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
156
157# 391 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
158
159# 397 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
160
161# 403 "/home/runner/work/MFC/MFC/src/common/include/omp_macros.fpp"
162! New line at end of file is required for FYPP
163# 3 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
164# 1 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 1
165# 1 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp" 1
166# 2 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
167# 3 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
168# 4 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
169# 5 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
170# 6 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
171
172# 8 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
173# 9 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
174# 10 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
175
176# 15 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
177# 16 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
178# 17 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
179
180# 24 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
181
182# 53 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
183
184# 65 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
185
186# 75 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
187
188# 105 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
189
190# 117 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
191
192# 127 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
193
194# 174 "/home/runner/work/MFC/MFC/src/common/include/shared_parallel_macros.fpp"
195! New line at end of file is required for FYPP
196# 2 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp" 2
197
198# 7 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
199
200# 17 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
201
202# 22 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
203
204# 27 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
205
206# 32 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
207
208# 37 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
209
210# 42 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
211
212# 47 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
213
214# 52 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
215
216# 57 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
217
218# 62 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
219
220# 73 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
221
222# 78 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
223
224# 83 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
225
226# 88 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
227
228# 103 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
229
230# 131 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
231
232# 160 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
233
234# 175 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
235
236# 193 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
237
238# 215 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
239
240# 244 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
241
242# 259 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
243
244# 269 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
245
246# 278 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
247
248# 294 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
249
250# 304 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
251
252# 311 "/home/runner/work/MFC/MFC/src/common/include/acc_macros.fpp"
253! New line at end of file is required for FYPP
254# 4 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp" 2
255
256! GPU parallel region (scalar reductions, maxval/minval)
257# 23 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
258
259! GPU parallel loop over threads (most common GPU macro)
260# 43 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
261
262! Required closing for GPU_PARALLEL_LOOP
263# 55 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
264
265! Mark routine for device compilation
266# 112 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
267
268! Declare device-resident data
269# 130 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
270
271! Inner loop within a GPU parallel region
272# 145 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
273
274! Scoped GPU data region
275# 164 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
276
277! Host code with device pointers (for MPI with GPU buffers)
278# 193 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
279
280! Allocate device memory (unscoped)
281# 207 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
282
283! Free device memory
284# 219 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
285
286! Atomic operation on device
287# 231 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
288
289! End atomic capture block
290# 242 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
291
292! Copy data between host and device
293# 254 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
294
295! Synchronization barrier
296# 266 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
297
298! Import GPU library module (openacc or omp_lib)
299# 275 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
300
301! Emit code only for AMD compiler
302# 282 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
303
304! Emit code for non-Cray compilers
305# 289 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
306
307! Emit code only for Cray compiler
308# 296 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
309
310! Emit code for non-NVIDIA compilers
311# 303 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
312
313# 305 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
314# 306 "/home/runner/work/MFC/MFC/src/common/include/parallel_macros.fpp"
315! New line at end of file is required for FYPP
316# 2 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp" 2
317
318# 14 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
319
320! Caution: This macro requires the use of a binding script to set CUDA_VISIBLE_DEVICES, such that we have one GPU device per MPI
321! rank. That's because for both cudaMemAdvise (preferred location) and cudaMemPrefetchAsync we use location = device_id = 0. For an
322! example see misc/nvidia_uvm/bind.sh.
323# 52 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
324
325! Allocate and create GPU device memory
326# 72 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
327
328! Free GPU device memory and deallocate
329# 80 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
330
331! Cray-specific GPU pointer setup for vector fields
332# 104 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
333
334! Cray-specific GPU pointer setup for scalar fields
335# 120 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
336
337! Cray-specific GPU pointer setup for acoustic source spatials
338# 145 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
339
340# 151 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
341
342# 158 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
343! New line at end of file is required for FYPP
344# 8 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp" 2
345
346module m_riemann_solver_hlld
347
348 use m_derived_types
349 use m_global_parameters
350 use m_variables_conversion
351 use m_riemann_state
352
353 implicit none
354
355contains
356
357 !> HLLD Riemann solver for MHD, Miyoshi & Kusano JCP (2005)
358 subroutine s_hlld_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, &
360 & flux_src_vf, 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), allocatable, dimension(:), intent(inout) :: dqL_prim_dx_vf, dqR_prim_dx_vf, dqL_prim_dy_vf, &
364 & dqR_prim_dy_vf, dqL_prim_dz_vf, dqR_prim_dz_vf
365
366 type(scalar_field), allocatable, dimension(:), intent(inout) :: qL_prim_vf, qR_prim_vf
367 type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
368 type(scalar_field), dimension(sys_size), intent(inout) :: flux_vf, flux_src_vf, flux_gsrc_vf
369 integer, intent(in) :: norm_dir
370 type(int_bounds_info), intent(in) :: ix, iy, iz
371
372 ! Local variables:
373
374# 40 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
375 real(wp), dimension(num_fluids) :: alpha_L, alpha_R, alpha_rho_L, alpha_rho_R
376# 42 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
377 type(riemann_states_vec3) :: vel
378 type(riemann_states) :: rho, pres, E, H_no_mag
379 type(riemann_states) :: gamma, pi_inf, qv
380 type(riemann_states) :: vel_rms
381 type(riemann_states_vec3) :: B
382 type(riemann_states) :: c, c_fast, pres_mag
383
384 ! HLLD speeds and intermediate state variables:
385 real(wp) :: s_L, s_R, s_M, s_starL, s_starR
386 real(wp) :: pTot_L, pTot_R, p_star, rhoL_star, rhoR_star, E_starL, E_starR
387 real(wp), dimension(7) :: U_L, U_R, U_starL, U_starR, U_doubleL, U_doubleR
388 real(wp), dimension(7) :: F_L, F_R, F_starL, F_starR, F_hlld
389
390 ! Indices for U and F: (rho, rho*vel(1), rho*vel(2), rho*vel(3), By, Bz, E) Note: vel and B are permutated, so vel(1) is the
391 ! normal velocity, and x is the normal direction Note: Bx is omitted as the magnetic flux is always zero in the normal
392 ! direction
393
394 real(wp) :: sqrt_rhoL_star, sqrt_rhoR_star, denom_ds, sign_Bx
395 real(wp) :: vL_star, vR_star, wL_star, wR_star
396 real(wp) :: v_double, w_double, By_double, Bz_double, E_doubleL, E_doubleR, E_double
397 integer :: i, j, k, l
398
399 call s_populate_riemann_states_variables_buffers(ql_prim_rsx_vf, dql_prim_dx_vf, dql_prim_dy_vf, dql_prim_dz_vf, &
400 & qr_prim_rsx_vf, dqr_prim_dx_vf, dqr_prim_dy_vf, dqr_prim_dz_vf, norm_dir, ix, iy, iz)
401
402 call s_initialize_riemann_solver(flux_src_vf, norm_dir)
403
404# 73 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
405# 74 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
406# 75 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
407 if (norm_dir == 1) then
408
409# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
410
411# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
412#if defined(MFC_OpenACC)
413# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
414!$acc parallel loop collapse(3) gang vector default(present) private(alpha_rho_L, alpha_rho_R, vel, alpha_L, alpha_R, rho, pres, E, H_no_mag, gamma, pi_inf, qv, vel_rms, B, c, c_fast, pres_mag, U_L, &
415# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
416!$acc& U_R, U_starL, U_starR, U_doubleL, U_doubleR, F_L, F_R, F_starL, F_starR, F_hlld, s_L, s_R, s_M, s_starL, s_starR, pTot_L, pTot_R, p_star, rhoL_star, rhoR_star, E_starL, E_starR, sqrt_rhoL_star, &
417# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
418!$acc& sqrt_rhoR_star, denom_ds, sign_Bx, vL_star, vR_star, wL_star, wR_star, v_double, w_double, By_double, Bz_double, E_doubleL, E_doubleR, E_double) copyin(norm_dir)
419# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
420#elif defined(MFC_OpenMP)
421# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
422
423# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
424
425# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
426
427# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
428!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(alpha_rho_L, &
429# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
430!$omp& alpha_rho_R, vel, alpha_L, alpha_R, rho, pres, E, H_no_mag, gamma, pi_inf, qv, vel_rms, B, c, c_fast, pres_mag, U_L, U_R, U_starL, U_starR, U_doubleL, U_doubleR, F_L, F_R, F_starL, F_starR, &
431# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
432!$omp& F_hlld, s_L, s_R, s_M, s_starL, s_starR, pTot_L, pTot_R, p_star, rhoL_star, rhoR_star, E_starL, E_starR, sqrt_rhoL_star, sqrt_rhoR_star, denom_ds, sign_Bx, vL_star, vR_star, wL_star, wR_star, &
433# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
434!$omp& v_double, w_double, By_double, Bz_double, E_doubleL, E_doubleR, E_double) map(to:norm_dir)
435# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
436#endif
437# 82 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
438 do l = is3%beg, is3%end
439 do k = is2%beg, is2%end
440 do j = is1%beg, is1%end
441 ! (1) Extract the left/right primitive states
442 do i = 1, eqn_idx%cont%end
443 alpha_rho_l(i) = ql_prim_rsx_vf(j, k, l, i)
444 alpha_rho_r(i) = qr_prim_rsx_vf(j + 1, k, l, i)
445 end do
446
447 ! NOTE: unlike HLL & HLLC, vel_L here is permutated by dir_idx for simpler logic
448 do i = 1, num_vels
449 vel%L(i) = ql_prim_rsx_vf(j, k, l, eqn_idx%cont%end + dir_idx(i))
450 vel%R(i) = qr_prim_rsx_vf(j + 1, k, l, eqn_idx%cont%end + dir_idx(i))
451 end do
452
453 vel_rms%L = sum(vel%L**2._wp)
454 vel_rms%R = sum(vel%R**2._wp)
455
456 do i = 1, num_fluids
457 alpha_l(i) = ql_prim_rsx_vf(j, k, l, eqn_idx%E + i)
458 alpha_r(i) = qr_prim_rsx_vf(j + 1, k, l, eqn_idx%E + i)
459 end do
460
461 pres%L = ql_prim_rsx_vf(j, k, l, eqn_idx%E)
462 pres%R = qr_prim_rsx_vf(j + 1, k, l, eqn_idx%E)
463
464 ! NOTE: unlike HLL, Bx, By, Bz are permutated by dir_idx for simpler logic
465 if (mhd) then
466 if (n == 0) then ! 1D: constant Bx; By, Bz as variables; only in x so not permutated
467 b%L = [bx0, ql_prim_rsx_vf(j, k, l, eqn_idx%B%beg), ql_prim_rsx_vf(j, k, l, &
468 & eqn_idx%B%beg + 1)]
469 b%R = [bx0, qr_prim_rsx_vf(j + 1, k, l, eqn_idx%B%beg), qr_prim_rsx_vf(j + 1, k, l, &
470 & eqn_idx%B%beg + 1)]
471 else ! 2D/3D: Bx, By, Bz as variables
472 b%L = [ql_prim_rsx_vf(j, k, l, eqn_idx%B%beg + dir_idx(1) - 1), ql_prim_rsx_vf(j, k, l, &
473 & eqn_idx%B%beg + dir_idx(2) - 1), ql_prim_rsx_vf(j, k, l, &
474 & eqn_idx%B%beg + dir_idx(3) - 1)]
475 b%R = [qr_prim_rsx_vf(j + 1, k, l, eqn_idx%B%beg + dir_idx(1) - 1), &
476 & qr_prim_rsx_vf(j + 1, k, l, eqn_idx%B%beg + dir_idx(2) - 1), &
477 & qr_prim_rsx_vf(j + 1, k, l, eqn_idx%B%beg + dir_idx(3) - 1)]
478 end if
479 end if
480
481 call s_compute_mixture_coefficients(alpha_rho_l, alpha_l, rho%L, gamma%L, pi_inf%L, qv%L)
482 call s_compute_mixture_coefficients(alpha_rho_r, alpha_r, rho%R, gamma%R, pi_inf%R, qv%R)
483
484 pres_mag%L = 0.5_wp*sum(b%L**2._wp)
485 pres_mag%R = 0.5_wp*sum(b%R**2._wp)
486 call s_compute_energy(pres%L, alpha_rho_l, alpha_l, vel_rms%L, e%L)
487 e%L = e%L + pres_mag%L
488 call s_compute_energy(pres%R, alpha_rho_r, alpha_r, vel_rms%R, e%R)
489 e%R = e%R + pres_mag%R ! includes magnetic energy
490 h_no_mag%L = (e%L + pres%L - pres_mag%L)/rho%L
491 ! stagnation enthalpy here excludes magnetic energy (only used to find speed of sound)
492 h_no_mag%R = (e%R + pres%R - pres_mag%R)/rho%R
493
494 ! (2) Compute fast wave speeds
495 call s_compute_speed_of_sound(pres%L, rho%L, gamma%L, pi_inf%L, alpha_l, c%L)
496 call s_compute_speed_of_sound(pres%R, rho%R, gamma%R, pi_inf%R, alpha_r, c%R)
497 call s_compute_fast_magnetosonic_speed(rho%L, c%L, b%L, norm_dir, c_fast%L, h_no_mag%L)
498 call s_compute_fast_magnetosonic_speed(rho%R, c%R, b%R, norm_dir, c_fast%R, h_no_mag%R)
499
500 ! (3) Compute contact speed s_M [Miyoshi Equ. (38)]
501 s_l = min(vel%L(1) - c_fast%L, vel%R(1) - c_fast%R)
502 s_r = max(vel%R(1) + c_fast%R, vel%L(1) + c_fast%L)
503
504 ptot_l = pres%L + pres_mag%L
505 ptot_r = pres%R + pres_mag%R
506
507 s_m = (((s_r - vel%R(1))*rho%R*vel%R(1) - (s_l - vel%L(1))*rho%L*vel%L(1) - ptot_r + ptot_l)/((s_r &
508 & - vel%R(1))*rho%R - (s_l - vel%L(1))*rho%L))
509
510 ! (4) Compute star state variables
511 rhol_star = rho%L*(s_l - vel%L(1))/(s_l - s_m)
512 rhor_star = rho%R*(s_r - vel%R(1))/(s_r - s_m)
513 p_star = ptot_l + rho%L*(s_l - vel%L(1))*(s_m - vel%L(1))/(s_l - s_m)
514 e_starl = ((s_l - vel%L(1))*e%L - ptot_l*vel%L(1) + p_star*s_m)/(s_l - s_m)
515 e_starr = ((s_r - vel%R(1))*e%R - ptot_r*vel%R(1) + p_star*s_m)/(s_r - s_m)
516
517 ! (5) Compute left/right state vectors and fluxes
518 u_l = [rho%L, rho%L*vel%L(1:3), b%L(2:3), e%L]
519 u_starl = [rhol_star, rhol_star*s_m, rhol_star*vel%L(2:3), b%L(2:3), e_starl]
520 u_r = [rho%R, rho%R*vel%R(1:3), b%R(2:3), e%R]
521 u_starr = [rhor_star, rhor_star*s_m, rhor_star*vel%R(2:3), b%R(2:3), e_starr]
522
523 ! Compute the left/right fluxes
524 f_l(1) = u_l(2)
525 f_l(2) = u_l(2)*vel%L(1) - b%L(1)*b%L(1) + ptot_l
526 f_l(3:4) = u_l(2)*vel%L(2:3) - b%L(1)*b%L(2:3)
527 f_l(5:6) = vel%L(1)*b%L(2:3) - vel%L(2:3)*b%L(1)
528 f_l(7) = (e%L + ptot_l)*vel%L(1) - b%L(1)*(vel%L(1)*b%L(1) + vel%L(2)*b%L(2) + vel%L(3)*b%L(3))
529
530 f_r(1) = u_r(2)
531 f_r(2) = u_r(2)*vel%R(1) - b%R(1)*b%R(1) + ptot_r
532 f_r(3:4) = u_r(2)*vel%R(2:3) - b%R(1)*b%R(2:3)
533 f_r(5:6) = vel%R(1)*b%R(2:3) - vel%R(2:3)*b%R(1)
534 f_r(7) = (e%R + ptot_r)*vel%R(1) - b%R(1)*(vel%R(1)*b%R(1) + vel%R(2)*b%R(2) + vel%R(3)*b%R(3))
535 ! HLLD star-state fluxes via HLL jump relation
536 f_starl = f_l + s_l*(u_starl - u_l)
537 f_starr = f_r + s_r*(u_starr - u_r)
538 ! Alfven wave speeds bounding the rotational discontinuities
539 s_starl = s_m - abs(b%L(1))/sqrt(rhol_star)
540 s_starr = s_m + abs(b%L(1))/sqrt(rhor_star)
541 ! HLLD double-star (intermediate) states across rotational discontinuities
542 sqrt_rhol_star = sqrt(rhol_star); sqrt_rhor_star = sqrt(rhor_star)
543 vl_star = vel%L(2); wl_star = vel%L(3)
544 vr_star = vel%R(2); wr_star = vel%R(3)
545
546 ! (6) Compute the double-star states [Miyoshi Eqns. (59)-(62)]
547 denom_ds = sqrt_rhol_star + sqrt_rhor_star
548 sign_bx = sign(1._wp, b%L(1))
549 v_double = (sqrt_rhol_star*vl_star + sqrt_rhor_star*vr_star + (b%R(2) - b%L(2))*sign_bx)/denom_ds
550 w_double = (sqrt_rhol_star*wl_star + sqrt_rhor_star*wr_star + (b%R(3) - b%L(3))*sign_bx)/denom_ds
551 by_double = (sqrt_rhol_star*b%R(2) + sqrt_rhor_star*b%L(2) + sqrt_rhol_star*sqrt_rhor_star*(vr_star &
552 & - vl_star)*sign_bx)/denom_ds
553 bz_double = (sqrt_rhol_star*b%R(3) + sqrt_rhor_star*b%L(3) + sqrt_rhol_star*sqrt_rhor_star*(wr_star &
554 & - wl_star)*sign_bx)/denom_ds
555
556 e_doublel = e_starl - sqrt_rhol_star*((vl_star*b%L(2) + wl_star*b%L(3)) - (v_double*by_double &
557 & + w_double*bz_double))*sign_bx
558 e_doubler = e_starr + sqrt_rhor_star*((vr_star*b%R(2) + wr_star*b%R(3)) - (v_double*by_double &
559 & + w_double*bz_double))*sign_bx
560 e_double = 0.5_wp*(e_doublel + e_doubler)
561
562 u_doublel = [rhol_star, rhol_star*s_m, rhol_star*v_double, rhol_star*w_double, by_double, bz_double, &
563 & e_double]
564 u_doubler = [rhor_star, rhor_star*s_m, rhor_star*v_double, rhor_star*w_double, by_double, bz_double, &
565 & e_double]
566
567 ! Select HLLD flux region
568 if (0.0_wp <= s_l) then
569 f_hlld = f_l
570 else if (0.0_wp <= s_starl) then
571 f_hlld = f_l + s_l*(u_starl - u_l)
572 else if (0.0_wp <= s_m) then
573 f_hlld = f_starl + s_starl*(u_doublel - u_starl)
574 else if (0.0_wp <= s_starr) then
575 f_hlld = f_starr + s_starr*(u_doubler - u_starr)
576 else if (0.0_wp <= s_r) then
577 f_hlld = f_r + s_r*(u_starr - u_r)
578 else
579 f_hlld = f_r
580 end if
581
582 ! (12) Write HLLD flux to output arrays
583 flux_rsx_vf(j, k, l, 1) = f_hlld(1) ! TODO multi-component
584 ! Momentum
585 flux_rsx_vf(j, k, l, eqn_idx%cont%end + dir_idx(1)) = f_hlld(2)
586 flux_rsx_vf(j, k, l, eqn_idx%cont%end + dir_idx(2)) = f_hlld(3)
587 flux_rsx_vf(j, k, l, eqn_idx%cont%end + dir_idx(3)) = f_hlld(4)
588 ! Magnetic field
589 if (n == 0) then
590 flux_rsx_vf(j, k, l, eqn_idx%B%beg) = f_hlld(5)
591 flux_rsx_vf(j, k, l, eqn_idx%B%beg + 1) = f_hlld(6)
592 else
593 flux_rsx_vf(j, k, l, eqn_idx%B%beg + dir_idx(1) - 1) = 0._wp
594 flux_rsx_vf(j, k, l, eqn_idx%B%beg + dir_idx(2) - 1) = f_hlld(5)
595 flux_rsx_vf(j, k, l, eqn_idx%B%beg + dir_idx(3) - 1) = f_hlld(6)
596 end if
597 ! Energy
598 flux_rsx_vf(j, k, l, eqn_idx%E) = f_hlld(7)
599 ! Volume fractions
600
601# 244 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
602#if defined(MFC_OpenACC)
603# 244 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
604!$acc loop seq
605# 244 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
606#elif defined(MFC_OpenMP)
607# 244 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
608
609# 244 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
610#endif
611 do i = eqn_idx%adv%beg, eqn_idx%adv%end
612 flux_rsx_vf(j, k, l, i) = 0._wp ! TODO multi-component (zero for now)
613 end do
614
615 flux_src_rsx_vf(j, k, l, eqn_idx%adv%beg) = 0._wp
616 end do
617 end do
618 end do
619
620# 253 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
621#if defined(MFC_OpenACC)
622# 253 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
623!$acc end parallel loop
624# 253 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
625#elif defined(MFC_OpenMP)
626# 253 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
627
628# 253 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
629!$omp end target teams loop
630# 253 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
631#endif
632 end if
633# 73 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
634# 74 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
635# 75 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
636 if (norm_dir == 2) then
637
638# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
639
640# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
641#if defined(MFC_OpenACC)
642# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
643!$acc parallel loop collapse(3) gang vector default(present) private(alpha_rho_L, alpha_rho_R, vel, alpha_L, alpha_R, rho, pres, E, H_no_mag, gamma, pi_inf, qv, vel_rms, B, c, c_fast, pres_mag, U_L, &
644# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
645!$acc& U_R, U_starL, U_starR, U_doubleL, U_doubleR, F_L, F_R, F_starL, F_starR, F_hlld, s_L, s_R, s_M, s_starL, s_starR, pTot_L, pTot_R, p_star, rhoL_star, rhoR_star, E_starL, E_starR, sqrt_rhoL_star, &
646# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
647!$acc& sqrt_rhoR_star, denom_ds, sign_Bx, vL_star, vR_star, wL_star, wR_star, v_double, w_double, By_double, Bz_double, E_doubleL, E_doubleR, E_double) copyin(norm_dir)
648# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
649#elif defined(MFC_OpenMP)
650# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
651
652# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
653
654# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
655
656# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
657!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(alpha_rho_L, &
658# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
659!$omp& alpha_rho_R, vel, alpha_L, alpha_R, rho, pres, E, H_no_mag, gamma, pi_inf, qv, vel_rms, B, c, c_fast, pres_mag, U_L, U_R, U_starL, U_starR, U_doubleL, U_doubleR, F_L, F_R, F_starL, F_starR, &
660# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
661!$omp& F_hlld, s_L, s_R, s_M, s_starL, s_starR, pTot_L, pTot_R, p_star, rhoL_star, rhoR_star, E_starL, E_starR, sqrt_rhoL_star, sqrt_rhoR_star, denom_ds, sign_Bx, vL_star, vR_star, wL_star, wR_star, &
662# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
663!$omp& v_double, w_double, By_double, Bz_double, E_doubleL, E_doubleR, E_double) map(to:norm_dir)
664# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
665#endif
666# 82 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
667 do l = is3%beg, is3%end
668 do k = is1%beg, is1%end
669 do j = is2%beg, is2%end
670 ! (1) Extract the left/right primitive states
671 do i = 1, eqn_idx%cont%end
672 alpha_rho_l(i) = ql_prim_rsx_vf(j, k, l, i)
673 alpha_rho_r(i) = qr_prim_rsx_vf(j, k + 1, l, i)
674 end do
675
676 ! NOTE: unlike HLL & HLLC, vel_L here is permutated by dir_idx for simpler logic
677 do i = 1, num_vels
678 vel%L(i) = ql_prim_rsx_vf(j, k, l, eqn_idx%cont%end + dir_idx(i))
679 vel%R(i) = qr_prim_rsx_vf(j, k + 1, l, eqn_idx%cont%end + dir_idx(i))
680 end do
681
682 vel_rms%L = sum(vel%L**2._wp)
683 vel_rms%R = sum(vel%R**2._wp)
684
685 do i = 1, num_fluids
686 alpha_l(i) = ql_prim_rsx_vf(j, k, l, eqn_idx%E + i)
687 alpha_r(i) = qr_prim_rsx_vf(j, k + 1, l, eqn_idx%E + i)
688 end do
689
690 pres%L = ql_prim_rsx_vf(j, k, l, eqn_idx%E)
691 pres%R = qr_prim_rsx_vf(j, k + 1, l, eqn_idx%E)
692
693 ! NOTE: unlike HLL, Bx, By, Bz are permutated by dir_idx for simpler logic
694 if (mhd) then
695 if (n == 0) then ! 1D: constant Bx; By, Bz as variables; only in x so not permutated
696 b%L = [bx0, ql_prim_rsx_vf(j, k, l, eqn_idx%B%beg), ql_prim_rsx_vf(j, k, l, &
697 & eqn_idx%B%beg + 1)]
698 b%R = [bx0, qr_prim_rsx_vf(j, k + 1, l, eqn_idx%B%beg), qr_prim_rsx_vf(j, k + 1, l, &
699 & eqn_idx%B%beg + 1)]
700 else ! 2D/3D: Bx, By, Bz as variables
701 b%L = [ql_prim_rsx_vf(j, k, l, eqn_idx%B%beg + dir_idx(1) - 1), ql_prim_rsx_vf(j, k, l, &
702 & eqn_idx%B%beg + dir_idx(2) - 1), ql_prim_rsx_vf(j, k, l, &
703 & eqn_idx%B%beg + dir_idx(3) - 1)]
704 b%R = [qr_prim_rsx_vf(j, k + 1, l, eqn_idx%B%beg + dir_idx(1) - 1), &
705 & qr_prim_rsx_vf(j, k + 1, l, eqn_idx%B%beg + dir_idx(2) - 1), &
706 & qr_prim_rsx_vf(j, k + 1, l, eqn_idx%B%beg + dir_idx(3) - 1)]
707 end if
708 end if
709
710 call s_compute_mixture_coefficients(alpha_rho_l, alpha_l, rho%L, gamma%L, pi_inf%L, qv%L)
711 call s_compute_mixture_coefficients(alpha_rho_r, alpha_r, rho%R, gamma%R, pi_inf%R, qv%R)
712
713 pres_mag%L = 0.5_wp*sum(b%L**2._wp)
714 pres_mag%R = 0.5_wp*sum(b%R**2._wp)
715 call s_compute_energy(pres%L, alpha_rho_l, alpha_l, vel_rms%L, e%L)
716 e%L = e%L + pres_mag%L
717 call s_compute_energy(pres%R, alpha_rho_r, alpha_r, vel_rms%R, e%R)
718 e%R = e%R + pres_mag%R ! includes magnetic energy
719 h_no_mag%L = (e%L + pres%L - pres_mag%L)/rho%L
720 ! stagnation enthalpy here excludes magnetic energy (only used to find speed of sound)
721 h_no_mag%R = (e%R + pres%R - pres_mag%R)/rho%R
722
723 ! (2) Compute fast wave speeds
724 call s_compute_speed_of_sound(pres%L, rho%L, gamma%L, pi_inf%L, alpha_l, c%L)
725 call s_compute_speed_of_sound(pres%R, rho%R, gamma%R, pi_inf%R, alpha_r, c%R)
726 call s_compute_fast_magnetosonic_speed(rho%L, c%L, b%L, norm_dir, c_fast%L, h_no_mag%L)
727 call s_compute_fast_magnetosonic_speed(rho%R, c%R, b%R, norm_dir, c_fast%R, h_no_mag%R)
728
729 ! (3) Compute contact speed s_M [Miyoshi Equ. (38)]
730 s_l = min(vel%L(1) - c_fast%L, vel%R(1) - c_fast%R)
731 s_r = max(vel%R(1) + c_fast%R, vel%L(1) + c_fast%L)
732
733 ptot_l = pres%L + pres_mag%L
734 ptot_r = pres%R + pres_mag%R
735
736 s_m = (((s_r - vel%R(1))*rho%R*vel%R(1) - (s_l - vel%L(1))*rho%L*vel%L(1) - ptot_r + ptot_l)/((s_r &
737 & - vel%R(1))*rho%R - (s_l - vel%L(1))*rho%L))
738
739 ! (4) Compute star state variables
740 rhol_star = rho%L*(s_l - vel%L(1))/(s_l - s_m)
741 rhor_star = rho%R*(s_r - vel%R(1))/(s_r - s_m)
742 p_star = ptot_l + rho%L*(s_l - vel%L(1))*(s_m - vel%L(1))/(s_l - s_m)
743 e_starl = ((s_l - vel%L(1))*e%L - ptot_l*vel%L(1) + p_star*s_m)/(s_l - s_m)
744 e_starr = ((s_r - vel%R(1))*e%R - ptot_r*vel%R(1) + p_star*s_m)/(s_r - s_m)
745
746 ! (5) Compute left/right state vectors and fluxes
747 u_l = [rho%L, rho%L*vel%L(1:3), b%L(2:3), e%L]
748 u_starl = [rhol_star, rhol_star*s_m, rhol_star*vel%L(2:3), b%L(2:3), e_starl]
749 u_r = [rho%R, rho%R*vel%R(1:3), b%R(2:3), e%R]
750 u_starr = [rhor_star, rhor_star*s_m, rhor_star*vel%R(2:3), b%R(2:3), e_starr]
751
752 ! Compute the left/right fluxes
753 f_l(1) = u_l(2)
754 f_l(2) = u_l(2)*vel%L(1) - b%L(1)*b%L(1) + ptot_l
755 f_l(3:4) = u_l(2)*vel%L(2:3) - b%L(1)*b%L(2:3)
756 f_l(5:6) = vel%L(1)*b%L(2:3) - vel%L(2:3)*b%L(1)
757 f_l(7) = (e%L + ptot_l)*vel%L(1) - b%L(1)*(vel%L(1)*b%L(1) + vel%L(2)*b%L(2) + vel%L(3)*b%L(3))
758
759 f_r(1) = u_r(2)
760 f_r(2) = u_r(2)*vel%R(1) - b%R(1)*b%R(1) + ptot_r
761 f_r(3:4) = u_r(2)*vel%R(2:3) - b%R(1)*b%R(2:3)
762 f_r(5:6) = vel%R(1)*b%R(2:3) - vel%R(2:3)*b%R(1)
763 f_r(7) = (e%R + ptot_r)*vel%R(1) - b%R(1)*(vel%R(1)*b%R(1) + vel%R(2)*b%R(2) + vel%R(3)*b%R(3))
764 ! HLLD star-state fluxes via HLL jump relation
765 f_starl = f_l + s_l*(u_starl - u_l)
766 f_starr = f_r + s_r*(u_starr - u_r)
767 ! Alfven wave speeds bounding the rotational discontinuities
768 s_starl = s_m - abs(b%L(1))/sqrt(rhol_star)
769 s_starr = s_m + abs(b%L(1))/sqrt(rhor_star)
770 ! HLLD double-star (intermediate) states across rotational discontinuities
771 sqrt_rhol_star = sqrt(rhol_star); sqrt_rhor_star = sqrt(rhor_star)
772 vl_star = vel%L(2); wl_star = vel%L(3)
773 vr_star = vel%R(2); wr_star = vel%R(3)
774
775 ! (6) Compute the double-star states [Miyoshi Eqns. (59)-(62)]
776 denom_ds = sqrt_rhol_star + sqrt_rhor_star
777 sign_bx = sign(1._wp, b%L(1))
778 v_double = (sqrt_rhol_star*vl_star + sqrt_rhor_star*vr_star + (b%R(2) - b%L(2))*sign_bx)/denom_ds
779 w_double = (sqrt_rhol_star*wl_star + sqrt_rhor_star*wr_star + (b%R(3) - b%L(3))*sign_bx)/denom_ds
780 by_double = (sqrt_rhol_star*b%R(2) + sqrt_rhor_star*b%L(2) + sqrt_rhol_star*sqrt_rhor_star*(vr_star &
781 & - vl_star)*sign_bx)/denom_ds
782 bz_double = (sqrt_rhol_star*b%R(3) + sqrt_rhor_star*b%L(3) + sqrt_rhol_star*sqrt_rhor_star*(wr_star &
783 & - wl_star)*sign_bx)/denom_ds
784
785 e_doublel = e_starl - sqrt_rhol_star*((vl_star*b%L(2) + wl_star*b%L(3)) - (v_double*by_double &
786 & + w_double*bz_double))*sign_bx
787 e_doubler = e_starr + sqrt_rhor_star*((vr_star*b%R(2) + wr_star*b%R(3)) - (v_double*by_double &
788 & + w_double*bz_double))*sign_bx
789 e_double = 0.5_wp*(e_doublel + e_doubler)
790
791 u_doublel = [rhol_star, rhol_star*s_m, rhol_star*v_double, rhol_star*w_double, by_double, bz_double, &
792 & e_double]
793 u_doubler = [rhor_star, rhor_star*s_m, rhor_star*v_double, rhor_star*w_double, by_double, bz_double, &
794 & e_double]
795
796 ! Select HLLD flux region
797 if (0.0_wp <= s_l) then
798 f_hlld = f_l
799 else if (0.0_wp <= s_starl) then
800 f_hlld = f_l + s_l*(u_starl - u_l)
801 else if (0.0_wp <= s_m) then
802 f_hlld = f_starl + s_starl*(u_doublel - u_starl)
803 else if (0.0_wp <= s_starr) then
804 f_hlld = f_starr + s_starr*(u_doubler - u_starr)
805 else if (0.0_wp <= s_r) then
806 f_hlld = f_r + s_r*(u_starr - u_r)
807 else
808 f_hlld = f_r
809 end if
810
811 ! (12) Write HLLD flux to output arrays
812 flux_rsx_vf(j, k, l, 1) = f_hlld(1) ! TODO multi-component
813 ! Momentum
814 flux_rsx_vf(j, k, l, eqn_idx%cont%end + dir_idx(1)) = f_hlld(2)
815 flux_rsx_vf(j, k, l, eqn_idx%cont%end + dir_idx(2)) = f_hlld(3)
816 flux_rsx_vf(j, k, l, eqn_idx%cont%end + dir_idx(3)) = f_hlld(4)
817 ! Magnetic field
818 if (n == 0) then
819 flux_rsx_vf(j, k, l, eqn_idx%B%beg) = f_hlld(5)
820 flux_rsx_vf(j, k, l, eqn_idx%B%beg + 1) = f_hlld(6)
821 else
822 flux_rsx_vf(j, k, l, eqn_idx%B%beg + dir_idx(1) - 1) = 0._wp
823 flux_rsx_vf(j, k, l, eqn_idx%B%beg + dir_idx(2) - 1) = f_hlld(5)
824 flux_rsx_vf(j, k, l, eqn_idx%B%beg + dir_idx(3) - 1) = f_hlld(6)
825 end if
826 ! Energy
827 flux_rsx_vf(j, k, l, eqn_idx%E) = f_hlld(7)
828 ! Volume fractions
829
830# 244 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
831#if defined(MFC_OpenACC)
832# 244 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
833!$acc loop seq
834# 244 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
835#elif defined(MFC_OpenMP)
836# 244 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
837
838# 244 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
839#endif
840 do i = eqn_idx%adv%beg, eqn_idx%adv%end
841 flux_rsx_vf(j, k, l, i) = 0._wp ! TODO multi-component (zero for now)
842 end do
843
844 flux_src_rsx_vf(j, k, l, eqn_idx%adv%beg) = 0._wp
845 end do
846 end do
847 end do
848
849# 253 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
850#if defined(MFC_OpenACC)
851# 253 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
852!$acc end parallel loop
853# 253 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
854#elif defined(MFC_OpenMP)
855# 253 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
856
857# 253 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
858!$omp end target teams loop
859# 253 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
860#endif
861 end if
862# 73 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
863# 74 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
864# 75 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
865 if (norm_dir == 3) then
866
867# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
868
869# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
870#if defined(MFC_OpenACC)
871# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
872!$acc parallel loop collapse(3) gang vector default(present) private(alpha_rho_L, alpha_rho_R, vel, alpha_L, alpha_R, rho, pres, E, H_no_mag, gamma, pi_inf, qv, vel_rms, B, c, c_fast, pres_mag, U_L, &
873# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
874!$acc& U_R, U_starL, U_starR, U_doubleL, U_doubleR, F_L, F_R, F_starL, F_starR, F_hlld, s_L, s_R, s_M, s_starL, s_starR, pTot_L, pTot_R, p_star, rhoL_star, rhoR_star, E_starL, E_starR, sqrt_rhoL_star, &
875# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
876!$acc& sqrt_rhoR_star, denom_ds, sign_Bx, vL_star, vR_star, wL_star, wR_star, v_double, w_double, By_double, Bz_double, E_doubleL, E_doubleR, E_double) copyin(norm_dir)
877# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
878#elif defined(MFC_OpenMP)
879# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
880
881# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
882
883# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
884
885# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
886!$omp target teams loop defaultmap(firstprivate:scalar) bind(teams,parallel) collapse(3) defaultmap(tofrom:aggregate) defaultmap(tofrom:allocatable) defaultmap(tofrom:pointer) private(alpha_rho_L, &
887# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
888!$omp& alpha_rho_R, vel, alpha_L, alpha_R, rho, pres, E, H_no_mag, gamma, pi_inf, qv, vel_rms, B, c, c_fast, pres_mag, U_L, U_R, U_starL, U_starR, U_doubleL, U_doubleR, F_L, F_R, F_starL, F_starR, &
889# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
890!$omp& F_hlld, s_L, s_R, s_M, s_starL, s_starR, pTot_L, pTot_R, p_star, rhoL_star, rhoR_star, E_starL, E_starR, sqrt_rhoL_star, sqrt_rhoR_star, denom_ds, sign_Bx, vL_star, vR_star, wL_star, wR_star, &
891# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
892!$omp& v_double, w_double, By_double, Bz_double, E_doubleL, E_doubleR, E_double) map(to:norm_dir)
893# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
894#endif
895# 82 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
896 do l = is1%beg, is1%end
897 do k = is2%beg, is2%end
898 do j = is3%beg, is3%end
899 ! (1) Extract the left/right primitive states
900 do i = 1, eqn_idx%cont%end
901 alpha_rho_l(i) = ql_prim_rsx_vf(j, k, l, i)
902 alpha_rho_r(i) = qr_prim_rsx_vf(j, k, l + 1, i)
903 end do
904
905 ! NOTE: unlike HLL & HLLC, vel_L here is permutated by dir_idx for simpler logic
906 do i = 1, num_vels
907 vel%L(i) = ql_prim_rsx_vf(j, k, l, eqn_idx%cont%end + dir_idx(i))
908 vel%R(i) = qr_prim_rsx_vf(j, k, l + 1, eqn_idx%cont%end + dir_idx(i))
909 end do
910
911 vel_rms%L = sum(vel%L**2._wp)
912 vel_rms%R = sum(vel%R**2._wp)
913
914 do i = 1, num_fluids
915 alpha_l(i) = ql_prim_rsx_vf(j, k, l, eqn_idx%E + i)
916 alpha_r(i) = qr_prim_rsx_vf(j, k, l + 1, eqn_idx%E + i)
917 end do
918
919 pres%L = ql_prim_rsx_vf(j, k, l, eqn_idx%E)
920 pres%R = qr_prim_rsx_vf(j, k, l + 1, eqn_idx%E)
921
922 ! NOTE: unlike HLL, Bx, By, Bz are permutated by dir_idx for simpler logic
923 if (mhd) then
924 if (n == 0) then ! 1D: constant Bx; By, Bz as variables; only in x so not permutated
925 b%L = [bx0, ql_prim_rsx_vf(j, k, l, eqn_idx%B%beg), ql_prim_rsx_vf(j, k, l, &
926 & eqn_idx%B%beg + 1)]
927 b%R = [bx0, qr_prim_rsx_vf(j, k, l + 1, eqn_idx%B%beg), qr_prim_rsx_vf(j, k, l + 1, &
928 & eqn_idx%B%beg + 1)]
929 else ! 2D/3D: Bx, By, Bz as variables
930 b%L = [ql_prim_rsx_vf(j, k, l, eqn_idx%B%beg + dir_idx(1) - 1), ql_prim_rsx_vf(j, k, l, &
931 & eqn_idx%B%beg + dir_idx(2) - 1), ql_prim_rsx_vf(j, k, l, &
932 & eqn_idx%B%beg + dir_idx(3) - 1)]
933 b%R = [qr_prim_rsx_vf(j, k, l + 1, eqn_idx%B%beg + dir_idx(1) - 1), &
934 & qr_prim_rsx_vf(j, k, l + 1, eqn_idx%B%beg + dir_idx(2) - 1), &
935 & qr_prim_rsx_vf(j, k, l + 1, eqn_idx%B%beg + dir_idx(3) - 1)]
936 end if
937 end if
938
939 call s_compute_mixture_coefficients(alpha_rho_l, alpha_l, rho%L, gamma%L, pi_inf%L, qv%L)
940 call s_compute_mixture_coefficients(alpha_rho_r, alpha_r, rho%R, gamma%R, pi_inf%R, qv%R)
941
942 pres_mag%L = 0.5_wp*sum(b%L**2._wp)
943 pres_mag%R = 0.5_wp*sum(b%R**2._wp)
944 call s_compute_energy(pres%L, alpha_rho_l, alpha_l, vel_rms%L, e%L)
945 e%L = e%L + pres_mag%L
946 call s_compute_energy(pres%R, alpha_rho_r, alpha_r, vel_rms%R, e%R)
947 e%R = e%R + pres_mag%R ! includes magnetic energy
948 h_no_mag%L = (e%L + pres%L - pres_mag%L)/rho%L
949 ! stagnation enthalpy here excludes magnetic energy (only used to find speed of sound)
950 h_no_mag%R = (e%R + pres%R - pres_mag%R)/rho%R
951
952 ! (2) Compute fast wave speeds
953 call s_compute_speed_of_sound(pres%L, rho%L, gamma%L, pi_inf%L, alpha_l, c%L)
954 call s_compute_speed_of_sound(pres%R, rho%R, gamma%R, pi_inf%R, alpha_r, c%R)
955 call s_compute_fast_magnetosonic_speed(rho%L, c%L, b%L, norm_dir, c_fast%L, h_no_mag%L)
956 call s_compute_fast_magnetosonic_speed(rho%R, c%R, b%R, norm_dir, c_fast%R, h_no_mag%R)
957
958 ! (3) Compute contact speed s_M [Miyoshi Equ. (38)]
959 s_l = min(vel%L(1) - c_fast%L, vel%R(1) - c_fast%R)
960 s_r = max(vel%R(1) + c_fast%R, vel%L(1) + c_fast%L)
961
962 ptot_l = pres%L + pres_mag%L
963 ptot_r = pres%R + pres_mag%R
964
965 s_m = (((s_r - vel%R(1))*rho%R*vel%R(1) - (s_l - vel%L(1))*rho%L*vel%L(1) - ptot_r + ptot_l)/((s_r &
966 & - vel%R(1))*rho%R - (s_l - vel%L(1))*rho%L))
967
968 ! (4) Compute star state variables
969 rhol_star = rho%L*(s_l - vel%L(1))/(s_l - s_m)
970 rhor_star = rho%R*(s_r - vel%R(1))/(s_r - s_m)
971 p_star = ptot_l + rho%L*(s_l - vel%L(1))*(s_m - vel%L(1))/(s_l - s_m)
972 e_starl = ((s_l - vel%L(1))*e%L - ptot_l*vel%L(1) + p_star*s_m)/(s_l - s_m)
973 e_starr = ((s_r - vel%R(1))*e%R - ptot_r*vel%R(1) + p_star*s_m)/(s_r - s_m)
974
975 ! (5) Compute left/right state vectors and fluxes
976 u_l = [rho%L, rho%L*vel%L(1:3), b%L(2:3), e%L]
977 u_starl = [rhol_star, rhol_star*s_m, rhol_star*vel%L(2:3), b%L(2:3), e_starl]
978 u_r = [rho%R, rho%R*vel%R(1:3), b%R(2:3), e%R]
979 u_starr = [rhor_star, rhor_star*s_m, rhor_star*vel%R(2:3), b%R(2:3), e_starr]
980
981 ! Compute the left/right fluxes
982 f_l(1) = u_l(2)
983 f_l(2) = u_l(2)*vel%L(1) - b%L(1)*b%L(1) + ptot_l
984 f_l(3:4) = u_l(2)*vel%L(2:3) - b%L(1)*b%L(2:3)
985 f_l(5:6) = vel%L(1)*b%L(2:3) - vel%L(2:3)*b%L(1)
986 f_l(7) = (e%L + ptot_l)*vel%L(1) - b%L(1)*(vel%L(1)*b%L(1) + vel%L(2)*b%L(2) + vel%L(3)*b%L(3))
987
988 f_r(1) = u_r(2)
989 f_r(2) = u_r(2)*vel%R(1) - b%R(1)*b%R(1) + ptot_r
990 f_r(3:4) = u_r(2)*vel%R(2:3) - b%R(1)*b%R(2:3)
991 f_r(5:6) = vel%R(1)*b%R(2:3) - vel%R(2:3)*b%R(1)
992 f_r(7) = (e%R + ptot_r)*vel%R(1) - b%R(1)*(vel%R(1)*b%R(1) + vel%R(2)*b%R(2) + vel%R(3)*b%R(3))
993 ! HLLD star-state fluxes via HLL jump relation
994 f_starl = f_l + s_l*(u_starl - u_l)
995 f_starr = f_r + s_r*(u_starr - u_r)
996 ! Alfven wave speeds bounding the rotational discontinuities
997 s_starl = s_m - abs(b%L(1))/sqrt(rhol_star)
998 s_starr = s_m + abs(b%L(1))/sqrt(rhor_star)
999 ! HLLD double-star (intermediate) states across rotational discontinuities
1000 sqrt_rhol_star = sqrt(rhol_star); sqrt_rhor_star = sqrt(rhor_star)
1001 vl_star = vel%L(2); wl_star = vel%L(3)
1002 vr_star = vel%R(2); wr_star = vel%R(3)
1003
1004 ! (6) Compute the double-star states [Miyoshi Eqns. (59)-(62)]
1005 denom_ds = sqrt_rhol_star + sqrt_rhor_star
1006 sign_bx = sign(1._wp, b%L(1))
1007 v_double = (sqrt_rhol_star*vl_star + sqrt_rhor_star*vr_star + (b%R(2) - b%L(2))*sign_bx)/denom_ds
1008 w_double = (sqrt_rhol_star*wl_star + sqrt_rhor_star*wr_star + (b%R(3) - b%L(3))*sign_bx)/denom_ds
1009 by_double = (sqrt_rhol_star*b%R(2) + sqrt_rhor_star*b%L(2) + sqrt_rhol_star*sqrt_rhor_star*(vr_star &
1010 & - vl_star)*sign_bx)/denom_ds
1011 bz_double = (sqrt_rhol_star*b%R(3) + sqrt_rhor_star*b%L(3) + sqrt_rhol_star*sqrt_rhor_star*(wr_star &
1012 & - wl_star)*sign_bx)/denom_ds
1013
1014 e_doublel = e_starl - sqrt_rhol_star*((vl_star*b%L(2) + wl_star*b%L(3)) - (v_double*by_double &
1015 & + w_double*bz_double))*sign_bx
1016 e_doubler = e_starr + sqrt_rhor_star*((vr_star*b%R(2) + wr_star*b%R(3)) - (v_double*by_double &
1017 & + w_double*bz_double))*sign_bx
1018 e_double = 0.5_wp*(e_doublel + e_doubler)
1019
1020 u_doublel = [rhol_star, rhol_star*s_m, rhol_star*v_double, rhol_star*w_double, by_double, bz_double, &
1021 & e_double]
1022 u_doubler = [rhor_star, rhor_star*s_m, rhor_star*v_double, rhor_star*w_double, by_double, bz_double, &
1023 & e_double]
1024
1025 ! Select HLLD flux region
1026 if (0.0_wp <= s_l) then
1027 f_hlld = f_l
1028 else if (0.0_wp <= s_starl) then
1029 f_hlld = f_l + s_l*(u_starl - u_l)
1030 else if (0.0_wp <= s_m) then
1031 f_hlld = f_starl + s_starl*(u_doublel - u_starl)
1032 else if (0.0_wp <= s_starr) then
1033 f_hlld = f_starr + s_starr*(u_doubler - u_starr)
1034 else if (0.0_wp <= s_r) then
1035 f_hlld = f_r + s_r*(u_starr - u_r)
1036 else
1037 f_hlld = f_r
1038 end if
1039
1040 ! (12) Write HLLD flux to output arrays
1041 flux_rsx_vf(j, k, l, 1) = f_hlld(1) ! TODO multi-component
1042 ! Momentum
1043 flux_rsx_vf(j, k, l, eqn_idx%cont%end + dir_idx(1)) = f_hlld(2)
1044 flux_rsx_vf(j, k, l, eqn_idx%cont%end + dir_idx(2)) = f_hlld(3)
1045 flux_rsx_vf(j, k, l, eqn_idx%cont%end + dir_idx(3)) = f_hlld(4)
1046 ! Magnetic field
1047 if (n == 0) then
1048 flux_rsx_vf(j, k, l, eqn_idx%B%beg) = f_hlld(5)
1049 flux_rsx_vf(j, k, l, eqn_idx%B%beg + 1) = f_hlld(6)
1050 else
1051 flux_rsx_vf(j, k, l, eqn_idx%B%beg + dir_idx(1) - 1) = 0._wp
1052 flux_rsx_vf(j, k, l, eqn_idx%B%beg + dir_idx(2) - 1) = f_hlld(5)
1053 flux_rsx_vf(j, k, l, eqn_idx%B%beg + dir_idx(3) - 1) = f_hlld(6)
1054 end if
1055 ! Energy
1056 flux_rsx_vf(j, k, l, eqn_idx%E) = f_hlld(7)
1057 ! Volume fractions
1058
1059# 244 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
1060#if defined(MFC_OpenACC)
1061# 244 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
1062!$acc loop seq
1063# 244 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
1064#elif defined(MFC_OpenMP)
1065# 244 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
1066
1067# 244 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
1068#endif
1069 do i = eqn_idx%adv%beg, eqn_idx%adv%end
1070 flux_rsx_vf(j, k, l, i) = 0._wp ! TODO multi-component (zero for now)
1071 end do
1072
1073 flux_src_rsx_vf(j, k, l, eqn_idx%adv%beg) = 0._wp
1074 end do
1075 end do
1076 end do
1077
1078# 253 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
1079#if defined(MFC_OpenACC)
1080# 253 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
1081!$acc end parallel loop
1082# 253 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
1083#elif defined(MFC_OpenMP)
1084# 253 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
1085
1086# 253 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
1087!$omp end target teams loop
1088# 253 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
1089#endif
1090 end if
1091# 256 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
1092
1093 call s_finalize_riemann_solver(flux_vf, flux_src_vf, flux_gsrc_vf, norm_dir)
1094
1095 end subroutine s_hlld_riemann_solver
1096
1097end module m_riemann_solver_hlld
m_derived_types
Shared derived types for field data, patch geometry, bubble dynamics, and MPI I/O structures.
Definition m_derived_types.fpp.f90:335
m_global_parameters
Global parameters for the computational domain, fluid properties, and simulation algorithm configurat...
Definition m_global_parameters.fpp.f90:346
m_global_parameters::dir_idx
integer, dimension(3) dir_idx
Definition m_global_parameters.fpp.f90:787
m_riemann_solver_hlld
HLLD approximate Riemann solver for MHD, Miyoshi & Kusano JCP (2005).
Definition m_riemann_solver_hlld.fpp.f90:346
m_riemann_solver_hlld::s_hlld_riemann_solver
subroutine s_hlld_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)
HLLD Riemann solver for MHD, Miyoshi & Kusano JCP (2005).
Definition m_riemann_solver_hlld.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:347
m_riemann_state::flux_src_rsx_vf
real(wp), dimension(:,:,:,:), allocatable flux_src_rsx_vf
Definition m_riemann_state.fpp.f90:360
m_riemann_state::is3
type(int_bounds_info) is3
Definition m_riemann_state.fpp.f90:489
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:360
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:2088
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:1006
m_riemann_state::is2
type(int_bounds_info) is2
Definition m_riemann_state.fpp.f90:489
m_riemann_state::s_finalize_riemann_solver
subroutine s_finalize_riemann_solver(flux_vf, flux_src_vf, flux_gsrc_vf, norm_dir)
Reshape and copy the Riemann-solver flux buffers back to the physical-space output arrays for the sel...
Definition m_riemann_state.fpp.f90:3311
m_riemann_state::is1
type(int_bounds_info) is1
Definition m_riemann_state.fpp.f90:489
m_variables_conversion
Conservative-to-primitive variable conversion, mixture property evaluation, and pressure computation.
Definition m_variables_conversion.fpp.f90:346
m_variables_conversion::s_compute_mixture_coefficients
subroutine, public s_compute_mixture_coefficients(alpha_rho_k, alpha_k, rho_k, gamma_k, pi_inf_k, qv_k)
Mixture coefficients of one state. Under bubbles_euler with num_fluids == 1 the sole advection slot a...
Definition m_variables_conversion.fpp.f90:2683
m_variables_conversion::s_compute_energy
subroutine, public s_compute_energy(pres, alpha_rho_k, alpha_k, vel_sum, e)
Total energy per unit volume, thermodynamic terms only. Callers add magnetic and elastic energy,...
Definition m_variables_conversion.fpp.f90:2852
m_variables_conversion::s_compute_speed_of_sound
subroutine, public s_compute_speed_of_sound(pres, rho, gamma, pi_inf, adv, c, alpha_rho)
Speed of sound of a thermodynamic state. Enthalpy is not an argument: for a real state H,...
Definition m_variables_conversion.fpp.f90:4023
m_variables_conversion::s_compute_fast_magnetosonic_speed
subroutine, public 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:4201
m_derived_types::int_bounds_info
Integer bounds for variables.
Definition m_derived_types.fpp.f90:429
m_derived_types::riemann_states_vec3
Left and right Riemann states for 3-component vectors.
Definition m_derived_types.fpp.f90:417
m_derived_types::riemann_states
Left and right Riemann states.
Definition m_derived_types.fpp.f90:411
m_derived_types::scalar_field
Derived type annexing a scalar field (SF).
Definition m_derived_types.fpp.f90:349