diff options
Diffstat (limited to 'src/EHFinder_SetMask.F90')
| -rw-r--r-- | src/EHFinder_SetMask.F90 | 86 |
1 files changed, 31 insertions, 55 deletions
diff --git a/src/EHFinder_SetMask.F90 b/src/EHFinder_SetMask.F90 index 566f654..eb13858 100644 --- a/src/EHFinder_SetMask.F90 +++ b/src/EHFinder_SetMask.F90 @@ -65,9 +65,8 @@ end subroutine EHFinder_MaskInit ! This routine will excise points and re-activate excised points -! after need. EHFinder_Setmask2 will then find the boundary of the +! as needed. EHFinder_Setmask2 will then find the boundary of the ! excsion region and make sure that the mask is correct there. - subroutine EHFinder_SetMask1(CCTK_ARGUMENTS) use EHFinder_mod @@ -218,9 +217,6 @@ subroutine EHFinder_SetMask1(CCTK_ARGUMENTS) call CCTK_WARN ( 0, 'Reduction of array fimax_loc failed' ) end if -! print*, imin_glob(1), imax_glob(1), fimin_glob(1), fimax_glob(1) -! print*, imin_glob(2), imax_glob(2), fimin_glob(2), fimax_glob(2) -! print*, imin_glob(3), imax_glob(3), fimin_glob(3), fimax_glob(3) end if ! Now check and see if any interior excised points need to be activated. @@ -307,11 +303,6 @@ subroutine EHFinder_SetMask1(CCTK_ARGUMENTS) end if end do -! print*,'New range' -! print*,imin_n(1),imax_n(1) -! print*,imin_n(2),imax_n(2) -! print*,imin_n(3),imax_n(3) - ! Use the new indices to actually activate points, taking care to ! activate points that are actually on the local grid. First do the ! faces of the rectangular box. @@ -322,11 +313,9 @@ subroutine EHFinder_SetMask1(CCTK_ARGUMENTS) j2 = min ( min ( imax_n(2), imax_glob(2) ) - cctk_lbnd(2), jyr ) k1 = max ( max ( imin_n(3), imin_glob(3) ) - cctk_lbnd(3), kzl ) k2 = min ( min ( imax_n(3), imax_glob(3) ) - cctk_lbnd(3), kzr ) -! print*,'x1 ',i1,j1,j2,k1,k2 if ( ( j1 .le. j2 ) .and. ( k1 .le. k2 ) ) then tm_mask(i1,j1:j2,k1:k2,l) = 0 ftmp(i1,j1:j2,k1:k2,l) = ftmp(i1+1,j1:j2,k1:k2,l) + delta -! print*,'done' end if end if end if @@ -337,11 +326,9 @@ subroutine EHFinder_SetMask1(CCTK_ARGUMENTS) j2 = min ( min ( imax_n(2), imax_glob(2) ) - cctk_lbnd(2), jyr ) k1 = max ( max ( imin_n(3), imin_glob(3) ) - cctk_lbnd(3), kzl ) k2 = min ( min ( imax_n(3), imax_glob(3) ) - cctk_lbnd(3), kzr ) -! print*,'x2 ',i2,j1,j2,k1,k2 if ( ( j1 .le. j2 ) .and. ( k1 .le. k2 ) ) then tm_mask(i2,j1:j2,k1:k2,l) = 0 ftmp(i2,j1:j2,k1:k2,l) = ftmp(i2-1,j1:j2,k1:k2,l) + delta -! print*,'done' end if end if end if @@ -352,11 +339,9 @@ subroutine EHFinder_SetMask1(CCTK_ARGUMENTS) i2 = min ( min ( imax_n(1), imax_glob(1) ) - cctk_lbnd(1), ixr ) k1 = max ( max ( imin_n(3), imin_glob(3) ) - cctk_lbnd(3), kzl ) k2 = min ( min ( imax_n(3), imax_glob(3) ) - cctk_lbnd(3), kzr ) -! print*,'y1 ',j1,i1,i2,k1,k2 if ( ( i1 .le. i2 ) .and. ( k1 .le. k2 ) ) then tm_mask(i1:i2,j1,k1:k2,l) = 0 ftmp(i1:i2,j1,k1:k2,l) = ftmp(i1:i2,j1+1,k1:k2,l) + delta -! print*,'done' end if end if end if @@ -367,11 +352,9 @@ subroutine EHFinder_SetMask1(CCTK_ARGUMENTS) i2 = min ( min ( imax_n(1), imax_glob(1) ) - cctk_lbnd(1), ixr ) k1 = max ( max ( imin_n(3), imin_glob(3) ) - cctk_lbnd(3), kzl ) k2 = min ( min ( imax_n(3), imax_glob(3) ) - cctk_lbnd(3), kzr ) -! print*,'y2 ',j2,i1,i2,k1,k2 if ( ( i1 .le. i2 ) .and. ( k1 .le. k2 ) ) then tm_mask(i1:i2,j2,k1:k2,l) = 0 ftmp(i1:i2,j2,k1:k2,l) = f(i1:i2,j2-1,k1:k2,l) + delta -! print*,'done' end if end if end if @@ -382,11 +365,9 @@ subroutine EHFinder_SetMask1(CCTK_ARGUMENTS) i2 = min ( min ( imax_n(1), imax_glob(1) ) - cctk_lbnd(1), ixr ) j1 = max ( max ( imin_n(2), imin_glob(2) ) - cctk_lbnd(2), jyl ) j2 = min ( min ( imax_n(2), imax_glob(2) ) - cctk_lbnd(2), jyr ) -! print*,'z1 ',k1,i1,i2,j1,j2 if ( ( i1 .le. i2 ) .and. ( j1 .le. j2 ) ) then tm_mask(i1:i2,j1:j2,k1,l) = 0 ftmp(i1:i2,j1:j2,k1,l) = ftmp(i1:i2,j1:j2,k1+1,l) + delta -! print*,'done' end if end if end if @@ -397,11 +378,9 @@ subroutine EHFinder_SetMask1(CCTK_ARGUMENTS) i2 = min ( min ( imax_n(1), imax_glob(1) ) - cctk_lbnd(1), ixr ) j1 = max ( max ( imin_n(2), imin_glob(2) ) - cctk_lbnd(2), jyl ) j2 = min ( min ( imax_n(2), imax_glob(2) ) - cctk_lbnd(2), jyr ) -! print*,'z2 ',k2,i1,i2,j1,j2 if ( ( i1 .le. i2 ) .and. ( j1 .le. j2 ) ) then tm_mask(i1:i2,j1:j2,k2,l) = 0 ftmp(i1:i2,j1:j2,k2,l) = ftmp(i1:i2,j1:j2,k2-1,l) + delta -! print*,'done' end if end if end if @@ -674,7 +653,6 @@ subroutine EHFinder_SetMask1(CCTK_ARGUMENTS) i2 = imax_n(1) - cctk_lbnd(1) j2 = imax_n(2) - cctk_lbnd(2) k2 = imax_n(3) - cctk_lbnd(3) -! print*,'Debug 3 : ',i2,j2,k2 if ( ( ixl .le. i2 ) .and. ( i2 .le. min(ixr,nx-1) ) .and. & ( jyl .le. j2 ) .and. ( j2 .le. min(jyr,ny-1) ) .and. & ( kzl .le. k2 ) .and. ( k2 .le. min(kzr,nz-1) ) ) then @@ -829,7 +807,6 @@ subroutine EHFinder_SetMask2(CCTK_ARGUMENTS) tm_mask(ixl:ixr,jyl:jyr,kzl:kzr,l) end where end if -! call CCTK_INFO('Mask modified') end do ! Indicate that it is not anymore the first time the mask is set. @@ -839,6 +816,10 @@ subroutine EHFinder_SetMask2(CCTK_ARGUMENTS) end subroutine EHFinder_SetMask2 +! This routine checks if any of the excision regions are to close together +! so that some points do not have enough active neigbouring points to be +! able to calculate derivatives. If this is the case these points are +! excised as well. subroutine EHFinder_SetMask3(CCTK_ARGUMENTS) use EHFinder_mod @@ -884,81 +865,76 @@ subroutine EHFinder_SetMask3(CCTK_ARGUMENTS) do j = jyl, jyr do i = ixl, ixr if ( eh_mask(i,j,k,l) .ge. 0 ) then - ! If we are on an excision boundary in the -x-direction... + +! If we are on an excision boundary in the -x-direction... if ( btest(eh_mask(i,j,k,l),0) ) then - ! If any of the two closest neighbours in the +x-direction is - ! excised we have to excise this point + +! If any of the two closest neighbours in the +x-direction is +! excised we have to excise this point if ( ( eh_mask(i+1,j,k,l) .eq. -1 ) .or. & ( eh_mask(i+2,j,k,l) .eq. -1 ) ) then tm_mask(i,j,k,l) = -1 mask_modified = .true. - print*, i, j, k, 'x-' - print*, eh_mask(i:i+2,j,k,l) end if end if - ! If we are on an excision boundary in the +x-direction... +! If we are on an excision boundary in the +x-direction... if ( btest(eh_mask(i,j,k,l),1) ) then - ! If any of the two closest neighbours in the -x-direction is - ! excised we have to excise this point + +! If any of the two closest neighbours in the -x-direction is +! excised we have to excise this point if ( ( eh_mask(i-1,j,k,l) .eq. -1 ) .or. & ( eh_mask(i-2,j,k,l) .eq. -1 ) ) then tm_mask(i,j,k,l) = -1 mask_modified = .true. - print*, i, j, k, 'x+' - print*, eh_mask(i:i-2,j,k,l) end if end if - ! If we are on an excision boundary in the -y-direction... +! If we are on an excision boundary in the -y-direction... if ( btest(eh_mask(i,j,k,l),2) ) then - ! If any of the two closest neighbours in the +y-direction is - ! excised we have to excise this point + +! If any of the two closest neighbours in the +y-direction is +! excised we have to excise this point if ( ( eh_mask(i,j+1,k,l) .eq. -1 ) .or. & ( eh_mask(i,j+2,k,l) .eq. -1 ) ) then tm_mask(i,j,k,l) = -1 mask_modified = .true. - print*, i, j, k, 'y-' - print*, eh_mask(i,j:j+2,k,l) end if end if - ! If we are on an excision boundary in the +y-direction... +! If we are on an excision boundary in the +y-direction... if ( btest(eh_mask(i,j,k,l),3) ) then - ! If any of the two closest neighbours in the -y-direction is - ! excised we have to excise this point + +! If any of the two closest neighbours in the -y-direction is +! excised we have to excise this point if ( ( eh_mask(i,j-1,k,l) .eq. -1 ) .or. & ( eh_mask(i,j-2,k,l) .eq. -1 ) ) then tm_mask(i,j,k,l) = -1 mask_modified = .true. - print*, i, j, k, 'y+' - print*, eh_mask(i,j:j-2,k,l) end if end if - ! If we are on an excision boundary in the -z-direction... +! If we are on an excision boundary in the -z-direction... if ( btest(eh_mask(i,j,k,l),4) ) then - ! If any of the two closest neighbours in the +z-direction is - ! excised we have to excise this point + +! If any of the two closest neighbours in the +z-direction is +! excised we have to excise this point if ( ( eh_mask(i,j,k+1,l) .eq. -1 ) .or. & ( eh_mask(i,j,k+2,l) .eq. -1 ) ) then tm_mask(i,j,k,l) = -1 mask_modified = .true. - print*, i, j, k, 'z-' - print*, eh_mask(i,j,k:k+2,l) end if end if - ! If we are on an excision boundary in the +z-direction... +! If we are on an excision boundary in the +z-direction... if ( btest(eh_mask(i,j,k,l),5) ) then - ! If any of the two closest neighbours in the -z-direction is - ! excised we have to excise this point + +! If any of the two closest neighbours in the -z-direction is +! excised we have to excise this point if ( ( eh_mask(i,j,k-1,l) .eq. -1 ) .or. & ( eh_mask(i,j,k-2,l) .eq. -1 ) ) then tm_mask(i,j,k,l) = -1 mask_modified = .true. - print*, i, j, k, 'z+' - print*, eh_mask(i,j,k:k-2,l) end if end if end if @@ -968,7 +944,7 @@ subroutine EHFinder_SetMask3(CCTK_ARGUMENTS) end do if ( mask_modified ) then - call CCTK_WARN ( 1, "Mask modified because it was not convex" ) + call CCTK_WARN ( 1, 'Mask modified because it was not convex' ) end if end if end do |