diff options
Diffstat (limited to 'src/GRHydro_HLLEPoly.F90')
| -rw-r--r-- | src/GRHydro_HLLEPoly.F90 | 83 |
1 files changed, 36 insertions, 47 deletions
diff --git a/src/GRHydro_HLLEPoly.F90 b/src/GRHydro_HLLEPoly.F90 index 084fb33..ed3c4f3 100644 --- a/src/GRHydro_HLLEPoly.F90 +++ b/src/GRHydro_HLLEPoly.F90 @@ -13,6 +13,7 @@ #include "cctk_Arguments.h" #include "cctk_Functions.h" +#include "GRHydro_Macros.h" #include "SpaceMask.h" /*@@ @@ -44,7 +45,7 @@ subroutine GRHydro_HLLE(CCTK_ARGUMENTS) CCTK_REAL, dimension(5) :: qdiff 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 CCTK_INT :: type_bits, trivial, not_trivial @@ -121,14 +122,8 @@ subroutine GRHydro_HLLE(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 @@ -160,8 +155,8 @@ subroutine GRHydro_HLLE(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 @@ -191,14 +186,14 @@ subroutine GRHydro_HLLE(CCTK_ARGUMENTS) velzminus(i+xoffset,j+yoffset,k+zoffset),& epsminus(i+xoffset,j+yoffset,k+zoffset),& w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),& - lamminus,psi4h*gxxh,psi4h*gxyh,psi4h*gxzh,psi4h*gyyh,& - psi4h*gyzh,psi4h*gzzh,& + lamminus,gxxh,gxyh,gxzh,gyyh,& + gyzh,gzzh,& usendh,avg_alp,avg_beta) call eigenvalues(GRHydro_eos_handle,rhoplus(i,j,k),& velxplus(i,j,k),velyplus(i,j,k),& velzplus(i,j,k),epsplus(i,j,k),w_lorentzplus(i,j,k),& - lamplus,psi4h*gxxh,psi4h*gxyh,psi4h*gxzh,psi4h*gyyh,& - psi4h*gyzh,psi4h*gzzh,& + lamplus,gxxh,gxyh,gxzh,gyyh,& + gyzh,gzzh,& usendh,avg_alp,avg_beta) call num_x_flux(consp(1),consp(2),consp(3),consp(4),consp(5),& fplus(1),fplus(2),fplus(3),fplus(4),& @@ -218,14 +213,14 @@ subroutine GRHydro_HLLE(CCTK_ARGUMENTS) velxminus(i+xoffset,j+yoffset,k+zoffset),& epsminus(i+xoffset,j+yoffset,k+zoffset),& w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),& - lamminus,psi4h*gyyh,psi4h*gyzh,psi4h*gxyh,psi4h*gzzh,& - psi4h*gxzh,psi4h*gxxh,& + lamminus,gyyh,gyzh,gxyh,gzzh,& + gxzh,gxxh,& usendh,avg_alp,avg_beta) call eigenvalues(GRHydro_eos_handle,rhoplus(i,j,k),& velyplus(i,j,k),velzplus(i,j,k),& velxplus(i,j,k),epsplus(i,j,k),w_lorentzplus(i,j,k),& - lamplus,psi4h*gyyh,psi4h*gyzh,psi4h*gxyh,psi4h*gzzh,& - psi4h*gxzh,psi4h*gxxh,& + lamplus,gyyh,gyzh,gxyh,gzzh,& + gxzh,gxxh,& usendh,avg_alp,avg_beta) call num_x_flux(consp(1),consp(3),consp(4),consp(2),consp(5),& fplus(1),fplus(3),fplus(4),fplus(2),& @@ -245,14 +240,14 @@ subroutine GRHydro_HLLE(CCTK_ARGUMENTS) velyminus(i+xoffset,j+yoffset,k+zoffset),& epsminus(i+xoffset,j+yoffset,k+zoffset),& w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),& - lamminus,psi4h*gzzh,psi4h*gxzh,psi4h*gyzh,& - psi4h*gxxh,psi4h*gxyh,psi4h*gyyh,& + lamminus,gzzh,gxzh,gyzh,& + gxxh,gxyh,gyyh,& usendh,avg_alp,avg_beta) call eigenvalues(GRHydro_eos_handle,rhoplus(i,j,k),& velzplus(i,j,k),velxplus(i,j,k),& velyplus(i,j,k),epsplus(i,j,k),w_lorentzplus(i,j,k),& - lamplus,psi4h*gzzh,psi4h*gxzh,psi4h*gyzh,& - psi4h*gxxh,psi4h*gxyh,psi4h*gyyh,& + lamplus,gzzh,gxzh,gyzh,& + gxxh,gxyh,gyyh,& usendh,avg_alp,avg_beta) call num_x_flux(consp(1),consp(4),consp(2),consp(3),consp(5),& fplus(1),fplus(4),fplus(2),fplus(3),& @@ -335,7 +330,7 @@ subroutine GRHydro_HLLE_Tracer(CCTK_ARGUMENTS) CCTK_REAL, dimension(number_of_tracers) :: qdiff 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 CCTK_INT :: type_bits, trivial, not_trivial @@ -404,18 +399,12 @@ subroutine GRHydro_HLLE_Tracer(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) + 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 @@ -441,14 +430,14 @@ subroutine GRHydro_HLLE_Tracer(CCTK_ARGUMENTS) velzminus(i+xoffset,j+yoffset,k+zoffset),& epsminus(i+xoffset,j+yoffset,k+zoffset),& w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),& - lamminus,psi4h*gxxh,psi4h*gxyh,psi4h*gxzh,psi4h*gyyh,& - psi4h*gyzh,psi4h*gzzh,& + lamminus,gxxh,gxyh,gxzh,gyyh,& + gyzh,gzzh,& usendh,avg_alp,avg_beta) call eigenvalues(GRHydro_eos_handle,rhoplus(i,j,k),& velxplus(i,j,k),velyplus(i,j,k),& velzplus(i,j,k),epsplus(i,j,k),w_lorentzplus(i,j,k),& - lamplus,psi4h*gxxh,psi4h*gxyh,psi4h*gxzh,psi4h*gyyh,& - psi4h*gyzh,psi4h*gzzh,& + lamplus,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,:) @@ -462,14 +451,14 @@ subroutine GRHydro_HLLE_Tracer(CCTK_ARGUMENTS) velxminus(i+xoffset,j+yoffset,k+zoffset),& epsminus(i+xoffset,j+yoffset,k+zoffset),& w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),& - lamminus,psi4h*gyyh,psi4h*gyzh,psi4h*gxyh,psi4h*gzzh,& - psi4h*gxzh,psi4h*gxxh,& + lamminus,gyyh,gyzh,gxyh,gzzh,& + gxzh,gxxh,& usendh,avg_alp,avg_beta) call eigenvalues(GRHydro_eos_handle,rhoplus(i,j,k),& velyplus(i,j,k),velzplus(i,j,k),& velxplus(i,j,k),epsplus(i,j,k),w_lorentzplus(i,j,k),& - lamplus,psi4h*gyyh,psi4h*gyzh,psi4h*gxyh,psi4h*gzzh,& - psi4h*gxzh,psi4h*gxxh,& + lamplus,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,:) @@ -483,14 +472,14 @@ subroutine GRHydro_HLLE_Tracer(CCTK_ARGUMENTS) velyminus(i+xoffset,j+yoffset,k+zoffset),& epsminus(i+xoffset,j+yoffset,k+zoffset),& w_lorentzminus(i+xoffset,j+yoffset,k+zoffset),& - lamminus,psi4h*gzzh,psi4h*gxzh,psi4h*gyzh,& - psi4h*gxxh,psi4h*gxyh,psi4h*gyyh,& + lamminus,gzzh,gxzh,gyzh,& + gxxh,gxyh,gyyh,& usendh,avg_alp,avg_beta) call eigenvalues(GRHydro_eos_handle,rhoplus(i,j,k),& velzplus(i,j,k),velxplus(i,j,k),& velyplus(i,j,k),epsplus(i,j,k),w_lorentzplus(i,j,k),& - lamplus,psi4h*gzzh,psi4h*gxzh,psi4h*gyzh,& - psi4h*gxxh,psi4h*gxyh,psi4h*gyyh,& + lamplus,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,:) |