From df8bb95c19fd2476e762daa0f4880fea475e3b34 Mon Sep 17 00:00:00 2001 From: pollney Date: Fri, 13 Dec 2002 15:39:01 +0000 Subject: Modified output so that data is appended after checkpoint, rather than deleting any existing files. * Starting a run from scratch will still wipe out any files that already exist in the output directory, as is usual for cactus 1d data. * If you do a horizon find on the timestep when you are also checkpointing, then you'll end up with a duplicate point in the output, as the horizon is found again at restart. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinder/trunk@320 89daf98e-ef62-4674-b946-b8ff9de2216c --- schedule.ccl | 5 + src/AHFinder_dat.F | 2 +- src/AHFinder_dis.F | 19 ++-- src/AHFinder_flow.F | 2 +- src/AHFinder_gau.F | 47 ++------- src/AHFinder_min.F | 2 +- src/AHFinder_output.F | 270 +++++++++++++++++++++++++++++++++++++++++--------- src/AHFinder_pow.F | 4 +- 8 files changed, 248 insertions(+), 103 deletions(-) diff --git a/schedule.ccl b/schedule.ccl index 4cbf7e7..1e17003 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -20,6 +20,11 @@ schedule AHFinder_ParamCheck at PARAMCHECK LANG: Fortran } "Check for physical or conformal metric" +schedule AHFinder_InitOutput at CCTK_INITIAL +{ + LANG: Fortran +} "Create output files, write headers" + if (ahf_active) { STORAGE: hole1_bounds,hole2_bounds,hole3_bounds,ahf_centroid diff --git a/src/AHFinder_dat.F b/src/AHFinder_dat.F index a7fc589..718d16e 100644 --- a/src/AHFinder_dat.F +++ b/src/AHFinder_dat.F @@ -59,7 +59,7 @@ CCTK_REAL, allocatable, dimension (:,:) :: hflowc,cflowc,nflowc CCTK_REAL, allocatable, dimension (:,:) :: hflows,cflows,nflows - character*200 filestr + character(len=200) :: filestr data ahf_ncall / 0 / data status_old / .false. / diff --git a/src/AHFinder_dis.F b/src/AHFinder_dis.F index 894271c..c18b0f8 100644 --- a/src/AHFinder_dis.F +++ b/src/AHFinder_dis.F @@ -40,7 +40,8 @@ CCTK_REAL, allocatable, dimension(:) :: xa,ya,za CCTK_REAL, allocatable, dimension(:) :: txx,tyy,tzz,txy,txz,tyz - character*200 disf + character(len=200) :: disf + logical disf_exists ! ************************** @@ -571,22 +572,16 @@ 100 continue - disf = filestr(1:nfile) // "/ahf_d12" // ".tl" - if (myproc.eq.0) then - - if (ahf_ncall.eq.1) then - open(1,file=disf,form='formatted',status='replace') - write(1,"(A18)") '" Horizon distance' + inquire(file=disf, exist=disf_exists) + if (disf_exists) then + open(1,file=disf,form='formatted',status='old', position='append') else - open(1,file=disf,form='formatted',status='old', - . position='append') + open(1,file=disf,form='formatted',status='replace') + write(1,"(A18)") '" Horizon distance' end if - write(1,"(2ES14.6)") cctk_time,d12 - close(1) - end if diff --git a/src/AHFinder_flow.F b/src/AHFinder_flow.F index 810019b..b5e4300 100644 --- a/src/AHFinder_flow.F +++ b/src/AHFinder_flow.F @@ -38,7 +38,7 @@ CCTK_REAL zero,one CCTK_REAL dd,aux - character*200 logf + character(len=200) :: logf ! Description of variables: ! diff --git a/src/AHFinder_gau.F b/src/AHFinder_gau.F index 28ab392..d66c214 100644 --- a/src/AHFinder_gau.F +++ b/src/AHFinder_gau.F @@ -66,6 +66,7 @@ CCTK_INT :: use_att, apower, nchars logical :: use_rot_att logical :: str_comp + logical :: gaussf_exists external str_comp character(len=32) :: tmpstr @@ -1143,10 +1144,9 @@ gaussf = filestr(1:nfile)//"/ahf.gauss" end if -! On first call replace file. + inquire(file=gaussf, exist=gaussf_exists) if (firstgau) then - if (find3) then if (mfind.eq.3) then firstgau = .false. @@ -1154,49 +1154,20 @@ else firstgau = .false. end if + end if - open(1,file=gaussf,form='formatted', - . status='replace') - - write(1,"(A20)") '# GAUSSIAN CURVATURE' - write(1,"(A1)") '#' - write(1,"(A35)") '# The data is written in a loop as:' - write(1,"(A1)") '#' - write(1,"(A17)") '# do i=1,ntheta' - write(1,"(A18)") '# do j=1,nphi' - write(1,"(A27)") '# write gaussian(i,j)' - write(1,"(A11)") '# end do' - write(1,"(A8)") '# end do' - write(1,"(A1)") '#' - write(1,"(A40)") '# theta and phi are subdivided uniformly' - write(1,"(A26)") '# according to symmetries:' - write(1,"(A1)") '#' - write(1,"(A36)") '# phi=[0,2 pi] (refx=refy=.false.)' - write(1,"(A44)") '# phi=[0,pi] (refx=.false., refy=.true.)' - write(1,"(A35)") '# phi=[0,pi/2] (refx=refy=.true.)' - write(1,"(A1)") '#' - write(1,"(A31)") '# theta=[0,pi] (refz=.false.)' - write(1,"(A30)") '# theta=[0,pi/2] (refz=.true.)' - write(1,"(A1)") '#' - - write(1,"(A9,L1)") '# refx = ',refx - write(1,"(A9,L1)") '# refy = ',refy - write(1,"(A9,L1)") '# refz = ',refz - write(1,"(A1)") '#' - write(1,"(A11,I7)") '# ntheta = ',(ntheta+1) - write(1,"(A11,I7)") '# nphi = ',(nphi+1) +! On first call replace file. - write(1,*) - write(1,"(A11,I4)") '# Time step',cctk_iteration - write(1,"(A6,ES11.3)") '# Time',cctk_time - write(1,"(A6,I4)") '# Call',ahf_ncall + if (.not.gaussf_exists) then + + call AHFinder_WriteGaussHeader(CCTK_PASS_FTOF, gaussf) ! On subsequent calls, append to file. else - open(1,file=gaussf,form='formatted', - . status='old',position='append') + open(1,file=gaussf,form='formatted', status='old', + . position='append') write(1,*) write(1,"(A11,I4)") '# Time step',cctk_iteration diff --git a/src/AHFinder_min.F b/src/AHFinder_min.F index b458d85..1c4a044 100644 --- a/src/AHFinder_min.F +++ b/src/AHFinder_min.F @@ -32,7 +32,7 @@ CCTK_REAL, dimension (NN) :: P CCTK_REAL, dimension (NN,NN) :: XI - character*200 logf + character(len=200) :: logf ! Description of variables: ! diff --git a/src/AHFinder_output.F b/src/AHFinder_output.F index 3f5b339..c9f4f8b 100644 --- a/src/AHFinder_output.F +++ b/src/AHFinder_output.F @@ -11,6 +11,167 @@ #include "cctk_Parameters.h" #include "cctk_Arguments.h" + subroutine AHFinder_InitOutput(CCTK_ARGUMENTS) + use AHFinder_dat + implicit none + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + character(len=3) :: ext + character(len=200) :: almf,areaf,massf,radf,circeqf,merip1f,merip2f + character(len=200) :: asymxf,asymyf,asymzf + character(len=200) :: disf + character(len=200) :: gaussf, gauss0f, gauss1f, gauss2f + + if (myproc.ne.0) then + return + end if + +! **************************** +! *** BUILD FILE NAMES *** +! **************************** + + call CCTK_FortranString(nfile,out_dir,filestr) + +c some compilers cannot trim an empty string, so we add at least one char + if (find3) then + ext = "_" // char(ichar('0') + mfind - 1) // "." + else + ext = "." + end if + + almf = filestr(1:nfile) // "/ahf_coeff" // trim(ext) // "alm" + areaf = filestr(1:nfile) // "/ahf_area" // trim(ext) // "tl" + massf = filestr(1:nfile) // "/ahf_mass" // trim(ext) // "tl" + radf = filestr(1:nfile) // "/ahf_rad" // trim(ext) // "tl" + + circeqf = filestr(1:nfile) // "/ahf_circ_eq" // trim(ext) // "tl" + merip1f = filestr(1:nfile) // "/ahf_meri_p1" // trim(ext) // "tl" + merip2f = filestr(1:nfile) // "/ahf_meri_p2" // trim(ext) // "tl" + + asymxf = filestr(1:nfile) // "/ahf_asymx" // trim(ext) // "tl" + asymyf = filestr(1:nfile) // "/ahf_asymy" // trim(ext) // "tl" + asymzf = filestr(1:nfile) // "/ahf_asymz" // trim(ext) // "tl" + + disf = filestr(1:nfile) // "/ahf_d12" // trim(ext) // "tl" + + gaussf = filestr(1:nfile) // "/ahf.gauss" + gauss0f = filestr(1:nfile) // "/ahf_0.gauss" + gauss1f = filestr(1:nfile) // "/ahf_1.gauss" + gauss2f = filestr(1:nfile) // "/ahf_2.gauss" + + + open(1,file=massf,form='formatted',status='replace') + write(1,"(A14)") '" Horizon mass' + close(1) + + open(1,file=radf,form='formatted',status='replace') + write(1,"(A26)") '" Horizon radius' + close(1) + + open(1,file=circeqf,form='formatted',status='replace') + write(1,"(A26)") '" Equatorial circumference' + close(1) + + open(1,file=merip1f,form='formatted',status='replace') + write(1,"(A27)") '" Length of meridian, phi=0' + close(1) + + open(1,file=merip1f,form='formatted',status='replace') + write(1,"(A27)") '" Length of meridian, phi=0' + close(1) + + open(1,file=merip2f,form='formatted',status='replace') + write(1,"(A30)") '" Length of meridian, phi=pi/2' + close(1) + + open(1,file=areaf,form='formatted',status='replace') + write(1,"(A14)") '" Horizon area' + close(1) + + open(1,file=asymxf,form='formatted',status='replace') + write(1,"(A29)") '" Asymmetries on x reflection' + close(1) + + open(1,file=asymyf,form='formatted',status='replace') + write(1,"(A29)") '" Asymmetries on y reflection' + close(1) + + open(1,file=asymzf,form='formatted',status='replace') + write(1,"(A29)") '" Asymmetries on z reflection' + close(1) + + open(1,file=almf,form='formatted',status='replace') + write(1,"(A21)") '# Radial coefficients' + write(1,"(A1)") '#' + close(1) + + open(1,file=disf,form='formatted',status='replace') + write(1,"(A18)") '" Horizon distance' + close(1) + + if (find3) then + call AHFinder_WriteGaussHeader(CCTK_PASS_FTOF, gauss0f) + call AHFinder_WriteGaussHeader(CCTK_PASS_FTOF, gauss1f) + call AHFinder_WriteGaussHeader(CCTK_PASS_FTOF, gauss2f) + else + call AHFinder_WriteGaussHeader(CCTK_PASS_FTOF, gaussf) + end if + + end subroutine AHFinder_InitOutput + + + + + subroutine AHFinder_WriteGaussHeader(CCTK_ARGUMENTS, filename) + use AHFinder_dat + implicit none + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + character(len=200) :: filename + + open(1,file=filename,form='formatted', status='replace') + + write(1,"(A20)") '# GAUSSIAN CURVATURE' + write(1,"(A1)") '#' + write(1,"(A35)") '# The data is written in a loop as:' + write(1,"(A1)") '#' + write(1,"(A17)") '# do i=1,ntheta' + write(1,"(A18)") '# do j=1,nphi' + write(1,"(A27)") '# write gaussian(i,j)' + write(1,"(A11)") '# end do' + write(1,"(A8)") '# end do' + write(1,"(A1)") '#' + write(1,"(A40)") '# theta and phi are subdivided uniformly' + write(1,"(A26)") '# according to symmetries:' + write(1,"(A1)") '#' + write(1,"(A36)") '# phi=[0,2 pi] (refx=refy=.false.)' + write(1,"(A44)") '# phi=[0,pi] (refx=.false., refy=.true.)' + write(1,"(A35)") '# phi=[0,pi/2] (refx=refy=.true.)' + write(1,"(A1)") '#' + write(1,"(A31)") '# theta=[0,pi] (refz=.false.)' + write(1,"(A30)") '# theta=[0,pi/2] (refz=.true.)' + write(1,"(A1)") '#' + + write(1,"(A9,L1)") '# refx = ',refx + write(1,"(A9,L1)") '# refy = ',refy + write(1,"(A9,L1)") '# refz = ',refz + write(1,"(A1)") '#' + write(1,"(A11,I7)") '# ntheta = ',(ntheta+1) + write(1,"(A11,I7)") '# nphi = ',(nphi+1) + + write(1,*) + write(1,"(A11,I4)") '# Time step',cctk_iteration + write(1,"(A6,ES11.3)") '# Time',cctk_time + write(1,"(A6,I4)") '# Call',ahf_ncall + + close(1) + + end subroutine AHFinder_WriteGaussHeader + + + subroutine AHFinder_Output(CCTK_ARGUMENTS, report, status, horizon, mtype, intarea_h) @@ -23,13 +184,16 @@ integer mtype logical report, status, horizon + logical massf_exists, radf_exists, circeqf_exists + logical merip1f_exists, merip2f_exists, areaf_exists + logical asymxf_exists, asymyf_exists, asymzf_exists, almf_exists CCTK_REAL intarea_h integer l, m, ierror - character*3 ext - character*200 almf,areaf,massf,radf,circeqf,merip1f,merip2f - character*200 asymxf,asymyf,asymzf - character*200 options + character(len=3) :: ext + character(len=200) :: almf,areaf,massf,radf,circeqf,merip1f,merip2f + character(len=200) :: asymxf,asymyf,asymzf + character(len=200) :: options ! **************************** @@ -97,12 +261,13 @@ c some compilers cannot trim an empty string, so we add at least one char if (myproc.eq.0) then ! Mass. - - if (ahf_ncall.eq.1) then - open(1,file=massf,form='formatted',status='replace') - write(1,"(A14)") '" Horizon mass' + + inquire(file=massf, exist=massf_exists) + if (massf_exists) then + open(1,file=massf,form='formatted',status='old',position='append') else - open(1,file=massf,form='formatted',status='old',position='append') + open(1,file=massf,form='formatted',status='replace') + write(1,"(A14)") '" Horizon mass' end if if (status.and.report) then @@ -116,11 +281,12 @@ c some compilers cannot trim an empty string, so we add at least one char ! Radius. - if (ahf_ncall.eq.1) then - open(1,file=radf,form='formatted',status='replace') - write(1,"(A16)") '" Horizon radius' + inquire(file=radf, exist=radf_exists) + if (radf_exists) then + open(1,file=radf,form='formatted',status='old',position='append') else - open(1,file=radf,form='formatted',status='old',position='append') + open(1,file=radf,form='formatted',status='replace') + write(1,"(A26)") '" Horizon radius' end if if (status.and.report) then @@ -134,11 +300,12 @@ c some compilers cannot trim an empty string, so we add at least one char ! Equatorial circumference. - if (ahf_ncall.eq.1) then - open(1,file=circeqf,form='formatted',status='replace') - write(1,"(A26)") '" Equatorial circumference' + inquire(file=circeqf, exist=circeqf_exists) + if (circeqf_exists) then + open(1,file=circeqf,form='formatted',status='old',position='append') else - open(1,file=circeqf,form='formatted',status='old',position='append') + open(1,file=circeqf,form='formatted',status='replace') + write(1,"(A26)") '" Equatorial circumference' end if if (status.and.report) then @@ -152,11 +319,12 @@ c some compilers cannot trim an empty string, so we add at least one char ! Meridians. - if (ahf_ncall.eq.1) then - open(1,file=merip1f,form='formatted',status='replace') - write(1,"(A27)") '" Length of meridian, phi=0' + inquire(file=merip1f, exist=merip1f_exists) + if (merip1f_exists) then + open(1,file=merip1f,form='formatted',status='old',position='append') else - open(1,file=merip1f,form='formatted',status='old',position='append') + open(1,file=merip1f,form='formatted',status='replace') + write(1,"(A27)") '" Length of meridian, phi=0' end if if (status.and.report) then @@ -168,11 +336,12 @@ c some compilers cannot trim an empty string, so we add at least one char close(1) - if (ahf_ncall.eq.1) then - open(1,file=merip2f,form='formatted',status='replace') - write(1,"(A30)") '" Length of meridian, phi=pi/2' + inquire(file=merip2f, exist=merip2f_exists) + if (merip2f_exists) then + open(1,file=merip2f,form='formatted',status='old', position='append') else - open(1,file=merip2f,form='formatted',status='old',position='append') + open(1,file=merip2f,form='formatted',status='replace') + write(1,"(A30)") '" Length of meridian, phi=pi/2' end if if (status.and.report) then @@ -185,12 +354,13 @@ c some compilers cannot trim an empty string, so we add at least one char close(1) ! Area. - - if (ahf_ncall.eq.1) then - open(1,file=areaf,form='formatted',status='replace') - write(1,"(A14)") '" Horizon area' + + inquire(file=areaf, exist=areaf_exists) + if (areaf_exists) then + open(1,file=areaf,form='formatted',status='old',position='append') else - open(1,file=areaf,form='formatted',status='old',position='append') + open(1,file=areaf,form='formatted',status='replace') + write(1,"(A14)") '" Horizon area' end if if (status.and.report) then @@ -206,11 +376,12 @@ c some compilers cannot trim an empty string, so we add at least one char ! coefficients cs with even m indicate asymmetries ! on refelction x<->-x. - if (ahf_ncall.eq.1) then - open(1,file=asymxf,form='formatted',status='replace') - write(1,"(A29)") '" Asymmetries on x reflection' + inquire(file=asymxf, exist=asymxf_exists) + if (asymxf_exists) then + open(1,file=asymxf,form='formatted',status='old', position='append') else - open(1,file=asymxf,form='formatted',status='old',position='append') + open(1,file=asymxf,form='formatted',status='replace') + write(1,"(A29)") '" Asymmetries on x reflection' end if out_asymx = 0.0D0 @@ -233,11 +404,12 @@ c some compilers cannot trim an empty string, so we add at least one char ! Asymmetries in y. Any cs coefficient different from ! zero indicates asymmetries on reflection y<->-y. - if (ahf_ncall.eq.1) then - open(1,file=asymyf,form='formatted',status='replace') - write(1,"(A29)") '" Asymmetries on y reflection' + inquire(file=asymyf, exist=asymyf_exists) + if (asymyf_exists) then + open(1,file=asymyf,form='formatted',status='old', position='append') else - open(1,file=asymyf,form='formatted',status='old',position='append') + open(1,file=asymyf,form='formatted',status='replace') + write(1,"(A29)") '" Asymmetries on y reflection' end if out_asymy = 0.0D0 @@ -257,11 +429,12 @@ c some compilers cannot trim an empty string, so we add at least one char ! or any coefficient cc and cs with odd (l-m) indicate ! asymmetries on reflection z<->-z. - if (ahf_ncall.eq.1) then - open(1,file=asymzf,form='formatted',status='replace') - write(1,"(A29)") '" Asymmetries on z reflection' + inquire(file=asymzf, exist=asymzf_exists) + if (asymzf_exists) then + open(1,file=asymzf,form='formatted',status='old',position='append') else - open(1,file=asymzf,form='formatted',status='old',position='append') + open(1,file=asymzf,form='formatted',status='replace') + write(1,"(A29)") '" Asymmetries on z reflection' end if out_asymz = 0.0D0 @@ -282,13 +455,14 @@ c some compilers cannot trim an empty string, so we add at least one char ! Expansion coefficients. - if (ahf_ncall.eq.1) then - open(1,file=almf,form='formatted',status='replace') - write(1,"(A21)") '# Radial coefficients' - write(1,"(A1)") '#' + inquire(file=almf, exist=almf_exists) + if (almf_exists) then + open(1,file=almf,form='formatted',status='old',position='append') + write(1,*) else - open(1,file=almf,form='formatted',status='old',position='append') - write(1,*) + open(1,file=almf,form='formatted',status='replace') + write(1,"(A21)") '# Radial coefficients' + write(1,"(A1)") '#' end if if (status.and.report) then diff --git a/src/AHFinder_pow.F b/src/AHFinder_pow.F index fe53079..a928f76 100644 --- a/src/AHFinder_pow.F +++ b/src/AHFinder_pow.F @@ -77,7 +77,7 @@ CCTK_REAL, dimension(1:N) :: P,PT,PTT,XIT CCTK_REAL, dimension(1:N,1:N) :: XI - character*200 logf + character(len=200) :: logf parameter (ZEPS=1.0D-10) @@ -691,7 +691,7 @@ CCTK_REAL, dimension(1:N) :: P - character*200 logf + character(len=200) :: logf save iter,logf -- cgit v1.2.3