From f6be10b987d8db12434c2607dcd6d434a277157b Mon Sep 17 00:00:00 2001 From: pollney Date: Thu, 26 Jun 2003 11:37:17 +0000 Subject: Update of macros for 4th order differencing. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinBase/ADMMacros/trunk@53 b1d164ef-f17a-46e7-89d4-021c7118ef4e --- src/macro/ADM_Derivative.h | 93 +++++++++++++ src/macro/ADM_Spacing.h | 34 +++++ src/macro/ADM_Spacing_declare.h | 13 ++ src/macro/DA_guts.h | 18 ++- src/macro/DDA_guts.h | 31 +++-- src/macro/DXDB_guts.h | 16 ++- src/macro/DXDCG_guts.h | 26 ++-- src/macro/DXDK_guts.h | 25 ++-- src/macro/DXXDG_guts.h | 50 +++++-- src/macro/DXYDG_guts.h | 55 +++++--- src/macro/DXZDG_guts.h | 56 +++++--- src/macro/DYDB_guts.h | 16 ++- src/macro/DYDCG_guts.h | 26 ++-- src/macro/DYDK_guts.h | 26 ++-- src/macro/DYYDG_guts.h | 49 +++++-- src/macro/DYZDG_guts.h | 54 +++++--- src/macro/DZDB_guts.h | 16 ++- src/macro/DZDCG_guts.h | 25 ++-- src/macro/DZDK_guts.h | 26 ++-- src/macro/DZZDG_guts.h | 48 +++++-- src/macro/HAMADM_guts.h | 2 - src/macro/LIEG_guts.h | 228 ++++++++++--------------------- src/macro/LIEK_guts.h | 290 ++++++++++------------------------------ 23 files changed, 659 insertions(+), 564 deletions(-) create mode 100644 src/macro/ADM_Derivative.h create mode 100644 src/macro/ADM_Spacing.h create mode 100644 src/macro/ADM_Spacing_declare.h diff --git a/src/macro/ADM_Derivative.h b/src/macro/ADM_Derivative.h new file mode 100644 index 0000000..af8cc91 --- /dev/null +++ b/src/macro/ADM_Derivative.h @@ -0,0 +1,93 @@ +/*@@ + @header ADM_Derivative.h + @date June 2002 + @author Denis Pollney + @desc + Derivative operators. + @enddesc +@@*/ + +#ifndef ADM_DERIVATIVE_H +#define ADM_DERIVATIVE_H + +#include "ADM_Spacing.h" + +/* + * 2nd order operators. + */ + +#define ADM_DX_2(var,i,j,k) I2DX*(var(i+1,j,k) - var(i-1,j,k)) +#define ADM_DY_2(var,i,j,k) I2DY*(var(i,j+1,k) - var(i,j-1,k)) +#define ADM_DZ_2(var,i,j,k) I2DZ*(var(i,j,k+1) - var(i,j,k-1)) + +#define ADM_DXX_2(var,i,j,k) IDXX*(var(i+1,j,k) - 2.d0*var(i,j,k) \ + + var(i-1,j,k)) +#define ADM_DYY_2(var,i,j,k) IDYY*(var(i,j+1,k) - 2.d0*var(i,j,k) \ + + var(i,j-1,k)) +#define ADM_DZZ_2(var,i,j,k) IDZZ*(var(i,j,k+1) - 2.d0*var(i,j,k) \ + + var(i,j,k-1)) + +#define ADM_DXY_2(var,i,j,k) IDXY*(var(i+1,j+1,k) - var(i+1,j-1,k) \ + - var(i-1,j+1,k) + var(i-1,j-1,k)) +#define ADM_DXZ_2(var,i,j,k) IDXZ*(var(i+1,j,k+1) - var(i+1,j,k-1) \ + - var(i-1,j,k+1) + var(i-1,j,k-1)) +#define ADM_DYZ_2(var,i,j,k) IDYZ*(var(i,j+1,k+1) - var(i,j+1,k-1) \ + - var(i,j-1,k+1) + var(i,j-1,k-1)) +/* + * 4th order operators. + */ + +#define ADM_DX_4(var,i,j,k) (I2DX/6.d0)*(-var(i+2,j,k) + 8.d0*var(i+1,j,k) \ + - 8.d0*var(i-1,j,k) + var(i-2,j,k)) +#define ADM_DY_4(var,i,j,k) (I2DY/6.d0)*(-var(i,j+2,k) + 8.d0*var(i,j+1,k) \ + - 8.d0*var(i,j-1,k) + var(i,j-2,k)) +#define ADM_DZ_4(var,i,j,k) (I2DZ/6.d0)*(-var(i,j,k+2) + 8.d0*var(i,j,k+1) \ + - 8.d0*var(i,j,k-1) + var(i,j,k-2)) + +#define ADM_DXX_4(var,i,j,k) (IDXX/12.d0)*(-var(i+2,j,k) + 16.d0*var(i+1,j,k)\ + - 30.d0*var(i,j,k) + 16.d0*var(i-1,j,k) - var(i-2,j,k)) +#define ADM_DYY_4(var,i,j,k) (IDYY/12.d0)*(-var(i,j+2,k) + 16.d0*var(i,j+1,k)\ + - 30.d0*var(i,j,k) + 16.d0*var(i,j-1,k) - var(i,j-2,k)) +#define ADM_DZZ_4(var,i,j,k) (IDZZ/12.d0)*(-var(i,j,k+2) + 16.d0*var(i,j,k+1)\ + - 30.d0*var(i,j,k) + 16.d0*var(i,j,k-1) - var(i,j,k-2)) + +#define ADM_DXY_4(var,i,j,k) (IDXY/36.d0)* \ + (var(i+2,j+2,k) - var(i+2,j-2,k) - var(i-2,j+2,k) + var(i-2,j-2,k) \ + + 8.d0*(-var(i+2,j+1,k) + var(i+2,j-1,k) - var(i+1,j+2,k) + var(i+1,j-2,k) \ + + var(i-2,j+1,k) - var(i-2,j-1,k) - var(i-1,j-2,k) + var(i-1,j+2,k)) \ + + 64.d0*(var(i+1,j+1,k) - var(i+1,j-1,k) - var(i-1,j+1,k) + var(i-1,j-1,k))) +#define ADM_DXZ_4(var,i,j,k) (IDXZ/36.d0)* \ + (var(i+2,j,k+2) - var(i+2,j,k-2) - var(i-2,j,k+2) + var(i-2,j,k-2) \ + + 8.d0*(-var(i+2,j,k+1) + var(i+2,j,k-1) - var(i+1,j,k+2) + var(i+1,j,k-2) \ + + var(i-2,j,k+1) - var(i-2,j,k-1) - var(i-1,j,k-2) + var(i-1,j,k+2)) \ + + 64.d0*(var(i+1,j,k+1) - var(i+1,j,k-1) - var(i-1,j,k+1) + var(i-1,j,k-1))) +#define ADM_DYZ_4(var,i,j,k) (IDYZ/36.d0)* \ + (var(i,j+2,k+2) - var(i,j+2,k-2) - var(i,j-2,k+2) + var(i,j-2,k-2) \ + + 8.d0*(-var(i,j+2,k+1) + var(i,j+2,k-1) - var(i,j+1,k+2) + var(i,j+1,k-2) \ + + var(i,j-2,k+1) - var(i,j-2,k-1) - var(i,j-1,k-2) + var(i,j-1,k+2)) \ + + 64.d0*(var(i,j+1,k+1) - var(i,j+1,k-1) - var(i,j-1,k+1) + var(i,j-1,k-1))) + + +#define ADM_ADV_DX_P1(var,i,j,k) betax(i,j,k)*IDX*(var(i+i,j,k) - var(i,j,k)) +#define ADM_ADV_DY_P1(var,i,j,k) betay(i,j,k)*IDY*(var(i,j+1,k) - var(i,j,k)) +#define ADM_ADV_DZ_P1(var,i,j,k) betaz(i,j,k)*IDZ*(var(i,j,k+1) - var(i,j,k)) + +#define ADM_ADV_DX_M1(var,i,j,k) betax(i,j,k)*IDX*(var(i,j,k) - var(i-1,j,k)) +#define ADM_ADV_DY_M1(var,i,j,k) betay(i,j,k)*IDY*(var(i,j,k) - var(i,j-1,k)) +#define ADM_ADV_DZ_M1(var,i,j,k) betaz(i,j,k)*IDZ*(var(i,j,k) - var(i,j,k-1)) + +#define ADM_ADV_DX_P2(var,i,j,k) betax(i,j,k)*I2DX*(-3.d0*var(i,j,k) \ + + 4.d0*var(i+1,j,k) - var(i+2,j,k)) +#define ADM_ADV_DY_P2(var,i,j,k) betay(i,j,k)*I2DY*(-3.d0*var(i,j,k) \ + + 4.d0*var(i,j+1,k) - var(i,j+2,k)) +#define ADM_ADV_DZ_P2(var,i,j,k) betaz(i,j,k)*I2DZ*(-3.d0*var(i,j,k) \ + + 4.d0*var(i,j,k+1) - var(i,j,k+2)) + +#define ADM_ADV_DX_M2(var,i,j,k) betax(i,j,k)*I2DX*(3.d0*var(i,j,k) \ + - 4.d0*var(i-1,j,k) + var(i-2,j,k)) +#define ADM_ADV_DY_M2(var,i,j,k) betay(i,j,k)*I2DY*(3.d0*var(i,j,k) \ + - 4.d0*var(i,j-1,k) + var(i,j-2,k)) +#define ADM_ADV_DZ_M2(var,i,j,k) betaz(i,j,k)*I2DZ*(3.d0*var(i,j,k) \ + - 4.d0*var(i,j,k-1) + var(i,j,k-2)) + +#endif diff --git a/src/macro/ADM_Spacing.h b/src/macro/ADM_Spacing.h new file mode 100644 index 0000000..234d554 --- /dev/null +++ b/src/macro/ADM_Spacing.h @@ -0,0 +1,34 @@ +/*@@ + @header ADM_Spacing.h + @date June 2002 + @author Denis Pollney + @desc + Calculate various spacing dependent scalars. + @enddesc +@@*/ + +#ifndef ADM_SPACING_H +#define ADM_SPACING_H + + dt=cctk_delta_time + dx=cctk_delta_space(1) + dy=cctk_delta_space(2) + dz=cctk_delta_space(3) + + idx=1.d0/dx + idy=1.d0/dy + idz=1.d0/dz + + i2dx=0.5d0*idx + i2dy=0.5d0*idy + i2dz=0.5d0*idz + + idxx=idx**2 + idyy=idy**2 + idzz=idz**2 + + idxy=i2dx*i2dy + idxz=i2dx*i2dz + idyz=i2dy*i2dz + +#endif diff --git a/src/macro/ADM_Spacing_declare.h b/src/macro/ADM_Spacing_declare.h new file mode 100644 index 0000000..21a1ee6 --- /dev/null +++ b/src/macro/ADM_Spacing_declare.h @@ -0,0 +1,13 @@ +/*@@ + @header ADM_Spacing.h + @date June 2002 + @author Denis Pollney + @desc + Calculate various spacing dependent scalars. + @enddesc +@@*/ + + CCTK_REAL :: dt, dx, dy, dz + CCTK_REAL :: idx, idy, idz, i2dx, i2dy, i2dz + CCTK_REAL :: idxx, idxy, idxz, idyy, idyz, idzz + diff --git a/src/macro/DA_guts.h b/src/macro/DA_guts.h index 2b3ea63..cceb5c8 100644 --- a/src/macro/DA_guts.h +++ b/src/macro/DA_guts.h @@ -12,13 +12,17 @@ #ifdef FCODE - DA_OO2DX = 1D0/(2D0*DA_DX) - DA_OO2DY = 1D0/(2D0*DA_DY) - DA_OO2DZ = 1D0/(2D0*DA_DZ) - - DA_DXDA = DA_OO2DX*(DA_A_IP - DA_A_IM) - DA_DYDA = DA_OO2DY*(DA_A_JP - DA_A_JM) - DA_DZDA = DA_OO2DZ*(DA_A_KP - DA_A_KM) +#include "ADM_Derivative.h" + + if (spatial_order.eq.2) then + DA_DXDA = ADM_DX_2(alp,i,j,k) + DA_DYDA = ADM_DY_2(alp,i,j,k) + DA_DZDA = ADM_DZ_2(alp,i,j,k) + else + DA_DXDA = ADM_DX_4(alp,i,j,k) + DA_DYDA = ADM_DY_4(alp,i,j,k) + DA_DZDA = ADM_DZ_4(alp,i,j,k) + end if #endif diff --git a/src/macro/DDA_guts.h b/src/macro/DDA_guts.h index de4c849..7f115f0 100644 --- a/src/macro/DDA_guts.h +++ b/src/macro/DDA_guts.h @@ -12,20 +12,23 @@ #ifdef FCODE - DDA_OODX2 = 1D0/(DDA_DX*DDA_DX) - DDA_OODY2 = 1D0/(DDA_DY*DDA_DY) - DDA_OODZ2 = 1D0/(DDA_DZ*DDA_DZ) - DDA_OO4DXDY = 1D0/(4D0*DDA_DX*DDA_DY) - DDA_OO4DXDZ = 1D0/(4D0*DDA_DX*DDA_DZ) - DDA_OO4DYDZ = 1D0/(4D0*DDA_DY*DDA_DZ) - - DDA_DXXDA = DDA_OODX2*(DDA_A_IP - 2D0*DDA_A + DDA_A_IM) - DDA_DYYDA = DDA_OODY2*(DDA_A_JP - 2D0*DDA_A + DDA_A_JM) - DDA_DZZDA = DDA_OODZ2*(DDA_A_KP - 2D0*DDA_A + DDA_A_KM) - - DDA_DXYDA = DDA_OO4DXDY*(DDA_A_IPJP-DDA_A_IPJM-DDA_A_IMJP+DDA_A_IMJM) - DDA_DXZDA = DDA_OO4DXDZ*(DDA_A_IPKP-DDA_A_IPKM-DDA_A_IMKP+DDA_A_IMKM) - DDA_DYZDA = DDA_OO4DYDZ*(DDA_A_JPKP-DDA_A_JPKM-DDA_A_JMKP+DDA_A_JMKM) +#include "ADM_Derivative.h" + + if (spatial_order.eq.2) then + DDA_DXXDA = ADM_DXX_2(alp,i,j,k) + DDA_DXYDA = ADM_DXY_2(alp,i,j,k) + DDA_DXZDA = ADM_DXZ_2(alp,i,j,k) + DDA_DYYDA = ADM_DYY_2(alp,i,j,k) + DDA_DYZDA = ADM_DYZ_2(alp,i,j,k) + DDA_DZZDA = ADM_DZZ_2(alp,i,j,k) + else + DDA_DXXDA = ADM_DXX_4(alp,i,j,k) + DDA_DXYDA = ADM_DXY_4(alp,i,j,k) + DDA_DXZDA = ADM_DXZ_4(alp,i,j,k) + DDA_DYYDA = ADM_DYY_4(alp,i,j,k) + DDA_DYZDA = ADM_DYZ_4(alp,i,j,k) + DDA_DZZDA = ADM_DZZ_4(alp,i,j,k) + end if #endif diff --git a/src/macro/DXDB_guts.h b/src/macro/DXDB_guts.h index c47e81d..45a773c 100644 --- a/src/macro/DXDB_guts.h +++ b/src/macro/DXDB_guts.h @@ -17,11 +17,17 @@ #ifdef FCODE - DXDB_OO2DX = 1.0D0/(2.0D0*DXDB_DX) - - DXDB_DXDBX = DXDB_OO2DX*(DXDB_BX_IP - DXDB_BX_IM) - DXDB_DXDBY = DXDB_OO2DX*(DXDB_BY_IP - DXDB_BY_IM) - DXDB_DXDBZ = DXDB_OO2DX*(DXDB_BZ_IP - DXDB_BZ_IM) +#include "ADM_Derivative.h" + + if (spatial_order.eq.2) then + DXDB_DXDBX = ADM_DX_2(betax,i,j,k) + DXDB_DXDBY = ADM_DX_2(betay,i,j,k) + DXDB_DXDBZ = ADM_DX_2(betaz,i,j,k) + else + DXDB_DXDBX = ADM_DX_4(betax,i,j,k) + DXDB_DXDBY = ADM_DX_4(betay,i,j,k) + DXDB_DXDBZ = ADM_DX_4(betaz,i,j,k) + end if #endif diff --git a/src/macro/DXDCG_guts.h b/src/macro/DXDCG_guts.h index 464fdb7..b922068 100644 --- a/src/macro/DXDCG_guts.h +++ b/src/macro/DXDCG_guts.h @@ -21,15 +21,23 @@ #ifdef FCODE - DXDCG_OO2DX = 1D0/(2D0*DXDCG_DX) - - DXDCG_DXDCGXX = DXDCG_OO2DX*(DXDCG_GXX_IP - DXDCG_GXX_IM) - DXDCG_DXDCGXY = DXDCG_OO2DX*(DXDCG_GXY_IP - DXDCG_GXY_IM) - DXDCG_DXDCGXZ = DXDCG_OO2DX*(DXDCG_GXZ_IP - DXDCG_GXZ_IM) - DXDCG_DXDCGYY = DXDCG_OO2DX*(DXDCG_GYY_IP - DXDCG_GYY_IM) - DXDCG_DXDCGYZ = DXDCG_OO2DX*(DXDCG_GYZ_IP - DXDCG_GYZ_IM) - DXDCG_DXDCGZZ = DXDCG_OO2DX*(DXDCG_GZZ_IP - DXDCG_GZZ_IM) - +#include "ADM_Derivative.h" + + if (spatial_order.eq.2) then + DXDCG_DXDCGXX = ADM_DX_2(gxx,i,j,k) + DXDCG_DXDCGXY = ADM_DX_2(gxy,i,j,k) + DXDCG_DXDCGXZ = ADM_DX_2(gxz,i,j,k) + DXDCG_DXDCGYY = ADM_DX_2(gyy,i,j,k) + DXDCG_DXDCGYZ = ADM_DX_2(gyz,i,j,k) + DXDCG_DXDCGZZ = ADM_DX_2(gzz,i,j,k) + else + DXDCG_DXDCGXX = ADM_DX_4(gxx,i,j,k) + DXDCG_DXDCGXY = ADM_DX_4(gxy,i,j,k) + DXDCG_DXDCGXZ = ADM_DX_4(gxz,i,j,k) + DXDCG_DXDCGYY = ADM_DX_4(gyy,i,j,k) + DXDCG_DXDCGYZ = ADM_DX_4(gyz,i,j,k) + DXDCG_DXDCGZZ = ADM_DX_4(gzz,i,j,k) + end if #endif #ifdef CCODE diff --git a/src/macro/DXDK_guts.h b/src/macro/DXDK_guts.h index 6a6df3f..fd8af92 100644 --- a/src/macro/DXDK_guts.h +++ b/src/macro/DXDK_guts.h @@ -13,14 +13,23 @@ #ifdef FCODE - DXDK_OO2DX = 1D0/(2D0*DXDK_DX) - - DXDK_DXDKXX = DXDK_OO2DX*(DXDK_KXX_IP - DXDK_KXX_IM) - DXDK_DXDKXY = DXDK_OO2DX*(DXDK_KXY_IP - DXDK_KXY_IM) - DXDK_DXDKXZ = DXDK_OO2DX*(DXDK_KXZ_IP - DXDK_KXZ_IM) - DXDK_DXDKYY = DXDK_OO2DX*(DXDK_KYY_IP - DXDK_KYY_IM) - DXDK_DXDKYZ = DXDK_OO2DX*(DXDK_KYZ_IP - DXDK_KYZ_IM) - DXDK_DXDKZZ = DXDK_OO2DX*(DXDK_KZZ_IP - DXDK_KZZ_IM) +#include "ADM_Derivative.h" + + if (spatial_order.eq.2) then + DXDK_DXDKXX = ADM_DX_2(kxx,i,j,k) + DXDK_DXDKXY = ADM_DX_2(kxy,i,j,k) + DXDK_DXDKXZ = ADM_DX_2(kxz,i,j,k) + DXDK_DXDKYY = ADM_DX_2(kyy,i,j,k) + DXDK_DXDKYZ = ADM_DX_2(kyz,i,j,k) + DXDK_DXDKZZ = ADM_DX_2(kzz,i,j,k) + else + DXDK_DXDKXX = ADM_DX_4(kxx,i,j,k) + DXDK_DXDKXY = ADM_DX_4(kxy,i,j,k) + DXDK_DXDKXZ = ADM_DX_4(kxz,i,j,k) + DXDK_DXDKYY = ADM_DX_4(kyy,i,j,k) + DXDK_DXDKYZ = ADM_DX_4(kyz,i,j,k) + DXDK_DXDKZZ = ADM_DX_4(kzz,i,j,k) + end if #endif diff --git a/src/macro/DXXDG_guts.h b/src/macro/DXXDG_guts.h index 997a62b..89dcedf 100644 --- a/src/macro/DXXDG_guts.h +++ b/src/macro/DXXDG_guts.h @@ -20,33 +20,55 @@ #ifdef FCODE - DXXDG_OODX2 = 1D0/(dx*dx) +#include "ADM_Derivative.h" /* Factor involving 2nd derivative of conformal factor */ IF (conformal_state .eq. 0) THEN DXXDG_FAC = 0 ELSE - DXXDG_FAC = DXDG_PSI4*(4*DXXDG_DXXDPSI_O_PSI + 12*DXDG_DXDPSI_O_PSI*DXDG_DXDPSI_O_PSI) + DXXDG_FAC = DXDG_PSI4*(4*DXXDG_DXXDPSI_O_PSI \ + + 12*DXDG_DXDPSI_O_PSI*DXDG_DXDPSI_O_PSI) ENDIF /* Now calculate the second derivatives */ - DXXDG_DXXDGXX = 2*DXDCG_DXDCGXX*DXDG_FAC+DXXDG_FAC*DXDG_GXX+DXDG_PSI4 \ - *DXXDG_OODX2*(DXDCG_GXX_IP-2*DXDG_GXX+DXDCG_GXX_IM) - DXXDG_DXXDGXY = 2*DXDCG_DXDCGXY*DXDG_FAC+DXXDG_FAC*DXDG_GXY+DXDG_PSI4 \ - *DXXDG_OODX2*(DXDCG_GXY_IP-2*DXDG_GXY+DXDCG_GXY_IM) + if (spatial_order.eq.2) then + DXXDG_DXXDGXX = 2*DXDCG_DXDCGXX*DXDG_FAC + DXXDG_FAC*DXDG_GXX \ + + DXDG_PSI4*ADM_DXX_2(gxx,i,j,k) - DXXDG_DXXDGXZ = 2*DXDCG_DXDCGXZ*DXDG_FAC+DXXDG_FAC*DXDG_GXZ+DXDG_PSI4 \ - *DXXDG_OODX2*(DXDCG_GXZ_IP-2*DXDG_GXZ+DXDCG_GXZ_IM) + DXXDG_DXXDGXY = 2*DXDCG_DXDCGXY*DXDG_FAC + DXXDG_FAC*DXDG_GXY \ + + DXDG_PSI4*ADM_DXX_2(gxy,i,j,k) - DXXDG_DXXDGYY = 2*DXDCG_DXDCGYY*DXDG_FAC+DXXDG_FAC*DXDG_GYY+DXDG_PSI4 \ - *DXXDG_OODX2*(DXDCG_GYY_IP-2*DXDG_GYY+DXDCG_GYY_IM) + DXXDG_DXXDGXZ = 2*DXDCG_DXDCGXZ*DXDG_FAC + DXXDG_FAC*DXDG_GXZ \ + + DXDG_PSI4*ADM_DXX_2(gxz,i,j,k) - DXXDG_DXXDGYZ = 2*DXDCG_DXDCGYZ*DXDG_FAC+DXXDG_FAC*DXDG_GYZ+DXDG_PSI4 \ - *DXXDG_OODX2*(DXDCG_GYZ_IP-2*DXDG_GYZ+DXDCG_GYZ_IM) + DXXDG_DXXDGYY = 2*DXDCG_DXDCGYY*DXDG_FAC + DXXDG_FAC*DXDG_GYY \ + + DXDG_PSI4*ADM_DXX_2(gyy,i,j,k) - DXXDG_DXXDGZZ = 2*DXDCG_DXDCGZZ*DXDG_FAC+DXXDG_FAC*DXDG_GZZ+DXDG_PSI4 \ - *DXXDG_OODX2*(DXDCG_GZZ_IP-2*DXDG_GZZ+DXDCG_GZZ_IM) + DXXDG_DXXDGYZ = 2*DXDCG_DXDCGYZ*DXDG_FAC + DXXDG_FAC*DXDG_GYZ \ + + DXDG_PSI4*ADM_DXX_2(gyz,i,j,k) + + DXXDG_DXXDGZZ = 2*DXDCG_DXDCGZZ*DXDG_FAC + DXXDG_FAC*DXDG_GZZ \ + + DXDG_PSI4*ADM_DXX_2(gzz,i,j,k) + else + DXXDG_DXXDGXX = 2*DXDCG_DXDCGXX*DXDG_FAC + DXXDG_FAC*DXDG_GXX \ + + DXDG_PSI4*ADM_DXX_4(gxx,i,j,k) + + DXXDG_DXXDGXY = 2*DXDCG_DXDCGXY*DXDG_FAC + DXXDG_FAC*DXDG_GXY \ + + DXDG_PSI4*ADM_DXX_4(gxy,i,j,k) + + DXXDG_DXXDGXZ = 2*DXDCG_DXDCGXZ*DXDG_FAC + DXXDG_FAC*DXDG_GXZ \ + + DXDG_PSI4*ADM_DXX_4(gxz,i,j,k) + + DXXDG_DXXDGYY = 2*DXDCG_DXDCGYY*DXDG_FAC + DXXDG_FAC*DXDG_GYY \ + + DXDG_PSI4*ADM_DXX_4(gyy,i,j,k) + + DXXDG_DXXDGYZ = 2*DXDCG_DXDCGYZ*DXDG_FAC + DXXDG_FAC*DXDG_GYZ \ + + DXDG_PSI4*ADM_DXX_4(gyz,i,j,k) + + DXXDG_DXXDGZZ = 2*DXDCG_DXDCGZZ*DXDG_FAC + DXXDG_FAC*DXDG_GZZ \ + + DXDG_PSI4*ADM_DXX_4(gzz,i,j,k) + end if #endif diff --git a/src/macro/DXYDG_guts.h b/src/macro/DXYDG_guts.h index de1744e..f43de6c 100644 --- a/src/macro/DXYDG_guts.h +++ b/src/macro/DXYDG_guts.h @@ -22,39 +22,54 @@ #ifdef FCODE - DXYDG_OO4DXDY = 1D0/(4D0*dx*dy) +#include "ADM_Derivative.h" /* Factor involving 2nd derivative of conformal factor */ IF (conformal_state .eq. 0) THEN DXYDG_FAC = 0 ELSE - DXYDG_FAC = DXDG_PSI4*(4*DXYDG_DXYDPSI_O_PSI + 12*DXDG_DXDPSI_O_PSI*DYDG_DYDPSI_O_PSI) + DXYDG_FAC = DXDG_PSI4*(4*DXYDG_DXYDPSI_O_PSI + \ + 12*DXDG_DXDPSI_O_PSI*DYDG_DYDPSI_O_PSI) ENDIF /* Now calculate the second deriatives */ - DXYDG_DXYDGXX = DYDCG_DYDCGXX*DXDG_FAC+DXDCG_DXDCGXX*DYDG_FAC+DXYDG_FAC*DXDG_GXX \ - +DXDG_PSI4*DXYDG_OO4DXDY* \ - (DXYDG_GXX_IPJP-DXYDG_GXX_IPJM-DXYDG_GXX_IMJP+DXYDG_GXX_IMJM); + if (spatial_order.eq.2) then + DXYDG_DXYDGXX = DYDCG_DYDCGXX*DXDG_FAC + DXDCG_DXDCGXX*DYDG_FAC \ + + DXYDG_FAC*DXDG_GXX + DXDG_PSI4*ADM_DXY_2(gxx,i,j,k) - DXYDG_DXYDGXY = DYDCG_DYDCGXY*DXDG_FAC+DXDCG_DXDCGXY*DYDG_FAC+DXYDG_FAC*DXDG_GXY \ - +DXDG_PSI4*DXYDG_OO4DXDY* \ - (DXYDG_GXY_IPJP-DXYDG_GXY_IPJM-DXYDG_GXY_IMJP+DXYDG_GXY_IMJM); + DXYDG_DXYDGXY = DYDCG_DYDCGXY*DXDG_FAC + DXDCG_DXDCGXY*DYDG_FAC \ + + DXYDG_FAC*DXDG_GXY + DXDG_PSI4*ADM_DXY_2(gxy,i,j,k) - DXYDG_DXYDGXZ = DYDCG_DYDCGXZ*DXDG_FAC+DXDCG_DXDCGXZ*DYDG_FAC+DXYDG_FAC*DXDG_GXZ \ - +DXDG_PSI4*DXYDG_OO4DXDY* \ - (DXYDG_GXZ_IPJP-DXYDG_GXZ_IPJM-DXYDG_GXZ_IMJP+DXYDG_GXZ_IMJM); + DXYDG_DXYDGXZ = DYDCG_DYDCGXZ*DXDG_FAC + DXDCG_DXDCGXZ*DYDG_FAC \ + + DXYDG_FAC*DXDG_GXZ + DXDG_PSI4*ADM_DXY_2(gxz,i,j,k) - DXYDG_DXYDGYY = DYDCG_DYDCGYY*DXDG_FAC+DXDCG_DXDCGYY*DYDG_FAC+DXYDG_FAC*DXDG_GYY \ - +DXDG_PSI4*DXYDG_OO4DXDY* \ - (DXYDG_GYY_IPJP-DXYDG_GYY_IPJM-DXYDG_GYY_IMJP+DXYDG_GYY_IMJM); + DXYDG_DXYDGYY = DYDCG_DYDCGYY*DXDG_FAC + DXDCG_DXDCGYY*DYDG_FAC \ + + DXYDG_FAC*DXDG_GYY + DXDG_PSI4*ADM_DXY_2(gyy,i,j,k) - DXYDG_DXYDGYZ = DYDCG_DYDCGYZ*DXDG_FAC+DXDCG_DXDCGYZ*DYDG_FAC+DXYDG_FAC*DXDG_GYZ \ - +DXDG_PSI4*DXYDG_OO4DXDY* \ - (DXYDG_GYZ_IPJP-DXYDG_GYZ_IPJM-DXYDG_GYZ_IMJP+DXYDG_GYZ_IMJM); + DXYDG_DXYDGYZ = DYDCG_DYDCGYZ*DXDG_FAC + DXDCG_DXDCGYZ*DYDG_FAC \ + + DXYDG_FAC*DXDG_GYZ + DXDG_PSI4*ADM_DXY_2(gyz,i,j,k) - DXYDG_DXYDGZZ = DYDCG_DYDCGZZ*DXDG_FAC+DXDCG_DXDCGZZ*DYDG_FAC+DXYDG_FAC*DXDG_GZZ \ - +DXDG_PSI4*DXYDG_OO4DXDY* \ - (DXYDG_GZZ_IPJP-DXYDG_GZZ_IPJM-DXYDG_GZZ_IMJP+DXYDG_GZZ_IMJM); + DXYDG_DXYDGZZ = DYDCG_DYDCGZZ*DXDG_FAC + DXDCG_DXDCGZZ*DYDG_FAC \ + + DXYDG_FAC*DXDG_GZZ + DXDG_PSI4*ADM_DXY_2(gzz,i,j,k) + else + DXYDG_DXYDGXX = DYDCG_DYDCGXX*DXDG_FAC + DXDCG_DXDCGXX*DYDG_FAC \ + + DXYDG_FAC*DXDG_GXX + DXDG_PSI4*ADM_DXY_4(gxx,i,j,k) + + DXYDG_DXYDGXY = DYDCG_DYDCGXY*DXDG_FAC + DXDCG_DXDCGXY*DYDG_FAC \ + + DXYDG_FAC*DXDG_GXY + DXDG_PSI4*ADM_DXY_4(gxy,i,j,k) + + DXYDG_DXYDGXZ = DYDCG_DYDCGXZ*DXDG_FAC + DXDCG_DXDCGXZ*DYDG_FAC \ + + DXYDG_FAC*DXDG_GXZ + DXDG_PSI4*ADM_DXY_4(gxz,i,j,k) + + DXYDG_DXYDGYY = DYDCG_DYDCGYY*DXDG_FAC + DXDCG_DXDCGYY*DYDG_FAC \ + + DXYDG_FAC*DXDG_GYY + DXDG_PSI4*ADM_DXY_4(gyy,i,j,k) + + DXYDG_DXYDGYZ = DYDCG_DYDCGYZ*DXDG_FAC + DXDCG_DXDCGYZ*DYDG_FAC \ + + DXYDG_FAC*DXDG_GYZ + DXDG_PSI4*ADM_DXY_4(gyz,i,j,k) + + DXYDG_DXYDGZZ = DYDCG_DYDCGZZ*DXDG_FAC + DXDCG_DXDCGZZ*DYDG_FAC \ + + DXYDG_FAC*DXDG_GZZ + DXDG_PSI4*ADM_DXY_4(gzz,i,j,k) + end if #endif diff --git a/src/macro/DXZDG_guts.h b/src/macro/DXZDG_guts.h index 4978507..cdee51c 100644 --- a/src/macro/DXZDG_guts.h +++ b/src/macro/DXZDG_guts.h @@ -22,40 +22,56 @@ #ifdef FCODE - DXZDG_OO4DXDZ = 1D0/(4D0*dx*dz) +#include "ADM_Derivative.h" /* Factor involving 2nd derivative of conformal factor */ IF (conformal_state .eq. 0) THEN DXZDG_FAC = 0 ELSE - DXZDG_FAC = DXDG_PSI4*(4*DXZDG_DXZDPSI_O_PSI + 12*DXDG_DXDPSI_O_PSI*DZDG_DZDPSI_O_PSI) + DXZDG_FAC = DXDG_PSI4*(4*DXZDG_DXZDPSI_O_PSI \ + + 12*DXDG_DXDPSI_O_PSI*DZDG_DZDPSI_O_PSI) ENDIF /* Now calculate the second deriatives */ - DXZDG_DXZDGXX = DZDCG_DZDCGXX*DXDG_FAC+DXDCG_DXDCGXX*DZDG_FAC+DXZDG_FAC*DXDG_GXX+ \ - DXDG_PSI4*DXZDG_OO4DXDZ*(DXZDG_GXX_IPKP-DXZDG_GXX_IPKM-DXZDG_GXX_IMKP+ \ - DXZDG_GXX_IMKM) + if (spatial_order.eq.2) then + DXZDG_DXZDGXX = DZDCG_DZDCGXX*DXDG_FAC + DXDCG_DXDCGXX*DZDG_FAC \ + + DXZDG_FAC*DXDG_GXX + DXDG_PSI4*ADM_DXZ_2(gxx,i,j,k) - DXZDG_DXZDGXY = DZDCG_DZDCGXY*DXDG_FAC+DXDCG_DXDCGXY*DZDG_FAC+DXZDG_FAC*DXDG_GXY+ \ - DXDG_PSI4*DXZDG_OO4DXDZ*(DXZDG_GXY_IPKP-DXZDG_GXY_IPKM-DXZDG_GXY_IMKP+ \ - DXZDG_GXY_IMKM) + DXZDG_DXZDGXY = DZDCG_DZDCGXY*DXDG_FAC + DXDCG_DXDCGXY*DZDG_FAC \ + + DXZDG_FAC*DXDG_GXY + DXDG_PSI4*ADM_DXZ_2(gxy,i,j,k) - DXZDG_DXZDGXZ = DZDCG_DZDCGXZ*DXDG_FAC+DXDCG_DXDCGXZ*DZDG_FAC+DXZDG_FAC*DXDG_GXZ+ \ - DXDG_PSI4*DXZDG_OO4DXDZ*(DXZDG_GXZ_IPKP-DXZDG_GXZ_IPKM-DXZDG_GXZ_IMKP+ \ - DXZDG_GXZ_IMKM) + DXZDG_DXZDGXZ = DZDCG_DZDCGXZ*DXDG_FAC + DXDCG_DXDCGXZ*DZDG_FAC \ + + DXZDG_FAC*DXDG_GXZ + DXDG_PSI4*ADM_DXZ_2(gxz,i,j,k) - DXZDG_DXZDGYY = DZDCG_DZDCGYY*DXDG_FAC+DXDCG_DXDCGYY*DZDG_FAC+DXZDG_FAC*DXDG_GYY+ \ - DXDG_PSI4*DXZDG_OO4DXDZ*(DXZDG_GYY_IPKP-DXZDG_GYY_IPKM-DXZDG_GYY_IMKP+ \ - DXZDG_GYY_IMKM) + DXZDG_DXZDGYY = DZDCG_DZDCGYY*DXDG_FAC + DXDCG_DXDCGYY*DZDG_FAC \ + + DXZDG_FAC*DXDG_GYY + DXDG_PSI4*ADM_DXZ_2(gyy,i,j,k) - DXZDG_DXZDGYZ = DZDCG_DZDCGYZ*DXDG_FAC+DXDCG_DXDCGYZ*DZDG_FAC+DXZDG_FAC*DXDG_GYZ+ \ - DXDG_PSI4*DXZDG_OO4DXDZ*(DXZDG_GYZ_IPKP-DXZDG_GYZ_IPKM-DXZDG_GYZ_IMKP+ \ - DXZDG_GYZ_IMKM) + DXZDG_DXZDGYZ = DZDCG_DZDCGYZ*DXDG_FAC + DXDCG_DXDCGYZ*DZDG_FAC \ + + DXZDG_FAC*DXDG_GYZ + DXDG_PSI4*ADM_DXZ_2(gyz,i,j,k) - DXZDG_DXZDGZZ = DZDCG_DZDCGZZ*DXDG_FAC+DXDCG_DXDCGZZ*DZDG_FAC+DXZDG_FAC*DXDG_GZZ+ \ - DXDG_PSI4*DXZDG_OO4DXDZ*(DXZDG_GZZ_IPKP-DXZDG_GZZ_IPKM-DXZDG_GZZ_IMKP+ \ - DXZDG_GZZ_IMKM) + DXZDG_DXZDGZZ = DZDCG_DZDCGZZ*DXDG_FAC + DXDCG_DXDCGZZ*DZDG_FAC \ + + DXZDG_FAC*DXDG_GZZ + DXDG_PSI4*ADM_DXZ_2(gzz,i,j,k) + else + + DXZDG_DXZDGXX = DZDCG_DZDCGXX*DXDG_FAC + DXDCG_DXDCGXX*DZDG_FAC \ + + DXZDG_FAC*DXDG_GXX + DXDG_PSI4*ADM_DXZ_4(gxx,i,j,k) + + DXZDG_DXZDGXY = DZDCG_DZDCGXY*DXDG_FAC + DXDCG_DXDCGXY*DZDG_FAC \ + + DXZDG_FAC*DXDG_GXY + DXDG_PSI4*ADM_DXZ_4(gxy,i,j,k) + + DXZDG_DXZDGXZ = DZDCG_DZDCGXZ*DXDG_FAC + DXDCG_DXDCGXZ*DZDG_FAC \ + + DXZDG_FAC*DXDG_GXZ + DXDG_PSI4*ADM_DXZ_4(gxz,i,j,k) + + DXZDG_DXZDGYY = DZDCG_DZDCGYY*DXDG_FAC + DXDCG_DXDCGYY*DZDG_FAC \ + + DXZDG_FAC*DXDG_GYY + DXDG_PSI4*ADM_DXZ_4(gyy,i,j,k) + + DXZDG_DXZDGYZ = DZDCG_DZDCGYZ*DXDG_FAC + DXDCG_DXDCGYZ*DZDG_FAC \ + + DXZDG_FAC*DXDG_GYZ + DXDG_PSI4*ADM_DXZ_4(gyz,i,j,k) + + DXZDG_DXZDGZZ = DZDCG_DZDCGZZ*DXDG_FAC + DXDCG_DXDCGZZ*DZDG_FAC \ + + DXZDG_FAC*DXDG_GZZ + DXDG_PSI4*ADM_DXZ_4(gzz,i,j,k) + end if #endif diff --git a/src/macro/DYDB_guts.h b/src/macro/DYDB_guts.h index 54cba6c..bd92f57 100644 --- a/src/macro/DYDB_guts.h +++ b/src/macro/DYDB_guts.h @@ -17,11 +17,17 @@ #ifdef FCODE - DYDB_OO2DY = 1.0D0/(2.0D0*DYDB_DY) - - DYDB_DYDBX = DYDB_OO2DY*(DYDB_BX_JP - DYDB_BX_JM) - DYDB_DYDBY = DYDB_OO2DY*(DYDB_BY_JP - DYDB_BY_JM) - DYDB_DYDBZ = DYDB_OO2DY*(DYDB_BZ_JP - DYDB_BZ_JM) +#include "ADM_Derivative.h" + + if (spatial_order.eq.2) then + DYDB_DYDBX = ADM_DY_2(betax,i,j,k) + DYDB_DYDBY = ADM_DY_2(betay,i,j,k) + DYDB_DYDBZ = ADM_DY_2(betaz,i,j,k) + else + DYDB_DYDBX = ADM_DY_4(betax,i,j,k) + DYDB_DYDBY = ADM_DY_4(betay,i,j,k) + DYDB_DYDBZ = ADM_DY_4(betaz,i,j,k) + end if #endif diff --git a/src/macro/DYDCG_guts.h b/src/macro/DYDCG_guts.h index 01bced2..3c328d7 100644 --- a/src/macro/DYDCG_guts.h +++ b/src/macro/DYDCG_guts.h @@ -21,15 +21,23 @@ #ifdef FCODE - DYDCG_OO2DY = 1D0/(2D0*DYDCG_DY) - - DYDCG_DYDCGXX = DYDCG_OO2DY*(DYDCG_GXX_JP - DYDCG_GXX_JM) - DYDCG_DYDCGXY = DYDCG_OO2DY*(DYDCG_GXY_JP - DYDCG_GXY_JM) - DYDCG_DYDCGXZ = DYDCG_OO2DY*(DYDCG_GXZ_JP - DYDCG_GXZ_JM) - DYDCG_DYDCGYY = DYDCG_OO2DY*(DYDCG_GYY_JP - DYDCG_GYY_JM) - DYDCG_DYDCGYZ = DYDCG_OO2DY*(DYDCG_GYZ_JP - DYDCG_GYZ_JM) - DYDCG_DYDCGZZ = DYDCG_OO2DY*(DYDCG_GZZ_JP - DYDCG_GZZ_JM) - +#include "ADM_Derivative.h" + + if (spatial_order.eq.2) then + DYDCG_DYDCGXX = ADM_DY_2(gxx,i,j,k) + DYDCG_DYDCGXY = ADM_DY_2(gxy,i,j,k) + DYDCG_DYDCGXZ = ADM_DY_2(gxz,i,j,k) + DYDCG_DYDCGYY = ADM_DY_2(gyy,i,j,k) + DYDCG_DYDCGYZ = ADM_DY_2(gyz,i,j,k) + DYDCG_DYDCGZZ = ADM_DY_2(gzz,i,j,k) + else + DYDCG_DYDCGXX = ADM_DY_4(gxx,i,j,k) + DYDCG_DYDCGXY = ADM_DY_4(gxy,i,j,k) + DYDCG_DYDCGXZ = ADM_DY_4(gxz,i,j,k) + DYDCG_DYDCGYY = ADM_DY_4(gyy,i,j,k) + DYDCG_DYDCGYZ = ADM_DY_4(gyz,i,j,k) + DYDCG_DYDCGZZ = ADM_DY_4(gzz,i,j,k) + end if #endif diff --git a/src/macro/DYDK_guts.h b/src/macro/DYDK_guts.h index 26364b3..7dbfc0c 100644 --- a/src/macro/DYDK_guts.h +++ b/src/macro/DYDK_guts.h @@ -13,15 +13,23 @@ #ifdef FCODE - DYDK_OO2DY = 1D0/(2D0*DYDK_DY) - - DYDK_DYDKXX = DYDK_OO2DY*(DYDK_KXX_JP - DYDK_KXX_JM) - DYDK_DYDKXY = DYDK_OO2DY*(DYDK_KXY_JP - DYDK_KXY_JM) - DYDK_DYDKXZ = DYDK_OO2DY*(DYDK_KXZ_JP - DYDK_KXZ_JM) - DYDK_DYDKYY = DYDK_OO2DY*(DYDK_KYY_JP - DYDK_KYY_JM) - DYDK_DYDKYZ = DYDK_OO2DY*(DYDK_KYZ_JP - DYDK_KYZ_JM) - DYDK_DYDKZZ = DYDK_OO2DY*(DYDK_KZZ_JP - DYDK_KZZ_JM) - +#include "ADM_Derivative.h" + + if (spatial_order.eq.2) then + DYDK_DYDKXX = ADM_DY_2(kxx,i,j,k) + DYDK_DYDKXY = ADM_DY_2(kxy,i,j,k) + DYDK_DYDKXZ = ADM_DY_2(kxz,i,j,k) + DYDK_DYDKYY = ADM_DY_2(kyy,i,j,k) + DYDK_DYDKYZ = ADM_DY_2(kyz,i,j,k) + DYDK_DYDKZZ = ADM_DY_2(kzz,i,j,k) + else + DYDK_DYDKXX = ADM_DY_4(kxx,i,j,k) + DYDK_DYDKXY = ADM_DY_4(kxy,i,j,k) + DYDK_DYDKXZ = ADM_DY_4(kxz,i,j,k) + DYDK_DYDKYY = ADM_DY_4(kyy,i,j,k) + DYDK_DYDKYZ = ADM_DY_4(kyz,i,j,k) + DYDK_DYDKZZ = ADM_DY_4(kzz,i,j,k) + end if #endif diff --git a/src/macro/DYYDG_guts.h b/src/macro/DYYDG_guts.h index 1896229..2dc495e 100644 --- a/src/macro/DYYDG_guts.h +++ b/src/macro/DYYDG_guts.h @@ -20,33 +20,54 @@ #ifdef FCODE - DYYDG_OODY2 = 1D0/(dy*dy) +#include "ADM_Derivative.h" /* Factor involving 2nd derivative of conformal factor */ IF (conformal_state .eq. 0) THEN DYYDG_FAC = 0 ELSE - DYYDG_FAC = DYDG_PSI4*(4*DYYDG_DYYDPSI_O_PSI + 12*DYDG_DYDPSI_O_PSI*DYDG_DYDPSI_O_PSI) + DYYDG_FAC = DYDG_PSI4*(4*DYYDG_DYYDPSI_O_PSI \ + + 12*DYDG_DYDPSI_O_PSI*DYDG_DYDPSI_O_PSI) ENDIF /* Now calculate the second deriatives */ - DYYDG_DYYDGXX = 2*DYDCG_DYDCGXX*DYDG_FAC+DYYDG_FAC*DYDG_GXX+DYDG_PSI4 \ - *DYYDG_OODY2*(DYDCG_GXX_JP-2*DYDG_GXX+DYDCG_GXX_JM) + if (spatial_order.eq.2) then + DYYDG_DYYDGXX = 2*DYDCG_DYDCGXX*DYDG_FAC + DYYDG_FAC*DYDG_GXX \ + + DYDG_PSI4*ADM_DYY_2(gxx,i,j,k) - DYYDG_DYYDGXY = 2*DYDCG_DYDCGXY*DYDG_FAC+DYYDG_FAC*DYDG_GXY+DYDG_PSI4 \ - *DYYDG_OODY2*(DYDCG_GXY_JP-2*DYDG_GXY+DYDCG_GXY_JM) + DYYDG_DYYDGXY = 2*DYDCG_DYDCGXY*DYDG_FAC + DYYDG_FAC*DYDG_GXY \ + + DYDG_PSI4*ADM_DYY_2(gxy,i,j,k) - DYYDG_DYYDGXZ = 2*DYDCG_DYDCGXZ*DYDG_FAC+DYYDG_FAC*DYDG_GXZ+DYDG_PSI4 \ - *DYYDG_OODY2*(DYDCG_GXZ_JP-2*DYDG_GXZ+DYDCG_GXZ_JM) + DYYDG_DYYDGXZ = 2*DYDCG_DYDCGXZ*DYDG_FAC + DYYDG_FAC*DYDG_GXZ \ + + DYDG_PSI4*ADM_DYY_2(gxz,i,j,k) - DYYDG_DYYDGYY = 2*DYDCG_DYDCGYY*DYDG_FAC+DYYDG_FAC*DYDG_GYY+DYDG_PSI4 \ - *DYYDG_OODY2*(DYDCG_GYY_JP-2*DYDG_GYY+DYDCG_GYY_JM) + DYYDG_DYYDGYY = 2*DYDCG_DYDCGYY*DYDG_FAC + DYYDG_FAC*DYDG_GYY \ + + DYDG_PSI4*ADM_DYY_2(gyy,i,j,k) - DYYDG_DYYDGYZ = 2*DYDCG_DYDCGYZ*DYDG_FAC+DYYDG_FAC*DYDG_GYZ+DYDG_PSI4 \ - *DYYDG_OODY2*(DYDCG_GYZ_JP-2*DYDG_GYZ+DYDCG_GYZ_JM) + DYYDG_DYYDGYZ = 2*DYDCG_DYDCGYZ*DYDG_FAC + DYYDG_FAC*DYDG_GYZ \ + + DYDG_PSI4*ADM_DYY_2(gyz,i,j,k) - DYYDG_DYYDGZZ = 2*DYDCG_DYDCGZZ*DYDG_FAC+DYYDG_FAC*DYDG_GZZ+DYDG_PSI4 \ - *DYYDG_OODY2*(DYDCG_GZZ_JP-2*DYDG_GZZ+DYDCG_GZZ_JM) + DYYDG_DYYDGZZ = 2*DYDCG_DYDCGZZ*DYDG_FAC + DYYDG_FAC*DYDG_GZZ \ + + DYDG_PSI4*ADM_DYY_2(gzz,i,j,k) + else + DYYDG_DYYDGXX = 2*DYDCG_DYDCGXX*DYDG_FAC + DYYDG_FAC*DYDG_GXX \ + + DYDG_PSI4*ADM_DYY_4(gxx,i,j,k) + + DYYDG_DYYDGXY = 2*DYDCG_DYDCGXY*DYDG_FAC + DYYDG_FAC*DYDG_GXY \ + + DYDG_PSI4*ADM_DYY_4(gxy,i,j,k) + + DYYDG_DYYDGXZ = 2*DYDCG_DYDCGXZ*DYDG_FAC + DYYDG_FAC*DYDG_GXZ \ + + DYDG_PSI4*ADM_DYY_4(gxz,i,j,k) + + DYYDG_DYYDGYY = 2*DYDCG_DYDCGYY*DYDG_FAC + DYYDG_FAC*DYDG_GYY \ + + DYDG_PSI4*ADM_DYY_4(gyy,i,j,k) + + DYYDG_DYYDGYZ = 2*DYDCG_DYDCGYZ*DYDG_FAC + DYYDG_FAC*DYDG_GYZ \ + + DYDG_PSI4*ADM_DYY_4(gyz,i,j,k) + + DYYDG_DYYDGZZ = 2*DYDCG_DYDCGZZ*DYDG_FAC + DYYDG_FAC*DYDG_GZZ \ + + DYDG_PSI4*ADM_DYY_4(gzz,i,j,k) + end if #endif diff --git a/src/macro/DYZDG_guts.h b/src/macro/DYZDG_guts.h index 80347ff..65a513f 100644 --- a/src/macro/DYZDG_guts.h +++ b/src/macro/DYZDG_guts.h @@ -22,40 +22,54 @@ #ifdef FCODE - DYZDG_OO4DYDZ = 1D0/(4D0*dy*dz) +#include "ADM_Derivative.h" /* Factor involving 2nd derivative of conformal factor */ IF (conformal_state .eq. 0) THEN DYZDG_FAC = 0 ELSE - DYZDG_FAC = DYDG_PSI4*(4*DYZDG_DYZDPSI_O_PSI + 12*DYDG_DYDPSI_O_PSI*DZDG_DZDPSI_O_PSI) + DYZDG_FAC = DYDG_PSI4*(4*DYZDG_DYZDPSI_O_PSI \ + + 12*DYDG_DYDPSI_O_PSI*DZDG_DZDPSI_O_PSI) ENDIF /* Now calculate the second deriatives */ - DYZDG_DYZDGXX = DZDCG_DZDCGXX*DYDG_FAC+DYDCG_DYDCGXX*DZDG_FAC+DYZDG_FAC*DYDG_GXX \ - +DYDG_PSI4*DYZDG_OO4DYDZ*(DYZDG_GXX_JPKP-DYZDG_GXX_JPKM-DYZDG_GXX_JMKP \ - +DYZDG_GXX_JMKM) + if (spatial_order.eq.2) then + DYZDG_DYZDGXX = DZDCG_DZDCGXX*DYDG_FAC + DYDCG_DYDCGXX*DZDG_FAC \ + + DYZDG_FAC*DYDG_GXX + DYDG_PSI4*ADM_DYZ_2(gxx,i,j,k) - DYZDG_DYZDGXY = DZDCG_DZDCGXY*DYDG_FAC+DYDCG_DYDCGXY*DZDG_FAC+DYZDG_FAC*DYDG_GXY \ - +DYDG_PSI4*DYZDG_OO4DYDZ*(DYZDG_GXY_JPKP-DYZDG_GXY_JPKM-DYZDG_GXY_JMKP \ - +DYZDG_GXY_JMKM) + DYZDG_DYZDGXY = DZDCG_DZDCGXY*DYDG_FAC + DYDCG_DYDCGXY*DZDG_FAC \ + + DYZDG_FAC*DYDG_GXY + DYDG_PSI4*ADM_DYZ_2(gxy,i,j,k) - DYZDG_DYZDGXZ = DZDCG_DZDCGXZ*DYDG_FAC+DYDCG_DYDCGXZ*DZDG_FAC+DYZDG_FAC*DYDG_GXZ \ - +DYDG_PSI4*DYZDG_OO4DYDZ*(DYZDG_GXZ_JPKP-DYZDG_GXZ_JPKM-DYZDG_GXZ_JMKP \ - +DYZDG_GXZ_JMKM) + DYZDG_DYZDGXZ = DZDCG_DZDCGXZ*DYDG_FAC + DYDCG_DYDCGXZ*DZDG_FAC \ + + DYZDG_FAC*DYDG_GXZ + DYDG_PSI4*ADM_DYZ_2(gxz,i,j,k) - DYZDG_DYZDGYY = DZDCG_DZDCGYY*DYDG_FAC+DYDCG_DYDCGYY*DZDG_FAC+DYZDG_FAC*DYDG_GYY \ - +DYDG_PSI4*DYZDG_OO4DYDZ*(DYZDG_GYY_JPKP-DYZDG_GYY_JPKM-DYZDG_GYY_JMKP \ - +DYZDG_GYY_JMKM) + DYZDG_DYZDGYY = DZDCG_DZDCGYY*DYDG_FAC + DYDCG_DYDCGYY*DZDG_FAC \ + + DYZDG_FAC*DYDG_GYY + DYDG_PSI4*ADM_DYZ_2(gyy,i,j,k) - DYZDG_DYZDGYZ = DZDCG_DZDCGYZ*DYDG_FAC+DYDCG_DYDCGYZ*DZDG_FAC+DYZDG_FAC*DYDG_GYZ \ - +DYDG_PSI4*DYZDG_OO4DYDZ*(DYZDG_GYZ_JPKP-DYZDG_GYZ_JPKM-DYZDG_GYZ_JMKP \ - +DYZDG_GYZ_JMKM) + DYZDG_DYZDGYZ = DZDCG_DZDCGYZ*DYDG_FAC + DYDCG_DYDCGYZ*DZDG_FAC \ + + DYZDG_FAC*DYDG_GYZ + DYDG_PSI4*ADM_DYZ_2(gyz,i,j,k) - DYZDG_DYZDGZZ = DZDCG_DZDCGZZ*DYDG_FAC+DYDCG_DYDCGZZ*DZDG_FAC+DYZDG_FAC*DYDG_GZZ \ - +DYDG_PSI4*DYZDG_OO4DYDZ*(DYZDG_GZZ_JPKP-DYZDG_GZZ_JPKM-DYZDG_GZZ_JMKP \ - +DYZDG_GZZ_JMKM) + DYZDG_DYZDGZZ = DZDCG_DZDCGZZ*DYDG_FAC + DYDCG_DYDCGZZ*DZDG_FAC \ + + DYZDG_FAC*DYDG_GZZ + DYDG_PSI4*ADM_DYZ_2(gzz,i,j,k) + else + DYZDG_DYZDGXX = DZDCG_DZDCGXX*DYDG_FAC + DYDCG_DYDCGXX*DZDG_FAC \ + + DYZDG_FAC*DYDG_GXX + DYDG_PSI4*ADM_DYZ_4(gxx,i,j,k) + DYZDG_DYZDGXY = DZDCG_DZDCGXY*DYDG_FAC + DYDCG_DYDCGXY*DZDG_FAC \ + + DYZDG_FAC*DYDG_GXY + DYDG_PSI4*ADM_DYZ_4(gxy,i,j,k) + + DYZDG_DYZDGXZ = DZDCG_DZDCGXZ*DYDG_FAC + DYDCG_DYDCGXZ*DZDG_FAC \ + + DYZDG_FAC*DYDG_GXZ + DYDG_PSI4*ADM_DYZ_4(gxz,i,j,k) + + DYZDG_DYZDGYY = DZDCG_DZDCGYY*DYDG_FAC + DYDCG_DYDCGYY*DZDG_FAC \ + + DYZDG_FAC*DYDG_GYY + DYDG_PSI4*ADM_DYZ_4(gyy,i,j,k) + + DYZDG_DYZDGYZ = DZDCG_DZDCGYZ*DYDG_FAC + DYDCG_DYDCGYZ*DZDG_FAC \ + + DYZDG_FAC*DYDG_GYZ + DYDG_PSI4*ADM_DYZ_4(gyz,i,j,k) + + DYZDG_DYZDGZZ = DZDCG_DZDCGZZ*DYDG_FAC + DYDCG_DYDCGZZ*DZDG_FAC \ + + DYZDG_FAC*DYDG_GZZ + DYDG_PSI4*ADM_DYZ_4(gzz,i,j,k) + end if #endif #ifdef CCODE diff --git a/src/macro/DZDB_guts.h b/src/macro/DZDB_guts.h index 5797db7..4b01986 100644 --- a/src/macro/DZDB_guts.h +++ b/src/macro/DZDB_guts.h @@ -17,12 +17,16 @@ #ifdef FCODE - DZDB_OO2DZ = 1.0D0/(2.0D0*DZDB_DZ) - - DZDB_DZDBX = DZDB_OO2DZ*(DZDB_BX_KP - DZDB_BX_KM) - DZDB_DZDBY = DZDB_OO2DZ*(DZDB_BY_KP - DZDB_BY_KM) - DZDB_DZDBZ = DZDB_OO2DZ*(DZDB_BZ_KP - DZDB_BZ_KM) - +#include "ADM_Derivative.h" + if (spatial_order.eq.2) then + DZDB_DZDBX = ADM_DZ_2(betax,i,j,k) + DZDB_DZDBY = ADM_DZ_2(betay,i,j,k) + DZDB_DZDBZ = ADM_DZ_2(betaz,i,j,k) + else + DZDB_DZDBX = ADM_DZ_4(betax,i,j,k) + DZDB_DZDBY = ADM_DZ_4(betay,i,j,k) + DZDB_DZDBZ = ADM_DZ_4(betaz,i,j,k) + end if #endif #ifdef CCODE diff --git a/src/macro/DZDCG_guts.h b/src/macro/DZDCG_guts.h index 82f3a09..08664d0 100644 --- a/src/macro/DZDCG_guts.h +++ b/src/macro/DZDCG_guts.h @@ -21,15 +21,22 @@ #ifdef FCODE - DZDCG_OO2DZ = 1D0/(2D0*DZDCG_DZ) - - DZDCG_DZDCGXX = DZDCG_OO2DZ*(DZDCG_GXX_KP - DZDCG_GXX_KM) - DZDCG_DZDCGXY = DZDCG_OO2DZ*(DZDCG_GXY_KP - DZDCG_GXY_KM) - DZDCG_DZDCGXZ = DZDCG_OO2DZ*(DZDCG_GXZ_KP - DZDCG_GXZ_KM) - DZDCG_DZDCGYY = DZDCG_OO2DZ*(DZDCG_GYY_KP - DZDCG_GYY_KM) - DZDCG_DZDCGYZ = DZDCG_OO2DZ*(DZDCG_GYZ_KP - DZDCG_GYZ_KM) - DZDCG_DZDCGZZ = DZDCG_OO2DZ*(DZDCG_GZZ_KP - DZDCG_GZZ_KM) - +#include "ADM_Derivative.h" + if (spatial_order.eq.2) then + DZDCG_DZDCGXX = ADM_DZ_2(gxx,i,j,k) + DZDCG_DZDCGXY = ADM_DZ_2(gxy,i,j,k) + DZDCG_DZDCGXZ = ADM_DZ_2(gxz,i,j,k) + DZDCG_DZDCGYY = ADM_DZ_2(gyy,i,j,k) + DZDCG_DZDCGYZ = ADM_DZ_2(gyz,i,j,k) + DZDCG_DZDCGZZ = ADM_DZ_2(gzz,i,j,k) + else + DZDCG_DZDCGXX = ADM_DZ_4(gxx,i,j,k) + DZDCG_DZDCGXY = ADM_DZ_4(gxy,i,j,k) + DZDCG_DZDCGXZ = ADM_DZ_4(gxz,i,j,k) + DZDCG_DZDCGYY = ADM_DZ_4(gyy,i,j,k) + DZDCG_DZDCGYZ = ADM_DZ_4(gyz,i,j,k) + DZDCG_DZDCGZZ = ADM_DZ_4(gzz,i,j,k) + end if #endif #ifdef CCODE diff --git a/src/macro/DZDK_guts.h b/src/macro/DZDK_guts.h index 1e205ab..cc8e5e6 100644 --- a/src/macro/DZDK_guts.h +++ b/src/macro/DZDK_guts.h @@ -13,15 +13,23 @@ #ifdef FCODE - DZDK_OO2DZ = 1D0/(2D0*DZDK_DZ) - - DZDK_DZDKXX = DZDK_OO2DZ*(DZDK_KXX_KP - DZDK_KXX_KM) - DZDK_DZDKXY = DZDK_OO2DZ*(DZDK_KXY_KP - DZDK_KXY_KM) - DZDK_DZDKXZ = DZDK_OO2DZ*(DZDK_KXZ_KP - DZDK_KXZ_KM) - DZDK_DZDKYY = DZDK_OO2DZ*(DZDK_KYY_KP - DZDK_KYY_KM) - DZDK_DZDKYZ = DZDK_OO2DZ*(DZDK_KYZ_KP - DZDK_KYZ_KM) - DZDK_DZDKZZ = DZDK_OO2DZ*(DZDK_KZZ_KP - DZDK_KZZ_KM) - +#include "ADM_Derivative.h" + + if (spatial_order.eq.2) then + DZDK_DZDKXX = ADM_DZ_2(kxx,i,j,k) + DZDK_DZDKXY = ADM_DZ_2(kxy,i,j,k) + DZDK_DZDKXZ = ADM_DZ_2(kxz,i,j,k) + DZDK_DZDKYY = ADM_DZ_2(kyy,i,j,k) + DZDK_DZDKYZ = ADM_DZ_2(kyz,i,j,k) + DZDK_DZDKZZ = ADM_DZ_2(kzz,i,j,k) + else + DZDK_DZDKXX = ADM_DZ_4(kxx,i,j,k) + DZDK_DZDKXY = ADM_DZ_4(kxy,i,j,k) + DZDK_DZDKXZ = ADM_DZ_4(kxz,i,j,k) + DZDK_DZDKYY = ADM_DZ_4(kyy,i,j,k) + DZDK_DZDKYZ = ADM_DZ_4(kyz,i,j,k) + DZDK_DZDKZZ = ADM_DZ_4(kzz,i,j,k) + end if #endif #ifdef CCODE diff --git a/src/macro/DZZDG_guts.h b/src/macro/DZZDG_guts.h index 601589b..853c12d 100644 --- a/src/macro/DZZDG_guts.h +++ b/src/macro/DZZDG_guts.h @@ -20,33 +20,53 @@ #ifdef FCODE - DZZDG_OODZ2 = 1D0/(dz*dz) +#include "ADM_Derivative.h" /* Factor involving 2nd derivative of conformal factor */ IF (conformal_state .eq. 0) THEN DZZDG_FAC = 0 ELSE - DZZDG_FAC = DZDG_PSI4*(4*DZZDG_DZZDPSI_O_PSI + 12*DZDG_DZDPSI_O_PSI*DZDG_DZDPSI_O_PSI) + DZZDG_FAC = DZDG_PSI4*(4*DZZDG_DZZDPSI_O_PSI \ + + 12*DZDG_DZDPSI_O_PSI*DZDG_DZDPSI_O_PSI) ENDIF - DZZDG_DZZDGXX = 2*DZDCG_DZDCGXX*DZDG_FAC+DZZDG_FAC*DZDG_GXX+DZDG_PSI4 \ - *DZZDG_OODZ2*(DZDCG_GXX_KP-2*DZDG_GXX+DZDCG_GXX_KM) + if (spatial_order.eq.2) then + DZZDG_DZZDGXX = 2*DZDCG_DZDCGXX*DZDG_FAC + DZZDG_FAC*DZDG_GXX \ + + DZDG_PSI4*ADM_DZZ_2(gxx,i,j,k) - DZZDG_DZZDGXY = 2*DZDCG_DZDCGXY*DZDG_FAC+DZZDG_FAC*DZDG_GXY+DZDG_PSI4 \ - *DZZDG_OODZ2*(DZDCG_GXY_KP-2*DZDG_GXY+DZDCG_GXY_KM) + DZZDG_DZZDGXY = 2*DZDCG_DZDCGXY*DZDG_FAC + DZZDG_FAC*DZDG_GXY \ + + DZDG_PSI4*ADM_DZZ_2(gxy,i,j,k) - DZZDG_DZZDGXZ = 2*DZDCG_DZDCGXZ*DZDG_FAC+DZZDG_FAC*DZDG_GXZ+DZDG_PSI4 \ - *DZZDG_OODZ2*(DZDCG_GXZ_KP-2*DZDG_GXZ+DZDCG_GXZ_KM) + DZZDG_DZZDGXZ = 2*DZDCG_DZDCGXZ*DZDG_FAC + DZZDG_FAC*DZDG_GXZ \ + + DZDG_PSI4*ADM_DZZ_2(gxz,i,j,k) - DZZDG_DZZDGYY = 2*DZDCG_DZDCGYY*DZDG_FAC+DZZDG_FAC*DZDG_GYY+DZDG_PSI4 \ - *DZZDG_OODZ2*(DZDCG_GYY_KP-2*DZDG_GYY+DZDCG_GYY_KM) + DZZDG_DZZDGYY = 2*DZDCG_DZDCGYY*DZDG_FAC + DZZDG_FAC*DZDG_GYY \ + + DZDG_PSI4*ADM_DZZ_2(gyy,i,j,k) - DZZDG_DZZDGYZ = 2*DZDCG_DZDCGYZ*DZDG_FAC+DZZDG_FAC*DZDG_GYZ+DZDG_PSI4 \ - *DZZDG_OODZ2*(DZDCG_GYZ_KP-2*DZDG_GYZ+DZDCG_GYZ_KM) + DZZDG_DZZDGYZ = 2*DZDCG_DZDCGYZ*DZDG_FAC + DZZDG_FAC*DZDG_GYZ \ + + DZDG_PSI4*ADM_DZZ_2(gyz,i,j,k) - DZZDG_DZZDGZZ = 2*DZDCG_DZDCGZZ*DZDG_FAC+DZZDG_FAC*DZDG_GZZ+DZDG_PSI4 \ - *DZZDG_OODZ2*(DZDCG_GZZ_KP-2*DZDG_GZZ+DZDCG_GZZ_KM) + DZZDG_DZZDGZZ = 2*DZDCG_DZDCGZZ*DZDG_FAC + DZZDG_FAC*DZDG_GZZ \ + + DZDG_PSI4*ADM_DZZ_2(gzz,i,j,k) + else + DZZDG_DZZDGXX = 2*DZDCG_DZDCGXX*DZDG_FAC + DZZDG_FAC*DZDG_GXX \ + + DZDG_PSI4*ADM_DZZ_4(gxx,i,j,k) + DZZDG_DZZDGXY = 2*DZDCG_DZDCGXY*DZDG_FAC + DZZDG_FAC*DZDG_GXY \ + + DZDG_PSI4*ADM_DZZ_4(gxy,i,j,k) + + DZZDG_DZZDGXZ = 2*DZDCG_DZDCGXZ*DZDG_FAC + DZZDG_FAC*DZDG_GXZ \ + + DZDG_PSI4*ADM_DZZ_4(gxz,i,j,k) + + DZZDG_DZZDGYY = 2*DZDCG_DZDCGYY*DZDG_FAC + DZZDG_FAC*DZDG_GYY \ + + DZDG_PSI4*ADM_DZZ_4(gyy,i,j,k) + + DZZDG_DZZDGYZ = 2*DZDCG_DZDCGYZ*DZDG_FAC + DZZDG_FAC*DZDG_GYZ \ + + DZDG_PSI4*ADM_DZZ_4(gyz,i,j,k) + + DZZDG_DZZDGZZ = 2*DZDCG_DZDCGZZ*DZDG_FAC + DZZDG_FAC*DZDG_GZZ \ + + DZDG_PSI4*ADM_DZZ_4(gzz,i,j,k) + end if #endif #ifdef CCODE diff --git a/src/macro/HAMADM_guts.h b/src/macro/HAMADM_guts.h index 523d0d5..9c87ccb 100644 --- a/src/macro/HAMADM_guts.h +++ b/src/macro/HAMADM_guts.h @@ -21,14 +21,12 @@ #ifdef FCODE HAMADM_HAMADM = TRRICCI_TRRICCI - TRKK_TRKK + TRK_TRK**2 - HAMADM_HAMADMABS = abs(TRRICCI_TRRICCI) + abs(TRKK_TRKK) + TRK_TRK**2 #endif #ifdef CCODE HAMADM_HAMADM = TRRICCI_TRRICCI - TRKK_TRKK + TRK_TRK**2 - HAMADM_HAMADMABS = fabs(TRRICCI_TRRICCI) + fabs(TRKK_TRKK) + TRK_TRK**2 #endif diff --git a/src/macro/LIEG_guts.h b/src/macro/LIEG_guts.h index f36c812..3afd945 100644 --- a/src/macro/LIEG_guts.h +++ b/src/macro/LIEG_guts.h @@ -21,6 +21,8 @@ #ifdef FCODE +#include "ADM_Derivative.h" + /* Advection. */ LIEG_LGXX = 0.0D0 @@ -43,63 +45,39 @@ else if (admmacros_advectionx .eq. 1) then /* UPWIND1 */ - LIEG_LGXX = LIEG_LGXX + LIEG_BX \ - *(gxx(i+1,j,k) - gxx(i,j,k))/dx - LIEG_LGYY = LIEG_LGYY + LIEG_BX \ - *(gyy(i+1,j,k) - gyy(i,j,k))/dx - LIEG_LGZZ = LIEG_LGZZ + LIEG_BX \ - *(gzz(i+1,j,k) - gzz(i,j,k))/dx - LIEG_LGXY = LIEG_LGXY + LIEG_BX \ - *(gxy(i+1,j,k) - gxy(i,j,k))/dx - LIEG_LGXZ = LIEG_LGXZ + LIEG_BX \ - *(gxz(i+1,j,k) - gxz(i,j,k))/dx - LIEG_LGYZ = LIEG_LGYZ + LIEG_BX \ - *(gyz(i+1,j,k) - gyz(i,j,k))/dx + LIEG_LGXX = LIEG_LGXX + ADM_ADV_DX_P1(gxx,i,j,k) + LIEG_LGYY = LIEG_LGYY + ADM_ADV_DX_P1(gyy,i,j,k) + LIEG_LGZZ = LIEG_LGZZ + ADM_ADV_DX_P1(gzz,i,j,k) + LIEG_LGXY = LIEG_LGXY + ADM_ADV_DX_P1(gxy,i,j,k) + LIEG_LGXZ = LIEG_LGXZ + ADM_ADV_DX_P1(gxz,i,j,k) + LIEG_LGYZ = LIEG_LGYZ + ADM_ADV_DX_P1(gyz,i,j,k) else if (admmacros_advectionx .eq. -1) then /* UPWIND1 */ - LIEG_LGXX = LIEG_LGXX + LIEG_BX \ - *(gxx(i,j,k) - gxx(i-1,j,k))/dx - LIEG_LGYY = LIEG_LGYY + LIEG_BX \ - *(gyy(i,j,k) - gyy(i-1,j,k))/dx - LIEG_LGZZ = LIEG_LGZZ + LIEG_BX \ - *(gzz(i,j,k) - gzz(i-1,j,k))/dx - LIEG_LGXY = LIEG_LGXY + LIEG_BX \ - *(gxy(i,j,k) - gxy(i-1,j,k))/dx - LIEG_LGXZ = LIEG_LGXZ + LIEG_BX \ - *(gxz(i,j,k) - gxz(i-1,j,k))/dx - LIEG_LGYZ = LIEG_LGYZ + LIEG_BX \ - *(gyz(i,j,k) - gyz(i-1,j,k))/dx + LIEG_LGXX = LIEG_LGXX + ADM_ADV_DX_M1(gxx,i,j,k) + LIEG_LGYY = LIEG_LGYY + ADM_ADV_DX_M1(gyy,i,j,k) + LIEG_LGZZ = LIEG_LGZZ + ADM_ADV_DX_M1(gzz,i,j,k) + LIEG_LGXY = LIEG_LGXY + ADM_ADV_DX_M1(gxy,i,j,k) + LIEG_LGXZ = LIEG_LGXZ + ADM_ADV_DX_M1(gxz,i,j,k) + LIEG_LGYZ = LIEG_LGYZ + ADM_ADV_DX_M1(gyz,i,j,k) else if (admmacros_advectionx .eq. 2) then /* UPWIND2 */ - LIEG_LGXX = LIEG_LGXX - 0.5D0*LIEG_BX/dx \ - *(3.0D0*gxx(i,j,k) - 4.0D0*gxx(i+1,j,k) + gxx(i+2,j,k)) - LIEG_LGYY = LIEG_LGYY - 0.5D0*LIEG_BX/dx \ - *(3.0D0*gyy(i,j,k) - 4.0D0*gyy(i+1,j,k) + gyy(i+2,j,k)) - LIEG_LGZZ = LIEG_LGZZ - 0.5D0*LIEG_BX/dx \ - *(3.0D0*gzz(i,j,k) - 4.0D0*gzz(i+1,j,k) + gzz(i+2,j,k)) - LIEG_LGXY = LIEG_LGXY - 0.5D0*LIEG_BX/dx \ - *(3.0D0*gxy(i,j,k) - 4.0D0*gxy(i+1,j,k) + gxy(i+2,j,k)) - LIEG_LGXZ = LIEG_LGXZ - 0.5D0*LIEG_BX/dx \ - *(3.0D0*gxz(i,j,k) - 4.0D0*gxz(i+1,j,k) + gxz(i+2,j,k)) - LIEG_LGYZ = LIEG_LGYZ - 0.5D0*LIEG_BX/dx \ - *(3.0D0*gyz(i,j,k) - 4.0D0*gyz(i+1,j,k) + gyz(i+2,j,k)) + LIEG_LGXX = LIEG_LGXX + ADM_ADV_DX_P2(gxx,i,j,k) + LIEG_LGYY = LIEG_LGYY + ADM_ADV_DX_P2(gyy,i,j,k) + LIEG_LGZZ = LIEG_LGZZ + ADM_ADV_DX_P2(gzz,i,j,k) + LIEG_LGXY = LIEG_LGXY + ADM_ADV_DX_P2(gxy,i,j,k) + LIEG_LGXZ = LIEG_LGXZ + ADM_ADV_DX_P2(gxz,i,j,k) + LIEG_LGYZ = LIEG_LGYZ + ADM_ADV_DX_P2(gyz,i,j,k) else if (admmacros_advectionx .eq. -2) then /* UPWIND2 */ - LIEG_LGXX = LIEG_LGXX + 0.5D0*LIEG_BX/dx \ - *(3.0D0*gxx(i,j,k) - 4.0D0*gxx(i-1,j,k) + gxx(i-2,j,k)) - LIEG_LGYY = LIEG_LGYY + 0.5D0*LIEG_BX/dx \ - *(3.0D0*gyy(i,j,k) - 4.0D0*gyy(i-1,j,k) + gyy(i-2,j,k)) - LIEG_LGZZ = LIEG_LGZZ + 0.5D0*LIEG_BX/dx \ - *(3.0D0*gzz(i,j,k) - 4.0D0*gzz(i-1,j,k) + gzz(i-2,j,k)) - LIEG_LGXY = LIEG_LGXY + 0.5D0*LIEG_BX/dx \ - *(3.0D0*gxy(i,j,k) - 4.0D0*gxy(i-1,j,k) + gxy(i-2,j,k)) - LIEG_LGXZ = LIEG_LGXZ + 0.5D0*LIEG_BX/dx \ - *(3.0D0*gxz(i,j,k) - 4.0D0*gxz(i-1,j,k) + gxz(i-2,j,k)) - LIEG_LGYZ = LIEG_LGYZ + 0.5D0*LIEG_BX/dx \ - *(3.0D0*gyz(i,j,k) - 4.0D0*gyz(i-1,j,k) + gyz(i-2,j,k)) + LIEG_LGXX = LIEG_LGXX + ADM_ADV_DX_M2(gxx,i,j,k) + LIEG_LGYY = LIEG_LGYY + ADM_ADV_DX_M2(gyy,i,j,k) + LIEG_LGZZ = LIEG_LGZZ + ADM_ADV_DX_M2(gzz,i,j,k) + LIEG_LGXY = LIEG_LGXY + ADM_ADV_DX_M2(gxy,i,j,k) + LIEG_LGXZ = LIEG_LGXZ + ADM_ADV_DX_M2(gxz,i,j,k) + LIEG_LGYZ = LIEG_LGYZ + ADM_ADV_DX_M2(gyz,i,j,k) end if @@ -117,63 +95,39 @@ else if (admmacros_advectiony .eq. 1) then /* UPWIND1 */ - LIEG_LGXX = LIEG_LGXX + LIEG_BY \ - *(gxx(i,j+1,k) - gxx(i,j,k))/dy - LIEG_LGYY = LIEG_LGYY + LIEG_BY \ - *(gyy(i,j+1,k) - gyy(i,j,k))/dy - LIEG_LGZZ = LIEG_LGZZ + LIEG_BY \ - *(gzz(i,j+1,k) - gzz(i,j,k))/dy - LIEG_LGXY = LIEG_LGXY + LIEG_BY \ - *(gxy(i,j+1,k) - gxy(i,j,k))/dy - LIEG_LGXZ = LIEG_LGXZ + LIEG_BY \ - *(gxz(i,j+1,k) - gxz(i,j,k))/dy - LIEG_LGYZ = LIEG_LGYZ + LIEG_BY \ - *(gyz(i,j+1,k) - gyz(i,j,k))/dy + LIEG_LGXX = LIEG_LGXX + ADM_ADV_DY_P1(gxx,i,j,k) + LIEG_LGYY = LIEG_LGYY + ADM_ADV_DY_P1(gyy,i,j,k) + LIEG_LGZZ = LIEG_LGZZ + ADM_ADV_DY_P1(gzz,i,j,k) + LIEG_LGXY = LIEG_LGXY + ADM_ADV_DY_P1(gxy,i,j,k) + LIEG_LGXZ = LIEG_LGXZ + ADM_ADV_DY_P1(gxz,i,j,k) + LIEG_LGYZ = LIEG_LGYZ + ADM_ADV_DY_P1(gyz,i,j,k) else if (admmacros_advectiony .eq. -1) then /* UPWIND1 */ - LIEG_LGXX = LIEG_LGXX + LIEG_BY \ - *(gxx(i,j,k) - gxx(i,j-1,k))/dy - LIEG_LGYY = LIEG_LGYY + LIEG_BY \ - *(gyy(i,j,k) - gyy(i,j-1,k))/dy - LIEG_LGZZ = LIEG_LGZZ + LIEG_BY \ - *(gzz(i,j,k) - gzz(i,j-1,k))/dy - LIEG_LGXY = LIEG_LGXY + LIEG_BY \ - *(gxy(i,j,k) - gxy(i,j-1,k))/dy - LIEG_LGXZ = LIEG_LGXZ + LIEG_BY \ - *(gxz(i,j,k) - gxz(i,j-1,k))/dy - LIEG_LGYZ = LIEG_LGYZ + LIEG_BY \ - *(gyz(i,j,k) - gyz(i,j-1,k))/dy + LIEG_LGXX = LIEG_LGXX + ADM_ADV_DY_M1(gxx,i,j,k) + LIEG_LGYY = LIEG_LGYY + ADM_ADV_DY_M1(gyy,i,j,k) + LIEG_LGZZ = LIEG_LGZZ + ADM_ADV_DY_M1(gzz,i,j,k) + LIEG_LGXY = LIEG_LGXY + ADM_ADV_DY_M1(gxy,i,j,k) + LIEG_LGXZ = LIEG_LGXZ + ADM_ADV_DY_M1(gxz,i,j,k) + LIEG_LGYZ = LIEG_LGYZ + ADM_ADV_DY_M1(gyz,i,j,k) else if (admmacros_advectiony .eq. 2) then /* UPWIND2 */ - LIEG_LGXX = LIEG_LGXX - 0.5D0*LIEG_BY/dy \ - *(3.0D0*gxx(i,j,k) - 4.0D0*gxx(i,j+1,k) + gxx(i,j+2,k)) - LIEG_LGYY = LIEG_LGYY - 0.5D0*LIEG_BY/dy \ - *(3.0D0*gyy(i,j,k) - 4.0D0*gyy(i,j+1,k) + gyy(i,j+2,k)) - LIEG_LGZZ = LIEG_LGZZ - 0.5D0*LIEG_BY/dy \ - *(3.0D0*gzz(i,j,k) - 4.0D0*gzz(i,j+1,k) + gzz(i,j+2,k)) - LIEG_LGXY = LIEG_LGXY - 0.5D0*LIEG_BY/dy \ - *(3.0D0*gxy(i,j,k) - 4.0D0*gxy(i,j+1,k) + gxy(i,j+2,k)) - LIEG_LGXZ = LIEG_LGXZ - 0.5D0*LIEG_BY/dy \ - *(3.0D0*gxz(i,j,k) - 4.0D0*gxz(i,j+1,k) + gxz(i,j+2,k)) - LIEG_LGYZ = LIEG_LGYZ - 0.5D0*LIEG_BY/dy \ - *(3.0D0*gyz(i,j,k) - 4.0D0*gyz(i,j+1,k) + gyz(i,j+2,k)) + LIEG_LGXX = LIEG_LGXX + ADM_ADV_DY_P2(gxx,i,j,k) + LIEG_LGYY = LIEG_LGYY + ADM_ADV_DY_P2(gyy,i,j,k) + LIEG_LGZZ = LIEG_LGZZ + ADM_ADV_DY_P2(gzz,i,j,k) + LIEG_LGXY = LIEG_LGXY + ADM_ADV_DY_P2(gxy,i,j,k) + LIEG_LGXZ = LIEG_LGXZ + ADM_ADV_DY_P2(gxz,i,j,k) + LIEG_LGYZ = LIEG_LGYZ + ADM_ADV_DY_P2(gyz,i,j,k) else if (admmacros_advectiony .eq. -2) then /* UPWIND2 */ - LIEG_LGXX = LIEG_LGXX + 0.5D0*LIEG_BY/dy \ - *(3.0D0*gxx(i,j,k) - 4.0D0*gxx(i,j-1,k) + gxx(i,j-2,k)) - LIEG_LGYY = LIEG_LGYY + 0.5D0*LIEG_BY/dy \ - *(3.0D0*gyy(i,j,k) - 4.0D0*gyy(i,j-1,k) + gyy(i,j-2,k)) - LIEG_LGZZ = LIEG_LGZZ + 0.5D0*LIEG_BY/dy \ - *(3.0D0*gzz(i,j,k) - 4.0D0*gzz(i,j-1,k) + gzz(i,j-2,k)) - LIEG_LGXY = LIEG_LGXY + 0.5D0*LIEG_BY/dy \ - *(3.0D0*gxy(i,j,k) - 4.0D0*gxy(i,j-1,k) + gxy(i,j-2,k)) - LIEG_LGXZ = LIEG_LGXZ + 0.5D0*LIEG_BY/dy \ - *(3.0D0*gxz(i,j,k) - 4.0D0*gxz(i,j-1,k) + gxz(i,j-2,k)) - LIEG_LGYZ = LIEG_LGYZ + 0.5D0*LIEG_BY/dy \ - *(3.0D0*gyz(i,j,k) - 4.0D0*gyz(i,j-1,k) + gyz(i,j-2,k)) + LIEG_LGXX = LIEG_LGXX + ADM_ADV_DY_M2(gxx,i,j,k) + LIEG_LGYY = LIEG_LGYY + ADM_ADV_DY_M2(gyy,i,j,k) + LIEG_LGZZ = LIEG_LGZZ + ADM_ADV_DY_M2(gzz,i,j,k) + LIEG_LGXY = LIEG_LGXY + ADM_ADV_DY_M2(gxy,i,j,k) + LIEG_LGXZ = LIEG_LGXZ + ADM_ADV_DY_M2(gxz,i,j,k) + LIEG_LGYZ = LIEG_LGYZ + ADM_ADV_DY_M2(gyz,i,j,k) end if @@ -191,73 +145,39 @@ else if ((admmacros_advectionz .eq. 1)) then /* UPWIND1 */ - LIEG_LGXX = LIEG_LGXX + LIEG_BZ \ - *(gxx(i,j,k+1) - gxx(i,j,k))/dz - LIEG_LGYY = LIEG_LGYY + LIEG_BZ \ - *(gyy(i,j,k+1) - gyy(i,j,k))/dz - LIEG_LGZZ = LIEG_LGZZ + LIEG_BZ \ - *(gzz(i,j,k+1) - gzz(i,j,k))/dz - LIEG_LGXY = LIEG_LGXY + LIEG_BZ \ - *(gxy(i,j,k+1) - gxy(i,j,k))/dz - LIEG_LGXZ = LIEG_LGXZ + LIEG_BZ \ - *(gxz(i,j,k+1) - gxz(i,j,k))/dz - LIEG_LGYZ = LIEG_LGYZ + LIEG_BZ \ - *(gyz(i,j,k+1) - gyz(i,j,k))/dz + LIEG_LGXX = LIEG_LGXX + ADM_ADV_DZ_P1(gxx,i,j,k) + LIEG_LGYY = LIEG_LGYY + ADM_ADV_DZ_P1(gyy,i,j,k) + LIEG_LGZZ = LIEG_LGZZ + ADM_ADV_DZ_P1(gzz,i,j,k) + LIEG_LGXY = LIEG_LGXY + ADM_ADV_DZ_P1(gxy,i,j,k) + LIEG_LGXZ = LIEG_LGXZ + ADM_ADV_DZ_P1(gxz,i,j,k) + LIEG_LGYZ = LIEG_LGYZ + ADM_ADV_DZ_P1(gyz,i,j,k) else if ((admmacros_advectionz .eq. -1)) then /* UPWIND1 */ - LIEG_LGXX = LIEG_LGXX + LIEG_BZ \ - *(gxx(i,j,k) - gxx(i,j,k-1))/dz - LIEG_LGYY = LIEG_LGYY + LIEG_BZ \ - *(gyy(i,j,k) - gyy(i,j,k-1))/dz - LIEG_LGZZ = LIEG_LGZZ + LIEG_BZ \ - *(gzz(i,j,k) - gzz(i,j,k-1))/dz - LIEG_LGXY = LIEG_LGXY + LIEG_BZ \ - *(gxy(i,j,k) - gxy(i,j,k-1))/dz - LIEG_LGXZ = LIEG_LGXZ + LIEG_BZ \ - *(gxz(i,j,k) - gxz(i,j,k-1))/dz - LIEG_LGYZ = LIEG_LGYZ + LIEG_BZ \ - *(gyz(i,j,k) - gyz(i,j,k-1))/dz + LIEG_LGXX = LIEG_LGXX + ADM_ADV_DZ_M1(gxx,i,j,k) + LIEG_LGYY = LIEG_LGYY + ADM_ADV_DZ_M1(gyy,i,j,k) + LIEG_LGZZ = LIEG_LGZZ + ADM_ADV_DZ_M1(gzz,i,j,k) + LIEG_LGXY = LIEG_LGXY + ADM_ADV_DZ_M1(gxy,i,j,k) + LIEG_LGXZ = LIEG_LGXZ + ADM_ADV_DZ_M1(gxz,i,j,k) + LIEG_LGYZ = LIEG_LGYZ + ADM_ADV_DZ_M1(gyz,i,j,k) else if (admmacros_advectionz .eq. 2) then /* UPWIND2 */ - LIEG_LGXX = LIEG_LGXX - 0.5D0*LIEG_BZ/dz \ - *(3.0D0*gxx(i,j,k) - 4.0D0*gxx(i,j,k+1) + gxx(i,j,k+2)) - - LIEG_LGYY = LIEG_LGYY - 0.5D0*LIEG_BZ/dz \ - *(3.0D0*gyy(i,j,k) - 4.0D0*gyy(i,j,k+1) + gyy(i,j,k+2)) - - LIEG_LGZZ = LIEG_LGZZ - 0.5D0*LIEG_BZ/dz \ - *(3.0D0*gzz(i,j,k) - 4.0D0*gzz(i,j,k+1) + gzz(i,j,k+2)) - - LIEG_LGXY = LIEG_LGXY - 0.5D0*LIEG_BZ/dz \ - *(3.0D0*gxy(i,j,k) - 4.0D0*gxy(i,j,k+1) + gxy(i,j,k+2)) - - LIEG_LGXZ = LIEG_LGXZ - 0.5D0*LIEG_BZ/dz \ - *(3.0D0*gxz(i,j,k) - 4.0D0*gxz(i,j,k+1) + gxz(i,j,k+2)) - - LIEG_LGYZ = LIEG_LGYZ - 0.5D0*LIEG_BZ/dz \ - *(3.0D0*gyz(i,j,k) - 4.0D0*gyz(i,j,k+1) + gyz(i,j,k+2)) + LIEG_LGXX = LIEG_LGXX + ADM_ADV_DZ_P2(gxx,i,j,k) + LIEG_LGYY = LIEG_LGYY + ADM_ADV_DZ_P2(gyy,i,j,k) + LIEG_LGZZ = LIEG_LGZZ + ADM_ADV_DZ_P2(gzz,i,j,k) + LIEG_LGXY = LIEG_LGXY + ADM_ADV_DZ_P2(gxy,i,j,k) + LIEG_LGXZ = LIEG_LGXZ + ADM_ADV_DZ_P2(gxz,i,j,k) + LIEG_LGYZ = LIEG_LGYZ + ADM_ADV_DZ_P2(gyz,i,j,k) else if (admmacros_advectionz .eq. -2) then /* UPWIND2 */ - LIEG_LGXX = LIEG_LGXX + 0.5D0*LIEG_BZ/dz \ - *(3.0D0*gxx(i,j,k) - 4.0D0*gxx(i,j,k-1) + gxx(i,j,k-2)) - - LIEG_LGYY = LIEG_LGYY + 0.5D0*LIEG_BZ/dz \ - *(3.0D0*gyy(i,j,k) - 4.0D0*gyy(i,j,k-1) + gyy(i,j,k-2)) - - LIEG_LGZZ = LIEG_LGZZ + 0.5D0*LIEG_BZ/dz \ - *(3.0D0*gzz(i,j,k) - 4.0D0*gzz(i,j,k-1) + gzz(i,j,k-2)) - - LIEG_LGXY = LIEG_LGXY + 0.5D0*LIEG_BZ/dz \ - *(3.0D0*gxy(i,j,k) - 4.0D0*gxy(i,j,k-1) + gxy(i,j,k-2)) - - LIEG_LGXZ = LIEG_LGXZ + 0.5D0*LIEG_BZ/dz \ - *(3.0D0*gxz(i,j,k) - 4.0D0*gxz(i,j,k-1) + gxz(i,j,k-2)) - - LIEG_LGYZ = LIEG_LGYZ + 0.5D0*LIEG_BZ/dz \ - *(3.0D0*gyz(i,j,k) - 4.0D0*gyz(i,j,k-1) + gyz(i,j,k-2)) + LIEG_LGXX = LIEG_LGXX + ADM_ADV_DZ_M2(gxx,i,j,k) + LIEG_LGYY = LIEG_LGYY + ADM_ADV_DZ_M2(gyy,i,j,k) + LIEG_LGZZ = LIEG_LGZZ + ADM_ADV_DZ_M2(gzz,i,j,k) + LIEG_LGXY = LIEG_LGXY + ADM_ADV_DZ_M2(gxy,i,j,k) + LIEG_LGXZ = LIEG_LGXZ + ADM_ADV_DZ_M2(gxz,i,j,k) + LIEG_LGYZ = LIEG_LGYZ + ADM_ADV_DZ_M2(gyz,i,j,k) end if diff --git a/src/macro/LIEK_guts.h b/src/macro/LIEK_guts.h index b18b3f9..bdb6141 100644 --- a/src/macro/LIEK_guts.h +++ b/src/macro/LIEK_guts.h @@ -22,6 +22,8 @@ #ifdef FCODE +#include "ADM_Derivative.h" + /* Advection. */ LIEK_LKXX = 0.0D0 @@ -34,96 +36,47 @@ if (admmacros_advectionx .eq. 0) then /* CENTER X */ LIEK_LKXX = LIEK_LKXX + DXDK_DXDKXX*LIEK_BX - LIEK_LKYY = LIEK_LKYY + DXDK_DXDKYY*LIEK_BX - LIEK_LKZZ = LIEK_LKZZ + DXDK_DXDKZZ*LIEK_BX - LIEK_LKXY = LIEK_LKXY + DXDK_DXDKXY*LIEK_BX - LIEK_LKXZ = LIEK_LKXZ + DXDK_DXDKXZ*LIEK_BX - LIEK_LKYZ = LIEK_LKYZ + DXDK_DXDKYZ*LIEK_BX else if (admmacros_advectionx .eq. 1) then /* UPWIND1 X */ - LIEK_LKXX = LIEK_LKXX + LIEK_BX \ - *(kxx(i+1,j,k) - kxx(i,j,k))/dx - - LIEK_LKYY = LIEK_LKYY + LIEK_BX \ - *(kyy(i+1,j,k) - kyy(i,j,k))/dx - - LIEK_LKZZ = LIEK_LKZZ + LIEK_BX \ - *(kzz(i+1,j,k) - kzz(i,j,k))/dx - - LIEK_LKXY = LIEK_LKXY + LIEK_BX \ - *(kxy(i+1,j,k) - kxy(i,j,k))/dx - - LIEK_LKXZ = LIEK_LKXZ + LIEK_BX \ - *(kxz(i+1,j,k) - kxz(i,j,k))/dx - - LIEK_LKYZ = LIEK_LKYZ + LIEK_BX \ - *(kyz(i+1,j,k) - kyz(i,j,k))/dx + LIEK_LKXX = LIEK_LKXX + ADM_ADV_DX_P1(kxx,i,j,k) + LIEK_LKYY = LIEK_LKYY + ADM_ADV_DX_P1(kyy,i,j,k) + LIEK_LKZZ = LIEK_LKZZ + ADM_ADV_DX_P1(kzz,i,j,k) + LIEK_LKXY = LIEK_LKXY + ADM_ADV_DX_P1(kxy,i,j,k) + LIEK_LKXZ = LIEK_LKXZ + ADM_ADV_DX_P1(kxz,i,j,k) + LIEK_LKYZ = LIEK_LKYZ + ADM_ADV_DX_P1(kyz,i,j,k) else if (admmacros_advectionx .eq. -1) then /* UPWIND1 X */ - LIEK_LKXX = LIEK_LKXX + LIEK_BX \ - *(kxx(i,j,k) - kxx(i-1,j,k))/dx - - LIEK_LKYY = LIEK_LKYY + LIEK_BX \ - *(kyy(i,j,k) - kyy(i-1,j,k))/dx - - LIEK_LKZZ = LIEK_LKZZ + LIEK_BX \ - *(kzz(i,j,k) - kzz(i-1,j,k))/dx - - LIEK_LKXY = LIEK_LKXY + LIEK_BX \ - *(kxy(i,j,k) - kxy(i-1,j,k))/dx - - LIEK_LKXZ = LIEK_LKXZ + LIEK_BX \ - *(kxz(i,j,k) - kxz(i-1,j,k))/dx - - LIEK_LKYZ = LIEK_LKYZ + LIEK_BX \ - *(kyz(i,j,k) - kyz(i-1,j,k))/dx + LIEK_LKXX = LIEK_LKXX + ADM_ADV_DX_M1(kxx,i,j,k) + LIEK_LKYY = LIEK_LKYY + ADM_ADV_DX_M1(kyy,i,j,k) + LIEK_LKZZ = LIEK_LKZZ + ADM_ADV_DX_M1(kzz,i,j,k) + LIEK_LKXY = LIEK_LKXY + ADM_ADV_DX_M1(kxy,i,j,k) + LIEK_LKXZ = LIEK_LKXZ + ADM_ADV_DX_M1(kxz,i,j,k) + LIEK_LKYZ = LIEK_LKYZ + ADM_ADV_DX_M1(kyz,i,j,k) else if (admmacros_advectionx .eq. 2) then /* UPWIND2 X */ - LIEK_LKXX = LIEK_LKXX - 0.5D0*LIEK_BX/dx \ - *(3.0D0*kxx(i,j,k) - 4.0D0*kxx(i+1,j,k) + kxx(i+2,j,k)) - - LIEK_LKYY = LIEK_LKYY - 0.5D0*LIEK_BX/dx \ - *(3.0D0*kyy(i,j,k) - 4.0D0*kyy(i+1,j,k) + kyy(i+2,j,k)) - - LIEK_LKZZ = LIEK_LKZZ - 0.5D0*LIEK_BX/dx \ - *(3.0D0*kzz(i,j,k) - 4.0D0*kzz(i+1,j,k) + kzz(i+2,j,k)) - - LIEK_LKXY = LIEK_LKXY - 0.5D0*LIEK_BX/dx \ - *(3.0D0*kxy(i,j,k) - 4.0D0*kxy(i+1,j,k) + kxy(i+2,j,k)) - - LIEK_LKXZ = LIEK_LKXZ - 0.5D0*LIEK_BX/dx \ - *(3.0D0*kxz(i,j,k) - 4.0D0*kxz(i+1,j,k) + kxz(i+2,j,k)) - - LIEK_LKYZ = LIEK_LKYZ - 0.5D0*LIEK_BX/dx \ - *(3.0D0*kyz(i,j,k) - 4.0D0*kyz(i+1,j,k) + kyz(i+2,j,k)) + LIEK_LKXX = LIEK_LKXX + ADM_ADV_DX_P2(kxx,i,j,k) + LIEK_LKYY = LIEK_LKYY + ADM_ADV_DX_P2(kyy,i,j,k) + LIEK_LKZZ = LIEK_LKZZ + ADM_ADV_DX_P2(kzz,i,j,k) + LIEK_LKXY = LIEK_LKXY + ADM_ADV_DX_P2(kxy,i,j,k) + LIEK_LKXZ = LIEK_LKXZ + ADM_ADV_DX_P2(kxz,i,j,k) + LIEK_LKYZ = LIEK_LKYZ + ADM_ADV_DX_P2(kyz,i,j,k) else if (admmacros_advectionx .eq. -2) then /* UPWIND2 X */ - LIEK_LKXX = LIEK_LKXX + 0.5D0*LIEK_BX/dx \ - *(3.0D0*kxx(i,j,k) - 4.0D0*kxx(i-1,j,k) + kxx(i-2,j,k)) - - LIEK_LKYY = LIEK_LKYY + 0.5D0*LIEK_BX/dx \ - *(3.0D0*kyy(i,j,k) - 4.0D0*kyy(i-1,j,k) + kyy(i-2,j,k)) - - LIEK_LKZZ = LIEK_LKZZ + 0.5D0*LIEK_BX/dx \ - *(3.0D0*kzz(i,j,k) - 4.0D0*kzz(i-1,j,k) + kzz(i-2,j,k)) - - LIEK_LKXY = LIEK_LKXY + 0.5D0*LIEK_BX/dx \ - *(3.0D0*kxy(i,j,k) - 4.0D0*kxy(i-1,j,k) + kxy(i-2,j,k)) - - LIEK_LKXZ = LIEK_LKXZ + 0.5D0*LIEK_BX/dx \ - *(3.0D0*kxz(i,j,k) - 4.0D0*kxz(i-1,j,k) + kxz(i-2,j,k)) - - LIEK_LKYZ = LIEK_LKYZ + 0.5D0*LIEK_BX/dx \ - *(3.0D0*kyz(i,j,k) - 4.0D0*kyz(i-1,j,k) + kyz(i-2,j,k)) + LIEK_LKXX = LIEK_LKXX + ADM_ADV_DX_M2(kxx,i,j,k) + LIEK_LKYY = LIEK_LKYY + ADM_ADV_DX_M2(kyy,i,j,k) + LIEK_LKZZ = LIEK_LKZZ + ADM_ADV_DX_M2(kzz,i,j,k) + LIEK_LKXY = LIEK_LKXY + ADM_ADV_DX_M2(kxy,i,j,k) + LIEK_LKXZ = LIEK_LKXZ + ADM_ADV_DX_M2(kxz,i,j,k) + LIEK_LKYZ = LIEK_LKYZ + ADM_ADV_DX_M2(kyz,i,j,k) end if @@ -131,96 +84,47 @@ if (admmacros_advectiony .eq. 0) then /* CENTER Y */ LIEK_LKXX = LIEK_LKXX + DYDK_DYDKXX*LIEK_BY - LIEK_LKYY = LIEK_LKYY + DYDK_DYDKYY*LIEK_BY - LIEK_LKZZ = LIEK_LKZZ + DYDK_DYDKZZ*LIEK_BY - LIEK_LKXY = LIEK_LKXY + DYDK_DYDKXY*LIEK_BY - LIEK_LKXZ = LIEK_LKXZ + DYDK_DYDKXZ*LIEK_BY - LIEK_LKYZ = LIEK_LKYZ + DYDK_DYDKYZ*LIEK_BY else if (admmacros_advectiony .eq. 1) then /* UPWIND1 Y */ - LIEK_LKXX = LIEK_LKXX + LIEK_BY \ - *(kxx(i,j+1,k) - kxx(i,j,k))/dy - - LIEK_LKYY = LIEK_LKYY + LIEK_BY \ - *(kyy(i,j+1,k) - kyy(i,j,k))/dy - - LIEK_LKZZ = LIEK_LKZZ + LIEK_BY \ - *(kzz(i,j+1,k) - kzz(i,j,k))/dy - - LIEK_LKXY = LIEK_LKXY + LIEK_BY \ - *(kxy(i,j+1,k) - kxy(i,j,k))/dy - - LIEK_LKXZ = LIEK_LKXZ + LIEK_BY \ - *(kxz(i,j+1,k) - kxz(i,j,k))/dy - - LIEK_LKYZ = LIEK_LKYZ + LIEK_BY \ - *(kyz(i,j+1,k) - kyz(i,j,k))/dy + LIEK_LKXX = LIEK_LKXX + ADM_ADV_DY_P1(kxx,i,j,k) + LIEK_LKYY = LIEK_LKYY + ADM_ADV_DY_P1(kyy,i,j,k) + LIEK_LKZZ = LIEK_LKZZ + ADM_ADV_DY_P1(kzz,i,j,k) + LIEK_LKXY = LIEK_LKXY + ADM_ADV_DY_P1(kxy,i,j,k) + LIEK_LKXZ = LIEK_LKXZ + ADM_ADV_DY_P1(kxz,i,j,k) + LIEK_LKYZ = LIEK_LKYZ + ADM_ADV_DY_P1(kyz,i,j,k) else if (admmacros_advectiony .eq. -1) then /* UPWIND1 Y */ - LIEK_LKXX = LIEK_LKXX + LIEK_BY \ - *(kxx(i,j,k) - kxx(i,j-1,k))/dy - - LIEK_LKYY = LIEK_LKYY + LIEK_BY \ - *(kyy(i,j,k) - kyy(i,j-1,k))/dy - - LIEK_LKZZ = LIEK_LKZZ + LIEK_BY \ - *(kzz(i,j,k) - kzz(i,j-1,k))/dy - - LIEK_LKXY = LIEK_LKXY + LIEK_BY \ - *(kxy(i,j,k) - kxy(i,j-1,k))/dy - - LIEK_LKXZ = LIEK_LKXZ + LIEK_BY \ - *(kxz(i,j,k) - kxz(i,j-1,k))/dy - - LIEK_LKYZ = LIEK_LKYZ + LIEK_BY \ - *(kyz(i,j,k) - kyz(i,j-1,k))/dy + LIEK_LKXX = LIEK_LKXX + ADM_ADV_DY_M1(kxx,i,j,k) + LIEK_LKYY = LIEK_LKYY + ADM_ADV_DY_M1(kyy,i,j,k) + LIEK_LKZZ = LIEK_LKZZ + ADM_ADV_DY_M1(kzz,i,j,k) + LIEK_LKXY = LIEK_LKXY + ADM_ADV_DY_M1(kxy,i,j,k) + LIEK_LKXZ = LIEK_LKXZ + ADM_ADV_DY_M1(kxz,i,j,k) + LIEK_LKYZ = LIEK_LKYZ + ADM_ADV_DY_M1(kyz,i,j,k) else if (admmacros_advectiony .eq. 2) then /* UPWIND2 Y */ - LIEK_LKXX = LIEK_LKXX - 0.5D0*LIEK_BY/dy \ - *(3.0D0*kxx(i,j,k) - 4.0D0*kxx(i,j+1,k) + kxx(i,j+2,k)) - - LIEK_LKYY = LIEK_LKYY - 0.5D0*LIEK_BY/dy \ - *(3.0D0*kyy(i,j,k) - 4.0D0*kyy(i,j+1,k) + kyy(i,j+2,k)) - - LIEK_LKZZ = LIEK_LKZZ - 0.5D0*LIEK_BY/dy \ - *(3.0D0*kzz(i,j,k) - 4.0D0*kzz(i,j+1,k) + kzz(i,j+2,k)) - - LIEK_LKXY = LIEK_LKXY - 0.5D0*LIEK_BY/dy \ - *(3.0D0*kxy(i,j,k) - 4.0D0*kxy(i,j+1,k) + kxy(i,j+2,k)) - - LIEK_LKXZ = LIEK_LKXZ - 0.5D0*LIEK_BY/dy \ - *(3.0D0*kxz(i,j,k) - 4.0D0*kxz(i,j+1,k) + kxz(i,j+2,k)) - - LIEK_LKYZ = LIEK_LKYZ - 0.5D0*LIEK_BY/dy \ - *(3.0D0*kyz(i,j,k) - 4.0D0*kyz(i,j+1,k) + kyz(i,j+2,k)) + LIEK_LKXX = LIEK_LKXX + ADM_ADV_DY_P2(kxx,i,j,k) + LIEK_LKYY = LIEK_LKYY + ADM_ADV_DY_P2(kyy,i,j,k) + LIEK_LKZZ = LIEK_LKZZ + ADM_ADV_DY_P2(kzz,i,j,k) + LIEK_LKXY = LIEK_LKXY + ADM_ADV_DY_P2(kxy,i,j,k) + LIEK_LKXZ = LIEK_LKXZ + ADM_ADV_DY_P2(kxz,i,j,k) + LIEK_LKYZ = LIEK_LKYZ + ADM_ADV_DY_P2(kyz,i,j,k) else if (admmacros_advectiony .eq. -2) then /* UPWIND2 Y */ - LIEK_LKXX = LIEK_LKXX + 0.5D0*LIEK_BY/dy \ - *(3.0D0*kxx(i,j,k) - 4.0D0*kxx(i,j-1,k) + kxx(i,j-2,k)) - - LIEK_LKYY = LIEK_LKYY + 0.5D0*LIEK_BY/dy \ - *(3.0D0*kyy(i,j,k) - 4.0D0*kyy(i,j-1,k) + kyy(i,j-2,k)) - - LIEK_LKZZ = LIEK_LKZZ + 0.5D0*LIEK_BY/dy \ - *(3.0D0*kzz(i,j,k) - 4.0D0*kzz(i,j-1,k) + kzz(i,j-2,k)) - - LIEK_LKXY = LIEK_LKXY + 0.5D0*LIEK_BY/dy \ - *(3.0D0*kxy(i,j,k) - 4.0D0*kxy(i,j-1,k) + kxy(i,j-2,k)) - - LIEK_LKXZ = LIEK_LKXZ + 0.5D0*LIEK_BY/dy \ - *(3.0D0*kxz(i,j,k) - 4.0D0*kxz(i,j-1,k) + kxz(i,j-2,k)) - - LIEK_LKYZ = LIEK_LKYZ + 0.5D0*LIEK_BY/dy \ - *(3.0D0*kyz(i,j,k) - 4.0D0*kyz(i,j-1,k) + kyz(i,j-2,k)) + LIEK_LKXX = LIEK_LKXX + ADM_ADV_DY_M2(kxx,i,j,k) + LIEK_LKYY = LIEK_LKYY + ADM_ADV_DY_M2(kyy,i,j,k) + LIEK_LKZZ = LIEK_LKZZ + ADM_ADV_DY_M2(kzz,i,j,k) + LIEK_LKXY = LIEK_LKXY + ADM_ADV_DY_M2(kxy,i,j,k) + LIEK_LKXZ = LIEK_LKXZ + ADM_ADV_DY_M2(kxz,i,j,k) + LIEK_LKYZ = LIEK_LKYZ + ADM_ADV_DY_M2(kyz,i,j,k) end if @@ -228,96 +132,50 @@ if (admmacros_advectionz .eq. 0) then /* CENTER Z */ LIEK_LKXX = LIEK_LKXX + DZDK_DZDKXX*LIEK_BZ - LIEK_LKYY = LIEK_LKYY + DZDK_DZDKYY*LIEK_BZ - LIEK_LKZZ = LIEK_LKZZ + DZDK_DZDKZZ*LIEK_BZ - LIEK_LKXY = LIEK_LKXY + DZDK_DZDKXY*LIEK_BZ - LIEK_LKXZ = LIEK_LKXZ + DZDK_DZDKXZ*LIEK_BZ - LIEK_LKYZ = LIEK_LKYZ + DZDK_DZDKYZ*LIEK_BZ else if (admmacros_advectionz .eq. 1) then /* UPWIND1 Z */ - LIEK_LKXX = LIEK_LKXX + LIEK_BZ \ - *(kxx(i,j,k+1) - kxx(i,j,k))/dz - - LIEK_LKYY = LIEK_LKYY + LIEK_BZ \ - *(kyy(i,j,k+1) - kyy(i,j,k))/dz - - LIEK_LKZZ = LIEK_LKZZ + LIEK_BZ \ - *(kzz(i,j,k+1) - kzz(i,j,k))/dz - - LIEK_LKXY = LIEK_LKXY + LIEK_BZ \ - *(kxy(i,j,k+1) - kxy(i,j,k))/dz - - LIEK_LKXZ = LIEK_LKXZ + LIEK_BZ \ - *(kxz(i,j,k+1) - kxz(i,j,k))/dz - - LIEK_LKYZ = LIEK_LKYZ + LIEK_BZ \ - *(kyz(i,j,k+1) - kyz(i,j,k))/dz + LIEK_LKXX = LIEK_LKXX + ADM_ADV_DY_P1(kxx,i,j,k) + LIEK_LKYY = LIEK_LKYY + ADM_ADV_DY_P1(kyy,i,j,k) + LIEK_LKZZ = LIEK_LKZZ + ADM_ADV_DY_P1(kzz,i,j,k) + LIEK_LKXY = LIEK_LKXY + ADM_ADV_DY_P1(kxy,i,j,k) + LIEK_LKXZ = LIEK_LKXZ + ADM_ADV_DY_P1(kxz,i,j,k) + LIEK_LKYZ = LIEK_LKYZ + ADM_ADV_DY_P1(kyz,i,j,k) else if (admmacros_advectionz .eq. -1) then /* UPWIND1 Z */ - LIEK_LKXX = LIEK_LKXX + LIEK_BZ \ - *(kxx(i,j,k) - kxx(i,j,k-1))/dz - - LIEK_LKYY = LIEK_LKYY + LIEK_BZ \ - *(kyy(i,j,k) - kyy(i,j,k-1))/dz - - LIEK_LKZZ = LIEK_LKZZ + LIEK_BZ \ - *(kzz(i,j,k) - kzz(i,j,k-1))/dz - - LIEK_LKXY = LIEK_LKXY + LIEK_BZ \ - *(kxy(i,j,k) - kxy(i,j,k-1))/dz - - LIEK_LKXZ = LIEK_LKXZ + LIEK_BZ \ - *(kxz(i,j,k) - kxz(i,j,k-1))/dz + LIEK_LKXX = LIEK_LKXX + ADM_ADV_DY_M1(kxx,i,j,k) + LIEK_LKYY = LIEK_LKYY + ADM_ADV_DY_M1(kyy,i,j,k) + LIEK_LKZZ = LIEK_LKZZ + ADM_ADV_DY_M1(kzz,i,j,k) + LIEK_LKXY = LIEK_LKXY + ADM_ADV_DY_M1(kxy,i,j,k) + LIEK_LKXZ = LIEK_LKXZ + ADM_ADV_DY_M1(kxz,i,j,k) + LIEK_LKYZ = LIEK_LKYZ + ADM_ADV_DY_M1(kyz,i,j,k) else if (admmacros_advectionz .eq. 2) then /* UPWIND2 Z */ - LIEK_LKXX = LIEK_LKXX - 0.5D0*LIEK_BZ/dz \ - *(3.0D0*kxx(i,j,k) - 4.0D0*kxx(i,j,k+1) + kxx(i,j,k+2)) - - LIEK_LKYY = LIEK_LKYY - 0.5D0*LIEK_BZ/dz \ - *(3.0D0*kyy(i,j,k) - 4.0D0*kyy(i,j,k+1) + kyy(i,j,k+2)) - - LIEK_LKZZ = LIEK_LKZZ - 0.5D0*LIEK_BZ/dz \ - *(3.0D0*kzz(i,j,k) - 4.0D0*kzz(i,j,k+1) + kzz(i,j,k+2)) - - LIEK_LKXY = LIEK_LKXY - 0.5D0*LIEK_BZ/dz \ - *(3.0D0*kxy(i,j,k) - 4.0D0*kxy(i,j,k+1) + kxy(i,j,k+2)) - - LIEK_LKXZ = LIEK_LKXZ - 0.5D0*LIEK_BZ/dz \ - *(3.0D0*kxz(i,j,k) - 4.0D0*kxz(i,j,k+1) + kxz(i,j,k+2)) - - LIEK_LKYZ = LIEK_LKYZ - 0.5D0*LIEK_BZ/dz \ - *(3.0D0*kyz(i,j,k) - 4.0D0*kyz(i,j,k+1) + kyz(i,j,k+2)) + LIEK_LKXX = LIEK_LKXX + ADM_ADV_DY_P2(kxx,i,j,k) + LIEK_LKYY = LIEK_LKYY + ADM_ADV_DY_P2(kyy,i,j,k) + LIEK_LKZZ = LIEK_LKZZ + ADM_ADV_DY_P2(kzz,i,j,k) + LIEK_LKXY = LIEK_LKXY + ADM_ADV_DY_P2(kxy,i,j,k) + LIEK_LKXZ = LIEK_LKXZ + ADM_ADV_DY_P2(kxz,i,j,k) + LIEK_LKYZ = LIEK_LKYZ + ADM_ADV_DY_P2(kyz,i,j,k) else if (admmacros_advectionz .eq. -2) then /* UPWIND2 Z */ - LIEK_LKXX = LIEK_LKXX + 0.5D0*LIEK_BZ/dz \ - *(3.0D0*kxx(i,j,k) - 4.0D0*kxx(i,j,k-1) + kxx(i,j,k-2)) - - LIEK_LKYY = LIEK_LKYY + 0.5D0*LIEK_BZ/dz \ - *(3.0D0*kyy(i,j,k) - 4.0D0*kyy(i,j,k-1) + kyy(i,j,k-2)) + LIEK_LKXX = LIEK_LKXX + ADM_ADV_DY_M2(kxx,i,j,k) + LIEK_LKYY = LIEK_LKYY + ADM_ADV_DY_M2(kyy,i,j,k) + LIEK_LKZZ = LIEK_LKZZ + ADM_ADV_DY_M2(kzz,i,j,k) + LIEK_LKXY = LIEK_LKXY + ADM_ADV_DY_M2(kxy,i,j,k) + LIEK_LKXZ = LIEK_LKXZ + ADM_ADV_DY_M2(kxz,i,j,k) + LIEK_LKYZ = LIEK_LKYZ + ADM_ADV_DY_M2(kyz,i,j,k) - LIEK_LKZZ = LIEK_LKZZ + 0.5D0*LIEK_BZ/dz \ - *(3.0D0*kzz(i,j,k) - 4.0D0*kzz(i,j,k-1) + kzz(i,j,k-2)) - - LIEK_LKXY = LIEK_LKXY + 0.5D0*LIEK_BZ/dz \ - *(3.0D0*kxy(i,j,k) - 4.0D0*kxy(i,j,k-1) + kxy(i,j,k-2)) - - LIEK_LKXZ = LIEK_LKXZ + 0.5D0*LIEK_BZ/dz \ - *(3.0D0*kxz(i,j,k) - 4.0D0*kxz(i,j,k-1) + kxz(i,j,k-2)) - - LIEK_LKYZ = LIEK_LKYZ + 0.5D0*LIEK_BZ/dz \ - *(3.0D0*kyz(i,j,k) - 4.0D0*kyz(i,j,k-1) + kyz(i,j,k-2)) end if - /* Extra terms in the Lie derivative. */ LIEK_LKXX = LIEK_LKXX + 2.0D0*(DXDB_DXDBX*LIEK_KXX \ -- cgit v1.2.3