diff options
Diffstat (limited to 'src/GRHydro_HLLE.F90')
| -rw-r--r-- | src/GRHydro_HLLE.F90 | 80 |
1 files changed, 33 insertions, 47 deletions
diff --git a/src/GRHydro_HLLE.F90 b/src/GRHydro_HLLE.F90 index d32f4a7..b0b5ba4 100644 --- a/src/GRHydro_HLLE.F90 +++ b/src/GRHydro_HLLE.F90 @@ -12,7 +12,7 @@ #include "cctk_Parameters.h" #include "cctk_Arguments.h" #include "cctk_Functions.h" - +#include "GRHydro_Macros.h" #include "SpaceMask.h" /*@@ @@ -47,7 +47,7 @@ subroutine GRHydro_HLLEGeneral(CCTK_ARGUMENTS) CCTK_REAL, dimension(6) :: prim_p, prim_m CCTK_REAL :: charmin, charmax, charpm,avg_alp,avg_det CCTK_REAL :: gxxh, gxyh, gxzh, gyyh, gyzh, gzzh, uxxh, uxyh, & - uxzh, uyyh, uyzh, uzzh, avg_beta, usendh, alp_l, alp_r, psi4h, & + uxzh, uyyh, uyzh, uzzh, avg_beta, usendh, alp_l, alp_r, & cs2_p, cs2_m, dpdeps_p, dpdeps_m CCTK_INT :: type_bits, trivial, not_trivial @@ -154,14 +154,7 @@ subroutine GRHydro_HLLEGeneral(CCTK_ARGUMENTS) gyzh = 0.5d0 * (gyz(i+xoffset,j+yoffset,k+zoffset) + gyz(i,j,k)) gzzh = 0.5d0 * (gzz(i+xoffset,j+yoffset,k+zoffset) + gzz(i,j,k)) - if (conformal_state .eq. 0) then - psi4h = 1.d0 - else - psi4h = (0.5d0*(psi(i,j,k)+psi(i+xoffset,j+yoffset,k+zoffset)))**4 - end if - - call SpatialDeterminant(psi4h*gxxh,psi4h*gxyh,psi4h*gxzh,& - psi4h*gyyh,psi4h*gyzh,psi4h*gzzh,avg_det) + avg_det = SPATIAL_DETERMINANT(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh) !!$ If the Riemann problem is trivial, just calculate the fluxes from the !!$ left state and skip to the next cell @@ -190,8 +183,8 @@ subroutine GRHydro_HLLEGeneral(CCTK_ARGUMENTS) else !!! The end of this branch is right at the bottom of the routine call UpperMetric(uxxh, uxyh, uxzh, uyyh, uyzh, uzzh, & - avg_det,psi4h*gxxh, psi4h*gxyh, psi4h*gxzh, & - psi4h*gyyh, psi4h*gyzh, psi4h*gzzh) + avg_det,gxxh,gxyh,gxzh, & + gyyh,gyzh, gzzh) if (flux_direction == 1) then usendh = uxxh @@ -217,14 +210,14 @@ subroutine GRHydro_HLLEGeneral(CCTK_ARGUMENTS) call eigenvalues_general(& prim_m(2),prim_m(3),prim_m(4),cs2_m, & lamminus,& - psi4h*gxxh,psi4h*gxyh,psi4h*gxzh,& - psi4h*gyyh,psi4h*gyzh,psi4h*gzzh,& + gxxh,gxyh,gxzh,& + gyyh,gyzh,gzzh,& usendh,avg_alp,avg_beta) call eigenvalues_general(& prim_p(2),prim_p(3),prim_p(4),cs2_p, & lamplus,& - psi4h*gxxh,psi4h*gxyh,psi4h*gxzh,& - psi4h*gyyh,psi4h*gyzh,psi4h*gzzh,& + gxxh,gxyh,gxzh,& + gyyh,gyzh,gzzh,& usendh,avg_alp,avg_beta) call num_x_flux(cons_p(1),cons_p(2),cons_p(3),cons_p(4),cons_p(5),& fplus(1),fplus(2),fplus(3),fplus(4),fplus(5),& @@ -238,14 +231,14 @@ subroutine GRHydro_HLLEGeneral(CCTK_ARGUMENTS) call eigenvalues_general(& prim_m(3),prim_m(4),prim_m(2),cs2_m, & lamminus,& - psi4h*gyyh,psi4h*gyzh,psi4h*gxyh,& - psi4h*gzzh,psi4h*gxzh,psi4h*gxxh,& + gyyh,gyzh,gxyh,& + gzzh,gxzh,gxxh,& usendh,avg_alp,avg_beta) call eigenvalues_general(& prim_p(3),prim_p(4),prim_p(2),cs2_p, & lamplus,& - psi4h*gyyh,psi4h*gyzh,psi4h*gxyh,& - psi4h*gzzh,psi4h*gxzh,psi4h*gxxh,& + gyyh,gyzh,gxyh,& + gzzh,gxzh,gxxh,& usendh,avg_alp,avg_beta) call num_x_flux(cons_p(1),cons_p(3),cons_p(4),cons_p(2),cons_p(5),& fplus(1),fplus(3),fplus(4),fplus(2),fplus(5),& @@ -259,14 +252,14 @@ subroutine GRHydro_HLLEGeneral(CCTK_ARGUMENTS) call eigenvalues_general(& prim_m(4),prim_m(2),prim_m(3),cs2_m, & lamminus,& - psi4h*gzzh,psi4h*gxzh,psi4h*gyzh,& - psi4h*gxxh,psi4h*gxyh,psi4h*gyyh,& + gzzh,gxzh,gyzh,& + gxxh,gxyh,gyyh,& usendh,avg_alp,avg_beta) call eigenvalues_general(& prim_p(4),prim_p(2),prim_p(3),cs2_p, & lamplus,& - psi4h*gzzh,psi4h*gxzh,psi4h*gyzh,& - psi4h*gxxh,psi4h*gxyh,psi4h*gyyh,& + gzzh,gxzh,gyzh,& + gxxh,gxyh,gyyh,& usendh,avg_alp,avg_beta) call num_x_flux(cons_p(1),cons_p(4),cons_p(2),cons_p(3),cons_p(5),& fplus(1),fplus(4),fplus(2),fplus(3),fplus(5), & @@ -358,7 +351,7 @@ subroutine GRHydro_HLLE_TracerGeneral(CCTK_ARGUMENTS) CCTK_REAL, dimension(6) :: prim_p, prim_m CCTK_REAL :: charmin, charmax, charpm,avg_alp,avg_det CCTK_REAL :: gxxh, gxyh, gxzh, gyyh, gyzh, gzzh, uxxh, uxyh, & - uxzh, uyyh, uyzh, uzzh, avg_beta, usendh, alp_l, alp_r, psi4h, & + uxzh, uyyh, uyzh, uzzh, avg_beta, usendh, alp_l, alp_r, & cs2_p, cs2_m, dpdeps_p, dpdeps_m integer tadmor @@ -438,21 +431,14 @@ subroutine GRHydro_HLLE_TracerGeneral(CCTK_ARGUMENTS) gyzh = 0.5d0 * (gyz(i+xoffset,j+yoffset,k+zoffset) + gyz(i,j,k)) gzzh = 0.5d0 * (gzz(i+xoffset,j+yoffset,k+zoffset) + gzz(i,j,k)) - if (conformal_state .eq. 0) then - psi4h = 1.d0 - else - psi4h = (0.5d0*(psi(i,j,k)+psi(i+xoffset,j+yoffset,k+zoffset)))**4 - end if - - call SpatialDeterminant(psi4h*gxxh,psi4h*gxyh,psi4h*gxzh,& - psi4h*gyyh,psi4h*gyzh,psi4h*gzzh,avg_det) + avg_det = SPATIAL_DETERMINANT(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh) !!$ If the Riemann problem is trivial, just calculate the fluxes from the !!$ left state and skip to the next cell call UpperMetric(uxxh, uxyh, uxzh, uyyh, uyzh, uzzh, & - avg_det,psi4h*gxxh, psi4h*gxyh, psi4h*gxzh, & - psi4h*gyyh, psi4h*gyzh, psi4h*gzzh) + avg_det,gxxh, gxyh, gxzh, & + gyyh, gyzh, gzzh) if (flux_direction == 1) then usendh = uxxh @@ -470,14 +456,14 @@ subroutine GRHydro_HLLE_TracerGeneral(CCTK_ARGUMENTS) call eigenvalues_general(& prim_m(2),prim_m(3),prim_m(4),cs2_m, & lamminus,& - psi4h*gxxh,psi4h*gxyh,psi4h*gxzh,& - psi4h*gyyh,psi4h*gyzh,psi4h*gzzh,& + gxxh,gxyh,gxzh,& + gyyh,gyzh,gzzh,& usendh,avg_alp,avg_beta) call eigenvalues_general(& prim_p(2),prim_p(3),prim_p(4),cs2_p, & lamplus,& - psi4h*gxxh,psi4h*gxyh,psi4h*gxzh,& - psi4h*gyyh,psi4h*gyzh,psi4h*gzzh,& + gxxh,gxyh,gxzh,& + gyyh,gyzh,gzzh,& usendh,avg_alp,avg_beta) fplus(:) = (velxplus(i,j,k) - avg_beta / avg_alp) * & cons_tracerplus(i,j,k,:) @@ -487,14 +473,14 @@ subroutine GRHydro_HLLE_TracerGeneral(CCTK_ARGUMENTS) call eigenvalues_general(& prim_m(3),prim_m(4),prim_m(2),cs2_m, & lamminus,& - psi4h*gyyh,psi4h*gyzh,psi4h*gxyh,& - psi4h*gzzh,psi4h*gxzh,psi4h*gxxh,& + gyyh,gyzh,gxyh,& + gzzh,gxzh,gxxh,& usendh,avg_alp,avg_beta) call eigenvalues_general(& prim_p(3),prim_p(4),prim_p(2),cs2_p, & lamplus,& - psi4h*gyyh,psi4h*gyzh,psi4h*gxyh,& - psi4h*gzzh,psi4h*gxzh,psi4h*gxxh,& + gyyh,gyzh,gxyh,& + gzzh,gxzh,gxxh,& usendh,avg_alp,avg_beta) fplus(:) = (velyplus(i,j,k) - avg_beta / avg_alp) * & cons_tracerplus(i,j,k,:) @@ -504,14 +490,14 @@ subroutine GRHydro_HLLE_TracerGeneral(CCTK_ARGUMENTS) call eigenvalues_general(& prim_m(4),prim_m(2),prim_m(3),cs2_m, & lamminus,& - psi4h*gzzh,psi4h*gxzh,psi4h*gyzh,& - psi4h*gxxh,psi4h*gxyh,psi4h*gyyh,& + gzzh,gxzh,gyzh,& + gxxh,gxyh,gyyh,& usendh,avg_alp,avg_beta) call eigenvalues_general(& prim_p(4),prim_p(2),prim_p(3),cs2_p, & lamplus,& - psi4h*gzzh,psi4h*gxzh,psi4h*gyzh,& - psi4h*gxxh,psi4h*gxyh,psi4h*gyyh,& + gzzh,gxzh,gyzh,& + gxxh,gxyh,gyyh,& usendh,avg_alp,avg_beta) fplus(:) = (velzplus(i,j,k) - avg_beta / avg_alp) * & cons_tracerplus(i,j,k,:) |