From d25673061d63055b0e2eef0175c73a84fda9d328 Mon Sep 17 00:00:00 2001 From: schnetter Date: Sat, 14 Feb 2004 15:04:53 +0000 Subject: Add parameters to set the origin of the excision region for "fixed excision". When checking for (real-valued) mask values, do not check for equality, but instead allow for a fudge factor. Check that the thorn is activated when any of its routines are called. When there is only a single point excised, do not abort with an internal error, but pretend the normal extends into the x direction. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/LegoExcision/trunk@42 f75ba9e5-694f-0410-ac2c-87ea7ce7132b --- param.ccl | 15 +++++++++++++++ src/extrapolate.F90 | 11 +++++++---- src/findboundary.F90 | 11 +++++++---- src/findnormals.F90 | 26 +++++++++++++++++++------- src/fixedsphere.F90 | 17 +++++++++++++---- 5 files changed, 61 insertions(+), 19 deletions(-) diff --git a/param.ccl b/param.ccl index 6769626..20176da 100644 --- a/param.ccl +++ b/param.ccl @@ -12,3 +12,18 @@ REAL fixed_size "The size (radius or corner location) of the fixed region" { 0:* :: "Positive" } 1.0 + +REAL fixed_origin_x "The x coordinate of the centre of the fixed region" +{ + *:* :: "" +} 0.0 + +REAL fixed_origin_y "The y coordinate of the centre of the fixed region" +{ + *:* :: "" +} 0.0 + +REAL fixed_origin_z "The z coordinate of the centre of the fixed region" +{ + *:* :: "" +} 0.0 diff --git a/src/extrapolate.F90 b/src/extrapolate.F90 index 19ca1b6..c09e35c 100644 --- a/src/extrapolate.F90 +++ b/src/extrapolate.F90 @@ -5,9 +5,8 @@ ! diry, and dirz for extrapolation. #include "cctk.h" -#include "cctk_Parameters.h" - #include "cctk_Functions.h" +#include "cctk_Parameters.h" #include "maskvalues.h" @@ -46,15 +45,19 @@ subroutine excision_extrapolate (ierr, var, oldvar, & integer i,j,k integer ii,jj,kk + if (CCTK_IsThornActive(CCTK_THORNSTRING) == 0) then + call CCTK_WARN (0, "The routine excision_extrapolate was called, but thorn " // CCTK_THORNSTRING // " is not active") + end if + do k=2,nk-1 do j=2,nj-1 do i=2,ni-1 - if (mask(i,j,k)==MASK_BOUNDARY) then + if (abs(mask(i,j,k)-MASK_BOUNDARY)<0.01) then ii = i + int(dirx(i,j,k)) jj = j + int(diry(i,j,k)) kk = k + int(dirz(i,j,k)) var(i,j,k) = oldvar(i,j,k) + var(ii,jj,kk) - oldvar(ii,jj,kk) - else if (mask(i,j,k)==MASK_EXCISED) then + else if (abs(mask(i,j,k)-MASK_EXCISED)<0.01) then var(i,j,k) = var0 end if end do diff --git a/src/findboundary.F90 b/src/findboundary.F90 index 454c91a..2eaa1e1 100644 --- a/src/findboundary.F90 +++ b/src/findboundary.F90 @@ -4,9 +4,8 @@ ! Return a mask where the outermost 0s have been replaced by 0.5s. #include "cctk.h" -#include "cctk_Parameters.h" - #include "cctk_Functions.h" +#include "cctk_Parameters.h" #include "maskvalues.h" @@ -34,6 +33,10 @@ subroutine excision_findboundary (ierr, mask, ni, nj, nk) integer i,j,k integer ii,jj,kk logical bnd + + if (CCTK_IsThornActive(CCTK_THORNSTRING) == 0) then + call CCTK_WARN (0, "The routine excision_findboundary was called, but thorn " // CCTK_THORNSTRING // " is not active") + end if ! Loop over grid points. @@ -43,7 +46,7 @@ subroutine excision_findboundary (ierr, mask, ni, nj, nk) ! Check if we are in an excised point. - if (mask(i,j,k)==MASK_EXCISED) then + if (abs(mask(i,j,k)-MASK_EXCISED)<0.01) then bnd = .false. @@ -53,7 +56,7 @@ subroutine excision_findboundary (ierr, mask, ni, nj, nk) do kk=k-1,k+1 do jj=j-1,j+1 do ii=i-1,i+1 - bnd = bnd.or.mask(ii,jj,kk)==MASK_ACTIVE + bnd = bnd.or.abs(mask(ii,jj,kk)-MASK_ACTIVE)<0.01 end do end do end do diff --git a/src/findnormals.F90 b/src/findnormals.F90 index ecd4d43..eb8e465 100644 --- a/src/findnormals.F90 +++ b/src/findnormals.F90 @@ -4,9 +4,8 @@ ! Return information about the normals in these locations. #include "cctk.h" -#include "cctk_Parameters.h" - #include "cctk_Functions.h" +#include "cctk_Parameters.h" #include "maskvalues.h" @@ -39,6 +38,10 @@ subroutine excision_findnormals (ierr, mask, dirx, diry, dirz, ni, nj, nk) CCTK_REAL sx,sy,sz,smag CCTK_REAL vx,vy,vz CCTK_REAL dot,dotmax + + if (CCTK_IsThornActive(CCTK_THORNSTRING) == 0) then + call CCTK_WARN (0, "The routine excision_findnormals was called, but thorn " // CCTK_THORNSTRING // " is not active") + end if ! Initialise direction arrays to zero. @@ -55,7 +58,7 @@ subroutine excision_findnormals (ierr, mask, dirx, diry, dirz, ni, nj, nk) ! Check if the point is on the boundary, if it ! is not forget about it. - if (mask(i,j,k)==MASK_BOUNDARY) then + if (abs(mask(i,j,k)-MASK_BOUNDARY)<0.01) then ! Initialise the vector (sx,sy,sz) to zero. @@ -72,7 +75,7 @@ subroutine excision_findnormals (ierr, mask, dirx, diry, dirz, ni, nj, nk) ! If a given neighbour is active, then add its ! relative coordinates to (sx,sy,sz). - if (mask(i+ii,j+jj,k+kk)==MASK_ACTIVE) then + if (abs(mask(i+ii,j+jj,k+kk)-MASK_ACTIVE)<0.01) then sx = sx + dble(ii) sy = sy + dble(jj) sz = sz + dble(kk) @@ -92,6 +95,11 @@ subroutine excision_findnormals (ierr, mask, dirx, diry, dirz, ni, nj, nk) sx = sx/smag sy = sy/smag sz = sz/smag + else + ! There is a tie. Break it. + sx = 1 + sy = 0 + sz = 0 end if ! Now we have a unit vector pointing in the average @@ -118,7 +126,7 @@ subroutine excision_findnormals (ierr, mask, dirx, diry, dirz, ni, nj, nk) ! If the point is not active, forget about it. - if (mask(i+ii,j+jj,k+kk)==MASK_ACTIVE) then + if (abs(mask(i+ii,j+jj,k+kk)-MASK_ACTIVE)<0.01) then ! Find normalized dot product. @@ -130,7 +138,7 @@ subroutine excision_findnormals (ierr, mask, dirx, diry, dirz, ni, nj, nk) / sqrt(vx**2 + vy**2 + vz**2) ! If the new product is larger than the -! largest so far, redifine the normal. +! largest so far, redefine the normal. if (dot > dotmax) then @@ -155,7 +163,11 @@ subroutine excision_findnormals (ierr, mask, dirx, diry, dirz, ni, nj, nk) jj = int(diry(i,j,k)) kk = int(dirz(i,j,k)) - if (mask(i+ii,j+jj,k+kk) /= MASK_ACTIVE) then + if (ii==0 .and. jj==0 .and. kk==0) then + call CCTK_WARN (0, "Could not decide for a normal direction") + end if + + if (abs(mask(i+ii,j+jj,k+kk)-MASK_ACTIVE)>0.01) then call CCTK_WARN (0, "Mask boundary layer is too thick") end if diff --git a/src/fixedsphere.F90 b/src/fixedsphere.F90 index 7d2294d..dafca7d 100644 --- a/src/fixedsphere.F90 +++ b/src/fixedsphere.F90 @@ -2,7 +2,6 @@ #include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" - #include "cctk_Functions.h" subroutine Lego_FixedSphere(CCTK_ARGUMENTS) @@ -12,20 +11,30 @@ subroutine Lego_FixedSphere(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS + + integer ierr + + if (CCTK_IsThornActive(CCTK_THORNSTRING) == 0) then + call CCTK_WARN (0, "The routine Lego_FixedSphere was called, but thorn " // CCTK_THORNSTRING // " is not active") + end if + + ! Only excise if there would be more than a single grid point excised + if (fixed_size < 2*minval(CCTK_DELTA_SPACE(:))) return if (CCTK_EQUALS(fixed_excision,"sphere")) then - where (r < fixed_size) + where ((x - fixed_origin_x)**2 + (y - fixed_origin_y)**2 + (z - fixed_origin_z)**2 < fixed_size**2) emask = 0.d0 elsewhere emask = 1.d0 end where else if (CCTK_EQUALS(fixed_excision,"cube")) then - where ((abs(x) < fixed_size).and.(abs(z) < fixed_size).and.& - (abs(y) < fixed_size)) + where (max(abs(x - fixed_origin_x), abs(y - fixed_origin_y), abs(z - fixed_origin_z)) < fixed_size) emask = 0.d0 elsewhere emask = 1.d0 end where end if + call excision_findboundary (ierr, emask, cctk_lsh(1), cctk_lsh(2), cctk_lsh(3)) + end subroutine Lego_FixedSphere -- cgit v1.2.3