c/*@@ c @file Extract.F c @date 28th October 1997 c @author Gab Allen c @desc Waveform extraction c c @enddesc c@@*/ #include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" c/*@@ c @routine Extract c @date 28th October 1997 c @author Gab Allen c @desc Entry point for waveform extraction routines c c @enddesc c @calls D3_extract c @calledby No idea c @history c c @endhistory c@@*/ SUBROUTINE Extract(CCTK_ARGUMENTS) USE D3_extract_int IMPLICIT NONE DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS c Non-Cactus input variables for D3_extract INTEGER :: & igrid,lmode,mmode,Psi_power=1 INTEGER :: & do_ADMmass(2) CCTK_REAL :: & orig(3),radius CCTK_REAL :: & x_1d(cctk_lsh(1)),y_1d(cctk_lsh(2)),z_1d(cctk_lsh(3)) c Output variables from D3_extract CCTK_REAL :: & dtaudt,ADMmass(2),mass,rsch,momentum(3),spin(3) CCTK_REAL,ALLOCATABLE :: & Qodd(:,:,:), Qeven(:,:,:) c Local variables INTEGER :: & ix,iy,iz,fn1,fn2,lmin,lmax,mmin,mmax,lstep,mstep, & il,im,it,ndet,idet,ioutput INTEGER,SAVE :: openfile = 1 INTEGER,PARAMETER :: & max_detectors = 9,number_timecoords = 2 CCTK_REAL,PARAMETER :: & two = 2.0D0 CCTK_REAL :: & time,r1,r2,dr,r_max,xmin,xmax,ymin,ymax,zmin,zmax CCTK_REAL :: & detector(max_detectors) CHARACTER*100 :: & out1,out2,massfile,rschfile,ADMmassfile1,ADMmassfile2, & momentumfile1,momentumfile2,momentumfile3,spinfile1,spinfile2,spinfile3 CHARACTER*3 :: & timestring LOGICAL :: use_proper,use_coordinate,do_output CCTK_REAL,DIMENSION(10),SAVE :: proper_time CCTK_REAL,SAVE :: last_time INTEGER, save :: open_file_level(10) = 0 INTEGER myproc,ierr CCTK_REAL dx,dy,dz character*80 infoline c ------------------------------------------------------------------ myproc = CCTK_MyProc(cctkGH) dx = cctk_delta_space(1) dy = cctk_delta_space(2) dz = cctk_delta_space(3) c ------------------------------------------------------------------ c c 0. Check to see if Extract should have been called c c ------------------------------------------------------------------ IF ((itout) .LE. 0) RETURN IF(MOD(int(cctk_iteration),int(itout)) .NE. 0) RETURN if (verbose == 1) then write(infoline,'(A24,G12.7)') & 'Calling Extract at time ',cctk_time call CCTK_INFO(infoline) end if c ------------------------------------------------------------------ c c 1. Initial Stuff c c ------------------------------------------------------------------ c Get the number of polar divisions IF (MOD(int(Nt),2) == 1 .OR. MOD(int(Np),2) == 1 .OR. & Nt < 0 .OR. Np < 0) THEN call CCTK_WARN(0,"Error in Nt or Np in Extract") END IF c Get the origin of spherical symmetry orig(1) = origin_x orig(2) = origin_y orig(3) = origin_z c Set the value of igrid IF (CCTK_EQUALS(domain,"octant")) THEN igrid = 1 ELSEIF (CCTK_EQUALS(domain,"bitant")) THEN igrid = 3 ELSE igrid = 0 END IF #ifdef THORN_CARTOON_2D IF (contains("cartoon_active","yes") .NE. 0) THEN igrid = 2 Np = 1 END IF #endif c Create 1D coordinate arrays do ix = 1, cctk_lsh(1) x_1d(ix) = x(ix,1,1) enddo do iy = 1, cctk_lsh(2) y_1d(iy) = y(1,iy,1) enddo do iz = 1, cctk_lsh(3) z_1d(iz) = z(1,1,iz) enddo c See if the ADM mass should be calculated and output IF (doADMmass == 1) THEN do_ADMmass(1) = 1 IF (use_conformal == 1) THEN do_ADMmass(2) = 1 ELSE do_ADMmass(2) = 0 ENDIF ELSE do_ADMmass(:) = 0 END IF c Set array containing gtt ... for the moment I am lazy and c just make it the lapse and do not include the shift! g00 = alp c What kind of timecoordinate do I want to use if (CCTK_EQUALS(timecoord,"both")) then use_proper = .true. use_coordinate = .true. if (cctk_iteration == 0) then proper_time = 0 last_time = cctk_time end if elseif (CCTK_EQUALS(timecoord,"proper")) then use_proper = .true. use_coordinate = .false. if (cctk_iteration == 0) then proper_time = 0 last_time = cctk_time end if elseif (CCTK_EQUALS(timecoord,"coord")) then use_proper = .false. use_coordinate = .true. else use_proper = .false. use_coordinate = .true. endif c ------------------------------------------------------------------ c c 2. Sort out which modes to extract c c ------------------------------------------------------------------ lmode = l_mode IF (all_modes == 1) THEN lmin = 2 lmax = lmode c Get m_mode in this case too, or the T3E has a junk value (PW) mmode = 0 ALLOCATE(Qodd(2,2:lmode,0:lmode),Qeven(2,2:lmode,0:lmode)) ELSE lmin = lmode lmax = lmode mmode = m_mode ALLOCATE(Qodd(2,2:lmode,0:mmode),Qeven(2,2:lmode,0:mmode)) ENDIF c Check do not have bad values for the modes IF (lmode < 2 .OR. mmode > lmode .OR. mmode < 0) THEN call CCTK_WARN(0,"Error in lmode,mmode in Extract") END IF IF (lmode > 9) THEN call CCTK_WARN(4,"Filenames will probably be crazy in Extract") END IF c If using an octant, do not do odd modes IF (igrid == 1 .OR. igrid == 2) THEN lstep = 2 ; mstep = 2 ELSE lstep = 1 ; mstep = 1 END IF c If we have cartoon, then jump past all the m-modes in loop IF (igrid == 2) mstep = lmode+1 c ------------------------------------------------------------------ c c 3. Find maximum radius of extraction on grid c c ------------------------------------------------------------------ call CCTK_CoordRange(ierr,cctkGH,xmin,xmax,-1,"x","cart3d") call CCTK_CoordRange(ierr,cctkGH,ymin,ymax,-1,"y","cart3d") call CCTK_CoordRange(ierr,cctkGH,zmin,zmax,-1,"z","cart3d") IF (igrid == 2) THEN r_max = MIN(xmax-two*Dx-orig(1), & zmax-two*Dz-orig(3)) ELSE IF (igrid == 1) THEN r_max = MIN(zmax-two*Dx-orig(1), & ymax-two*Dy-orig(2), & zmax-two*Dz-orig(3)) ELSE IF (igrid == 3) THEN r_max = MIN(xmax-two*Dx-orig(1), & ymax-two*Dy-orig(2), & zmax-two*Dz-orig(3), & DABS(xmin+two*Dx-orig(1)), & DABS(ymin+two*Dy-orig(2))) ELSE r_max = MIN(xmax-two*Dx-orig(1), & ymax-two*Dy-orig(2), & zmax-two*Dz-orig(3), & DABS(xmin+two*Dx-orig(1)), & DABS(ymin+two*Dy-orig(2)), & DABS(zmin+two*Dz-orig(3))) END IF c ------------------------------------------------------------------ c c 4. Extract Cauchy initial data for a linear wave equation c c ------------------------------------------------------------------ test_Cauchy: IF (Cauchy == 1) THEN it = Cauchy_timestep c check_timestep: IF (cctk_iteration == it) THEN check_timestep: IF (open_file_level(1) .eq. 0) THEN c Open output files write (*,*) "Open output file on level ", 1 open_file_level(1)=1; test_myproc1: IF (CCTK_MyProc(cctkGH) == 0 ) THEN out1 = "rsch_ini.rl" call Extract_open(cctkGH,0,out1,77) out1 = "mass_ini.rl" call Extract_open(cctkGH,0,out1,78) IF (do_ADMmass(1) == 1) THEN out1 = "ADMmass_ini.rl" call Extract_open(cctkGH,0,out1,79) ENDIF IF (do_ADMmass(2) == 1) THEN out1 = "ADMmassc_ini.rl" call Extract_open(cctkGH,0,out1,80) ENDIF IF (do_momentum == 1) THEN out1 = "momentum_x_ini.rl" call Extract_open(cctkGH,0,out1,781) out1 = "momentum_y_ini.rl" call Extract_open(cctkGH,0,out1,782) out1 = "momentum_z_ini.rl" call Extract_open(cctkGH,0,out1,783) ENDIF IF (do_spin == 1) THEN out1 = "spin_x_ini.rl" call Extract_open(cctkGH,0,out1,784) out1 = "spin_y_ini.rl" call Extract_open(cctkGH,0,out1,785) out1 = "spin_z_ini.rl" call Extract_open(cctkGH,0,out1,786) ENDIF loop_l1: DO il = lmin,lmax,lstep IF (all_modes == 0) THEN mmin = mmode ; mmax = mmode ELSE mmin = 0 ; mmax = il END IF loop_m1: DO im = mmin,mmax,lstep fn1 = il*10+im out1 = "Qeven_ini_"// & CHAR(il+48)//CHAR(im+48)//".rl" call Extract_open(cctkGH,0,out1,fn1) c Only print odd-parity if full grid IF (igrid == 0) THEN fn2 = 100+il*10+im out2 = "Qodd_ini_"// & CHAR(il+48)//CHAR(im+48)//".rl" call Extract_open(cctkGH,0,out2,fn2) END IF END DO loop_m1 END DO loop_l1 END IF test_myproc1 c Find range of extraction radii r1 = Cauchy_r1 r2 = r_max dr = Cauchy_dr c Do extraction at each radius IF (verbose == 1) THEN WRITE(*,*) "Extracting Cauchy initial data" WRITE(*,*) " r = ",r1," to ",r2," step ",dr END IF radius = r1 extract_at_each_radius: DO WHILE (radius < r2) CALL D3_extract(cctkGH,do_ADMmass,do_momentum,do_spin, & igrid,orig,myproc,interpolation_order,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) IF (verbose == 1) THEN WRITE(*,*) "Extracted at r =",radius WRITE(*,*) " Sch radius/mass =",rsch,mass ENDIF c Write to file at each radius test_myproc2: IF (CCTK_MyProc(cctkGH) == 0) THEN WRITE( 77,*) radius,rsch WRITE( 78,*) radius,mass IF (do_ADMmass(1) == 1) WRITE( 79,*) radius,ADMmass(1) IF (do_ADMmass(2) == 1) WRITE( 80,*) radius,ADMmass(2) IF (do_momentum == 1) then WRITE(781,*) radius,momentum(1) WRITE(782,*) radius,momentum(2) WRITE(783,*) radius,momentum(3) end if IF (do_spin == 1) THEN WRITE(784,*) radius,spin(1) WRITE(785,*) radius,spin(2) WRITE(786,*) radius,spin(3) END IF loop_l2: DO il = lmin,lmax,lstep IF (all_modes == 0) THEN mmin = mmode ; mmax = mmode ELSE mmin = 0 ; mmax = il END IF loop_m2: DO im = mmin,mmax,mstep fn1 = il*10+im WRITE(fn1,*) rsch,Qeven(:,il,im) c Only print odd-parity if full grid IF (igrid == 0) THEN fn2 = 100+il*10+im WRITE(fn2,*) rsch,Qeven(:,il,im) END IF END DO loop_m2 END DO loop_l2 END IF test_myproc2 radius = radius+dr END DO extract_at_each_radius c Now close all the files test_myproc3: IF (myproc == 0) THEN CLOSE( 77) ! Schwarzschild radius CLOSE( 78) ! Schwarzschild mass IF (do_ADMmass(1) == 1) CLOSE( 79) ! ADM mass IF (do_ADMmass(2) == 1) CLOSE( 80) ! ADM mass IF (do_momentum == 1) THEN CLOSE(781) ! momentum CLOSE(782) ! momentum CLOSE(783) ! momentum end if IF (do_spin == 1) THEN CLOSE(784) ! spin CLOSE(785) ! spin CLOSE(786) ! spin END IF loop_l3: DO il = lmin,lmax,lstep IF (all_modes == 0) THEN mmin = mmode ; mmax = mmode ELSE mmin = 0 ; mmax = il END IF loop_m3: DO im = mmin,mmax,mstep fn1 = il*10+im CLOSE(fn1) ! Qeven_ini_lm.dat IF (igrid == 0) THEN fn2 = 100+il*10+im CLOSE(fn2) ! Qodd_ini_lm.dat END IF END DO loop_m3 END DO loop_l3 END IF test_myproc3 END IF check_timestep END IF test_Cauchy c ------------------------------------------------------------------ c c 5. Extract waveforms at detectors c c ------------------------------------------------------------------ c Cannot use the conformal equation for ADM mass now do_ADMmass(2) = 0 ndet = num_detectors test_detectors: IF (ndet > 0) THEN IF (ndet > 9) THEN call CCTK_WARN(0,"Too many detectors in Extract") END IF detector(1) = detector1 detector(2) = detector2 detector(3) = detector3 detector(4) = detector4 detector(5) = detector5 detector(6) = detector6 detector(7) = detector7 detector(8) = detector8 detector(9) = detector9 DO idet = ndet+1, max_detectors detector(idet) = 0.0D0 END DO DO idet = 1, ndet IF (detector(idet) > r_max) THEN IF (openfile == 1 .OR. cctk_iteration == it) THEN call CCTK_WARN(8,"Removing detectors outside coordinate range") END IF ndet =idet-1 GOTO 404 ENDIF END DO 404 CONTINUE IF (verbose == 1) THEN IF (openfile == 1) THEN WRITE(*,*) "Extracting at ",ndet," detectors " DO idet=1,ndet WRITE(*,*) " r = ",detector(idet) ENDDO END IF END IF detector_loop: DO idet = 1,ndet radius = detector(idet) IF (verbose == 1) THEN WRITE(*,*) "Calling extract at radius ... ",radius END IF CALL D3_extract(cctkGH,do_ADMmass,do_momentum,do_spin, & igrid,orig,myproc,interpolation_order,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) IF (verbose == 1) THEN WRITE(*,*) "Detector ",idet," ..." WRITE(*,*) " Sch radius/mass =",rsch,mass END IF c Output to files test_myproc4: IF (myproc == 0) THEN DO ioutput=1,number_timecoords do_output = .false. timestring = "ERR" time = 0D0 IF (ioutput == 1 .AND. use_coordinate) THEN timestring = ".tl" do_output = .true. time = cctk_time ELSEIF (ioutput == 2 .AND. use_proper) THEN timestring = ".ul" do_output = .true. time = proper_time(idet) + (cctk_time-last_time)*dtaudt proper_time(idet) = time if (idet == ndet) last_time = cctk_time ENDIF IF (do_output) THEN c Output extracted radius and mass rschfile = "rsch_"//"R"// & CHAR(idet+48)//timestring massfile = "mass_"//"R"// & CHAR(idet+48)//timestring ADMmassfile1 = "ADMmass_"//"R"// & CHAR(idet+48)//timestring ADMmassfile2 = "ADMmassc_"//"R"// & CHAR(idet+48)//timestring momentumfile1 = "momentum_x_"//"R"// & CHAR(idet+48)//timestring momentumfile2 = "momentum_y_"//"R"// & CHAR(idet+48)//timestring momentumfile3 = "momentum_z_"//"R"// & CHAR(idet+48)//timestring spinfile1 = "spin_x_"//"R"// & CHAR(idet+48)//timestring spinfile2 = "spin_y_"//"R"// & CHAR(idet+48)//timestring spinfile3 = "spin_z_"//"R"// & CHAR(idet+48)//timestring call Extract_write(cctkGH,openfile,rschfile,time,rsch) call Extract_write(cctkGH,openfile,massfile,time,mass) IF (do_ADMmass(1) == 1) THEN call Extract_write(cctkGH,openfile,ADMmassfile1,time,ADMmass(1)) END IF IF (do_ADMmass(2) == 1) THEN call Extract_write(cctkGH,openfile,ADMmassfile2,time,ADMmass(2)) END IF IF (do_momentum == 1) THEN call Extract_write(cctkGH,openfile,momentumfile1,time,momentum(1)) call Extract_write(cctkGH,openfile,momentumfile2,time,momentum(2)) call Extract_write(cctkGH,openfile,momentumfile3,time,momentum(3)) END IF IF (do_spin == 1) THEN call Extract_write(cctkGH,openfile,spinfile1,cctk_time,spin(1)) call Extract_write(cctkGH,openfile,spinfile2,cctk_time,spin(2)) call Extract_write(cctkGH,openfile,spinfile3,cctk_time,spin(3)) END IF c Output gauge invariant variables loop_l5: DO il = lmin,lmax,lstep IF (all_modes == 0) THEN mmin = mmode ; mmax = mmode ELSE mmin = 0 ; mmax = il END IF loop_m5: DO im = mmin,mmax,mstep out1 = "Qeven_"//"R" & //CHAR(idet+48)//"_"//CHAR(il+48)// & CHAR(im+48)//timestring out2 = "Qodd_"//"R" & //CHAR(idet+48)//"_"//CHAR(il+48)// & CHAR(im+48)//timestring c Write even parity waveforms call Extract_write2(cctkGH,openfile,out1,time,Qeven(:,il,im)) c Only write odd parity waveforms if full grid IF (igrid == 0) THEN call Extract_write2(cctkGH,openfile,out2,time,Qodd(:,il,im)) END IF END DO loop_m5 END DO loop_l5 ENDIF END DO END IF test_myproc4 END DO detector_loop END IF test_detectors DEALLOCATE(Qodd,Qeven) openfile = 0 20 format(e24.15,e24.15,e24.15) 21 format(e24.15,e24.15) END SUBROUTINE Extract SUBROUTINE Extract_open(cctkGH,openfile,filename,filehandle) implicit none DECLARE_CCTK_PARAMETERS integer filehandle,openfile CCTK_INT nchar character*(*) filename CCTK_POINTER cctkGH character*200 fullname character*200 filestr c Get the output directory sorted out call CCTK_FortranString(nchar,outdir,filestr) fullname = filestr(1:nchar)//"/"//filename if (openfile==1) then OPEN(UNIT= filehandle,FILE=fullname,STATUS="unknown") write(filehandle,101) char(34),filename 101 format(a1,a14) call Extract_Advertise(cctkGH,fullname) else OPEN(UNIT=filehandle,FILE=fullname,STATUS="old", . POSITION="append") end if END SUBROUTINE Extract_open SUBROUTINE Extract_write(cctkGH, openfile,filename,value1,value2) implicit none DECLARE_CCTK_PARAMETERS CCTK_POINTER cctkGH CCTK_REAL value1,value2 integer openfile CCTK_INT nchar character*(*) filename character*200 fullname character*200 filestr character*80 infoline c Get the output directory sorted out call CCTK_FortranString(nchar,outdir,filestr) fullname = filestr(1:nchar)//"/"//filename if (openfile==1) then OPEN(UNIT= 444,FILE=fullname,STATUS="unknown") write(444,101) char(34),filename 101 format(a1,a14) call Extract_Advertise(cctkGH,fullname) else OPEN(UNIT=444,FILE=fullname,STATUS="old",POSITION="append") end if write(444,*) value1,value2 close(444) END SUBROUTINE Extract_write SUBROUTINE Extract_write2(cctkGH, openfile,filename,value1,value2) implicit none DECLARE_CCTK_PARAMETERS CCTK_POINTER cctkGH CCTK_REAL value1 CCTK_REAL value2(2) integer openfile CCTK_INT nchar character*(*) filename character*200 fullname character*200 filestr character*80 infoline c Get the output directory sorted out call CCTK_FortranString(nchar,outdir,filestr) fullname = filestr(1:nchar)//"/"//filename if (openfile==1) then OPEN(UNIT= 444,FILE=fullname,STATUS="unknown") write(444,101) char(34),filename 101 format(a1,a14) call Extract_Advertise(cctkGH,fullname) else OPEN(UNIT=444,FILE=fullname,STATUS="old",POSITION="append") end if write(444,*) value1,value2(1),value2(2) close(444) END SUBROUTINE Extract_write2