Skip to content

Commit 3fff3b1

Browse files
Malmahrouqi3Malmahrouqi3
committed
s_hll_riemann_solver passed Hypoelasticity tests
1 parent 062d51d commit 3fff3b1

File tree

1 file changed

+41
-56
lines changed

1 file changed

+41
-56
lines changed

‎src/simulation/m_riemann_solvers.fpp‎

Lines changed: 41 additions & 56 deletions
Original file line numberDiff line numberDiff line change
@@ -474,31 +474,19 @@ contains
474474
!$acc loop seq
475475
do i = 1, 2
476476
Re_L(i) = dflt_real
477-
478-
if (Re_size(i) > 0) Re_L(i) = 0._wp
479-
477+
Re_R(i) = dflt_real
478+
if (Re_size(i) > 0) then
479+
Re_L(i) = 0._wp
480+
Re_R(i) = 0._wp
481+
end if
480482
!$acc loop seq
481483
do q = 1, Re_size(i)
482484
Re_L(i) = alpha_L(Re_idx(i, q))/Res(i, q) &
483485
+ Re_L(i)
484-
end do
485-
486-
Re_L(i) = 1._wp/max(Re_L(i), sgm_eps)
487-
488-
end do
489-
490-
!$acc loop seq
491-
do i = 1, 2
492-
Re_R(i) = dflt_real
493-
494-
if (Re_size(i) > 0) Re_R(i) = 0._wp
495-
496-
!$acc loop seq
497-
do q = 1, Re_size(i)
498486
Re_R(i) = alpha_R(Re_idx(i, q))/Res(i, q) &
499487
+ Re_R(i)
500488
end do
501-
489+
Re_L(i) = 1._wp/max(Re_L(i), sgm_eps)
502490
Re_R(i) = 1._wp/max(Re_R(i), sgm_eps)
503491
end do
504492
end if
@@ -552,44 +540,41 @@ contains
552540
E_R = rho_R*E_R + 5e-1*rho_R*vel_R_rms
553541
H_L = (E_L + pres_L)/rho_L
554542
H_R = (E_R + pres_R)/rho_R
555-
elseif (mhd .and. relativity) then
556-
Ga%L = 1._wp/sqrt(1._wp - vel_L_rms)
557-
Ga%R = 1._wp/sqrt(1._wp - vel_R_rms)
558-
vdotB%L = vel_L(1)*B%L(1) + vel_L(2)*B%L(2) + vel_L(3)*B%L(3)
559-
vdotB%R = vel_R(1)*B%R(1) + vel_R(2)*B%R(2) + vel_R(3)*B%R(3)
560-
561-
b4%L(1) = B%L(1)/Ga%L + Ga%L*vel_L(1)*vdotB%L
562-
b4%L(2) = B%L(2)/Ga%L + Ga%L*vel_L(2)*vdotB%L
563-
b4%L(3) = B%L(3)/Ga%L + Ga%L*vel_L(3)*vdotB%L
564-
b4%R(1) = B%R(1)/Ga%R + Ga%R*vel_R(1)*vdotB%R
565-
b4%R(2) = B%R(2)/Ga%R + Ga%R*vel_R(2)*vdotB%R
566-
b4%R(3) = B%R(3)/Ga%R + Ga%R*vel_R(3)*vdotB%R
567-
B2%L = B%L(1)**2._wp + B%L(2)**2._wp + B%L(3)**2._wp
568-
B2%R = B%R(1)**2._wp + B%R(2)**2._wp + B%R(3)**2._wp
569-
570-
pres_mag%L = 0.5_wp*(B2%L/Ga%L**2._wp + vdotB%L**2._wp)
571-
pres_mag%R = 0.5_wp*(B2%R/Ga%R**2._wp + vdotB%R**2._wp)
572-
573-
! Hard-coded EOS
574-
H_L = 1._wp + (gamma_L + 1)*pres_L/rho_L
575-
H_R = 1._wp + (gamma_R + 1)*pres_R/rho_R
576-
577-
cm%L(1) = (rho_L*H_L*Ga%L**2 + B2%L)*vel_L(1) - vdotB%L*B%L(1)
578-
cm%L(2) = (rho_L*H_L*Ga%L**2 + B2%L)*vel_L(2) - vdotB%L*B%L(2)
579-
cm%L(3) = (rho_L*H_L*Ga%L**2 + B2%L)*vel_L(3) - vdotB%L*B%L(3)
580-
cm%R(1) = (rho_R*H_R*Ga%R**2 + B2%R)*vel_R(1) - vdotB%R*B%R(1)
581-
cm%R(2) = (rho_R*H_R*Ga%R**2 + B2%R)*vel_R(2) - vdotB%R*B%R(2)
582-
cm%R(3) = (rho_R*H_R*Ga%R**2 + B2%R)*vel_R(3) - vdotB%R*B%R(3)
583-
584-
E_L = rho_L*H_L*Ga%L**2 - pres_L + 0.5_wp*(B2%L + vel_L_rms*B2%L - vdotB%L**2._wp) - rho_L*Ga%L
585-
E_R = rho_R*H_R*Ga%R**2 - pres_R + 0.5_wp*(B2%R + vel_R_rms*B2%R - vdotB%R**2._wp) - rho_R*Ga%R
586-
elseif (mhd .and. .not. relativity) then
587-
pres_mag%L = 0.5_wp*(B%L(1)**2._wp + B%L(2)**2._wp + B%L(3)**2._wp)
588-
pres_mag%R = 0.5_wp*(B%R(1)**2._wp + B%R(2)**2._wp + B%R(3)**2._wp)
589-
E_L = gamma_L*pres_L + pi_inf_L + 0.5_wp*rho_L*vel_L_rms + qv_L + pres_mag%L
590-
E_R = gamma_R*pres_R + pi_inf_R + 0.5_wp*rho_R*vel_R_rms + qv_R + pres_mag%R ! includes magnetic energy
591-
H_L = (E_L + pres_L - pres_mag%L)/rho_L
592-
H_R = (E_R + pres_R - pres_mag%R)/rho_R ! stagnation enthalpy here excludes magnetic energy (only used to find speed of sound)
543+
elseif (mhd) then
544+
if (relativity) then
545+
Ga%L = 1._wp/sqrt(1._wp - vel_L_rms)
546+
Ga%R = 1._wp/sqrt(1._wp - vel_R_rms)
547+
vdotB%L = vel_L(1)*B%L(1) + vel_L(2)*B%L(2) + vel_L(3)*B%L(3)
548+
vdotB%R = vel_R(1)*B%R(1) + vel_R(2)*B%R(2) + vel_R(3)*B%R(3)
549+
550+
!acc loop seq
551+
do i = 1, 3
552+
b4%L(i) = B%L(i)/Ga%L + Ga%L*vel_L(i)*vdotB%L
553+
b4%R(i) = B%R(i)/Ga%R + Ga%R*vel_R(i)*vdotB%R
554+
end do
555+
556+
B2%L = sum(B%L**2._wp)
557+
B2%R = sum(B%R**2._wp)
558+
pres_mag%L = 0.5_wp*(B2%L/Ga%L**2._wp + vdotB%L**2._wp)
559+
pres_mag%R = 0.5_wp*(B2%R/Ga%R**2._wp + vdotB%R**2._wp)
560+
! Hard-coded EOS
561+
H_L = 1._wp + (gamma_L + 1)*pres_L/rho_L
562+
H_R = 1._wp + (gamma_R + 1)*pres_R/rho_R
563+
!acc loop seq
564+
do i = 1, 3
565+
cm%L(i) = (rho_L*H_L*Ga%L**2 + B2%L)*vel_L(i) - vdotB%L*B%L(i)
566+
cm%R(i) = (rho_R*H_R*Ga%R**2 + B2%R)*vel_R(i) - vdotB%R*B%R(i)
567+
end do
568+
E_L = rho_L*H_L*Ga%L**2 - pres_L + 0.5_wp*(B2%L + vel_L_rms*B2%L - vdotB%L**2._wp) - rho_L*Ga%L
569+
E_R = rho_R*H_R*Ga%R**2 - pres_R + 0.5_wp*(B2%R + vel_R_rms*B2%R - vdotB%R**2._wp) - rho_R*Ga%R
570+
elseif (.not. relativity) then
571+
pres_mag%L = 0.5_wp*sum(B%L**2._wp)
572+
pres_mag%R = 0.5_wp*sum(B%R**2._wp)
573+
E_L = gamma_L*pres_L + pi_inf_L + 0.5_wp*rho_L*vel_L_rms + qv_L + pres_mag%L
574+
E_R = gamma_R*pres_R + pi_inf_R + 0.5_wp*rho_R*vel_R_rms + qv_R + pres_mag%R ! includes magnetic energy
575+
H_L = (E_L + pres_L - pres_mag%L)/rho_L
576+
H_R = (E_R + pres_R - pres_mag%R)/rho_R ! stagnation enthalpy here excludes magnetic energy (only used to find speed of sound)
577+
end if
593578
else
594579
E_L = gamma_L*pres_L + pi_inf_L + 5e-1*rho_L*vel_L_rms + qv_L
595580
E_R = gamma_R*pres_R + pi_inf_R + 5e-1*rho_R*vel_R_rms + qv_R

0 commit comments

Comments
 (0)