diff options
| -rw-r--r-- | param.ccl | 215 | ||||
| -rw-r--r-- | src/AHFinder.F | 142 |
2 files changed, 241 insertions, 116 deletions
@@ -118,6 +118,11 @@ REAL dhole3_zmax "Upper z bound of excised region 2" private: + +############################ +### General parameters ### +############################ + BOOLEAN ahf_find3 "Searching for 3 surfaces?" { } "no" @@ -133,22 +138,25 @@ INT ahf_findevery "How often to look for horizons" STEERABLE = ALWAYS INT ahf_findafter "After how many iterations look for horizons" STEERABLE = ALWAYS { - : :: "By default set to 0, Any positive integer if start later in an evolution" + : :: "By default set to 0, any positive integer if start later in an evolution" } 0 REAL ahf_findaftertime "After how much time look for horizons." { - : :: "Any positive real number. Default is 0.0. Overides ahf_findafter" + : :: "Any positive real number. If non-zero overides ahf_findafter" } 0.0 REAL trapped_surface_delta "find (expansion = delta) surface" { - : :: "Just a real number. Range ??" + : :: "Just a real number" } 0.0 -# Parameters for surface expansion. -BOOLEAN ahf_phi "Expand also in phi? I.e seach for non-axisymmetric surface" STEERABLE = ALWAYS +######################################## +### Parameters for surface expansion ### +######################################## + +BOOLEAN ahf_phi "Expand also in phi? (seach for non-axisymmetric surface)" STEERABLE = ALWAYS { } "no" @@ -162,23 +170,90 @@ BOOLEAN ahf_wander "Allow the center to wander?" INT ahf_lmax "Maximum number of terms in theta expansion" STEERABLE = ALWAYS { -0:20 :: "Range from 0 to 20. 8 is the default" +0:20 :: "Range from 0 to 20" } 8 -# Parameters for surface integrals. +REAL ahf_xc "x-coordinate of center of expansion" +{ + : :: +} 0.0 + +REAL ahf_yc "y-coordinate of center of expansion" +{ + : :: +} 0.0 + +REAL ahf_zc "z-coordinate of center of expansion" +{ + : :: +} 0.0 + +REAL ahf_xc_0 "x-coordinate of center of expansion for first surface with find3" +{ + : :: +} 0.0 + +REAL ahf_yc_0 "y-coordinate of center of expansion for first surface with find3" +{ + : :: +} 0.0 + +REAL ahf_zc_0 "z-coordinate of center of expansion for first surface with find3" +{ + : :: +} 0.0 + +REAL ahf_xc_1 "x-coordinate of center of expansion for second surface with find3" +{ + : :: +} 0.0 + +REAL ahf_yc_1 "y-coordinate of center of expansion for second surface with find3" +{ + : :: +} 0.0 + +REAL ahf_zc_1 "z-coordinate of center of expansion for second surface with find3" +{ + : :: +} 0.0 + +REAL ahf_xc_2 "x-coordinate of center of expansion for third surface with find3" +{ + : :: +} 0.0 + +REAL ahf_yc_2 "y-coordinate of center of expansion for third surface with find3" +{ + : :: +} 0.0 + +REAL ahf_zc_2 "z-coordinate of center of expansion for third surface with find3" +{ + : :: +} 0.0 + + +######################################## +### Parameters for surface integrals ### +######################################## INT ahf_ntheta "Number of subdivisions in theta" { - : :: "Any sensible integer. The default of 100 is useful" + 1: :: "Any sensible integer" } 100 INT ahf_nphi "Number of subdivisions in phi" { - : :: "Any sensible integer. The default of 100 is useful" + 1: :: "Any sensible integer" } 100 -# Parameters to indicate symmetries. This symmetries refer to the -# surface itself, and NOT the grid. + +######################################### +### Parameters to indicate symmetries ### +######################################### + +# IMPORTANT: This symmetries refer to the surface itself, and NOT the grid. # Notice that octant = "yes" enforces reflection symmetries on all # three coordinate planes, while octant = "high" enforces also the @@ -196,6 +271,10 @@ BOOLEAN ahf_refz "Reflection symmetry z->-z?" STEERABLE = ALWAYS { } "no" +BOOLEAN ahf_cartoon "Cartoon mode?" STEERABLE = ALWAYS +{ +} "no" + KEYWORD ahf_octant "Octant symmetry?" STEERABLE = ALWAYS { "no" :: "No octant symmetry" @@ -203,7 +282,10 @@ KEYWORD ahf_octant "Octant symmetry?" STEERABLE = ALWAYS "high" :: "Octant symmetry + symmetry of rotation of pi/2 around z axis" } "no" -# Parameters for minimization algorithm. + +############################################# +### Parameters for minimization algorithm ### +############################################# BOOLEAN ahf_minarea "Minimize area instead of expansion?" { @@ -211,15 +293,18 @@ BOOLEAN ahf_minarea "Minimize area instead of expansion?" INT ahf_maxiter "Maximum number of iterations of POWELL" { - : :: "Any sensible integer value. Default is 10" + : :: "Any sensible integer value" } 10 REAL ahf_tol "Tolerance for minimization routines" { - : :: "A sensible real number." +0: :: "A sensible positive number" } 0.1 -#Parameters for initial guess. + +#################################### +### Parameters for initial guess ### +#################################### BOOLEAN ahf_sloppyguess "Use sphere as initial guess?" { @@ -301,87 +386,31 @@ REAL ahf_l18_guess "Manual guess for l=18 component" : :: } 0.0 -REAL ahf_r0 "Radius of initial sphere (0. forces largest sphere)" -{ - : :: -} 0.0 - -REAL ahf_r0_0 "Radius of first initial sphere for find3 (largest for 0.)" -{ - : :: -} 0.0 - -REAL ahf_r0_1 "Radius of second initial sphere for find3 (largest for 0.)" -{ - : :: -} 0.0 - -REAL ahf_r0_2 "Radius of third initial sphere for find3 (largest for 0.)" -{ - : :: -} 0.0 - -REAL ahf_xc "x-coordinate of center of expansion" -{ - : :: -} 0.0 - -REAL ahf_yc "y-coordinate of center of expansion" -{ - : :: -} 0.0 - -REAL ahf_zc "z-coordinate of center of expansion" -{ - : :: -} 0.0 - -REAL ahf_xc_0 "x-coordinate of center of expansion for first surface with find3" -{ - : :: -} 0.0 - -REAL ahf_yc_0 "y-coordinate of center of expansion for first surface with find3" -{ - : :: -} 0.0 - -REAL ahf_zc_0 "z-coordinate of center of expansion for first surface with find3" -{ - : :: -} 0.0 - -REAL ahf_xc_1 "x-coordinate of center of expansion for second surface with find3" +REAL ahf_r0 "Radius of initial sphere (0 forces largest sphere)" { : :: } 0.0 -REAL ahf_yc_1 "y-coordinate of center of expansion for second surface with find3" +REAL ahf_r0_0 "Radius of first initial sphere for find3 (0 forces largest sphere)" { : :: } 0.0 -REAL ahf_zc_1 "z-coordinate of center of expansion for second surface with find3" +REAL ahf_r0_1 "Radius of second initial sphere for find3 (0 forces largest sphere)" { : :: } 0.0 -REAL ahf_xc_2 "x-coordinate of center of expansion for third surface with find3" +REAL ahf_r0_2 "Radius of third initial sphere for find3 (0 forces largest sphere)" { : :: } 0.0 -REAL ahf_yc_2 "y-coordinate of center of expansion for third surface with find3" -{ - : :: -} 0.0 -REAL ahf_zc_2 "z-coordinate of center of expansion for third surface with find3" -{ - : :: -} 0.0 -# Parameters for flow algorithm. +##################################### +### Parameters for flow algorithm ### +##################################### BOOLEAN ahf_flow "Use flow instead of minimization?" STEERABLE = ALWAYS { @@ -432,7 +461,10 @@ REAL ahf_minchange "Minimum relative difference between 1 big and 2 small steps" : :: } 0.01 -# Parameters for output. + +############################# +### Parameters for output ### +############################# # IMPORTANT: 3D output for grid functions is still not implemented! @@ -473,7 +505,10 @@ BOOLEAN ahf_gaussout "Output gaussian curvature of horizon?" STEERABLE = ALWAYS } "yes" -# Parameters for mask. + +########################### +### Parameters for mask ### +########################### KEYWORD ahf_mask "Use mask?" { @@ -491,8 +526,6 @@ REAL ahf_maskshrink "Shrink factor for mask" 0.0:1.0 :: "Must be positive and not larger than 1" } 0.9 -# Parameters for shift. - REAL ahf_shiftcoeff "Coefficient for shift" { : :: "Anything goes" @@ -500,12 +533,14 @@ REAL ahf_shiftcoeff "Coefficient for shift" -#################################### -### PARAMETERS FROM OTHER THORNS ### -#################################### +############################################### +### PARAMETERS SHARED FROM OTHER THORNS ### +############################################### -### FROM IO ### +############### +### FROM IO ### +############### shares: IO @@ -514,7 +549,9 @@ STRING outdir "" } -### FROM GRID ### +################# +### FROM GRID ### +################# shares: grid @@ -527,7 +564,9 @@ KEYWORD bitant_plane "" } -### FROM EINSTEIN ### +##################### +### FROM EINSTEIN ### +##################### shares: einstein diff --git a/src/AHFinder.F b/src/AHFinder.F index 0f38823..c050dcf 100644 --- a/src/AHFinder.F +++ b/src/AHFinder.F @@ -409,19 +409,19 @@ if ((xc.gt.xmx).or.(xc.lt.xmn)) then write(*,*) - write(*,*) 'xc is out of the grid, giving up ...' + write(*,*) 'xc is out of the grid, giving up.' goto 5 end if if ((yc.gt.ymx).or.(yc.lt.ymn)) then write(*,*) - write(*,*) 'yc is out of the grid, giving up ...' + write(*,*) 'yc is out of the grid, giving up.' goto 5 end if if ((zc.gt.zmx).or.(zc.lt.zmn)) then write(*,*) - write(*,*) 'zc is out of the grid, giving up ...' + write(*,*) 'zc is out of the grid, giving up.' goto 5 end if @@ -455,24 +455,30 @@ nonaxi = .false. end if -! Find out reflection symmetries. +! Find out reflection symmetries depending on +! the values of the parameters. if ((CCTK_Equals(ahf_octant,"yes").eq.1).or. . (CCTK_Equals(ahf_octant,"high").eq.1)) then + refx = .true. refy = .true. refz = .true. + else + if (ahf_refx.ne.0) then refx = .true. else refx = .false. end if + if (ahf_refy.ne.0) then refy = .true. else refy = .false. end if + if (ahf_refz.ne.0) then refz = .true. else @@ -486,16 +492,26 @@ refy = .true. end if -! Force grid symmetries. +! Force grid symmetries. If the grid has some +! symmetry (octant, quadrant, bitant), and the +! horizon is centered on one of the symmetry +! planes, then the grid symmetries are forced +! on it regardless of the values of the symmetry +! parameters. if (CCTK_Equals(domain,"octant").eq.1) then + if (xc.eq.zero) refx = .true. if (yc.eq.zero) refy = .true. if (zc.eq.zero) refz = .true. + else if (CCTK_Equals(domain,"quadrant").eq.1) then + if (xc.eq.zero) refx = .true. if (yc.eq.zero) refy = .true. + else if (CCTK_Equals(domain,"bitant").eq.1) then + if (CCTK_Equals(bitant_plane,"xy").eq.1) then if (zc.eq.zero) refz = .true. else if (CCTK_Equals(bitant_plane,"xz").eq.1) then @@ -503,9 +519,37 @@ else if (xc.eq.zero) refx = .true. end if + + end if + +! Check if we are using cartoon. If we are, then the +! corresponding symmetries are forced on the horizon. + + if (ahf_cartoon.ne.0) then + cartoon = .true. + else + cartoon = .false. + end if + + if (cartoon) then + nonaxi = .false. + refx = .true. + refy = .true. + end if + +! If desired, write values of symmetry flags to screen. + + if (veryver) then + write(*,*) + write(*,*) 'Symmetries used for horizon:' + write(*,*) 'refx = ',refx + write(*,*) 'refy = ',refy + write(*,*) 'refz = ',refz + write(*,*) end if -! Find {stepx,stepy,stepz}. +! Find the values of {stepx,stepy,stepz} depending +! on the symmetries. stepx = 0 stepy = 0 @@ -517,22 +561,6 @@ if (CCTK_Equals(ahf_octant,"high").eq.1) stepx=3 -! Check if we are using cartoon. - - cartoon = .false. - -#ifdef THORN_cartoon_2d - if (cartoon_active.ne.0) then - cartoon = .true. - end if -#endif - - if (cartoon) then - nonaxi = .false. - refx = .true. - refy = .true. - end if - ! ************************ ! *** FIND {lmax} *** @@ -1024,15 +1052,10 @@ . status='old',position='append') end if -! Surface information. Notice that here I define -! the surface "mass" as: -! -! M = sqrt( A / (16 pi) ) = 0.141047396 sqrt(A) - ! Find inverse of area. Be carefull not to divide by zero! ! Also, a value of 1.0D10 indicates an error in the integration ! routine, and all integrals end up with the same value. -! So dont divide by it or I will be confused later. +! So do not divide by it or I will be confused later. if ((intarea_h.ne.zero).and.(intarea.lt.1.0D10)) then aux = one/intarea_h @@ -1040,6 +1063,11 @@ aux = one end if +! Surface information. Notice that here I define the +! surface "mass" as: +! +! M = sqrt( A / (16 pi) ) = 0.141047396 sqrt(A) + write(*,*) write(*,*) 'AHFinder: Surface found.' write(*,*) @@ -1532,6 +1560,59 @@ ! *** MASK *** ! **************** +! Initialize the sides of the excised region to the center +! of the expansion. + + if (ahf_ncall.eq.1) then + + if (find3) then + + if (mfind.eq.0) then + + dhole1_xmin = xc + dhole1_ymin = yc + dhole1_zmin = zc + + dhole1_xmax = xc + dhole1_ymax = yc + dhole1_zmax = zc + + else if (mfind.eq.1) then + + dhole2_xmin = xc + dhole2_ymin = yc + dhole2_zmin = zc + + dhole2_xmax = xc + dhole2_ymax = yc + dhole2_zmax = zc + + else + + dhole3_xmin = xc + dhole3_ymin = yc + dhole3_zmin = zc + + dhole3_xmax = xc + dhole3_ymax = yc + dhole3_zmax = zc + + end if + + else + + dhole1_xmin = xc + dhole1_ymin = yc + dhole1_zmin = zc + + dhole1_xmax = xc + dhole1_ymax = yc + dhole1_zmax = zc + + end if + + end if + ! The mask should only be different from 1 inside a horizon. ! I therefore need to check first if I do have a horizon, and ! also if the mask is desired. @@ -1571,6 +1652,11 @@ if (dhole2_ymax.lt.(yc+0.5D0*rhor)) dhole2_ymax = yc+0.5D0*rhor if (dhole2_zmax.lt.(zc+0.5D0*rhor)) dhole2_zmax = zc+0.5D0*rhor + print *, dhole2_xmin,dhole2_xmax + print *, dhole2_ymin,dhole2_ymax + print *, dhole2_zmin,dhole2_zmax + print * + else ! Third hole. |