diff options
author | allen <allen@b7a48df3-cbbf-4440-997f-b4b717c9f7fc> | 2002-05-27 18:26:50 +0000 |
---|---|---|
committer | allen <allen@b7a48df3-cbbf-4440-997f-b4b717c9f7fc> | 2002-05-27 18:26:50 +0000 |
commit | 812c1f169ce7678111dd0c21c878f38ced014b67 (patch) | |
tree | 720c52f8a10cbf6cfc592aec7547ecf4f0107c56 /src | |
parent | 0e9d7de52d12fc6f3f2d7179aec4187bae08e90a (diff) |
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
Diffstat (limited to 'src')
-rw-r--r-- | src/ADMConstraints.F | 240 | ||||
-rw-r--r-- | src/InitSymBound.F | 25 | ||||
-rw-r--r-- | src/ParamCheck.c | 6 |
3 files changed, 137 insertions, 134 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 diff --git a/src/InitSymBound.F b/src/InitSymBound.F index fdefad6..00cfefa 100644 --- a/src/InitSymBound.F +++ b/src/InitSymBound.F @@ -5,6 +5,7 @@ @desc Sets the symmetries for the Constraint grid functions @enddesc + @version $Header$ @@*/ #include "cctk.h" @@ -31,29 +32,41 @@ DECLARE_CCTK_ARGUMENTS - INTEGER, PARAMETER :: one = 1 - INTEGER, DIMENSION(3) :: sym - INTEGER :: ierr + INTEGER one + PARAMETER (one=1) + + INTEGER sym(3) + INTEGER ierr -c GROUP: constraints sym(1) = one sym(2) = one sym(3) = one call SetCartSymVN(ierr, cctkGH, sym, 'admconstraints::ham') + if (ierr .lt. 0) then + call CCTK_WARN(0,"Error setting symmetries for ADMConstraints::ham") + end if sym(1) = -one sym(2) = one sym(3) = one call SetCartSymVN(ierr, cctkGH, sym, 'admconstraints::momx') + if (ierr .lt. 0) then + call CCTK_WARN(0,"Error setting symmetries for ADMConstraints::momx") + end if sym(1) = one sym(2) = -one sym(3) = one call SetCartSymVN(ierr, cctkGH, sym, 'admconstraints::momy') + if (ierr .lt. 0) then + call CCTK_WARN(0,"Error setting symmetries for ADMConstraints::momy") + end if sym(1) = one sym(2) = one sym(3) = -one call SetCartSymVN(ierr, cctkGH, sym, 'admconstraints::momz') + if (ierr .lt. 0) then + call CCTK_WARN(0,"Error setting symmetries for ADMConstraints::momz") + end if return - - end subroutine ADMConstraint_InitSymbound + end diff --git a/src/ParamCheck.c b/src/ParamCheck.c index 45c15e9..0ca88e8 100644 --- a/src/ParamCheck.c +++ b/src/ParamCheck.c @@ -75,6 +75,12 @@ void ADMConstraints_ParamCheck(CCTK_ARGUMENTS) CCTK_PARAMWARN("ADMConstraints can currently only work with a physical metric or a static conformal metric with second derivatives"); } + CCTK_INFO("ADMConstraints contains ifdefs for Cartoon2D"); + if (cartoon && !CCTK_IsThornActive("Cartoon2D")) + { + CCTK_WARN(0," ... using Cartoon2D but it is not active"); + } + } /*@@ |