Skip to content

Refactor m_riemann_solvers Module (HLLD Solver) #909

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 6 commits into from
Jul 3, 2025
Merged

Refactor m_riemann_solvers Module (HLLD Solver) #909

Changes from 4 commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
37650e2
HLLD refactor
Malmahrouqi3 Jun 29, 2025
3657f95
typo in HLLD U_doubleR
Malmahrouqi3 Jun 30, 2025
e3a22de
Corrected F_starR & F_starL
Malmahrouqi3 Jun 30, 2025
ee298e6
lint & format
Malmahrouqi3 Jun 30, 2025
63e818c
Merge branch 'master' into refactor-HLLD
sbryngelson Jul 1, 2025
93ad7f0
Merge branch 'master' into refactor-HLLD
sbryngelson Jul 3, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
127 changes: 35 additions & 92 deletions src/simulation/m_riemann_solvers.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3092,19 +3092,15 @@
! NOTE: unlike HLL, Bx, By, Bz are permutated by dir_idx for simpler logic
if (mhd) then
if (n == 0) then ! 1D: constant Bx; By, Bz as variables; only in x so not permutated
B%L(1) = Bx0
B%R(1) = Bx0
B%L(2) = qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg)
B%R(2) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, B_idx%beg)
B%L(3) = qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg + 1)
B%R(3) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, B_idx%beg + 1)
B%L = [Bx0, qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg), qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg + 1)]
B%R = [Bx0, qR_prim_rs${XYZ}$_vf(j + 1, k, l, B_idx%beg), qR_prim_rs${XYZ}$_vf(j + 1, k, l, B_idx%beg + 1)]

Check warning on line 3096 in src/simulation/m_riemann_solvers.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_riemann_solvers.fpp#L3095-L3096 

Added lines #L3095 - L3096 were not covered by tests
else ! 2D/3D: Bx, By, Bz as variables
B%L(1) = qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg + dir_idx(1) - 1)
B%R(1) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, B_idx%beg + dir_idx(1) - 1)
B%L(2) = qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg + dir_idx(2) - 1)
B%R(2) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, B_idx%beg + dir_idx(2) - 1)
B%L(3) = qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg + dir_idx(3) - 1)
B%R(3) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, B_idx%beg + dir_idx(3) - 1)
B%L = [qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg + dir_idx(1) - 1), &
qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg + dir_idx(2) - 1), &
qL_prim_rs${XYZ}$_vf(j, k, l, B_idx%beg + dir_idx(3) - 1)]
B%R = [qR_prim_rs${XYZ}$_vf(j + 1, k, l, B_idx%beg + dir_idx(1) - 1), &
qR_prim_rs${XYZ}$_vf(j + 1, k, l, B_idx%beg + dir_idx(2) - 1), &
qR_prim_rs${XYZ}$_vf(j + 1, k, l, B_idx%beg + dir_idx(3) - 1)]

Check warning on line 3103 in src/simulation/m_riemann_solvers.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_riemann_solvers.fpp#L3098-L3103 

Added lines #L3098 - L3103 were not covered by tests
end if
end if
Expand Down Expand Up @@ -3155,74 +3151,38 @@
E_starL = ((s_L - vel%L(1))*E%L - pTot_L*vel%L(1) + p_star*s_M)/(s_L - s_M)
E_starR = ((s_R - vel%R(1))*E%R - pTot_R*vel%R(1) + p_star*s_M)/(s_R - s_M)
! (5) Compute the left/right conserved state vectors
U_L(1) = rho%L
U_L(2) = rho%L*vel%L(1)
U_L(3) = rho%L*vel%L(2)
U_L(4) = rho%L*vel%L(3)
U_L(5) = B%L(2)
U_L(6) = B%L(3)
U_L(7) = E%L
U_R(1) = rho%R
U_R(2) = rho%R*vel%R(1)
U_R(3) = rho%R*vel%R(2)
U_R(4) = rho%R*vel%R(3)
U_R(5) = B%R(2)
U_R(6) = B%R(3)
U_R(7) = E%R
! (6) Compute the left/right star state vectors
U_starL(1) = rhoL_star
U_starL(2) = rhoL_star*s_M
U_starL(3) = rhoL_star*vel%L(2)
U_starL(4) = rhoL_star*vel%L(3)
U_starL(5) = B%L(2)
U_starL(6) = B%L(3)
U_starL(7) = E_starL
U_starR(1) = rhoR_star
U_starR(2) = rhoR_star*s_M
U_starR(3) = rhoR_star*vel%R(2)
U_starR(4) = rhoR_star*vel%R(3)
U_starR(5) = B%R(2)
U_starR(6) = B%R(3)
U_starR(7) = E_starR
! (7) Compute the left/right fluxes
F_L(1) = rho%L*vel%L(1)
F_L(2) = rho%L*vel%L(1)*vel%L(1) - B%L(1)*B%L(1) + pTot_L
F_L(3) = rho%L*vel%L(1)*vel%L(2) - B%L(1)*B%L(2)
F_L(4) = rho%L*vel%L(1)*vel%L(3) - B%L(1)*B%L(3)
F_L(5) = vel%L(1)*B%L(2) - vel%L(2)*B%L(1)
F_L(6) = vel%L(1)*B%L(3) - vel%L(3)*B%L(1)
! (5) Compute left/right state vectors and fluxes
U_L = [rho%L, rho%L*vel%L(1:3), B%L(2:3), E%L]
U_starL = [rhoL_star, rhoL_star*s_M, rhoL_star*vel%L(2:3), B%L(2:3), E_starL]
U_R = [rho%R, rho%R*vel%R(1:3), B%R(2:3), E%R]
U_starR = [rhoR_star, rhoR_star*s_M, rhoR_star*vel%R(2:3), B%R(2:3), E_starR]

