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