diff options
Diffstat (limited to 'src/DistortedBHIVP.F.sav')
-rw-r--r-- | src/DistortedBHIVP.F.sav | 681 |
1 files changed, 0 insertions, 681 deletions
diff --git a/src/DistortedBHIVP.F.sav b/src/DistortedBHIVP.F.sav deleted file mode 100644 index 815d25d..0000000 --- a/src/DistortedBHIVP.F.sav +++ /dev/null @@ -1,681 +0,0 @@ -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 :: amp,eta0,sigma,c,etamax,deta,dq,dphi - integer :: ne,nq,np,n - 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 - amp = distortedbh_amp - eta0 = distortedbh_eta0 - c = distortedbh_c - sigma = distortedbh_sigma - n = distortedbh_n - etamax = distortedbh_etamax - - 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 = distortedbh_ne + 2 - nq = distortedbh_nq + 2 - np = distortedbh_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 - do j = 1,nq - qgrd(j) = (j-1.5)*dq - do i=1,ne - etagrd(i) = (i-2)*deta - enddo - enddo - 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 - do k = 1,np - do j = 1,nq - do i = 1,ne - psisph(i,j,k) = 2.*cosh(0.5*etagrd(i)) - psiprim(i,j,k) = 0. - enddo - enddo - enddo -c -c Initialize stencil coefficients: -c - do k = 1,np - do j = 1,nq - do i = 1,ne - ac(i,j,k) = 0. - ae(i,j,k) = 0. - aw(i,j,k) = 0. - an(i,j,k) = 0. - as(i,j,k) = 0. - aq(i,j,k) = 0. - ab(i,j,k) = 0. - rhs(i,j,k) = 0. - enddo - enddo - enddo -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 - do k = 1,nz - do j = 1,ny - do i = 1,nx - eta(i,j,k) = 0.5d0*dlog(x(i,j,k)**2+y(i,j,k)**2+z(i,j,k)**2) - abseta(i,j,k) = abs(eta(i,j,k)) - if(eta(i,j,k) .lt. 0)then - sign_eta(i,j,k) = -1.0d0 - else - sign_eta(i,j,k) = 1.0d0 - endif - q(i,j,k) = datan2(sqrt(x(i,j,k)**2+y(i,j,k)**2),z(i,j,k)) - phi(i,j,k) = datan2(y(i,j,k),x(i,j,k)) - enddo - enddo - enddo - - call CCTK_GetInterpHandle (handle, "simple_local") - - npoints = nx*ny*nz - - call CCTK_Interp (cctkGH,ierror,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) - - do k=1,nz - do j=1,ny - do i=1,nx - psi(i,j,k) = psi3d(i,j,k)*exp(-0.5*eta(i,j,k)) - - detapsi3d(i,j,k) = sign_eta(i,j,k)*detapsi3d(i,j,k) - -c psix = \partial psi / \partial x / psi -#include "CactusEinstein/DistortedBHIVP/src/psi_1st_deriv.x" - - detaqpsi3d(i,j,k) = sign_eta(i,j,k)*detaqpsi3d(i,j,k) - detaphipsi3d(i,j,k) = sign_eta(i,j,k)*detaphipsi3d(i,j,k) - -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 - eta(i,j,k) = 0.5d0*dlog(x(i,j,k)**2+y(i,j,k)**2+z(i,j,k)**2) - q(i,j,k) = datan2(sqrt(x(i,j,k)**2+y(i,j,k)**2),z(i,j,k)) - phi(i,j,k) = datan2(y(i,j,k),x(i,j,k)) -#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 |