diff options
Diffstat (limited to 'src/GRHydro_Prim2Con.F90')
-rw-r--r-- | src/GRHydro_Prim2Con.F90 | 131 |
1 files changed, 24 insertions, 107 deletions
diff --git a/src/GRHydro_Prim2Con.F90 b/src/GRHydro_Prim2Con.F90 index 7951846..268e8e4 100644 --- a/src/GRHydro_Prim2Con.F90 +++ b/src/GRHydro_Prim2Con.F90 @@ -11,6 +11,7 @@ #include "cctk_Arguments.h" #include "cctk_Parameters.h" #include "cctk_Functions.h" +#include "GRHydro_Macros.h" #define velx(i,j,k) vel(i,j,k,1) #define vely(i,j,k) vel(i,j,k,2) @@ -63,25 +64,8 @@ subroutine primitive2conservative(CCTK_ARGUMENTS) gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset)) gzzr = 0.5d0 * (gzz(i,j,k) + gzz(i+xoffset,j+yoffset,k+zoffset)) - if (.not.(conformal_state .eq. 0)) then - psi4r = (0.5d0*(psi(i,j,k)+psi(i+xoffset,j+yoffset,k+zoffset)))**4 - psi4l = (0.5d0*(psi(i,j,k)+psi(i-xoffset,j-yoffset,k-zoffset)))**4 - gxxr = gxxr * psi4r - gxyr = gxyr * psi4r - gxzr = gxzr * psi4r - gyyr = gyyr * psi4r - gyzr = gyzr * psi4r - gzzr = gzzr * psi4r - gxxl = gxxl * psi4l - gxyl = gxyl * psi4l - gxzl = gxzl * psi4l - gyyl = gyyl * psi4l - gyzl = gyzl * psi4l - gzzl = gzzl * psi4l - end if - - call SpatialDeterminant(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl,avg_detl) - call SpatialDeterminant(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr,avg_detr) + avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl) + avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr) call prim2con(GRHydro_eos_handle, gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, & avg_detl,densminus(i,j,k),sxminus(i,j,k),& @@ -201,23 +185,17 @@ subroutine Primitive2ConservativeCells(CCTK_ARGUMENTS) do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1 do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1 - if (conformal_state .eq. 0) then - psi4pt = 1d0 - call SpatialDeterminant(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),& - gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),det) - else - psi4pt = psi(i,j,k)**4 - call SpatialDeterminant(psi4pt*gxx(i,j,k),psi4pt*gxy(i,j,k),& - psi4pt*gxz(i,j,k),psi4pt*gyy(i,j,k),psi4pt*gyz(i,j,k),& - psi4pt*gzz(i,j,k),det) - end if - call prim2con(GRHydro_eos_handle,psi4pt*gxx(i,j,k),& - psi4pt*gxy(i,j,k),psi4pt*gxz(i,j,k),& - psi4pt*gyy(i,j,k),psi4pt*gyz(i,j,k),psi4pt*gzz(i,j,k),& - det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& - tau(i,j,k),rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),& - eps(i,j,k),press(i,j,k),w_lorentz(i,j,k)) - + psi4pt = 1.0d0 + det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k), \ + gyy(i,j,k),gyz(i,j,k),gzz(i,j,k)) + + call prim2con(GRHydro_eos_handle,psi4pt*gxx(i,j,k),& + psi4pt*gxy(i,j,k),psi4pt*gxz(i,j,k),& + psi4pt*gyy(i,j,k),psi4pt*gyz(i,j,k),psi4pt*gzz(i,j,k),& + det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + tau(i,j,k),rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),& + eps(i,j,k),press(i,j,k),w_lorentz(i,j,k)) + end do end do end do @@ -272,25 +250,8 @@ subroutine Prim2ConservativePolytype(CCTK_ARGUMENTS) gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset)) gzzr = 0.5d0 * (gzz(i,j,k) + gzz(i+xoffset,j+yoffset,k+zoffset)) - if (.not.(conformal_state .eq. 0)) then - psi4r = (0.5d0*(psi(i,j,k)+psi(i+xoffset,j+yoffset,k+zoffset)))**4 - psi4l = (0.5d0*(psi(i,j,k)+psi(i-xoffset,j-yoffset,k-zoffset)))**4 - gxxr = gxxr * psi4r - gxyr = gxyr * psi4r - gxzr = gxzr * psi4r - gyyr = gyyr * psi4r - gyzr = gyzr * psi4r - gzzr = gzzr * psi4r - gxxl = gxxl * psi4l - gxyl = gxyl * psi4l - gxzl = gxzl * psi4l - gyyl = gyyl * psi4l - gyzl = gyzl * psi4l - gzzl = gzzl * psi4l - end if - - call SpatialDeterminant(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl,avg_detl) - call SpatialDeterminant(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr,avg_detr) + avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl) + avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr) call prim2conpolytype(GRHydro_eos_handle, gxxl,gxyl,gxzl,& gyyl,gyzl,gzzl, & @@ -437,8 +398,9 @@ subroutine Primitive2ConservativePolyCells(CCTK_ARGUMENTS) do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1 do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1 - call SpatialDeterminant(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),& - gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),det) + det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),\ + gyy(i,j,k),gyz(i,j,k),gzz(i,j,k)) + call prim2conpolytype(GRHydro_eos_handle,gxx(i,j,k),gxy(i,j,k),& gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),& det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& @@ -494,26 +456,8 @@ subroutine Prim2ConservativeTracer(CCTK_ARGUMENTS) gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset)) gzzr = 0.5d0 * (gzz(i,j,k) + gzz(i+xoffset,j+yoffset,k+zoffset)) - if (.not.(conformal_state .eq. 0)) then - psi4r = (0.5d0*(psi(i,j,k)+psi(i+xoffset,j+yoffset,k+zoffset)))**4 - psi4l = (0.5d0*(psi(i,j,k)+psi(i-xoffset,j-yoffset,k-zoffset)))**4 - gxxr = gxxr * psi4r - gxyr = gxyr * psi4r - gxzr = gxzr * psi4r - gyyr = gyyr * psi4r - gyzr = gyzr * psi4r - gzzr = gzzr * psi4r - gxxl = gxxl * psi4l - gxyl = gxyl * psi4l - gxzl = gxzl * psi4l - gyyl = gyyl * psi4l - gyzl = gyzl * psi4l - gzzl = gzzl * psi4l - end if - - call SpatialDeterminant(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl,avg_detl) - call SpatialDeterminant(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr,avg_detr) - + avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl) + avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr) cons_tracerplus(i,j,k,:) = tracerplus(i,j,k,:) * & sqrt(avg_detr) * rhoplus(i,j,k) / & sqrt(1.d0 - & @@ -591,25 +535,8 @@ subroutine primitive2conservativegeneral(CCTK_ARGUMENTS) gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset)) gzzr = 0.5d0 * (gzz(i,j,k) + gzz(i+xoffset,j+yoffset,k+zoffset)) - if (.not.(conformal_state .eq. 0)) then - psi4r = (0.5d0*(psi(i,j,k)+psi(i+xoffset,j+yoffset,k+zoffset)))**4 - psi4l = (0.5d0*(psi(i,j,k)+psi(i-xoffset,j-yoffset,k-zoffset)))**4 - gxxr = gxxr * psi4r - gxyr = gxyr * psi4r - gxzr = gxzr * psi4r - gyyr = gyyr * psi4r - gyzr = gyzr * psi4r - gzzr = gzzr * psi4r - gxxl = gxxl * psi4l - gxyl = gxyl * psi4l - gxzl = gxzl * psi4l - gyyl = gyyl * psi4l - gyzl = gyzl * psi4l - gzzl = gzzl * psi4l - end if - - call SpatialDeterminant(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl,avg_detl) - call SpatialDeterminant(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr,avg_detr) + avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl) + avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr) call prim2conpolytype(-1, gxxl,gxyl,gxzl,& gyyl,gyzl,gzzl, & @@ -674,17 +601,7 @@ subroutine Primitive2ConservativeCellsGeneral(CCTK_ARGUMENTS) gyzpt = gyz(i,j,k) gzzpt = gzz(i,j,k) - if (.not.(conformal_state .eq. 0)) then - psi4pt = psi(i,j,k)**4 - gxxpt = gxxpt * psi4pt - gxypt = gxypt * psi4pt - gxzpt = gxzpt * psi4pt - gyypt = gyypt * psi4pt - gyzpt = gyzpt * psi4pt - gzzpt = gzzpt * psi4pt - end if - - call SpatialDeterminant(gxxpt,gxypt,gxzpt,gyypt,gyzpt,gzzpt,det) + det = SPATIAL_DETERMINANT(gxxpt,gxypt,gxzpt,gyypt,gyzpt,gzzpt) w_lorentz(i,j,k) = 1.d0 / sqrt(1.d0 - & (gxxpt * velx(i,j,k)**2 + & |