diff options
author | tradke <tradke@0a4070d5-58f5-498f-b6c0-2693e757fa0f> | 2000-11-21 10:59:18 +0000 |
---|---|---|
committer | tradke <tradke@0a4070d5-58f5-498f-b6c0-2693e757fa0f> | 2000-11-21 10:59:18 +0000 |
commit | fa69a46020cef9636482c9d45f40f7115eaff740 (patch) | |
tree | bf72fecdb478c01120a3d763cbbf5267a579ed5f | |
parent | a1c689dac4969b9dd76742267db33c12d318e265 (diff) |
Don't change parameters in place but assign their values to local variables
instead. This fixes problems with recovery where the modified parameter values
caused different results.
git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/IDAxiBrillBH/trunk@18 0a4070d5-58f5-498f-b6c0-2693e757fa0f
-rw-r--r-- | src/IDAxiBrillBH.F | 99 |
1 files changed, 52 insertions, 47 deletions
diff --git a/src/IDAxiBrillBH.F b/src/IDAxiBrillBH.F index 2c7f913..5535734 100644 --- a/src/IDAxiBrillBH.F +++ b/src/IDAxiBrillBH.F @@ -68,6 +68,7 @@ c Perhaps this and others should go into cctk.h integer :: nx,ny,nz integer i,j,k,nquads integer npoints,handle,ierror + integer neb, nqb pi = 4.0d0*atan(1.0d0) @@ -86,13 +87,17 @@ c Brill wave parameters c Solve on this sized cartesian grid c 2D grid size NE x NQ c Add 2 zones for boundaries... - ne = ne+2 - nq = nq+2 +c 21/11/00 TR: dont change parameters in place +c but keep a copy in local variables +c Otherwise the changed parameters cause trouble +c after recovery. + neb = ne+2 + nqb = nq+2 ! do I need to call free? - allocate(cc(ne,nq),ce(ne,nq),cw(ne,nq),cn(ne,nq),cs(ne,nq), - $ rhs(ne,nq),psi2d(ne,nq),detapsi2d(ne,nq),dqpsi2d(ne,nq), - $ detaetapsi2d(ne,nq),detaqpsi2d(ne,nq),dqqpsi2d(ne,nq), - $ etagrd(ne),qgrd(nq)) + allocate(cc(neb,nqb),ce(neb,nqb),cw(neb,nqb),cn(neb,nqb),cs(neb,nqb), + $ rhs(neb,nqb),psi2d(neb,nqb),detapsi2d(neb,nqb),dqpsi2d(neb,nqb), + $ detaetapsi2d(neb,nqb),detaqpsi2d(neb,nqb),dqqpsi2d(neb,nqb), + $ etagrd(neb),qgrd(nqb)) allocate(eta(nx,ny,nz),abseta(nx,ny,nz),sign_eta(nx,ny,nz), $ q(nx,ny,nz),phi(nx,ny,nz), $ psi2dv(nx,ny,nz),detapsi2dv(nx,ny,nz),dqpsi2dv(nx,ny,nz), @@ -103,37 +108,37 @@ c Initialize some arrays detapsi2d = 0.0d0 nquads = 2 - dq = nquads*0.5*pi/(nq-2) - deta = etamax/(ne-3) + dq = nquads*0.5*pi/(nqb-2) + deta = etamax/(neb-3) - do j=1,nq + do j=1,nqb qgrd(j) = (j-1.5)*dq - do i=1,ne + do i=1,neb etagrd(i) = (i-2)*deta #include "CactusEinstein/IDAxiBrillBH/src/bhbrill.x" enddo enddo c Boundary conditions - do j=1,nq + do j=1,nqb ce(2,j)=ce(2,j)+cw(2,j) cw(2,j)=0.0 - cw(ne-1,j)=cw(ne-1,j)+ce(ne-1,j) - cc(ne-1,j)=cc(ne-1,j)-deta*ce(ne-1,j) - ce(ne-1,j)=0.0 + cw(neb-1,j)=cw(neb-1,j)+ce(neb-1,j) + cc(neb-1,j)=cc(neb-1,j)-deta*ce(neb-1,j) + ce(neb-1,j)=0.0 enddo - do i=1,ne + do i=1,neb cc(i,2)=cc(i,2)+cs(i,2) cs(i,2)=0.0 - cc(i,nq-1)=cc(i,nq-1)+cn(i,nq-1) - cn(i,nq-1)=0.0 + cc(i,nqb-1)=cc(i,nqb-1)+cn(i,nqb-1) + cn(i,nqb-1)=0.0 enddo c Do the solve print *, " Calling axisymmetric solver" - call mgparm (levels,5,id5,id9,idi,idg,ne,nq) - call mg5 (ne,2,ne-1,nq,2,nq-1, + call mgparm (levels,5,id5,id9,idi,idg,neb,nqb) + call mg5 (neb,2,neb-1,nqb,2,nqb-1, $ cc,cn,cs,cw,ce,psi2d,rhs, $ id5,id9,idi,idg,1,axibheps,rmax,ier) print *, " Solve complete" @@ -147,8 +152,8 @@ c be called if the IVP solve is not successful print *,'psi2d = ',maxval(psi2d),' ',minval(psi2d) ep2 = 0.0 - do j=2,nq-1 - do i=2,ne-1 + do j=2,nqb-1 + do i=2,neb-1 ep1 = rhs(i,j)-psi2d(i,j)*cc(i,j)-psi2d(i,j+1)*cn(i,j)-psi2d(i,j-1)*cs(i,j)- & psi2d(i+1,j)*ce(i,j)-psi2d(i-1,j)*cw(i,j) ep2 = max(abs(ep1),ep2) @@ -157,17 +162,17 @@ c be called if the IVP solve is not successful print *,'Resulting eps =',ep2 ! what a pain all this is.... - do j=1,nq + do j=1,nqb psi2d(1,j)=psi2d(3,j) - psi2d(ne,j)=-deta*psi2d(ne-1,j)+psi2d(ne-2,j) + psi2d(neb,j)=-deta*psi2d(neb-1,j)+psi2d(neb-2,j) enddo - do i=1,ne + do i=1,neb psi2d(i,1)=psi2d(i,2) - psi2d(i,nq)=psi2d(i,nq-1) + psi2d(i,nqb)=psi2d(i,nqb-1) enddo c goto 111 - do j=2,nq-1 - do i=2,ne-1 + do j=2,nqb-1 + do i=2,neb-1 dqpsi2d(i,j)=0.5*(psi2d(i,j+1)-psi2d(i,j-1))/dq dqqpsi2d(i,j)=(psi2d(i,j+1)+psi2d(i,j-1)-2.*psi2d(i,j))/dq**2 detapsi2d(i,j)=sinh(0.5*etagrd(i))+0.5*(psi2d(i+1,j)-psi2d(i-1,j))/deta @@ -175,46 +180,46 @@ c goto 111 $ (psi2d(i+1,j)+psi2d(i-1,j)-2.*psi2d(i,j))/deta**2 enddo enddo - do j=1,nq + do j=1,nqb detapsi2d(1,j)=-detapsi2d(3,j) - detapsi2d(ne,j)=detapsi2d(ne-2,j) ! simplified + detapsi2d(neb,j)=detapsi2d(neb-2,j) ! simplified detaetapsi2d(1,j)=detaetapsi2d(3,j) - detaetapsi2d(ne,j)=detaetapsi2d(ne-2,j) ! simplified... + detaetapsi2d(neb,j)=detaetapsi2d(neb-2,j) ! simplified... dqqpsi2d(1,j)=dqqpsi2d(3,j) - dqqpsi2d(ne,j)=dqqpsi2d(ne-2,j) ! simplified + dqqpsi2d(neb,j)=dqqpsi2d(neb-2,j) ! simplified dqpsi2d(1,j)=dqpsi2d(3,j) - dqpsi2d(ne,j)=-dq*dqpsi2d(ne-1,j)+dqpsi2d(ne-2,j) + dqpsi2d(neb,j)=-dq*dqpsi2d(neb-1,j)+dqpsi2d(neb-2,j) enddo - do i=1,ne + do i=1,neb detapsi2d(i,1)=detapsi2d(i,2) - detapsi2d(i,nq)=detapsi2d(i,nq-1) + detapsi2d(i,nqb)=detapsi2d(i,nqb-1) detaetapsi2d(i,1)=detaetapsi2d(i,2) - detaetapsi2d(i,nq)=detaetapsi2d(i,nq-1) + detaetapsi2d(i,nqb)=detaetapsi2d(i,nqb-1) dqqpsi2d(i,1)=dqqpsi2d(i,2) - dqqpsi2d(i,nq)=dqqpsi2d(i,nq-1) + dqqpsi2d(i,nqb)=dqqpsi2d(i,nqb-1) dqpsi2d(i,1)=-dqpsi2d(i,2) - dqpsi2d(i,nq)=-dqpsi2d(i,nq-1) + dqpsi2d(i,nqb)=-dqpsi2d(i,nqb-1) enddo - do j=2,nq-1 - do i=2,ne-1 + do j=2,nqb-1 + do i=2,neb-1 detaqpsi2d(i,j)=0.5*(detapsi2d(i,j+1)-detapsi2d(i,j-1))/dq enddo enddo - do j=1,nq + do j=1,nqb detaqpsi2d(1,j)=-detaqpsi2d(3,j) - detaqpsi2d(ne,j)=detaqpsi2d(ne-2,j) ! simplified + detaqpsi2d(neb,j)=detaqpsi2d(neb-2,j) ! simplified enddo do i=1,ne detaqpsi2d(i,1)=-detaqpsi2d(i,2) - detaqpsi2d(i,nq)=-detaqpsi2d(i,nq-1) + detaqpsi2d(i,nqb)=-detaqpsi2d(i,nqb-1) enddo - do j=1,nq + do j=1,nqb psi2d(:,j)=psi2d(:,j)+2.*cosh(0.5*etagrd) enddo @@ -252,7 +257,7 @@ c phi(i,j,k)= datan2(y(i,j,k),x(i,j,k)) npoints = nx*ny*nz call CCTK_Interp (ierror,cctkGH,handle,npoints,2,6,6, - $ ne,nq,abseta,q, + $ neb,nqb,abseta,q, $ CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL, $ etagrd(1),qgrd(1),deta,dq, $ psi2d,detapsi2d,dqpsi2d,detaetapsi2d,detaqpsi2d,dqqpsi2d, @@ -313,13 +318,13 @@ c Curvature 111 continue c Set ADM mass - i = ne-15 + i = neb-15 adm = 0.0 - do j=2,nq-1 + do j=2,nqb-1 adm=adm+(psi2d(i,j)-(psi2d(i+1,j)-psi2d(i-1,j))/deta)*exp(0.5* $ etagrd(i)) enddo - adm=adm/(nq-2) + adm=adm/(nqb-2) print *,'ADM mass: ',adm if (CCTK_Equals(initial_lapse,"schwarz")==1) then write (*,*)"Initial with schwarzschild-like lapse" |