diff options
Diffstat (limited to 'src/GRHydro_HLLEM.F90')
| -rw-r--r-- | src/GRHydro_HLLEM.F90 | 25 |
1 files changed, 20 insertions, 5 deletions
diff --git a/src/GRHydro_HLLEM.F90 b/src/GRHydro_HLLEM.F90 index 1e274eb..3c694ca 100644 --- a/src/GRHydro_HLLEM.F90 +++ b/src/GRHydro_HLLEM.F90 @@ -55,6 +55,7 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) CCTK_REAL :: pressstarp,pressstarm,velxlowp,velxlowm,velylowp,velylowm,velzlowp,velzlowm CCTK_REAL :: psidcp, psidcm, psidcf, psidcdiff, psidcfp, psidcfm + CCTK_REAL :: charmax_dc, charmin_dc, charpm_dc CCTK_INT :: type_bits, trivial, not_trivial @@ -440,13 +441,27 @@ subroutine GRHydro_HLLEM(CCTK_ARGUMENTS) end if case(1) !!$ Wavespeeds for psidc are +/-c, not Fast Magnetosonic? - psidcf = 0.5d0 * (1.d0 * psidcfp - (-1.d0) * psidcfm + & - 1.d0 * (-1.d0) * psidcdiff) + !!$ psidcf = 0.5d0 * (1.d0 * psidcfp - (-1.d0) * psidcfm + & + !!$ 1.d0 * (-1.d0) * psidcdiff) + + !!$ The fastest speed for psidc comes from the condition + !!$ that the normal vector to the characteristic hypersurface + !!$ be spacelike (Eq. 60 of Anton et al.) + + charmax_dc = sqrt(usendh) - avg_beta/avg_alp + charmin_dc = -1.d0*sqrt(usendh) - avg_beta/avg_alp + charpm_dc = charmax_dc - charmin_dc + + psidcf = (charmax_dc * psidcfp - charmin_dc * psidcfm + & + charmax_dc * charmin_dc * psidcdiff) / charpm_dc + if(decouple_normal_Bfield .ne. 0) then ! same expression for HLLE and LLF !!$ B^i field decouples from the others and has same propagation - !!$ speed as divergence - f1(5+flux_direction) = 0.5d0 * (1.d0 * fplus(5+flux_direction) - (-1.d0) * fminus(5+flux_direction) + & - 1.d0 * (-1.d0) * qdiff(5+flux_direction)) + !!$ speed as divergence -null direction, + !!$ \pm sqrt(g^{xx}} - beta^x/alpha + f1(5+flux_direction) = (charmax_dc * fplus(5+flux_direction) & + - charmin_dc * fminus(5+flux_direction) + & + charmax_dc * charmin_dc * qdiff(5+flux_direction)) / charpm_dc end if case(2) charmax = setcharmax |