aboutsummaryrefslogtreecommitdiff
path: root/src/GRHydro_Prim2Con.F90
diff options
context:
space:
mode:
Diffstat (limited to 'src/GRHydro_Prim2Con.F90')
-rw-r--r--src/GRHydro_Prim2Con.F90131
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 + &