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