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# 167 "/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# 167 "/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# 167 "/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.
311# 55 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
312
313! Allocate and create GPU device memory
314# 75 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
315
316! Free GPU device memory and deallocate
317# 83 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
318
319! Cray-specific GPU pointer setup for vector fields
320# 107 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
321
322! Cray-specific GPU pointer setup for scalar fields
323# 123 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
324
325! Cray-specific GPU pointer setup for acoustic source spatials
326# 148 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
327
328# 154 "/home/runner/work/MFC/MFC/src/common/include/macros.fpp"
329
330# 161 "/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) 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, &
403# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
404!$acc& U_R, U_starL, U_starR, U_doubleL, U_doubleR, F_L, F_R, F_starL, F_starR, F_hlld, s_L, s_R, s_M, s_starL, s_starR, pTot_L, pTot_R, p_star, rhoL_star, rhoR_star, E_starL, E_starR, sqrt_rhoL_star, &
405# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
406!$acc& sqrt_rhoR_star, denom_ds, sign_Bx, vL_star, vR_star, wL_star, wR_star, v_double, w_double, By_double, Bz_double, E_doubleL, E_doubleR, E_double) copyin(norm_dir)
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) private(alpha_rho_L, &
417# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
418!$omp& alpha_rho_R, vel, alpha_L, alpha_R, rho, pres, E, H_no_mag, gamma, pi_inf, qv, vel_rms, B, c, c_fast, pres_mag, U_L, U_R, U_starL, U_starR, U_doubleL, U_doubleR, F_L, F_R, F_starL, F_starR, &
419# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
420!$omp& F_hlld, s_L, s_R, s_M, s_starL, s_starR, pTot_L, pTot_R, p_star, rhoL_star, rhoR_star, E_starL, E_starR, sqrt_rhoL_star, sqrt_rhoR_star, denom_ds, sign_Bx, vL_star, vR_star, wL_star, wR_star, &
421# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
422!$omp& v_double, w_double, By_double, Bz_double, E_doubleL, E_doubleR, E_double) map(to:norm_dir)
423# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
424#endif
425# 82 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
426 do l = is3%beg, is3%end
427 do k = is2%beg, is2%end
428 do j = is1%beg, is1%end
429 ! (1) Extract the left/right primitive states
430 do i = 1, eqn_idx%cont%end
431 alpha_rho_l(i) = ql_prim_rsx_vf(j, k, l, i)
432 alpha_rho_r(i) = qr_prim_rsx_vf(j + 1, k, l, i)
433 end do
434
435 ! NOTE: unlike HLL & HLLC, vel_L here is permutated by dir_idx for simpler logic
436 do i = 1, num_vels
437 vel%L(i) = ql_prim_rsx_vf(j, k, l, eqn_idx%cont%end + dir_idx(i))
438 vel%R(i) = qr_prim_rsx_vf(j + 1, k, l, eqn_idx%cont%end + dir_idx(i))
439 end do
440
441 vel_rms%L = sum(vel%L**2._wp)
442 vel_rms%R = sum(vel%R**2._wp)
443
444 do i = 1, num_fluids
445 alpha_l(i) = ql_prim_rsx_vf(j, k, l, eqn_idx%E + i)
446 alpha_r(i) = qr_prim_rsx_vf(j + 1, k, l, eqn_idx%E + i)
447 end do
448
449 pres%L = ql_prim_rsx_vf(j, k, l, eqn_idx%E)
450 pres%R = qr_prim_rsx_vf(j + 1, k, l, eqn_idx%E)
451
452 ! NOTE: unlike HLL, Bx, By, Bz are permutated by dir_idx for simpler logic
453 if (mhd) then
454 if (n == 0) then ! 1D: constant Bx; By, Bz as variables; only in x so not permutated
455 b%L = [bx0, ql_prim_rsx_vf(j, k, l, eqn_idx%B%beg), ql_prim_rsx_vf(j, k, l, &
456 & eqn_idx%B%beg + 1)]
457 b%R = [bx0, qr_prim_rsx_vf(j + 1, k, l, eqn_idx%B%beg), qr_prim_rsx_vf(j + 1, k, l, &
458 & eqn_idx%B%beg + 1)]
459 else ! 2D/3D: Bx, By, Bz as variables
460 b%L = [ql_prim_rsx_vf(j, k, l, eqn_idx%B%beg + dir_idx(1) - 1), ql_prim_rsx_vf(j, k, l, &
461 & eqn_idx%B%beg + dir_idx(2) - 1), ql_prim_rsx_vf(j, k, l, &
462 & eqn_idx%B%beg + dir_idx(3) - 1)]
463 b%R = [qr_prim_rsx_vf(j + 1, k, l, eqn_idx%B%beg + dir_idx(1) - 1), &
464 & qr_prim_rsx_vf(j + 1, k, l, eqn_idx%B%beg + dir_idx(2) - 1), &
465 & qr_prim_rsx_vf(j + 1, k, l, eqn_idx%B%beg + dir_idx(3) - 1)]
466 end if
467 end if
468
469 ! Sum properties of all fluid components
470 rho%L = 0._wp; gamma%L = 0._wp; pi_inf%L = 0._wp; qv%L = 0._wp
471 rho%R = 0._wp; gamma%R = 0._wp; pi_inf%R = 0._wp; qv%R = 0._wp
472
473# 128 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
474#if defined(MFC_OpenACC)
475# 128 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
476!$acc loop seq
477# 128 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
478#elif defined(MFC_OpenMP)
479# 128 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
480
481# 128 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
482#endif
483 do i = 1, num_fluids
484 rho%L = rho%L + alpha_rho_l(i)
485 gamma%L = gamma%L + alpha_l(i)*gammas(i)
486 pi_inf%L = pi_inf%L + alpha_l(i)*pi_infs(i)
487 qv%L = qv%L + alpha_rho_l(i)*qvs(i)
488
489 rho%R = rho%R + alpha_rho_r(i)
490 gamma%R = gamma%R + alpha_r(i)*gammas(i)
491 pi_inf%R = pi_inf%R + alpha_r(i)*pi_infs(i)
492 qv%R = qv%R + alpha_rho_r(i)*qvs(i)
493 end do
494
495 pres_mag%L = 0.5_wp*sum(b%L**2._wp)
496 pres_mag%R = 0.5_wp*sum(b%R**2._wp)
497 e%L = gamma%L*pres%L + pi_inf%L + 0.5_wp*rho%L*vel_rms%L + qv%L + pres_mag%L
498 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
499 h_no_mag%L = (e%L + pres%L - pres_mag%L)/rho%L
500 ! stagnation enthalpy here excludes magnetic energy (only used to find speed of sound)
501 h_no_mag%R = (e%R + pres%R - pres_mag%R)/rho%R
502
503 ! (2) Compute fast wave speeds
504 call s_compute_speed_of_sound(pres%L, rho%L, gamma%L, pi_inf%L, h_no_mag%L, alpha_l, vel_rms%L, &
505 & 0._wp, c%L, qv%L)
506 call s_compute_speed_of_sound(pres%R, rho%R, gamma%R, pi_inf%R, h_no_mag%R, alpha_r, vel_rms%R, &
507 & 0._wp, c%R, qv%R)
508 call s_compute_fast_magnetosonic_speed(rho%L, c%L, b%L, norm_dir, c_fast%L, h_no_mag%L)
509 call s_compute_fast_magnetosonic_speed(rho%R, c%R, b%R, norm_dir, c_fast%R, h_no_mag%R)
510
511 ! (3) Compute contact speed s_M [Miyoshi Equ. (38)]
512 s_l = min(vel%L(1) - c_fast%L, vel%R(1) - c_fast%R)
513 s_r = max(vel%R(1) + c_fast%R, vel%L(1) + c_fast%L)
514
515 ptot_l = pres%L + pres_mag%L
516 ptot_r = pres%R + pres_mag%R
517
518 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 &
519 & - vel%R(1))*rho%R - (s_l - vel%L(1))*rho%L))
520
521 ! (4) Compute star state variables
522 rhol_star = rho%L*(s_l - vel%L(1))/(s_l - s_m)
523 rhor_star = rho%R*(s_r - vel%R(1))/(s_r - s_m)
524 p_star = ptot_l + rho%L*(s_l - vel%L(1))*(s_m - vel%L(1))/(s_l - s_m)
525 e_starl = ((s_l - vel%L(1))*e%L - ptot_l*vel%L(1) + p_star*s_m)/(s_l - s_m)
526 e_starr = ((s_r - vel%R(1))*e%R - ptot_r*vel%R(1) + p_star*s_m)/(s_r - s_m)
527
528 ! (5) Compute left/right state vectors and fluxes
529 u_l = [rho%L, rho%L*vel%L(1:3), b%L(2:3), e%L]
530 u_starl = [rhol_star, rhol_star*s_m, rhol_star*vel%L(2:3), b%L(2:3), e_starl]
531 u_r = [rho%R, rho%R*vel%R(1:3), b%R(2:3), e%R]
532 u_starr = [rhor_star, rhor_star*s_m, rhor_star*vel%R(2:3), b%R(2:3), e_starr]
533
534 ! Compute the left/right fluxes
535 f_l(1) = u_l(2)
536 f_l(2) = u_l(2)*vel%L(1) - b%L(1)*b%L(1) + ptot_l
537 f_l(3:4) = u_l(2)*vel%L(2:3) - b%L(1)*b%L(2:3)
538 f_l(5:6) = vel%L(1)*b%L(2:3) - vel%L(2:3)*b%L(1)
539 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))
540
541 f_r(1) = u_r(2)
542 f_r(2) = u_r(2)*vel%R(1) - b%R(1)*b%R(1) + ptot_r
543 f_r(3:4) = u_r(2)*vel%R(2:3) - b%R(1)*b%R(2:3)
544 f_r(5:6) = vel%R(1)*b%R(2:3) - vel%R(2:3)*b%R(1)
545 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))
546 ! HLLD star-state fluxes via HLL jump relation
547 f_starl = f_l + s_l*(u_starl - u_l)
548 f_starr = f_r + s_r*(u_starr - u_r)
549 ! Alfven wave speeds bounding the rotational discontinuities
550 s_starl = s_m - abs(b%L(1))/sqrt(rhol_star)
551 s_starr = s_m + abs(b%L(1))/sqrt(rhor_star)
552 ! HLLD double-star (intermediate) states across rotational discontinuities
553 sqrt_rhol_star = sqrt(rhol_star); sqrt_rhor_star = sqrt(rhor_star)
554 vl_star = vel%L(2); wl_star = vel%L(3)
555 vr_star = vel%R(2); wr_star = vel%R(3)
556
557 ! (6) Compute the double-star states [Miyoshi Eqns. (59)-(62)]
558 denom_ds = sqrt_rhol_star + sqrt_rhor_star
559 sign_bx = sign(1._wp, b%L(1))
560 v_double = (sqrt_rhol_star*vl_star + sqrt_rhor_star*vr_star + (b%R(2) - b%L(2))*sign_bx)/denom_ds
561 w_double = (sqrt_rhol_star*wl_star + sqrt_rhor_star*wr_star + (b%R(3) - b%L(3))*sign_bx)/denom_ds
562 by_double = (sqrt_rhol_star*b%R(2) + sqrt_rhor_star*b%L(2) + sqrt_rhol_star*sqrt_rhor_star*(vr_star &
563 & - vl_star)*sign_bx)/denom_ds
564 bz_double = (sqrt_rhol_star*b%R(3) + sqrt_rhor_star*b%L(3) + sqrt_rhol_star*sqrt_rhor_star*(wr_star &
565 & - wl_star)*sign_bx)/denom_ds
566
567 e_doublel = e_starl - sqrt_rhol_star*((vl_star*b%L(2) + wl_star*b%L(3)) - (v_double*by_double &
568 & + w_double*bz_double))*sign_bx
569 e_doubler = e_starr + sqrt_rhor_star*((vr_star*b%R(2) + wr_star*b%R(3)) - (v_double*by_double &
570 & + w_double*bz_double))*sign_bx
571 e_double = 0.5_wp*(e_doublel + e_doubler)
572
573 u_doublel = [rhol_star, rhol_star*s_m, rhol_star*v_double, rhol_star*w_double, by_double, bz_double, &
574 & e_double]
575 u_doubler = [rhor_star, rhor_star*s_m, rhor_star*v_double, rhor_star*w_double, by_double, bz_double, &
576 & e_double]
577
578 ! Select HLLD flux region
579 if (0.0_wp <= s_l) then
580 f_hlld = f_l
581 else if (0.0_wp <= s_starl) then
582 f_hlld = f_l + s_l*(u_starl - u_l)
583 else if (0.0_wp <= s_m) then
584 f_hlld = f_starl + s_starl*(u_doublel - u_starl)
585 else if (0.0_wp <= s_starr) then
586 f_hlld = f_starr + s_starr*(u_doubler - u_starr)
587 else if (0.0_wp <= s_r) then
588 f_hlld = f_r + s_r*(u_starr - u_r)
589 else
590 f_hlld = f_r
591 end if
592
593 ! (12) Write HLLD flux to output arrays
594 flux_rsx_vf(j, k, l, 1) = f_hlld(1) ! TODO multi-component
595 ! Momentum
596 flux_rsx_vf(j, k, l, eqn_idx%cont%end + dir_idx(1)) = f_hlld(2)
597 flux_rsx_vf(j, k, l, eqn_idx%cont%end + dir_idx(2)) = f_hlld(3)
598 flux_rsx_vf(j, k, l, eqn_idx%cont%end + dir_idx(3)) = f_hlld(4)
599 ! Magnetic field
600 if (n == 0) then
601 flux_rsx_vf(j, k, l, eqn_idx%B%beg) = f_hlld(5)
602 flux_rsx_vf(j, k, l, eqn_idx%B%beg + 1) = f_hlld(6)
603 else
604 flux_rsx_vf(j, k, l, eqn_idx%B%beg + dir_idx(1) - 1) = 0._wp
605 flux_rsx_vf(j, k, l, eqn_idx%B%beg + dir_idx(2) - 1) = f_hlld(5)
606 flux_rsx_vf(j, k, l, eqn_idx%B%beg + dir_idx(3) - 1) = f_hlld(6)
607 end if
608 ! Energy
609 flux_rsx_vf(j, k, l, eqn_idx%E) = f_hlld(7)
610 ! Volume fractions
611
612# 257 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
613#if defined(MFC_OpenACC)
614# 257 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
615!$acc loop seq
616# 257 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
617#elif defined(MFC_OpenMP)
618# 257 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
619
620# 257 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
621#endif
622 do i = eqn_idx%adv%beg, eqn_idx%adv%end
623 flux_rsx_vf(j, k, l, i) = 0._wp ! TODO multi-component (zero for now)
624 end do
625
626 flux_src_rsx_vf(j, k, l, eqn_idx%adv%beg) = 0._wp
627 end do
628 end do
629 end do
630
631# 266 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
632#if defined(MFC_OpenACC)
633# 266 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
634!$acc end parallel loop
635# 266 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
636#elif defined(MFC_OpenMP)
637# 266 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
638
639# 266 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
640!$omp end target teams loop
641# 266 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
642#endif
643 end if
644# 73 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
645# 74 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
646# 75 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
647 if (norm_dir == 2) then
648
649# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
650
651# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
652#if defined(MFC_OpenACC)
653# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
654!$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, &
655# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
656!$acc& U_R, U_starL, U_starR, U_doubleL, U_doubleR, F_L, F_R, F_starL, F_starR, F_hlld, s_L, s_R, s_M, s_starL, s_starR, pTot_L, pTot_R, p_star, rhoL_star, rhoR_star, E_starL, E_starR, sqrt_rhoL_star, &
657# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
658!$acc& sqrt_rhoR_star, denom_ds, sign_Bx, vL_star, vR_star, wL_star, wR_star, v_double, w_double, By_double, Bz_double, E_doubleL, E_doubleR, E_double) copyin(norm_dir)
659# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
660#elif defined(MFC_OpenMP)
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
667# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
668!$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, &
669# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
670!$omp& alpha_rho_R, vel, alpha_L, alpha_R, rho, pres, E, H_no_mag, gamma, pi_inf, qv, vel_rms, B, c, c_fast, pres_mag, U_L, U_R, U_starL, U_starR, U_doubleL, U_doubleR, F_L, F_R, F_starL, F_starR, &
671# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
672!$omp& F_hlld, s_L, s_R, s_M, s_starL, s_starR, pTot_L, pTot_R, p_star, rhoL_star, rhoR_star, E_starL, E_starR, sqrt_rhoL_star, sqrt_rhoR_star, denom_ds, sign_Bx, vL_star, vR_star, wL_star, wR_star, &
673# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
674!$omp& v_double, w_double, By_double, Bz_double, E_doubleL, E_doubleR, E_double) map(to:norm_dir)
675# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
676#endif
677# 82 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
678 do l = is3%beg, is3%end
679 do k = is1%beg, is1%end
680 do j = is2%beg, is2%end
681 ! (1) Extract the left/right primitive states
682 do i = 1, eqn_idx%cont%end
683 alpha_rho_l(i) = ql_prim_rsx_vf(j, k, l, i)
684 alpha_rho_r(i) = qr_prim_rsx_vf(j, k + 1, l, i)
685 end do
686
687 ! NOTE: unlike HLL & HLLC, vel_L here is permutated by dir_idx for simpler logic
688 do i = 1, num_vels
689 vel%L(i) = ql_prim_rsx_vf(j, k, l, eqn_idx%cont%end + dir_idx(i))
690 vel%R(i) = qr_prim_rsx_vf(j, k + 1, l, eqn_idx%cont%end + dir_idx(i))
691 end do
692
693 vel_rms%L = sum(vel%L**2._wp)
694 vel_rms%R = sum(vel%R**2._wp)
695
696 do i = 1, num_fluids
697 alpha_l(i) = ql_prim_rsx_vf(j, k, l, eqn_idx%E + i)
698 alpha_r(i) = qr_prim_rsx_vf(j, k + 1, l, eqn_idx%E + i)
699 end do
700
701 pres%L = ql_prim_rsx_vf(j, k, l, eqn_idx%E)
702 pres%R = qr_prim_rsx_vf(j, k + 1, l, eqn_idx%E)
703
704 ! NOTE: unlike HLL, Bx, By, Bz are permutated by dir_idx for simpler logic
705 if (mhd) then
706 if (n == 0) then ! 1D: constant Bx; By, Bz as variables; only in x so not permutated
707 b%L = [bx0, ql_prim_rsx_vf(j, k, l, eqn_idx%B%beg), ql_prim_rsx_vf(j, k, l, &
708 & eqn_idx%B%beg + 1)]
709 b%R = [bx0, qr_prim_rsx_vf(j, k + 1, l, eqn_idx%B%beg), qr_prim_rsx_vf(j, k + 1, l, &
710 & eqn_idx%B%beg + 1)]
711 else ! 2D/3D: Bx, By, Bz as variables
712 b%L = [ql_prim_rsx_vf(j, k, l, eqn_idx%B%beg + dir_idx(1) - 1), ql_prim_rsx_vf(j, k, l, &
713 & eqn_idx%B%beg + dir_idx(2) - 1), ql_prim_rsx_vf(j, k, l, &
714 & eqn_idx%B%beg + dir_idx(3) - 1)]
715 b%R = [qr_prim_rsx_vf(j, k + 1, l, eqn_idx%B%beg + dir_idx(1) - 1), &
716 & qr_prim_rsx_vf(j, k + 1, l, eqn_idx%B%beg + dir_idx(2) - 1), &
717 & qr_prim_rsx_vf(j, k + 1, l, eqn_idx%B%beg + dir_idx(3) - 1)]
718 end if
719 end if
720
721 ! Sum properties of all fluid components
722 rho%L = 0._wp; gamma%L = 0._wp; pi_inf%L = 0._wp; qv%L = 0._wp
723 rho%R = 0._wp; gamma%R = 0._wp; pi_inf%R = 0._wp; qv%R = 0._wp
724
725# 128 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
726#if defined(MFC_OpenACC)
727# 128 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
728!$acc loop seq
729# 128 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
730#elif defined(MFC_OpenMP)
731# 128 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
732
733# 128 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
734#endif
735 do i = 1, num_fluids
736 rho%L = rho%L + alpha_rho_l(i)
737 gamma%L = gamma%L + alpha_l(i)*gammas(i)
738 pi_inf%L = pi_inf%L + alpha_l(i)*pi_infs(i)
739 qv%L = qv%L + alpha_rho_l(i)*qvs(i)
740
741 rho%R = rho%R + alpha_rho_r(i)
742 gamma%R = gamma%R + alpha_r(i)*gammas(i)
743 pi_inf%R = pi_inf%R + alpha_r(i)*pi_infs(i)
744 qv%R = qv%R + alpha_rho_r(i)*qvs(i)
745 end do
746
747 pres_mag%L = 0.5_wp*sum(b%L**2._wp)
748 pres_mag%R = 0.5_wp*sum(b%R**2._wp)
749 e%L = gamma%L*pres%L + pi_inf%L + 0.5_wp*rho%L*vel_rms%L + qv%L + pres_mag%L
750 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
751 h_no_mag%L = (e%L + pres%L - pres_mag%L)/rho%L
752 ! stagnation enthalpy here excludes magnetic energy (only used to find speed of sound)
753 h_no_mag%R = (e%R + pres%R - pres_mag%R)/rho%R
754
755 ! (2) Compute fast wave speeds
756 call s_compute_speed_of_sound(pres%L, rho%L, gamma%L, pi_inf%L, h_no_mag%L, alpha_l, vel_rms%L, &
757 & 0._wp, c%L, qv%L)
758 call s_compute_speed_of_sound(pres%R, rho%R, gamma%R, pi_inf%R, h_no_mag%R, alpha_r, vel_rms%R, &
759 & 0._wp, c%R, qv%R)
760 call s_compute_fast_magnetosonic_speed(rho%L, c%L, b%L, norm_dir, c_fast%L, h_no_mag%L)
761 call s_compute_fast_magnetosonic_speed(rho%R, c%R, b%R, norm_dir, c_fast%R, h_no_mag%R)
762
763 ! (3) Compute contact speed s_M [Miyoshi Equ. (38)]
764 s_l = min(vel%L(1) - c_fast%L, vel%R(1) - c_fast%R)
765 s_r = max(vel%R(1) + c_fast%R, vel%L(1) + c_fast%L)
766
767 ptot_l = pres%L + pres_mag%L
768 ptot_r = pres%R + pres_mag%R
769
770 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 &
771 & - vel%R(1))*rho%R - (s_l - vel%L(1))*rho%L))
772
773 ! (4) Compute star state variables
774 rhol_star = rho%L*(s_l - vel%L(1))/(s_l - s_m)
775 rhor_star = rho%R*(s_r - vel%R(1))/(s_r - s_m)
776 p_star = ptot_l + rho%L*(s_l - vel%L(1))*(s_m - vel%L(1))/(s_l - s_m)
777 e_starl = ((s_l - vel%L(1))*e%L - ptot_l*vel%L(1) + p_star*s_m)/(s_l - s_m)
778 e_starr = ((s_r - vel%R(1))*e%R - ptot_r*vel%R(1) + p_star*s_m)/(s_r - s_m)
779
780 ! (5) Compute left/right state vectors and fluxes
781 u_l = [rho%L, rho%L*vel%L(1:3), b%L(2:3), e%L]
782 u_starl = [rhol_star, rhol_star*s_m, rhol_star*vel%L(2:3), b%L(2:3), e_starl]
783 u_r = [rho%R, rho%R*vel%R(1:3), b%R(2:3), e%R]
784 u_starr = [rhor_star, rhor_star*s_m, rhor_star*vel%R(2:3), b%R(2:3), e_starr]
785
786 ! Compute the left/right fluxes
787 f_l(1) = u_l(2)
788 f_l(2) = u_l(2)*vel%L(1) - b%L(1)*b%L(1) + ptot_l
789 f_l(3:4) = u_l(2)*vel%L(2:3) - b%L(1)*b%L(2:3)
790 f_l(5:6) = vel%L(1)*b%L(2:3) - vel%L(2:3)*b%L(1)
791 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))
792
793 f_r(1) = u_r(2)
794 f_r(2) = u_r(2)*vel%R(1) - b%R(1)*b%R(1) + ptot_r
795 f_r(3:4) = u_r(2)*vel%R(2:3) - b%R(1)*b%R(2:3)
796 f_r(5:6) = vel%R(1)*b%R(2:3) - vel%R(2:3)*b%R(1)
797 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))
798 ! HLLD star-state fluxes via HLL jump relation
799 f_starl = f_l + s_l*(u_starl - u_l)
800 f_starr = f_r + s_r*(u_starr - u_r)
801 ! Alfven wave speeds bounding the rotational discontinuities
802 s_starl = s_m - abs(b%L(1))/sqrt(rhol_star)
803 s_starr = s_m + abs(b%L(1))/sqrt(rhor_star)
804 ! HLLD double-star (intermediate) states across rotational discontinuities
805 sqrt_rhol_star = sqrt(rhol_star); sqrt_rhor_star = sqrt(rhor_star)
806 vl_star = vel%L(2); wl_star = vel%L(3)
807 vr_star = vel%R(2); wr_star = vel%R(3)
808
809 ! (6) Compute the double-star states [Miyoshi Eqns. (59)-(62)]
810 denom_ds = sqrt_rhol_star + sqrt_rhor_star
811 sign_bx = sign(1._wp, b%L(1))
812 v_double = (sqrt_rhol_star*vl_star + sqrt_rhor_star*vr_star + (b%R(2) - b%L(2))*sign_bx)/denom_ds
813 w_double = (sqrt_rhol_star*wl_star + sqrt_rhor_star*wr_star + (b%R(3) - b%L(3))*sign_bx)/denom_ds
814 by_double = (sqrt_rhol_star*b%R(2) + sqrt_rhor_star*b%L(2) + sqrt_rhol_star*sqrt_rhor_star*(vr_star &
815 & - vl_star)*sign_bx)/denom_ds
816 bz_double = (sqrt_rhol_star*b%R(3) + sqrt_rhor_star*b%L(3) + sqrt_rhol_star*sqrt_rhor_star*(wr_star &
817 & - wl_star)*sign_bx)/denom_ds
818
819 e_doublel = e_starl - sqrt_rhol_star*((vl_star*b%L(2) + wl_star*b%L(3)) - (v_double*by_double &
820 & + w_double*bz_double))*sign_bx
821 e_doubler = e_starr + sqrt_rhor_star*((vr_star*b%R(2) + wr_star*b%R(3)) - (v_double*by_double &
822 & + w_double*bz_double))*sign_bx
823 e_double = 0.5_wp*(e_doublel + e_doubler)
824
825 u_doublel = [rhol_star, rhol_star*s_m, rhol_star*v_double, rhol_star*w_double, by_double, bz_double, &
826 & e_double]
827 u_doubler = [rhor_star, rhor_star*s_m, rhor_star*v_double, rhor_star*w_double, by_double, bz_double, &
828 & e_double]
829
830 ! Select HLLD flux region
831 if (0.0_wp <= s_l) then
832 f_hlld = f_l
833 else if (0.0_wp <= s_starl) then
834 f_hlld = f_l + s_l*(u_starl - u_l)
835 else if (0.0_wp <= s_m) then
836 f_hlld = f_starl + s_starl*(u_doublel - u_starl)
837 else if (0.0_wp <= s_starr) then
838 f_hlld = f_starr + s_starr*(u_doubler - u_starr)
839 else if (0.0_wp <= s_r) then
840 f_hlld = f_r + s_r*(u_starr - u_r)
841 else
842 f_hlld = f_r
843 end if
844
845 ! (12) Write HLLD flux to output arrays
846 flux_rsx_vf(j, k, l, 1) = f_hlld(1) ! TODO multi-component
847 ! Momentum
848 flux_rsx_vf(j, k, l, eqn_idx%cont%end + dir_idx(1)) = f_hlld(2)
849 flux_rsx_vf(j, k, l, eqn_idx%cont%end + dir_idx(2)) = f_hlld(3)
850 flux_rsx_vf(j, k, l, eqn_idx%cont%end + dir_idx(3)) = f_hlld(4)
851 ! Magnetic field
852 if (n == 0) then
853 flux_rsx_vf(j, k, l, eqn_idx%B%beg) = f_hlld(5)
854 flux_rsx_vf(j, k, l, eqn_idx%B%beg + 1) = f_hlld(6)
855 else
856 flux_rsx_vf(j, k, l, eqn_idx%B%beg + dir_idx(1) - 1) = 0._wp
857 flux_rsx_vf(j, k, l, eqn_idx%B%beg + dir_idx(2) - 1) = f_hlld(5)
858 flux_rsx_vf(j, k, l, eqn_idx%B%beg + dir_idx(3) - 1) = f_hlld(6)
859 end if
860 ! Energy
861 flux_rsx_vf(j, k, l, eqn_idx%E) = f_hlld(7)
862 ! Volume fractions
863
864# 257 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
865#if defined(MFC_OpenACC)
866# 257 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
867!$acc loop seq
868# 257 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
869#elif defined(MFC_OpenMP)
870# 257 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
871
872# 257 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
873#endif
874 do i = eqn_idx%adv%beg, eqn_idx%adv%end
875 flux_rsx_vf(j, k, l, i) = 0._wp ! TODO multi-component (zero for now)
876 end do
877
878 flux_src_rsx_vf(j, k, l, eqn_idx%adv%beg) = 0._wp
879 end do
880 end do
881 end do
882
883# 266 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
884#if defined(MFC_OpenACC)
885# 266 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
886!$acc end parallel loop
887# 266 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
888#elif defined(MFC_OpenMP)
889# 266 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
890
891# 266 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
892!$omp end target teams loop
893# 266 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
894#endif
895 end if
896# 73 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
897# 74 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
898# 75 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
899 if (norm_dir == 3) then
900
901# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
902
903# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
904#if defined(MFC_OpenACC)
905# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
906!$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, &
907# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
908!$acc& U_R, U_starL, U_starR, U_doubleL, U_doubleR, F_L, F_R, F_starL, F_starR, F_hlld, s_L, s_R, s_M, s_starL, s_starR, pTot_L, pTot_R, p_star, rhoL_star, rhoR_star, E_starL, E_starR, sqrt_rhoL_star, &
909# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
910!$acc& sqrt_rhoR_star, denom_ds, sign_Bx, vL_star, vR_star, wL_star, wR_star, v_double, w_double, By_double, Bz_double, E_doubleL, E_doubleR, E_double) copyin(norm_dir)
911# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
912#elif defined(MFC_OpenMP)
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
917# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
918
919# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
920!$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, &
921# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
922!$omp& alpha_rho_R, vel, alpha_L, alpha_R, rho, pres, E, H_no_mag, gamma, pi_inf, qv, vel_rms, B, c, c_fast, pres_mag, U_L, U_R, U_starL, U_starR, U_doubleL, U_doubleR, F_L, F_R, F_starL, F_starR, &
923# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
924!$omp& F_hlld, s_L, s_R, s_M, s_starL, s_starR, pTot_L, pTot_R, p_star, rhoL_star, rhoR_star, E_starL, E_starR, sqrt_rhoL_star, sqrt_rhoR_star, denom_ds, sign_Bx, vL_star, vR_star, wL_star, wR_star, &
925# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
926!$omp& v_double, w_double, By_double, Bz_double, E_doubleL, E_doubleR, E_double) map(to:norm_dir)
927# 76 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
928#endif
929# 82 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
930 do l = is1%beg, is1%end
931 do k = is2%beg, is2%end
932 do j = is3%beg, is3%end
933 ! (1) Extract the left/right primitive states
934 do i = 1, eqn_idx%cont%end
935 alpha_rho_l(i) = ql_prim_rsx_vf(j, k, l, i)
936 alpha_rho_r(i) = qr_prim_rsx_vf(j, k, l + 1, i)
937 end do
938
939 ! NOTE: unlike HLL & HLLC, vel_L here is permutated by dir_idx for simpler logic
940 do i = 1, num_vels
941 vel%L(i) = ql_prim_rsx_vf(j, k, l, eqn_idx%cont%end + dir_idx(i))
942 vel%R(i) = qr_prim_rsx_vf(j, k, l + 1, eqn_idx%cont%end + dir_idx(i))
943 end do
944
945 vel_rms%L = sum(vel%L**2._wp)
946 vel_rms%R = sum(vel%R**2._wp)
947
948 do i = 1, num_fluids
949 alpha_l(i) = ql_prim_rsx_vf(j, k, l, eqn_idx%E + i)
950 alpha_r(i) = qr_prim_rsx_vf(j, k, l + 1, eqn_idx%E + i)
951 end do
952
953 pres%L = ql_prim_rsx_vf(j, k, l, eqn_idx%E)
954 pres%R = qr_prim_rsx_vf(j, k, l + 1, eqn_idx%E)
955
956 ! NOTE: unlike HLL, Bx, By, Bz are permutated by dir_idx for simpler logic
957 if (mhd) then
958 if (n == 0) then ! 1D: constant Bx; By, Bz as variables; only in x so not permutated
959 b%L = [bx0, ql_prim_rsx_vf(j, k, l, eqn_idx%B%beg), ql_prim_rsx_vf(j, k, l, &
960 & eqn_idx%B%beg + 1)]
961 b%R = [bx0, qr_prim_rsx_vf(j, k, l + 1, eqn_idx%B%beg), qr_prim_rsx_vf(j, k, l + 1, &
962 & eqn_idx%B%beg + 1)]
963 else ! 2D/3D: Bx, By, Bz as variables
964 b%L = [ql_prim_rsx_vf(j, k, l, eqn_idx%B%beg + dir_idx(1) - 1), ql_prim_rsx_vf(j, k, l, &
965 & eqn_idx%B%beg + dir_idx(2) - 1), ql_prim_rsx_vf(j, k, l, &
966 & eqn_idx%B%beg + dir_idx(3) - 1)]
967 b%R = [qr_prim_rsx_vf(j, k, l + 1, eqn_idx%B%beg + dir_idx(1) - 1), &
968 & qr_prim_rsx_vf(j, k, l + 1, eqn_idx%B%beg + dir_idx(2) - 1), &
969 & qr_prim_rsx_vf(j, k, l + 1, eqn_idx%B%beg + dir_idx(3) - 1)]
970 end if
971 end if
972
973 ! Sum properties of all fluid components
974 rho%L = 0._wp; gamma%L = 0._wp; pi_inf%L = 0._wp; qv%L = 0._wp
975 rho%R = 0._wp; gamma%R = 0._wp; pi_inf%R = 0._wp; qv%R = 0._wp
976
977# 128 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
978#if defined(MFC_OpenACC)
979# 128 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
980!$acc loop seq
981# 128 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
982#elif defined(MFC_OpenMP)
983# 128 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
984
985# 128 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
986#endif
987 do i = 1, num_fluids
988 rho%L = rho%L + alpha_rho_l(i)
989 gamma%L = gamma%L + alpha_l(i)*gammas(i)
990 pi_inf%L = pi_inf%L + alpha_l(i)*pi_infs(i)
991 qv%L = qv%L + alpha_rho_l(i)*qvs(i)
992
993 rho%R = rho%R + alpha_rho_r(i)
994 gamma%R = gamma%R + alpha_r(i)*gammas(i)
995 pi_inf%R = pi_inf%R + alpha_r(i)*pi_infs(i)
996 qv%R = qv%R + alpha_rho_r(i)*qvs(i)
997 end do
998
999 pres_mag%L = 0.5_wp*sum(b%L**2._wp)
1000 pres_mag%R = 0.5_wp*sum(b%R**2._wp)
1001 e%L = gamma%L*pres%L + pi_inf%L + 0.5_wp*rho%L*vel_rms%L + qv%L + pres_mag%L
1002 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
1003 h_no_mag%L = (e%L + pres%L - pres_mag%L)/rho%L
1004 ! stagnation enthalpy here excludes magnetic energy (only used to find speed of sound)
1005 h_no_mag%R = (e%R + pres%R - pres_mag%R)/rho%R
1006
1007 ! (2) Compute fast wave speeds
1008 call s_compute_speed_of_sound(pres%L, rho%L, gamma%L, pi_inf%L, h_no_mag%L, alpha_l, vel_rms%L, &
1009 & 0._wp, c%L, qv%L)
1010 call s_compute_speed_of_sound(pres%R, rho%R, gamma%R, pi_inf%R, h_no_mag%R, alpha_r, vel_rms%R, &
1011 & 0._wp, c%R, qv%R)
1012 call s_compute_fast_magnetosonic_speed(rho%L, c%L, b%L, norm_dir, c_fast%L, h_no_mag%L)
1013 call s_compute_fast_magnetosonic_speed(rho%R, c%R, b%R, norm_dir, c_fast%R, h_no_mag%R)
1014
1015 ! (3) Compute contact speed s_M [Miyoshi Equ. (38)]
1016 s_l = min(vel%L(1) - c_fast%L, vel%R(1) - c_fast%R)
1017 s_r = max(vel%R(1) + c_fast%R, vel%L(1) + c_fast%L)
1018
1019 ptot_l = pres%L + pres_mag%L
1020 ptot_r = pres%R + pres_mag%R
1021
1022 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 &
1023 & - vel%R(1))*rho%R - (s_l - vel%L(1))*rho%L))
1024
1025 ! (4) Compute star state variables
1026 rhol_star = rho%L*(s_l - vel%L(1))/(s_l - s_m)
1027 rhor_star = rho%R*(s_r - vel%R(1))/(s_r - s_m)
1028 p_star = ptot_l + rho%L*(s_l - vel%L(1))*(s_m - vel%L(1))/(s_l - s_m)
1029 e_starl = ((s_l - vel%L(1))*e%L - ptot_l*vel%L(1) + p_star*s_m)/(s_l - s_m)
1030 e_starr = ((s_r - vel%R(1))*e%R - ptot_r*vel%R(1) + p_star*s_m)/(s_r - s_m)
1031
1032 ! (5) Compute left/right state vectors and fluxes
1033 u_l = [rho%L, rho%L*vel%L(1:3), b%L(2:3), e%L]
1034 u_starl = [rhol_star, rhol_star*s_m, rhol_star*vel%L(2:3), b%L(2:3), e_starl]
1035 u_r = [rho%R, rho%R*vel%R(1:3), b%R(2:3), e%R]
1036 u_starr = [rhor_star, rhor_star*s_m, rhor_star*vel%R(2:3), b%R(2:3), e_starr]
1037
1038 ! Compute the left/right fluxes
1039 f_l(1) = u_l(2)
1040 f_l(2) = u_l(2)*vel%L(1) - b%L(1)*b%L(1) + ptot_l
1041 f_l(3:4) = u_l(2)*vel%L(2:3) - b%L(1)*b%L(2:3)
1042 f_l(5:6) = vel%L(1)*b%L(2:3) - vel%L(2:3)*b%L(1)
1043 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))
1044
1045 f_r(1) = u_r(2)
1046 f_r(2) = u_r(2)*vel%R(1) - b%R(1)*b%R(1) + ptot_r
1047 f_r(3:4) = u_r(2)*vel%R(2:3) - b%R(1)*b%R(2:3)
1048 f_r(5:6) = vel%R(1)*b%R(2:3) - vel%R(2:3)*b%R(1)
1049 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))
1050 ! HLLD star-state fluxes via HLL jump relation
1051 f_starl = f_l + s_l*(u_starl - u_l)
1052 f_starr = f_r + s_r*(u_starr - u_r)
1053 ! Alfven wave speeds bounding the rotational discontinuities
1054 s_starl = s_m - abs(b%L(1))/sqrt(rhol_star)
1055 s_starr = s_m + abs(b%L(1))/sqrt(rhor_star)
1056 ! HLLD double-star (intermediate) states across rotational discontinuities
1057 sqrt_rhol_star = sqrt(rhol_star); sqrt_rhor_star = sqrt(rhor_star)
1058 vl_star = vel%L(2); wl_star = vel%L(3)
1059 vr_star = vel%R(2); wr_star = vel%R(3)
1060
1061 ! (6) Compute the double-star states [Miyoshi Eqns. (59)-(62)]
1062 denom_ds = sqrt_rhol_star + sqrt_rhor_star
1063 sign_bx = sign(1._wp, b%L(1))
1064 v_double = (sqrt_rhol_star*vl_star + sqrt_rhor_star*vr_star + (b%R(2) - b%L(2))*sign_bx)/denom_ds
1065 w_double = (sqrt_rhol_star*wl_star + sqrt_rhor_star*wr_star + (b%R(3) - b%L(3))*sign_bx)/denom_ds
1066 by_double = (sqrt_rhol_star*b%R(2) + sqrt_rhor_star*b%L(2) + sqrt_rhol_star*sqrt_rhor_star*(vr_star &
1067 & - vl_star)*sign_bx)/denom_ds
1068 bz_double = (sqrt_rhol_star*b%R(3) + sqrt_rhor_star*b%L(3) + sqrt_rhol_star*sqrt_rhor_star*(wr_star &
1069 & - wl_star)*sign_bx)/denom_ds
1070
1071 e_doublel = e_starl - sqrt_rhol_star*((vl_star*b%L(2) + wl_star*b%L(3)) - (v_double*by_double &
1072 & + w_double*bz_double))*sign_bx
1073 e_doubler = e_starr + sqrt_rhor_star*((vr_star*b%R(2) + wr_star*b%R(3)) - (v_double*by_double &
1074 & + w_double*bz_double))*sign_bx
1075 e_double = 0.5_wp*(e_doublel + e_doubler)
1076
1077 u_doublel = [rhol_star, rhol_star*s_m, rhol_star*v_double, rhol_star*w_double, by_double, bz_double, &
1078 & e_double]
1079 u_doubler = [rhor_star, rhor_star*s_m, rhor_star*v_double, rhor_star*w_double, by_double, bz_double, &
1080 & e_double]
1081
1082 ! Select HLLD flux region
1083 if (0.0_wp <= s_l) then
1084 f_hlld = f_l
1085 else if (0.0_wp <= s_starl) then
1086 f_hlld = f_l + s_l*(u_starl - u_l)
1087 else if (0.0_wp <= s_m) then
1088 f_hlld = f_starl + s_starl*(u_doublel - u_starl)
1089 else if (0.0_wp <= s_starr) then
1090 f_hlld = f_starr + s_starr*(u_doubler - u_starr)
1091 else if (0.0_wp <= s_r) then
1092 f_hlld = f_r + s_r*(u_starr - u_r)
1093 else
1094 f_hlld = f_r
1095 end if
1096
1097 ! (12) Write HLLD flux to output arrays
1098 flux_rsx_vf(j, k, l, 1) = f_hlld(1) ! TODO multi-component
1099 ! Momentum
1100 flux_rsx_vf(j, k, l, eqn_idx%cont%end + dir_idx(1)) = f_hlld(2)
1101 flux_rsx_vf(j, k, l, eqn_idx%cont%end + dir_idx(2)) = f_hlld(3)
1102 flux_rsx_vf(j, k, l, eqn_idx%cont%end + dir_idx(3)) = f_hlld(4)
1103 ! Magnetic field
1104 if (n == 0) then
1105 flux_rsx_vf(j, k, l, eqn_idx%B%beg) = f_hlld(5)
1106 flux_rsx_vf(j, k, l, eqn_idx%B%beg + 1) = f_hlld(6)
1107 else
1108 flux_rsx_vf(j, k, l, eqn_idx%B%beg + dir_idx(1) - 1) = 0._wp
1109 flux_rsx_vf(j, k, l, eqn_idx%B%beg + dir_idx(2) - 1) = f_hlld(5)
1110 flux_rsx_vf(j, k, l, eqn_idx%B%beg + dir_idx(3) - 1) = f_hlld(6)
1111 end if
1112 ! Energy
1113 flux_rsx_vf(j, k, l, eqn_idx%E) = f_hlld(7)
1114 ! Volume fractions
1115
1116# 257 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
1117#if defined(MFC_OpenACC)
1118# 257 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
1119!$acc loop seq
1120# 257 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
1121#elif defined(MFC_OpenMP)
1122# 257 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
1123
1124# 257 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
1125#endif
1126 do i = eqn_idx%adv%beg, eqn_idx%adv%end
1127 flux_rsx_vf(j, k, l, i) = 0._wp ! TODO multi-component (zero for now)
1128 end do
1129
1130 flux_src_rsx_vf(j, k, l, eqn_idx%adv%beg) = 0._wp
1131 end do
1132 end do
1133 end do
1134
1135# 266 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
1136#if defined(MFC_OpenACC)
1137# 266 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
1138!$acc end parallel loop
1139# 266 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
1140#elif defined(MFC_OpenMP)
1141# 266 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
1142
1143# 266 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
1144!$omp end target teams loop
1145# 266 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
1146#endif
1147 end if
1148# 269 "/home/runner/work/MFC/MFC/src/simulation/m_riemann_solver_hlld.fpp"
1149
1150 call s_finalize_riemann_solver(flux_vf, flux_src_vf, flux_gsrc_vf, norm_dir)
1151
1152 end subroutine s_hlld_riemann_solver
1153
1154end 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:744
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:477
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:1630
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:548
m_riemann_state::is2
type(int_bounds_info) is2
Definition m_riemann_state.fpp.f90:477
m_riemann_state::s_finalize_riemann_solver
subroutine s_finalize_riemann_solver(flux_vf, flux_src_vf, flux_gsrc_vf, norm_dir)
Reshape and copy the Riemann-solver flux buffers back to the physical-space output arrays for the sel...
Definition m_riemann_state.fpp.f90:2923
m_riemann_state::is1
type(int_bounds_info) is1
Definition m_riemann_state.fpp.f90:477
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_speed_of_sound
subroutine, public 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:2532
m_variables_conversion::s_compute_fast_magnetosonic_speed
subroutine, public s_compute_fast_magnetosonic_speed(rho, c, b, norm, c_fast, h)
Compute the fast magnetosonic wave speed from the sound speed, density, and magnetic field components...
Definition m_variables_conversion.fpp.f90:2614
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