Check warning on line 3158 in src/simulation/m_riemann_solvers.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_riemann_solvers.fpp#L3155-L3158 

Added lines #L3155 - L3158 were not covered by tests
! Compute the left/right fluxes
F_L(1) = U_L(2)
F_L(2) = U_L(2)*vel%L(1) - B%L(1)*B%L(1) + pTot_L
F_L(3:4) = U_L(2)*vel%L(2:3) - B%L(1)*B%L(2:3)
F_L(5:6) = vel%L(1)*B%L(2:3) - vel%L(2:3)*B%L(1)

Check warning on line 3164 in src/simulation/m_riemann_solvers.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_riemann_solvers.fpp#L3161-L3164 

Added lines #L3161 - L3164 were not covered by tests
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))
F_R(1) = rho%R*vel%R(1)
F_R(2) = rho%R*vel%R(1)*vel%R(1) - B%R(1)*B%R(1) + pTot_R
F_R(3) = rho%R*vel%R(1)*vel%R(2) - B%R(1)*B%R(2)
F_R(4) = rho%R*vel%R(1)*vel%R(3) - B%R(1)*B%R(3)
F_R(5) = vel%R(1)*B%R(2) - vel%R(2)*B%R(1)
F_R(6) = vel%R(1)*B%R(3) - vel%R(3)*B%R(1)
F_R(1) = U_R(2)
F_R(2) = U_R(2)*vel%R(1) - B%R(1)*B%R(1) + pTot_R
F_R(3:4) = U_R(2)*vel%R(2:3) - B%R(1)*B%R(2:3)
F_R(5:6) = vel%R(1)*B%R(2:3) - vel%R(2:3)*B%R(1)

Check warning on line 3170 in src/simulation/m_riemann_solvers.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_riemann_solvers.fpp#L3167-L3170 

