Skip to content

Commit 83a0315

Browse files
sbryngelsonsbryngelson
committed
src: split s_hlld_riemann_solver into m_riemann_solver_hlld
NFC, pure motion (calibration commit for the solver extractions). The solver body accesses flux_rsx_vf, flux_src_rsx_vf, is1/2/3 by host association and calls the per-sweep lifecycle helpers, both now use-associated from m_riemann_state. Core re-exports s_hlld_riemann_solver through its existing public list. Statement-multiset union vs the pre-split file differs only by module boilerplate; GPU-directive census md5 unchanged; declare-target scoping verified for all files.
1 parent b1243d9 commit 83a0315

3 files changed

Lines changed: 277 additions & 255 deletions

File tree

‎docs/module_categories.json‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,7 @@
77
"m_weno",
88
"m_riemann_solvers",
99
"m_riemann_state",
10+
"m_riemann_solver_hlld",
1011
"m_muscl",
1112
"m_variables_conversion",
1213
"m_thinc"

‎src/simulation/m_riemann_solver_hlld.fpp‎

Lines changed: 275 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,275 @@
1+
!>
2+
!! @file
3+
!! @brief Contains module m_riemann_solver_hlld
4+
5+
!> @brief HLLD approximate Riemann solver for MHD, Miyoshi & Kusano JCP (2005)
6+
#:include 'case.fpp'
7+
#:include 'macros.fpp'
8+
9+
module m_riemann_solver_hlld
10+
11+
use m_derived_types
12+
use m_global_parameters
13+
use m_variables_conversion
14+
use m_riemann_state
15+
16+
implicit none
17+
18+
contains
19+
20+
!> HLLD Riemann solver for MHD, Miyoshi & Kusano JCP (2005)
21+
subroutine s_hlld_riemann_solver(qL_prim_rsx_vf, dqL_prim_dx_vf, dqL_prim_dy_vf, &
22+
23+
& 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, &
24+
& flux_vf, flux_src_vf, flux_gsrc_vf, norm_dir, ix, iy, iz)
25+
26+
real(wp), dimension(idwbuff(1)%beg:,idwbuff(2)%beg:,idwbuff(3)%beg:,1:), intent(inout) :: qL_prim_rsx_vf, qR_prim_rsx_vf
27+
type(scalar_field), allocatable, dimension(:), intent(inout) :: dqL_prim_dx_vf, dqR_prim_dx_vf, dqL_prim_dy_vf, &
28+
& dqR_prim_dy_vf, dqL_prim_dz_vf, dqR_prim_dz_vf
29+
30+
type(scalar_field), allocatable, dimension(:), intent(inout) :: qL_prim_vf, qR_prim_vf
31+
type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf
32+
type(scalar_field), dimension(sys_size), intent(inout) :: flux_vf, flux_src_vf, flux_gsrc_vf
33+
integer, intent(in) :: norm_dir
34+
type(int_bounds_info), intent(in) :: ix, iy, iz
35+
36+
! Local variables:
37+
38+
#:if not MFC_CASE_OPTIMIZATION and USING_AMD
39+
real(wp), dimension(3) :: alpha_L, alpha_R, alpha_rho_L, alpha_rho_R
40+
#:else
41+
real(wp), dimension(num_fluids) :: alpha_L, alpha_R, alpha_rho_L, alpha_rho_R
42+
#:endif
43+
type(riemann_states_vec3) :: vel
44+
type(riemann_states) :: rho, pres, E, H_no_mag
45+
type(riemann_states) :: gamma, pi_inf, qv
46+
type(riemann_states) :: vel_rms
47+
type(riemann_states_vec3) :: B
48+
type(riemann_states) :: c, c_fast, pres_mag
49+
50+
! HLLD speeds and intermediate state variables:
51+
real(wp) :: s_L, s_R, s_M, s_starL, s_starR
52+
real(wp) :: pTot_L, pTot_R, p_star, rhoL_star, rhoR_star, E_starL, E_starR
53+
real(wp), dimension(7) :: U_L, U_R, U_starL, U_starR, U_doubleL, U_doubleR
54+
real(wp), dimension(7) :: F_L, F_R, F_starL, F_starR, F_hlld
55+
56+
! 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
57+
! normal velocity, and x is the normal direction Note: Bx is omitted as the magnetic flux is always zero in the normal
58+
! direction
59+
60+
real(wp) :: sqrt_rhoL_star, sqrt_rhoR_star, denom_ds, sign_Bx
61+
real(wp) :: vL_star, vR_star, wL_star, wR_star
62+
real(wp) :: v_double, w_double, By_double, Bz_double, E_doubleL, E_doubleR, E_double
63+
integer :: i, j, k, l
64+
65+
call s_populate_riemann_states_variables_buffers(qL_prim_rsx_vf, dqL_prim_dx_vf, dqL_prim_dy_vf, dqL_prim_dz_vf, &
66+
& qR_prim_rsx_vf, dqR_prim_dx_vf, dqR_prim_dy_vf, dqR_prim_dz_vf, norm_dir, ix, iy, iz)
67+
68+
call s_initialize_riemann_solver(flux_src_vf, norm_dir)
69+
70+
#:for NORM_DIR, XYZ, STENCIL_VAR, COORDS, X_BND, Y_BND, Z_BND in &
71+
[(1, 'x', 'j', '{STENCIL_IDX}, k, l', 'is1', 'is2', 'is3'), &
72+
(2, 'y', 'k', 'j, {STENCIL_IDX}, l', 'is2', 'is1', 'is3'), &
73+
(3, 'z', 'l', 'j, k, {STENCIL_IDX}', 'is3', 'is2', 'is1')]
74+
#:set SV = STENCIL_VAR
75+
#:set SF = lambda offs: COORDS.format(STENCIL_IDX = SV + offs)
76+
if (norm_dir == ${NORM_DIR}$) then
77+
$:GPU_PARALLEL_LOOP(collapse=3, private='[alpha_rho_L, alpha_rho_R, vel, alpha_L, alpha_R, rho, pres, E, &
78+
& H_no_mag, gamma, pi_inf, qv, vel_rms, B, c, c_fast, pres_mag, U_L, U_R, U_starL, U_starR, &
79+
& U_doubleL, U_doubleR, F_L, F_R, F_starL, F_starR, F_hlld, s_L, s_R, s_M, s_starL, s_starR, &
80+
& pTot_L, pTot_R, p_star, rhoL_star, rhoR_star, E_starL, E_starR, sqrt_rhoL_star, &
81+
& sqrt_rhoR_star, denom_ds, sign_Bx, vL_star, vR_star, wL_star, wR_star, v_double, w_double, &
82+
& By_double, Bz_double, E_doubleL, E_doubleR, E_double]', copyin='[norm_dir]')
83+
do l = ${Z_BND}$%beg, ${Z_BND}$%end
84+
do k = ${Y_BND}$%beg, ${Y_BND}$%end
85+
do j = ${X_BND}$%beg, ${X_BND}$%end
86+
! (1) Extract the left/right primitive states
87+
do i = 1, eqn_idx%cont%end
88+
alpha_rho_L(i) = qL_prim_rsx_vf(${SF('')}$, i)
89+
alpha_rho_R(i) = qR_prim_rsx_vf(${SF(' + 1')}$, i)
90+
end do
91+
92+
! NOTE: unlike HLL & HLLC, vel_L here is permutated by dir_idx for simpler logic
93+
do i = 1, num_vels
94+
vel%L(i) = qL_prim_rsx_vf(${SF('')}$, eqn_idx%cont%end + dir_idx(i))
95+
vel%R(i) = qR_prim_rsx_vf(${SF(' + 1')}$, eqn_idx%cont%end + dir_idx(i))
96+
end do
97+
98+
vel_rms%L = sum(vel%L**2._wp)
99+
vel_rms%R = sum(vel%R**2._wp)
100+
101+
do i = 1, num_fluids
102+
alpha_L(i) = qL_prim_rsx_vf(${SF('')}$, eqn_idx%E + i)
103+
alpha_R(i) = qR_prim_rsx_vf(${SF(' + 1')}$, eqn_idx%E + i)
104+
end do
105+
106+
pres%L = qL_prim_rsx_vf(${SF('')}$, eqn_idx%E)
107+
pres%R = qR_prim_rsx_vf(${SF(' + 1')}$, eqn_idx%E)
108+
109+
! NOTE: unlike HLL, Bx, By, Bz are permutated by dir_idx for simpler logic
110+
if (mhd) then
111+
if (n == 0) then ! 1D: constant Bx; By, Bz as variables; only in x so not permutated
112+
B%L = [Bx0, qL_prim_rsx_vf(${SF('')}$, eqn_idx%B%beg), qL_prim_rsx_vf(${SF('')}$, &
113+
& eqn_idx%B%beg + 1)]
114+
B%R = [Bx0, qR_prim_rsx_vf(${SF(' + 1')}$, eqn_idx%B%beg), qR_prim_rsx_vf(${SF(' + 1')}$, &
115+
& eqn_idx%B%beg + 1)]
116+
else ! 2D/3D: Bx, By, Bz as variables
117+
B%L = [qL_prim_rsx_vf(${SF('')}$, eqn_idx%B%beg + dir_idx(1) - 1), qL_prim_rsx_vf(${SF('')}$, &
118+
& eqn_idx%B%beg + dir_idx(2) - 1), qL_prim_rsx_vf(${SF('')}$, &
119+
& eqn_idx%B%beg + dir_idx(3) - 1)]
120+
B%R = [qR_prim_rsx_vf(${SF(' + 1')}$, eqn_idx%B%beg + dir_idx(1) - 1), &
121+
& qR_prim_rsx_vf(${SF(' + 1')}$, eqn_idx%B%beg + dir_idx(2) - 1), &
122+
& qR_prim_rsx_vf(${SF(' + 1')}$, eqn_idx%B%beg + dir_idx(3) - 1)]
123+
end if
124+
end if
125+
126+
! Sum properties of all fluid components
127+
rho%L = 0._wp; gamma%L = 0._wp; pi_inf%L = 0._wp; qv%L = 0._wp
128+
rho%R = 0._wp; gamma%R = 0._wp; pi_inf%R = 0._wp; qv%R = 0._wp
129+
$:GPU_LOOP(parallelism='[seq]')
130+
do i = 1, num_fluids
131+
rho%L = rho%L + alpha_rho_L(i)
132+
gamma%L = gamma%L + alpha_L(i)*gammas(i)
133+
pi_inf%L = pi_inf%L + alpha_L(i)*pi_infs(i)
134+
qv%L = qv%L + alpha_rho_L(i)*qvs(i)
135+
136+
rho%R = rho%R + alpha_rho_R(i)
137+
gamma%R = gamma%R + alpha_R(i)*gammas(i)
138+
pi_inf%R = pi_inf%R + alpha_R(i)*pi_infs(i)
139+
qv%R = qv%R + alpha_rho_R(i)*qvs(i)
140+
end do
141+
142+
pres_mag%L = 0.5_wp*sum(B%L**2._wp)
143+
pres_mag%R = 0.5_wp*sum(B%R**2._wp)
144+
E%L = gamma%L*pres%L + pi_inf%L + 0.5_wp*rho%L*vel_rms%L + qv%L + pres_mag%L
145+
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
146+
H_no_mag%L = (E%L + pres%L - pres_mag%L)/rho%L
147+
! stagnation enthalpy here excludes magnetic energy (only used to find speed of sound)
148+
H_no_mag%R = (E%R + pres%R - pres_mag%R)/rho%R
149+
150+
! (2) Compute fast wave speeds
151+
call s_compute_speed_of_sound(pres%L, rho%L, gamma%L, pi_inf%L, H_no_mag%L, alpha_L, vel_rms%L, &
152+
& 0._wp, c%L, qv%L)
153+
call s_compute_speed_of_sound(pres%R, rho%R, gamma%R, pi_inf%R, H_no_mag%R, alpha_R, vel_rms%R, &
154+
& 0._wp, c%R, qv%R)
155+
call s_compute_fast_magnetosonic_speed(rho%L, c%L, B%L, norm_dir, c_fast%L, H_no_mag%L)
156+
call s_compute_fast_magnetosonic_speed(rho%R, c%R, B%R, norm_dir, c_fast%R, H_no_mag%R)
157+
158+
! (3) Compute contact speed s_M [Miyoshi Equ. (38)]
159+
s_L = min(vel%L(1) - c_fast%L, vel%R(1) - c_fast%R)
160+
s_R = max(vel%R(1) + c_fast%R, vel%L(1) + c_fast%L)
161+
162+
pTot_L = pres%L + pres_mag%L
163+
pTot_R = pres%R + pres_mag%R
164+
165+
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 &
166+
& - vel%R(1))*rho%R - (s_L - vel%L(1))*rho%L))
167+
168+
! (4) Compute star state variables
169+
rhoL_star = rho%L*(s_L - vel%L(1))/(s_L - s_M)
170+
rhoR_star = rho%R*(s_R - vel%R(1))/(s_R - s_M)
171+
p_star = pTot_L + rho%L*(s_L - vel%L(1))*(s_M - vel%L(1))/(s_L - s_M)
172+
E_starL = ((s_L - vel%L(1))*E%L - pTot_L*vel%L(1) + p_star*s_M)/(s_L - s_M)
173+
E_starR = ((s_R - vel%R(1))*E%R - pTot_R*vel%R(1) + p_star*s_M)/(s_R - s_M)
174+
175+
! (5) Compute left/right state vectors and fluxes
176+
U_L = [rho%L, rho%L*vel%L(1:3), B%L(2:3), E%L]
177+
U_starL = [rhoL_star, rhoL_star*s_M, rhoL_star*vel%L(2:3), B%L(2:3), E_starL]
178+
U_R = [rho%R, rho%R*vel%R(1:3), B%R(2:3), E%R]
179+
U_starR = [rhoR_star, rhoR_star*s_M, rhoR_star*vel%R(2:3), B%R(2:3), E_starR]
180+
181+
! Compute the left/right fluxes
182+
F_L(1) = U_L(2)
183+
F_L(2) = U_L(2)*vel%L(1) - B%L(1)*B%L(1) + pTot_L
184+
F_L(3:4) = U_L(2)*vel%L(2:3) - B%L(1)*B%L(2:3)
185+
F_L(5:6) = vel%L(1)*B%L(2:3) - vel%L(2:3)*B%L(1)
186+
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))
187+
188+
F_R(1) = U_R(2)
189+
F_R(2) = U_R(2)*vel%R(1) - B%R(1)*B%R(1) + pTot_R
190+
F_R(3:4) = U_R(2)*vel%R(2:3) - B%R(1)*B%R(2:3)
191+
F_R(5:6) = vel%R(1)*B%R(2:3) - vel%R(2:3)*B%R(1)
192+
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))
193+
! HLLD star-state fluxes via HLL jump relation
194+
F_starL = F_L + s_L*(U_starL - U_L)
195+
F_starR = F_R + s_R*(U_starR - U_R)
196+
! Alfven wave speeds bounding the rotational discontinuities
197+
s_starL = s_M - abs(B%L(1))/sqrt(rhoL_star)
198+
s_starR = s_M + abs(B%L(1))/sqrt(rhoR_star)
199+
! HLLD double-star (intermediate) states across rotational discontinuities
200+
sqrt_rhoL_star = sqrt(rhoL_star); sqrt_rhoR_star = sqrt(rhoR_star)
201+
vL_star = vel%L(2); wL_star = vel%L(3)
202+
vR_star = vel%R(2); wR_star = vel%R(3)
203+
204+
! (6) Compute the double-star states [Miyoshi Eqns. (59)-(62)]
205+
denom_ds = sqrt_rhoL_star + sqrt_rhoR_star
206+
sign_Bx = sign(1._wp, B%L(1))
207+
v_double = (sqrt_rhoL_star*vL_star + sqrt_rhoR_star*vR_star + (B%R(2) - B%L(2))*sign_Bx)/denom_ds
208+
w_double = (sqrt_rhoL_star*wL_star + sqrt_rhoR_star*wR_star + (B%R(3) - B%L(3))*sign_Bx)/denom_ds
209+
By_double = (sqrt_rhoL_star*B%R(2) + sqrt_rhoR_star*B%L(2) + sqrt_rhoL_star*sqrt_rhoR_star*(vR_star &
210+
& - vL_star)*sign_Bx)/denom_ds
211+
Bz_double = (sqrt_rhoL_star*B%R(3) + sqrt_rhoR_star*B%L(3) + sqrt_rhoL_star*sqrt_rhoR_star*(wR_star &
212+
& - wL_star)*sign_Bx)/denom_ds
213+
214+
E_doubleL = E_starL - sqrt_rhoL_star*((vL_star*B%L(2) + wL_star*B%L(3)) - (v_double*By_double &
215+
& + w_double*Bz_double))*sign_Bx
216+
E_doubleR = E_starR + sqrt_rhoR_star*((vR_star*B%R(2) + wR_star*B%R(3)) - (v_double*By_double &
217+
& + w_double*Bz_double))*sign_Bx
218+
E_double = 0.5_wp*(E_doubleL + E_doubleR)
219+
220+
U_doubleL = [rhoL_star, rhoL_star*s_M, rhoL_star*v_double, rhoL_star*w_double, By_double, Bz_double, &
221+
& E_double]
222+
U_doubleR = [rhoR_star, rhoR_star*s_M, rhoR_star*v_double, rhoR_star*w_double, By_double, Bz_double, &
223+
& E_double]
224+
225+
! Select HLLD flux region
226+
if (0.0_wp <= s_L) then
227+
F_hlld = F_L
228+
else if (0.0_wp <= s_starL) then
229+
F_hlld = F_L + s_L*(U_starL - U_L)
230+
else if (0.0_wp <= s_M) then
231+
F_hlld = F_starL + s_starL*(U_doubleL - U_starL)
232+
else if (0.0_wp <= s_starR) then
233+
F_hlld = F_starR + s_starR*(U_doubleR - U_starR)
234+
else if (0.0_wp <= s_R) then
235+
F_hlld = F_R + s_R*(U_starR - U_R)
236+
else
237+
F_hlld = F_R
238+
end if
239+
240+
! (12) Write HLLD flux to output arrays
241+
flux_rsx_vf(${SF('')}$, 1) = F_hlld(1) ! TODO multi-component
242+
! Momentum
243+
flux_rsx_vf(${SF('')}$, eqn_idx%cont%end + dir_idx(1)) = F_hlld(2)
244+
flux_rsx_vf(${SF('')}$, eqn_idx%cont%end + dir_idx(2)) = F_hlld(3)
245+
flux_rsx_vf(${SF('')}$, eqn_idx%cont%end + dir_idx(3)) = F_hlld(4)
246+
! Magnetic field
247+
if (n == 0) then
248+
flux_rsx_vf(${SF('')}$, eqn_idx%B%beg) = F_hlld(5)
249+
flux_rsx_vf(${SF('')}$, eqn_idx%B%beg + 1) = F_hlld(6)
250+
else
251+
flux_rsx_vf(${SF('')}$, eqn_idx%B%beg + dir_idx(1) - 1) = 0._wp
252+
flux_rsx_vf(${SF('')}$, eqn_idx%B%beg + dir_idx(2) - 1) = F_hlld(5)
253+
flux_rsx_vf(${SF('')}$, eqn_idx%B%beg + dir_idx(3) - 1) = F_hlld(6)
254+
end if
255+
! Energy
256+
flux_rsx_vf(${SF('')}$, eqn_idx%E) = F_hlld(7)
257+
! Volume fractions
258+
$:GPU_LOOP(parallelism='[seq]')
259+
do i = eqn_idx%adv%beg, eqn_idx%adv%end
260+
flux_rsx_vf(${SF('')}$, i) = 0._wp ! TODO multi-component (zero for now)
261+
end do
262+
263+
flux_src_rsx_vf(${SF('')}$, eqn_idx%adv%beg) = 0._wp
264+
end do
265+
end do
266+
end do
267+
$:END_GPU_PARALLEL_LOOP()
268+
end if
269+
#:endfor
270+
271+
call s_finalize_riemann_solver(flux_vf, flux_src_vf, flux_gsrc_vf, norm_dir)
272+
273+
end subroutine s_hlld_riemann_solver
274+
275+
end module m_riemann_solver_hlld

0 commit comments

Comments
 (0)