diff options
| -rw-r--r-- | param.ccl | 9 | ||||
| -rw-r--r-- | src/AHFinder_mask.F | 153 |
2 files changed, 117 insertions, 45 deletions
@@ -427,14 +427,17 @@ KEYWORD ahf_mask "Use mask?" "weak" :: "Mask set for both definite and probable horizons" } "off" -BOOLEAN ahf_maskcube "Do we want a cubical mask?" +KEYWORD ahf_masktype "Type of mask" { -} "no" +"lego" :: "Mask is a lego sphere" +"cube" :: "Mask is a cube" +"poly" :: "Mask is a polyhedra" +} "cube" REAL ahf_maskshrink "Shrink factor for mask" { 0.0:1.0 :: "Must be positive and not larger than 1" -} 0.9 +} 0.8 REAL ahf_shiftcoeff "Coefficient for shift" { diff --git a/src/AHFinder_mask.F b/src/AHFinder_mask.F index d3b881f..8a5900e 100644 --- a/src/AHFinder_mask.F +++ b/src/AHFinder_mask.F @@ -27,65 +27,72 @@ DECLARE_CCTK_FARGUMENTS DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS integer i,j,k integer reduce_handle CCTK_REAL zero,one - CCTK_REAL aux,rhor,rhortemp + CCTK_REAL rhor,rhortemp + CCTK_REAL buffer,aux + CCTK_REAL xa,ya,za -! ***************************** -! *** DEFINE {zero,one} *** -! ***************************** +! ******************* +! *** NUMBERS *** +! ******************* zero = 0.0D0 one = 1.0D0 -! ************************* -! *** START ROUTINE *** -! ************************* +! ********************************* +! *** FIND HORIZON FUNCTION *** +! ********************************* -! If inside of horizon (minus a bufferzone) set the mask to 0. -! To find out how big is the buffer zone I check the value -! of the parameter ahf_maskshrink. This parameter says -! by which factor of the monopole term to shrink the mask -! with respect to the real horizon. Here I use the fact that -! the function ahfgrid is zero at the horizon, and is linear -! on the monopole term (which measures the overall scale). -! -! Notice that if I want a cubical mask, I first find out the radius -! of the largest sphere that fits inside the horizon and only after -! that I set the mask. + call AHFinder_fun(CCTK_FARGUMENTS) - rhor = 1.0D10 - call AHFinder_fun(CCTK_FARGUMENTS) +! *********************************** +! *** FIND BUFFER ZONE FACTOR *** +! *********************************** + +! The mask will always be set inside the horizon minus +! a safety buffer zone. To find out how big is the buffer +! zone I check the value of the parameter "ahf_maskshrink". +! This parameter says by which factor of the monopole term +! to shrink the mask with respect to the real horizon. +! Here I use the fact that the function ahfgrid is zero +! at the horizon, and is linear on the monopole term +! (which measures the overall scale). + + buffer = (ahf_maskshrink-one)*c0(0) + + +! ***************************************** +! *** FIND RADIUS OF LARGEST SPHERE *** +! ***************************************** + +! Find the radius of the largest sphere that +! fits inside the horizon. + + rhor = 1.0D10 do k=1,nz do j=1,ny do i=1,nx - if (ahfgrid(i,j,k).ge.(ahf_maskshrink-one)*c0(0)) then + if (ahfgrid(i,j,k).ge.buffer) then aux = dsqrt((x(i,j,k)-xc)**2 + (y(i,j,k)-yc)**2 . + (z(i,j,k)-zc)**2) if (rhor.gt.aux) rhor = aux - else if (ahf_maskcube.eq.0) then - if (ahfgrid(i,j,k).lt.(0.8D0*ahf_maskshrink-one)*c0(0)) then - ahmask(i,j,k) = zero - else - ahmask(i,j,k) = (ahfgrid(i,j,k) - . - (0.8D0*ahf_maskshrink-one)*c0(0)) - . / (0.2D0*ahf_maskshrink*c0(0)) - end if end if end do end do end do -! Now find the minumum of rhor across processors. +! Now find the minimum of rhor across processors. call CCTK_ReductionArrayHandle(reduce_handle,"minimum") call CCTK_ReduceLocalScalar(ierr ,cctkGH,-1,reduce_handle, @@ -97,25 +104,87 @@ end if +! *********************** +! *** LEGO SPHERE *** +! *********************** + +! If we want a lego shere mask life is simple. +! Just set the mask inside of horizon (minus a +! buffer zone) equal to 0. + + if (CCTK_EQUALS(ahf_masktype,'lego')) then + + do k=1,nz + do j=1,ny + do i=1,nx + + if (ahfgrid(i,j,k).lt.buffer) then + ahmask(i,j,k) = zero + end if + + end do + end do + end do + + ! ************************ ! *** CUBICAL MASK *** ! ************************ -! If we want a cubical mask, set the mask based -! on the value of rhor. +! If we want a cubical mask, set the mask based on the value +! of rhor. Notice that if we want the corners of the cube to +! touch the sphere of radius rhor, we must make the side of +! the cube equal to rhor/sqrt(3), which is about 0.57*rhor. + + else if (CCTK_EQUALS(ahf_masktype,'cube')) then + + do k=1,nz + do j=1,ny + do i=1,nx + + aux = 0.57D0*rhor + + xa = dabs(x(i,j,k)-xc) + ya = dabs(y(i,j,k)-yc) + za = dabs(z(i,j,k)-zc) + + if ((xa.lt.aux).and.(ya.lt.aux).and.(za.lt.aux)) then + ahmask(i,j,k) = zero + end if + + end do + end do + end do - if (ahf_maskcube.ne.0) then + +! ************************** +! *** POLYHEDRA MASK *** +! ************************** + +! Here we want to excise a polyhedra. But not just any +! polyhedra, a polyhedra defined by taking a cube and +! cutting away the corners all the way up to the middle +! of the edges. This results in a figure made up of +! 6 squares and 8 equilateral triangles. + +! Notice that for such a figure, if we want the corners +! to touch a sphere of radius rhor, then the sides have +! to be rhor/sqrt(2) which is roughly 0.7*rhor. + + else if (CCTK_EQUALS(ahf_masktype,'poly')) then do k=1,nz do j=1,ny do i=1,nx - if ((x(i,j,k).gt.(xc-0.5D0*rhor)).and. - . (y(i,j,k).gt.(yc-0.5D0*rhor)).and. - . (z(i,j,k).gt.(zc-0.5D0*rhor)).and. - . (x(i,j,k).lt.(xc+0.5D0*rhor)).and. - . (y(i,j,k).lt.(yc+0.5D0*rhor)).and. - . (z(i,j,k).lt.(zc+0.5D0*rhor))) then + aux = 0.7D0*rhor + + xa = dabs(x(i,j,k)-xc) + ya = dabs(y(i,j,k)-yc) + za = dabs(z(i,j,k)-zc) + + if ((xa.lt.aux).and.(ya.lt.aux).and.(za.lt.aux).and. + . ((xa+ya+za).lt.2.0D0*aux)) then ahmask(i,j,k) = zero end if @@ -126,9 +195,9 @@ end if -! **************************************************** -! *** APPLY SYMMETRY BOUNDARIES AND SYNCHRONIZE *** -! **************************************************** +! ***************************************************** +! *** APPLY SYMMETRY BOUNDARIES AND SYNCHRONIZE *** +! ***************************************************** call CCTK_SyncGroup(cctkGH,"ahfinder::ahfmask") call CartSymGN(ierr,cctkGH,"ahfinder::ahfmask") |