diff options
Diffstat (limited to 'src/GRHydro_HLLEM.F90')
-rw-r--r-- | src/GRHydro_HLLEM.F90 | 34 |
1 files changed, 25 insertions, 9 deletions
diff --git a/src/GRHydro_HLLEM.F90 b/src/GRHydro_HLLEM.F90 index 200575f..7e9c541 100644 --- a/src/GRHydro_HLLEM.F90 +++ b/src/GRHydro_HLLEM.F90 @@ -33,6 +33,7 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) USE GRHydro_EigenproblemM + USE GRHydro_Scalars implicit none @@ -43,7 +44,7 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) CCTK_REAL, dimension(8) :: cons_p,cons_m,fplus,fminus,f1,qdiff CCTK_REAL, dimension(10) :: prim_p, prim_m CCTK_REAL, dimension(5) :: lamminus,lamplus - CCTK_REAL :: charmin, charmax, charpm,avg_alp,avg_det, sdet + CCTK_REAL :: charmin, charmax, charpm,chartop,avg_alp,avg_det, sdet CCTK_REAL :: gxxh, gxyh, gxzh, gyyh, gyzh, gzzh, uxxh, uxyh, & uxzh, uyyh, uyzh, uzzh, avg_beta, usendh, alp_l, alp_r, & cs2_p, cs2_m, dpdeps_p, dpdeps_m @@ -404,6 +405,8 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) charmax = max(0.d0, lamplus(1), lamplus(2), lamplus(3), & lamplus(4),lamplus(5), lamminus(1),lamminus(2),lamminus(3),& lamminus(4),lamminus(5)) + + chartop = max(-charmin,charmax) charpm = charmax - charmin @@ -413,8 +416,12 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) qdiff(m) = cons_m(m) - cons_p(m) - f1(m) = (charmax * fplus(m) - charmin * fminus(m) + & - charmax * charmin * qdiff(m)) / charpm + if (HLLE) then + f1(m) = (charmax * fplus(m) - charmin * fminus(m) + & + charmax * charmin * qdiff(m)) / charpm + else if (LLF) then + f1(m) = 0.5d0 * (fplus(m) + fminus(m) - chartop * qdiff(m)) + end if end do @@ -422,17 +429,26 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) psidcdiff = psidcm - psidcp select case(whichpsidcspeed) case(0) - psidcf = (charmax * psidcfp - charmin * psidcfm + & - charmax * charmin * psidcdiff) / charpm + if (HLLE) then + psidcf = (charmax * psidcfp - charmin * psidcfm + & + charmax * charmin * psidcdiff) / charpm + else if (LLF) then + psidcf = 0.5d0 * (psidcfp + psidcfm - chartop * psidcdiff) + end if case(1) !!$ Wavespeeds for psidc are +/-c, not Fast Magnetosonic? - psidcf = (1.d0 * psidcfp - (-1.d0) * psidcfm + & - 1.d0 * (-1.d0) * psidcdiff) / 2.d0 + psidcf = 0.5d0 * (1.d0 * psidcfp - (-1.d0) * psidcfm + & + 1.d0 * (-1.d0) * psidcdiff) case(2) charmax = setcharmax charmin = setcharmin - psidcf = (charmax * psidcfp - charmin * psidcfm + & - charmax * charmin * psidcdiff) / charpm + if (HLLE) then + psidcf = (charmax * psidcfp - charmin * psidcfm + & + charmax * charmin * psidcdiff) / charpm + else if (LLF) then + chartop = max(-charmin,charmax) + psidcf = 0.5d0 * (psidcfp + psidcfm - chartop * psidcdiff) + end if end select endif |