CactusBase is of course a totally inappropriate place to put the

constraints, so I'll try again ...



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

4 files changed, 199 insertions, 0 deletions

diff --git a/src/ADMConstraints.F b/src/ADMConstraints.F
new file mode 100644
index 0000000..d91cdf7
--- /dev/null
+++ b/src/ADMConstraints.F
@@ -0,0 +1,188 @@
+cc/*@@
+c  @file      ADMConstraints.F
+c  @date      Aug 98
+c  @desc 
+c     Calculate the ADM Constraints for output:
+c    
+c     Hamiltonian Constraint is:
+c
+c       R - K^i_j K^j_i + trK^2 - 16 Pi rho
+c
+c     Momentum Constraints are:
+c
+c       Del_j K_i^j - Del_i trK - 8 Pi j_i
+c
+c  @enddesc 
+c@@*/
+
+c --------------------------------------------------------------
+
+c CCTK Header files
+#include "cctk.h"
+#include "declare_parameters.h"
+#include "declare_arguments.h"
+
+c Using macro definitions from Einstein
+#include "../../packages/CactusEinstein/Einstein/src/Einstein.h"
+
+c --------------------------------------------------------------
+
+c/*@@
+c  @routine    ADMConstraints
+c  @date       Aug 98
+c  @desc 
+c    Calculate the ADM Constraints
+c  @enddesc 
+c  @calls     
+c  @calledby   
+c  @history 
+c
+c  @endhistory 
+c@@*/
+
+      subroutine ADMConstraints(CCTK_FARGUMENTS)
+      
+      USE ADM_Scalars
+
+      implicit none
+
+c     Declare variables in the dynamic argument list
+      DECLARE_CCTK_FARGUMENTS
+
+c     Declare parameters
+      DECLARE_PARAMETERS
+
+c     Perhaps this and others should go into cctk.h
+      INTEGER CCTK_Equals
+
+
+
+      integer :: i,j,k,i1,i2,j1,j2,k1,k2
+      REAL :: dx,dy,dz
+
+      REAL :: pi,det,uxx,uxy,uxz,uyy,uyz,uzz
+
+c     Temporaries for the Stress-Energy tensor
+      REAL :: Ttt,Ttx,Tty,Ttz,Txx,Txy,Txz,Tyy,Tyz,Tzz
+
+
+#include "CalcTmunu_temps.inc"
+
+c     Macros from Standard Einstein
+#include "../../packages/CactusEinstein/Einstein/src/macro/HAMADM_declare.h"
+#include "../../packages/CactusEinstein/Einstein/src/macro/MOMXADM_declare.h"
+#include "../../packages/CactusEinstein/Einstein/src/macro/MOMYADM_declare.h"
+#include "../../packages/CactusEinstein/Einstein/src/macro/MOMZADM_declare.h"
+#include "../../packages/CactusEinstein/Einstein/src/macro/DETG_declare.h"
+#include "../../packages/CactusEinstein/Einstein/src/macro/UPPERMET_declare.h"
+
+      dx = delta_space(1)
+      dy = delta_space(2)
+      dz = delta_space(3)
+
+      i1 = 2
+      i2 = sh(1)-1
+      j1 = 2
+      j2 = sh(2)-1
+      k1 = 2
+      k2 = sh(2)-1
+  
+      do i=i1,i2
+        do j=j1,j2
+          do k=k1,k2
+            ham(i,j,k)  = 0.0D0
+            momx(i,j,k) = 0.0D0
+            momy(i,j,k) = 0.0D0
+            momz(i,j,k) = 0.0D0
+          end do
+        end do
+      end do
+
+      Pi =  ACOS(-1D0)
+
+      do k = k1, k2
+         do j = j1, j2
+            do i = k1, k2
+               
+c              Calculate the stress energy tensor at this point
+c              ------------------------------------------------
+
+               Ttt = 0.0D0; Ttx = 0.0D0; Tty = 0.0D0; Ttz = 0.0D0
+               Txx = 0.0D0; Txy = 0.0D0; Txz = 0.0D0; Tyy = 0.0D0
+               Tyz = 0.0D0; Tzz = 0.0D0
+               
+c              This may be needed for CalcTmunu etc
+#include "../../packages/CactusEinstein/Einstein/src/macro/DETG_guts.h"
+               det = DETG_DETCG
+#include "../../packages/CactusEinstein/Einstein/src/macro/UPPERMET_guts.h"
+               uxx = UPPERMET_UXX; uxy = UPPERMET_UXY; uxz = UPPERMET_UXZ
+               uyy = UPPERMET_UYY; uyz = UPPERMET_UYZ; uzz = UPPERMET_UZZ
+#include "CalcTmunu.inc"
+              
+
+c              Calculate the hamiltonian constraint
+c              ------------------------------------
+
+#include "../../packages/CactusEinstein/Einstein/src/macro/HAMADM_guts.h"
+
+               ham(i,j,k) = HAMADM_HAMADM - 16D0*Pi*Ttt/alp(i,j,k)**2
+
+#ifdef THORN_MAHC
+               if ((some_hydro .eq. MAHCHYDRO).or.
+     .              (some_hydro .eq. MAHCHYDRODUST)) then
+                  ham(i,j,k) = ADM_ham(i,j,k) - 16.0d0 * pi *
+     .                 ( (rho(i,j,k) + rho(i,j,k)*eps(i,j,k) + press(i,j,k)) *
+     .                 w_lorentz(i,j,k)**2 - press(i,j,k))
+               endif
+#endif
+
+
+c              Calculate the Momentum constraints
+c              ----------------------------------
+
+#include "../../packages/CactusEinstein/Einstein/src/macro/MOMXADM_guts.h"
+#include "../../packages/CactusEinstein/Einstein/src/macro/MOMYADM_guts.h"
+#include "../../packages/CactusEinstein/Einstein/src/macro/MOMZADM_guts.h"
+
+               momx(i,j,k) = MOMXADM_MOMXADM - 8D0*Pi*Ttx/alp(i,j,k)**2
+               momy(i,j,k) = MOMYADM_MOMYADM - 8D0*Pi*Tty/alp(i,j,k)**2
+               momz(i,j,k) = MOMZADM_MOMZADM - 8D0*Pi*Ttz/alp(i,j,k)**2
+
+
+#ifdef THORN_MAHC
+               if ((some_hydro .eq. MAHCHYDRO).or.
+     .              (some_hydro .eq. MAHCHYDRODUST)) then
+                  momx(i,j,k) = ADM_momx(i,j,k) -
+     .                 8.0d0 * pi * sx(i,j,k) / sqrt(det)
+                  momy(i,j,k) = ADM_momy(i,j,k) -
+     .                 8.0d0 * pi * sy(i,j,k) / sqrt(det)
+                  momz(i,j,k) = ADM_momz(i,j,k) -
+     .                 8.0d0 * pi * sz(i,j,k) / sqrt(det)
+               endif
+#endif
+               
+            enddo
+         enddo
+      enddo
+
+#include "../../packages/CactusEinstein/Einstein/src/macro/DETG_undefine.h"
+#include "../../packages/CactusEinstein/Einstein/src/macro/HAMADM_undefine.h"
+#include "../../packages/CactusEinstein/Einstein/src/macro/MOMXADM_undefine.h"
+#include "../../packages/CactusEinstein/Einstein/src/macro/MOMYADM_undefine.h"
+#include "../../packages/CactusEinstein/Einstein/src/macro/MOMZADM_undefine.h"
+#include "../../packages/CactusEinstein/Einstein/src/macro/UPPERMET_undefine.h"
+
+
+c     Synchronize and apply boundaries conditions
+
+c      call onebound(ARENA, ADM_ham,  1d0, 1d0, 1d0)
+c      call onebound(ARENA, ADM_momx,-1d0, 1d0, 1d0)
+c      call onebound(ARENA, ADM_momy, 1d0,-1d0, 1d0)
+c      call onebound(ARENA, ADM_momz, 1d0,-1d0, 1d0)
+
+      if (constraint_communication.eq.1) then
+         call CCTK_SyncGroup(GH,"adm::ADM_constraints")
+      end if
+
+      end subroutine ADMConstraints
+
diff --git a/src/CalcTmunu.inc b/src/CalcTmunu.inc
new file mode 100644
index 0000000..9b32a17
--- /dev/null
+++ b/src/CalcTmunu.inc
@@ -0,0 +1 @@
+/*Stand in for .inc files Ed Evans*/
diff --git a/src/CalcTmunu_temps.inc b/src/CalcTmunu_temps.inc
new file mode 100644
index 0000000..9b32a17
--- /dev/null
+++ b/src/CalcTmunu_temps.inc
@@ -0,0 +1 @@
+/*Stand in for .inc files Ed Evans*/
diff --git a/src/make.code.defn b/src/make.code.defn
new file mode 100644
index 0000000..b5f7acc
--- /dev/null
+++ b/src/make.code.defn
@@ -0,0 +1,9 @@
+# Main make.code.defn file for thorn ADMConstraints
+# $Header$
+
+# Source files in this directory
+SRCS =  ADMConstraints.F
+
+# Subdirectories containing source files
+SUBDIRS = 
+