diff options
Diffstat (limited to 'src/DistortedBHIVP.F')
| -rw-r--r-- | src/DistortedBHIVP.F | 665 |
1 files changed, 665 insertions, 0 deletions
diff --git a/src/DistortedBHIVP.F b/src/DistortedBHIVP.F new file mode 100644 index 0000000..9682cf4 --- /dev/null +++ b/src/DistortedBHIVP.F @@ -0,0 +1,665 @@ +c /*@@ +c @file DistortedBHIVP.F +c @date +c @author +c @desc +c +c @enddesc +c@@ */ + +#include "cctk.h" +#include "cctk_parameters.h" +#include "cctk_arguments.h" + +c /*@@ +c @routine DistortedBHIVP +c @date +c @author +c @desc +c +c @enddesc +c @calls +c @calledby +c @history +c +c @endhistory +c@@ */ + +c Need include file from Einstein +#include "CactusEinstein/Einstein/src/Einstein.h" + + subroutine DistortedBHIVP(CCTK_FARGUMENTS) + + implicit none + + DECLARE_CCTK_FARGUMENTS + DECLARE_CCTK_PARAMETERS + +c Perhaps this and others should go into cctk.h + integer CCTK_Equals + + real*8 :: deta,dq,dphi + real*8, allocatable :: ac(:,:,:),ae(:,:,:),aw(:,:,:),an(:,:,:), + $ as(:,:,:),aq(:,:,:),ab(:,:,:),rhs(:,:,:), + $ qf(:,:,:),qfetaeta(:,:,:),qfqq(:,:,:),qfphi(:,:,:), + $ qfphiphi(:,:,:), + $ psisph(:,:,:),psiprim(:,:,:),detapsisph(:,:,:), + $ dqpsisph(:,:,:),dphipsisph(:,:,:),detaetapsisph(:,:,:), + $ detaqpsisph(:,:,:),detaphipsisph(:,:,:),dqqpsisph(:,:,:), + $ dqphipsisph(:,:,:),dphiphipsisph(:,:,:), + $ psi3d(:,:,:),detapsi3d(:,:,:),dqpsi3d(:,:,:), + $ dphipsi3d(:,:,:),detaetapsi3d(:,:,:),detaqpsi3d(:,:,:), + $ detaphipsi3d(:,:,:),dqqpsi3d(:,:,:),dqphipsi3d(:,:,:), + $ dphiphipsi3d(:,:,:) + real*8, allocatable :: etagrd(:),qgrd(:),phigrd(:) + real*8, allocatable :: eta(:,:,:),abseta(:,:,:),sign_eta(:,:,:), + $ q(:,:,:),phi(:,:,:) + real*8 o1,o2,o3,o4,o5,o6,o7,o8,o9,o10,o11,o12,o13,o14,o15,o16,o17, + $ o18,o19,o20,o21,o22,o23,o24,o25,o26,o27,o28,o29,o30,o31,o32, + $ o33,o34,o35,o36,o37,o38,o39,o40,o41,o42,o43,o44,o45,o46,o47, + $ o48,o49,o50,o51,o52,o53,o54,o55,o56,o57,o58,o59,o60,o61,o62, + $ o63,o64,o65,o66,o67 + real*8 rmax,adm + real*8,parameter :: dbh_eps = 1.0d-9 + real*8 r,pi + integer :: nx,ny,nz + integer i,j,k,ier,nquads,nocts,order + integer npoints,handle,ierror + + conformal_state = CONFORMAL_METRIC + + pi = 4.0d0*atan(1.0d0) + +c Set up the grid spacings + nx = cctk_lsh(1) + ny = cctk_lsh(2) + nz = cctk_lsh(3) + +c Distorted Schwarzchild BH parameters + + print *,"Brill wave + Distorted BH solve" + write(*,123)amp,eta0,c,sigma,n + print*,'etamax=',etamax + 123 format(1x, 'Pars: amp',f8.5,' eta0',f8.5,' c',f8.5,' sigma',f8.5,' n',i3) + +c Sovle on this sized cartesian grid +c 3D grid size NE x NT x NP +c Add 2 zones for eta coordinate and 4 for theta +c and phi coordenate. + ne = ne + 2 + nq = nq + 2 + np = np + 2 +c + allocate(ac(ne,nq,np),ae(ne,nq,np),aw(ne,nq,np),an(ne,nq,np), + $ as(ne,nq,np),aq(ne,nq,np),ab(ne,nq,np),rhs(ne,nq,np), + $ qf(ne,nq,np),qfetaeta(ne,nq,np),qfqq(ne,nq,np), + $ qfphi(ne,nq,np),qfphiphi(ne,nq,np), + $ psisph(ne,nq,np),psiprim(ne,nq,np),detapsisph(ne,nq,np), + $ dqpsisph(ne,nq,np),dphipsisph(ne,nq,np), + $ detaetapsisph(ne,nq,np),detaqpsisph(ne,nq,np), + $ detaphipsisph(ne,nq,np),dqqpsisph(ne,nq,np), + $ dqphipsisph(ne,nq,np),dphiphipsisph(ne,nq,np), + $ psi3d(nx,ny,nz),detapsi3d(nx,ny,nz),dqpsi3d(nx,ny,nz), + $ dphipsi3d(nx,ny,nz),detaetapsi3d(nx,ny,nz), + $ detaqpsi3d(nx,ny,nz),detaphipsi3d(nx,ny,nz), + $ dqqpsi3d(nx,ny,nz),dqphipsi3d(nx,ny,nz), + $ dphiphipsi3d(nx,ny,nz)) + allocate(etagrd(ne),qgrd(nq),phigrd(np)) + allocate(eta(nx,ny,nz),abseta(nx,ny,nz),sign_eta(nx,ny,nz), + $ q(nx,ny,nz),phi(nx,ny,nz)) +c +c Initialize some array +c + nocts = 4 + nquads = 2 + dphi = nocts*0.5*pi/(np-2) + dq = nquads*0.5*pi/(nq-2) + deta = etamax/(ne-3) + + do k = 1,np + phigrd(k) = (k-1.5)*dphi + enddo + do j = 1,nq + qgrd(j) = (j-1.5)*dq + enddo + do i=1,ne + etagrd(i) = (i-2)*deta + enddo +c +c Initialize q-function and its derivatives: should be generalized +c + do k = 1,np + do j = 1,nq + do i = 1,ne +#include "CactusEinstein/DistortedBHIVP/src/qfunc.x" + enddo + enddo + enddo + +c +c Initialize psi to the Schwarzschild solution: +c + psiprim = 0. + + do k = 1,np + do j = 1,nq + do i = 1,ne + psisph(i,j,k) = 2.*cosh(0.5*etagrd(i)) + enddo + enddo + enddo +c +c Initialize stencil coefficients: +c + ac = 0. + ae = 0. + aw = 0. + an = 0. + as = 0. + aq = 0. + ab = 0. + rhs = 0. +c + + do k = 2,np-1 + do j = 2,nq-1 + do i = 2,ne-1 + ac(i,j,k) = -2./deta**2-2./dq**2-2.*exp(2.* + & qf(i,j,k))/(dphi**2* + & sin(qgrd(j))**2)+0.25* + & (qfetaeta(i,j,k)+qfqq(i,j,k)+2.* + & exp(2.*qf(i,j,k))*qfphiphi(i,j,k)/ + & sin(qgrd(j))**2+3.*exp(2.* + & qf(i,j,k))*qfphi(i,j,k)**2/ + & sin(qgrd(j))**2-1.) + ae(i,j,k) = 1./deta**2 + aw(i,j,k) = 1./deta**2 + an(i,j,k) = 1./dq**2+0.5/(dq*tan(qgrd(j))) + + as(i,j,k) = 1./dq**2-0.5/(dq*tan(qgrd(j))) + + aq(i,j,k) = exp(2.*qf(i,j,k))/(dphi**2* + & sin(qgrd(j))**2)+exp(2.* + & qf(i,j,k))*qfphi(i,j,k)/(dphi* + & sin(qgrd(j))**2) + + ab(i,j,k) = exp(2.*qf(i,j,k))/(dphi**2* + & sin(qgrd(j))**2)-exp(2.* + & qf(i,j,k))*qfphi(i,j,k)/(dphi* + & sin(qgrd(j))**2) + + rhs(i,j,k) = -0.25*(qfetaeta(i,j,k)+ + & qfqq(i,j,k)+2.*exp(2.*qf(i,j,k))* + & qfphiphi(i,j,k)/sin(qgrd(j))**2+3.* + & exp(2.*qf(i,j,k))*qfphi(i,j,k)**2/ + & sin(qgrd(j))**2)*psisph(i,j,k) + enddo + enddo + enddo +c +c Apply boundary conditions to the faces of the cube: +c +c i=2: + do k = 3,np-2 + do j = 3,nq-2 + ae(2,j,k) = ae(2,j,k) + aw(2,j,k) + aw(2,j,k) = 0. +c +c i=ne-1: + ac(ne-1,j,k) = ac(ne-1,j,k)+4.*ae(ne-1,j,k)/(3.+deta) + aw(ne-1,j,k) = aw(ne-1,j,k) - ae(ne-1,j,k)/(3.+deta) + ae(ne-1,j,k) = 0. + enddo + enddo +c +c j=2: + do k = 3,np-2 + do i = 3,ne-2 + ac(i,2,k) = ac(i,2,k) + as(i,2,k) + as(i,2,k) = 0. +c +c j=nq-1: + ac(i,nq-1,k) = ac(i,nq-1,k) + an(i,nq-1,k) + an(i,nq-1,k) = 0. + enddo + enddo +c +c k=2: + do j = 3,nq-2 + do i = 3,ne-2 + ac(i,j,2) = ac(i,j,2) + ab(i,j,2) + ab(i,j,2) = 0. +c +c k=np-1: + ac(i,j,np-1) = ac(i,j,np-1) + aq(i,j,np-1) + aq(i,j,np-1) = 0. + enddo + enddo +c +c Apply boundary conditions to the edges of the cube: +c +c i=2, j=2: + do k = 3,np-2 + ae(2,2,k) = ae(2,2,k) + aw(2,2,k) + ac(2,2,k) = ac(2,2,k) + as(2,2,k) + aw(2,2,k) = 0. + as(2,2,k) = 0. +c +c i=ne-1, j=2: + aw(ne-1,2,k) = aw(ne-1,2,k) - ae(ne-1,2,k)/(3.+deta) + ac(ne-1,2,k) = ac(ne-1,2,k) + as(ne-1,2,k) + + & 4.*ae(ne-1,2,k)/(3.+deta) + ae(ne-1,2,k) = 0. + as(ne-1,2,k) = 0. +c +c i=2, j=nq-1: + ae(2,nq-1,k) = ae(2,nq-1,k) + aw(2,nq-1,k) + ac(2,nq-1,k) = ac(2,nq-1,k) + an(2,nq-1,k) + aw(2,nq-1,k) = 0. + an(2,nq-1,k) = 0. +c +c i=ne-1, j=nq-1: + aw(ne-1,nq-1,k) = aw(ne-1,nq-1,k) - ae(ne-1,nq-1,k)/ + & (3.+deta) + ac(ne-1,nq-1,k) = ac(ne-1,nq-1,k) + an(ne-1,nq-1,k) + + & 4.*ae(ne-1,nq-1,k)/(3.+deta) + ae(ne-1,nq-1,k) = 0. + an(ne-1,nq-1,k) = 0. + enddo +c +c i=2, k=2: + do j = 3,nq-2 + ae(2,j,2) = ae(2,j,2) + aw(2,j,2) + ac(2,j,2) = ac(2,j,2) + ab(2,j,2) + aw(2,j,2) = 0. + ab(2,j,2) = 0. +c +c i=ne-1, k=2: + aw(ne-1,j,2) = aw(ne-1,j,2) - ae(ne-1,j,2)/(3.+deta) + ac(ne-1,j,2) = ac(ne-1,j,2) + ab(ne-1,j,2) + + & 4.*ae(ne-1,j,2)/(3.+deta) + ae(ne-1,j,2) = 0. + ab(ne-1,j,2) = 0. +c +c i=2, k=np-1: + ae(2,j,np-1) = ae(2,j,np-1) + aw(2,j,np-1) + ac(2,j,np-1) = ac(2,j,np-1) + aq(2,j,np-1) + aw(2,j,np-1) = 0. + aq(2,j,np-1) = 0. +c +c i=ne-1, k=np-1: + aw(ne-1,j,np-1) = aw(ne-1,j,np-1) - ae(ne-1,j,np-1)/ + & (3.+deta) + ac(ne-1,j,np-1) = ac(ne-1,j,np-1) + aq(ne-1,j,np-1) + + & 4.*ae(ne-1,j,np-1)/(3.+deta) + ae(ne-1,j,np-1) = 0. + aq(ne-1,j,np-1) = 0. + enddo +c +c j=2, k=2: + do i = 3,ne-2 + ac(i,2,2) = ac(i,2,2) + as(i,2,2) + ab(i,2,2) + as(i,2,2) = 0. + ab(i,2,2) = 0. +c +c j=nq-1, k=2: + ac(i,nq-1,2) = ac(i,nq-1,2) + an(i,nq-1,2) + + & ab(i,nq-1,2) + an(i,nq-1,2) = 0. + ab(i,nq-1,2) = 0. +c +c j=2, k=np-1: + ac(i,2,np-1) = ac(i,2,np-1) + as(i,2,np-1) + + & aq(i,2,np-1) + as(i,2,np-1) = 0. + aq(i,2,np-1) = 0. +c +c j=nq-1, k=np-1: + ac(i,nq-1,np-1) = ac(i,nq-1,np-1) + an(i,nq-1,np-1) + + & aq(i,nq-1,np-1) + an(i,nq-1,np-1) = 0. + aq(i,nq-1,np-1) = 0. + enddo +c +c Apply boundary conditions to the corners of the cube: +c +c i=2, j=2, k=2: + ae(2,2,2) = ae(2,2,2) + aw(2,2,2) + ac(2,2,2) = ac(2,2,2) + as(2,2,2) + ab(2,2,2) + aw(2,2,2) = 0. + as(2,2,2) = 0. + ab(2,2,2) = 0. +c +c i=ne-1, j=2, k=2: + aw(ne-1,2,2) = aw(ne-1,2,2) - ae(ne-1,2,2)/(3.+deta) + ac(ne-1,2,2) = ac(ne-1,2,2) + as(ne-1,2,2) + ab(ne-1,2,2) + + & 4.*ae(ne-1,2,2)/(3.+deta) + ae(ne-1,2,2) = 0. + as(ne-1,2,2) = 0. + ab(ne-1,2,2) = 0. +c +c i=2, j=nq-1, k=2: + ae(2,nq-1,2) = ae(2,nq-1,2) + aw(2,nq-1,2) + ac(2,nq-1,2) = ac(2,nq-1,2) + an(2,nq-1,2) + ab(2,nq-1,2) + aw(2,nq-1,2) = 0. + an(2,nq-1,2) = 0. + ab(2,nq-1,2) = 0. +c +c i=2, j=2, k=np-1: + ae(2,2,np-1) = ae(2,2,np-1) + aw(2,2,np-1) + ac(2,2,np-1) = ac(2,2,np-1) + as(2,2,np-1) + aq(2,2,np-1) + aw(2,2,np-1) = 0. + as(2,2,np-1) = 0. + aq(2,2,np-1) = 0. +c +c i=ne-1, j=nq-1, k=2: + aw(ne-1,nq-1,2) = aw(ne-1,nq-1,2) - ae(ne-1,nq-1,2)/(3.+deta) + ac(ne-1,nq-1,2) = ac(ne-1,nq-1,2) + an(ne-1,nq-1,2) + ab(ne-1,nq-1,2) + + & 4.*ae(ne-1,nq-1,2)/(3.+deta) + ae(ne-1,nq-1,2) = 0. + an(ne-1,nq-1,2) = 0. + ab(ne-1,nq-1,2) = 0. +c +c i=ne-1, j=2, k=np-1: + aw(ne-1,2,np-1) = aw(ne-1,2,np-1) - ae(ne-1,2,np-1)/(3.+deta) + ac(ne-1,2,np-1) = ac(ne-1,2,np-1) + as(ne-1,2,np-1) + aq(ne-1,2,np-1) + + & 4.*ae(ne-1,2,np-1)/(3.+deta) + ae(ne-1,2,np-1) = 0. + as(ne-1,2,np-1) = 0. + aq(ne-1,2,np-1) = 0. +c +c i=2, j=nq-1, k=np-1: + ae(2,nq-1,np-1) = ae(2,nq-1,np-1) + aw(2,nq-1,np-1) + ac(2,nq-1,np-1) = ac(2,nq-1,np-1) + an(2,nq-1,np-1) + aq(2,nq-1,np-1) + aw(2,nq-1,np-1) = 0. + an(2,nq-1,np-1) = 0. + aq(2,nq-1,np-1) = 0. +c +c i=ne-1, j=nq-1, k=np-1: + aw(ne-1,nq-1,np-1) = aw(ne-1,nq-1,np-1) - ae(ne-1,nq-1,np-1)/(3.+deta) + ac(ne-1,nq-1,np-1) = ac(ne-1,nq-1,np-1) + an(ne-1,nq-1,np-1) + + & aq(ne-1,nq-1,np-1) + 4.*ae(ne-1,nq-1,np-1)/(3.+deta) + ae(ne-1,nq-1,np-1) = 0. + an(ne-1,nq-1,np-1) = 0. + aq(ne-1,nq-1,np-1) = 0. +c +c Solve for psi: +c + call bicgst3d(ac,ae,aw,an,as,aq,ab,psiprim,rhs,dbh_eps,rmax,ier, + $ ne,nq,np) +c + if (rmax.gt.dbh_eps) then + write(*,*) '***WARNING: bicgst3d did not converge.' + endif + if (ier.eq.-1) then + write(*,*) '***WARNING: ier=-1' + endif + print *,'psiprim = ',maxval(psiprim),' ',minval(psiprim) +c +c Now, apply boundary conditions to psiprim: +c + do k = 1,np + do j = 1,nq + psiprim(1,j,k) = psiprim(3,j,k) + psiprim(ne,j,k) = (4.*psiprim(ne-1,j,k)-psiprim(ne-2,j,k))/ + $ (3.+deta) + enddo + enddo + do k = 1,np + do i = 1,ne + psiprim(i,1,k) = psiprim(i,2,k) + psiprim(i,nq,k) = psiprim(i,nq-1,k) + enddo + enddo + do j = 1,nq + do i = 1,ne + psiprim(i,j,1) = psiprim(i,j,2) + psiprim(i,j,np) = psiprim(i,j,np-1) + enddo + enddo +c +c Here, compute the derivatives of the spherical conformal factor +c +c goto 110 + + do k = 1,np + do j = 1,nq + do i = 2,ne-1 + detapsisph(i,j,k)=0.5*(psiprim(i+1,j,k)-psiprim(i-1,j,k)) + $ /deta + sinh(0.5*etagrd(i)) + enddo + detapsisph(1,j,k) = -detapsisph(3,j,k) + enddo + enddo +c + do k = 1,np + do j = 2,nq-1 + do i = 1,ne + dqpsisph(i,j,k)=0.5*(psiprim(i,j+1,k)-psiprim(i,j-1,k))/ + $ dq + enddo + enddo + do i = 1,ne + dqpsisph(i,1,k) = -dqpsisph(i,2,k) + dqpsisph(i,nq,k) = -dqpsisph(i,nq-1,k) + enddo + enddo +c + do k = 2,np-1 + do j = 1,nq + do i = 1,ne + dphipsisph(i,j,k)=0.5*(psiprim(i,j,k+1)-psiprim(i,j,k-1)) + $ /dphi + enddo + enddo + enddo + do j = 1,nq + do i = 1,ne + dphipsisph(i,j,1) = -dphipsisph(i,j,2) + dphipsisph(i,j,np) = -dphipsisph(i,j,np-1) + enddo + enddo +c + do k = 1,np + do j = 1,nq + do i = 2,ne-1 + detaetapsisph(i,j,k)=(psiprim(i+1,j,k)-2.*psiprim(i,j,k)+ + & psiprim(i-1,j,k))/deta**2+sqrt(0.25)* + & cosh(0.5*etagrd(i)) + enddo + detaetapsisph(1,j,k) = detaetapsisph(3,j,k) + enddo + enddo +c + do k = 1,np + do j = 2,nq-1 + do i = 1,ne + detaqpsisph(i,j,k)=0.5*(detapsisph(i,j+1,k)-detapsisph(i, + $ j-1,k))/dq + enddo + enddo + do i = 1,ne + detaqpsisph(i,1,k) = -detaqpsisph(i,2,k) + detaqpsisph(i,nq,k) = -detaqpsisph(i,nq-1,k) + enddo + enddo +c + do k = 2,np-1 + do j = 1,nq + do i = 1,ne + detaphipsisph(i,j,k)=0.5*(detapsisph(i,j,k+1)-detapsisph( + $ i,j,k-1))/dphi + enddo + enddo + enddo + do j = 1,nq + do i = 1,ne + detaphipsisph(i,j,1) = -detaphipsisph(i,j,2) + detaphipsisph(i,j,np) = -detaphipsisph(i,j,np-1) + enddo + enddo +c + do k = 1,np + do j = 2,nq-1 + do i = 1,ne + dqqpsisph(i,j,k)=0.5*(dqpsisph(i,j+1,k)-dqpsisph(i,j-1,k))/ + $ dq + enddo + enddo + do i = 1,ne + dqqpsisph(i,1,k) = dqqpsisph(i,2,k) + dqqpsisph(i,nq,k) = dqqpsisph(i,nq-1,k) + enddo + enddo +c + do k = 2,np-1 + do j = 1,nq + do i = 1,ne + dqphipsisph(i,j,k)=0.5*(dqpsisph(i,j,k+1)-dqpsisph(i,j,k-1))/ + $ dphi + enddo + enddo + enddo + do j = 1,nq + do i = 1,ne + dqphipsisph(i,j,1) = -dqphipsisph(i,j,2) + dqphipsisph(i,j,np) = -dqphipsisph(i,j,np-1) + enddo + enddo +c + do k = 2,np-1 + do j = 1,nq + do i = 1,ne + dphiphipsisph(i,j,k)=0.5*(dphipsisph(i,j,k+1)- + $ dphipsisph(i,j,k-1))/dphi + enddo + enddo + enddo + do j = 1,nq + do i = 1,ne + dphiphipsisph(i,j,1) = dphiphipsisph(i,j,2) + dphiphipsisph(i,j,np) = dphiphipsisph(i,j,np-1) + enddo + enddo +c + do k = 1,np + do j = 1,nq + psisph(:,j,k)=psiprim(:,j,k)+2.0*cosh(0.5*etagrd) + enddo + enddo + +c +c Now compute on the Cartesian coordinate. +c +c Compute eta,q,phi at the each points of cartesian grid + + eta = 0.5d0*dlog(x**2+y**2+z**2) + abseta = abs(eta) + q = datan2(sqrt(x**2+y**2),z) + phi = datan2(y,x) + + do k = 1,nz + do j = 1,ny + do i = 1,nx + + if(eta(i,j,k) .lt. 0)then + sign_eta(i,j,k) = -1.0d0 + else + sign_eta(i,j,k) = 1.0d0 + endif + + enddo + enddo + enddo + + call CCTK_GetInterpHandle (handle, "simple_local") + + npoints = nx*ny*nz + + call CCTK_Interp (ierror,cctkGH,handle,npoints,3,10,10, + $ ne,nq,np,abseta,q,phi, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ etagrd(1),qgrd(1),phigrd(1)-pi,deta,dq,dphi, + $ psisph,detapsisph,dqpsisph,dphipsisph,detaetapsisph, + $ detaqpsisph,detaphipsisph,dqqpsisph,dqphipsisph,dphiphipsisph, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ CCTK_VARIABLE_REAL, + $ psi3d,detapsi3d,dqpsi3d,dphipsi3d,detaetapsi3d,detaqpsi3d, + $ detaphipsi3d,dqqpsi3d,dqphipsi3d,dphiphipsi3d, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, + $ CCTK_VARIABLE_REAL) + + psi = psi3d*exp(-0.5*eta) + detapsi3d = sign_eta*detapsi3d + detaqpsi3d = sign_eta*detaqpsi3d + detaphipsi3d = sign_eta*detaphipsi3d + + do k=1,nz + do j=1,ny + do i=1,nx +c psix = \partial psi / \partial x / psi +#include "CactusEinstein/DistortedBHIVP/src/psi_1st_deriv.x" +c psixx = \partial^2\psi / \partial x^2 / psi +#include "CactusEinstein/DistortedBHIVP/src/psi_2nd_deriv.x" + enddo + enddo + enddo + +c Conformal metric +c gxx = ... + +c Derivatives of the metric +c dxgxx = 1/2 \partial gxx / \partial x +c + do k=1,nz + do j=1,ny + do i=1,nx +#include "CactusEinstein/DistortedBHIVP/src/gij.x" + enddo + enddo + enddo + +c Courvature + kxx = 0.0d0 + kxy = 0.0d0 + kxz = 0.0d0 + kyy = 0.0d0 + kyz = 0.0d0 + kzz = 0.0d0 + + 110 continue +c Set ADM mass + i = ne-15 + adm = 0.0 + do k=2,np-1 + do j=2,nq-1 + adm=adm+(psisph(i,j,k)-(psisph(i+1,j,k)-psisph(i-1,j,k))/ + $ deta)*exp(0.5*etagrd(i)) + enddo + enddo + adm=adm/(nq-2)/(np-2) + print *,'ADM mass: ',adm + if (CCTK_Equals(initial_lapse,"schwarz")==1) then + write (*,*)"Initial with schwarzschild-like lapse" + write (*,*)"using alp = (2.*r - adm)/(2.*r+adm)." + alp = (2.*r - adm)/(2.*r+adm) + endif + + deallocate(ac,ae,aw,an,as,aq,ab,rhs, + $ qf,qfetaeta,qfqq,qfphi,qfphiphi, + $ psisph,psiprim,detapsisph,dqpsisph,dphipsisph,detaetapsisph, + $ detaqpsisph,detaphipsisph,dqqpsisph,dqphipsisph, + $ dphiphipsisph, + $ psi3d,detapsi3d,dqpsi3d, + $ dphipsi3d,detaetapsi3d,detaqpsi3d, + $ detaphipsi3d,dqqpsi3d,dqphipsi3d, + $ dphiphipsi3d) + deallocate(etagrd,qgrd,phigrd) + deallocate(eta,abseta,sign_eta,q,phi) + + return + end |