Look at admmacros_advection? variables for which way to difference beta.

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinBase/ADMMacros/trunk@39 b1d164ef-f17a-46e7-89d4-021c7118ef4e

6 files changed, 274 insertions, 10 deletions

diff --git a/src/macro/DXDB_declare.h b/src/macro/DXDB_declare.h
index e102004..54aaaf6 100644
--- a/src/macro/DXDB_declare.h
+++ b/src/macro/DXDB_declare.h
@@ -14,18 +14,36 @@
 #ifdef FCODE
 
 /* Input variables */
+#undef  DXDB_BX
+#define DXDB_BX betax(i,j,k)
+#undef  DXDB_BY
+#define DXDB_BY betay(i,j,k)
+#undef  DXDB_BZ
+#define DXDB_BZ betaz(i,j,k)
 #undef  DXDB_BX_IP 
 #define DXDB_BX_IP betax(i+1,j,k)
+#undef  DXDB_BX_IPP
+#define DXDB_BX_IPP betax(i+2,j,k)
 #undef  DXDB_BX_IM
 #define DXDB_BX_IM betax(i-1,j,k)
+#undef  DXDB_BX_IMM
+#define DXDB_BX_IMM betax(i-2,j,k)
 #undef  DXDB_BY_IP 
 #define DXDB_BY_IP betay(i+1,j,k)
+#undef  DXDB_BY_IPP 
+#define DXDB_BY_IPP betay(i+2,j,k)
 #undef  DXDB_BY_IM 
 #define DXDB_BY_IM betay(i-1,j,k)
+#undef  DXDB_BY_IMM
+#define DXDB_BY_IMM betay(i-2,j,k)
 #undef  DXDB_BZ_IP 
 #define DXDB_BZ_IP betaz(i+1,j,k)
+#undef  DXDB_BZ_IPP 
+#define DXDB_BZ_IPP betaz(i+2,j,k)
 #undef  DXDB_BZ_IM
 #define DXDB_BZ_IM betaz(i-1,j,k)
+#undef  DXDB_BZ_IMM
+#define DXDB_BZ_IMM betaz(i-2,j,k)
 
 /* Output variables */ 
 #undef  DXDB_DXDBX
@@ -55,6 +73,12 @@
 #ifdef CCODE
 
 /* Input variables */
+#undef  DXDB_BX
+#define DXDB_BX betax[ijk]
+#undef  DXDB_BY
+#define DXDB_BY betay[ijk]
+#undef  DXDB_BZ
+#define DXDB_BZ betaz[ijk]
 #undef  DXDB_BX_IP 
 #define DXDB_BX_IP betax[ di+ijk]
 #undef  DXDB_BX_IM
diff --git a/src/macro/DXDB_guts.h b/src/macro/DXDB_guts.h
index 46a7f9f..3b1531b 100644
--- a/src/macro/DXDB_guts.h
+++ b/src/macro/DXDB_guts.h
@@ -18,10 +18,78 @@
 #ifdef FCODE 
 
       DXDB_OO2DX = 1D0/(2D0*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)
