diff options
| -rw-r--r-- | src/AHFinder_mask.F | 58 |
1 files changed, 44 insertions, 14 deletions
diff --git a/src/AHFinder_mask.F b/src/AHFinder_mask.F index 21fac55..0050f2b 100644 --- a/src/AHFinder_mask.F +++ b/src/AHFinder_mask.F @@ -14,10 +14,8 @@ subroutine AHFinder_mask(CCTK_FARGUMENTS,rhor) ! Set up the horizon mask. All this routine does is -! set up the value of the grid function "ahmask" to -! one for points outside the horizon, and to zero -! for points inside the horizon. This is only necessary -! if we want to use black hole excision. +! set up the value of the grid function "ahmask" to zero +! for points inside the horizon. ! ! It also finds out the radius of the largest sphere that ! fits inside the horizon and it centered in (xc,yc,zc). @@ -48,28 +46,31 @@ ! *** START ROUTINE *** ! ************************* - call AHFinder_fun(CCTK_FARGUMENTS) - -! If outside of horizon set mask to 1. For this +! If inside of horizon (minus a bufferzone) set the mask to 0. +! To find out how big is the buffer zone inside theFor this ! I check the parameter ahf_maskshrink. This parameter -! says by which factor of the monopole term to shrink -! the mass 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). +! 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. rhor = 1.0D10 + call AHFinder_fun(CCTK_FARGUMENTS) + do k=1,nz do j=1,ny do i=1,nx if (ahfgrid(i,j,k).ge.(ahf_maskshrink-one)*c0(0)) then - ahmask(i,j,k) = one aux = dsqrt((x(i,j,k)-xc)**2 + (y(i,j,k)-yc)**2 . + (z(i,j,k)-xc)**2) if (rhor.gt.aux) rhor = aux - else + else if (ahf_maskcube.eq.0) then ahmask(i,j,k) = zero end if @@ -78,6 +79,35 @@ end do +! ************************ +! *** CUBICAL MASK *** +! ************************ + +! If we want a cubical mask, set the mask based +! on the value of rhor. + + if (ahf_maskcube.ne.0) 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 + ahmask(i,j,k) = zero + end if + + end do + end do + end do + + end if + + ! *************** ! *** END *** ! *************** |