diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/GRHydro_Boundaries.F90 | 27 | ||||
| -rw-r--r-- | src/GRHydro_Con2Prim.F90 | 125 | ||||
| -rw-r--r-- | src/GRHydro_EoSChangeGamma.F90 | 9 | ||||
| -rw-r--r-- | src/GRHydro_FluxSplit.F90 | 154 | ||||
| -rw-r--r-- | src/GRHydro_HLLE.F90 | 80 | ||||
| -rw-r--r-- | src/GRHydro_HLLEPoly.F90 | 83 | ||||
| -rw-r--r-- | src/GRHydro_Marquina.F90 | 31 | ||||
| -rw-r--r-- | src/GRHydro_PPM.F90 | 5 | ||||
| -rw-r--r-- | src/GRHydro_ParamCheck.F90 | 7 | ||||
| -rw-r--r-- | src/GRHydro_Prim2Con.F90 | 131 | ||||
| -rw-r--r-- | src/GRHydro_Reconstruct.F90 | 6 | ||||
| -rw-r--r-- | src/GRHydro_ReconstructPoly.F90 | 6 | ||||
| -rw-r--r-- | src/GRHydro_RiemannSolve.F90 | 15 | ||||
| -rw-r--r-- | src/GRHydro_RoeSolver.F90 | 29 | ||||
| -rw-r--r-- | src/GRHydro_Source.F90 | 174 | ||||
| -rw-r--r-- | src/GRHydro_Tmunu.F90 | 56 | ||||
| -rw-r--r-- | src/GRHydro_UpdateMask.F90 | 30 | ||||
| -rw-r--r-- | src/Utils.F90 | 36 |
18 files changed, 273 insertions, 731 deletions
diff --git a/src/GRHydro_Boundaries.F90 b/src/GRHydro_Boundaries.F90 index 5e8c093..b4dcafb 100644 --- a/src/GRHydro_Boundaries.F90 +++ b/src/GRHydro_Boundaries.F90 @@ -11,6 +11,7 @@ #include "cctk_Parameters.h" #include "cctk_Arguments.h" #include "cctk_Functions.h" +#include "GRHydro_Macros.h" #include "util_Table.h" @@ -246,16 +247,9 @@ subroutine GRHydro_OutflowBoundaries(CCTK_ARGUMENTS) press(i,j,k) = press(i,j,k-1) w_lorentz(i,j,k) = w_lorentz(i,j,k-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 + 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),& @@ -293,16 +287,9 @@ subroutine GRHydro_OutflowBoundaries(CCTK_ARGUMENTS) press(i,j,k) = press(i,j,k+1) w_lorentz(i,j,k) = w_lorentz(i,j,k+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 + 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),& diff --git a/src/GRHydro_Con2Prim.F90 b/src/GRHydro_Con2Prim.F90 index 160493c..24811c5 100644 --- a/src/GRHydro_Con2Prim.F90 +++ b/src/GRHydro_Con2Prim.F90 @@ -13,6 +13,7 @@ #include "cctk_Functions.h" #include "SpaceMask.h" #include "GRHydro_Interfaces.h" +#include "GRHydro_Macros.h" /*@@ @routine Conservative2Primitive @@ -26,7 +27,7 @@ @calledby @history Trimmed and altered from the GR3D routines, original author Mark Miller. - 2007?: Bruno escluded the points in the atmosphere and excision region from the computation. + 2007?: Bruno excluded the points in the atmosphere and excision region from the computation. Aug. 2008: Luca added a check on whether a failure at a given point may be disregarded, because that point will then be restricted from a finer level. This should be completely safe only if *regridding happens at times when all levels are evolved.* @@ -115,21 +116,11 @@ subroutine Conservative2Primitive(CCTK_ARGUMENTS) epsnegative = .false. - if (conformal_state .eq. 0) then - 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) - call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det,& - gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),& - gyz(i,j,k),gzz(i,j,k)) - 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) - call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det,& - 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)) - end if + 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 UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det,& + gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),& + gyz(i,j,k),gzz(i,j,k)) !!$ if (det < 0.e0) then !!$ call CCTK_WARN(0, "nan produced (det negative)") @@ -637,28 +628,12 @@ subroutine Conservative2PrimitiveBounds(CCTK_ARGUMENTS) epsnegative = .false. - if (conformal_state .eq. 0) then - call SpatialDeterminant(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl,avg_detl) - call SpatialDeterminant(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr,avg_detr) - call UpperMetric(uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl,& - gxxl,gxyl,gxzl,gyyl,gyzl,gzzl) - call UpperMetric(uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr,& - gxxr,gxyr,gxzr,gyyr,gyzr,gzzr) - else - 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 - - call SpatialDeterminant(psi4l*gxxl,psi4l*gxyl,psi4l*gxzl,& - psi4l*gyyl,psi4l*gyzl,psi4l*gzzl,avg_detl) - call SpatialDeterminant(psi4r*gxxr,psi4r*gxyr,psi4r*gxzr,& - psi4r*gyyr,psi4r*gyzr,psi4r*gzzr,avg_detr) - call UpperMetric(uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl,& - psi4l*gxxl,psi4l*gxyl,psi4l*gxzl,psi4l*gyyl,psi4l*gyzl,& - psi4l*gzzl) - call UpperMetric(uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr,& - psi4r*gxxr,psi4r*gxyr,psi4r*gxzr,psi4r*gyyr,psi4r*gyzr,& - psi4r*gzzr) - end if + avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl) + avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr) + call UpperMetric(uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl,& + gxxl,gxyl,gxzl,gyyl,gyzl,gzzl) + call UpperMetric(uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr,& + gxxr,gxyr,gxzr,gyyr,gyzr,gzzr) if (evolve_tracer .ne. 0) then do itracer=1,number_of_tracers @@ -837,21 +812,11 @@ subroutine Conservative2PrimitivePolytype(CCTK_ARGUMENTS) if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. & GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0)) cycle - if (conformal_state .eq. 0) then - 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) - call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det,& - gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),& - gyz(i,j,k),gzz(i,j,k)) - 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) - call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det,& - 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)) - end if + 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 UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det,& + gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),& + gyz(i,j,k),gzz(i,j,k)) if (evolve_tracer .ne. 0) then do itracer=1,number_of_tracers @@ -1252,28 +1217,12 @@ subroutine Con2PrimBoundsPolytype(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 (conformal_state .eq. 0) then - call SpatialDeterminant(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl,avg_detl) - call SpatialDeterminant(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr,avg_detr) - call UpperMetric(uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl,& - gxxl,gxyl,gxzl,gyyl,gyzl,gzzl) - call UpperMetric(uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr,& - gxxr,gxyr,gxzr,gyyr,gyzr,gzzr) - else - 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 - - call SpatialDeterminant(psi4l*gxxl,psi4l*gxyl,psi4l*gxzl,& - psi4l*gyyl,psi4l*gyzl,psi4l*gzzl,avg_detl) - call SpatialDeterminant(psi4r*gxxr,psi4r*gxyr,psi4r*gxzr,& - psi4r*gyyr,psi4r*gyzr,psi4r*gzzr,avg_detr) - call UpperMetric(uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl,& - psi4l*gxxl,psi4l*gxyl,psi4l*gxzl,psi4l*gyyl,psi4l*gyzl,& - psi4l*gzzl) - call UpperMetric(uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr,& - psi4r*gxxr,psi4r*gxyr,psi4r*gxzr,psi4r*gyyr,psi4r*gyzr,& - psi4r*gzzr) - end if + avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl) + avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr) + call UpperMetric(uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl,& + gxxl,gxyl,gxzl,gyyl,gyzl,gzzl) + call UpperMetric(uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr,& + gxxr,gxyr,gxzr,gyyr,gyzr,gzzr) call Con2Prim_ptPolytype(GRHydro_eos_handle, densminus(i,j,k),& sxminus(i,j,k),syminus(i,j,k),szminus(i,j,k),& @@ -1350,28 +1299,12 @@ subroutine Con2PrimBoundsTracer(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 (conformal_state .eq. 0) then - call SpatialDeterminant(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl,avg_detl) - call SpatialDeterminant(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr,avg_detr) - call UpperMetric(uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl,& - gxxl,gxyl,gxzl,gyyl,gyzl,gzzl) - call UpperMetric(uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr,& + avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl) + avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr) + call UpperMetric(uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl,& + gxxl,gxyl,gxzl,gyyl,gyzl,gzzl) + call UpperMetric(uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr,& gxxr,gxyr,gxzr,gyyr,gyzr,gzzr) - else - 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 - - call SpatialDeterminant(psi4l*gxxl,psi4l*gxyl,psi4l*gxzl,& - psi4l*gyyl,psi4l*gyzl,psi4l*gzzl,avg_detl) - call SpatialDeterminant(psi4r*gxxr,psi4r*gxyr,psi4r*gxzr,& - psi4r*gyyr,psi4r*gyzr,psi4r*gzzr,avg_detr) - call UpperMetric(uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl,& - psi4l*gxxl,psi4l*gxyl,psi4l*gxzl,psi4l*gyyl,psi4l*gyzl,& - psi4l*gzzl) - call UpperMetric(uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr,& - psi4r*gxxr,psi4r*gxyr,psi4r*gxzr,psi4r*gyyr,psi4r*gyzr,& - psi4r*gzzr) - end if do itracer=1,number_of_tracers call Con2Prim_ptBoundsTracer(cons_tracerplus(i,j,k,itracer), & diff --git a/src/GRHydro_EoSChangeGamma.F90 b/src/GRHydro_EoSChangeGamma.F90 index c075aae..6a502f4 100644 --- a/src/GRHydro_EoSChangeGamma.F90 +++ b/src/GRHydro_EoSChangeGamma.F90 @@ -12,6 +12,7 @@ #include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.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) @@ -97,8 +98,8 @@ subroutine GRHydro_EoSChangeGamma(CCTK_ARGUMENTS) do j = 1, cctk_lsh(2) do i = 1, cctk_lsh(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_polytrope_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),& @@ -313,8 +314,8 @@ subroutine GRHydro_EoSChangeGammaK_Shibata(CCTK_ARGUMENTS) do j = 1, cctk_lsh(2) do i = 1, cctk_lsh(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 prim2con(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),& diff --git a/src/GRHydro_FluxSplit.F90 b/src/GRHydro_FluxSplit.F90 index c5fe6d7..d54e446 100644 --- a/src/GRHydro_FluxSplit.F90 +++ b/src/GRHydro_FluxSplit.F90 @@ -11,6 +11,7 @@ #include "cctk_Arguments.h" #include "cctk_Parameters.h" #include "cctk_Functions.h" +#include "GRHydro_Macros.h" /*@@ @routine GRHydro_FSAlpha @@ -66,22 +67,11 @@ subroutine GRHydro_FSAlpha(CCTK_ARGUMENTS) do j = 1, ny do i = 1, nx - if (conformal_state .eq. 0) then - 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) - call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det,& - gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),& - gyz(i,j,k),gzz(i,j,k)) - 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) - call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det,& - 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)) - end if + 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 UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det,& + gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),& + gyz(i,j,k),gzz(i,j,k)) if (shift_state .eq. 0) then beta = 0.d0 @@ -123,22 +113,11 @@ subroutine GRHydro_FSAlpha(CCTK_ARGUMENTS) do j = 1, ny do i = 1, nx - if (conformal_state .eq. 0) then - 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) - call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det,& - gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),& - gyz(i,j,k),gzz(i,j,k)) - 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) - call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det,& - 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)) - end if + 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 UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det,& + gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),& + gyz(i,j,k),gzz(i,j,k)) if (shift_state .eq. 0) then beta = 0.d0 @@ -180,22 +159,11 @@ subroutine GRHydro_FSAlpha(CCTK_ARGUMENTS) do j = 1, ny do i = 1, nx - if (conformal_state .eq. 0) then - 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) - call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det,& - gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),& - gyz(i,j,k),gzz(i,j,k)) - 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) - call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det,& - 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)) - end if + 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 UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det,& + gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),& + gyz(i,j,k),gzz(i,j,k)) if (shift_state .eq. 0) then beta = 0.d0 @@ -347,28 +315,14 @@ subroutine GRHydro_SplitFlux(CCTK_ARGUMENTS) else dummy = betax(:,j,k) end if - if (conformal_state .eq. 0) then - do i = 1, cctk_lsh(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(i)) - call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det(i),& - gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),& - gyz(i,j,k),gzz(i,j,k)) - upper(i) = uxx - end do - else - do i = 1, cctk_lsh(1) - 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(i)) - call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det(i),& - 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)) - upper(i) = uxx - end do - end if + do i = 1, cctk_lsh(1) + 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 UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det(i),& + gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),& + gyz(i,j,k),gzz(i,j,k)) + upper(i) = uxx + end do call GRHydro_SplitFlux_1D(GRHydro_eos_handle, nx,& fs_alpha1, fs_alpha2, fs_alpha3, fs_alpha4, fs_alpha5, & @@ -398,28 +352,14 @@ subroutine GRHydro_SplitFlux(CCTK_ARGUMENTS) else dummy = betay(i,:,k) end if - if (conformal_state .eq. 0) then - do j = 1, cctk_lsh(2) - 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(j)) - call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det(j),& - gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),& - gyz(i,j,k),gzz(i,j,k)) - upper(j) = uyy - end do - else - do j = 1, cctk_lsh(2) - 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(j)) - call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det(j),& - 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)) - upper(j) = uyy - end do - end if + do j = 1, cctk_lsh(2) + 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 UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det(j),& + gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),& + gyz(i,j,k),gzz(i,j,k)) + upper(j) = uyy + end do call GRHydro_SplitFlux_1D(GRHydro_eos_handle, ny,& fs_alpha1, fs_alpha2, fs_alpha3, fs_alpha4, fs_alpha5, & @@ -450,28 +390,14 @@ subroutine GRHydro_SplitFlux(CCTK_ARGUMENTS) dummy = betaz(i,j,:) end if - if (conformal_state .eq. 0) then - do k = 1, cctk_lsh(3) - 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(k)) - call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det(k),& - gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),& - gyz(i,j,k),gzz(i,j,k)) - upper(k) = uzz - end do - else - do k = 1, cctk_lsh(3) - 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(k)) - call UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det(k),& - 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)) - upper(k) = uzz - end do - end if + do k = 1, cctk_lsh(3) + 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 UpperMetric(uxx,uxy,uxz,uyy,uyz,uzz,det(k),& + gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),& + gyz(i,j,k),gzz(i,j,k)) + upper(k) = uzz + end do call GRHydro_SplitFlux_1D(GRHydro_eos_handle, nz,& fs_alpha1, fs_alpha2, fs_alpha3, fs_alpha4, fs_alpha5, & 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,:) 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,:) diff --git a/src/GRHydro_Marquina.F90 b/src/GRHydro_Marquina.F90 index 9d1e946..f3c2030 100644 --- a/src/GRHydro_Marquina.F90 +++ b/src/GRHydro_Marquina.F90 @@ -15,6 +15,7 @@ #include "cctk_Parameters.h" #include "cctk_Functions.h" +#include "GRHydro_Macros.h" #include "SpaceMask.h" /*@@ @@ -139,21 +140,9 @@ subroutine GRHydro_Marquina(CCTK_ARGUMENTS) gzzh = 0.5d0 * (gzz(i+xoffset,j+yoffset,k+zoffset) + & gzz(i,j,k)) - if (conformal_state > 0) then - - psi4h = (0.5d0*(psi(i,j,k)+psi(i+xoffset,j+yoffset,k+zoffset)))**4 - gxxh = gxxh * psi4h - gxyh = gxyh * psi4h - gxzh = gxzh * psi4h - gyyh = gyyh * psi4h - gyzh = gyzh * psi4h - gzzh = gzzh * psi4h - - end if - !!$ routine to calculate the determinant of the metric - call SpatialDeterminant(gxxh,gxyh,gxzh,gyyh,gyzh,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 @@ -491,22 +480,10 @@ subroutine GRHydro_MarquinaGeneral(CCTK_ARGUMENTS) gzzh = 0.5d0 * (gzz(i+xoffset,j+yoffset,k+zoffset) + & gzz(i,j,k)) - if (conformal_state > 0) then - - psi4h = (0.5d0*(psi(i,j,k)+psi(i+xoffset,j+yoffset,k+zoffset)))**4 - gxxh = gxxh * psi4h - gxyh = gxyh * psi4h - gxzh = gxzh * psi4h - gyyh = gyyh * psi4h - gyzh = gyzh * psi4h - gzzh = gzzh * psi4h - - end if - !!$ routine to calculate the determinant of the metric - call SpatialDeterminant(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,avg_det) - + avg_det = SPATIAL_DETERMINANT(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh) + !!$ If the Riemann problem is trivial the flux is already correct if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, trivial)) then diff --git a/src/GRHydro_PPM.F90 b/src/GRHydro_PPM.F90 index cc42787..ea1607f 100644 --- a/src/GRHydro_PPM.F90 +++ b/src/GRHydro_PPM.F90 @@ -10,6 +10,7 @@ #include "cctk.h" #include "cctk_Arguments.h" #include "cctk_Parameters.h" +#include "GRHydro_Macros.h" subroutine PPM_TVD(origm, orig, origp, bextm, bextp) CCTK_REAL :: origm, orig, origp, bextm, bextp @@ -217,8 +218,8 @@ trivial_rp = .true. agyy = 0.5d0*( psi4(i)*gyy(i) + psi4(i+1)*gyy(i+1) ) agyz = 0.5d0*( psi4(i)*gyz(i) + psi4(i+1)*gyz(i+1) ) agzz = 0.5d0*( psi4(i)*gzz(i) + psi4(i+1)*gzz(i+1) ) - call SpatialDeterminant(agxx, agxy, agxz,& - agyy, agyz, agzz, det) + det = SPATIAL_DETERMINANT(agxx, agxy, agxz, \ + agyy, agyz, agzz) call UpperMetric (uxx, uxy, uxz, uyy, uyz, uzz, & det, agxx, agxy, agxz, agyy, agyz, agzz) call eigenvalues(handle,& diff --git a/src/GRHydro_ParamCheck.F90 b/src/GRHydro_ParamCheck.F90 index f311659..f192071 100644 --- a/src/GRHydro_ParamCheck.F90 +++ b/src/GRHydro_ParamCheck.F90 @@ -67,12 +67,7 @@ subroutine GRHydro_ParamCheck(CCTK_ARGUMENTS) if (.not.(CCTK_EQUALS(metric_type, "Physical").or.& CCTK_EQUALS(metric_type, "Static Conformal"))) then - call CCTK_PARAMWARN("GRHydro only knows how to deal with physical metric type at the moment. Complain to the maintainers...") - end if - - if (.not.CCTK_EQUALS(metric_type, "Static Conformal")) then - call CCTK_INFO("GRHydro will use the physical metric and set conformal_state=0.") - conformal_state = 0 + call CCTK_WARN(0,"GRHydro only knows how to deal with physical metric type at the moment. Complain to the maintainers...") end if 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 + & diff --git a/src/GRHydro_Reconstruct.F90 b/src/GRHydro_Reconstruct.F90 index d7af59a..3d4b843 100644 --- a/src/GRHydro_Reconstruct.F90 +++ b/src/GRHydro_Reconstruct.F90 @@ -102,11 +102,7 @@ subroutine Reconstruction(CCTK_ARGUMENTS) !!$ Currently only option is reconstruction on primitive variables. !!$ Should change this. - if (conformal_state .eq. 0) then - psi4 = 1.d0 - else - psi4 = psi**4 - end if + psi4 = 1.d0 if (shift_state .ne. 0) then lbetax = betax lbetay = betay diff --git a/src/GRHydro_ReconstructPoly.F90 b/src/GRHydro_ReconstructPoly.F90 index 0467c27..26b2acd 100644 --- a/src/GRHydro_ReconstructPoly.F90 +++ b/src/GRHydro_ReconstructPoly.F90 @@ -102,11 +102,7 @@ subroutine ReconstructionPolytype(CCTK_ARGUMENTS) !!$ Currently only option is reconstruction on primitive variables. !!$ Should change this. - if (conformal_state .eq. 0) then - psi4 = 1.d0 - else - psi4 = psi**4 - end if + psi4 = 1.0d0 if (shift_state .ne. 0) then lbetax = betax lbetay = betay diff --git a/src/GRHydro_RiemannSolve.F90 b/src/GRHydro_RiemannSolve.F90 index 76739a3..dd01bff 100644 --- a/src/GRHydro_RiemannSolve.F90 +++ b/src/GRHydro_RiemannSolve.F90 @@ -11,6 +11,7 @@ #include "cctk_Parameters.h" #include "cctk_Arguments.h" #include "cctk_Functions.h" +#include "GRHydro_Macros.h" /*@@ @routine RiemannSolve @@ -251,19 +252,7 @@ subroutine RiemannSolveGeneral(CCTK_ARGUMENTS) gzzh = 0.5d0 * (gzz(i+xoffset,j+yoffset,k+zoffset) + & gzz(i,j,k)) - if (conformal_state > 0) then - - psi4h = (0.5d0*(psi(i,j,k)+psi(i+xoffset,j+yoffset,k+zoffset)))**4 - gxxh = gxxh * psi4h - gxyh = gxyh * psi4h - gxzh = gxzh * psi4h - gyyh = gyyh * psi4h - gyzh = gyzh * psi4h - gzzh = gzzh * psi4h - - end if - - call SpatialDeterminant(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,avg_det) + avg_det = SPATIAL_DETERMINANT(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh) if (flux_direction == 1) then diff --git a/src/GRHydro_RoeSolver.F90 b/src/GRHydro_RoeSolver.F90 index 7d77599..329c081 100644 --- a/src/GRHydro_RoeSolver.F90 +++ b/src/GRHydro_RoeSolver.F90 @@ -11,6 +11,7 @@ #include "cctk_Arguments.h" #include "cctk_Parameters.h" #include "cctk_Functions.h" +#include "GRHydro_Macros.h" #include "SpaceMask.h" @@ -138,19 +139,7 @@ subroutine GRHydro_RoeSolve(CCTK_ARGUMENTS) gzzh = 0.5d0 * (gzz(i+xoffset,j+yoffset,k+zoffset) + & gzz(i,j,k)) - if (conformal_state > 0) then - - psi4h = (0.5d0*(psi(i,j,k)+psi(i+xoffset,j+yoffset,k+zoffset)))**4 - gxxh = gxxh * psi4h - gxyh = gxyh * psi4h - gxzh = gxzh * psi4h - gyyh = gyyh * psi4h - gyzh = gyzh * psi4h - gzzh = gzzh * psi4h - - end if - - call SpatialDeterminant(gxxh,gxyh,gxzh,gyyh,gyzh,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 @@ -493,20 +482,8 @@ subroutine GRHydro_RoeSolveGeneral(CCTK_ARGUMENTS) gyz(i,j,k)) gzzh = 0.5d0 * (gzz(i+xoffset,j+yoffset,k+zoffset) + & gzz(i,j,k)) - - if (conformal_state > 0) then - - psi4h = (0.5d0*(psi(i,j,k)+psi(i+xoffset,j+yoffset,k+zoffset)))**4 - gxxh = gxxh * psi4h - gxyh = gxyh * psi4h - gxzh = gxzh * psi4h - gyyh = gyyh * psi4h - gyzh = gyzh * psi4h - gzzh = gzzh * psi4h - - end if - call SpatialDeterminant(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh,avg_det) + avg_det = SPATIAL_DETERMINANT(gxxh,gxyh,gxzh,gyyh,gyzh,gzzh) !!$ If the Riemann problem is trivial the flux is already correct diff --git a/src/GRHydro_Source.F90 b/src/GRHydro_Source.F90 index b2833a8..e8c00ec 100644 --- a/src/GRHydro_Source.F90 +++ b/src/GRHydro_Source.F90 @@ -25,6 +25,7 @@ #include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.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) @@ -153,31 +154,15 @@ subroutine SourceTerms(CCTK_ARGUMENTS) !!$ Set the metric terms. - if (conformal_state .eq. 0) then + localgxx = gxx(i,j,k) + localgxy = gxy(i,j,k) + localgxz = gxz(i,j,k) + localgyy = gyy(i,j,k) + localgyz = gyz(i,j,k) + localgzz = gzz(i,j,k) - localgxx = gxx(i,j,k) - localgxy = gxy(i,j,k) - localgxz = gxz(i,j,k) - localgyy = gyy(i,j,k) - localgyz = gyz(i,j,k) - localgzz = gzz(i,j,k) - - else - - psi4pt = psi(i,j,k)**4 - - localgxx = psi4pt*gxx(i,j,k) - localgxy = psi4pt*gxy(i,j,k) - localgxz = psi4pt*gxz(i,j,k) - localgyy = psi4pt*gyy(i,j,k) - localgyz = psi4pt*gyz(i,j,k) - localgzz = psi4pt*gzz(i,j,k) - - end if - - - call SpatialDeterminant(localgxx, localgxy, localgxz,& - localgyy, localgyz, localgzz, det) + det = SPATIAL_DETERMINANT(localgxx, localgxy, localgxz,\ + localgyy, localgyz, localgzz) sqrtdet = sqrt(det) call UpperMetric(uxx, uxy, uxz, uyy, uyz, uzz, det, localgxx,& localgxy, localgxz, localgyy, localgyz, localgzz) @@ -294,111 +279,50 @@ subroutine SourceTerms(CCTK_ARGUMENTS) end if - if (conformal_state .eq. 0) then - - if (local_spatial_order .eq. 2) then - - dx_gxx = DIFF_X_2(gxx) - dx_gxy = DIFF_X_2(gxy) - dx_gxz = DIFF_X_2(gxz) - dx_gyy = DIFF_X_2(gyy) - dx_gyz = DIFF_X_2(gyz) - dx_gzz = DIFF_X_2(gzz) - dy_gxx = DIFF_Y_2(gxx) - dy_gxy = DIFF_Y_2(gxy) - dy_gxz = DIFF_Y_2(gxz) - dy_gyy = DIFF_Y_2(gyy) - dy_gyz = DIFF_Y_2(gyz) - dy_gzz = DIFF_Y_2(gzz) - dz_gxx = DIFF_Z_2(gxx) - dz_gxy = DIFF_Z_2(gxy) - dz_gxz = DIFF_Z_2(gxz) - dz_gyy = DIFF_Z_2(gyy) - dz_gyz = DIFF_Z_2(gyz) - dz_gzz = DIFF_Z_2(gzz) - - else - - dx_gxx = DIFF_X_4(gxx) - dx_gxy = DIFF_X_4(gxy) - dx_gxz = DIFF_X_4(gxz) - dx_gyy = DIFF_X_4(gyy) - dx_gyz = DIFF_X_4(gyz) - dx_gzz = DIFF_X_4(gzz) - dy_gxx = DIFF_Y_4(gxx) - dy_gxy = DIFF_Y_4(gxy) - dy_gxz = DIFF_Y_4(gxz) - dy_gyy = DIFF_Y_4(gyy) - dy_gyz = DIFF_Y_4(gyz) - dy_gzz = DIFF_Y_4(gzz) - dz_gxx = DIFF_Z_4(gxx) - dz_gxy = DIFF_Z_4(gxy) - dz_gxz = DIFF_Z_4(gxz) - dz_gyy = DIFF_Z_4(gyy) - dz_gyz = DIFF_Z_4(gyz) - dz_gzz = DIFF_Z_4(gzz) - - end if - - else if (conformal_state > 1) then - -!!$ Using conformal factor, and the derivatives are -!!$ already calculated. - - psi4pt = psi(i,j,k)**4 - dx_psi4 = 4.d0*psi4pt*psix(i,j,k) - dy_psi4 = 4.d0*psi4pt*psiy(i,j,k) - dz_psi4 = 4.d0*psi4pt*psiz(i,j,k) - - if (local_spatial_order .eq. 2) then - - dx_gxx = DIFF_X_2(gxx)*psi4pt + dx_psi4*gxx(i,j,k) - dx_gxy = DIFF_X_2(gxy)*psi4pt + dx_psi4*gxy(i,j,k) - dx_gxz = DIFF_X_2(gxz)*psi4pt + dx_psi4*gxz(i,j,k) - dx_gyy = DIFF_X_2(gyy)*psi4pt + dx_psi4*gyy(i,j,k) - dx_gyz = DIFF_X_2(gyz)*psi4pt + dx_psi4*gyz(i,j,k) - dx_gzz = DIFF_X_2(gzz)*psi4pt + dx_psi4*gzz(i,j,k) - dy_gxx = DIFF_Y_2(gxx)*psi4pt + dy_psi4*gxx(i,j,k) - dy_gxy = DIFF_Y_2(gxy)*psi4pt + dy_psi4*gxy(i,j,k) - dy_gxz = DIFF_Y_2(gxz)*psi4pt + dy_psi4*gxz(i,j,k) - dy_gyy = DIFF_Y_2(gyy)*psi4pt + dy_psi4*gyy(i,j,k) - dy_gyz = DIFF_Y_2(gyz)*psi4pt + dy_psi4*gyz(i,j,k) - dy_gzz = DIFF_Y_2(gzz)*psi4pt + dy_psi4*gzz(i,j,k) - dz_gxx = DIFF_Z_2(gxx)*psi4pt + dz_psi4*gxx(i,j,k) - dz_gxy = DIFF_Z_2(gxy)*psi4pt + dz_psi4*gxy(i,j,k) - dz_gxz = DIFF_Z_2(gxz)*psi4pt + dz_psi4*gxz(i,j,k) - dz_gyy = DIFF_Z_2(gyy)*psi4pt + dz_psi4*gyy(i,j,k) - dz_gyz = DIFF_Z_2(gyz)*psi4pt + dz_psi4*gyz(i,j,k) - dz_gzz = DIFF_Z_2(gzz)*psi4pt + dz_psi4*gzz(i,j,k) + if (local_spatial_order .eq. 2) then - else + dx_gxx = DIFF_X_2(gxx) + dx_gxy = DIFF_X_2(gxy) + dx_gxz = DIFF_X_2(gxz) + dx_gyy = DIFF_X_2(gyy) + dx_gyz = DIFF_X_2(gyz) + dx_gzz = DIFF_X_2(gzz) + dy_gxx = DIFF_Y_2(gxx) + dy_gxy = DIFF_Y_2(gxy) + dy_gxz = DIFF_Y_2(gxz) + dy_gyy = DIFF_Y_2(gyy) + dy_gyz = DIFF_Y_2(gyz) + dy_gzz = DIFF_Y_2(gzz) + dz_gxx = DIFF_Z_2(gxx) + dz_gxy = DIFF_Z_2(gxy) + dz_gxz = DIFF_Z_2(gxz) + dz_gyy = DIFF_Z_2(gyy) + dz_gyz = DIFF_Z_2(gyz) + dz_gzz = DIFF_Z_2(gzz) + + else - dx_gxx = DIFF_X_4(gxx)*psi4pt + dx_psi4*gxx(i,j,k) - dx_gxy = DIFF_X_4(gxy)*psi4pt + dx_psi4*gxy(i,j,k) - dx_gxz = DIFF_X_4(gxz)*psi4pt + dx_psi4*gxz(i,j,k) - dx_gyy = DIFF_X_4(gyy)*psi4pt + dx_psi4*gyy(i,j,k) - dx_gyz = DIFF_X_4(gyz)*psi4pt + dx_psi4*gyz(i,j,k) - dx_gzz = DIFF_X_4(gzz)*psi4pt + dx_psi4*gzz(i,j,k) - dy_gxx = DIFF_Y_4(gxx)*psi4pt + dy_psi4*gxx(i,j,k) - dy_gxy = DIFF_Y_4(gxy)*psi4pt + dy_psi4*gxy(i,j,k) - dy_gxz = DIFF_Y_4(gxz)*psi4pt + dy_psi4*gxz(i,j,k) - dy_gyy = DIFF_Y_4(gyy)*psi4pt + dy_psi4*gyy(i,j,k) - dy_gyz = DIFF_Y_4(gyz)*psi4pt + dy_psi4*gyz(i,j,k) - dy_gzz = DIFF_Y_4(gzz)*psi4pt + dy_psi4*gzz(i,j,k) - dz_gxx = DIFF_Z_4(gxx)*psi4pt + dz_psi4*gxx(i,j,k) - dz_gxy = DIFF_Z_4(gxy)*psi4pt + dz_psi4*gxy(i,j,k) - dz_gxz = DIFF_Z_4(gxz)*psi4pt + dz_psi4*gxz(i,j,k) - dz_gyy = DIFF_Z_4(gyy)*psi4pt + dz_psi4*gyy(i,j,k) - dz_gyz = DIFF_Z_4(gyz)*psi4pt + dz_psi4*gyz(i,j,k) - dz_gzz = DIFF_Z_4(gzz)*psi4pt + dz_psi4*gzz(i,j,k) + dx_gxx = DIFF_X_4(gxx) + dx_gxy = DIFF_X_4(gxy) + dx_gxz = DIFF_X_4(gxz) + dx_gyy = DIFF_X_4(gyy) + dx_gyz = DIFF_X_4(gyz) + dx_gzz = DIFF_X_4(gzz) + dy_gxx = DIFF_Y_4(gxx) + dy_gxy = DIFF_Y_4(gxy) + dy_gxz = DIFF_Y_4(gxz) + dy_gyy = DIFF_Y_4(gyy) + dy_gyz = DIFF_Y_4(gyz) + dy_gzz = DIFF_Y_4(gzz) + dz_gxx = DIFF_Z_4(gxx) + dz_gxy = DIFF_Z_4(gxy) + dz_gxz = DIFF_Z_4(gxz) + dz_gyy = DIFF_Z_4(gyy) + dz_gyz = DIFF_Z_4(gyz) + dz_gzz = DIFF_Z_4(gzz) - end if - - else -!!$ Have to finite difference the conformal factor - call CCTK_WARN(0,"At this point in the source terms we should be differencing the conformal factor. Ian has been too lazy to put this bit of code in yet.") end if - + !!$ Contract the shift with the extrinsic curvature shiftshiftk = shiftx*shiftx*kxx(i,j,k) + & diff --git a/src/GRHydro_Tmunu.F90 b/src/GRHydro_Tmunu.F90 index 6ba34e3..7939a07 100644 --- a/src/GRHydro_Tmunu.F90 +++ b/src/GRHydro_Tmunu.F90 @@ -86,32 +86,20 @@ do j = 1, cctk_lsh(2) do i = 1, cctk_lsh(1) - if (conformal_state .eq. 0) then - velxlow = gxx(i,j,k)*velx(i,j,k) + gxy(i,j,k)*vely(i,j,k) + gxz(i,j,k)*velz(i,j,k) - velylow = gxy(i,j,k)*velx(i,j,k) + gyy(i,j,k)*vely(i,j,k) + gyz(i,j,k)*velz(i,j,k) - velzlow = gxz(i,j,k)*velx(i,j,k) + gyz(i,j,k)*vely(i,j,k) + gzz(i,j,k)*velz(i,j,k) - else - velxlow = psi(i,j,k)**4 * (gxx(i,j,k)*velx(i,j,k) + gxy(i,j,k)*vely(i,j,k) + gxz(i,j,k)*velz(i,j,k) ) - velylow = psi(i,j,k)**4 * (gxy(i,j,k)*velx(i,j,k) + gyy(i,j,k)*vely(i,j,k) + gyz(i,j,k)*velz(i,j,k) ) - velzlow = psi(i,j,k)**4 * (gxz(i,j,k)*velx(i,j,k) + gyz(i,j,k)*vely(i,j,k) + gzz(i,j,k)*velz(i,j,k) ) - end if + velxlow = gxx(i,j,k)*velx(i,j,k) + gxy(i,j,k)*vely(i,j,k) + gxz(i,j,k)*velz(i,j,k) + velylow = gxy(i,j,k)*velx(i,j,k) + gyy(i,j,k)*vely(i,j,k) + gyz(i,j,k)*velz(i,j,k) + velzlow = gxz(i,j,k)*velx(i,j,k) + gyz(i,j,k)*vely(i,j,k) + gzz(i,j,k)*velz(i,j,k) !!$ Calculate lower components and square of shift vector. if (shift_state .ne. 0) then - if (conformal_state .eq. 0) then - betaxlow = gxx(i,j,k)*betax(i,j,k) + gxy(i,j,k)*betay(i,j,k) + gxz(i,j,k)*betaz(i,j,k) - betaylow = gxy(i,j,k)*betax(i,j,k) + gyy(i,j,k)*betay(i,j,k) + gyz(i,j,k)*betaz(i,j,k) - betazlow = gxz(i,j,k)*betax(i,j,k) + gyz(i,j,k)*betay(i,j,k) + gzz(i,j,k)*betaz(i,j,k) - else - betaxlow = psi(i,j,k)**4 * (gxx(i,j,k)*betax(i,j,k) + gxy(i,j,k)*betay(i,j,k) + gxz(i,j,k)*betaz(i,j,k) ) - betaylow = psi(i,j,k)**4 * (gxy(i,j,k)*betax(i,j,k) + gyy(i,j,k)*betay(i,j,k) + gyz(i,j,k)*betaz(i,j,k) ) - betazlow = psi(i,j,k)**4 * (gxz(i,j,k)*betax(i,j,k) + gyz(i,j,k)*betay(i,j,k) + gzz(i,j,k)*betaz(i,j,k) ) - end if - - beta2 = betax(i,j,k)*betaxlow + betay(i,j,k)*betaylow + betaz(i,j,k)*betazlow + betaxlow = gxx(i,j,k)*betax(i,j,k) + gxy(i,j,k)*betay(i,j,k) + gxz(i,j,k)*betaz(i,j,k) + betaylow = gxy(i,j,k)*betax(i,j,k) + gyy(i,j,k)*betay(i,j,k) + gyz(i,j,k)*betaz(i,j,k) + betazlow = gxz(i,j,k)*betax(i,j,k) + gyz(i,j,k)*betay(i,j,k) + gzz(i,j,k)*betaz(i,j,k) + beta2 = betax(i,j,k)*betaxlow + betay(i,j,k)*betaylow + betaz(i,j,k)*betazlow + else betaxlow = 0.0D0 @@ -143,27 +131,13 @@ eTty(i,j,k) = eTty(i,j,k) + rhoenthalpy*utlow*uylow + press(i,j,k)*betaylow eTtz(i,j,k) = eTtz(i,j,k) + rhoenthalpy*utlow*uzlow + press(i,j,k)*betazlow - if (conformal_state .eq. 0) then - - eTxx(i,j,k) = eTxx(i,j,k) + rhoenthalpy*uxlow**2 + press(i,j,k)*gxx(i,j,k) - eTyy(i,j,k) = eTyy(i,j,k) + rhoenthalpy*uylow**2 + press(i,j,k)*gyy(i,j,k) - eTzz(i,j,k) = eTzz(i,j,k) + rhoenthalpy*uzlow**2 + press(i,j,k)*gzz(i,j,k) - - eTxy(i,j,k) = eTxy(i,j,k) + rhoenthalpy*uxlow*uylow + press(i,j,k)*gxy(i,j,k) - eTxz(i,j,k) = eTxz(i,j,k) + rhoenthalpy*uxlow*uzlow + press(i,j,k)*gxz(i,j,k) - eTyz(i,j,k) = eTyz(i,j,k) + rhoenthalpy*uylow*uzlow + press(i,j,k)*gyz(i,j,k) - - else - - eTxx(i,j,k) = eTxx(i,j,k) + rhoenthalpy*uxlow**2 + psi(i,j,k)**4 * press(i,j,k)*gxx(i,j,k) - eTyy(i,j,k) = eTyy(i,j,k) + rhoenthalpy*uylow**2 + psi(i,j,k)**4 * press(i,j,k)*gyy(i,j,k) - eTzz(i,j,k) = eTzz(i,j,k) + rhoenthalpy*uzlow**2 + psi(i,j,k)**4 * press(i,j,k)*gzz(i,j,k) - - eTxy(i,j,k) = eTxy(i,j,k) + rhoenthalpy*uxlow*uylow + psi(i,j,k)**4 * press(i,j,k)*gxy(i,j,k) - eTxz(i,j,k) = eTxz(i,j,k) + rhoenthalpy*uxlow*uzlow + psi(i,j,k)**4 * press(i,j,k)*gxz(i,j,k) - eTyz(i,j,k) = eTyz(i,j,k) + rhoenthalpy*uylow*uzlow + psi(i,j,k)**4 * press(i,j,k)*gyz(i,j,k) - - end if + eTxx(i,j,k) = eTxx(i,j,k) + rhoenthalpy*uxlow**2 + press(i,j,k)*gxx(i,j,k) + eTyy(i,j,k) = eTyy(i,j,k) + rhoenthalpy*uylow**2 + press(i,j,k)*gyy(i,j,k) + eTzz(i,j,k) = eTzz(i,j,k) + rhoenthalpy*uzlow**2 + press(i,j,k)*gzz(i,j,k) + + eTxy(i,j,k) = eTxy(i,j,k) + rhoenthalpy*uxlow*uylow + press(i,j,k)*gxy(i,j,k) + eTxz(i,j,k) = eTxz(i,j,k) + rhoenthalpy*uxlow*uzlow + press(i,j,k)*gxz(i,j,k) + eTyz(i,j,k) = eTyz(i,j,k) + rhoenthalpy*uylow*uzlow + press(i,j,k)*gyz(i,j,k) end do end do diff --git a/src/GRHydro_UpdateMask.F90 b/src/GRHydro_UpdateMask.F90 index 3753797..d7cdff0 100644 --- a/src/GRHydro_UpdateMask.F90 +++ b/src/GRHydro_UpdateMask.F90 @@ -11,6 +11,7 @@ #include "cctk_Parameters.h" #include "cctk_Arguments.h" +#include "GRHydro_Macros.h" #include "SpaceMask.h" #define velx(i,j,k) vel(i,j,k,1) @@ -209,27 +210,14 @@ subroutine GRHydro_AtmosphereReset(CCTK_ARGUMENTS) velx(i,j,k) = 0.0d0 vely(i,j,k) = 0.0d0 velz(i,j,k) = 0.0d0 - if (conformal_state .eq. 0) then - 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) - call prim2conpolytype(GRHydro_polytrope_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), scon(i,j,k,1), scon(i,j,k,2), scon(i,j,k,3), & - 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)) - 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) - call prim2conpolytype(GRHydro_polytrope_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), scon(i,j,k,1), scon(i,j,k,2), scon(i,j,k,3), & - 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 if + 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_polytrope_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), scon(i,j,k,1), scon(i,j,k,2), scon(i,j,k,3), & + 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)) if (wk_atmosphere .eq. 0) then atmosphere_mask(i, j, k) = 0 SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits,\ diff --git a/src/Utils.F90 b/src/Utils.F90 index 440c945..0ef9af2 100644 --- a/src/Utils.F90 +++ b/src/Utils.F90 @@ -11,6 +11,7 @@ #include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" +#include "GRHydro_Macros.h" /*@@ @routine GRHydro_RefinementLevel @@ -42,6 +43,7 @@ end subroutine GRHydro_RefinementLevel @author @desc Calculates the determinant of the spatial metric. + PLEASE USE THE MACRO SPATIAL_DETERMINANT from now on!!! @enddesc @calls @calledby @@ -59,7 +61,7 @@ subroutine SpatialDeterminant(gxx,gxy,gxz,gyy,gyz,gzz,det) det = -gxz**2*gyy + 2*gxy*gxz*gyz - gxx*gyz**2 - gxy**2*gzz + gxx*gyy*gzz - + !!$ Why is this weird order necessary? Search me. It just seemed !!$ to make a really odd NaN go away. @@ -142,30 +144,14 @@ subroutine SetMetricTemps(CCTK_ARGUMENTS) do j = 1, ny do i = 1, nx - if (conformal_state .eq. 0) then - 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),& - GRHydro_det(i,j,k)) - call UpperMetric(& - GRHydro_uxx(i,j,k),GRHydro_uxy(i,j,k),GRHydro_uxz(i,j,k),& - GRHydro_uyy(i,j,k),GRHydro_uyz(i,j,k),GRHydro_uzz(i,j,k),& - GRHydro_det(i,j,k),& - gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),& - gyy(i,j,k),gyz(i,j,k),gzz(i,j,k)) - 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),& - GRHydro_det(i,j,k)) - call UpperMetric(& - GRHydro_uxx(i,j,k),GRHydro_uxy(i,j,k),GRHydro_uxz(i,j,k),& - GRHydro_uyy(i,j,k),GRHydro_uyz(i,j,k),GRHydro_uzz(i,j,k),& - GRHydro_det(i,j,k),& - 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)) - end if + GRHydro_det(i,j,k) = 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 UpperMetric(& + GRHydro_uxx(i,j,k),GRHydro_uxy(i,j,k),GRHydro_uxz(i,j,k),& + GRHydro_uyy(i,j,k),GRHydro_uyz(i,j,k),GRHydro_uzz(i,j,k),& + GRHydro_det(i,j,k),& + gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),& + gyy(i,j,k),gyz(i,j,k),gzz(i,j,k)) end do end do |