+
+      if (admmacros_advectionx .eq. 0) then /* CENTER */
+
+        DXDB_DXDBX = DXDB_OO2DX*(DXDB_BX_IP - DXDB_BX_IM)
+
+      else if (admmacros_advectionx .eq. 1) then /* UPWIND1 */
+
+        DXDB_DXDBX = (DXDB_BX_IP - DXDB_BX)/DXDB_DX
+
+      else if (admmacros_advectionx .eq. -1) then /* UPWIND1 */
+
+        DXDB_DXDBX = (DXDB_BX - DXDB_BX_IM)/DXDB_DX
+
+      else if (admmacros_advectionx .eq. 2) then /* UPWIND2 */
+     
+        DXDB_DXDBX = - 0.5D0/DXDB_DX \
+          *(3.0D0*DXDB_BX - 4.0D0*DXDB_BX_IP + DXDB_BX_IPP)     
+
+      else if (admmacros_advectionx .eq. -2) then /* UPWIND2 */
+
+         DXDB_DXDBX =  0.5D0/DXDB_DX \
+          *(3.0D0*DXDB_BX - 4.0D0*DXDB_BX_IM + DXDB_BX_IMM)
+ 
+      end if
+
+      if (admmacros_advectiony .eq. 0) then /* CENTER */
+
+        DXDB_DXDBY = DXDB_OO2DX*(DXDB_BY_IP - DXDB_BY_IM)
+
+      else if (admmacros_advectiony .eq. 1) then /* UPWIND1 */
+
+        DXDB_DXDBY = (DXDB_BY_IP - DXDB_BY)/DXDB_DX
+
+      else if (admmacros_advectiony .eq. -1) then /* UPWIND1 */
+
+        DXDB_DXDBY = (DXDB_BY - DXDB_BY_IM)/DXDB_DX
+
+      else if (admmacros_advectiony .eq. 2) then /* UPWIND2 */
+     
+        DXDB_DXDBY = - 0.5D0/DXDB_DX \
+          *(3.0D0*DXDB_BY - 4.0D0*DXDB_BY_IP + DXDB_BY_IPP)     
+
+      else if (admmacros_advectiony .eq. -2) then /* UPWIND2 */
+
+         DXDB_DXDBY =  0.5D0/DXDB_DX \
+          *(3.0D0*DXDB_BY - 4.0D0*DXDB_BY_IM + DXDB_BY_IMM)
+
+      end if
+
+      if (admmacros_advectionz .eq. 0) then /* CENTER */
+
+        DXDB_DXDBZ = DXDB_OO2DX*(DXDB_BZ_IP - DXDB_BZ_IM)
+
+      else if (admmacros_advectionz .eq. 1) then /* UPWIND1 */
+
+        DXDB_DXDBZ = (DXDB_BZ_IP - DXDB_BZ)/DXDB_DX
+
+      else if (admmacros_advectionz .eq. -1) then /* UPWIND1 */
+
+        DXDB_DXDBZ = (DXDB_BZ - DXDB_BZ_IM)/DXDB_DX
+
+      else if (admmacros_advectionz .eq. 2) then /* UPWIND2 */
+     
+        DXDB_DXDBZ = - 0.5D0/DXDB_DX \
+          *(3.0D0*DXDB_BZ - 4.0D0*DXDB_BZ_IP + DXDB_BZ_IPP)     
+
+      else if (admmacros_advectionz .eq. -2) then /* UPWIND2 */
+
+         DXDB_DXDBY =  0.5D0/DXDB_DX \
+          *(3.0D0*DXDB_BZ - 4.0D0*DXDB_BZ_IM + DXDB_BZ_IMM)
+
+      end if
 
 #endif
 
diff --git a/src/macro/DYDB_declare.h b/src/macro/DYDB_declare.h
index 3bf2eac..3b73b7d 100644
--- a/src/macro/DYDB_declare.h
+++ b/src/macro/DYDB_declare.h
@@ -14,18 +14,36 @@
 #ifdef FCODE
 
 /* Input variables */
+#undef  DYDB_BX
+#define DYDB_BX betax(i,j,k)
+#undef  DYDB_BY
+#define DYDB_BY betay(i,j,k)
+#undef  DYDB_BZ
+#define DYDB_BZ betaz(i,j,k)
 #undef  DYDB_BX_JP 
 #define DYDB_BX_JP betax(i,j+1,k)
+#undef  DYDB_BX_JPP 
+#define DYDB_BX_JPP betax(i,j+2,k)
 #undef  DYDB_BX_JM
 #define DYDB_BX_JM betax(i,j-1,k)
+#undef  DYDB_BX_JMM
+#define DYDB_BX_JMM betax(i,j-2,k)
 #undef  DYDB_BY_JP 
 #define DYDB_BY_JP betay(i,j+1,k)
+#undef  DYDB_BY_JPP 
+#define DYDB_BY_JPP betay(i,j+2,k)
 #undef  DYDB_BY_JM 
 #define DYDB_BY_JM betay(i,j-1,k)
+#undef  DYDB_BY_JMM 
+#define DYDB_BY_JMM betay(i,j-2,k)
 #undef  DYDB_BZ_JP 
 #define DYDB_BZ_JP betaz(i,j+1,k)
+#undef  DYDB_BZ_JPP 
+#define DYDB_BZ_JPP betaz(i,j+2,k)
 #undef  DYDB_BZ_JM
 #define DYDB_BZ_JM betaz(i,j-1,k)
+#undef  DYDB_BZ_JMM
+#define DYDB_BZ_JMM betaz(i,j-2,k)
 
 /* Output variables */ 
 #undef  DYDB_DYDBX
diff --git a/src/macro/DYDB_guts.h b/src/macro/DYDB_guts.h
index 24013d3..015c8d4 100644
--- a/src/macro/DYDB_guts.h
+++ b/src/macro/DYDB_guts.h
@@ -19,9 +19,77 @@
 
       DYDB_OO2DY = 1D0/(2D0*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)
