diff options
author | allen <allen@5301f0c2-dbc4-4cee-b2f5-8d7afba4d129> | 2000-01-11 09:11:58 +0000 |
---|---|---|
committer | allen <allen@5301f0c2-dbc4-4cee-b2f5-8d7afba4d129> | 2000-01-11 09:11:58 +0000 |
commit | 2c5a9898aade0ef52e6208bbb3fca110e66b742d (patch) | |
tree | 773c841581385c5de9fde6ab17dad0d96dcf9846 /src/Extract.F | |
parent | 9e60c939734aa6bdcdd6b9927d74043fed5b3b15 (diff) |
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
Diffstat (limited to 'src/Extract.F')
-rw-r--r-- | src/Extract.F | 266 |
1 files changed, 86 insertions, 180 deletions
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 |