Added lines #L3167 - L3170 were not covered by tests
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))
! (8) Compute the left/right star fluxes (note array operations)
! Compute the star flux using HLL relation
F_starL = F_L + s_L*(U_starL - U_L)
F_starR = F_R + s_R*(U_starR - U_R)
! (9) Compute the rotational (Alfvén) speeds
! Compute the rotational (Alfvén) speeds
s_starL = s_M - abs(B%L(1))/sqrt(rhoL_star)
s_starR = s_M + abs(B%L(1))/sqrt(rhoR_star)
! Compute the double–star states [Miyoshi Eqns. (59)-(62)]
sqrt_rhoL_star = sqrt(rhoL_star); sqrt_rhoR_star = sqrt(rhoR_star)
vL_star = vel%L(2); wL_star = vel%L(3)
vR_star = vel%R(2); wR_star = vel%R(3)
! (10) Compute the double–star states [Miyoshi Eqns. (59)-(62)]
sqrt_rhoL_star = sqrt(rhoL_star)
sqrt_rhoR_star = sqrt(rhoR_star)
! (6) Compute the double–star states [Miyoshi Eqns. (59)-(62)]
denom_ds = sqrt_rhoL_star + sqrt_rhoR_star
sign_Bx = sign(1._wp, B%L(1))
vL_star = vel%L(2)
wL_star = vel%L(3)
vR_star = vel%R(2)
wR_star = vel%R(3)
v_double = (sqrt_rhoL_star*vL_star + sqrt_rhoR_star*vR_star + (B%R(2) - B%L(2))*sign_Bx)/denom_ds
w_double = (sqrt_rhoL_star*wL_star + sqrt_rhoR_star*wR_star + (B%R(3) - B%L(3))*sign_Bx)/denom_ds
By_double = (sqrt_rhoL_star*B%R(2) + sqrt_rhoR_star*B%L(2) + sqrt_rhoL_star*sqrt_rhoR_star*(vR_star - vL_star)*sign_Bx)/denom_ds
Expand All @@ -3232,21 +3192,8 @@
E_doubleR = E_starR + sqrt_rhoR_star*((vR_star*B%R(2) + wR_star*B%R(3)) - (v_double*By_double + w_double*Bz_double))*sign_Bx
E_double = 0.5_wp*(E_doubleL + E_doubleR)
U_doubleL(1) = rhoL_star
U_doubleL(2) = rhoL_star*s_M
U_doubleL(3) = rhoL_star*v_double
U_doubleL(4) = rhoL_star*w_double
U_doubleL(5) = By_double
U_doubleL(6) = Bz_double
U_doubleL(7) = E_double
U_doubleR(1) = rhoR_star
U_doubleR(2) = rhoR_star*s_M
U_doubleR(3) = rhoR_star*v_double
U_doubleR(4) = rhoR_star*w_double
U_doubleR(5) = By_double
U_doubleR(6) = Bz_double
U_doubleR(7) = E_double
U_doubleL = [rhoL_star, rhoL_star*s_M, rhoL_star*v_double, rhoL_star*w_double, By_double, Bz_double, E_double]
U_doubleR = [rhoR_star, rhoR_star*s_M, rhoR_star*v_double, rhoR_star*w_double, By_double, Bz_double, E_double]

Check warning on line 3196 in src/simulation/m_riemann_solvers.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_riemann_solvers.fpp#L3195-L3196 

Added lines #L3195 - L3196 were not covered by tests
! (11) Choose HLLD flux based on wave-speed regions
if (0.0_wp <= s_L) then
Expand All @@ -3267,16 +3214,12 @@
! Mass
flux_rs${XYZ}$_vf(j, k, l, 1) = F_hlld(1) ! TODO multi-component
! Momentum
flux_rs${XYZ}$_vf(j, k, l, contxe + dir_idx(1)) = F_hlld(2)
flux_rs${XYZ}$_vf(j, k, l, contxe + dir_idx(2)) = F_hlld(3)
flux_rs${XYZ}$_vf(j, k, l, contxe + dir_idx(3)) = F_hlld(4)
flux_rs${XYZ}$_vf(j, k, l, [contxe + dir_idx(1), contxe + dir_idx(2), contxe + dir_idx(3)]) = F_hlld([2, 3, 4])

Check warning on line 3217 in src/simulation/m_riemann_solvers.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_riemann_solvers.fpp#L3217 

Added line #L3217 was not covered by tests
! Magnetic field
if (n == 0) then
flux_rs${XYZ}$_vf(j, k, l, B_idx%beg) = F_hlld(5)
flux_rs${XYZ}$_vf(j, k, l, B_idx%beg + 1) = F_hlld(6)
flux_rs${XYZ}$_vf(j, k, l, [B_idx%beg, B_idx%beg + 1]) = F_hlld([5, 6])

Check warning on line 3220 in src/simulation/m_riemann_solvers.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_riemann_solvers.fpp#L3220 

Added line #L3220 was not covered by tests
else
flux_rs${XYZ}$_vf(j, k, l, B_idx%beg + dir_idx(2) - 1) = F_hlld(5)
flux_rs${XYZ}$_vf(j, k, l, B_idx%beg + dir_idx(3) - 1) = F_hlld(6)
flux_rs${XYZ}$_vf(j, k, l, [B_idx%beg + dir_idx(2) - 1, B_idx%beg + dir_idx(3) - 1]) = F_hlld([5, 6])

Check warning on line 3222 in src/simulation/m_riemann_solvers.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_riemann_solvers.fpp#L3222 

Added line #L3222 was not covered by tests
end if
! Energy
flux_rs${XYZ}$_vf(j, k, l, E_idx) = F_hlld(7)
Expand Down
Loading