From eefab44a3f055102e6ac7d5f81e82501692a71f0 Mon Sep 17 00:00:00 2001 From: diener Date: Wed, 23 Feb 2005 16:55:52 +0000 Subject: Added a parameter to only start setting the excision mask after a certain time. The default value is set so that the normal behaviour is not changed. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinder/trunk@376 89daf98e-ef62-4674-b946-b8ff9de2216c --- param.ccl | 5 + src/AHFinder_mask.F | 272 ++++++++++++++++++++++++++-------------------------- 2 files changed, 142 insertions(+), 135 deletions(-) diff --git a/param.ccl b/param.ccl index debc0b0..10161ff 100644 --- a/param.ccl +++ b/param.ccl @@ -453,6 +453,11 @@ KEYWORD ahf_masktype "Type of mask" STEERABLE = ALWAYS "poly" :: "Mask is a polyhedra" } "cube" +REAL ahf_mask_time "Time after which to start setting the mask" +{ +: :: "Anything goes. Negative number means setting the mask as soon as possible" +} -1.0 + BOOLEAN ahf_mask_0 "Mask for first horizon with find3?" { } "yes" diff --git a/src/AHFinder_mask.F b/src/AHFinder_mask.F index de8f347..4153a27 100644 --- a/src/AHFinder_mask.F +++ b/src/AHFinder_mask.F @@ -46,195 +46,197 @@ zero = 0.0D0 one = 1.0D0 + if ( (ahf_mask_time .lt. 0.0) .or. (cctk_time .ge. ahf_mask_time) ) then -! ********************************* -! *** FIND HORIZON FUNCTION *** -! ********************************* +! ********************************* +! *** FIND HORIZON FUNCTION *** +! ********************************* - call AHFinder_fun(CCTK_ARGUMENTS) + call AHFinder_fun(CCTK_ARGUMENTS) -! *********************************** -! *** FIND BUFFER ZONE FACTOR *** -! *********************************** +! *********************************** +! *** 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. +! 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). What I then do is -! set up the value of "buffer" to a negative number that -! roughly measures the desired coordinate distance from -! the mask to the horizon. The reason why it is negative -! is because the horizon function is negative inside the -! horizon, and I just want to compare its value with the -! value of "buffer" to decide when a given point is masked. +! 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). What I then do is +! set up the value of "buffer" to a negative number that +! roughly measures the desired coordinate distance from +! the mask to the horizon. The reason why it is negative +! is because the horizon function is negative inside the +! horizon, and I just want to compare its value with the +! value of "buffer" to decide when a given point is masked. - buffer = - (one-ahf_maskshrink)*c0(0) + buffer = - (one-ahf_maskshrink)*c0(0) -! Now I make sure that there are at least "ahf_maskbuffer" -! grid points in the buffer zone. The final buffer zone will -! be whatever is largest between the number of grid points -! specified by "ahf_maskbuffer" and the shrink factor -! specified by "ahf_maskshrink". +! Now I make sure that there are at least "ahf_maskbuffer" +! grid points in the buffer zone. The final buffer zone will +! be whatever is largest between the number of grid points +! specified by "ahf_maskbuffer" and the shrink factor +! specified by "ahf_maskshrink". - dd = max(dx,dy,dz) - aux = ahf_maskbuffer*dd + dd = max(dx,dy,dz) + aux = ahf_maskbuffer*dd - if (abs(buffer).lt.aux) buffer = - aux + if (abs(buffer).lt.aux) buffer = - aux -! ***************************************** -! *** FIND RADIUS OF LARGEST SPHERE *** -! ***************************************** +! ***************************************** +! *** FIND RADIUS OF LARGEST SPHERE *** +! ***************************************** -! Find the radius of the largest sphere that -! fits inside the horizon. +! Find the radius of the largest sphere that +! fits inside the horizon. - rhor = 1.0D10 + rhor = 1.0D10 - do k=1,nz - do j=1,ny - do i=1,nx + do k=1,nz + do j=1,ny + do i=1,nx - if (ahfgrid(i,j,k).ge.buffer) then - aux = sqrt((x(i,j,k)-xc)**2 + (y(i,j,k)-yc)**2 - . + (z(i,j,k)-zc)**2) - if (rhor.gt.aux) rhor = aux - end if + if (ahfgrid(i,j,k).ge.buffer) then + aux = sqrt((x(i,j,k)-xc)**2 + (y(i,j,k)-yc)**2 + . + (z(i,j,k)-zc)**2) + if (rhor.gt.aux) rhor = aux + end if - end do - end do - end do + end do + end do + end do -! Now find the minimum of rhor across processors. +! Now find the minimum of rhor across processors. - reduce_handle = -1 - call CCTK_ReductionArrayHandle(reduce_handle,"minimum") - if (reduce_handle.lt.0) then - call CCTK_WARN(1,"Cannot get handle for minimum reduction ! Forgot to activate an implementation providing reduction operators ??") - end if + reduce_handle = -1 + call CCTK_ReductionArrayHandle(reduce_handle,"minimum") + if (reduce_handle.lt.0) then + call CCTK_WARN(1,"Cannot get handle for minimum reduction ! Forgot to activate an implementation providing reduction operators ??") + end if - call CCTK_ReduceLocalScalar(ierr ,cctkGH,-1,reduce_handle, - . rhor,rhortemp,CCTK_VARIABLE_REAL) - rhor = rhortemp + call CCTK_ReduceLocalScalar(ierr ,cctkGH,-1,reduce_handle, + . rhor,rhortemp,CCTK_VARIABLE_REAL) + rhor = rhortemp - if (ierr.ne.0) then - call CCTK_WARN(1,"AHFinder_mask:Reduction of rhor failed!") - end if + if (ierr.ne.0) then + call CCTK_WARN(1,"AHFinder_mask:Reduction of rhor failed!") + end if -! *********************** -! *** LEGO SPHERE *** -! *********************** +! *********************** +! *** 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 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 (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 + if (ahfgrid(i,j,k).lt.buffer) then + ahmask(i,j,k) = zero + end if - end do - end do - end do + end do + end do + end do -! ************************ -! *** CUBICAL MASK *** -! ************************ +! ************************ +! *** CUBICAL MASK *** +! ************************ -! 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. +! 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 + else if (CCTK_EQUALS(ahf_masktype,'cube')) then - do k=1,nz - do j=1,ny - do i=1,nx + do k=1,nz + do j=1,ny + do i=1,nx - aux = 0.57D0*rhor + aux = 0.57D0*rhor - xa = abs(x(i,j,k)-xc) - ya = abs(y(i,j,k)-yc) - za = abs(z(i,j,k)-zc) + xa = abs(x(i,j,k)-xc) + ya = abs(y(i,j,k)-yc) + za = abs(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 + 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 + end do + end do + end do -! ************************** -! *** POLYHEDRA MASK *** -! ************************** +! ************************** +! *** 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. +! 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. +! 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 + else if (CCTK_EQUALS(ahf_masktype,'poly')) then - do k=1,nz - do j=1,ny - do i=1,nx + do k=1,nz + do j=1,ny + do i=1,nx - aux = 0.7D0*rhor + aux = 0.7D0*rhor - xa = abs(x(i,j,k)-xc) - ya = abs(y(i,j,k)-yc) - za = abs(z(i,j,k)-zc) + xa = abs(x(i,j,k)-xc) + ya = abs(y(i,j,k)-yc) + za = abs(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 + 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 - end do - end do - end do + end do + end do + end do - end if + end if -! ***************************************************** -! *** APPLY SYMMETRY BOUNDARIES AND SYNCHRONIZE *** -! ***************************************************** +! ***************************************************** +! *** APPLY SYMMETRY BOUNDARIES AND SYNCHRONIZE *** +! ***************************************************** - call CCTK_SyncGroup(ierr,cctkGH,"ahfinder::ahfmask") - call CartSymGN(ierr,cctkGH,"ahfinder::ahfmask") + call CCTK_SyncGroup(ierr,cctkGH,"ahfinder::ahfmask") + call CartSymGN(ierr,cctkGH,"ahfinder::ahfmask") -! If desired, copy ahf_mask to spacemask mask (emask). +! If desired, copy ahf_mask to spacemask mask (emask). - if (use_mask.eq.1) then - emask = ahmask - end if + if (use_mask.eq.1) then + emask = ahmask + end if -! *************** -! *** END *** -! *************** +! *************** +! *** END *** +! *************** + end if end subroutine AHFinder_mask -- cgit v1.2.3