From ba7e3f62e276db1cac91bfcc74add1306501237a Mon Sep 17 00:00:00 2001 From: miguel Date: Fri, 16 Feb 2001 11:33:36 +0000 Subject: Cleaning up and adding comments to make it more readable. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/ADMConstraints/trunk@43 b7a48df3-cbbf-4440-997f-b4b717c9f7fc --- src/ADMConstraints.F | 106 ++++++++++++++++++++++++++++----------------------- 1 file changed, 59 insertions(+), 47 deletions(-) (limited to 'src') diff --git a/src/ADMConstraints.F b/src/ADMConstraints.F index 1a999dd..7c61d21 100644 --- a/src/ADMConstraints.F +++ b/src/ADMConstraints.F @@ -47,7 +47,8 @@ c@@*/ DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS - integer :: i,j,k,i1,i2,j1,j2,k1,k2 + integer :: i,j,k + integer :: nx,ny,nz c Stencil width used for calculating constraints (for outer boundary condition) integer, dimension(3),parameter :: sw = 1 @@ -71,33 +72,31 @@ c Macros from Standard Einstein #include "CactusEinstein/Einstein/src/macro/DETG_declare.h" #include "CactusEinstein/Einstein/src/macro/UPPERMET_declare.h" +c Grid parameters. + dx = cctk_delta_space(1) dy = cctk_delta_space(2) dz = cctk_delta_space(3) - i1 = 2 - i2 = cctk_lsh(1)-1 - j1 = 2 - j2 = cctk_lsh(2)-1 - k1 = 2 - k2 = cctk_lsh(3)-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 + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + +c Fill with zeros. + + ham = 0.0D0 + + momx = 0.0D0 + momy = 0.0D0 + momz = 0.0D0 + +c Calculate constraints. - Pi = ACOS(-1D0) + pi = acos(-1.0D0) - do k = k1,k2 - do j = j1,j2 - do i = i1,i2 + do k=2,nz-1 + do j=2,ny-1 + do i=2,nx-1 ialp = 1.0D0/alp(i,j,k) ialp2 = ialp**2 @@ -105,20 +104,25 @@ c Macros from Standard Einstein c Calculate the stress energy tensor at this point c ------------------------------------------------ +c Fill with zeroes. + 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 -c -------------------------------- + + Txx = 0.0D0; Txy = 0.0D0; Txz = 0.0D0 + Tyy = 0.0D0; Tyz = 0.0D0; Tzz = 0.0D0 + +c Inverse metric may be needed for CalcTmunu. #include "CactusEinstein/Einstein/src/macro/DETG_guts.h" #include "CactusEinstein/Einstein/src/macro/UPPERMET_guts.h" det = DETG_DETCG + uxx = UPPERMET_UXX; uxy = UPPERMET_UXY; uxz = UPPERMET_UXZ uyy = UPPERMET_UYY; uyz = UPPERMET_UYZ; uzz = UPPERMET_UZZ +c Call macro for stress energy tensor. + #include "CalcTmunu.inc" c Calculate the hamiltonian constraint @@ -154,13 +158,15 @@ c = (T_00 - 2 beta^i T_{i0} + beta^i beta^j T_{ij})/alpha^2 ham(i,j,k) = HAMADM_HAMADM - 16.0D0*pi*m_rho +c From thorn MAHC (this overwrites all we just did). + #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)) * + 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 + end if #endif c Calculate the Momentum constraints @@ -206,37 +212,43 @@ c = - (T_{i0} - beta^j T_{ij})/alpha momy(i,j,k) = MOMYADM_MOMYADM - 8.0D0*pi*m_sy momz(i,j,k) = MOMZADM_MOMZADM - 8.0D0*pi*m_sz +c From thorn MAHC (this overwrites all we just did). + #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 + 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) + end if #endif - - enddo - enddo - enddo + + end do + end do + end do #include "CactusEinstein/Einstein/src/macro/DETG_undefine.h" +#include "CactusEinstein/Einstein/src/macro/UPPERMET_undefine.h" #include "CactusEinstein/Einstein/src/macro/HAMADM_undefine.h" #include "CactusEinstein/Einstein/src/macro/MOMXADM_undefine.h" #include "CactusEinstein/Einstein/src/macro/MOMYADM_undefine.h" #include "CactusEinstein/Einstein/src/macro/MOMZADM_undefine.h" -#include "CactusEinstein/Einstein/src/macro/UPPERMET_undefine.h" -c Synchronize and apply flat boundary conditions +c Apply symmetry boundary conditions. call CartSymGN(ierr,cctkGH,"admconstraints::admconstraints") +c Apply flat boundary conditions at outer boundaries. + if (CCTK_Equals(bound,"flat") == 1) then call BndFlatGN(ierr,cctkGH,sw,"admconstraints::admconstraints") end if +c Synchronize. + if (constraint_communication.eq.1) then call CCTK_SyncGroup(cctkGH,"admconstraints::admconstraints") end if -- cgit v1.2.3