+      if (admmacros_advectionx .eq. 0) then /* CENTER */
+
+        DYDB_DYDBX = DYDB_OO2DY*(DYDB_BX_JP - DYDB_BX_JM)
+
+      else if (admmacros_advectionx .eq. 1) then /* UPWIND1 */
+
+        DYDB_DYDBX = (DYDB_BX_JP - DYDB_BX)/DYDB_DY
+
+      else if (admmacros_advectionx .eq. -1) then /* UPWIND1 */
+
+        DYDB_DYDBX = (DYDB_BX - DYDB_BX_JM)/DYDB_DY
+
+      else if (admmacros_advectionx .eq. 2) then /* UPWIND2 */
+     
+        DYDB_DYDBX = - 0.5D0/DYDB_DY \
+          *(3.0D0*DYDB_BX - 4.0D0*DYDB_BX_JP + DYDB_BX_JPP)     
+
+      else if (admmacros_advectionx .eq. -2) then /* UPWIND2 */
+
+         DYDB_DYDBX =  0.5D0/DYDB_DY \
+          *(3.0D0*DYDB_BX - 4.0D0*DYDB_BX_JM + DYDB_BX_JMM)
+ 
+      end if
+
+      if (admmacros_advectiony .eq. 0) then /* CENTER */
+
+        DYDB_DYDBY = DYDB_OO2DY*(DYDB_BY_JP - DYDB_BY_JM)
+
+      else if (admmacros_advectiony .eq. 1) then /* UPWIND1 */
+
+        DYDB_DYDBY = (DYDB_BY_JP - DYDB_BY)/DYDB_DY
+
+      else if (admmacros_advectiony .eq. -1) then /* UPWIND1 */
+
+        DYDB_DYDBY = (DYDB_BY - DYDB_BY_JM)/DYDB_DY
+
+      else if (admmacros_advectiony .eq. 2) then /* UPWIND2 */
+     
+        DYDB_DYDBY = - 0.5D0/DYDB_DY \
+          *(3.0D0*DYDB_BY - 4.0D0*DYDB_BY_JP + DYDB_BY_JPP)     
+
+      else if (admmacros_advectiony .eq. -2) then /* UPWIND2 */
+
+         DYDB_DYDBY =  0.5D0/DYDB_DY \
+          *(3.0D0*DYDB_BY - 4.0D0*DYDB_BY_JM + DYDB_BY_JMM)
+
+      end if
+
+      if (admmacros_advectionz .eq. 0) then /* CENTER */
+
+        DYDB_DYDBZ = DYDB_OO2DY*(DYDB_BZ_JP - DYDB_BZ_JM)
+
+      else if (admmacros_advectionz .eq. 1) then /* UPWIND1 */
+
+        DYDB_DYDBZ = (DYDB_BZ_JP - DYDB_BZ)/DYDB_DY
+
+      else if (admmacros_advectionz .eq. -1) then /* UPWIND1 */
+
+        DYDB_DYDBZ = (DYDB_BZ - DYDB_BZ_JM)/DYDB_DY
+
+      else if (admmacros_advectionz .eq. 2) then /* UPWIND2 */
+     
+        DYDB_DYDBZ = - 0.5D0/DYDB_DY \
+          *(3.0D0*DYDB_BZ - 4.0D0*DYDB_BZ_JP + DYDB_BZ_JPP)     
+
+      else if (admmacros_advectionz .eq. -2) then /* UPWIND2 */
+
+         DYDB_DYDBY =  0.5D0/DYDB_DY \
+          *(3.0D0*DYDB_BZ - 4.0D0*DYDB_BZ_JM + DYDB_BZ_JMM)
+
+      end if
 
 #endif
 
diff --git a/src/macro/DZDB_declare.h b/src/macro/DZDB_declare.h
index eceeec4..6e08b3a 100644
--- a/src/macro/DZDB_declare.h
+++ b/src/macro/DZDB_declare.h
@@ -14,18 +14,36 @@
 #ifdef FCODE
 
 /* Input variables */
+#undef  DZDB_BX
+#define DZDB_BX betax(i,j,k)
+#undef  DZDB_BY
+#define DZDB_BY betay(i,j,k)
+#undef  DZDB_BZ
+#define DZDB_BZ betaz(i,j,k)
 #undef  DZDB_BX_KP 
 #define DZDB_BX_KP betax(i,j,k+1)
+#undef  DZDB_BX_KPP 
+#define DZDB_BX_KPP betax(i,j,k+2)
 #undef  DZDB_BX_KM
 #define DZDB_BX_KM betax(i,j,k-1)
+#undef  DZDB_BX_KMM
+#define DZDB_BX_KMM betax(i,j,k-2)
 #undef  DZDB_BY_KP 
 #define DZDB_BY_KP betay(i,j,k+1)
+#undef  DZDB_BY_KPP 
+#define DZDB_BY_KPP betay(i,j,k+2)
 #undef  DZDB_BY_KM 
 #define DZDB_BY_KM betay(i,j,k-1)
