Use mask without using excision calls, this may need some more tweaking

depending on how excision labels points in the mask, but this will be
sorted out in the near future.

Converted to F77, then realised that this might be a bad idea since people
might be including F90 in the CalcTmunu stuff, so left the files labeled
as F90 for now.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/ADMConstraints/trunk@69 b7a48df3-cbbf-4440-997f-b4b717c9f7fc

1 files changed, 112 insertions, 128 deletions

diff --git a/src/ADMConstraints.F b/src/ADMConstraints.F
index 8eaba09..7715ce9 100644
--- a/src/ADMConstraints.F
+++ b/src/ADMConstraints.F
@@ -1,5 +1,5 @@
 /*@@
-  @file      ADMConstraints.F
+  @file      ADMConstraints.F77
   @date      August 98
   @desc 
              Calculate the ADM Constraints for output:
@@ -34,17 +34,13 @@
       DECLARE_CCTK_PARAMETERS
       DECLARE_CCTK_FUNCTIONS
 
-      integer :: i,j,k
-      integer :: nx,ny,nz
-
-#ifdef EXCISION_LEGOEXCISION 
-      CCTK_REAL, allocatable, dimension (:,:,:) :: dirx,diry,dirz,aux
-#endif
+      integer i,j,k
+      integer nx,ny,nz
 
 c     Stencil width used for calculating constraints
 c     (for outer boundary condition)
 
-      integer, dimension(3),parameter :: sw = 1
+      integer sw(3)
 
 c     Return code from Cactus sync routine and boundary conditions.
 
@@ -52,14 +48,14 @@ c     Return code from Cactus sync routine and boundary conditions.
 
 c     Various real variables.
 
-      CCTK_REAL :: dx,dy,dz
-      CCTK_REAL :: m_rho,m_sx,m_sy,m_sz
-      CCTK_REAL :: pi,ialp,ialp2
-      CCTK_REAL :: det,uxx,uyy,uzz,uxy,uxz,uyz
+      CCTK_REAL dx,dy,dz
+      CCTK_REAL m_rho,m_sx,m_sy,m_sz
+      CCTK_REAL pi,ialp,ialp2
+      CCTK_REAL det,uxx,uyy,uzz,uxy,uxz,uyz
 
 c     Temporaries for the Stress-Energy tensor.
 
-      CCTK_REAL :: Ttt,Ttx,Tty,Ttz,Txx,Txy,Txz,Tyy,Tyz,Tzz
+      CCTK_REAL Ttt,Ttx,Tty,Ttz,Txx,Txy,Txz,Tyy,Tyz,Tzz
 
 c     Matter declarations.
 
@@ -76,6 +72,10 @@ c     Macros from Standard Einstein.
 
 c --------------------------------------------------------------
 
+      sw(1) = 1
+      sw(2) = 1
+      sw(3) = 1
+
 c     Grid parameters.
 
       dx = cctk_delta_space(1)
@@ -88,11 +88,19 @@ c     Grid parameters.
 
 c     Fill with zeros.
 
-      ham = 0.0D0
+      do k=1,nz
+         do j=1,ny
+            do i=1,nx
+
+	      ham(i,j,k) = 0.0D0
+	
+      	      momx(i,j,k) = 0.0D0
+     	      momy(i,j,k) = 0.0D0
+              momz(i,j,k) = 0.0D0	
 
-      momx = 0.0D0
-      momy = 0.0D0
-      momz = 0.0D0
+            end do
+         end do
+      end do
 
 c     Calculate constraints.
 
@@ -101,114 +109,125 @@ c     Calculate constraints.
       do k=2,nz-1
          do j=2,ny-1
             do i=2,nx-1
+               
+               if ((use_mask .eq. 0) .or. (emask(i,j,k) .eq. 1)) then
 
-               ialp  = 1.0D0/alp(i,j,k)
-               ialp2 = ialp**2
+                  ialp  = 1.0D0/alp(i,j,k)
+                  ialp2 = ialp**2
+                  
+c                 Calculate the stress energy tensor at this point
+c                 ------------------------------------------------
 
-c              Calculate the stress energy tensor at this point
-c              ------------------------------------------------
-
-c              This may be needed for CalcTmunu
+c                 This may be needed for CalcTmunu
 
 #include "CactusEinstein/ADMMacros/src/macro/DETG_guts.h"
 
-               det = DETG_DETCG
+                  det = DETG_DETCG
 
 #include "CactusEinstein/ADMMacros/src/macro/UPPERMET_guts.h"
 
-               uxx = UPPERMET_UXX; uxy = UPPERMET_UXY; uxz = UPPERMET_UXZ
-               uyy = UPPERMET_UYY; uyz = UPPERMET_UYZ; uzz = UPPERMET_UZZ
-
-c              Initialize stress-energy tensor components.
+                  uxx = UPPERMET_UXX; uxy = UPPERMET_UXY; uxz = UPPERMET_UXZ
+                  uyy = UPPERMET_UYY; uyz = UPPERMET_UYZ; uzz = UPPERMET_UZZ
 
-               Ttt = 0.0D0
+c                 Initialize stress-energy tensor components.
 
-               Ttx = 0.0D0; Tty = 0.0D0; Ttz = 0.0D0
+                  Ttt = 0.0D0
+                  
+                  Ttx = 0.0D0; Tty = 0.0D0; Ttz = 0.0D0
+                  
+                  Txx = 0.0D0; Tyy = 0.0D0; Tzz = 0.0D0
+                  Txy = 0.0D0; Txz = 0.0D0; Tyz = 0.0D0
 
-               Txx = 0.0D0; Tyy = 0.0D0; Tzz = 0.0D0
-               Txy = 0.0D0; Txz = 0.0D0; Tyz = 0.0D0
-
-c              Include macro for stress energy tensor.
+c                 Include macro for stress energy tensor.
 
 #include "CalcTmunu.inc"
 
-c              Calculate the hamiltonian constraint
-c              ------------------------------------
+c                 Calculate the hamiltonian constraint
+c                 ------------------------------------
 
-c              Geometric piece.
+c                 Geometric piece.
 
 #include "CactusEinstein/ADMMacros/src/macro/HAMADM_guts.h"
 
-c              Add matter terms: - 16*pi*rho
+c                 Add matter terms: - 16*pi*rho
 c
-c              with rho defined as:
+c                 with rho defined as:
 c
-c              rho = n_a n_b T^{ab} = n^a n^b T_{ab}
-c                  = (T_00 - 2 beta^i T_{i0} + beta^i beta^j T_{ij})/alpha^2
+c                 rho = n_a n_b T^{ab} = n^a n^b T_{ab}
+c                     = (T_00 - 2 beta^i T_{i0} + beta^i beta^j T_{ij})/alpha^2
 
-               m_rho = ialp2*Ttt
+                  m_rho = ialp2*Ttt
 
-               if (shift_state .eq. 1) then
+                  if (shift_state .eq. 1) then
 
-                  m_rho = m_rho + ialp2
-     &               *(betax(i,j,k)**2*Txx
-     &               + betay(i,j,k)**2*Tyy
-     &               + betaz(i,j,k)**2*Tzz
-     &               +(betax(i,j,k)*betay(i,j,k)*Txy
-     &               + betax(i,j,k)*betaz(i,j,k)*Txz
-     &               + betay(i,j,k)*betaz(i,j,k)*Tyz)*2.0D0
-     &               -(betax(i,j,k)*Ttx
-     &               + betay(i,j,k)*Tty
-     &               + betaz(i,j,k)*Ttz)*2.0D0)
+                     m_rho = m_rho + ialp2
+     &                  *(betax(i,j,k)**2*Txx
+     &                  + betay(i,j,k)**2*Tyy
+     &                  + betaz(i,j,k)**2*Tzz
+     &                  +(betax(i,j,k)*betay(i,j,k)*Txy
+     &                  + betax(i,j,k)*betaz(i,j,k)*Txz
+     &                  + betay(i,j,k)*betaz(i,j,k)*Tyz)*2.0D0
+     &                  -(betax(i,j,k)*Ttx
+     &                  + betay(i,j,k)*Tty
+     &                  + betaz(i,j,k)*Ttz)*2.0D0)
 
-               end if
+                  end if
 
-               ham(i,j,k) = HAMADM_HAMADM - 16.0D0*pi*m_rho
+                  ham(i,j,k) = HAMADM_HAMADM - 16.0D0*pi*m_rho
 
 
-c              Calculate the Momentum constraints
-c              ----------------------------------
+c                 Calculate the Momentum constraints
+c                 ----------------------------------
 
-c              Geometric piece.
+c                 Geometric piece.
 
 #include "CactusEinstein/ADMMacros/src/macro/MOMXADM_guts.h"
 #include "CactusEinstein/ADMMacros/src/macro/MOMYADM_guts.h"
 #include "CactusEinstein/ADMMacros/src/macro/MOMZADM_guts.h"
 
-c              Add matter terms: - 8*pi*S_i
+c                 Add matter terms: - 8*pi*S_i
 c
-c              with S_i defined as:
+c                 with S_i defined as:
 c
-c              S_i = - g_{ia} n_b T^{ab} = - g_i^a n^b T_{ab} 
-c                  = - (T_{i0} - beta^j T_{ij})/alpha
-
-               m_sx = - ialp*Ttx
-               m_sy = - ialp*Tty
-               m_sz = - ialp*Ttz
-
-               if (shift_state .eq. 1) then
-
-                  m_sx = m_sx + ialp
-     &                 *(betax(i,j,k)*Txx
-     &                 + betay(i,j,k)*Txy
-     &                 + betaz(i,j,k)*Txz)
-
-                  m_sy = m_sy + ialp
-     &                 *(betax(i,j,k)*Txy
-     &                 + betay(i,j,k)*Tyy
-     &                 + betaz(i,j,k)*Tyz)
-
-                  m_sz = m_sz + ialp
-     &                 *(betax(i,j,k)*Txz
-     &                 + betay(i,j,k)*Tyz
-     &                 + betaz(i,j,k)*Tzz)
+c                 S_i = - g_{ia} n_b T^{ab} = - g_i^a n^b T_{ab} 
+c                     = - (T_{i0} - beta^j T_{ij})/alpha
+
+                  m_sx = - ialp*Ttx
+                  m_sy = - ialp*Tty
+                  m_sz = - ialp*Ttz
+
+                  if (shift_state .eq. 1) then
+
+                     m_sx = m_sx + ialp
+     &                    *(betax(i,j,k)*Txx
+     &                    + betay(i,j,k)*Txy
+     &                    + betaz(i,j,k)*Txz)
+
+                     m_sy = m_sy + ialp
+     &                    *(betax(i,j,k)*Txy
+     &                    + betay(i,j,k)*Tyy
+     &                    + betaz(i,j,k)*Tyz)
+                     
+                     m_sz = m_sz + ialp
+     &                    *(betax(i,j,k)*Txz
+     &                    + betay(i,j,k)*Tyz
+     &                    + betaz(i,j,k)*Tzz)
+                     
+                  end if
+                  
+                  momx(i,j,k) = MOMXADM_MOMXADM - 8.0D0*pi*m_sx
+                  momy(i,j,k) = MOMYADM_MOMYADM - 8.0D0*pi*m_sy
+                  momz(i,j,k) = MOMZADM_MOMZADM - 8.0D0*pi*m_sz
+
+               else
+
+                  ham(i,j,k) = excised_value
+                  momx(i,j,k) = excised_value
+                  momy(i,j,k) = excised_value
+                  momz(i,j,k) = excised_value
 
                end if
 
-               momx(i,j,k) = MOMXADM_MOMXADM - 8.0D0*pi*m_sx
-               momy(i,j,k) = MOMYADM_MOMYADM - 8.0D0*pi*m_sy
-               momz(i,j,k) = MOMZADM_MOMZADM - 8.0D0*pi*m_sz
-
             end do
          end do 
       end do
@@ -220,42 +239,6 @@ c                  = - (T_{i0} - beta^j T_{ij})/alpha
 #include "CactusEinstein/ADMMacros/src/macro/MOMYADM_undefine.h"
 #include "CactusEinstein/ADMMacros/src/macro/MOMZADM_undefine.h"
 
-c     LegoExcision (must be done before symmetries are applied).
-
-      if (excise==1) then
-
-#ifdef EXCISION_LEGOEXCISION 
-
-         allocate(dirx(nx,ny,nz),diry(nx,ny,nz),dirz(nx,ny,nz))
-         allocate(aux(nx,ny,nz))
-
-         aux = 0.0D0
-
-         call excision_findboundary(ierr,emask,nx,ny,nz)
-
-         call CCTK_SyncGroup(ierr,cctkGH,"einstein::mask")
-         call CartSymGN(ierr,cctkGH,"einstein::mask")
-
-         call excision_findnormals(ierr,emask,dirx,diry,dirz,nx,ny,nz)
-
-         call excision_extrapolate(ierr,ham ,aux,emask,
-     .      dirx,diry,dirz,nx,ny,nz,0.0D0)
-         call excision_extrapolate(ierr,momx,aux,emask,
-     .      dirx,diry,dirz,nx,ny,nz,0.0D0)
-         call excision_extrapolate(ierr,momy,aux,emask,
-     .      dirx,diry,dirz,nx,ny,nz,0.0D0)
-         call excision_extrapolate(ierr,momz,aux,emask,
-     .      dirx,diry,dirz,nx,ny,nz,0.0D0)
-
-         deallocate(dirx,diry,dirz)
-         deallocate(aux)
-
-#else
-         call CCTK_WARN(0,"You have not compiled with LegoExcision")
-#endif
-
-      end if
-
 c     Apply symmetry boundary conditions.
 
       call CartSymGN(ierr,cctkGH,"admconstraints::hamiltonian")
@@ -280,15 +263,16 @@ c     Cartoon.
       if (cartoon==1) then
 
 #ifdef BETATHORNS_CARTOON2D
-         call BndCartoon2DGN(ierr,cctkGH,TENSORTYPE_SCALAR,"admconstraints::hamiltonian")
-         call BndCartoon2DGN(ierr,cctkGH,TENSORTYPE_U,"admconstraints::momentum")
+         if (CCTK_IsThornActive("Cartoon2D")==1) then
+            call BndCartoon2DGN(ierr,cctkGH,TENSORTYPE_SCALAR,"admconstraints::hamiltonian")
+            call BndCartoon2DGN(ierr,cctkGH,TENSORTYPE_U,"admconstraints::momentum")
+         end if
 #else
          call CCTK_WARN(0,"You have not compiled with Cartoon2D")
 #endif
 
       end if
 
-c     End
-
-      end subroutine ADMConstraints
+      return
+      end