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

1 files changed, 72 insertions, 4 deletions

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