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