diff options
author | miguel <miguel@89daf98e-ef62-4674-b946-b8ff9de2216c> | 1999-10-29 16:14:09 +0000 |
---|---|---|
committer | miguel <miguel@89daf98e-ef62-4674-b946-b8ff9de2216c> | 1999-10-29 16:14:09 +0000 |
commit | b5c4ea91e453f1587044f3380824d35e663e1db0 (patch) | |
tree | c92ef48b92b18bd808941ac066a1be3762020532 /src/AHFinder.F | |
parent | 22d15bf63112317c25624309ea2d75a3dad10740 (diff) |
Some cleaning up. Gerd: Sorry, but c0(0) was initialized inside AHfinder_flow,
so there must be another problem.
Miguel
git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinder/trunk@21 89daf98e-ef62-4674-b946-b8ff9de2216c
Diffstat (limited to 'src/AHFinder.F')
-rw-r--r-- | src/AHFinder.F | 204 |
1 files changed, 90 insertions, 114 deletions
diff --git a/src/AHFinder.F b/src/AHFinder.F index 679715e..5dba527 100644 --- a/src/AHFinder.F +++ b/src/AHFinder.F @@ -50,8 +50,6 @@ c Note thadmcevol_onpunc_0.250_128_neu2_drop/at including cactus.h will also inc #include "CactusEinstein/Einstein/src/Einstein.h" - - ! Description of variables: ! ! status Surface found? @@ -101,6 +99,7 @@ c Note thadmcevol_onpunc_0.250_128_neu2_drop/at including cactus.h will also inc dy = cctk_delta_space(2) dz = cctk_delta_space(3) + ! ******************************** ! *** INITIALIZE {status} *** ! ******************************** @@ -154,19 +153,12 @@ c Note thadmcevol_onpunc_0.250_128_neu2_drop/at including cactus.h will also inc if (mod(cctk_iteration-ahfafter,ahfso).ne.0) return -! BoxInBox may turn of AHF this way during runtime (BB) +! Sorry for this hack, but the call order for BoxInBox +! is screwy in cactus 3.2 -! if (contains("ahf_active","no").ne.0) then - if (ahf_active.eq.0) then - return - endif - -! Sorry for this hack, but the call order for BoxInBox is screwy in -! cactus 3.2 #ifdef THORN_BOXINBOX -! if ((contains("boxinbox","yes").ne.0).and.(ahf_ncall.eq.0)) then if ((boxinbox.ne.0).and.(ahf_ncall.eq.0)) then - return; + return endif #endif @@ -214,15 +206,12 @@ c Note thadmcevol_onpunc_0.250_128_neu2_drop/at including cactus.h will also inc ! *** FIND OUT IF WE WANT TO FIND TRAPPED SURFACE *** ! ******************************************************* -! find_trapped_surface = contains("ahf_trapped_surface","yes").ne.0 if (ahf_trapped_surface.ne.0) then find_trapped_surface = .true. else find_trapped_surface = .false. end if -! trapped_surface_delta = getr8("trapped_surface_delta") - ! ******************************************** ! *** GET THE NAME OF OUTPUT DIRECTORY *** @@ -230,70 +219,64 @@ c Note thadmcevol_onpunc_0.250_128_neu2_drop/at including cactus.h will also inc call CCTK_FortranString(nfile,outdir,filestr) -cGAB : should there still be an if statement here - - if (.true.) then - - if (find3) then + if (find3) then - if (mfind.eq.0) then + if (mfind.eq.0) then - logf = filestr(1:nfile)//"ahf_logfile_0" - almf = filestr(1:nfile)//"ahf_coeff_0.alm" + logf = filestr(1:nfile)//"ahf_logfile_0" + almf = filestr(1:nfile)//"ahf_coeff_0.alm" - areaf = filestr(1:nfile)//"ahf_area_0.tl" - massf = filestr(1:nfile)//"ahf_mass_0.tl" + areaf = filestr(1:nfile)//"ahf_area_0.tl" + massf = filestr(1:nfile)//"ahf_mass_0.tl" - circeqf = filestr(1:nfile)//"ahf_circ_eq_0.tl" - merip1f = filestr(1:nfile)//"ahf_meri_p1_0.tl" - merip2f = filestr(1:nfile)//"ahf_meri_p2_0.tl" + circeqf = filestr(1:nfile)//"ahf_circ_eq_0.tl" + merip1f = filestr(1:nfile)//"ahf_meri_p1_0.tl" + merip2f = filestr(1:nfile)//"ahf_meri_p2_0.tl" - else if (mfind.eq.1) then + else if (mfind.eq.1) then - logf = filestr(1:nfile)//"ahf_logfile_1" - almf = filestr(1:nfile)//"ahf_coeff_1.alm" + logf = filestr(1:nfile)//"ahf_logfile_1" + almf = filestr(1:nfile)//"ahf_coeff_1.alm" - areaf = filestr(1:nfile)//"ahf_area_1.tl" - massf = filestr(1:nfile)//"ahf_mass_1.tl" + areaf = filestr(1:nfile)//"ahf_area_1.tl" + massf = filestr(1:nfile)//"ahf_mass_1.tl" - circeqf = filestr(1:nfile)//"ahf_circ_eq_1.tl" - merip1f = filestr(1:nfile)//"ahf_meri_p1_1.tl" - merip2f = filestr(1:nfile)//"ahf_meri_p2_1.tl" - - else + circeqf = filestr(1:nfile)//"ahf_circ_eq_1.tl" + merip1f = filestr(1:nfile)//"ahf_meri_p1_1.tl" + merip2f = filestr(1:nfile)//"ahf_meri_p2_1.tl" - logf = filestr(1:nfile)//"ahf_logfile_2" - almf = filestr(1:nfile)//"ahf_coeff_2.alm" + else - areaf = filestr(1:nfile)//"ahf_area_2.tl" - massf = filestr(1:nfile)//"ahf_mass_2.tl" + logf = filestr(1:nfile)//"ahf_logfile_2" + almf = filestr(1:nfile)//"ahf_coeff_2.alm" - circeqf = filestr(1:nfile)//"ahf_circ_eq_2.tl" - merip1f = filestr(1:nfile)//"ahf_meri_p1_2.tl" - merip2f = filestr(1:nfile)//"ahf_meri_p2_2.tl" + areaf = filestr(1:nfile)//"ahf_area_2.tl" + massf = filestr(1:nfile)//"ahf_mass_2.tl" - end if + circeqf = filestr(1:nfile)//"ahf_circ_eq_2.tl" + merip1f = filestr(1:nfile)//"ahf_meri_p1_2.tl" + merip2f = filestr(1:nfile)//"ahf_meri_p2_2.tl" - else + end if - logf = filestr(1:nfile)//"ahf_logfile" - almf = filestr(1:nfile)//"ahf_coeff.alm" + else + + logf = filestr(1:nfile)//"ahf_logfile" + almf = filestr(1:nfile)//"ahf_coeff.alm" - areaf = filestr(1:nfile)//"ahf_area.tl" - massf = filestr(1:nfile)//"ahf_mass.tl" + areaf = filestr(1:nfile)//"ahf_area.tl" + massf = filestr(1:nfile)//"ahf_mass.tl" - circeqf = filestr(1:nfile)//"ahf_circ_eq.tl" - merip1f = filestr(1:nfile)//"ahf_meri_p1.tl" - merip2f = filestr(1:nfile)//"ahf_meri_p2.tl" + circeqf = filestr(1:nfile)//"ahf_circ_eq.tl" + merip1f = filestr(1:nfile)//"ahf_meri_p1.tl" + merip2f = filestr(1:nfile)//"ahf_meri_p2.tl" - end if end if + ! ************************************************** ! *** TRANSFORM CONFORMAL TO PHYSICAL METRIC *** ! ************************************************** - -! call ConfToPhys(nx,ny,nz,psi,gxx,gxy,gxz,gyy,gyz,gzz) call CCTK_WARN(4,"Converting metric: conformal -> physical") @@ -306,54 +289,51 @@ cGAB : should there still be an if statement here conformal_state = NOCONFORMAL_METRIC + ! *************************************** ! *** FIND OUT THE TYPE OF OUTPUT *** ! *************************************** -! logfile = contains("ahf_logfile","yes").ne.0 if (ahf_logfile.ne.0) then logfile = .true. else logfile = .false. end if -! verbose = contains("ahf_verbose","yes").ne.0 if (ahf_verbose.ne.0) then verbose = .true. else verbose = .false. end if - -! veryver = contains("ahf_veryverbose","yes").ne.0 + if (ahf_veryverbose.ne.0) then veryver = .true. else veryver = .false. end if -! guessverbose = contains("ahf_guessverbose","yes").ne.0 if (ahf_guessverbose.ne.0) then guessverbose = .true. else guessverbose = .false. end if - if (logfile.and.(myproc.eq.0)) then if (logfile) then if (myproc.eq.0) then open(1,file=logf,form='formatted',status='replace') write(1,*) - write(1,*) 'LOG FILE FOR MINIMIZATION APPARENT HORIZON FINDER' + write(1,*) 'LOG FILE FOR APPARENT HORIZON FINDER' write(1,*) - write(1,*) 'Note: During an evolution, This file will only' + write(1,*) 'Note: During an evolution, this file will only' write(1,*) 'contain information about the last time that the' write(1,*) 'finder was called.' write(1,*) close(1) end if end if - endif + end if + ! ******************************************* ! *** FIND {xmn,xmx,ymn,ymx,zmn,zmx} *** @@ -416,13 +396,12 @@ cGAB : should there still be an if statement here ! Find out if we want to move the center. -! offset = contains("ahf_offset","yes").ne.0 if (ahf_offset.ne.0) then offset = .true. else offset = .false. end if -! wander = contains("ahf_wander","yes").ne.0 + if (ahf_wander.ne.0) then wander = .true. else @@ -439,13 +418,11 @@ cGAB : should there still be an if statement here else nonaxi = .false. end if -! nonaxi = contains("ahf_phi","yes").ne.0 ! Find out reflection symmetries. if ((CCTK_Equals(ahf_octant,"yes").eq.1).or. - . (CCTK_Equals(ahf_octant,"high").eq.1)) then - + . (CCTK_Equals(ahf_octant,"high").eq.1)) then refx = .true. refy = .true. refz = .true. @@ -468,7 +445,6 @@ cGAB : should there still be an if statement here end if -! if (contains("ahf_phi","no").ne.0) then if (ahf_phi.eq.0) then refx = .true. refy = .true. @@ -504,7 +480,6 @@ cGAB : should there still be an if statement here cartoon = .false. #ifdef THORN_cartoon_2d -! cartoon = contains("cartoon_active","yes").ne.0 if (cartoon_active.ne.0) then cartoon = .true. else @@ -523,7 +498,6 @@ cGAB : should there still be an if statement here ! *** FIND {lmax} *** ! ************************ -! lmax = geti("ahf_lmax") lmax = ahf_lmax if (lmax.ge.20) then @@ -702,63 +676,55 @@ cGAB : should there still be an if statement here ! initializing flags - if (ahf_sloppyguess.ne.0) then - sloppy = .true. - else - sloppy = .false. - end if + sloppy = (ahf_sloppyguess.eq.1) -! sloppy = contains("ahf_sloppyguess","yes").ne.0 if (cartoon) sloppy = .true. - if (ahf_flow.ne.0) then - flow = .true. - else - flow = .false. - end if -! flow = contains("ahf_flow","yes").ne.0 - if (ahf_inner.ne.0) then - inner = .true. - else - inner = .false. - end if -! inner = contains("ahf_inner","yes").ne.0 + flow = (ahf_flow.eq.1) + + inner = (ahf_inner.eq.1) + if (ahf_guessold.ne.0) then guessold = .true. else guessold = .false. end if -! guessold = contains("ahf_guessold","yes").ne.0 + if (ahf_guess_absmin.ne.0) then guess_absmin = .true. else guess_absmin = .false. end if -! guess_absmin = contains("ahf_guess_absmin","yes").ne.0 + if (ahf_manual_guess.ne.0) then manual_guess = .true. else manual_guess = .false. end if -! manual_guess = contains("ahf_manual_guess","yes").ne.0 + if (ahf_minarea.ne.0) then minarea = .true. else minarea = .false. end if -! minarea = contains("ahf_minarea","yes").ne.0 -! get initial guess and jump into finder routines +! Get initial guess and jump into finder routines if ((.not.flow).and.guessold.and.(ahf_ncall.gt.1)) then + +! Minimization algorithm and old guess. + call AHFinder_initguess(CCTK_FARGUMENTS,rmx) call AHFinder_min(CCTK_FARGUMENTS,NN,status,logf) + else if ((.not.flow).and.manual_guess) then -! if (contains("grid","octant").eq.0) then + +! Minimization algorithm and manual guess. + if (CCTK_Equals(domain,"octant").eq.0) then write(*,*) 'AHFinder: if you want to use ahf_manual_guess' - write(*,*) ' without grid=octant, you will have to make' - write(*,*) ' it work!' + write(*,*) 'without grid=octant, you will have to make' + write(*,*) 'it work!' STOP endif @@ -772,23 +738,34 @@ cGAB : should there still be an if statement here c0_temp(14) = ahf_l14_guess c0_temp(16) = ahf_l16_guess c0_temp(18) = ahf_l18_guess - do i = 2,lmax,2 + + do i=2,lmax,2 c0(i) = c0_temp(i) enddo + call AHFinder_initguess(CCTK_FARGUMENTS,rmx) call AHFinder_min(CCTK_FARGUMENTS,NN,status,logf) + else if (.not.flow) then + +! Straight minimization algorithm. + call AHFinder_initguess(CCTK_FARGUMENTS,rmx) call AHFinder_min(CCTK_FARGUMENTS,NN,status,logf) + else - c0(0) =0.8D0*rmx + +! Flow algorithm. + if ((ahf_ncall.eq.1).and.(mfind.eq.0)) then allocate(hflow0(0:lmax),cflow0(0:lmax),nflow0(0:lmax)) allocate(hflowc(lmax,lmax),cflowc(lmax,lmax),nflowc(lmax,lmax)) allocate(hflows(lmax,lmax),cflows(lmax,lmax),nflows(lmax,lmax)) end if + call AHFinder_flow(CCTK_FARGUMENTS,rmx,status,logf) flow = .false. + end if @@ -826,8 +803,6 @@ cGAB : should there still be an if statement here call AHFinder_exp(CCTK_FARGUMENTS) call AHFinder_int(CCTK_FARGUMENTS) -! gauss = contains("ahf_gaussout","yes").ne.0 - if (ahf_gaussout.ne.0) then call AHFinder_gau(CCTK_FARGUMENTS,circ_eq,meri_p1,meri_p2) end if @@ -901,7 +876,9 @@ cGAB : should there still be an if statement here end if ! Messages to screen. + if (myproc.eq.0) then + ! Open logfile. if (logfile) then @@ -1206,6 +1183,7 @@ cGAB : should there still be an if statement here end if end if + ! **************************** ! *** WRITE DATA FILES *** ! **************************** @@ -1348,12 +1326,13 @@ c Expansion coefficients. if (status) write(1,"(2ES14.6)") cctk_time,meri_p2 close(1) end if -! modified + endif -! end modified -! ******************************************************** + + +! ****************************************************** ! *** PREPARING {ahfgrid3,ahf_exp3} FOR OUTPUT *** -! ******************************************************** +! ****************************************************** if (find3) then call AHFinder_find3(CCTK_FARGUMENTS,mtype) @@ -1418,25 +1397,21 @@ c Expansion coefficients. ! *** MASK *** ! **************** -! if (horizon.and.CONTAINS("ahf_mask","yes")) then if (horizon.and.(ahf_mask.ne.0)) then call AHFinder_mask(CCTK_FARGUMENTS) -! else if ((mtype.eq.1).and.CONTAINS("ahf_weakmask","yes")) then - else if ((mtype.eq.1).and.(ahf_weakmask.ne.0)) then + else if ((mtype.eq.1).and.(ahf_weakmask.ne.0)) then call AHFinder_mask(CCTK_FARGUMENTS) else ahmask = one end if -! call synconefunc(ahmask) call CCTK_SyncGroup(cctkGH,"ahfinder::ahfmask") + ! ************************************************** ! *** TRANSFORM PHYSICAL TO CONFORMAL METRIC *** ! ************************************************** -! call PhysToConf(nx,ny,nz,psi,gxx,gxy,gxz,gyy,gyz,gzz) - call CCTK_WARN(4,"Converting metric: physical -> conformal") gxx=gxx/psi**4 @@ -1448,6 +1423,7 @@ c Expansion coefficients. conformal_state = CONFORMAL_METRIC + ! ********************************************** ! *** UPDATING COUNTER AND LOOP IF FIND3 *** ! ********************************************** |