From 605e01b46d8b0ee092f9a172a07e5874ea91c864 Mon Sep 17 00:00:00 2001 From: ryoji Date: Mon, 14 Feb 2000 16:05:57 +0000 Subject: remove from here. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/DistortedBHIVP/trunk@11 971fb155-194f-0410-9daf-e2eca44e59f5 --- src/DistortedBHIVP.F.sav | 681 ----------------------------------------------- src/interp3.F | 264 ------------------ 2 files changed, 945 deletions(-) delete mode 100644 src/DistortedBHIVP.F.sav delete mode 100644 src/interp3.F 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 diff --git a/src/interp3.F b/src/interp3.F deleted file mode 100644 index 0bbdf90..0000000 --- a/src/interp3.F +++ /dev/null @@ -1,264 +0,0 @@ -c /*@@ -c @routine interp3 -c @date Fri Feb 14 08:46:53 1997 -c @author Ryoji Takahashi -c @desc -c Interpolates from 3D data var with coordinates x, y, and z and -c sizes nx , ny, and nz onto 1D data out with position outx, outy -c ,and outz nout points. -c

-c This has linear, quadratic and cubic interpolators in it. -c Or will one day very soon. -c @enddesc -c @calls -c @calledby numerical_nonaxi -c@@*/ - - subroutine interp3(dat,x,y,z,nx,ny,nz,out,outx,outy,outz,nout - $ ,order) - implicit none - integer nx,ny,nz,nout - real*8 dat(nx,ny,nz), x(nx), y(ny), z(nz) - real*8 out(nout),outx(nout),outy(nout),outz(nout) - integer order -c Interpolation goes from ibelow to ibelow+[1,2,3] depending on order - integer i,j,k,l,m,ibelow,jbelow,kbelow,pt - real*8 xsym,ysym,zsym,findx,findy,findz,frac - real*8 ydir(order+1) - real*8 zdir(order+1) - real*8 f(4), fp(4,4), fl(4) - real*8 dbh_linear, dbh_quad, dbh_cubic - real*8 dx, dy, dz, PI - integer twobhjsad - - pi = 4.0d0*atan(1.0d0) - - dx = x(2) - x(1) - dy = y(2) - y(1) - dz = z(2) - z(1) - -c Loop over all out points - do pt=1,nout - zsym = 1.0D0 - ysym = 1.0D0 - xsym = 1.0D0 -c Check bounds - findx = outx(pt) - if (findx .lt. x(1)) then - write (*,*) "Below inner bound at ",pt,outx(pt) - STOP - endif - if (findx .gt. x(nx)) then - write (*,*) "Above x bounds at ",pt,outx(pt),x(nx) - STOP - endif - findy = outy(pt) - if (findy .lt. y(1)) then - write (*,*) "Below y inner bound at ",pt,outy(pt) - STOP - endif - if (findy .gt. y(ny)) then - write (*,*) "Below y inner bound at ",pt,outy(pt) - STOP - endif - findz = outz(pt)+pi - if (findz .lt. z(1)) then - write (*,*) "Below z inner bound at ",pt,outz(pt) - STOP - endif - if (findz .gt. z(nz)) then - write (*,*) "Below z inner bound at ",pt,outz(pt) - STOP - endif - -c Locate ourselves in i,j space -c do i=1,nx -c if (x(i) .lt. findx) then -c ibelow = i -c endif -c enddo -c Assume a regular grid - ibelow = (findx-x(1))/dx+1 - - if (order .eq. 3 .and. ibelow .gt. 1) then - ibelow = ibelow - 1 - endif - if (ibelow + order .gt. nx) then - ibelow = nx - order - endif - -c do i=1,ny -c if (y(i) .lt. findy) then -c jbelow = i -c endif -c enddo -c Assume a regular grid - jbelow = (findy-y(1))/dy+1 - if (order .eq. 3 .and. jbelow .gt. 1) then - jbelow = jbelow - 1 - endif - if (jbelow + order .gt. ny) then - jbelow = ny - order - endif - -c do i=1,nz -c if (z(i) .lt. findz) then -c kbelow = i -c endif -c enddo -c Assume a regular grid - kbelow = (findz-z(1))/dz+1 - if (order .eq. 3 .and. kbelow .gt. 1) then - kbelow = kbelow - 1 - endif - if (kbelow + order .gt. nz) then - kbelow = nz - order - endif - -c write (*,*) "PT :",findx,findy -c write (*,*) "SYM:",sym -c write (*,*) "BOUND X ",ibelow,x(ibelow),x(ibelow+1) -c write (*,*) "BOUND Y ",jbelow,y(jbelow),y(jbelow+1) -c write (*,*) "BOUND z ",kbelow,z(kbelow),z(kbelow+1) - -c So do the interpolation - if (order .eq. 1) then -c Interp in the x direction -! frac = (findx-x(ibelow))/(x(ibelow+1)-x(ibelow)) - do l = 1,2 - do m = 1,2 - f(1) = dat(ibelow,jbelow+l-1,kbelow+m-1) - f(2) = dat(ibelow+1,jbelow+l-1,kbelow+m-1) - fp(l,m) = dbh_linear(f,x(ibelow),dx,findx) - enddo - enddo -c Now take our 2x2 plane and interp to the center of both -c in the y direction -! frac = (findy-y(jbelow))/(y(jbelow+1)-y(jbelow)) - do m = 1,2 - f(1) = fp(1,m) - f(2) = fp(2,m) - fl(m) = dbh_linear(f,y(jbelow),dy,findy) - enddo -c And finally, interp in the z direction - out(pt) = xsym * ysym * zsym * - $ dbh_linear(fl,z(kbelow),dz,findz) - else if (order .eq. 2) then -c Load up for calls to poly2inter - do l=1,3 - do m=1,3 - f(1) = dat(ibelow,jbelow+l-1,kbelow+m-1) - f(2) = dat(ibelow+1,jbelow+l-1,kbelow+m-1) - f(3) = dat(ibelow+2,jbelow+l-1,kbelow+m-1) - fp(l,m) = dbh_quad(f,x(ibelow),dx,findx) - enddo - enddo -c Now take our 2x2 plane and interp to the center of both -c in the y direction - do m=1,3 - f(1) = fp(1,m) - f(2) = fp(2,m) - f(3) = fp(3,m) - fl(m) = dbh_quad(f,y(jbelow),dy,findy) - enddo -c And finally, interp in the z direction - out(pt) = xsym * ysym * zsym * - $ dbh_quad(fl,z(kbelow),dz,findz) - else if (order .eq. 3) then -c Load up for calls to cubic - do l=1,4 - do m=1,4 - f(1) = dat(ibelow,jbelow+l-1,kbelow+m-1) - f(2) = dat(ibelow+1,jbelow+l-1,kbelow+m-1) - f(3) = dat(ibelow+2,jbelow+l-1,kbelow+m-1) - f(4) = dat(ibelow+3,jbelow+l-1,kbelow+m-1) - fp(l,m) = dbh_cubic(f,x(ibelow),dx,findx) - enddo - enddo -c Now take our 2x2 plane and interp to the center of both -c in the y direction - do m=1,4 - f(1) = fp(1,m) - f(2) = fp(2,m) - f(3) = fp(3,m) - f(4) = fp(4,m) - fl(m) = dbh_cubic(f,y(jbelow),dy,findz) - enddo -c And finally, interp in the z direction - out(pt) = xsym * ysym * zsym * - $ dbh_cubic(fl,z(kbelow),dz,findz) - else - write (*,*) "ORDER set wrong in interp3d",order - stop - endif - - enddo - return - end - - real*8 function dbh_linear(f, x0, dx, findx) - implicit none - real*8 f(2),x0,dx,findx - real*8 frac - - frac = (findx-x0)/dx - dbh_linear = (frac)*f(2) + (1.0-frac)*f(1) - - return - end - - real*8 function dbh_quad(f, x0, dx, findx) - implicit none - real*8 f(3),x0, dx, findx, dbh_quad - real*8 f0,f1,f2 - real*8 a,b,c, dx2, x02, o2dx2 -c Mathematica tells us -c - List(List(Rule(c,(2*dx**2*f0 + 3*dx*f0*x0 - 4*dx*f1*x0 -c - + dx*f2*x0 + f0*x0**2 - 2*f1*x0**2 + f2*x0**2) -c - /(2*dx**2)), Rule(b,(-3*dx*f0 + 4*dx*f1 - dx*f2 - 2*f0*x0 -c - + 4*f1*x0 - 2*f2*x0)/(2*dx**2)),Rule(a,(f0 - 2*f1 + -c - f2)/(2*dx**2)))) - - f0 = f(1) - f1 = f(2) - f2 = f(3) - dx2 = dx**2 - x02 = x0**2 - o2dx2 = 1.0D0/(2.0D0*dx2) - - c = (2.0D0*dx2*f0 + dx*x0*(3.0D0*f0 - 4.0D0*f1 + f2) + - $ x02*(f0 - 2.0D0*f1 + f2))*o2dx2 - b = (dx * (-3.0D0*f0 + 4.0D0*f1 - f2) + x0 * (- 2.0D0*f0 + - $ 4.0D0*f1 - 2.0D0*f2))*o2dx2 - - a = (f0 - 2.0D0*f1 + f2)*o2dx2 - - dbh_quad = a*findx**2 + b*findx + c - - end - - - real*8 function dbh_cubic(f, x0, dx, findx) - implicit none - real*8 a,b,c,d,dbh_cubic - real*8 f(4),x0,dx,findx - - a = -(f(1)-3.0*f(2)+3.0*f(3)-f(4)) / (6.0*(dx**3)) - - b = (f(1)-2.0*f(2)+f(3))/(2.0*(dx**2)) + - $ (f(1)-3.0*f(2)+3.0*f(3)-f(4))*(dx+x0)/(2.0*(dx**3)) - - c = ((dx**2)*(-11.0*f(1) + 18.0*f(2) - 9.0*f(3) + 2.0*f(4)) + - $ dx*x0* (-12.0*f(1) + 30.0*f(2) - 24.0*f(3) + 6.0*f(4)) + - $ (x0**2)*( -3.0*f(1) + 9.0*f(2) - 9.0*f(3) + 3.0*f(4))) / - $ (6.0*(dx**3)) - - d = ((dx**3)* ( 6.0*f(1) ) + - $ (dx**2)*x0*( 11.0*f(1) - 18.0*f(2) + 9.0*f(3) - 2.0*f(4)) + - $ (x0**2)*dx*( 6.0*f(1) - 15.0*f(2) + 12.0*f(3) - 3.0*f(4)) + - $ (x0**3)* ( 1.0*f(1) - 3.0*f(2) + 3.0*f(3) - 1.0*f(4)))/ - $ (6.0*(dx**3)) - - dbh_cubic = ((a*findx + b)*findx + c)*findx + d - return - end -- cgit v1.2.3