From 7a442cbe3519c63b16625669adf8762b8dc13f46 Mon Sep 17 00:00:00 2001 From: allen Date: Sat, 3 Apr 1999 23:53:33 +0000 Subject: 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 --- src/ADMConstraints.F | 188 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 188 insertions(+) create mode 100644 src/ADMConstraints.F (limited to 'src/ADMConstraints.F') 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 + -- cgit v1.2.3