diff options
Diffstat (limited to 'src/GRHydro_Con2Prim.F90')
| -rw-r--r-- | src/GRHydro_Con2Prim.F90 | 125 |
1 files changed, 29 insertions, 96 deletions
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), & |