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