aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/ADMConstraints.F240
-rw-r--r--src/InitSymBound.F25
-rw-r--r--src/ParamCheck.c6
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");
+ }
+
}
/*@@
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-15 01:51:57 +0000