+#undef  DZDB_BY_KMM 
+#define DZDB_BY_KMM betay(i,j,k-2)
 #undef  DZDB_BZ_KP 
 #define DZDB_BZ_KP betaz(i,j,k+1)
+#undef  DZDB_BZ_KPP 
+#define DZDB_BZ_KPP betaz(i,j,k+2)
 #undef  DZDB_BZ_KM
 #define DZDB_BZ_KM betaz(i,j,k-1)
+#undef  DZDB_BZ_KMM
+#define DZDB_BZ_KMM betaz(i,j,k-2)
 
 /* Output variables */ 
 #undef  DZDB_DZDBX
diff --git a/src/macro/DZDB_guts.h b/src/macro/DZDB_guts.h
index fbc1240..adcb1d1 100644
--- a/src/macro/DZDB_guts.h
+++ b/src/macro/DZDB_guts.h
@@ -19,9 +19,77 @@
 
       DZDB_OO2DZ = 1D0/(2D0*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)
+      if (admmacros_advectionx .eq. 0) then /* CENTER */
+
+        DZDB_DZDBX = DZDB_OO2DZ*(DZDB_BX_KP - DZDB_BX_KM)
+
+      else if (admmacros_advectionx .eq. 1) then /* UPWIND1 */
+
+        DZDB_DZDBX = (DZDB_BX_KP - DZDB_BX)/DZDB_DX
+
+      else if (admmacros_advectionx .eq. -1) then /* UPWIND1 */
+
+        DZDB_DZDBX = (DZDB_BX - DZDB_BX_KM)/DZDB_DX
+
+      else if (admmacros_advectionx .eq. 2) then /* UPWIND2 */
+     
+        DZDB_DZDBX = - 0.5D0/DZDB_DX \
+          *(3.0D0*DZDB_BX - 4.0D0*DZDB_BX_KP + DZDB_BX_KPP)     
+
+      else if (admmacros_advectionx .eq. -2) then /* UPWIND2 */
+
+         DZDB_DZDBX =  0.5D0/DZDB_DX \
+          *(3.0D0*DZDB_BX - 4.0D0*DZDB_BX_KM + DZDB_BX_KMM)
+ 
+      end if
+
+      if (admmacros_advectiony .eq. 0) then /* CENTER */
+
+        DZDB_DZDBY = DZDB_OO2DZ*(DZDB_BY_KP - DZDB_BY_KM)
+
+      else if (admmacros_advectiony .eq. 1) then /* UPWIND1 */
+
+        DZDB_DZDBY = (DZDB_BY_KP - DZDB_BY)/DZDB_DX
+
+      else if (admmacros_advectiony .eq. -1) then /* UPWIND1 */
+
+        DZDB_DZDBY = (DZDB_BY - DZDB_BY_KM)/DZDB_DX
+
+      else if (admmacros_advectiony .eq. 2) then /* UPWIND2 */
+     
+        DZDB_DZDBY = - 0.5D0/DZDB_DX \
+          *(3.0D0*DZDB_BY - 4.0D0*DZDB_BY_KP + DZDB_BY_KPP)     
+
+      else if (admmacros_advectiony .eq. -2) then /* UPWIND2 */
+
+         DZDB_DZDBY =  0.5D0/DZDB_DX \
+          *(3.0D0*DZDB_BY - 4.0D0*DZDB_BY_KM + DZDB_BY_KMM)
+
+      end if
+
+      if (admmacros_advectionz .eq. 0) then /* CENTER */
+
+        DZDB_DZDBZ = DZDB_OO2DZ*(DZDB_BZ_KP - DZDB_BZ_KM)
+
+      else if (admmacros_advectionz .eq. 1) then /* UPWIND1 */
+
+        DZDB_DZDBZ = (DZDB_BZ_KP - DZDB_BZ)/DZDB_DX
+
+      else if (admmacros_advectionz .eq. -1) then /* UPWIND1 */
+
+        DZDB_DZDBZ = (DZDB_BZ - DZDB_BZ_KM)/DZDB_DX
+
+      else if (admmacros_advectionz .eq. 2) then /* UPWIND2 */
+     
+        DZDB_DZDBZ = - 0.5D0/DZDB_DX \
+          *(3.0D0*DZDB_BZ - 4.0D0*DZDB_BZ_KP + DZDB_BZ_KPP)     
+
+      else if (admmacros_advectionz .eq. -2) then /* UPWIND2 */
+
+         DZDB_DZDBY =  0.5D0/DZDB_DX \
+          *(3.0D0*DZDB_BZ - 4.0D0*DZDB_BZ_KM + DZDB_BZ_KMM)
+
+      end if
 
 #endif