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

4 files changed, 102 insertions, 198 deletions

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 
 
 
-