From 2c5a9898aade0ef52e6208bbb3fca110e66b742d Mon Sep 17 00:00:00 2001 From: allen Date: Tue, 11 Jan 2000 09:11:58 +0000 Subject: Fixed bug in extract on multiprocessors, and tidied up especially reorganising the Extract.F file git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/Extract/trunk@9 5301f0c2-dbc4-4cee-b2f5-8d7afba4d129 --- src/D3_extract.F | 10 +-- src/D3_to_D2.F | 11 +-- src/Extract.F | 266 ++++++++++++++++++------------------------------------- src/energy.f | 13 +-- 4 files changed, 102 insertions(+), 198 deletions(-) (limited to 'src') diff --git a/src/D3_extract.F b/src/D3_extract.F index f5f6850..854ae1d 100644 --- a/src/D3_extract.F +++ b/src/D3_extract.F @@ -169,15 +169,11 @@ c ------------------------------------------------------------------ & ADMmass_int1,ADMmass_int2,momentum_int1,momentum_int2, & momentum_int3,spin_int1,spin_int2,spin_int3,Extract_temp3d) -c ------------------------------------------------------------------ -c Now this is done the grid on processor 0 has correct values and the -c grid on all other processors have junk. This means that processor -c 0 needs to calculate the wave forms and the other processors have -c to sit and do nothing. So we shall put this in a conditional -c to only happen on processor 0. -c PW Jan 31 98 if (myproc .eq. 0) then + + print *,psis(5,5) + c ------------------------------------------------------------------ c c 3. Transform tensor components to polar coordinates on 2-surface diff --git a/src/D3_to_D2.F b/src/D3_to_D2.F index 05e1f57..74c3cdb 100644 --- a/src/D3_to_D2.F +++ b/src/D3_to_D2.F @@ -113,11 +113,12 @@ c Choose the interpolator c ----------------------- call CCTK_GetInterpHandle (handle, "simple") -c Find the origin of the spatial coordinates -c ------------------------------------------ - call CCTK_CoordRange(ierror,cctkGH,xmin,xmax,"x") - call CCTK_CoordRange(ierror,cctkGH,ymin,ymax,"y") - call CCTK_CoordRange(ierror,cctkGH,zmin,zmax,"z") + +c Find the local origin of the spatial coordinates +c ------------------------------------------------ + call CCTK_CoordLocalRange(ierror,cctkGH,xmin,xmax,"x") + call CCTK_CoordLocalRange(ierror,cctkGH,ymin,ymax,"y") + call CCTK_CoordLocalRange(ierror,cctkGH,zmin,zmax,"z") c Project un-physical metric and conformal factor onto sphere diff --git a/src/Extract.F b/src/Extract.F index 1aa8534..761af55 100644 --- a/src/Extract.F +++ b/src/Extract.F @@ -61,7 +61,7 @@ c Local variables INTEGER :: & ix,iy,iz,fn1,fn2,lmin,lmax,mmin,mmax,lstep,mstep, & il,im,it,ndet,idet,nx_parfile,nchar,iconv,ioutput - INTEGER,SAVE :: openfiles = 1 + INTEGER,SAVE :: openfile = 1 INTEGER,PARAMETER :: & max_detectors = 9,number_timecoords = 2 CCTK_REAL,PARAMETER :: @@ -91,10 +91,6 @@ c ------------------------------------------------------------------ c Get the output directory sorted out call CCTK_FortranString(nchar,outdir,filestr) - if (firsttime) then - call CCTK_mkdir(ierr,filestr); - firsttime = .FALSE. - end if c ------------------------------------------------------------------ @@ -307,33 +303,35 @@ c Open output files test_myproc1: IF (CCTK_MyProc(cctkGH) == 0 ) THEN out1 = filestr(1:nchar)//"/rsch_ini.rl" - OPEN(UNIT= 77,file = out1,STATUS="unknown") + call Extract_open(0,out1,77) + out1 = filestr(1:nchar)//"/mass_ini.rl" - OPEN(UNIT= 78,file = out1,STATUS="unknown") + call Extract_open(0,out1,78) + IF (do_ADMmass(1) == 1) THEN out1 = filestr(1:nchar)//"/ADMmass_ini.rl" - OPEN(UNIT= 79,file = out1,STATUS="unknown") + call Extract_open(0,out1,79) ENDIF IF (do_ADMmass(2) == 1) THEN out1 = filestr(1:nchar)//"/ADMmassc_ini.rl" - OPEN(UNIT= 80,file = out1,STATUS="unknown") + call Extract_open(0,out1,80) ENDIF IF (do_momentum == 1) THEN out1 = filestr(1:nchar)//"/momentum_x_ini.rl" - OPEN(UNIT=781,file = out1,STATUS="unknown") + call Extract_open(0,out1,781) out1 = filestr(1:nchar)//"/momentum_y_ini.rl" - OPEN(UNIT=782,file = out1,STATUS="unknown") + call Extract_open(0,out1,782) out1 = filestr(1:nchar)//"/momentum_z_ini.rl" - OPEN(UNIT=783,file = out1,STATUS="unknown") + call Extract_open(0,out1,783) ENDIF IF (do_spin == 1) THEN out1 = filestr(1:nchar)//"/spin_x_ini.rl" - OPEN(UNIT=784,file = out1,STATUS="unknown") + call Extract_open(0,out1,784) out1 = filestr(1:nchar)//"/spin_y_ini.rl" - OPEN(UNIT=785,file = out1,STATUS="unknown") + call Extract_open(0,out1,785) out1 = filestr(1:nchar)//"/spin_z_ini.rl" - OPEN(UNIT=786,file = out1,STATUS="unknown") + call Extract_open(0,out1,786) ENDIF loop_l1: DO il = lmin,lmax,lstep @@ -346,18 +344,17 @@ c Open output files loop_m1: DO im = mmin,mmax,lstep - fn1 = il*10+im out1 = filestr(1:nchar)//"/Qeven_ini_"// & CHAR(il+48)//CHAR(im+48)//".rl" - OPEN(UNIT=fn1,FILE=out1,STATUS="unknown") + call Extract_open(0,out1,fn1) c Only print odd-parity if full grid IF (igrid == 0) THEN fn2 = 100+il*10+im out2 = filestr(1:nchar)//"/Qodd_ini_"// & CHAR(il+48)//CHAR(im+48)//".rl" - OPEN(UNIT=fn2,FILE=out2,STATUS="unknown") + call Extract_open(0,out2,fn2) END IF END DO loop_m1 @@ -374,17 +371,11 @@ c Find range of extraction radii dr = Cauchy_dr - c Do extraction at each radius IF (verbose == 1) THEN - WRITE(*,*) WRITE(*,*) "Extracting Cauchy initial data" - WRITE(*,*) "------------------------------" - WRITE(*,*) "From r =",r1 - WRITE(*,*) "To r =",r2 - WRITE(*,*) "In steps of ",dr - WRITE(*,*) + WRITE(*,*) " r = ",r1," to ",r2," step ",dr END IF radius = r1 @@ -400,25 +391,9 @@ c Do extraction at each radius IF (verbose == 1) THEN WRITE(*,*) "Extracted at r =",radius - WRITE(*,*) "Schwarzschild radius =",rsch - WRITE(*,*) "Schwarzschild mass =",mass - if (do_ADMmass(1) == 1) then - WRITE(*,*) "ADM mass =",ADMmass(1) - end if - if (do_ADMmass(2) == 1) then - WRITE(*,*) "ADM mass (conformal) =",ADMmass(2) - end if - if (do_momentum == 1) then - WRITE(*,*) "momentum = (",momentum(1),",", - & momentum(2),",",momentum(3),")" - end if - if (do_spin == 1) then - WRITE(*,*) "spin = (",spin(1),",", - & spin(2),",",spin(3),")" - end if + WRITE(*,*) " Sch radius/mass =",rsch,mass ENDIF - - + c Write to file at each radius test_myproc2: IF (CCTK_MyProc(cctkGH) == 0) THEN @@ -552,7 +527,7 @@ c Cannot use the conformal equation for ADM mass now DO idet = 1, ndet IF (detector(idet) > r_max) THEN - IF (openfiles == 1 .OR. cctk_iteration == it) THEN + IF (openfile == 1 .OR. cctk_iteration == it) THEN call CCTK_WARN(8,"Removing detectors outside coordinate range") c IF (iconv == 0) STOP END IF @@ -563,13 +538,10 @@ c IF (iconv == 0) STOP 404 CONTINUE IF (verbose == 1) THEN - IF (openfiles == 1) THEN - WRITE(*,*) - WRITE(*,*) "Extracting at detectors" - WRITE(*,*) "-----------------------" - WRITE(*,*) "Using ",ndet," detectors at" + IF (openfile == 1) THEN + WRITE(*,*) "Extracting at ",ndet," detectors " DO idet=1,ndet - WRITE(*,*) "r = ",detector(idet) + WRITE(*,*) " r = ",detector(idet) ENDDO END IF END IF @@ -583,41 +555,15 @@ c IF (iconv == 0) STOP END IF CALL D3_extract(cctkGH,do_ADMmass,do_momentum,do_spin, - & igrid,orig,myproc,Nt, - & Np,all_modes,lmode,mmode,x_1d,y_1d,z_1d, - & Dx,Dy,Dz,Psi_power,Psi,g00, - & gxx,gxy,gxz,gyy,gyz,gzz,kxx,kxy,kxz,kyy,kyz,kzz, - & radius,ADMmass,momentum,spin,mass,rsch, + & igrid,orig,myproc,Nt,Np,all_modes,lmode,mmode,x_1d,y_1d,z_1d, + & Dx,Dy,Dz,Psi_power,Psi,g00,gxx,gxy,gxz,gyy,gyz,gzz,kxx,kxy,kxz, + & kyy,kyz,kzz,radius,ADMmass,momentum,spin,mass,rsch, & Qodd,Qeven,temp3d,dtaudt) - test_verbose: - & IF (verbose == 1) THEN - - WRITE(*,*) + IF (verbose == 1) THEN WRITE(*,*) "Detector ",idet," ..." - WRITE(*,*) "Schwarzschild radius =",rsch - WRITE(*,*) "Schwarzschild mass =",mass - - loop_l4: DO il = lmin,lmax,lstep - - IF (all_modes == 0) THEN - mmin = mmode ; mmax = mmode - ELSE - mmin = 0 ; mmax = il - END IF - - loop_m4: DO im = mmin,mmax,mstep - - WRITE(*,*) filestr(1:nchar)//"/Qeven_"// - & CHAR(il+48)//CHAR(im+48),Qeven(:,il,im) - WRITE(*,*) filestr(1:nchar)//"/Qodd_"// - & CHAR(il+48)//CHAR(im+48),Qodd(:,il,im) - - END DO loop_m4 - - END DO loop_l4 - - END IF test_verbose + WRITE(*,*) " Sch radius/mass =",rsch,mass + END IF c Output to files @@ -664,85 +610,25 @@ c Output extracted radius and mass spinfile3 = filestr(1:nchar)//"/spin_z_"//"R"// & CHAR(idet+48)//timestring - IF (openfiles == 1) THEN - OPEN(UNIT= 77,FILE=rschfile,STATUS="unknown") - OPEN(UNIT= 78,FILE=massfile,STATUS="unknown") - IF (do_ADMmass(1) == 1) - & OPEN(UNIT= 79,FILE=ADMmassfile1,STATUS="unknown") - IF (do_ADMmass(2) == 1) - & OPEN(UNIT= 80,FILE=ADMmassfile2,STATUS="unknown") - IF (do_momentum == 1) THEN - OPEN(UNIT=781,FILE=momentumfile1,STATUS="unknown") - OPEN(UNIT=782,FILE=momentumfile2,STATUS="unknown") - OPEN(UNIT=783,FILE=momentumfile3,STATUS="unknown") - END IF - - IF (do_spin == 1) THEN - OPEN(UNIT=784,FILE=spinfile1,STATUS="unknown") - OPEN(UNIT=785,FILE=spinfile2,STATUS="unknown") - OPEN(UNIT=786,FILE=spinfile3,STATUS="unknown") - END IF - ELSE - OPEN(UNIT= 77,FILE=rschfile,STATUS="old", - & POSITION="append") - OPEN(UNIT= 78,FILE=massfile,STATUS="old", - & POSITION="append") - IF (do_ADMmass(1) == 1) - & OPEN(UNIT= 79,FILE=ADMmassfile1,STATUS="old", - & POSITION="append") - IF (do_ADMmass(2) == 1) - & OPEN(UNIT= 80,FILE=ADMmassfile2,STATUS="old", - & POSITION="append") - IF (do_momentum == 1) THEN - OPEN(UNIT=781,FILE=momentumfile1,STATUS="old", - & POSITION="append") - OPEN(UNIT=782,FILE=momentumfile2,STATUS="old", - & POSITION="append") - OPEN(UNIT=783,FILE=momentumfile3,STATUS="old", - & POSITION="append") - END IF - - IF (do_spin == 1)THEN - OPEN(UNIT=784,FILE=spinfile1,STATUS="old", - & POSITION="append") - OPEN(UNIT=785,FILE=spinfile2,STATUS="old", - & POSITION="append") - OPEN(UNIT=786,FILE=spinfile3,STATUS="old", - & POSITION="append") - END IF - - ENDIF - WRITE( 77,21) time,rsch - WRITE( 78,21) time,mass - IF (do_ADMmass(1) == 1) WRITE( 79,21) time,ADMmass(1) - IF (do_ADMmass(2) == 1) WRITE( 80,21) time,ADMmass(2) - IF (do_momentum == 1) THEN - WRITE(781,*) cctk_time,momentum(1) - WRITE(782,*) cctk_time,momentum(2) - WRITE(783,*) cctk_time,momentum(3) + call Extract_write(openfile,rschfile,time,rsch) + call Extract_write(openfile,massfile,time,mass) + IF (do_ADMmass(1) == 1) THEN + call Extract_write(openfile,ADMmassfile1,time,ADMmass(1)) END IF - IF (do_spin == 1) THEN - WRITE(784,*) cctk_time,spin(1) - WRITE(785,*) cctk_time,spin(2) - WRITE(786,*) cctk_time,spin(3) + IF (do_ADMmass(2) == 1) THEN + call Extract_write(openfile,ADMmassfile2,time,ADMmass(2)) END IF - - CLOSE( 77) - CLOSE( 78) - IF (do_ADMmass(1) == 1) CLOSE( 79) - IF (do_ADMmass(2) == 1) CLOSE( 80) - IF (do_momentum == 1) THEN - CLOSE(781) - CLOSE(782) - CLOSE(783) - END IF + call Extract_write(openfile,momentumfile1,time,momentum(1)) + call Extract_write(openfile,momentumfile2,time,momentum(2)) + call Extract_write(openfile,momentumfile3,time,momentum(3)) + END IF IF (do_spin == 1) THEN - CLOSE(782) - CLOSE(783) - CLOSE(784) + call Extract_write(openfile,spinfile1,cctk_time,spin(1)) + call Extract_write(openfile,spinfile2,cctk_time,spin(2)) + call Extract_write(openfile,spinfile3,cctk_time,spin(3)) END IF - + c Output gauge invariant variables @@ -756,9 +642,6 @@ c Output gauge invariant variables loop_m5: DO im = mmin,mmax,mstep - fn1 = 11 - fn2 = 12 - out1 = filestr(1:nchar)//"/Qeven_"//"R" & //CHAR(idet+48)//"_"//CHAR(il+48)// & CHAR(im+48)//timestring @@ -767,29 +650,11 @@ c Output gauge invariant variables & CHAR(im+48)//timestring c Write even parity waveforms - IF (openfiles == 1) THEN - OPEN(UNIT=fn1,FILE=out1,STATUS="unknown") - ELSE - OPEN(UNIT=fn1,FILE=out1,STATUS="old", - & POSITION="append") - ENDIF - - WRITE(fn1,20) time,Qeven(:,il,im) - - CLOSE(fn1) + call Extract_write(openfile,out1,time,Qeven(:,il,im)) c Only write odd parity waveforms if full grid IF (igrid == 0) THEN - IF (openfiles == 1) THEN - OPEN(UNIT=fn2,FILE=out2,STATUS="unknown") - ELSE - OPEN(UNIT=fn2,FILE=out2,STATUS="old", - & POSITION="append") - ENDIF - - WRITE(fn2,20) time,Qodd(:,il,im) - - CLOSE(fn2) + call Extract_write(openfile,out2,time,Qodd(:,il,im)) END IF END DO loop_m5 @@ -806,9 +671,50 @@ c Only write odd parity waveforms if full grid DEALLOCATE(Qodd,Qeven) - openfiles = 0 + openfile = 0 20 format(e24.15,e24.15,e24.15) 21 format(e24.15,e24.15) END SUBROUTINE Extract + + + SUBROUTINE Extract_open(openfile,filename,filehandle) + + implicit none + + DECLARE_CCTK_PARAMETERS + + integer filehandle,openfile + character*(*) filename + + if (openfile==1) then + OPEN(UNIT= filehandle,FILE=filename,STATUS="unknown") + else + OPEN(UNIT=filehandle,FILE=filename,STATUS="old",POSITION="append") + end if + + END SUBROUTINE Extract_open + + + SUBROUTINE Extract_write(openfile,filename,value1,value2) + + implicit none + + DECLARE_CCTK_PARAMETERS + + CCTK_REAL value1,value2 + integer openfile + character*(*) filename + + if (openfile==1) then + OPEN(UNIT= 444,FILE=filename,STATUS="unknown") + else + OPEN(UNIT=444,FILE=filename,STATUS="old",POSITION="append") + end if + + write(444,*) value1,value2 + + close(444) + + END SUBROUTINE Extract_write diff --git a/src/energy.f b/src/energy.f index 19ef5e2..2466e83 100644 --- a/src/energy.f +++ b/src/energy.f @@ -1,3 +1,4 @@ + PROGRAM energy c ----------------------------------------------------------------- @@ -17,7 +18,7 @@ c ----------------------------------------------------------------- CHARACTER*50 filename REAL*8 t(nmax),Qplus(nmax),dEdt(nmax),E,ddQ(nmax),Qnew(nmax) REAL*8 Dtnew,fac,time,t1,t2,dQdt_l,dQdt_r,dtmin,dt - + INTEGER :: ioerror WRITE(6,*) WRITE(6,*) "Starting PROGRAM energy ..." WRITE(6,*) "---------------------------" @@ -30,7 +31,8 @@ c ----------------------------------------------------------------- OPEN(UNIT=12,FILE=filename,STATUS="old") DO i=1,nmax - READ(12,*,ERR=101) t(i),Qplus(i) + READ(12,*,ERR=101, IOSTAT=ioerror) t(i),Qplus(i) + if(ioerror /= 0) goto 101 itotal = i END DO WRITE(6,*) "Didn't read all of data file" @@ -101,12 +103,12 @@ c Calculate energy flux dEdt(Nt) = fac*((Qnew(Nt)-Qnew(Nt-1))/Dtnew)**2 c Write results to file - OPEN(UNIT=101,FILE="dEdt",STATUS="unknown") + OPEN(UNIT=21,FILE="dEdt",STATUS="unknown") DO i = 1,Nt time = t1+DBLE(i-1)*Dtnew - WRITE(101,*) time,dEdt(i) + WRITE(21,*) time,dEdt(i) END DO - CLOSE(101) + CLOSE(21) c Use Simpsons rule to calculate integral IF ( MOD(Nt,2) /= 0) THEN @@ -202,4 +204,3 @@ c Write out total energy END - -- cgit v1.2.3