From ade581e1acd584ea6e54ff7a2aad11bc23477f9d Mon Sep 17 00:00:00 2001 From: diener Date: Wed, 13 Aug 2003 17:11:22 +0000 Subject: Major modification of almost all source files to use vector groups for various variables. Not tested in all detail, but the standard features seem to work. The version before all these changes was tagged with PRE_MULTI. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/EHFinder/trunk@126 2a26948c-0e4f-0410-aee8-f1d3e353619c --- src/EHFinder_Check.F90 | 140 ++-- src/EHFinder_FindSurface.F90 | 161 ++-- src/EHFinder_Generator_Sources.F90 | 299 ++++---- src/EHFinder_Generator_Sources2.F90 | 217 +++--- src/EHFinder_Init.F90 | 165 ++-- src/EHFinder_Integrate2.F90 | 140 +--- src/EHFinder_ReParametrize.F90 | 906 +++++++++++----------- src/EHFinder_SetMask.F90 | 1452 ++++++++++++++++++----------------- src/EHFinder_SetSym.F90 | 100 ++- src/EHFinder_Sources.F90 | 192 ++--- src/EHFinder_mod.F90 | 3 +- src/MoL_Init.F90 | 41 +- 12 files changed, 1906 insertions(+), 1910 deletions(-) (limited to 'src') diff --git a/src/EHFinder_Check.F90 b/src/EHFinder_Check.F90 index 0ef8f28..4a6e07a 100644 --- a/src/EHFinder_Check.F90 +++ b/src/EHFinder_Check.F90 @@ -14,7 +14,7 @@ subroutine EHFinder_ReParametrize_Check(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS - CCTK_INT :: i, j, k + CCTK_INT :: i, j, k, l CCTK_REAL :: fmin_extremum, fmin_extremum_loc, fmax_loc, fmin_bound_loc character(len=128) :: info_message logical :: check_reparam @@ -42,79 +42,81 @@ subroutine EHFinder_ReParametrize_Check(CCTK_ARGUMENTS) #include "include/physical_part2.h" if ( check_reparam ) then - fmin_bound_loc = huge -! print*,cctk_lbnd -! print*,cctk_ubnd -! print*,cctk_gsh - if ( cctk_lbnd(1) .eq. 0 .and. .not.symx ) then - fmin_bound_loc = min ( fmin_bound_loc, minval ( f(1,:,:) ) ) - end if - if ( cctk_ubnd(1) .eq. cctk_gsh(1)-1 ) then - fmin_bound_loc = min ( fmin_bound_loc, minval ( f(nx,:,:) ) ) - end if - if ( cctk_lbnd(2) .eq. 0 .and. .not.symy ) then - fmin_bound_loc = min ( fmin_bound_loc, minval ( f(:,1,:) ) ) - end if - if ( cctk_ubnd(2) .eq. cctk_gsh(2)-1 ) then - fmin_bound_loc = min ( fmin_bound_loc, minval ( f(:,ny,:) ) ) - end if - if ( cctk_lbnd(3) .eq. 0 .and. .not.symz ) then - fmin_bound_loc = min ( fmin_bound_loc, minval ( f(:,:,1) ) ) - end if - if ( cctk_ubnd(3) .eq. cctk_gsh(3)-1 ) then - fmin_bound_loc = min ( fmin_bound_loc, minval ( f(:,:,nz) ) ) - end if -! print*,fmin_bound_loc + do l = 1, eh_number_level_sets + fmin_bound_loc = huge + if ( cctk_lbnd(1) .eq. 0 .and. .not.symx ) then + fmin_bound_loc = min ( fmin_bound_loc, minval ( f(1,:,:,l) ) ) + end if + if ( cctk_ubnd(1) .eq. cctk_gsh(1)-1 ) then + fmin_bound_loc = min ( fmin_bound_loc, minval ( f(nx,:,:,l) ) ) + end if + if ( cctk_lbnd(2) .eq. 0 .and. .not.symy ) then + fmin_bound_loc = min ( fmin_bound_loc, minval ( f(:,1,:,l) ) ) + end if + if ( cctk_ubnd(2) .eq. cctk_gsh(2)-1 ) then + fmin_bound_loc = min ( fmin_bound_loc, minval ( f(:,ny,:,l) ) ) + end if + if ( cctk_lbnd(3) .eq. 0 .and. .not.symz ) then + fmin_bound_loc = min ( fmin_bound_loc, minval ( f(:,:,1,l) ) ) + end if + if ( cctk_ubnd(3) .eq. cctk_gsh(3)-1 ) then + fmin_bound_loc = min ( fmin_bound_loc, minval ( f(:,:,nz,l) ) ) + end if + ! print*,fmin_bound_loc + + fmax_loc = zero - fmax_loc = zero - - fmin_extremum_loc = ex_value - do k = kzl, kzr - do j = jyl, jyr - do i = ixl, ixr - fmax_loc = max ( f(i,j,k), fmax_loc ) - if ( f(i,j,k) .lt. zero ) then - if ( ( (f(i,j,k)-f(i-1,j,k))*(f(i+1,j,k)-f(i,j,k)).le.zero ).and. & - ( (f(i,j,k)-f(i,j-1,k))*(f(i,j+1,k)-f(i,j,k)).le.zero ).and. & - ( (f(i,j,k)-f(i,j,k-1))*(f(i,j,k+1)-f(i,j,k)).le.zero ) ) then - fmin_extremum_loc = max ( f(i,j,k), fmin_extremum_loc ) + fmin_extremum_loc = ex_value + do k = kzl, kzr + do j = jyl, jyr + do i = ixl, ixr + fmax_loc = max ( f(i,j,k,l), fmax_loc ) + if ( f(i,j,k,l) .lt. zero ) then + if ( ( (f(i,j,k,l)-f(i-1,j,k,l))*(f(i+1,j,k,l)-f(i,j,k,l)).le.zero ).and. & + ( (f(i,j,k,l)-f(i,j-1,k,l))*(f(i,j+1,k,l)-f(i,j,k,l)).le.zero ).and. & + ( (f(i,j,k,l)-f(i,j,k-1,l))*(f(i,j,k+1,l)-f(i,j,k,l)).le.zero ) ) then + fmin_extremum_loc = max ( f(i,j,k,l), fmin_extremum_loc ) + end if end if - end if + end do end do end do - end do - - call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, max_handle, & - fmin_extremum_loc, fmin_extremum, CCTK_VARIABLE_REAL ) - - if ( ierr .ne. 0 ) then - call CCTK_WARN(0,'Reduction of fmax_extremum failed') - end if - - call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, min_handle, & - fmin_bound_loc, fmin_bound, CCTK_VARIABLE_REAL ) - - if ( ierr .ne. 0 ) then - call CCTK_WARN(0,'Reduction of fmax failed') - end if - - write(info_message,'(a30,f12.5)') & - 'Minimum f at the extrema is : ',fmin_extremum - call CCTK_INFO(info_message) - - write(info_message,'(a30,f12.5)') & - 'Limit for centered differences is : ', & - fmin_bound - three*maxval ( cctk_delta_space ) - call CCTK_INFO(info_message) - - if ( reparam_undo .gt. 0 ) then - if ( fmin_extremum .gt. -two * minval(cctk_delta_space) ) then - call CCTK_INFO('Re-parametrization undone due to imminent pinch-off') - f = fbak - eh_mask = eh_mask_bak - reparam_undone = .true. + + call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, max_handle, & + fmin_extremum_loc, fmin_extremum, CCTK_VARIABLE_REAL ) + + if ( ierr .ne. 0 ) then + call CCTK_WARN(0,'Reduction of fmax_extremum failed') end if - end if + + call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, min_handle, & + fmin_bound_loc, fmin_bound, CCTK_VARIABLE_REAL ) + + if ( ierr .ne. 0 ) then + call CCTK_WARN(0,'Reduction of fmax failed') + end if + + write(info_message,'(a39,i3,a6,f12.5)') & + 'Minimum f at the extrema for level set ', l, ' is : ',fmin_extremum + call CCTK_INFO(info_message) + + write(info_message,'(a45,i3,a6,f12.5)') & + 'Limit for centered differences for level set ', l, ' is : ', & + fmin_bound - three*maxval ( cctk_delta_space ) + call CCTK_INFO(info_message) + + if ( reparam_undo .gt. 0 ) then + if ( fmin_extremum .gt. -two * minval(cctk_delta_space) ) then + write(info_message,'(a45,i3,a6,f12.5)') & + 'Re-parametrization of level set ', l, & + 'undone due to imminent pinch-off' + call CCTK_INFO(info_message) + f(:,:,:,l) = fbak(:,:,:,l) + eh_mask(:,:,:,l) = eh_mask_bak(:,:,:,l) + reparam_undone(l) = .true. + end if + end if + end do end if end subroutine EHFinder_ReParametrize_Check diff --git a/src/EHFinder_FindSurface.F90 b/src/EHFinder_FindSurface.F90 index 95e0e63..2537408 100644 --- a/src/EHFinder_FindSurface.F90 +++ b/src/EHFinder_FindSurface.F90 @@ -15,15 +15,16 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS - CCTK_INT :: i, j, k, im, jm, sn + CCTK_INT :: i, j, k, l, im, jm, sn CCTK_REAL, parameter :: eps = 1.0d-10 CCTK_INT :: interp_handle, table_handle, status, coord_system_handle CCTK_INT :: interp_handle2 - character(len=30) :: info_message + character(len=80) :: info_message character(len=200) :: surface_interp CCTK_INT :: surface_interp_len character(len=7) :: surface_order + character(len=15) :: vname CCTK_INT, dimension(4) :: bbox CCTK_INT, dimension(2) :: gsh, lsh, lbnd, ubnd, nghost, maxl @@ -54,7 +55,9 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) ! set to -1. find_surface_status = 0 sn = surface_counter - integrate_counter + 1 - write(info_message,'(a26,i3)') 'Finding points on surface ',sn + l = levelset_counter + write(info_message,'(a26,i3,a14,i3)') 'Finding points on surface ', sn, & + ' in level set ',l call CCTK_INFO ( info_message ) ! Find information about the local and global properties of the 2D @@ -103,13 +106,13 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) do k = kzl, kzr do j = jyl, jyr do i = ixl, ixr - if ( eh_mask(i,j,k) .eq. 0 ) then - if ( ( ( f(i,j,k) * f(i-1,j,k) .lt. zero ) .or. & - ( f(i,j,k) * f(i+1,j,k) .lt. zero ) .or. & - ( f(i,j,k) * f(i,j-1,k) .lt. zero ) .or. & - ( f(i,j,k) * f(i,j+1,k) .lt. zero ) .or. & - ( f(i,j,k) * f(i,j,k-1) .lt. zero ) .or. & - ( f(i,j,k) * f(i,j,k+1) .lt. zero ) ) .and. & + if ( eh_mask(i,j,k,l) .eq. 0 ) then + if ( ( ( f(i,j,k,l) * f(i-1,j,k,l) .lt. zero ) .or. & + ( f(i,j,k,l) * f(i+1,j,k,l) .lt. zero ) .or. & + ( f(i,j,k,l) * f(i,j-1,k,l) .lt. zero ) .or. & + ( f(i,j,k,l) * f(i,j+1,k,l) .lt. zero ) .or. & + ( f(i,j,k,l) * f(i,j,k-1,l) .lt. zero ) .or. & + ( f(i,j,k,l) * f(i,j,k+1,l) .lt. zero ) ) .and. & ( ( sc(i,j,k) .eq. sn ) .or. & ( sc(i-1,j,k) .eq. sn ) .or. & ( sc(i+1,j,k) .eq. sn ) .or. & @@ -352,16 +355,6 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) theta_sym_factor = two * pi / ( utheta - ltheta ) phi_sym_factor = two * pi / ( uphi - lphi ) -! print* -! print*,ltheta/pi,utheta/pi,lphi/pi,uphi/pi -! print*,sym_factor -! print*,center -! print*,symx,symy,symz -! ! For now symmetries are not handled properly -! if ( symx ) center(1) = zero -! if ( symy ) center(2) = zero -! if ( symz ) center(3) = zero - ! Find dtheta and dphi and initialise the theta and phi grid arrays. ! Here i + lbnd(1) - 1 is the global index for theta and ! j + lbnd(2) - 1 is the global index for phi. @@ -384,13 +377,6 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) costheta = cos(ctheta) sinphi = sin(cphi) cosphi = cos(cphi) -! do j = 1, lsh(2) -! do i = 1, lsh(1) -! sinthetacosphi(i,j) = sin(ctheta(i,j)) * cos(cphi(i,j)) -! sinthetasinphi(i,j) = sin(ctheta(i,j)) * sin(cphi(i,j)) -! costheta(i,j) = cos(ctheta(i,j)) -! end do -! end do ! rsurf is the radial grid array (as function of theta and phi) and is ! initialised to zero. The grid array drsurf contains the changes to rsurf @@ -442,7 +428,8 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) out_array(1) = CCTK_PointerTo(f_interp) ! Get the variable index for f - call CCTK_VarIndex ( in_array, "ehfinder::f" ) + write(vname,'(a12,i1,a1)') 'ehfinder::f[', l-1, ']' + call CCTK_VarIndex ( in_array, vname ) ! Loop to find starting point in all directions at the same time. For ! simplicity all ghost zones are found on all processors @@ -492,9 +479,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) maxf_loc, maxf, CCTK_VARIABLE_REAL ) maxl = maxloc ( abs ( f_interp ) ) -! print*,maxl,rsurf(maxl(1),maxl(2)),& -! drsurf(maxl(1),maxl(2)),& -! f_interp(maxl(1),maxl(2)) + if ( maxdr .lt. min_delta ) exit find_starting_points do j = 1, lsh(2) @@ -559,14 +544,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) maxf_loc, maxf, CCTK_VARIABLE_REAL ) maxl = maxloc ( abs ( f_interp ) ) -! print*,k,maxf -! print*,maxf,maxl,f_interp(maxl(1),maxl(2)),rsurf(maxl(1),maxl(2)), & -! dfdx_interp(maxl(1),maxl(2))*sintheta(maxl(1),maxl(2))* & -! cosphi(maxl(1),maxl(2))+ & -! dfdy_interp(maxl(1),maxl(2))*sintheta(maxl(1),maxl(2))* & -! sinphi(maxl(1),maxl(2))+ & -! dfdz_interp(maxl(1),maxl(2))*costheta(maxl(1),maxl(2)) -! print* + if ( maxf .gt. eps ) then do j = 1, lsh(2) do i = 1, lsh(1) @@ -694,6 +672,42 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) end subroutine EHFinder_FindSurface +subroutine EHFinder_LevelSetLoopInit(CCTK_ARGUMENTS) + + use EHFinder_mod + + implicit none + + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + + more_levelsets = 1 + levelset_counter = 1 + +end subroutine EHFinder_LevelSetLoopInit + + +subroutine EHFinder_UpdateCounter(CCTK_ARGUMENTS) + + use EHFinder_mod + + implicit none + + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + + if ( levelset_counter .lt. eh_number_level_sets ) then + levelset_counter = levelset_counter + 1 + more_levelsets = 1 + else + more_levelsets = 0 + end if + +end subroutine EHFinder_UpdateCounter + + subroutine EHFinder_CountSurfacesInit(CCTK_ARGUMENTS) use EHFinder_mod @@ -732,12 +746,12 @@ subroutine EHFinder_CountSurfaces(CCTK_ARGUMENTS) ! Find the location of the minimum among active points not already found ! to be inside a surface. - minl = minloc ( f(ixl:ixr,jyl:jyr,kzl:kzr), & + minl = minloc ( f(ixl:ixr,jyl:jyr,kzl:kzr,levelset_counter), & ( sc(ixl:ixr,jyl:jyr,kzl:kzr) .eq. 0 ) .and. & - ( eh_mask(ixl:ixr,jyl:jyr,kzl:kzr) .ge. 0 ) ) + ( eh_mask(ixl:ixr,jyl:jyr,kzl:kzr,levelset_counter) .ge. 0 ) ) ! Find the value of f at that location. - minf_loc = f(minl(1)+ixl-1,minl(2)+jyl-1,minl(3)+kzl-1) + minf_loc = f(minl(1)+ixl-1,minl(2)+jyl-1,minl(3)+kzl-1,levelset_counter) ! Find the minimum over all processors. call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, min_handle, & @@ -810,8 +824,8 @@ subroutine EHFinder_MarkSurfaces(CCTK_ARGUMENTS) do j = jyl, jyr do i = ixl, ixr marked = .false. - if ( ( eh_mask(i,j,k) .ge. 0 ) .and. & - ( f(i,j,k) .lt. 0 ) .and. & + if ( ( eh_mask(i,j,k,levelset_counter) .ge. 0 ) .and. & + ( f(i,j,k,levelset_counter) .lt. 0 ) .and. & ( sc(i,j,k) .ne. surface_counter ) ) then ! Find out if we are on the boundary @@ -857,58 +871,6 @@ subroutine EHFinder_MarkSurfaces(CCTK_ARGUMENTS) marked = .true. end if -! if ( any ( sc(i+i1:i+i2,j+j1:j+j2,k+k1:k+k2) & -! .eq. surface_counter ) ) then -! marked = .true. -! end if - -! ! First look in the -x direction. -! if ( .not. btest(eh_mask(i,j,k),0) ) then -! if ( sc(i-1,j,k) .eq. surface_counter ) then -! marked = .true. -! end if -! if ( .not. btest(eh_mask(i-1,j,k),2) ) then -! if ( sc(i-1,j-1,k) .eq. surface_counter ) then -! marked = .true. -! end if -! if ( .not. btest(eh_mask(i-1,j,k),3) ) then -! if ( sc(i-1,j+1,k) .eq. surface_counter ) then -! marked = .true. -! end if -! if ( .not. btest(eh_mask(i-1,j,k),4) ) then -! if ( sc(i-1,j,k-1) .eq. surface_counter ) then -! marked = .true. -! end if -! if ( .not. btest(eh_mask(i-1,j,k),5) ) then -! if ( sc(i-1,j,k+1) .eq. surface_counter ) then -! marked = .true. -! end if -! end if -! if ( .not. btest(eh_mask(i,j,k),1) ) then -! if ( sc(i+1,j,k) .eq. surface_counter ) then -! marked = .true. -! end if -! end if -! if ( .not. btest(eh_mask(i,j,k),2) ) then -! if ( sc(i,j-1,k) .eq. surface_counter ) then -! marked = .true. -! end if -! end if -! if ( .not. btest(eh_mask(i,j,k),3) ) then -! if ( sc(i,j+1,k) .eq. surface_counter ) then -! marked = .true. -! end if -! end if -! if ( .not. btest(eh_mask(i,j,k),4) ) then -! if ( sc(i,j,k-1) .eq. surface_counter ) then -! marked = .true. -! end if -! end if -! if ( .not. btest(eh_mask(i,j,k),5) ) then -! if ( sc(i,j,k+1) .eq. surface_counter ) then -! marked = .true. -! end if -! end if if ( marked ) then sc(i,j,k) = surface_counter n_loc = n_loc + 1 @@ -958,13 +920,14 @@ subroutine EHFinder_InfoSurfaces(CCTK_ARGUMENTS) ! Write out how many surfaces were found. if ( surface_counter .eq. 1 ) then - write(info_message,'(a10,i3,a8)') 'There are ',surface_counter,' surface' + write(info_message,'(a9,i3,a22,i3)') 'There is ',surface_counter, & + ' surface in level set ',levelset_counter else - write(info_message,'(a10,i3,a9)') 'There are ',surface_counter,' surfaces' + write(info_message,'(a10,i3,a23,i3)') 'There are ',surface_counter, & + ' surfaces in level set ',levelset_counter end if call CCTK_INFO(info_message) integrate_counter = surface_counter -! print*,'EHFinder_InfoSurfaces : ',integrate_counter end subroutine EHFinder_InfoSurfaces diff --git a/src/EHFinder_Generator_Sources.F90 b/src/EHFinder_Generator_Sources.F90 index a89a548..bb8e350 100644 --- a/src/EHFinder_Generator_Sources.F90 +++ b/src/EHFinder_Generator_Sources.F90 @@ -15,12 +15,13 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS - CCTK_INT :: i + CCTK_INT :: i, l CCTK_INT :: interp_handle, table_handle, status, coord_system_handle character(len=200) :: gen_interp CCTK_INT :: gen_interp_len character(len=7) :: gen_order + character(len=15) :: vname CCTK_INT, dimension(1) :: lsh CCTK_POINTER, dimension(3) :: interp_coords @@ -70,16 +71,11 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS) endif ! Find out how many interpolation points are located on this processor. - call CCTK_GrouplshGN ( status, cctkGH, 1, lsh, "ehfinder::generators" ) + call CCTK_GrouplshGN ( status, cctkGH, 1, lsh, "ehfinder::xg" ) if ( status .lt. 0 ) then call CCTK_WARN ( 0, "cannot get local size for surface arrays" ) end if - ! Set the pointers to the points to be interpolated to. - interp_coords(1) = CCTK_PointerTo(xg) - interp_coords(2) = CCTK_PointerTo(yg) - interp_coords(3) = CCTK_PointerTo(zg) - ! Set the pointers to the output arrays. out_arrays(1) = CCTK_PointerTo(alpg) out_arrays(2) = CCTK_PointerTo(betaxg) @@ -96,147 +92,158 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS) out_arrays(13) = CCTK_PointerTo(dfzg) out_arrays(14) = CCTK_PointerTo(psig) - ! Set the indices to the input grid functions. - call CCTK_VarIndex ( in_arrays(1), "admbase::alp" ) - call CCTK_VarIndex ( in_arrays(2), "admbase::betax" ) - call CCTK_VarIndex ( in_arrays(3), "admbase::betay" ) - call CCTK_VarIndex ( in_arrays(4), "admbase::betaz" ) - call CCTK_VarIndex ( in_arrays(5), "admbase::gxx" ) - call CCTK_VarIndex ( in_arrays(6), "admbase::gxy" ) - call CCTK_VarIndex ( in_arrays(7), "admbase::gxz" ) - call CCTK_VarIndex ( in_arrays(8), "admbase::gyy" ) - call CCTK_VarIndex ( in_arrays(9), "admbase::gyz" ) - call CCTK_VarIndex ( in_arrays(10), "admbase::gzz" ) - call CCTK_VarIndex ( in_arrays(11), "ehfinder::f" ) - call CCTK_VarIndex ( in_arrays(12), "staticconformal::psi" ) - - ! Check the metric type. At present physical and static_conformal are - ! supported. - if ( CCTK_EQUALS ( metric_type, 'physical' ) ) then - - ! Set the operand indices table entry, corresponding - ! to interpolation of admbase::alp (1), admbase::shift (3), - ! admbase::metric(6) and the first derivatives of ehfinder::f (3) for - ! a total of 13 output arrays. - call Util_TableSetIntArray ( status, table_handle, 13, & - op_indices(1:13), "operand_indices" ) - if ( status .lt. 0 ) then - call CCTK_WARN ( 0, "Cannot set operand indices array in parameter table" ) - endif - - ! Set the corresponding table entry for the operation codes. - call Util_TableSetIntArray ( status, table_handle, 13, & - op_codes(1:13), "operation_codes" ) - if ( status .lt. 0 ) then - call CCTK_WARN ( 0, "Cannot set operation codes array in parameter table" ) - endif - - ! Call the interpolator. - call CCTK_InterpGridArrays ( status, cctkGH, 3, interp_handle, & - table_handle, coord_system_handle, & - lsh(1), CCTK_VARIABLE_REAL, & - interp_coords, 11, in_arrays(1:11), & - 13, out_types(1:13), out_arrays(1:13) ) - - if ( status .lt. 0 ) then - call CCTK_INFO ( 'Interpolation failed.' ) - end if - - ! For each point on this processor calculate the right hand side of the - ! characteristic evolution equation. - do i = 1, lsh(1) - - ! calculate the square of the lapse. - alp2 = alpg(i)**2 - - ! Calculate the inverse of the 3-metric. - guxx = gyyg(i) * gzzg(i) - gyzg(i)**2 - guxy = gxzg(i) * gyzg(i) - gxyg(i) * gzzg(i) - guxz = gxyg(i) * gyzg(i) - gxzg(i) * gyyg(i) - - idetg = one / ( gxxg(i) * guxx + gxyg(i) * guxy + gxzg(i) * guxz ) - - guxx = idetg * guxx - guxy = idetg * guxy - guxz = idetg * guxz - - guyy = ( gxxg(i) * gzzg(i) - gxzg(i)**2 ) * idetg - guyz = ( gxyg(i) * gxzg(i) - gxxg(i) * gyzg(i) ) * idetg - guzz = ( gxxg(i) * gyyg(i) - gxyg(i)**2 ) * idetg + ! Loop over the level sets + do l = 1, eh_number_level_sets + + ! Set the pointers to the points to be interpolated to. + interp_coords(1) = CCTK_PointerTo(xg(:,l)) + interp_coords(2) = CCTK_PointerTo(yg(:,l)) + interp_coords(3) = CCTK_PointerTo(zg(:,l)) + + ! Set the indices to the input grid functions. + call CCTK_VarIndex ( in_arrays(1), "admbase::alp" ) + call CCTK_VarIndex ( in_arrays(2), "admbase::betax" ) + call CCTK_VarIndex ( in_arrays(3), "admbase::betay" ) + call CCTK_VarIndex ( in_arrays(4), "admbase::betaz" ) + call CCTK_VarIndex ( in_arrays(5), "admbase::gxx" ) + call CCTK_VarIndex ( in_arrays(6), "admbase::gxy" ) + call CCTK_VarIndex ( in_arrays(7), "admbase::gxz" ) + call CCTK_VarIndex ( in_arrays(8), "admbase::gyy" ) + call CCTK_VarIndex ( in_arrays(9), "admbase::gyz" ) + call CCTK_VarIndex ( in_arrays(10), "admbase::gzz" ) + write(vname,'(a12,i1,a1)') 'ehfinder::f[', l-1,']' + call CCTK_VarIndex ( in_arrays(11), vname ) + call CCTK_VarIndex ( in_arrays(12), "staticconformal::psi" ) + + ! Check the metric type. At present physical and static_conformal are + ! supported. + if ( CCTK_EQUALS ( metric_type, 'physical' ) ) then + + ! Set the operand indices table entry, corresponding + ! to interpolation of admbase::alp (1), admbase::shift (3), + ! admbase::metric(6) and the first derivatives of ehfinder::f (3) for + ! a total of 13 output arrays. + call Util_TableSetIntArray ( status, table_handle, 13, & + op_indices(1:13), "operand_indices" ) + if ( status .lt. 0 ) then + call CCTK_WARN ( 0, "Cannot set operand indices array in parameter table" ) + endif + + ! Set the corresponding table entry for the operation codes. + call Util_TableSetIntArray ( status, table_handle, 13, & + op_codes(1:13), "operation_codes" ) + if ( status .lt. 0 ) then + call CCTK_WARN ( 0, "Cannot set operation codes array in parameter table" ) + endif + + ! Call the interpolator. + call CCTK_InterpGridArrays ( status, cctkGH, 3, interp_handle, & + table_handle, coord_system_handle, & + lsh(1), CCTK_VARIABLE_REAL, & + interp_coords, 11, in_arrays(1:11), & + 13, out_types(1:13), out_arrays(1:13) ) + + if ( status .lt. 0 ) then + call CCTK_INFO ( 'Interpolation failed.' ) + end if + + ! For each point on this processor calculate the right hand side of the + ! characteristic evolution equation. + do i = 1, lsh(1) + + ! calculate the square of the lapse. + alp2 = alpg(i)**2 + + ! Calculate the inverse of the 3-metric. + guxx = gyyg(i) * gzzg(i) - gyzg(i)**2 + guxy = gxzg(i) * gyzg(i) - gxyg(i) * gzzg(i) + guxz = gxyg(i) * gyzg(i) - gxzg(i) * gyyg(i) + + idetg = one / ( gxxg(i) * guxx + gxyg(i) * guxy + gxzg(i) * guxz ) + + guxx = idetg * guxx + guxy = idetg * guxy + guxz = idetg * guxz + + guyy = ( gxxg(i) * gzzg(i) - gxzg(i)**2 ) * idetg + guyz = ( gxyg(i) * gxzg(i) - gxxg(i) * gyzg(i) ) * idetg + guzz = ( gxxg(i) * gyyg(i) - gxyg(i)**2 ) * idetg - ! Raise the index of the partial derivatives of f. - dfux = guxx * dfxg(i) + guxy * dfyg(i) + guxz * dfzg(i) - dfuy = guxy * dfxg(i) + guyy * dfyg(i) + guyz * dfzg(i) - dfuz = guxz * dfxg(i) + guyz * dfyg(i) + guzz * dfzg(i) - - ! Calculate the overall multiplication factor. - factor = alp2 / sqrt ( alp2 * ( dfux * dfxg(i) + & - dfuy * dfyg(i) + & - dfuz * dfzg(i) ) ) - - ! Finally obtain dx^i/dt. - dxg(i) = - betaxg(i) + ssign * factor * dfux - dyg(i) = - betayg(i) + ssign * factor * dfuy - dzg(i) = - betazg(i) + ssign * factor * dfuz - end do - else if ( CCTK_EQUALS ( metric_type, 'static conformal' ) ) then - call Util_TableSetIntArray ( status, table_handle, 14, & - op_indices, "operand_indices" ) - if ( status .lt. 0 ) then - call CCTK_WARN ( 0, "Cannot set operand indices array in parameter table" ) - endif - - call Util_TableSetIntArray ( status, table_handle, 14, & - op_codes, "operation_codes" ) - if ( status .lt. 0 ) then - call CCTK_WARN ( 0, "Cannot set operation codes array in parameter table" ) - endif - - call CCTK_InterpGridArrays ( status, cctkGH, 3, interp_handle, & - table_handle, coord_system_handle, & - lsh(1), CCTK_VARIABLE_REAL, & - interp_coords, 12, in_arrays, & - 14, out_types, out_arrays ) - - if ( status .lt. 0 ) then - call CCTK_INFO ( 'Interpolation failed.' ) + ! Raise the index of the partial derivatives of f. + dfux = guxx * dfxg(i) + guxy * dfyg(i) + guxz * dfzg(i) + dfuy = guxy * dfxg(i) + guyy * dfyg(i) + guyz * dfzg(i) + dfuz = guxz * dfxg(i) + guyz * dfyg(i) + guzz * dfzg(i) + + ! Calculate the overall multiplication factor. + factor = alp2 / sqrt ( alp2 * ( dfux * dfxg(i) + & + dfuy * dfyg(i) + & + dfuz * dfzg(i) ) ) + + ! Finally obtain dx^i/dt. + dxg(i,l) = - betaxg(i) + ssign * factor * dfux + dyg(i,l) = - betayg(i) + ssign * factor * dfuy + dzg(i,l) = - betazg(i) + ssign * factor * dfuz + end do + else if ( CCTK_EQUALS ( metric_type, 'static conformal' ) ) then + call Util_TableSetIntArray ( status, table_handle, 14, & + op_indices, "operand_indices" ) + if ( status .lt. 0 ) then + call CCTK_WARN ( 0, "Cannot set operand indices array in parameter table" ) + endif + + call Util_TableSetIntArray ( status, table_handle, 14, & + op_codes, "operation_codes" ) + if ( status .lt. 0 ) then + call CCTK_WARN ( 0, "Cannot set operation codes array in parameter table" ) + endif + + call CCTK_InterpGridArrays ( status, cctkGH, 3, interp_handle, & + table_handle, coord_system_handle, & + lsh(1), CCTK_VARIABLE_REAL, & + interp_coords, 12, in_arrays, & + 14, out_types, out_arrays ) + + if ( status .lt. 0 ) then + call CCTK_INFO ( 'Interpolation failed.' ) + end if + + do i = 1, lsh(1) + alp2 = alpg(i)**2 + +! The inverse of psi^4 + psi4 = one / psig(i)**4 + + guxx = gyyg(i) * gzzg(i) - gyzg(i)**2 + guxy = gxzg(i) * gyzg(i) - gxyg(i) * gzzg(i) + guxz = gxyg(i) * gyzg(i) - gxzg(i) * gyyg(i) + +! The determinant divided by psi^4. + idetg = psi4 / ( gxxg(i) * guxx + & + gxyg(i) * guxy + & + gxzg(i) * guxz ) + +! The inverse metric. Since the determinant is already divided +! by psi^4, this gives the inverse of the physical metric. + guxx = idetg * guxx + guxy = idetg * guxy + guxz = idetg * guxz + + guyy = ( gxxg(i) * gzzg(i) - gxzg(i)**2 ) * idetg + guyz = ( gxyg(i) * gxzg(i) - gxxg(i) * gyzg(i) ) * idetg + guzz = ( gxxg(i) * gyyg(i) - gxyg(i)**2 ) * idetg + + dfux = guxx * dfxg(i) + guxy * dfyg(i) + guxz * dfzg(i) + dfuy = guxy * dfxg(i) + guyy * dfyg(i) + guyz * dfzg(i) + dfuz = guxz * dfxg(i) + guyz * dfyg(i) + guzz * dfzg(i) + factor = alp2 / sqrt ( alp2 * ( dfux * dfxg(i) + & + dfuy * dfyg(i) + & + dfuz * dfzg(i) ) ) + dxg(i,l) = - betaxg(i) + ssign * factor * dfux + dyg(i,l) = - betayg(i) + ssign * factor * dfuy + dzg(i,l) = - betazg(i) + ssign * factor * dfuz + end do end if + end do - do i = 1, lsh(1) - alp2 = alpg(i)**2 - -! The inverse of psi^4 - psi4 = one / psig(i)**4 - - guxx = gyyg(i) * gzzg(i) - gyzg(i)**2 - guxy = gxzg(i) * gyzg(i) - gxyg(i) * gzzg(i) - guxz = gxyg(i) * gyzg(i) - gxzg(i) * gyyg(i) - -! The determinant divided by psi^4. - idetg = psi4 / ( gxxg(i) * guxx + & - gxyg(i) * guxy + & - gxzg(i) * guxz ) - -! The inverse metric. Since the determinant is already divided -! by psi^4, this gives the inverse of the physical metric. - guxx = idetg * guxx - guxy = idetg * guxy - guxz = idetg * guxz - - guyy = ( gxxg(i) * gzzg(i) - gxzg(i)**2 ) * idetg - guyz = ( gxyg(i) * gxzg(i) - gxxg(i) * gyzg(i) ) * idetg - guzz = ( gxxg(i) * gyyg(i) - gxyg(i)**2 ) * idetg - - dfux = guxx * dfxg(i) + guxy * dfyg(i) + guxz * dfzg(i) - dfuy = guxy * dfxg(i) + guyy * dfyg(i) + guyz * dfzg(i) - dfuz = guxz * dfxg(i) + guyz * dfyg(i) + guzz * dfzg(i) - factor = alp2 / sqrt ( alp2 * ( dfux * dfxg(i) + & - dfuy * dfyg(i) + & - dfuz * dfzg(i) ) ) - dxg(i) = - betaxg(i) + ssign * factor * dfux - dyg(i) = - betayg(i) + ssign * factor * dfuy - dzg(i) = - betazg(i) + ssign * factor * dfuz - end do - end if return end subroutine EHFinder_Generator_Sources diff --git a/src/EHFinder_Generator_Sources2.F90 b/src/EHFinder_Generator_Sources2.F90 index b50f9f1..b557381 100644 --- a/src/EHFinder_Generator_Sources2.F90 +++ b/src/EHFinder_Generator_Sources2.F90 @@ -15,7 +15,7 @@ subroutine EHFinder_Generator_Sources2(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS - CCTK_INT :: i, j, k + CCTK_INT :: i, j, k, l CCTK_INT :: interp_handle, table_handle, status, coord_system_handle character(len=200) :: gen_interp @@ -68,21 +68,11 @@ subroutine EHFinder_Generator_Sources2(CCTK_ARGUMENTS) #include "include/physical_part.h" ! Find out how many interpolation points are located on this processor. - call CCTK_GrouplshGN ( status, cctkGH, 1, lsh, "ehfinder::generators" ) + call CCTK_GrouplshGN ( status, cctkGH, 1, lsh, "ehfinder::xg" ) if ( status .lt. 0 ) then call CCTK_WARN ( 0, "cannot get local size for surface arrays" ) end if - ! Set the pointers to the points to be interpolated to. - interp_coords(1) = CCTK_PointerTo(xg) - interp_coords(2) = CCTK_PointerTo(yg) - interp_coords(3) = CCTK_PointerTo(zg) - - ! Set the pointers to the output arrays. - out_arrays(1) = CCTK_PointerTo(dxg) - out_arrays(2) = CCTK_PointerTo(dyg) - out_arrays(3) = CCTK_PointerTo(dzg) - ! Set the indices to the input grid functions. call CCTK_VarIndex ( in_arrays(1), "ehfinder::xgf" ) call CCTK_VarIndex ( in_arrays(2), "ehfinder::ygf" ) @@ -103,107 +93,132 @@ subroutine EHFinder_Generator_Sources2(CCTK_ARGUMENTS) call CCTK_WARN ( 0, "Cannot set operation codes array in parameter table" ) endif - ! Check the metric type. At present physical and static_conformal are - ! supported. - if ( CCTK_EQUALS ( metric_type, 'physical' ) ) then - - do k = kzl, kzr - do j = jyl, jyr - do i = ixl, ixr - if ( eh_mask(i,j,k) .ge. 0 ) then - - ! calculate the square of the lapse. - alp2 = alp(i,j,k)**2 - - ! Calculate the inverse of the 3-metric. - guxx = gyy(i,j,k) * gzz(i,j,k) - gyz(i,j,k)**2 - guxy = gxz(i,j,k) * gyz(i,j,k) - gxy(i,j,k) * gzz(i,j,k) - guxz = gxy(i,j,k) * gyz(i,j,k) - gxz(i,j,k) * gyy(i,j,k) - - idetg = one / ( gxx(i,j,k) * guxx + & - gxy(i,j,k) * guxy + & - gxz(i,j,k) * guxz ) - - guxx = idetg * guxx - guxy = idetg * guxy - guxz = idetg * guxz - - guyy = ( gxx(i,j,k) * gzz(i,j,k) - gxz(i,j,k)**2 ) * idetg - guyz = ( gxy(i,j,k) * gxz(i,j,k) - gxx(i,j,k) * gyz(i,j,k) ) * idetg - guzz = ( gxx(i,j,k) * gyy(i,j,k) - gxy(i,j,k)**2 ) * idetg - - ! Raise the index of the partial derivatives of f. - dfux = guxx * dfx(i,j,k) + guxy * dfy(i,j,k) + guxz * dfz(i,j,k) - dfuy = guxy * dfx(i,j,k) + guyy * dfy(i,j,k) + guyz * dfz(i,j,k) - dfuz = guxz * dfx(i,j,k) + guyz * dfy(i,j,k) + guzz * dfz(i,j,k) - - ! Calculate the overall multiplication factor. - factor = alp2 / sqrt ( alp2 * ( dfux * dfx(i,j,k) + & - dfuy * dfy(i,j,k) + & - dfuz * dfz(i,j,k) ) ) - - ! Finally obtain dx^i/dt. - xgf(i,j,k) = - betax(i,j,k) + ssign * factor * dfux - ygf(i,j,k) = - betay(i,j,k) + ssign * factor * dfuy - zgf(i,j,k) = - betaz(i,j,k) + ssign * factor * dfuz - else - xgf(i,j,k) = zero - ygf(i,j,k) = zero - zgf(i,j,k) = zero - end if + ! Loop over the level sets + do l = 1, eh_number_level_sets + + ! Set the pointers to the points to be interpolated to. + interp_coords(1) = CCTK_PointerTo(xg(:,l)) + interp_coords(2) = CCTK_PointerTo(yg(:,l)) + interp_coords(3) = CCTK_PointerTo(zg(:,l)) + + ! Set the pointers to the output arrays. + out_arrays(1) = CCTK_PointerTo(dxg(:,l)) + out_arrays(2) = CCTK_PointerTo(dyg(:,l)) + out_arrays(3) = CCTK_PointerTo(dzg(:,l)) + + ! Check the metric type. At present physical and static_conformal are + ! supported. + if ( CCTK_EQUALS ( metric_type, 'physical' ) ) then + + do k = kzl, kzr + do j = jyl, jyr + do i = ixl, ixr + if ( eh_mask(i,j,k,l) .ge. 0 ) then + + ! calculate the square of the lapse. + alp2 = alp(i,j,k)**2 + + ! Calculate the inverse of the 3-metric. + guxx = gyy(i,j,k) * gzz(i,j,k) - gyz(i,j,k)**2 + guxy = gxz(i,j,k) * gyz(i,j,k) - gxy(i,j,k) * gzz(i,j,k) + guxz = gxy(i,j,k) * gyz(i,j,k) - gxz(i,j,k) * gyy(i,j,k) + + idetg = one / ( gxx(i,j,k) * guxx + & + gxy(i,j,k) * guxy + & + gxz(i,j,k) * guxz ) + + guxx = idetg * guxx + guxy = idetg * guxy + guxz = idetg * guxz + + guyy = ( gxx(i,j,k) * gzz(i,j,k) - gxz(i,j,k)**2 ) * idetg + guyz = ( gxy(i,j,k) * gxz(i,j,k) - gxx(i,j,k) * gyz(i,j,k) ) * idetg + guzz = ( gxx(i,j,k) * gyy(i,j,k) - gxy(i,j,k)**2 ) * idetg + + ! Raise the index of the partial derivatives of f. + dfux = guxx * dfx(i,j,k,l) + guxy * dfy(i,j,k,l) + guxz * dfz(i,j,k,l) + dfuy = guxy * dfx(i,j,k,l) + guyy * dfy(i,j,k,l) + guyz * dfz(i,j,k,l) + dfuz = guxz * dfx(i,j,k,l) + guyz * dfy(i,j,k,l) + guzz * dfz(i,j,k,l) + + ! Calculate the overall multiplication factor. + factor = alp2 / sqrt ( alp2 * ( dfux * dfx(i,j,k,l) + & + dfuy * dfy(i,j,k,l) + & + dfuz * dfz(i,j,k,l) ) ) + + ! Finally obtain dx^i/dt. + xgf(i,j,k) = - betax(i,j,k) + ssign * factor * dfux + ygf(i,j,k) = - betay(i,j,k) + ssign * factor * dfuy + zgf(i,j,k) = - betaz(i,j,k) + ssign * factor * dfuz + else + xgf(i,j,k) = zero + ygf(i,j,k) = zero + zgf(i,j,k) = zero + end if + end do end do end do - end do - else if ( CCTK_EQUALS ( metric_type, 'static conformal' ) ) then + else if ( CCTK_EQUALS ( metric_type, 'static conformal' ) ) then - do i = 1, lsh(1) - alp2 = alpg(i)**2 + do k = kzl, kzr + do j = jyl, jyr + do i = ixl, ixr + if ( eh_mask(i,j,k,l) .ge. 0 ) then + alp2 = alp(i,j,k)**2 -! The inverse of psi^4 - psi4 = one / psig(i)**4 + ! The inverse of psi^4 + psi4 = one / psi(i,j,k)**4 - guxx = gyyg(i) * gzzg(i) - gyzg(i)**2 - guxy = gxzg(i) * gyzg(i) - gxyg(i) * gzzg(i) - guxz = gxyg(i) * gyzg(i) - gxzg(i) * gyyg(i) + guxx = gyy(i,j,k) * gzz(i,j,k) - gyz(i,j,k)**2 + guxy = gxz(i,j,k) * gyz(i,j,k) - gxy(i,j,k) * gzz(i,j,k) + guxz = gxy(i,j,k) * gyz(i,j,k) - gxz(i,j,k) * gyy(i,j,k) -! The determinant divided by psi^4. - idetg = psi4 / ( gxxg(i) * guxx + & - gxyg(i) * guxy + & - gxzg(i) * guxz ) +! The determinant divided by psi^4. + idetg = psi4 / ( gxx(i,j,k) * guxx + & + gxy(i,j,k) * guxy + & + gxz(i,j,k) * guxz ) -! The inverse metric. Since the determinant is already divided -! by psi^4, this gives the inverse of the physical metric. - guxx = idetg * guxx - guxy = idetg * guxy - guxz = idetg * guxz +! The inverse metric. Since the determinant is already divided +! by psi^4, this gives the inverse of the physical metric. + guxx = idetg * guxx + guxy = idetg * guxy + guxz = idetg * guxz - guyy = ( gxxg(i) * gzzg(i) - gxzg(i)**2 ) * idetg - guyz = ( gxyg(i) * gxzg(i) - gxxg(i) * gyzg(i) ) * idetg - guzz = ( gxxg(i) * gyyg(i) - gxyg(i)**2 ) * idetg + guyy = ( gxx(i,j,k) * gzz(i,j,k) - gxz(i,j,k)**2 ) * idetg + guyz = ( gxy(i,j,k) * gxz(i,j,k) - gxx(i,j,k) * gyz(i,j,k) ) * idetg + guzz = ( gxx(i,j,k) * gyy(i,j,k) - gxy(i,j,k)**2 ) * idetg - dfux = guxx * dfxg(i) + guxy * dfyg(i) + guxz * dfzg(i) - dfuy = guxy * dfxg(i) + guyy * dfyg(i) + guyz * dfzg(i) - dfuz = guxz * dfxg(i) + guyz * dfyg(i) + guzz * dfzg(i) - factor = alp2 / sqrt ( alp2 * ( dfux * dfxg(i) + & - dfuy * dfyg(i) + & - dfuz * dfzg(i) ) ) - dxg(i) = - betaxg(i) + ssign * factor * dfux - dyg(i) = - betayg(i) + ssign * factor * dfuy - dzg(i) = - betazg(i) + ssign * factor * dfuz - end do - end if + dfux = guxx * dfx(i,j,k,l) + guxy * dfy(i,j,k,l) + guxz * dfz(i,j,k,l) + dfuy = guxy * dfx(i,j,k,l) + guyy * dfy(i,j,k,l) + guyz * dfz(i,j,k,l) + dfuz = guxz * dfx(i,j,k,l) + guyz * dfy(i,j,k,l) + guzz * dfz(i,j,k,l) + factor = alp2 / sqrt ( alp2 * ( dfux * dfx(i,j,k,l) + & + dfuy * dfy(i,j,k,l) + & + dfuz * dfz(i,j,k,l) ) ) + xgf(i,j,k) = - betax(i,j,k) + ssign * factor * dfux + ygf(i,j,k) = - betay(i,j,k) + ssign * factor * dfuy + zgf(i,j,k) = - betaz(i,j,k) + ssign * factor * dfuz + else + xgf(i,j,k) = zero + ygf(i,j,k) = zero + zgf(i,j,k) = zero + end if + end do + end do + end do + end if - ! Call the interpolator. - call CCTK_InterpGridArrays ( status, cctkGH, 3, interp_handle, & - table_handle, coord_system_handle, & - lsh(1), CCTK_VARIABLE_REAL, & - interp_coords, 3, in_arrays, & - 3, out_types, out_arrays ) + ! Call the interpolator. + call CCTK_InterpGridArrays ( status, cctkGH, 3, interp_handle, & + table_handle, coord_system_handle, & + lsh(1), CCTK_VARIABLE_REAL, & + interp_coords, 3, in_arrays, & + 3, out_types, out_arrays ) - if ( status .lt. 0 ) then - call CCTK_INFO ( 'Interpolation failed.' ) - end if + if ( status .lt. 0 ) then + call CCTK_INFO ( 'Interpolation failed.' ) + end if + + end do return end subroutine EHFinder_Generator_Sources2 diff --git a/src/EHFinder_Init.F90 b/src/EHFinder_Init.F90 index 2bf3d7d..8fc9587 100644 --- a/src/EHFinder_Init.F90 +++ b/src/EHFinder_Init.F90 @@ -15,7 +15,7 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS - CCTK_INT :: i, j, k, status + CCTK_INT :: i, j, k, l, status CCTK_REAL, dimension(3) :: xp, xpt CCTK_REAL, dimension(3,3) :: txyz CCTK_REAL :: cosa, sina, cosb, sinb, cosc, sinc @@ -35,13 +35,25 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS) cctk_time = last_time + ! Allocate the logical array containing the flag determining if the + ! corresponding level set should be re-initialized. + + if ( allocated(reparam_this_level_set) ) then + deallocate ( reparam_this_level_set ) + end if + allocate ( reparam_this_level_set(eh_number_level_sets) ) + if ( allocated(reparam_undone) ) then + deallocate ( reparam_undone ) + end if + allocate ( reparam_undone(eh_number_level_sets) ) + if ( evolve_generators .gt. 0 ) then - call CCTK_GrouplbndGN ( status, cctkGH, 1, lbnd, "ehfinder::generators" ) + call CCTK_GrouplbndGN ( status, cctkGH, 1, lbnd, "ehfinder::xg" ) if ( status .lt. 0 ) then call CCTK_WARN ( 0, "cannot get lower bounds for generator arrays" ) end if - call CCTK_GrouplshGN ( status, cctkGH, 1, lsh, "ehfinder::generators" ) + call CCTK_GrouplshGN ( status, cctkGH, 1, lsh, "ehfinder::xg" ) if ( status .lt. 0 ) then call CCTK_WARN ( 0, "cannot get local size for generator arrays" ) end if @@ -75,90 +87,93 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS) end if end if - ! If a sphere is requested... - if ( CCTK_EQUALS( initial_f, 'sphere' ) ) then - - ! Set up a sphere of radius initial_rad and translated - ! (translate_x,translate_y,translate_z) away from the origin. - f = sqrt( ( x - translate_x )**2 + & - ( y - translate_y )**2 + & - ( z - translate_z )**2 ) - initial_rad - if ( evolve_generators .gt. 0 ) then - if ( CCTK_EQUALS( generator_distribution, 'line' ) ) then - do i = 1, lsh(1) - theta = thetamin + dtheta * ( i + lbnd(1) - 1 ) - xg(i) = initial_rad * sin(theta) + translate_x - yg(i) = translate_y - zg(i) = initial_rad * cos(theta) + translate_z - end do + do l = 1, eh_number_level_sets + + ! If a sphere is requested... + if ( CCTK_EQUALS( initial_f(l), 'sphere' ) ) then + + ! Set up a sphere of radius initial_rad and translated + ! (translate_x,translate_y,translate_z) away from the origin. + f(:,:,:,l) = sqrt( ( x - translate_x )**2 + & + ( y - translate_y )**2 + & + ( z - translate_z )**2 ) - initial_rad(l) + if ( evolve_generators .gt. 0 ) then + if ( CCTK_EQUALS( generator_distribution, 'line' ) ) then + do i = 1, lsh(1) + theta = thetamin + dtheta * ( i + lbnd(1) - 1 ) + xg(i,l) = initial_rad(l) * sin(theta) + translate_x + yg(i,l) = translate_y + zg(i,l) = initial_rad(l) * cos(theta) + translate_z + end do + end if end if end if - end if - ! If an ellipsoid is requested... - if ( CCTK_EQUALS( initial_f, 'ellipsoid' ) ) then - - ! Calculate sines and cosines of the rotation parameters. - cosa = cos(rotation_alpha) - sina = sin(rotation_alpha) - cosb = cos(rotation_beta) - sinb = sin(rotation_beta) - cosc = cos(rotation_gamma) - sinc = sin(rotation_gamma) - - ! Set up the rotation matrix. The order is alpha around the z-axis, - ! beta around the y-axis and finally gamma around the x-axis. - txyz(1,1) = cosa * cosb - txyz(1,2) = sina * cosb - txyz(1,3) = -sinb - txyz(2,1) = cosa * sinb * sinc - sina * cosc - txyz(2,2) = sina * sinb * sinc + cosa * cosc - txyz(2,3) = cosb * sinc - txyz(3,1) = cosa * sinb * cosc + sina * sinc - txyz(3,2) = sina * sinb * cosc - cosa * sinc - txyz(3,3) = cosb * cosc -! print*,txyz + ! If an ellipsoid is requested... + if ( CCTK_EQUALS( initial_f(l), 'ellipsoid' ) ) then + + ! Calculate sines and cosines of the rotation parameters. + cosa = cos(rotation_alpha) + sina = sin(rotation_alpha) + cosb = cos(rotation_beta) + sinb = sin(rotation_beta) + cosc = cos(rotation_gamma) + sinc = sin(rotation_gamma) + + ! Set up the rotation matrix. The order is alpha around the z-axis, + ! beta around the y-axis and finally gamma around the x-axis. + txyz(1,1) = cosa * cosb + txyz(1,2) = sina * cosb + txyz(1,3) = -sinb + txyz(2,1) = cosa * sinb * sinc - sina * cosc + txyz(2,2) = sina * sinb * sinc + cosa * cosc + txyz(2,3) = cosb * sinc + txyz(3,1) = cosa * sinb * cosc + sina * sinc + txyz(3,2) = sina * sinb * cosc - cosa * sinc + txyz(3,3) = cosb * cosc +! print*,txyz ! Apply the rotations and translation for all points on the grid. ! Even though at first glance it looks like the translation is done ! first, the opposite is actually true. - do k = 1, nz - do j = 1, ny - do i = 1, nx - xp(1) = x(i,j,k) - translate_x - xp(2) = y(i,j,k) - translate_y - xp(3) = z(i,j,k) - translate_z - xpt = matmul ( txyz, xp ) - f(i,j,k) = sqrt( xpt(1)**2 / initial_a**2 + & - xpt(2)**2 / initial_b**2 + & - xpt(3)**2 / initial_c**2) - 1.0 + do k = 1, nz + do j = 1, ny + do i = 1, nx + xp(1) = x(i,j,k) - translate_x + xp(2) = y(i,j,k) - translate_y + xp(3) = z(i,j,k) - translate_z + xpt = matmul ( txyz, xp ) + f(i,j,k,l) = sqrt( xpt(1)**2 / initial_a**2 + & + xpt(2)**2 / initial_b**2 + & + xpt(3)**2 / initial_c**2) - 1.0 + end do end do end do - end do - - if ( evolve_generators .gt. 0 ) then - if ( CCTK_EQUALS( generator_distribution, 'line' ) ) then - do i = 1, lsh(1) - theta = thetamin + dtheta * ( i + lbnd(1) - 1 ) - r_el = sqrt ( one / ( sin(theta)**2 / initial_a**2 + & - cos(theta)**2 / initial_c**2 ) ) - xp(1) = r_el * sin(theta) + translate_x - xp(2) = translate_y - xp(3) = r_el * cos(theta) + translate_z - xpt = matmul ( txyz, xp ) - xg(i) = xpt(1) - yg(i) = xpt(2) - zg(i) = xpt(3) - end do + + if ( evolve_generators .gt. 0 ) then + if ( CCTK_EQUALS( generator_distribution, 'line' ) ) then + do i = 1, lsh(1) + theta = thetamin + dtheta * ( i + lbnd(1) - 1 ) + r_el = sqrt ( one / ( sin(theta)**2 / initial_a**2 + & + cos(theta)**2 / initial_c**2 ) ) + xp(1) = r_el * sin(theta) + translate_x + xp(2) = translate_y + xp(3) = r_el * cos(theta) + translate_z + xpt = matmul ( txyz, xp ) + xg(i,l) = xpt(1) + yg(i,l) = xpt(2) + zg(i,l) = xpt(3) + end do + end if end if end if - end if - ! if an ovaloid of Cassini is requested... - if ( CCTK_EQUALS( initial_f, 'cassini' ) ) then - f = (x**2+y**2+z**2)**2 + cas_a**4 - & - 2*cas_a**2*(x**2 - (y**2+z**2)) - cas_b**4 - end if + ! if an ovaloid of Cassini is requested... + if ( CCTK_EQUALS( initial_f, 'cassini' ) ) then + f(:,:,:,l) = (x**2+y**2+z**2)**2 + cas_a**4 - & + 2*cas_a**2*(x**2 - (y**2+z**2)) - cas_b**4 + end if + end do ! Initialise the internal mask. eh_mask = 0 diff --git a/src/EHFinder_Integrate2.F90 b/src/EHFinder_Integrate2.F90 index 6433533..d822050 100644 --- a/src/EHFinder_Integrate2.F90 +++ b/src/EHFinder_Integrate2.F90 @@ -15,7 +15,7 @@ subroutine EHFinder_FindSurfaceElement(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS - CCTK_INT :: i, j, k, im, jm + CCTK_INT :: i, j, k, l, im, jm CCTK_INT :: interp_handle, table_handle, coord_system_handle character(len=200) :: area_interp @@ -44,32 +44,29 @@ subroutine EHFinder_FindSurfaceElement(CCTK_ARGUMENTS) return endif + l = levelset_counter + call CCTK_INFO ( 'Finding surface element' ) call CCTK_GroupbboxGN ( ierr, cctkGH, 4, bbox, "ehfinder::surface_arrays" ) if ( ierr .lt. 0 ) then call CCTK_WARN ( 0, "cannot get bounding box for surface arrays" ) end if -! print*,bbox call CCTK_GroupgshGN ( ierr, cctkGH, 2, gsh, "ehfinder::surface_arrays" ) if ( ierr .lt. 0 ) then call CCTK_WARN ( 0, "cannot get global size for surface arrays" ) end if -! print*,gsh call CCTK_GrouplbndGN ( ierr, cctkGH, 2, lbnd, "ehfinder::surface_arrays" ) if ( ierr .lt. 0 ) then call CCTK_WARN ( 0, "cannot get lower bounds for surface arrays" ) end if -! print*,lbnd call CCTK_GroupubndGN ( ierr, cctkGH, 2, ubnd, "ehfinder::surface_arrays" ) if ( ierr .lt. 0 ) then call CCTK_WARN ( 0, "cannot get upper bounds for surface arrays" ) end if -! print*,ubnd call CCTK_GrouplshGN ( ierr, cctkGH, 2, lsh, "ehfinder::surface_arrays" ) if ( ierr .lt. 0 ) then call CCTK_WARN ( 0, "cannot get local size for surface arrays" ) end if -! print*,lsh call CCTK_GroupnghostzonesGN ( ierr, cctkGH, 2, nghost, "ehfinder::surface_arrays" ) if ( ierr .lt. 0 ) then call CCTK_WARN ( 0, "cannot get local size for surface arrays" ) @@ -115,7 +112,6 @@ subroutine EHFinder_FindSurfaceElement(CCTK_ARGUMENTS) ! find the cartesian coordinates for the interpolation points -! print*,center do j = 1, lsh(2) do i = 1, lsh(1) interp_x(i,j) = center(1) + rsurf(i,j) * sintheta(i,j) * cosphi(i,j) @@ -155,24 +151,16 @@ subroutine EHFinder_FindSurfaceElement(CCTK_ARGUMENTS) do j = 1, lsh(2) do i = 1, lsh(1) -! print*,i,j,ctheta(i,j),cphi(i,j) dxdth = ( drdtheta(i,j) * sintheta(i,j) + & rsurf(i,j) * costheta(i,j) ) * cosphi(i,j) -! print*,dxdth,rsurf(i,j)*costheta(i,j)*cosphi(i,j) dydth = ( drdtheta(i,j) * sintheta(i,j) + & rsurf(i,j) * costheta(i,j) ) * sinphi(i,j) -! print*,dydth,rsurf(i,j)*costheta(i,j)*sinphi(i,j) dzdth = drdtheta(i,j) * costheta(i,j) - rsurf(i,j) * sintheta(i,j) -! print*,dzdth,-rsurf(i,j)*sintheta(i,j) dxdph = ( drdphi(i,j) * cosphi(i,j) - & rsurf(i,j) * sinphi(i,j) ) * sintheta(i,j) -! print*,dxdph,-rsurf(i,j)*sintheta(i,j)*sinphi(i,j) dydph = ( drdphi(i,j) * sinphi(i,j) + & rsurf(i,j) * cosphi(i,j) ) * sintheta(i,j) -! print*,dydph,rsurf(i,j)*sintheta(i,j)*cosphi(i,j) dzdph = drdphi(i,j) * costheta(i,j) -! print*,dzdph,0 -! print* gtt(i,j) = dxdth**2 * gxxi(i,j) + dydth**2 * gyyi(i,j) + & dzdth**2 * gzzi(i,j) + & @@ -195,14 +183,6 @@ subroutine EHFinder_FindSurfaceElement(CCTK_ARGUMENTS) dlphi(i,j) = dphi * sqrt ( gpp(i,j) ) end do end do -! print*,da -! do j = 1, lsh(2) -! print*,'phi = ',cphi(1,j) -! print*,da(:,j) -! print*,gtt(:,j) -! print*,gtp(:,j) -! print*,gpp(:,j) -! end do end subroutine EHFinder_FindSurfaceElement @@ -217,7 +197,7 @@ subroutine EHFinder_IntegrateArea(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS - CCTK_INT :: int_var, sn + CCTK_INT :: int_var, sn, l CCTK_REAL :: eh_area_tmp character(len=30) :: info_message CCTK_INT, dimension(1) :: lbnd, ubnd, lsh @@ -227,6 +207,8 @@ subroutine EHFinder_IntegrateArea(CCTK_ARGUMENTS) return endif + l = levelset_counter + call CCTK_INFO ( 'Integrating area' ) sn = surface_counter - integrate_counter @@ -236,9 +218,6 @@ subroutine EHFinder_IntegrateArea(CCTK_ARGUMENTS) if ( int_var .lt. 0 ) then call CCTK_WARN ( 0, 'Could not get index to grid array int_tmp' ) end if -! call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, & -! eh_area2(sn), 1, int_var ) -! write(info_message,'(a15,f10.5)') 'Horizon area = ',eh_area2(sn) call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, & eh_area_tmp, 1, int_var ) @@ -259,12 +238,9 @@ subroutine EHFinder_IntegrateArea(CCTK_ARGUMENTS) call CCTK_WARN ( 0, "cannot get local size for area array" ) end if -! if ( ( sn .ge. lbnd(1) + lsh(1) ) .and. ( sn .le. ubnd(1) + lsh(1) ) ) then if ( ( sn .ge. lbnd(1) + 1 ) .and. ( sn .le. ubnd(1) + 1 ) ) then - eh_area2(sn-lbnd(1)) = eh_area_tmp + eh_area2(sn-lbnd(1),l) = eh_area_tmp end if - print*,'Debug information : ',sn,lbnd,ubnd,lsh -! print*,eh_area2 end subroutine EHFinder_IntegrateArea @@ -279,7 +255,7 @@ subroutine EHFinder_IntegrateCentroid(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS - CCTK_INT :: int_var, sn, k + CCTK_INT :: int_var, sn, k, l character(len=50) :: info_message CCTK_INT, dimension(1) :: lbnd, ubnd, lsh CCTK_REAL :: centroid_x, centroid_y, centroid_z @@ -289,6 +265,8 @@ subroutine EHFinder_IntegrateCentroid(CCTK_ARGUMENTS) return endif + l = levelset_counter + call CCTK_INFO ( 'Integrating centroid' ) sn = surface_counter - integrate_counter @@ -302,22 +280,6 @@ subroutine EHFinder_IntegrateCentroid(CCTK_ARGUMENTS) call CCTK_WARN ( 0, 'Could not get index to grid array int_tmp' ) end if -! int_tmp = sym_factor * interp_x * weights * da -! call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, & -! eh_centroid2_x(sn), 1, int_var ) -! -! int_tmp = sym_factor * interp_y * weights * da -! call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, & -! eh_centroid2_y(sn), 1, int_var ) -! -! int_tmp = sym_factor * interp_z * weights * da -! call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, & -! eh_centroid2_z(sn), 1, int_var ) - -! if ( s_symx ) eh_centroid2_x(sn) = zero -! if ( s_symy ) eh_centroid2_y(sn) = zero -! if ( s_symz ) eh_centroid2_z(sn) = zero - int_tmp = sym_factor * interp_x * weights * da call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, & centroid_x, 1, int_var ) @@ -334,41 +296,30 @@ subroutine EHFinder_IntegrateCentroid(CCTK_ARGUMENTS) if ( s_symy ) centroid_y = zero if ( s_symz ) centroid_z = zero - call CCTK_GrouplbndGN ( ierr, cctkGH, 1, lbnd, "ehfinder::eh_centroid2" ) + call CCTK_GrouplbndGN ( ierr, cctkGH, 1, lbnd, "ehfinder::eh_centroid2_x" ) if ( ierr .lt. 0 ) then call CCTK_WARN ( 0, "cannot get lower bounds for centroid array" ) end if - call CCTK_GroupubndGN ( ierr, cctkGH, 1, ubnd, "ehfinder::eh_centroid2" ) + call CCTK_GroupubndGN ( ierr, cctkGH, 1, ubnd, "ehfinder::eh_centroid2_x" ) if ( ierr .lt. 0 ) then call CCTK_WARN ( 0, "cannot get upper bounds for centroid array" ) end if - call CCTK_GrouplshGN ( ierr, cctkGH, 1, lsh, "ehfinder::eh_centroid2" ) + call CCTK_GrouplshGN ( ierr, cctkGH, 1, lsh, "ehfinder::eh_centroid2_x" ) if ( ierr .lt. 0 ) then call CCTK_WARN ( 0, "cannot get local size for centroid array" ) end if -! eh_centroid2_x(sn) = eh_centroid2_x(sn) / eh_area2(sn) -! eh_centroid2_y(sn) = eh_centroid2_y(sn) / eh_area2(sn) -! eh_centroid2_z(sn) = eh_centroid2_z(sn) / eh_area2(sn) - if ( ( sn .ge. lbnd(1) + 1 ) .and. ( sn .le. ubnd(1) + 1 ) ) then k = sn - lbnd(1) - eh_centroid2_x(k) = centroid_x / eh_area2(k) - eh_centroid2_y(k) = centroid_y / eh_area2(k) - eh_centroid2_z(k) = centroid_z / eh_area2(k) -! print*,centroid_x,centroid_y,centroid_z -! print*,eh_area2(k) + eh_centroid2_x(k,l) = centroid_x / eh_area2(k,l) + eh_centroid2_y(k,l) = centroid_y / eh_area2(k,l) + eh_centroid2_z(k,l) = centroid_z / eh_area2(k,l) end if -! write(info_message,'(a19,3(f10.5))') 'Horizon centroid = ', & -! eh_centroid2_x(sn), & -! eh_centroid2_y(sn), & -! eh_centroid2_z(sn) - write(info_message,'(a19,3(f10.5))') 'Horizon centroid = ', & - eh_centroid2_x(sn), & - eh_centroid2_y(sn), & - eh_centroid2_z(sn) + eh_centroid2_x(sn,l), & + eh_centroid2_y(sn,l), & + eh_centroid2_z(sn,l) call CCTK_INFO(info_message) end subroutine EHFinder_IntegrateCentroid @@ -384,7 +335,7 @@ subroutine EHFinder_IntegrateEquatorial(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS - CCTK_INT :: int_var, sn, k + CCTK_INT :: int_var, sn, k, l character(len=50) :: info_message CCTK_INT, dimension(1) :: lbnd1, ubnd1, lsh1 CCTK_INT, dimension(2) :: lbnd, lsh, nghost @@ -397,6 +348,8 @@ subroutine EHFinder_IntegrateEquatorial(CCTK_ARGUMENTS) return endif + l = levelset_counter + sn = surface_counter - integrate_counter ! Find out the lower bounds of the distributed integration grid arrays. @@ -423,11 +376,6 @@ subroutine EHFinder_IntegrateEquatorial(CCTK_ARGUMENTS) call CCTK_WARN ( 0, "cannot get ghost zones for surface arrays" ) end if -! print*,lbnd -! print*,lsh -! print*,bbox -! print*,nghost - ithl = 1; ithr = lsh(1) jphl = 1; jphr = lsh(2) @@ -436,9 +384,6 @@ subroutine EHFinder_IntegrateEquatorial(CCTK_ARGUMENTS) if ( bbox(3) .eq. 0 ) jphl = jphl + nghost(2) if ( bbox(4) .eq. 0 ) jphr = jphr - nghost(2) -! print*,ithl,ithr -! print*,jphl,jphr - call CCTK_INFO ( 'Integrating equatorial circumference' ) int_tmp = phi_sym_factor * phiweights * dlphi @@ -450,46 +395,35 @@ subroutine EHFinder_IntegrateEquatorial(CCTK_ARGUMENTS) eq_circ_loc = zero end if -! print*,itha,githeta,eq_circ_loc -! print* -! print*,int_tmp(itha,jphl:jphr) -! print* -! print*,phi_sym_factor -! print* -! print*,phiweights(itha,jphl:jphr) -! print* -! print*,dlphi(itha,jphl:jphr) - call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, sum_handle, & eq_circ_loc, eq_circ, CCTK_VARIABLE_REAL ) - call CCTK_GrouplbndGN ( ierr, cctkGH, 1, lbnd1, "ehfinder::eh_circumference2" ) + call CCTK_GrouplbndGN ( ierr, cctkGH, 1, lbnd1, "ehfinder::eh_circ_eq2" ) if ( ierr .lt. 0 ) then call CCTK_WARN ( 0, "cannot get lower bounds for area array" ) end if - call CCTK_GroupubndGN ( ierr, cctkGH, 1, ubnd1, "ehfinder::eh_circumference2" ) + call CCTK_GroupubndGN ( ierr, cctkGH, 1, ubnd1, "ehfinder::eh_circ_eq2" ) if ( ierr .lt. 0 ) then call CCTK_WARN ( 0, "cannot get upper bounds for area array" ) end if - call CCTK_GrouplshGN ( ierr, cctkGH, 1, lsh1, "ehfinder::eh_circumference2" ) + call CCTK_GrouplshGN ( ierr, cctkGH, 1, lsh1, "ehfinder::eh_circ_eq2" ) if ( ierr .lt. 0 ) then call CCTK_WARN ( 0, "cannot get local size for area array" ) end if if ( ( sn .ge. lbnd1(1) + 1 ) .and. ( sn .le. ubnd1(1) + 1 ) ) then k = sn - lbnd1(1) - eh_circ_eq2(k) = eq_circ + eh_circ_eq2(k,l) = eq_circ end if -! eh_circ_eq2(sn) = eq_circ - write(info_message,'(a35,f10.5)') 'Horizon equatorial circumference = ',eq_circ + write(info_message,'(a35,f10.5)') 'Horizon equatorial circumference = ', & + eq_circ call CCTK_INFO(info_message) call CCTK_INFO ( 'Integrating polar circumference' ) int_tmp = theta_sym_factor * thetaweights * dltheta -! print*,jphl+lbnd(2),jphr+lbnd(2),gjphi if ( ( jphl+lbnd(2) .le. gjphi ) .and. ( gjphi .le. jphr+lbnd(2) ) ) then jpha = gjphi - lbnd(2) pol_circ_loc = sum ( int_tmp(ithl:ithr,jpha) ) @@ -497,23 +431,14 @@ subroutine EHFinder_IntegrateEquatorial(CCTK_ARGUMENTS) pol_circ_loc = zero end if -! print*,jpha,gjphi,pol_circ_loc -! print* -! print*,int_tmp(ithl:ithr,jpha) -! print* -! print*,theta_sym_factor -! print*,thetaweights(ithl:ithr,jpha) -! print* -! print*,dltheta(ithl:ithr,jpha) call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, sum_handle, & pol_circ_loc, pol_circ, CCTK_VARIABLE_REAL ) if ( ( sn .ge. lbnd1(1) + 1 ) .and. ( sn .le. ubnd1(1) + 1 ) ) then k = sn - lbnd1(1) - eh_circ_pol2(k) = pol_circ + eh_circ_pol2(k,l) = pol_circ end if -! eh_circ_pol2(sn) = pol_circ write(info_message,'(a30,f10.5)') 'Horizon polar circumference = ',pol_circ call CCTK_INFO(info_message) @@ -611,11 +536,7 @@ subroutine EHFinder_CopyArea(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS -! print*,Xeh_area0 -! print*,Xeh_area20 - ! If finding of surface failed set area to zero. - print*,'Debug2 : ',find_surface_status if ( find_surface_status .lt. 0 ) then eh_area = zero return @@ -636,9 +557,6 @@ subroutine EHFinder_CopyCentroid(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS -! print*,Xeh_centroid0 -! print*,Xeh_centroid20 - ! If finding of surface failed set centroids to zero. if ( find_surface_status .lt. 0 ) then eh_centroid_x = zero diff --git a/src/EHFinder_ReParametrize.F90 b/src/EHFinder_ReParametrize.F90 index 768a5e1..7559395 100644 --- a/src/EHFinder_ReParametrize.F90 +++ b/src/EHFinder_ReParametrize.F90 @@ -17,8 +17,8 @@ subroutine EHFinder_ReParametrizeControl(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS - CCTK_REAL :: fmin, fmax, fminloc, fmaxloc - CCTK_INT :: i, j, k + CCTK_REAL, dimension(eh_number_level_sets) :: fmin, fmax, fminloc, fmaxloc + CCTK_INT :: i, j, k, l ! print*,'EHFinder_ReParametrize1 Entered' ! Initialize the control variables. The default means no re-parametrization. @@ -88,27 +88,41 @@ subroutine EHFinder_ReParametrizeControl(CCTK_ARGUMENTS) ftmp = f fbak = f eh_mask_bak = eh_mask - fminloc = minval(f) - fmaxloc = maxval(f) - rep_mask = 0 +! rep_mask = 0 + do l = 1, eh_number_level_sets + fminloc(l) = minval(f(:,:,:,l)) + fmaxloc(l) = maxval(f(:,:,:,l)) + end do - call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, max_handle, & - fmaxloc, fmax, CCTK_VARIABLE_REAL ) + call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, max_handle, & + fmaxloc, fmax, eh_number_level_sets, & + CCTK_VARIABLE_REAL ) if ( ierr .ne. 0 ) then call CCTK_WARN(0,'Reduction of fmax failed') end if - call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, min_handle, & - fminloc, fmin, CCTK_VARIABLE_REAL ) + call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, min_handle, & + fminloc, fmin, eh_number_level_sets, & + CCTK_VARIABLE_REAL ) if ( ierr .ne. 0 ) then call CCTK_WARN(0,'Reduction of fmin failed') end if + reparam_this_level_set = .true. + ! If fmin and fmax have the same sign, there is no surface present and -! re-parametrization will not be performed. - if ( fmin * fmax .gt. zero ) then +! re-parametrization will not be performed for the given surface. + do l = 1, eh_number_level_sets + if ( fmin(l) * fmax(l) .gt. zero ) then + reparam_this_level_set(l) = .false. + end if + end do + +! If none of the surfaces should be reparametrized set the control variables +! to zero. + if ( all ( .not. reparam_this_level_set ) ) then re_param_control_approx = 0 re_param_control_pde = 0 - call CCTK_INFO ( 'No zero-point of the level set function. No re-parametrization performed.' ) + call CCTK_INFO ( 'No zero-points of the level set functions. No re-parametrization performed.' ) end if end if ! print*,'EHFinder_ReParametrizeControl Exited' @@ -116,151 +130,151 @@ subroutine EHFinder_ReParametrizeControl(CCTK_ARGUMENTS) end subroutine EHFinder_ReParametrizeControl -subroutine EHFinder_ReParametrizeRK2_1(CCTK_ARGUMENTS) - - use EHFinder_mod - - implicit none - - DECLARE_CCTK_PARAMETERS - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_FUNCTIONS - - CCTK_INT :: i, j, k - CCTK_REAL :: idx, idy, idz - CCTK_REAL :: al, ar, bl, br, cl, cr - CCTK_REAL :: alminus, alplus, blminus, blplus, clminus, clplus - CCTK_REAL :: arminus, arplus, brminus, brplus, crminus, crplus - -! print*,'EHFinder_ReParametrizeRK2_1 Entered' - h = hfac*minval(cctk_delta_space) - - if ( CCTK_EQUALS ( pde_differences, 'centered' ) ) then - -# include "include/centered_second.h" - - end if - - if ( CCTK_EQUALS ( pde_differences, 'upwind' ) ) then - -# include "include/upwind_first.h" - - end if - - if ( CCTK_EQUALS ( pde_differences, 'upwind2' ) ) then - -# include "include/upwind_second.h" - - end if - - where ( eh_mask .ge. 0 ) - sftmp = sqrt( dfx**2 + dfy**2 + dfz**2 ) - sftmp = - half * h * f / sqrt(f**2+one) * ( sftmp - one ) - f = f + sftmp - elsewhere - sftmp = one - end where -! print*,'EHFinder_ReParametrizeRK2_1 Exited' - -end subroutine EHFinder_ReParametrizeRK2_1 - - -subroutine EHFinder_ReParametrizeRK2_2(CCTK_ARGUMENTS) - - - use EHFinder_mod - - implicit none - - DECLARE_CCTK_PARAMETERS - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_FUNCTIONS - - CCTK_INT :: i, j, k - CCTK_REAL :: idx, idy, idz, maxf, maxdfloc, maxdf, mindfloc, mindf - CCTK_REAL :: al, ar, bl, br, cl, cr - CCTK_REAL :: alminus, alplus, blminus, blplus, clminus, clplus - CCTK_REAL :: arminus, arplus, brminus, brplus, crminus, crplus - character(len=128) :: info_message - -! print*,'EHFinder_ReParametrizeRK2_2 Entered' - - if ( CCTK_EQUALS ( pde_differences, 'centered' ) ) then - -# include "include/centered_second.h" - - end if - - if ( CCTK_EQUALS ( pde_differences, 'upwind' ) ) then - -# include "include/upwind_first.h" - - end if - - if ( CCTK_EQUALS ( pde_differences, 'upwind2' ) ) then - -# include "include/upwind_second.h" - - end if - - where ( eh_mask .ge. 0 ) - - sftmp = sqrt( dfx**2 + dfy**2 + dfz**2 ) - sftmp = - h * f / sqrt(f**2+one) * ( sftmp - one ) - f = ftmp + sftmp - elsewhere - sftmp = one - end where - - maxdfloc = zero - mindfloc = 1.0d23 - do k = kzl, kzr - do j = jyl, jyr - do i =ixl, ixr - if ( eh_mask(i,j,k) .ge. 0 ) then - maxdfloc = max ( maxdfloc, abs(sftmp(i,j,k)) ) - mindfloc = min ( mindfloc, abs(sftmp(i,j,k)) ) - end if - end do - end do - end do - - call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, max_handle, & - maxdfloc, maxdf, CCTK_VARIABLE_REAL ) - if ( ierr .ne. 0 ) then - call CCTK_WARN(0,'Reduction of maxdf failed') - end if - call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, min_handle, & - mindfloc, mindf, CCTK_VARIABLE_REAL ) - if ( ierr .ne. 0 ) then - call CCTK_WARN(0,'Reduction of mindf failed') - end if - - if ( maxdf .lt. h*minval(cctk_delta_space)**2 ) re_param_control_pde = 0 - - ncalls = ncalls + 1 - - if ( re_param_control_pde .eq. 0 ) then - write(info_message,'(a35,i5,a12)') & - 'PDE re-parametrization complete in ',ncalls,' iterations.' - call CCTK_INFO(info_message) - end if - - if ( ncalls .gt. re_param_max_iter ) then - call CCTK_WARN(1,'Re-parametrization failed to converge') - re_param_control_pde = 0 - end if - - ftmp = f -! print*,ncalls,maxdf,mindf -! print*,mx1,my1,mz2,mx2,my2,mz2 -! print*,f(mx1,my1,mz1),f(mx2,my2,mz2),eh_mask(mx1,my1,mz1),eh_mask(mx2,my2,mz2) -! print*,f(27,53,27),dfsq(27,53,27),ftmp(27,53,27) -! -! call CCTK_INFO('I was called') -! print*,'EHFinder_ReParametrizeRK2_2 Exited' - -end subroutine EHFinder_ReParametrizeRK2_2 +!subroutine EHFinder_ReParametrizeRK2_1(CCTK_ARGUMENTS) +! +! use EHFinder_mod +! +! implicit none +! +! DECLARE_CCTK_PARAMETERS +! DECLARE_CCTK_ARGUMENTS +! DECLARE_CCTK_FUNCTIONS +! +! CCTK_INT :: i, j, k +! CCTK_REAL :: idx, idy, idz +! CCTK_REAL :: al, ar, bl, br, cl, cr +! CCTK_REAL :: alminus, alplus, blminus, blplus, clminus, clplus +! CCTK_REAL :: arminus, arplus, brminus, brplus, crminus, crplus +! +!! print*,'EHFinder_ReParametrizeRK2_1 Entered' +! h = hfac*minval(cctk_delta_space) +! +! if ( CCTK_EQUALS ( pde_differences, 'centered' ) ) then +! +!# include "include/centered_second.h" +! +! end if +! +! if ( CCTK_EQUALS ( pde_differences, 'upwind' ) ) then +! +!# include "include/upwind_first.h" +! +! end if +! +! if ( CCTK_EQUALS ( pde_differences, 'upwind2' ) ) then +! +!# include "include/upwind_second.h" +! +! end if +! +! where ( eh_mask .ge. 0 ) +! sftmp = sqrt( dfx**2 + dfy**2 + dfz**2 ) +! sftmp = - half * h * f / sqrt(f**2+one) * ( sftmp - one ) +! f = f + sftmp +! elsewhere +! sftmp = one +! end where +!! print*,'EHFinder_ReParametrizeRK2_1 Exited' +! +!end subroutine EHFinder_ReParametrizeRK2_1 + + +!subroutine EHFinder_ReParametrizeRK2_2(CCTK_ARGUMENTS) +! +! +! use EHFinder_mod +! +! implicit none +! +! DECLARE_CCTK_PARAMETERS +! DECLARE_CCTK_ARGUMENTS +! DECLARE_CCTK_FUNCTIONS +! +! CCTK_INT :: i, j, k +! CCTK_REAL :: idx, idy, idz, maxf, maxdfloc, maxdf, mindfloc, mindf +! CCTK_REAL :: al, ar, bl, br, cl, cr +! CCTK_REAL :: alminus, alplus, blminus, blplus, clminus, clplus +! CCTK_REAL :: arminus, arplus, brminus, brplus, crminus, crplus +! character(len=128) :: info_message +! +!! print*,'EHFinder_ReParametrizeRK2_2 Entered' +! +! if ( CCTK_EQUALS ( pde_differences, 'centered' ) ) then +! +!# include "include/centered_second.h" +! +! end if +! +! if ( CCTK_EQUALS ( pde_differences, 'upwind' ) ) then +! +!# include "include/upwind_first.h" +! +! end if +! +! if ( CCTK_EQUALS ( pde_differences, 'upwind2' ) ) then +! +!# include "include/upwind_second.h" +! +! end if +! +! where ( eh_mask .ge. 0 ) +! +! sftmp = sqrt( dfx**2 + dfy**2 + dfz**2 ) +! sftmp = - h * f / sqrt(f**2+one) * ( sftmp - one ) +! f = ftmp + sftmp +! elsewhere +! sftmp = one +! end where +! +! maxdfloc = zero +! mindfloc = 1.0d23 +! do k = kzl, kzr +! do j = jyl, jyr +! do i =ixl, ixr +! if ( eh_mask(i,j,k) .ge. 0 ) then +! maxdfloc = max ( maxdfloc, abs(sftmp(i,j,k)) ) +! mindfloc = min ( mindfloc, abs(sftmp(i,j,k)) ) +! end if +! end do +! end do +! end do +! +! call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, max_handle, & +! maxdfloc, maxdf, CCTK_VARIABLE_REAL ) +! if ( ierr .ne. 0 ) then +! call CCTK_WARN(0,'Reduction of maxdf failed') +! end if +! call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, min_handle, & +! mindfloc, mindf, CCTK_VARIABLE_REAL ) +! if ( ierr .ne. 0 ) then +! call CCTK_WARN(0,'Reduction of mindf failed') +! end if +! +! if ( maxdf .lt. h*minval(cctk_delta_space)**2 ) re_param_control_pde = 0 +! +! ncalls = ncalls + 1 +! +! if ( re_param_control_pde .eq. 0 ) then +! write(info_message,'(a35,i5,a12)') & +! 'PDE re-parametrization complete in ',ncalls,' iterations.' +! call CCTK_INFO(info_message) +! end if +! +! if ( ncalls .gt. re_param_max_iter ) then +! call CCTK_WARN(1,'Re-parametrization failed to converge') +! re_param_control_pde = 0 +! end if +! +! ftmp = f +!! print*,ncalls,maxdf,mindf +!! print*,mx1,my1,mz2,mx2,my2,mz2 +!! print*,f(mx1,my1,mz1),f(mx2,my2,mz2),eh_mask(mx1,my1,mz1),eh_mask(mx2,my2,mz2) +!! print*,f(27,53,27),dfsq(27,53,27),ftmp(27,53,27) +!! +!! call CCTK_INFO('I was called') +!! print*,'EHFinder_ReParametrizeRK2_2 Exited' +! +!end subroutine EHFinder_ReParametrizeRK2_2 subroutine EHFinder_ReParametrizeEuler(CCTK_ARGUMENTS) @@ -273,8 +287,10 @@ subroutine EHFinder_ReParametrizeEuler(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS - CCTK_INT :: i, j, k - CCTK_REAL :: idx, idy, idz, maxf, maxdfloc, maxdf, mindfloc, mindf + CCTK_INT :: i, j, k, l + CCTK_REAL :: idx, idy, idz, maxf + CCTK_REAL, dimension(eh_number_level_sets) :: maxdfloc, maxdf, & + mindfloc, mindf CCTK_REAL :: al, ar, bl, br, cl, cr CCTK_REAL :: alminus, alplus, blminus, blplus, clminus, clplus CCTK_REAL :: arminus, arplus, brminus, brplus, crminus, crplus @@ -303,43 +319,53 @@ subroutine EHFinder_ReParametrizeEuler(CCTK_ARGUMENTS) where ( eh_mask .ge. 0 ) sftmp = sqrt( dfx**2 + dfy**2 + dfz**2 ) -! where ( f .gt. rdel ) -! sftmp = - h * ( sftmp - one ) * & -! ( rdel + ( f - rdel ) / sqrt ( one + ( f - rdel )**2 ) ) -! elsewhere -! sftmp = - h * f * ( sftmp - one ) -! end where sftmp = - h * f / sqrt(f**2+one) * ( sftmp - one ) f = f + sftmp elsewhere sftmp = one end where - - maxdfloc = zero - mindfloc = 1.0d23 - do k = kzl, kzr - do j = jyl, jyr - do i = ixl, ixr - if ( eh_mask(i,j,k) .ge. 0 ) then - maxdfloc = max ( maxdfloc, abs(sftmp(i,j,k)) ) - mindfloc = min ( mindfloc, abs(sftmp(i,j,k)) ) - end if - end do - end do +! print*,f(2:4,2:4,2:4,5) +! print*,sftmp(3,3,3,5) +! stop + +! do l = 1, eh_number_level_sets +! maxdfloc = zero +! mindfloc = 1.0d23 +! do k = kzl, kzr +! do j = jyl, jyr +! do i = ixl, ixr +! if ( eh_mask(i,j,k,l) .ge. 0 ) then +! maxdfloc = max ( maxdfloc, abs(sftmp(i,j,k,l)) ) +! mindfloc = min ( mindfloc, abs(sftmp(i,j,k,l)) ) +! end if +! end do +! end do +! end do +! end do + + do l = 1, eh_number_level_sets + maxdfloc(l) = maxval ( abs(sftmp(ixl:ixr,jyl:jyr,kzl:kzr,l)), & + mask = eh_mask(ixl:ixr,jyl:jyr,kzl:kzr,l) .ge. 0) + mindfloc(l) = minval ( abs(sftmp(ixl:ixr,jyl:jyr,kzl:kzr,l)), & + mask = eh_mask(ixl:ixr,jyl:jyr,kzl:kzr,l) .ge. 0) end do - call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, max_handle, & - maxdfloc, maxdf, CCTK_VARIABLE_REAL ) + call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, max_handle, & + maxdfloc, maxdf, eh_number_level_sets, & + CCTK_VARIABLE_REAL ) if ( ierr .ne. 0 ) then call CCTK_WARN(0,'Reduction of maxdf failed') end if - call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, min_handle, & - mindfloc, mindf, CCTK_VARIABLE_REAL ) + call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, min_handle, & + mindfloc, mindf, eh_number_level_sets, & + CCTK_VARIABLE_REAL ) if ( ierr .ne. 0 ) then call CCTK_WARN(0,'Reduction of mindf failed') end if - if ( maxdf .lt. h*minval(cctk_delta_space)**2 ) re_param_control_pde = 0 + if ( all ( maxdf .lt. h*minval(cctk_delta_space)**2 ) ) then + re_param_control_pde = 0 + endif ncalls = ncalls + 1 @@ -354,264 +380,266 @@ subroutine EHFinder_ReParametrizeEuler(CCTK_ARGUMENTS) re_param_control_pde = 0 end if -! print*,ncalls,maxdf,mindf +! print*,ncalls +! print*,maxdf +! print*,mindf ! print*,'EHFinder_ReParametrizeEuler Exited' end subroutine EHFinder_ReParametrizeEuler -subroutine EHFinder_ReParametrize5(CCTK_ARGUMENTS) - - use EHFinder_mod - - implicit none - - DECLARE_CCTK_PARAMETERS - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_FUNCTIONS - - CCTK_INT :: i, j, k - CCTK_INT :: rep_countloc, rep_totalloc - CCTK_REAL :: idx2, idy2, idz2, idx, idy, idz - CCTK_REAL :: d2fdx2, d2fdy2, d2fdz2, d2fdxdy, d2fdxdz, d2fdydz - CCTK_REAL :: n_x, n_y, n_z - CCTK_REAL :: a, b, c, lambda - - interface - CCTK_REAL function solve_quadratic ( a, b, c ) - CCTK_REAL, intent(in) :: a, b, c - end function solve_quadratic - end interface - - if ( re_param_control_approx .eq. 1 ) then - if ( reparametrize_every_approx .gt. 0 ) then - if ( mod(cctk_iteration,reparametrize_every_approx) .eq. 0 ) then - call CCTK_ReductionArrayHandle ( sum_handle, 'sum' ) - if ( sum_handle .lt. 0 ) then - call CCTK_WARN(0,'Could not obtain a handle for sum reduction') - end if - -# include "include/centered_second.h" - - ncalls = 0 - - rep_totalloc = count ( eh_mask(ixl:ixr,jyl:jyr,kzl:kzr) .ge. 0 ) - - rep_countloc = 0 - - ftmp = f - - ixr = min ( ixr, nx-1 ) - jyr = min ( jyr, ny-1 ) - kzr = min ( kzr, nz-1 ) - - do k = kzl, kzr - do j = jyl, jyr - do i = ixl, ixr - if ( eh_mask(i,j,k) .eq. 0 ) then - if ( ( f(i,j,k)*f(i-1,j,k) .le. zero ) .or. & - ( f(i,j,k)*f(i+1,j,k) .le. zero ) .or. & - ( f(i,j,k)*f(i,j-1,k) .le. zero ) .or. & - ( f(i,j,k)*f(i,j+1,k) .le. zero ) .or. & - ( f(i,j,k)*f(i,j,k-1) .le. zero ) .or. & - ( f(i,j,k)*f(i,j,k+1) .le. zero ) ) then - - rep_mask(i,j,k) = 1 - rep_countloc = rep_countloc + 1 - - d2fdx2 = idx2 * ( f(i-1,j,k) - two * f(i,j,k) + f(i+1,j,k) ) - d2fdy2 = idy2 * ( f(i,j-1,k) - two * f(i,j,k) + f(i,j+1,k) ) - d2fdz2 = idz2 * ( f(i,j,k-1) - two * f(i,j,k) + f(i,j,k+1) ) - - d2fdxdy = idx * idy * ( f(i-1,j-1,k) - f(i+1,j-1,k) - & - f(i-1,j+1,k) + f(i+1,j+1,k) ) - d2fdxdz = idx * idz * ( f(i-1,j,k-1) - f(i+1,j,k-1) - & - f(i-1,j,k+1) + f(i+1,j,k+1) ) - d2fdydz = idy * idz * ( f(i,j-1,k-1) - f(i,j+1,k-1) - & - f(i,j-1,k+1) + f(i,j+1,k+1) ) - - c = f(i,j,k) - b = sqrt ( dfx(i,j,k)**2 + dfy(i,j,k)**2 + dfz(i,j,k)**2 ) - n_x = dfx(i,j,k) / b - n_y = dfy(i,j,k) / b - n_z = dfz(i,j,k) / b - a = half*( n_x**2*d2fdx2 + n_y**2*d2fdy2 + n_z**2*d2fdz2 ) +& - n_x*n_y*d2fdxdy + n_x*n_z*d2fdxdz + n_y*n_z*d2fdydz - - if ( a .ne. 0 ) then - lambda = solve_quadratic ( a, b, c ) - else - lambda = abs ( c / b ) - endif - if ( lambda .eq. huge ) then - lambda = abs ( c / b ) - end if - - if ( lambda .gt. two*delta ) then - print*,i,j,k - print*,f(i-1:i+1,j-1:j+1,k-1:k+1) - print*,eh_mask(i-1:i+1,j-1:j+1,k-1:k+1) - print*,a,b,c,lambda - call CCTK_WARN(0,'Re-parametrization failed') - end if - - ftmp(i,j,k) = sign ( lambda, f(i,j,k) ) -! print*,i,j,k -! print*,lambda,f(i,j,k), ftmp(i,j,k) -! print*,f(i-1,j,k),f(i+1,j,k) -! print*,f(i,j-1,k),f(i,j+1,k) -! print*,f(i,j,k-1),f(i,j,k+1) - end if - end if - end do - end do - end do - - where ( eh_mask .ge. 0 ) f = ftmp - -! print*,rep_countloc,rep_totalloc - call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, sum_handle, & - rep_countloc, rep_count, CCTK_VARIABLE_INT ) - if ( ierr .ne. 0 ) then - call CCTK_WARN(0,'Sum of rep_count failed') - end if - - call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, sum_handle, & - rep_totalloc, rep_total, CCTK_VARIABLE_INT ) - if ( ierr .ne. 0 ) then - call CCTK_WARN(0,'Sum of rep_total failed') - end if -! print*,'First step' -! print*,rep_count,rep_total -! call CCTK_WARN(0,'Finished') - -! call CartSymGN(ierr,cctkGH,'ehfinder::level_set') -! call CartSymGN(ierr,cctkGH,'ehfinder::rep_mask') - - end if - end if - - end if - -end subroutine EHFinder_ReParametrize5 - - -subroutine EHFinder_ReParametrize6(CCTK_ARGUMENTS) - - use EHFinder_mod - - implicit none - - DECLARE_CCTK_PARAMETERS - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_FUNCTIONS - - CCTK_INT :: i, j, k, l, nneigh, rep_loc1, rep_loc2 - CCTK_INT, dimension(1) :: lmax - CCTK_REAL, dimension(3) :: n0, n1, vec - CCTK_REAL, dimension(0:5) :: dp1, dp2, fnew - CCTK_REAL :: dfs1, dfs2 - -# include "include/physical_part.h" - - ncalls = ncalls + 1 - - rep_loc1 = 0 - - do k = kzl, kzr - do j = jyl, jyr - do i = ixl, ixr - if ( ( eh_mask(i,j,k) .ge. 0 ) .and. ( rep_mask(i,j,k) .eq. 0 ) ) then - nneigh = 0 - fnew = zero - dfs1 = one / sqrt ( dfx(i,j,k)**2 + dfy(i,j,k)**2 + dfz(i,j,k)**2 ) - n0(1) = dfx(i,j,k) * dfs1 - n0(2) = dfy(i,j,k) * dfs1 - n0(3) = dfz(i,j,k) * dfs1 - do l = 0, 5 - if ( .not. btest(eh_mask(i,j,k),l) ) then - if ( rep_mask(i+ix(l),j+jy(l),k+kz(l)) .eq. 1 ) then - vec(1) = -ix(l) * cctk_delta_space(1) - vec(2) = -jy(l) * cctk_delta_space(2) - vec(3) = -kz(l) * cctk_delta_space(3) - dfs2 = one / sqrt ( dfx(i+ix(l),j+jy(l),k+kz(l))**2 + & - dfy(i+ix(l),j+jy(l),k+kz(l))**2 + & - dfz(i+ix(l),j+jy(l),k+kz(l))**2 ) - n1(1) = dfx(i+ix(l),j+jy(l),k+kz(l)) * dfs2 - n1(2) = dfy(i+ix(l),j+jy(l),k+kz(l)) * dfs2 - n1(3) = dfz(i+ix(l),j+jy(l),k+kz(l)) * dfs2 - dp1(nneigh) = vec(1) * n0(1) + vec(2) * n0(2) + vec(3) * n0(3) - dp2(nneigh) = vec(1) * n1(1) + vec(2) * n1(2) + vec(3) * n1(3) - fnew(nneigh) = f(i+ix(l),j+jy(l),k+kz(l)) + & - half * ( dp1(nneigh) + dp2(nneigh) ) - nneigh = nneigh + 1 - - end if - end if - end do - - if ( nneigh .gt. 0 ) then - f(i,j,k) = sum(fnew(0:nneigh-1)) / nneigh - rep_mask(i,j,k) = 2 - rep_loc1 = rep_loc1 + 1 - end if - end if - end do - end do - end do - - call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, sum_handle, & - rep_loc1, rep_loc2, CCTK_VARIABLE_INT ) - if ( ierr .ne. 0 ) then - call CCTK_WARN(0,'Sum of rep_count failed') - end if -! print*,rep_loc1,rep_loc2 - - rep_count = rep_count + rep_loc2 -! print*,ncalls,rep_count,rep_total - - if ( rep_count .ge. rep_total ) then - re_param_control_approx = 0 - call CCTK_INFO('Approximation re-parametrization complete') - end if - -! call CCTK_WARN(0,'Testing') -! print*,'after' -! print*,rep_mask(24,21,1:4),eh_mask(24,21,1:4) -! print*,f(24,21,1:4) -! print* -end subroutine EHFinder_ReParametrize6 - - -subroutine EHFinder_ReParametrize7(CCTK_ARGUMENTS) - - use EHFinder_mod - - implicit none - - DECLARE_CCTK_PARAMETERS - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_FUNCTIONS - - where ( rep_mask .eq. 2 ) rep_mask = 1 - -end subroutine EHFinder_ReParametrize7 - - -CCTK_REAL function solve_quadratic ( a, b, c ) - - use EHFinder_Constants - - CCTK_REAL, intent(in) :: a, b, c - CCTK_REAL :: d2, q - - d2 = b**2 - four * a * c - - if ( d2 .ge. zero ) then - q = - half * ( b + sign(one,b) * sqrt(d2) ) - - solve_quadratic = min ( abs(q/a), abs(c/q) ) - else - solve_quadratic = huge - end if - -end function solve_quadratic +!subroutine EHFinder_ReParametrize5(CCTK_ARGUMENTS) +! +! use EHFinder_mod +! +! implicit none +! +! DECLARE_CCTK_PARAMETERS +! DECLARE_CCTK_ARGUMENTS +! DECLARE_CCTK_FUNCTIONS +! +! CCTK_INT :: i, j, k +! CCTK_INT :: rep_countloc, rep_totalloc +! CCTK_REAL :: idx2, idy2, idz2, idx, idy, idz +! CCTK_REAL :: d2fdx2, d2fdy2, d2fdz2, d2fdxdy, d2fdxdz, d2fdydz +! CCTK_REAL :: n_x, n_y, n_z +! CCTK_REAL :: a, b, c, lambda +! +! interface +! CCTK_REAL function solve_quadratic ( a, b, c ) +! CCTK_REAL, intent(in) :: a, b, c +! end function solve_quadratic +! end interface +! +! if ( re_param_control_approx .eq. 1 ) then +! if ( reparametrize_every_approx .gt. 0 ) then +! if ( mod(cctk_iteration,reparametrize_every_approx) .eq. 0 ) then +! call CCTK_ReductionArrayHandle ( sum_handle, 'sum' ) +! if ( sum_handle .lt. 0 ) then +! call CCTK_WARN(0,'Could not obtain a handle for sum reduction') +! end if +! +!# include "include/centered_second.h" +! +! ncalls = 0 +! +! rep_totalloc = count ( eh_mask(ixl:ixr,jyl:jyr,kzl:kzr) .ge. 0 ) +! +! rep_countloc = 0 +! +! ftmp = f +! +! ixr = min ( ixr, nx-1 ) +! jyr = min ( jyr, ny-1 ) +! kzr = min ( kzr, nz-1 ) +! +! do k = kzl, kzr +! do j = jyl, jyr +! do i = ixl, ixr +! if ( eh_mask(i,j,k) .eq. 0 ) then +! if ( ( f(i,j,k)*f(i-1,j,k) .le. zero ) .or. & +! ( f(i,j,k)*f(i+1,j,k) .le. zero ) .or. & +! ( f(i,j,k)*f(i,j-1,k) .le. zero ) .or. & +! ( f(i,j,k)*f(i,j+1,k) .le. zero ) .or. & +! ( f(i,j,k)*f(i,j,k-1) .le. zero ) .or. & +! ( f(i,j,k)*f(i,j,k+1) .le. zero ) ) then +! +! rep_mask(i,j,k) = 1 +! rep_countloc = rep_countloc + 1 +! +! d2fdx2 = idx2 * ( f(i-1,j,k) - two * f(i,j,k) + f(i+1,j,k) ) +! d2fdy2 = idy2 * ( f(i,j-1,k) - two * f(i,j,k) + f(i,j+1,k) ) +! d2fdz2 = idz2 * ( f(i,j,k-1) - two * f(i,j,k) + f(i,j,k+1) ) +! +! d2fdxdy = idx * idy * ( f(i-1,j-1,k) - f(i+1,j-1,k) - & +! f(i-1,j+1,k) + f(i+1,j+1,k) ) +! d2fdxdz = idx * idz * ( f(i-1,j,k-1) - f(i+1,j,k-1) - & +! f(i-1,j,k+1) + f(i+1,j,k+1) ) +! d2fdydz = idy * idz * ( f(i,j-1,k-1) - f(i,j+1,k-1) - & +! f(i,j-1,k+1) + f(i,j+1,k+1) ) +! +! c = f(i,j,k) +! b = sqrt ( dfx(i,j,k)**2 + dfy(i,j,k)**2 + dfz(i,j,k)**2 ) +! n_x = dfx(i,j,k) / b +! n_y = dfy(i,j,k) / b +! n_z = dfz(i,j,k) / b +! a = half*( n_x**2*d2fdx2 + n_y**2*d2fdy2 + n_z**2*d2fdz2 ) +& +! n_x*n_y*d2fdxdy + n_x*n_z*d2fdxdz + n_y*n_z*d2fdydz +! +! if ( a .ne. 0 ) then +! lambda = solve_quadratic ( a, b, c ) +! else +! lambda = abs ( c / b ) +! endif +! if ( lambda .eq. huge ) then +! lambda = abs ( c / b ) +! end if +! +! if ( lambda .gt. two*delta ) then +! print*,i,j,k +! print*,f(i-1:i+1,j-1:j+1,k-1:k+1) +! print*,eh_mask(i-1:i+1,j-1:j+1,k-1:k+1) +! print*,a,b,c,lambda +! call CCTK_WARN(0,'Re-parametrization failed') +! end if +! +! ftmp(i,j,k) = sign ( lambda, f(i,j,k) ) +!! print*,i,j,k +!! print*,lambda,f(i,j,k), ftmp(i,j,k) +!! print*,f(i-1,j,k),f(i+1,j,k) +!! print*,f(i,j-1,k),f(i,j+1,k) +!! print*,f(i,j,k-1),f(i,j,k+1) +! end if +! end if +! end do +! end do +! end do +! +! where ( eh_mask .ge. 0 ) f = ftmp +! +!! print*,rep_countloc,rep_totalloc +! call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, sum_handle, & +! rep_countloc, rep_count, CCTK_VARIABLE_INT ) +! if ( ierr .ne. 0 ) then +! call CCTK_WARN(0,'Sum of rep_count failed') +! end if +! +! call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, sum_handle, & +! rep_totalloc, rep_total, CCTK_VARIABLE_INT ) +! if ( ierr .ne. 0 ) then +! call CCTK_WARN(0,'Sum of rep_total failed') +! end if +!! print*,'First step' +!! print*,rep_count,rep_total +!! call CCTK_WARN(0,'Finished') +! +!! call CartSymGN(ierr,cctkGH,'ehfinder::level_set') +!! call CartSymGN(ierr,cctkGH,'ehfinder::rep_mask') +! +! end if +! end if +! +! end if +! +!end subroutine EHFinder_ReParametrize5 + + +!subroutine EHFinder_ReParametrize6(CCTK_ARGUMENTS) +! +! use EHFinder_mod +! +! implicit none +! +! DECLARE_CCTK_PARAMETERS +! DECLARE_CCTK_ARGUMENTS +! DECLARE_CCTK_FUNCTIONS +! +! CCTK_INT :: i, j, k, l, nneigh, rep_loc1, rep_loc2 +! CCTK_INT, dimension(1) :: lmax +! CCTK_REAL, dimension(3) :: n0, n1, vec +! CCTK_REAL, dimension(0:5) :: dp1, dp2, fnew +! CCTK_REAL :: dfs1, dfs2 +! +!# include "include/physical_part.h" +! +! ncalls = ncalls + 1 +! +! rep_loc1 = 0 +! +! do k = kzl, kzr +! do j = jyl, jyr +! do i = ixl, ixr +! if ( ( eh_mask(i,j,k) .ge. 0 ) .and. ( rep_mask(i,j,k) .eq. 0 ) ) then +! nneigh = 0 +! fnew = zero +! dfs1 = one / sqrt ( dfx(i,j,k)**2 + dfy(i,j,k)**2 + dfz(i,j,k)**2 ) +! n0(1) = dfx(i,j,k) * dfs1 +! n0(2) = dfy(i,j,k) * dfs1 +! n0(3) = dfz(i,j,k) * dfs1 +! do l = 0, 5 +! if ( .not. btest(eh_mask(i,j,k),l) ) then +! if ( rep_mask(i+ix(l),j+jy(l),k+kz(l)) .eq. 1 ) then +! vec(1) = -ix(l) * cctk_delta_space(1) +! vec(2) = -jy(l) * cctk_delta_space(2) +! vec(3) = -kz(l) * cctk_delta_space(3) +! dfs2 = one / sqrt ( dfx(i+ix(l),j+jy(l),k+kz(l))**2 + & +! dfy(i+ix(l),j+jy(l),k+kz(l))**2 + & +! dfz(i+ix(l),j+jy(l),k+kz(l))**2 ) +! n1(1) = dfx(i+ix(l),j+jy(l),k+kz(l)) * dfs2 +! n1(2) = dfy(i+ix(l),j+jy(l),k+kz(l)) * dfs2 +! n1(3) = dfz(i+ix(l),j+jy(l),k+kz(l)) * dfs2 +! dp1(nneigh) = vec(1) * n0(1) + vec(2) * n0(2) + vec(3) * n0(3) +! dp2(nneigh) = vec(1) * n1(1) + vec(2) * n1(2) + vec(3) * n1(3) +! fnew(nneigh) = f(i+ix(l),j+jy(l),k+kz(l)) + & +! half * ( dp1(nneigh) + dp2(nneigh) ) +! nneigh = nneigh + 1 +! +! end if +! end if +! end do +! +! if ( nneigh .gt. 0 ) then +! f(i,j,k) = sum(fnew(0:nneigh-1)) / nneigh +! rep_mask(i,j,k) = 2 +! rep_loc1 = rep_loc1 + 1 +! end if +! end if +! end do +! end do +! end do +! +! call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, sum_handle, & +! rep_loc1, rep_loc2, CCTK_VARIABLE_INT ) +! if ( ierr .ne. 0 ) then +! call CCTK_WARN(0,'Sum of rep_count failed') +! end if +!! print*,rep_loc1,rep_loc2 +! +! rep_count = rep_count + rep_loc2 +!! print*,ncalls,rep_count,rep_total +! +! if ( rep_count .ge. rep_total ) then +! re_param_control_approx = 0 +! call CCTK_INFO('Approximation re-parametrization complete') +! end if +! +!! call CCTK_WARN(0,'Testing') +!! print*,'after' +!! print*,rep_mask(24,21,1:4),eh_mask(24,21,1:4) +!! print*,f(24,21,1:4) +!! print* +!end subroutine EHFinder_ReParametrize6 + + +!subroutine EHFinder_ReParametrize7(CCTK_ARGUMENTS) +! +! use EHFinder_mod +! +! implicit none +! +! DECLARE_CCTK_PARAMETERS +! DECLARE_CCTK_ARGUMENTS +! DECLARE_CCTK_FUNCTIONS +! +! where ( rep_mask .eq. 2 ) rep_mask = 1 +! +!end subroutine EHFinder_ReParametrize7 + + +!CCTK_REAL function solve_quadratic ( a, b, c ) +! +! use EHFinder_Constants +! +! CCTK_REAL, intent(in) :: a, b, c +! CCTK_REAL :: d2, q +! +! d2 = b**2 - four * a * c +! +! if ( d2 .ge. zero ) then +! q = - half * ( b + sign(one,b) * sqrt(d2) ) +! +! solve_quadratic = min ( abs(q/a), abs(c/q) ) +! else +! solve_quadratic = huge +! end if +! +!end function solve_quadratic diff --git a/src/EHFinder_SetMask.F90 b/src/EHFinder_SetMask.F90 index a55c410..40f057d 100644 --- a/src/EHFinder_SetMask.F90 +++ b/src/EHFinder_SetMask.F90 @@ -42,22 +42,22 @@ subroutine EHFinder_MaskInit(CCTK_ARGUMENTS) ! eh_mask appropriately. The array ll is defined in EHFinder_mod.F90 ! and contains the values ( 1, 2, 4, 8, 16, 32 ). if ( ( cctk_bbox(1) .eq. 1 ) .and. ( ixl .eq. 1 ) ) then - eh_mask(1,:,:) = eh_mask(1,:,:) + ll(0) + eh_mask(1,:,:,:) = eh_mask(1,:,:,:) + ll(0) end if if ( ( cctk_bbox(2) .eq. 1 ) .and. ( ixr .eq. nx ) ) then - eh_mask(nx,:,:) = eh_mask(nx,:,:) + ll(1) + eh_mask(nx,:,:,:) = eh_mask(nx,:,:,:) + ll(1) end if if ( ( cctk_bbox(3) .eq. 1 ) .and. ( jyl .eq. 1 ) ) then - eh_mask(:,1,:) = eh_mask(:,1,:) + ll(2) + eh_mask(:,1,:,:) = eh_mask(:,1,:,:) + ll(2) end if if ( ( cctk_bbox(4) .eq. 1 ) .and. ( jyr .eq. ny ) ) then - eh_mask(:,ny,:) = eh_mask(:,ny,:) + ll(3) + eh_mask(:,ny,:,:) = eh_mask(:,ny,:,:) + ll(3) end if if ( ( cctk_bbox(5) .eq. 1 ) .and. ( kzl .eq. 1 ) ) then - eh_mask(:,:,1) = eh_mask(:,:,1) + ll(4) + eh_mask(:,:,1,:) = eh_mask(:,:,1,:) + ll(4) end if if ( ( cctk_bbox(6) .eq. 1 ) .and. ( kzr .eq. nz ) ) then - eh_mask(:,:,nz) = eh_mask(:,:,nz) + ll(5) + eh_mask(:,:,nz,:) = eh_mask(:,:,nz,:) + ll(5) end if return @@ -78,7 +78,7 @@ subroutine EHFinder_SetMask1(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS - CCTK_INT :: i, j, k, i1, i2, j1, j2, k1, k2 + CCTK_INT :: i, j, k, l, i1, i2, j1, j2, k1, k2 logical :: active CCTK_INT, dimension(3) :: imin_loc, imax_loc, imin_n, imax_n @@ -109,641 +109,638 @@ subroutine EHFinder_SetMask1(CCTK_ARGUMENTS) end if ! If the reparametrization was not undone... - if ( active .and. .not.reparam_undone ) then + if ( active .and. .not. all(reparam_undone) ) then -! Get the minimum and maximum index excluding ghost and symmetry cells. +! Get the minimum and maximum index excluding ghost and symmetry cells. # include "include/physical_part.h" -! Store the current mask and level set function - tm_mask(ixl:ixr,jyl:jyr,kzl:kzr) = eh_mask(ixl:ixr,jyl:jyr,kzl:kzr) - ftmp(ixl:ixr,jyl:jyr,kzl:kzr) = f(ixl:ixr,jyl:jyr,kzl:kzr) - -! Next we try to locate the minimum and maximum of the global indeces of -! the currently active cells giving us the smallest rectangular box -! containing all active cells. - - if ( use_outer_excision .gt. 0 ) then -! First initialize some variables. -! fimin_loc, and fimax_loc are 3 element arrays that will contain the -! minimum value of f on the boundaries of this rectangular box. They -! will be used to decide if the box should be changed. They are -! initialized to the negative of the value used in excised cells. - -! The 3 element arrays imin_loc and imax_loc will contain the min and -! max global indices of the ractangular box. They are initialised to -! the maximum global index and 0, respectively. - - fimin_loc = -ex_value; fimax_loc = -ex_value - imin_loc = cctk_gsh; imax_loc = 0 - -! Find all cells with f-shell_width*delta then activate the point by setting -! the temporary mask to zero. The new value of f will be the value -! if f in its neighbour point - delta. Do this for all directions. - - if ( ( eh_mask(i,j,k) .eq. -1 ) .and. ( f(i,j,k) .lt. 0 ) ) then - if ( eh_mask(i-1,j,k) .ge. 0 ) then - if ( f(i-1,j,k) - delta .gt. -shell_width * delta ) then - tm_mask(i,j,k) = 0 - ftmp(i,j,k) = f(i-1,j,k) - delta - end if - end if +! Now check and see if any interior excised points need to be activated. + if ( use_inner_excision .gt. 0 ) then + do k = kzl, kzr + do j = jyl, jyr + do i = ixl, ixr + +! If an interior excised point has a non-excised neighbour where +! f-delta>-shell_width*delta then activate the point by setting +! the temporary mask to zero. The new value of f will be the value +! if f in its neighbour point - delta. Do this for all directions. + + if ( ( eh_mask(i,j,k,l) .eq. -1 ) .and. & + ( f(i,j,k,l) .lt. 0 ) ) then + if ( eh_mask(i-1,j,k,l) .ge. 0 ) then + if ( f(i-1,j,k,l) - delta .gt. -shell_width * delta ) then + tm_mask(i,j,k,l) = 0 + ftmp(i,j,k,l) = f(i-1,j,k,l) - delta + end if + end if + + if ( eh_mask(i+1,j,k,l) .ge. 0 ) then + if ( f(i+1,j,k,l) - delta .gt. -shell_width * delta ) then + tm_mask(i,j,k,l) = 0 + ftmp(i,j,k,l) = f(i+1,j,k,l) - delta + end if + end if + + if ( eh_mask(i,j-1,k,l) .ge. 0 ) then + if ( f(i,j-1,k,l) - delta .gt. -shell_width * delta ) then + tm_mask(i,j,k,l) = 0 + ftmp(i,j,k,l) = f(i,j-1,k,l) - delta + end if + end if + + if ( eh_mask(i,j+1,k,l) .ge. 0 ) then + if ( f(i,j+1,k,l) - delta .gt. -shell_width * delta ) then + tm_mask(i,j,k,l) = 0 + ftmp(i,j,k,l) = f(i,j+1,k,l) - delta + end if + end if + + if ( eh_mask(i,j,k-1,l) .ge. 0 ) then + if ( f(i,j,k-1,l) - delta .gt. -shell_width * delta ) then + tm_mask(i,j,k,l) = 0 + ftmp(i,j,k,l) = f(i,j,k-1,l) - delta + end if + end if + + if ( eh_mask(i,j,k+1,l) .ge. 0 ) then + if ( f(i,j,k+1,l) - delta .gt. -shell_width * delta ) then + tm_mask(i,j,k,l) = 0 + ftmp(i,j,k,l) = f(i,j,k+1,l) - delta + end if + end if - if ( eh_mask(i+1,j,k) .ge. 0 ) then - if ( f(i+1,j,k) - delta .gt. -shell_width * delta ) then - tm_mask(i,j,k) = 0 - ftmp(i,j,k) = f(i+1,j,k) - delta end if - end if + end do + end do + end do + end if - if ( eh_mask(i,j-1,k) .ge. 0 ) then - if ( f(i,j-1,k) - delta .gt. -shell_width * delta ) then - tm_mask(i,j,k) = 0 - ftmp(i,j,k) = f(i,j-1,k) - delta - end if - end if +! Check and see if the boundary of the exterior excision region should +! be changed and if so find the indices describing the new excision +! region. + if ( use_outer_excision .gt. 0 ) then + do i = 1, 3 + imin_n(i) = imin_glob(i) + imax_n(i) = imax_glob(i) +! If the minimum value of f on a face on the box plus delta is less +! than shell_width * delta then the active region has to be increased +! in size in the corresponding region, i.e. inactive cells has to be +! activated. + if ( fimin_glob(i) + delta .lt. shell_width * delta ) then + if ( imin_glob(i) .gt. 1 ) then + imin_n(i) = imin_glob(i) - 1 + endif + end if + if ( fimax_glob(i) + delta .lt. shell_width * delta ) then + if ( imax_glob(i) .lt. cctk_gsh(i) ) then + imax_n(i) = imax_glob(i) + 1 + endif + end if + end do - if ( eh_mask(i,j+1,k) .ge. 0 ) then - if ( f(i,j+1,k) - delta .gt. -shell_width * delta ) then - tm_mask(i,j,k) = 0 - ftmp(i,j,k) = f(i,j+1,k) - delta - end if +! 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. + if ( imin_n(1) .ne. imin_glob(1) ) then + i1 = imin_n(1) - cctk_lbnd(1) + if ( ( max(ixl,2) .le. i1 ) .and. ( i1 .le. ixr ) ) then + 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 ) + 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 - - if ( eh_mask(i,j,k-1) .ge. 0 ) then - if ( f(i,j,k-1) - delta .gt. -shell_width * delta ) then - tm_mask(i,j,k) = 0 - ftmp(i,j,k) = f(i,j,k-1) - delta - end if + end if + end if + if ( imax_n(1) .ne. imax_glob(1) ) then + i2 = imax_n(1) - cctk_lbnd(1) + if ( ( ixl .lt. i2 ) .and. ( i2 .le. min(ixr,nx-1) ) ) then + 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 ) + 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 - - if ( eh_mask(i,j,k+1) .ge. 0 ) then - if ( f(i,j,k+1) - delta .gt. -shell_width * delta ) then - tm_mask(i,j,k) = 0 - ftmp(i,j,k) = f(i,j,k+1) - delta - end if + end if + end if + if ( imin_n(2) .ne. imin_glob(2) ) then + j1 = imin_n(2) - cctk_lbnd(2) + if ( ( max(jyl,2) .le. j1 ) .and. ( j1 .lt. jyr ) ) then + i1 = max ( max ( imin_n(1), imin_glob(1) ) - cctk_lbnd(1), ixl ) + 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 do - end do - end do - end if + end if + if ( imax_n(2) .ne. imax_glob(2) ) then + j2 = imax_n(2) - cctk_lbnd(2) + if ( ( jyl .le. j2 ) .and. ( j2 .le. min(jyr,ny-1) ) ) then + i1 = max ( max ( imin_n(1), imin_glob(1) ) - cctk_lbnd(1), ixl ) + 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 + if ( imin_n(3) .ne. imin_glob(3) ) then + k1 = imin_n(3) - cctk_lbnd(3) + if ( ( max(kzl,2) .le. k1 ) .and. ( k1 .le. kzr ) ) then + i1 = max ( max ( imin_n(1), imin_glob(1) ) - cctk_lbnd(1), ixl ) + 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 + if ( imax_n(3) .ne. imax_glob(3) ) then + k2 = imax_n(3) - cctk_lbnd(3) + if ( ( kzl .lt. k2 ) .and. ( k2 .le. min(kzr,nz-1) ) ) then + i1 = max ( max ( imin_n(1), imin_glob(1) ) - cctk_lbnd(1), ixl ) + 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 -! Check and see if the boundary of the exterior excision region should -! be changed and if so find the indices describing the new excision region. - if ( use_outer_excision .gt. 0 ) then - do i = 1, 3 - imin_n(i) = imin_glob(i) - imax_n(i) = imax_glob(i) -! If the minimum value of f on a face on the box plus delta is less -! than shell_width * delta then the active region has to be increased -! in size in the corresponding region, i.e. inactive cells has to be -! activated. - if ( fimin_glob(i) + delta .lt. shell_width * delta ) then - if ( imin_glob(i) .gt. 1 ) then - imin_n(i) = imin_glob(i) - 1 - endif - end if - if ( fimax_glob(i) + delta .lt. shell_width * delta ) then - if ( imax_glob(i) .lt. cctk_gsh(i) ) then - imax_n(i) = imax_glob(i) + 1 - endif - 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. - if ( imin_n(1) .ne. imin_glob(1) ) then - i1 = imin_n(1) - cctk_lbnd(1) - if ( ( max(ixl,2) .le. i1 ) .and. ( i1 .le. ixr ) ) then - 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 ) - 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) = 0 - ftmp(i1,j1:j2,k1:k2) = f(i1+1,j1:j2,k1:k2) + delta -! print*,'done' +! Then do the edges of the box. + if ( ( imin_n(1) .ne. imin_glob(1) ) .and. & + ( imin_n(2) .ne. imin_glob(2) ) ) then + i1 = imin_n(1) - cctk_lbnd(1) + j1 = imin_n(2) - cctk_lbnd(2) + if ( ( max(ixl,2) .le. i1 ) .and. ( i1 .le. ixr ) .and. & + ( max(jyl,2) .le. j1 ) .and. ( j1 .le. jyr ) ) then + 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 ) + if ( k1 .le. k2 ) then + tm_mask(i1,j1,k1:k2,l) = 0 + ftmp(i1,j1,k1:k2,l) = f(i1+1,j1+1,k1:k2,l) + sqrt(two)*delta + end if + end if end if - end if - end if - if ( imax_n(1) .ne. imax_glob(1) ) then - i2 = imax_n(1) - cctk_lbnd(1) - if ( ( ixl .le. i2 ) .and. ( i2 .le. min(ixr,nx-1) ) ) then - 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 ) - 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) = 0 - ftmp(i2,j1:j2,k1:k2) = f(i2-1,j1:j2,k1:k2) + delta -! print*,'done' + if ( ( imin_n(1) .ne. imin_glob(1) ) .and. & + ( imax_n(2) .ne. imax_glob(2) ) ) then + i1 = imin_n(1) - cctk_lbnd(1) + j2 = imax_n(2) - cctk_lbnd(2) + if ( ( max(ixl,2) .le. i1 ) .and. ( i1 .le. ixr ) .and. & + ( jyl .le. j2 ) .and. ( j2 .le. min(jyr,ny-1) ) ) then + 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 ) + if ( k1 .le. k2 ) then + tm_mask(i1,j2,k1:k2,l) = 0 + ftmp(i1,j2,k1:k2,l) = f(i1+1,j2-1,k1:k2,l) + sqrt(two)*delta + end if + end if end if - end if - end if - if ( imin_n(2) .ne. imin_glob(2) ) then - j1 = imin_n(2) - cctk_lbnd(2) - if ( ( max(jyl,2) .le. j1 ) .and. ( j1 .le. jyr ) ) then - i1 = max ( max ( imin_n(1), imin_glob(1) ) - cctk_lbnd(1), ixl ) - 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) = 0 - ftmp(i1:i2,j1,k1:k2) = f(i1:i2,j1+1,k1:k2) + delta -! print*,'done' + if ( ( imin_n(1) .ne. imin_glob(1) ) .and. & + ( imin_n(3) .ne. imin_glob(3) ) ) then + i1 = imin_n(1) - cctk_lbnd(1) + k1 = imin_n(3) - cctk_lbnd(3) + if ( ( max(ixl,2) .le. i1 ) .and. ( i1 .le. ixr ) .and. & + ( max(kzl,2) .le. k1 ) .and. ( k1 .le. kzr ) ) then + 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 ) + if ( j1 .le. j2 ) then + tm_mask(i1,j1:j2,k1,l) = 0 + ftmp(i1,j1:j2,k1,l) = f(i1+1,j1:j2,k1+1,l) + sqrt(two)*delta + end if + end if end if - end if - end if - if ( imax_n(2) .ne. imax_glob(2) ) then - j2 = imax_n(2) - cctk_lbnd(2) - if ( ( jyl .le. j2 ) .and. ( j2 .le. min(jyr,ny-1) ) ) then - i1 = max ( max ( imin_n(1), imin_glob(1) ) - cctk_lbnd(1), ixl ) - 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) = 0 - ftmp(i1:i2,j2,k1:k2) = f(i1:i2,j2-1,k1:k2) + delta -! print*,'done' + if ( ( imin_n(1) .ne. imin_glob(1) ) .and. & + ( imax_n(3) .ne. imax_glob(3) ) ) then + i1 = imin_n(1) - cctk_lbnd(1) + k2 = imax_n(3) - cctk_lbnd(3) + if ( ( max(ixl,2) .le. i1 ) .and. ( i1 .le. ixr ) .and. & + ( kzl .le. k2 ) .and. ( k2 .le. min(kzr,nz-1) ) ) then + 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 ) + if ( j1 .le. j2 ) then + tm_mask(i1,j1:j2,k2,l) = 0 + ftmp(i1,j1:j2,k2,l) = f(i1+1,j1:j2,k2-1,l) + sqrt(two)*delta + end if + end if end if - end if - end if - if ( imin_n(3) .ne. imin_glob(3) ) then - k1 = imin_n(3) - cctk_lbnd(3) - print*,'Debug 1 : ', imin_n(3), imin_glob(3), k1, kzl, nz - if ( ( max(kzl,2) .le. k1 ) .and. ( k1 .le. kzr ) ) then - i1 = max ( max ( imin_n(1), imin_glob(1) ) - cctk_lbnd(1), ixl ) - 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) = 0 - ftmp(i1:i2,j1:j2,k1) = f(i1:i2,j1:j2,k1+1) + delta -! print*,'done' + if ( ( imax_n(1) .ne. imax_glob(1) ) .and. & + ( imin_n(2) .ne. imin_glob(2) ) ) then + i2 = imax_n(1) - cctk_lbnd(1) + j1 = imin_n(2) - cctk_lbnd(2) + if ( ( ixl .le. i2 ) .and. ( i2 .le. min(ixr,nx-1) ) .and. & + ( max(jyl,2) .le. j1 ) .and. ( j1 .le. jyr ) ) then + 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 ) + if ( k1 .le. k2 ) then + tm_mask(i2,j1,k1:k2,l) = 0 + ftmp(i2,j1,k1:k2,l) = f(i2-1,j1+1,k1:k2,l) + sqrt(two)*delta + end if + end if end if - end if - end if - if ( imax_n(3) .ne. imax_glob(3) ) then - k2 = imax_n(3) - cctk_lbnd(3) - print*, 'Debug 2 : ', imax_n(3), imax_glob(3), k2, kzl, nz-1 - if ( ( kzl .le. k2 ) .and. ( k2 .le. min(kzr,nz-1) ) ) then - i1 = max ( max ( imin_n(1), imin_glob(1) ) - cctk_lbnd(1), ixl ) - 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) = 0 - ftmp(i1:i2,j1:j2,k2) = f(i1:i2,j1:j2,k2-1) + delta -! print*,'done' + if ( ( imax_n(1) .ne. imax_glob(1) ) .and. & + ( imax_n(2) .ne. imax_glob(2) ) ) then + i2 = imax_n(1) - cctk_lbnd(1) + j2 = imax_n(2) - cctk_lbnd(2) + if ( ( ixl .le. i2 ) .and. ( i2 .le. min(ixr,nx-1) ) .and. & + ( jyl .le. j2 ) .and. ( j2 .le. min(jyr,ny-1) ) ) then + 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 ) + if ( k1 .le. k2 ) then + tm_mask(i2,j2,k1:k2,l) = 0 + ftmp(i2,j2,k1:k2,l) = f(i2-1,j2-1,k1:k2,l) + sqrt(two)*delta + end if + end if end if - end if - end if - -! Then do the edges of the box. - if ( ( imin_n(1) .ne. imin_glob(1) ) .and. & - ( imin_n(2) .ne. imin_glob(2) ) ) then - i1 = imin_n(1) - cctk_lbnd(1) - j1 = imin_n(2) - cctk_lbnd(2) - if ( ( max(ixl,2) .le. i1 ) .and. ( i1 .le. ixr ) .and. & - ( max(jyl,2) .le. j1 ) .and. ( j1 .le. jyr ) ) then - 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 ) - if ( k1 .le. k2 ) then - tm_mask(i1,j1,k1:k2) = 0 - ftmp(i1,j1,k1:k2) = f(i1+1,j1+1,k1:k2) + sqrt(two)*delta + if ( ( imax_n(1) .ne. imax_glob(1) ) .and. & + ( imin_n(3) .ne. imin_glob(3) ) ) then + i2 = imax_n(1) - cctk_lbnd(1) + k1 = imin_n(3) - cctk_lbnd(3) + if ( ( ixl .le. i2 ) .and. ( i2 .le. min(ixr,nx-1) ) .and. & + ( max(kzl,2) .le. k1 ) .and. ( k1 .le. kzr ) ) then + 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 ) + if ( j1 .le. j2 ) then + tm_mask(i2,j1:j2,k1,l) = 0 + ftmp(i2,j1:j2,k1,l) = f(i2-1,j1:j2,k1+1,l) + sqrt(two)*delta + end if + end if end if - end if - end if - if ( ( imin_n(1) .ne. imin_glob(1) ) .and. & - ( imax_n(2) .ne. imax_glob(2) ) ) then - i1 = imin_n(1) - cctk_lbnd(1) - j2 = imax_n(2) - cctk_lbnd(2) - if ( ( max(ixl,2) .le. i1 ) .and. ( i1 .le. ixr ) .and. & - ( jyl .le. j2 ) .and. ( j2 .le. min(jyr,ny-1) ) ) then - 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 ) - if ( k1 .le. k2 ) then - tm_mask(i1,j2,k1:k2) = 0 - ftmp(i1,j2,k1:k2) = f(i1+1,j2-1,k1:k2) + sqrt(two)*delta + if ( ( imax_n(1) .ne. imax_glob(1) ) .and. & + ( imax_n(3) .ne. imax_glob(3) ) ) then + i2 = imax_n(1) - cctk_lbnd(1) + k2 = imax_n(3) - cctk_lbnd(3) + if ( ( ixl .le. i2 ) .and. ( i2 .le. min(ixr,nx-1) ) .and. & + ( kzl .le. k2 ) .and. ( k2 .le. min(kzr,nz-1) ) ) then + 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 ) + if ( j1 .le. j2 ) then + tm_mask(i2,j1:j2,k2,l) = 0 + ftmp(i2,j1:j2,k2,l) = f(i2-1,j1:j2,k2-1,l) + sqrt(two)*delta + end if + end if end if - end if - end if - if ( ( imin_n(1) .ne. imin_glob(1) ) .and. & - ( imin_n(3) .ne. imin_glob(3) ) ) then - i1 = imin_n(1) - cctk_lbnd(1) - k1 = imin_n(3) - cctk_lbnd(3) - if ( ( max(ixl,2) .le. i1 ) .and. ( i1 .le. ixr ) .and. & - ( max(kzl,2) .le. k1 ) .and. ( k1 .le. kzr ) ) then - 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 ) - if ( j1 .le. j2 ) then - tm_mask(i1,j1:j2,k1) = 0 - ftmp(i1,j1:j2,k1) = f(i1+1,j1:j2,k1+1) + sqrt(two)*delta + if ( ( imin_n(2) .ne. imin_glob(2) ) .and. & + ( imin_n(3) .ne. imin_glob(3) ) ) then + j1 = imin_n(2) - cctk_lbnd(2) + k1 = imin_n(3) - cctk_lbnd(3) + if ( ( max(jyl,2) .le. j1 ) .and. ( j1 .le. jyr ) .and. & + ( max(kzl,2) .le. k1 ) .and. ( k1 .le. kzr ) ) then + i1 = max ( max ( imin_n(1), imin_glob(1) ) - cctk_lbnd(1), ixl ) + i2 = min ( min ( imax_n(1), imax_glob(1) ) - cctk_lbnd(1), ixr ) + if ( i1 .le. i2 ) then + tm_mask(i1:i2,j1,k1,l) = 0 + ftmp(i1:i2,j1,k1,l) = f(i1:i2,j1+1,k1+1,l) + sqrt(two)*delta + end if + end if end if - end if - end if - if ( ( imin_n(1) .ne. imin_glob(1) ) .and. & - ( imax_n(3) .ne. imax_glob(3) ) ) then - i1 = imin_n(1) - cctk_lbnd(1) - k2 = imax_n(3) - cctk_lbnd(3) - if ( ( max(ixl,2) .le. i1 ) .and. ( i1 .le. ixr ) .and. & - ( kzl .le. k2 ) .and. ( k2 .le. min(kzr,nz-1) ) ) then - 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 ) - if ( j1 .le. j2 ) then - tm_mask(i1,j1:j2,k2) = 0 - ftmp(i1,j1:j2,k2) = f(i1+1,j1:j2,k2-1) + sqrt(two)*delta + if ( ( imin_n(2) .ne. imin_glob(2) ) .and. & + ( imax_n(3) .ne. imax_glob(3) ) ) then + j1 = imin_n(2) - cctk_lbnd(2) + k2 = imax_n(3) - cctk_lbnd(3) + if ( ( max(jyl,2) .le. j1 ) .and. ( j1 .le. jyr ) .and. & + ( kzl .le. k2 ) .and. ( k2 .le. min(kzr,nz-1) ) ) then + i1 = max ( max ( imin_n(1), imin_glob(1) ) - cctk_lbnd(1), ixl ) + i2 = min ( min ( imax_n(1), imax_glob(1) ) - cctk_lbnd(1), ixr ) + if ( i1 .le. i2 ) then + tm_mask(i1:i2,j1,k2,l) = 0 + ftmp(i1:i2,j1,k2,l) = f(i1:i2,j1+1,k2-1,l) + sqrt(two)*delta + end if + end if end if - end if - end if - if ( ( imax_n(1) .ne. imax_glob(1) ) .and. & - ( imin_n(2) .ne. imin_glob(2) ) ) then - i2 = imax_n(1) - cctk_lbnd(1) - j1 = imin_n(2) - cctk_lbnd(2) - if ( ( ixl .le. i2 ) .and. ( i2 .le. min(ixr,nx-1) ) .and. & - ( max(jyl,2) .le. j1 ) .and. ( j1 .le. jyr ) ) then - 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 ) - if ( k1 .le. k2 ) then - tm_mask(i2,j1,k1:k2) = 0 - ftmp(i2,j1,k1:k2) = f(i2-1,j1+1,k1:k2) + sqrt(two)*delta + if ( ( imax_n(2) .ne. imax_glob(2) ) .and. & + ( imin_n(3) .ne. imin_glob(3) ) ) then + j2 = imax_n(2) - cctk_lbnd(2) + k1 = imin_n(3) - cctk_lbnd(3) + if ( ( jyl .le. j2 ) .and. ( j2 .le. min(jyr,ny-1) ) .and. & + ( max(kzl,2) .le. k1 ) .and. ( k1 .le. kzr ) ) then + i1 = max ( max ( imin_n(1), imin_glob(1) ) - cctk_lbnd(1), ixl ) + i2 = min ( min ( imax_n(1), imax_glob(1) ) - cctk_lbnd(1), ixr ) + if ( i1 .le. i2 ) then + tm_mask(i1:i2,j2,k1,l) = 0 + ftmp(i1:i2,j2,k1,l) = f(i1:i2,j2-1,k1+1,l) + sqrt(two)*delta + end if + end if end if - end if - end if - if ( ( imax_n(1) .ne. imax_glob(1) ) .and. & - ( imax_n(2) .ne. imax_glob(2) ) ) then - i2 = imax_n(1) - cctk_lbnd(1) - j2 = imax_n(2) - cctk_lbnd(2) - if ( ( ixl .le. i2 ) .and. ( i2 .le. min(ixr,nx-1) ) .and. & - ( jyl .le. j2 ) .and. ( j2 .le. min(jyr,ny-1) ) ) then - 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 ) - if ( k1 .le. k2 ) then - tm_mask(i2,j2,k1:k2) = 0 - ftmp(i2,j2,k1:k2) = f(i2-1,j2-1,k1:k2) + sqrt(two)*delta + if ( ( imax_n(2) .ne. imax_glob(2) ) .and. & + ( imax_n(3) .ne. imax_glob(3) ) ) then + j2 = imax_n(2) - cctk_lbnd(2) + k2 = imax_n(3) - cctk_lbnd(3) + if ( ( jyl .le. j2 ) .and. ( j2 .le. min(jyr,ny-1) ) .and. & + ( kzl .le. k2 ) .and. ( k2 .le. min(kzr,nz-1) ) ) then + i1 = max ( max ( imin_n(1), imin_glob(1) ) - cctk_lbnd(1), ixl ) + i2 = min ( min ( imax_n(1), imax_glob(1) ) - cctk_lbnd(1), ixr ) + if ( i1 .le. i2 ) then + tm_mask(i1:i2,j2,k2,l) = 0 + ftmp(i1:i2,j2,k2,l) = f(i1:i2,j2-1,k2-1,l) + sqrt(two)*delta + end if + end if end if - end if - end if - if ( ( imax_n(1) .ne. imax_glob(1) ) .and. & - ( imin_n(3) .ne. imin_glob(3) ) ) then - i2 = imax_n(1) - cctk_lbnd(1) - k1 = imin_n(3) - cctk_lbnd(3) - if ( ( ixl .le. i2 ) .and. ( i2 .le. min(ixr,nx-1) ) .and. & - ( max(kzl,2) .le. k1 ) .and. ( k1 .le. kzr ) ) then - 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 ) - if ( j1 .le. j2 ) then - tm_mask(i2,j1:j2,k1) = 0 - ftmp(i2,j1:j2,k1) = f(i2-1,j1:j2,k1+1) + sqrt(two)*delta + +! And finally do the corners. + if ( ( imin_n(1) .ne. imin_glob(1) ) .and. & + ( imin_n(2) .ne. imin_glob(2) ) .and. & + ( imin_n(3) .ne. imin_glob(3) ) ) then + i1 = imin_n(1) - cctk_lbnd(1) + j1 = imin_n(2) - cctk_lbnd(2) + k1 = imin_n(3) - cctk_lbnd(3) + if ( ( max(ixl,2) .le. i1 ) .and. ( i1 .le. ixr ) .and. & + ( max(jyl,2) .le. j1 ) .and. ( j1 .le. jyr ) .and. & + ( max(kzl,2) .le. k1 ) .and. ( k1 .le. kzr ) ) then + tm_mask(i1,j1,k1,l) = 0 + ftmp(i1,j1,k1,l) = f(i1+1,j1+1,k1+1,l) + sqrt(three)*delta + end if end if - end if - end if - if ( ( imax_n(1) .ne. imax_glob(1) ) .and. & - ( imax_n(3) .ne. imax_glob(3) ) ) then - i2 = imax_n(1) - cctk_lbnd(1) - k2 = imax_n(3) - cctk_lbnd(3) - if ( ( ixl .le. i2 ) .and. ( i2 .le. min(ixr,nx-1) ) .and. & - ( kzl .le. k2 ) .and. ( k2 .le. min(kzr,nz-1) ) ) then - 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 ) - if ( j1 .le. j2 ) then - tm_mask(i2,j1:j2,k2) = 0 - ftmp(i2,j1:j2,k2) = f(i2-1,j1:j2,k2-1) + sqrt(two)*delta + if ( ( imin_n(1) .ne. imin_glob(1) ) .and. & + ( imin_n(2) .ne. imin_glob(2) ) .and. & + ( imax_n(3) .ne. imax_glob(3) ) ) then + i1 = imin_n(1) - cctk_lbnd(1) + j1 = imin_n(2) - cctk_lbnd(2) + k2 = imax_n(3) - cctk_lbnd(3) + if ( ( max(ixl,2) .le. i1 ) .and. ( i1 .le. ixr ) .and. & + ( max(jyl,2) .le. j1 ) .and. ( j1 .le. jyr ) .and. & + ( kzl .le. k2 ) .and. ( k2 .le. min(kzr,nz-1) ) ) then + tm_mask(i1,j1,k2,l) = 0 + ftmp(i1,j1,k2,l) = f(i1+1,j1+1,k2-1,l) + sqrt(three)*delta + end if end if - end if - end if - if ( ( imin_n(2) .ne. imin_glob(2) ) .and. & - ( imin_n(3) .ne. imin_glob(3) ) ) then - j1 = imin_n(2) - cctk_lbnd(2) - k1 = imin_n(3) - cctk_lbnd(3) - if ( ( max(jyl,2) .le. j1 ) .and. ( j1 .le. jyr ) .and. & - ( max(kzl,2) .le. k1 ) .and. ( k1 .le. kzr ) ) then - i1 = max ( max ( imin_n(1), imin_glob(1) ) - cctk_lbnd(1), ixl ) - i2 = min ( min ( imax_n(1), imax_glob(1) ) - cctk_lbnd(1), ixr ) - if ( i1 .le. i2 ) then - tm_mask(i1:i2,j1,k1) = 0 - ftmp(i1:i2,j1,k1) = f(i1:i2,j1+1,k1+1) + sqrt(two)*delta + if ( ( imin_n(1) .ne. imin_glob(1) ) .and. & + ( imax_n(2) .ne. imax_glob(2) ) .and. & + ( imin_n(3) .ne. imin_glob(3) ) ) then + i1 = imin_n(1) - cctk_lbnd(1) + j2 = imax_n(2) - cctk_lbnd(2) + k1 = imin_n(3) - cctk_lbnd(3) + if ( ( max(ixl,2) .le. i1 ) .and. ( i1 .le. ixr ) .and. & + ( jyl .le. j2 ) .and. ( j2 .le. min(jyr,ny-1) ) .and. & + ( max(kzl,2) .le. k1 ) .and. ( k1 .le. kzr ) ) then + tm_mask(i1,j2,k1,l) = 0 + ftmp(i1,j2,k1,l) = f(i1+1,j2-1,k1+1,l) + sqrt(three)*delta + end if end if - end if - end if - if ( ( imin_n(2) .ne. imin_glob(2) ) .and. & - ( imax_n(3) .ne. imax_glob(3) ) ) then - j1 = imin_n(2) - cctk_lbnd(2) - k2 = imax_n(3) - cctk_lbnd(3) - if ( ( max(jyl,2) .le. j1 ) .and. ( j1 .le. jyr ) .and. & - ( kzl .le. k2 ) .and. ( k2 .le. min(kzr,nz-1) ) ) then - i1 = max ( max ( imin_n(1), imin_glob(1) ) - cctk_lbnd(1), ixl ) - i2 = min ( min ( imax_n(1), imax_glob(1) ) - cctk_lbnd(1), ixr ) - if ( i1 .le. i2 ) then - tm_mask(i1:i2,j1,k2) = 0 - ftmp(i1:i2,j1,k2) = f(i1:i2,j1+1,k2-1) + sqrt(two)*delta + if ( ( imin_n(1) .ne. imin_glob(1) ) .and. & + ( imax_n(2) .ne. imax_glob(2) ) .and. & + ( imax_n(3) .ne. imax_glob(3) ) ) then + i1 = imin_n(1) - cctk_lbnd(1) + j2 = imax_n(2) - cctk_lbnd(2) + k2 = imax_n(3) - cctk_lbnd(3) + if ( ( max(ixl,2) .le. i1 ) .and. ( i1 .le. ixr ) .and. & + ( jyl .le. j2 ) .and. ( j2 .le. min(jyr,ny-1) ) .and. & + ( kzl .le. k2 ) .and. ( k2 .le. min(kzr,nz-1) ) ) then + tm_mask(i1,j2,k2,l) = 0 + ftmp(i1,j2,k2,l) = f(i1+1,j2-1,k2-1,l) + sqrt(three)*delta + end if end if - end if - end if - if ( ( imax_n(2) .ne. imax_glob(2) ) .and. & - ( imin_n(3) .ne. imin_glob(3) ) ) then - j2 = imax_n(2) - cctk_lbnd(2) - k1 = imin_n(3) - cctk_lbnd(3) - if ( ( jyl .le. j2 ) .and. ( j2 .le. min(jyr,ny-1) ) .and. & - ( max(kzl,2) .le. k1 ) .and. ( k1 .le. kzr ) ) then - i1 = max ( max ( imin_n(1), imin_glob(1) ) - cctk_lbnd(1), ixl ) - i2 = min ( min ( imax_n(1), imax_glob(1) ) - cctk_lbnd(1), ixr ) - if ( i1 .le. i2 ) then - tm_mask(i1:i2,j2,k1) = 0 - ftmp(i1:i2,j2,k1) = f(i1:i2,j2-1,k1+1) + sqrt(two)*delta + if ( ( imax_n(1) .ne. imax_glob(1) ) .and. & + ( imin_n(2) .ne. imin_glob(2) ) .and. & + ( imin_n(3) .ne. imin_glob(3) ) ) then + i2 = imax_n(1) - cctk_lbnd(1) + j1 = imin_n(2) - cctk_lbnd(2) + k1 = imin_n(3) - cctk_lbnd(3) + if ( ( ixl .le. i2 ) .and. ( i2 .le. min(ixr,nx-1) ) .and. & + ( max(jyl,2) .le. j1 ) .and. ( j1 .le. jyr ) .and. & + ( max(kzl,2) .le. k1 ) .and. ( k1 .le. kzr ) ) then + tm_mask(i2,j1,k1,l) = 0 + ftmp(i2,j1,k1,l) = f(i2-1,j1+1,k1+1,l) + sqrt(three)*delta + end if end if - end if - end if - if ( ( imax_n(2) .ne. imax_glob(2) ) .and. & - ( imax_n(3) .ne. imax_glob(3) ) ) then - j2 = imax_n(2) - cctk_lbnd(2) - k2 = imax_n(3) - cctk_lbnd(3) - if ( ( jyl .le. j2 ) .and. ( j2 .le. min(jyr,ny-1) ) .and. & - ( kzl .le. k2 ) .and. ( k2 .le. min(kzr,nz-1) ) ) then - i1 = max ( max ( imin_n(1), imin_glob(1) ) - cctk_lbnd(1), ixl ) - i2 = min ( min ( imax_n(1), imax_glob(1) ) - cctk_lbnd(1), ixr ) - if ( i1 .le. i2 ) then - tm_mask(i1:i2,j2,k2) = 0 - ftmp(i1:i2,j2,k2) = f(i1:i2,j2-1,k2-1) + sqrt(two)*delta + if ( ( imax_n(1) .ne. imax_glob(1) ) .and. & + ( imin_n(2) .ne. imin_glob(2) ) .and. & + ( imax_n(3) .ne. imax_glob(3) ) ) then + i2 = imax_n(1) - cctk_lbnd(1) + j1 = imin_n(2) - cctk_lbnd(2) + k2 = imax_n(3) - cctk_lbnd(3) + if ( ( ixl .le. i2 ) .and. ( i2 .le. min(ixr,nx-1) ) .and. & + ( max(jyl,2) .le. j1 ) .and. ( j1 .le. jyr ) .and. & + ( kzl .le. k2 ) .and. ( k2 .le. min(kzr,nz-1) ) ) then + tm_mask(i2,j1,k2,l) = 0 + ftmp(i2,j1,k2,l) = f(i2-1,j1+1,k2-1,l) + sqrt(three)*delta + end if + end if + if ( ( imax_n(1) .ne. imax_glob(1) ) .and. & + ( imax_n(2) .ne. imax_glob(2) ) .and. & + ( imin_n(3) .ne. imin_glob(3) ) ) then + i2 = imax_n(1) - cctk_lbnd(1) + j2 = imax_n(2) - cctk_lbnd(2) + k1 = imin_n(3) - cctk_lbnd(3) + if ( ( ixl .le. i2 ) .and. ( i2 .le. min(ixr,nx-1) ) .and. & + ( jyl .le. j2 ) .and. ( j2 .le. min(jyr,ny-1) ) .and. & + ( max(kzl,2) .le. k1 ) .and. ( k1 .le. kzr ) ) then + tm_mask(i2,j2,k1,l) = 0 + ftmp(i2,j2,k1,l) = f(i2-1,j2-1,k1+1,l) + sqrt(three)*delta + end if + end if + if ( ( imax_n(1) .ne. imax_glob(1) ) .and. & + ( imax_n(2) .ne. imax_glob(2) ) .and. & + ( imax_n(3) .ne. imax_glob(3) ) ) then + 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 + tm_mask(i2,j2,k2,l) = 0 + ftmp(i2,j2,k2,l) = f(i2-1,j2-1,k2-1,l) + sqrt(three)*delta + end if end if end if - end if -! And finally do the corners. - if ( ( imin_n(1) .ne. imin_glob(1) ) .and. & - ( imin_n(2) .ne. imin_glob(2) ) .and. & - ( imin_n(3) .ne. imin_glob(3) ) ) then - i1 = imin_n(1) - cctk_lbnd(1) - j1 = imin_n(2) - cctk_lbnd(2) - k1 = imin_n(3) - cctk_lbnd(3) - if ( ( max(ixl,2) .le. i1 ) .and. ( i1 .le. ixr ) .and. & - ( max(jyl,2) .le. j1 ) .and. ( j1 .le. jyr ) .and. & - ( max(kzl,2) .le. k1 ) .and. ( k1 .le. kzr ) ) then - tm_mask(i1,j1,k1) = 0 - ftmp(i1,j1,k1) = f(i1+1,j1+1,k1+1) + sqrt(three)*delta - end if - end if - if ( ( imin_n(1) .ne. imin_glob(1) ) .and. & - ( imin_n(2) .ne. imin_glob(2) ) .and. & - ( imax_n(3) .ne. imax_glob(3) ) ) then - i1 = imin_n(1) - cctk_lbnd(1) - j1 = imin_n(2) - cctk_lbnd(2) - k2 = imax_n(3) - cctk_lbnd(3) - if ( ( max(ixl,2) .le. i1 ) .and. ( i1 .le. ixr ) .and. & - ( max(jyl,2) .le. j1 ) .and. ( j1 .le. jyr ) .and. & - ( kzl .le. k2 ) .and. ( k2 .le. min(kzr,nz-1) ) ) then - tm_mask(i1,j1,k2) = 0 - ftmp(i1,j1,k2) = f(i1+1,j1+1,k2-1) + sqrt(three)*delta - end if - end if - if ( ( imin_n(1) .ne. imin_glob(1) ) .and. & - ( imax_n(2) .ne. imax_glob(2) ) .and. & - ( imin_n(3) .ne. imin_glob(3) ) ) then - i1 = imin_n(1) - cctk_lbnd(1) - j2 = imax_n(2) - cctk_lbnd(2) - k1 = imin_n(3) - cctk_lbnd(3) - if ( ( max(ixl,2) .le. i1 ) .and. ( i1 .le. ixr ) .and. & - ( jyl .le. j2 ) .and. ( j2 .le. min(jyr,ny-1) ) .and. & - ( max(kzl,2) .le. k1 ) .and. ( k1 .le. kzr ) ) then - tm_mask(i1,j2,k1) = 0 - ftmp(i1,j2,k1) = f(i1+1,j2-1,k1+1) + sqrt(three)*delta - end if - end if - if ( ( imin_n(1) .ne. imin_glob(1) ) .and. & - ( imax_n(2) .ne. imax_glob(2) ) .and. & - ( imax_n(3) .ne. imax_glob(3) ) ) then - i1 = imin_n(1) - cctk_lbnd(1) - j2 = imax_n(2) - cctk_lbnd(2) - k2 = imax_n(3) - cctk_lbnd(3) - if ( ( max(ixl,2) .le. i1 ) .and. ( i1 .le. ixr ) .and. & - ( jyl .le. j2 ) .and. ( j2 .le. min(jyr,ny-1) ) .and. & - ( kzl .le. k2 ) .and. ( k2 .le. min(kzr,nz-1) ) ) then - tm_mask(i1,j2,k2) = 0 - ftmp(i1,j2,k2) = f(i1+1,j2-1,k2-1) + sqrt(three)*delta +! Copy the modified mask and level set function into the proper place. + eh_mask(ixl:ixr,jyl:jyr,kzl:kzr,l) = tm_mask(ixl:ixr,jyl:jyr,kzl:kzr,l) + f(ixl:ixr,jyl:jyr,kzl:kzr,l) = ftmp(ixl:ixr,jyl:jyr,kzl:kzr,l) + +! Now check if any interior points are far enough away from the f=0 +! surface and if so excise them. + if ( use_inner_excision .gt. 0 ) then + where ( ( eh_mask(ixl:ixr,jyl:jyr,kzl:kzr,l) .ge. 0 ) .and. & + ( f(ixl:ixr,jyl:jyr,kzl:kzr,l) .lt. -shell_width*delta ) ) + eh_mask(ixl:ixr,jyl:jyr,kzl:kzr,l) = -1 + f(ixl:ixr,jyl:jyr,kzl:kzr,l) = ex_value + end where end if - end if - if ( ( imax_n(1) .ne. imax_glob(1) ) .and. & - ( imin_n(2) .ne. imin_glob(2) ) .and. & - ( imin_n(3) .ne. imin_glob(3) ) ) then - i2 = imax_n(1) - cctk_lbnd(1) - j1 = imin_n(2) - cctk_lbnd(2) - k1 = imin_n(3) - cctk_lbnd(3) - if ( ( ixl .le. i2 ) .and. ( i2 .le. min(ixr,nx-1) ) .and. & - ( max(jyl,2) .le. j1 ) .and. ( j1 .le. jyr ) .and. & - ( max(kzl,2) .le. k1 ) .and. ( k1 .le. kzr ) ) then - tm_mask(i2,j1,k1) = 0 - ftmp(i2,j1,k1) = f(i2-1,j1+1,k1+1) + sqrt(three)*delta - end if - end if - if ( ( imax_n(1) .ne. imax_glob(1) ) .and. & - ( imin_n(2) .ne. imin_glob(2) ) .and. & - ( imax_n(3) .ne. imax_glob(3) ) ) then - i2 = imax_n(1) - cctk_lbnd(1) - j1 = imin_n(2) - cctk_lbnd(2) - k2 = imax_n(3) - cctk_lbnd(3) - if ( ( ixl .le. i2 ) .and. ( i2 .le. min(ixr,nx-1) ) .and. & - ( max(jyl,2) .le. j1 ) .and. ( j1 .le. jyr ) .and. & - ( kzl .le. k2 ) .and. ( k2 .le. min(kzr,nz-1) ) ) then - tm_mask(i2,j1,k2) = 0 - ftmp(i2,j1,k2) = f(i2-1,j1+1,k2-1) + sqrt(three)*delta - end if - end if - if ( ( imax_n(1) .ne. imax_glob(1) ) .and. & - ( imax_n(2) .ne. imax_glob(2) ) .and. & - ( imin_n(3) .ne. imin_glob(3) ) ) then - i2 = imax_n(1) - cctk_lbnd(1) - j2 = imax_n(2) - cctk_lbnd(2) - k1 = imin_n(3) - cctk_lbnd(3) - if ( ( ixl .le. i2 ) .and. ( i2 .le. min(ixr,nx-1) ) .and. & - ( jyl .le. j2 ) .and. ( j2 .le. min(jyr,ny-1) ) .and. & - ( max(kzl,2) .le. k1 ) .and. ( k1 .le. kzr ) ) then - tm_mask(i2,j2,k1) = 0 - ftmp(i2,j2,k1) = f(i2-1,j2-1,k1+1) + sqrt(three)*delta - end if - end if - if ( ( imax_n(1) .ne. imax_glob(1) ) .and. & - ( imax_n(2) .ne. imax_glob(2) ) .and. & - ( imax_n(3) .ne. imax_glob(3) ) ) then - 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 - tm_mask(i2,j2,k2) = 0 - ftmp(i2,j2,k2) = f(i2-1,j2-1,k2-1) + sqrt(three)*delta - end if - end if - end if -! print*,'Debug 4' -! print*,ftmp(imin_n(1),imin_n(2),imin_n(3)) -! print*,ftmp(imin_n(1),imin_n(2),imax_n(3)) -! print*,ftmp(imin_n(1),imax_n(2),imin_n(3)) -! print*,ftmp(imin_n(1),imax_n(2),imax_n(3)) -! print*,ftmp(imax_n(1),imin_n(2),imin_n(3)) -! print*,ftmp(imax_n(1),imin_n(2),imax_n(3)) -! print*,ftmp(imax_n(1),imax_n(2),imin_n(3)) -! print*,ftmp(imax_n(1),imax_n(2),imax_n(3)) - -! Copy the modified mask and level set function into the proper place. - eh_mask(ixl:ixr,jyl:jyr,kzl:kzr) = tm_mask(ixl:ixr,jyl:jyr,kzl:kzr) - f(ixl:ixr,jyl:jyr,kzl:kzr) = ftmp(ixl:ixr,jyl:jyr,kzl:kzr) - -! Now check if any interior points are far enough away from the f=0 -! surface and if so excise them. - if ( use_inner_excision .gt. 0 ) then - where ( ( eh_mask(ixl:ixr,jyl:jyr,kzl:kzr) .ge. 0 ) .and. & - ( f(ixl:ixr,jyl:jyr,kzl:kzr) .lt. -shell_width*delta ) ) - eh_mask(ixl:ixr,jyl:jyr,kzl:kzr) = -1 - f(ixl:ixr,jyl:jyr,kzl:kzr) = ex_value - end where - end if - -! Now check if any remaining active points where actually excised in the -! numerical run generating the metric data. - if ( use_mask .gt. 0 ) then - where ( emask(ixl:ixr,jyl:jyr,kzl:kzr) .lt. three*quarter ) - eh_mask(ixl:ixr,jyl:jyr,kzl:kzr) = -1 - f(ixl:ixr,jyl:jyr,kzl:kzr) = ex_value - end where - end if +! Now check if any remaining active points where actually excised in the +! numerical run generating the metric data. + + if ( use_mask .gt. 0 ) then + where ( emask(ixl:ixr,jyl:jyr,kzl:kzr) .lt. three*quarter ) + eh_mask(ixl:ixr,jyl:jyr,kzl:kzr,l) = -1 + f(ixl:ixr,jyl:jyr,kzl:kzr,l) = ex_value + end where + end if -! Make sure to mark all points outside of the rectangular box as excised -! points and set the value of f to -ex_value. - if ( use_outer_excision .gt. 0 ) then - do k = kzl, kzr - do j = jyl, jyr - do i = ixl, ixr - if ( ( ( i + cctk_lbnd(1) .lt. imin_n(1) ) .or. & - ( i + cctk_lbnd(1) .gt. imax_n(1) ) .or. & - ( j + cctk_lbnd(2) .lt. imin_n(2) ) .or. & - ( j + cctk_lbnd(2) .gt. imax_n(2) ) .or. & - ( k + cctk_lbnd(3) .lt. imin_n(3) ) .or. & - ( k + cctk_lbnd(3) .gt. imax_n(3) ) ) .and. & - ( eh_mask(i,j,k) .ge. 0 ) ) then - eh_mask(i,j,k) = -1 - f(i,j,k) = -ex_value - end if +! Make sure to mark all points outside of the rectangular box as excised +! points and set the value of f to -ex_value. + if ( use_outer_excision .gt. 0 ) then + do k = kzl, kzr + do j = jyl, jyr + do i = ixl, ixr + if ( ( ( i + cctk_lbnd(1) .lt. imin_n(1) ) .or. & + ( i + cctk_lbnd(1) .gt. imax_n(1) ) .or. & + ( j + cctk_lbnd(2) .lt. imin_n(2) ) .or. & + ( j + cctk_lbnd(2) .gt. imax_n(2) ) .or. & + ( k + cctk_lbnd(3) .lt. imin_n(3) ) .or. & + ( k + cctk_lbnd(3) .gt. imax_n(3) ) ) .and. & + ( eh_mask(i,j,k,l) .ge. 0 ) ) then + eh_mask(i,j,k,l) = -1 + f(i,j,k,l) = -ex_value + end if + end do + end do end do - end do - end do - end if + end if + end if not_undone + end do loop_over_l end if end subroutine EHFinder_SetMask1 @@ -760,7 +757,7 @@ subroutine EHFinder_SetMask2(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS - CCTK_INT :: i, j, k + CCTK_INT :: i, j, k, l logical :: active active = .false. @@ -784,7 +781,7 @@ subroutine EHFinder_SetMask2(CCTK_ARGUMENTS) end if ! If the reparametrization was not undone... - if ( active .and. .not.reparam_undone ) then + if ( active .and. .not. all(reparam_undone) ) then ! Get the minimum and maximum index excluding ghost and symmetry cells. # include "include/physical_part.h" @@ -797,61 +794,65 @@ subroutine EHFinder_SetMask2(CCTK_ARGUMENTS) if ( kzl .eq. 1 ) kzl = 2 if ( kzr .eq. nz ) kzr = nz - 1 -! Initialize the temporary mask to zero. - tm_mask(ixl:ixr,jyl:jyr,kzl:kzr) = 0 + do l = 1, eh_number_level_sets +! Initialize the temporary mask to zero. + tm_mask(ixl:ixr,jyl:jyr,kzl:kzr,l) = 0 -! We loop over all points and adjust the mask if the point is -! on the excision boundary. - do k = kzl, kzr - do j = jyl, jyr - do i = ixl, ixr + if ( .not. reparam_undone(l) ) then -! If the point is active..... - if ( eh_mask(i,j,k) .ge. 0 ) then +! We loop over all points and adjust the mask if the point is +! on the excision boundary. + do k = kzl, kzr + do j = jyl, jyr + do i = ixl, ixr -! If the neighbour in the -x-directions is excised.... - if ( eh_mask(i-1,j,k) .eq. -1 ) then - tm_mask(i,j,k) = tm_mask(i,j,k) + ll(0) - end if +! If the point is active..... + if ( eh_mask(i,j,k,l) .ge. 0 ) then -! If the neighbour in the +x-directions is excised.... - if ( eh_mask(i+1,j,k) .eq. -1 ) then - tm_mask(i,j,k) = tm_mask(i,j,k) + ll(1) - end if +! If the neighbour in the -x-directions is excised.... + if ( eh_mask(i-1,j,k,l) .eq. -1 ) then + tm_mask(i,j,k,l) = tm_mask(i,j,k,l) + ll(0) + end if -! If the neighbour in the -y-directions is excised.... - if ( eh_mask(i,j-1,k) .eq. -1 ) then - tm_mask(i,j,k) = tm_mask(i,j,k) + ll(2) - end if +! If the neighbour in the +x-directions is excised.... + if ( eh_mask(i+1,j,k,l) .eq. -1 ) then + tm_mask(i,j,k,l) = tm_mask(i,j,k,l) + ll(1) + end if -! If the neighbour in the +y-directions is excised.... - if ( eh_mask(i,j+1,k) .eq. -1 ) then - tm_mask(i,j,k) = tm_mask(i,j,k) + ll(3) - end if +! If the neighbour in the -y-directions is excised.... + if ( eh_mask(i,j-1,k,l) .eq. -1 ) then + tm_mask(i,j,k,l) = tm_mask(i,j,k,l) + ll(2) + end if -! If the neighbour in the -z-directions is excised.... - if ( eh_mask(i,j,k-1) .eq. -1 ) then - tm_mask(i,j,k) = tm_mask(i,j,k) + ll(4) - end if +! If the neighbour in the +y-directions is excised.... + if ( eh_mask(i,j+1,k,l) .eq. -1 ) then + tm_mask(i,j,k,l) = tm_mask(i,j,k,l) + ll(3) + end if -! If the neighbour in the +z-directions is excised.... - if ( eh_mask(i,j,k+1) .eq. -1 ) then - tm_mask(i,j,k) = tm_mask(i,j,k) + ll(5) - end if +! If the neighbour in the -z-directions is excised.... + if ( eh_mask(i,j,k-1,l) .eq. -1 ) then + tm_mask(i,j,k,l) = tm_mask(i,j,k,l) + ll(4) + end if - end if +! If the neighbour in the +z-directions is excised.... + if ( eh_mask(i,j,k+1,l) .eq. -1 ) then + tm_mask(i,j,k,l) = tm_mask(i,j,k,l) + ll(5) + end if - end do - end do - end do + end if -! Copy the changes back into eh_mask - where ( eh_mask(ixl:ixr,jyl:jyr,kzl:kzr) .ge. 0 ) - eh_mask(ixl:ixr,jyl:jyr,kzl:kzr) = & - tm_mask(ixl:ixr,jyl:jyr,kzl:kzr) - end where + end do + end do + end do - call CCTK_INFO('Mask modified') +! Copy the changes back into eh_mask + where ( eh_mask(ixl:ixr,jyl:jyr,kzl:kzr,l) .ge. 0 ) + eh_mask(ixl:ixr,jyl:jyr,kzl:kzr,l) = & + 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. mask_first = .false. @@ -870,7 +871,7 @@ subroutine EHFinder_SetMask3(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS - CCTK_INT :: i, j, k + CCTK_INT :: i, j, k, l logical active, mask_modified active = .false. @@ -894,7 +895,7 @@ subroutine EHFinder_SetMask3(CCTK_ARGUMENTS) end if ! If the reparametrization was not undone... - if ( active .and. .not.reparam_undone ) then + if ( active .and. .not. all(reparam_undone) ) then mask_modified = .false. # include "include/physical_part.h" @@ -907,98 +908,103 @@ subroutine EHFinder_SetMask3(CCTK_ARGUMENTS) if ( kzl .eq. 1 ) kzl = 2 if ( kzr .eq. nz ) kzr = nz - 1 - tm_mask = eh_mask - - do k = kzl, kzr - do j = jyl, jyr - do i = ixl, ixr - if ( eh_mask(i,j,k) .ge. 0 ) then - ! If we are on an excision boundary in the -x-direction... - if ( btest(eh_mask(i,j,k),0) ) then - ! 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) .eq. -1 ) .or. & - ( eh_mask(i+2,j,k) .eq. -1 ) ) then - tm_mask(i,j,k) = -1 - mask_modified = .true. - print*, i, j, k, 'x-' - print*, eh_mask(i:i+2,j,k) - end if - end if + do l = 1, eh_number_level_sets + tm_mask(:,:,:,l) = eh_mask(:,:,:,l) + + if ( .not. reparam_undone(l) ) then + + do k = kzl, kzr + 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 ( 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 ( ( 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 ( btest(eh_mask(i,j,k),1) ) then - ! 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) .eq. -1 ) .or. & - ( eh_mask(i-2,j,k) .eq. -1 ) ) then - tm_mask(i,j,k) = -1 - mask_modified = .true. - print*, i, j, k, 'x+' - print*, eh_mask(i:i-2,j,k) - end if - end if + ! 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 ( ( 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 ( btest(eh_mask(i,j,k),2) ) then - ! 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) .eq. -1 ) .or. & - ( eh_mask(i,j+2,k) .eq. -1 ) ) then - tm_mask(i,j,k) = -1 - mask_modified = .true. - print*, i, j, k, 'y-' - print*, eh_mask(i,j:j+2,k) - end if - end if + ! 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 ( ( 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 ( btest(eh_mask(i,j,k),3) ) then - ! 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) .eq. -1 ) .or. & - ( eh_mask(i,j-2,k) .eq. -1 ) ) then - tm_mask(i,j,k) = -1 - mask_modified = .true. - print*, i, j, k, 'y+' - print*, eh_mask(i,j:j-2,k) - end if - end if + ! 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 ( ( 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 ( btest(eh_mask(i,j,k),4) ) then - ! 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) .eq. -1 ) .or. & - ( eh_mask(i,j,k+2) .eq. -1 ) ) then - tm_mask(i,j,k) = -1 - mask_modified = .true. - print*, i, j, k, 'z-' - print*, eh_mask(i,j,k:k+2) - end if - end if + ! 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 ( ( 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 ( btest(eh_mask(i,j,k),5) ) then - ! 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) .eq. -1 ) .or. & - ( eh_mask(i,j,k-2) .eq. -1 ) ) then - tm_mask(i,j,k) = -1 - mask_modified = .true. - print*, i, j, k, 'z+' - print*, eh_mask(i,j,k:k-2) + ! 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 ( ( 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 - end if - end if + end do + end do end do - end do - end do - if ( mask_modified ) then - call CCTK_WARN ( 1, "Mask modified because it was not convex" ) - end if + if ( mask_modified ) then + call CCTK_WARN ( 1, "Mask modified because it was not convex" ) + end if + end if + end do eh_mask = tm_mask end if diff --git a/src/EHFinder_SetSym.F90 b/src/EHFinder_SetSym.F90 index 817ec2e..b04672d 100644 --- a/src/EHFinder_SetSym.F90 +++ b/src/EHFinder_SetSym.F90 @@ -15,15 +15,24 @@ subroutine EHFinder_SetSym(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS CCTK_INT, dimension(3) :: sym + CCTK_INT :: l + character(len=14) :: fname sym = 1 - call SetCartSymVN(ierr,cctkGH,sym,'ehfinder::f') - if ( ierr .gt. 0 ) then - call CCTK_WARN(1,'Failed to register symmetry for the level set function') - end if - - call SetCartSymVN(ierr,cctkGH,sym,'ehfinder::sc') + do l = 0, eh_number_level_sets - 1 + write(fname,'(a12,i1,a1)') 'ehfinder::f[', l, ']' + call SetCartSymVN(ierr,cctkGH,sym,fname) + if ( ierr .gt. 0 ) then + call CCTK_WARN(1,'Failed to register symmetry for the level set function') + end if +! call SetCartSymVN(ierr,cctkGH,sym,'ehfinder::f[2]') +! if ( ierr .gt. 0 ) then +! call CCTK_WARN(1,'Failed to register symmetry for the level set function') +! end if + end do + + call SetCartSymGN(ierr,cctkGH,sym,'ehfinder::surface_index') if ( ierr .gt. 0 ) then call CCTK_WARN(1,'Failed to register symmetry for sc') end if @@ -47,21 +56,27 @@ subroutine EHFinder_ApplySymAll(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS - call CartSymVN(ierr,cctkGH,'ehfinder::f') + CCTK_INT :: l + character(len=14) :: fname + + do l = 0, eh_number_level_sets - 1 + write(fname,'(a12,i1,a1)') 'ehfinder::f[', l, ']' + call CartSymVN(ierr,cctkGH,fname) + end do # include "include/physical_part.h" if ( symx ) then - eh_mask(ngx:1:-1,:,:) = eh_mask(ngx+1:ngx+ngx,:,:) - rep_mask(ngx:1:-1,:,:) = rep_mask(ngx+1:ngx+ngx,:,:) + eh_mask(ngx:1:-1,:,:,:) = eh_mask(ngx+1:ngx+ngx,:,:,:) +! rep_mask(ngx:1:-1,:,:,:) = rep_mask(ngx+1:ngx+ngx,:,:,:) end if if ( symy ) then - eh_mask(:,ngy:1:-1,:) = eh_mask(:,ngy+1:ngy+ngy,:) - rep_mask(:,ngy:1:-1,:) = rep_mask(:,ngy+1:ngy+ngy,:) + eh_mask(:,ngy:1:-1,:,:) = eh_mask(:,ngy+1:ngy+ngy,:,:) +! rep_mask(:,ngy:1:-1,:,:) = rep_mask(:,ngy+1:ngy+ngy,:,:) end if if ( symz ) then - eh_mask(:,:,ngz:1:-1) = eh_mask(:,:,ngz+1:ngz+ngz) - rep_mask(:,:,ngz:1:-1) = rep_mask(:,:,ngz+1:ngz+ngz) + eh_mask(:,:,ngz:1:-1,:) = eh_mask(:,:,ngz+1:ngz+ngz,:) +! rep_mask(:,:,ngz:1:-1,:) = rep_mask(:,:,ngz+1:ngz+ngz,:) end if return end subroutine EHFinder_ApplySymAll @@ -77,7 +92,14 @@ subroutine EHFinder_ApplySymF(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS - call CartSymVN(ierr,cctkGH,'ehfinder::f') + CCTK_INT :: l + character(len=14) :: fname + + do l = 0, eh_number_level_sets - 1 + write(fname,'(a12,i1,a1)') 'ehfinder::f[', l, ']' + call CartSymVN(ierr,cctkGH,fname) + end do + return end subroutine EHFinder_ApplySymF @@ -94,15 +116,15 @@ subroutine EHFinder_ApplySymRep(CCTK_ARGUMENTS) # include "include/physical_part.h" - if ( symx ) then - rep_mask(ngx:1:-1,:,:) = rep_mask(ngx+1:ngx+ngx,:,:) - end if - if ( symy ) then - rep_mask(:,ngy:1:-1,:) = rep_mask(:,ngy+1:ngy+ngy,:) - end if - if ( symz ) then - rep_mask(:,:,ngz:1:-1) = rep_mask(:,:,ngz+1:ngz+ngz) - end if +! if ( symx ) then +! rep_mask(ngx:1:-1,:,:,:) = rep_mask(ngx+1:ngx+ngx,:,:,:) +! end if +! if ( symy ) then +! rep_mask(:,ngy:1:-1,:,:) = rep_mask(:,ngy+1:ngy+ngy,:,:) +! end if +! if ( symz ) then +! rep_mask(:,:,ngz:1:-1,:) = rep_mask(:,:,ngz+1:ngz+ngz,:) +! end if return end subroutine EHFinder_ApplySymRep @@ -117,19 +139,25 @@ subroutine EHFinder_ApplySymFRep(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS + CCTK_INT :: l + character(len=14) :: fname + # include "include/physical_part.h" - call CartSymVN(ierr,cctkGH,'ehfinder::f') + do l = 0, eh_number_level_sets - 1 + write(fname,'(a12,i1,a1)') 'ehfinder::f[', l, ']' + call CartSymVN(ierr,cctkGH,fname) + end do - if ( symx ) then - rep_mask(ngx:1:-1,:,:) = rep_mask(ngx+1:ngx+ngx,:,:) - end if - if ( symy ) then - rep_mask(:,ngy:1:-1,:) = rep_mask(:,ngy+1:ngy+ngy,:) - end if - if ( symz ) then - rep_mask(:,:,ngz:1:-1) = rep_mask(:,:,ngz+1:ngz+ngz) - end if +! if ( symx ) then +! rep_mask(ngx:1:-1,:,:,:) = rep_mask(ngx+1:ngx+ngx,:,:,:) +! end if +! if ( symy ) then +! rep_mask(:,ngy:1:-1,:,:) = rep_mask(:,ngy+1:ngy+ngy,:,:) +! end if +! if ( symz ) then +! rep_mask(:,:,ngz:1:-1,:) = rep_mask(:,:,ngz+1:ngz+ngz,:) +! end if return end subroutine EHFinder_ApplySymFRep @@ -147,13 +175,13 @@ subroutine EHFinder_ApplySymMask(CCTK_ARGUMENTS) # include "include/physical_part.h" if ( symx ) then - eh_mask(ngx:1:-1,:,:) = eh_mask(ngx+1:ngx+ngx,:,:) + eh_mask(ngx:1:-1,:,:,:) = eh_mask(ngx+1:ngx+ngx,:,:,:) end if if ( symy ) then - eh_mask(:,ngy:1:-1,:) = eh_mask(:,ngy+1:ngy+ngy,:) + eh_mask(:,ngy:1:-1,:,:) = eh_mask(:,ngy+1:ngy+ngy,:,:) end if if ( symz ) then - eh_mask(:,:,ngz:1:-1) = eh_mask(:,:,ngz+1:ngz+ngz) + eh_mask(:,:,ngz:1:-1,:) = eh_mask(:,:,ngz+1:ngz+ngz,:) end if return end subroutine EHFinder_ApplySymMask diff --git a/src/EHFinder_Sources.F90 b/src/EHFinder_Sources.F90 index cdf2095..65e0488 100644 --- a/src/EHFinder_Sources.F90 +++ b/src/EHFinder_Sources.F90 @@ -15,7 +15,7 @@ subroutine EHFinder_Sources(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS - CCTK_INT :: i, j, k + CCTK_INT :: i, j, k, l CCTK_REAL :: idx, idy, idz, mdelta CCTK_REAL :: a, b, c CCTK_REAL :: g3xx, g3xy, g3xz, g3yy, g3yz, g3zz @@ -41,132 +41,142 @@ subroutine EHFinder_Sources(CCTK_ARGUMENTS) if ( CCTK_EQUALS ( surface_direction, 'inward' ) ) ssign = -one if ( CCTK_EQUALS ( upwind_type, 'intrinsic' ) ) then - do k = kzl, kzr - do j = jyl, jyr - do i = ixl, ixr - if ( f(i,j,k) .gt. fmin_bound - three*mdelta ) then + do l = 1, eh_number_level_sets + do k = kzl, kzr + do j = jyl, jyr + do i = ixl, ixr + if ( f(i,j,k,l) .gt. fmin_bound - three*mdelta ) then #include "include/centered_second2.h" - else + else #include "include/upwind_second2.h" - end if + end if + end do end do end do end do end if if ( CCTK_EQUALS ( upwind_type, 'shift' ) ) then - do k = kzl, kzr - do j = jyl, jyr - do i = ixl, ixr + do l = 1, eh_number_level_sets + do k = kzl, kzr + do j = jyl, jyr + do i = ixl, ixr #include "include/upwind_shift_second2.h" + end do end do end do end do end if if ( CCTK_EQUALS ( upwind_type, 'characteristic' ) ) then - do k = kzl, kzr - do j = jyl, jyr - do i = ixl, ixr + do l = 1, eh_number_level_sets + do k = kzl, kzr + do j = jyl, jyr + do i = ixl, ixr #include "include/metric.h" #include "include/upwind_characteristic_second2.h" + end do end do end do end do end if if ( CCTK_EQUALS ( metric_type, 'physical' ) ) then - do k = kzl, kzr - do j = jyl, jyr - do i = ixl, ixr - if ( eh_mask(i,j,k) .ge. 0 ) then - alp2 = alp(i,j,k)**2 - - tmp1 = gyy(i,j,k)*gzz(i,j,k) - gyz(i,j,k)**2 - tmp2 = gxz(i,j,k)*gyz(i,j,k) - gxy(i,j,k)*gzz(i,j,k) - tmp3 = gxy(i,j,k)*gyz(i,j,k) - gxz(i,j,k)*gyy(i,j,k) - - idetg = one / ( gxx(i,j,k)*tmp1 + gxy(i,j,k)*tmp2 + gxz(i,j,k)*tmp3 ) + do l = 1, eh_number_level_sets + do k = kzl, kzr + do j = jyl, jyr + do i = ixl, ixr + if ( eh_mask(i,j,k,l) .ge. 0 ) then + alp2 = alp(i,j,k)**2 - g3xx = tmp1 * idetg - g3xy = tmp2 * idetg - g3xz = tmp3 * idetg - g3yy = ( gxx(i,j,k)*gzz(i,j,k) - gxz(i,j,k)**2 ) * idetg - g3yz = ( gxy(i,j,k)*gxz(i,j,k) - gxx(i,j,k)*gyz(i,j,k) ) * idetg - g3zz = ( gxx(i,j,k)*gyy(i,j,k) - gxy(i,j,k)**2 ) * idetg - - tmp1 = betax(i,j,k) * dfx(i,j,k) + & - betay(i,j,k) * dfy(i,j,k) + & - betaz(i,j,k) * dfz(i,j,k) - - tmp2 = g3xx * dfx(i,j,k)**2 + & - g3yy * dfy(i,j,k)**2 + & - g3zz * dfz(i,j,k)**2 + & - two * ( g3xy * dfx(i,j,k) * dfy(i,j,k) + & - g3xz * dfx(i,j,k) * dfz(i,j,k) + & - g3yz * dfy(i,j,k) * dfz(i,j,k) ) - if ( tmp2 .ge. zero ) then - sf(i,j,k) = tmp1 - ssign * sqrt ( alp2 * tmp2 ) + tmp1 = gyy(i,j,k)*gzz(i,j,k) - gyz(i,j,k)**2 + tmp2 = gxz(i,j,k)*gyz(i,j,k) - gxy(i,j,k)*gzz(i,j,k) + tmp3 = gxy(i,j,k)*gyz(i,j,k) - gxz(i,j,k)*gyy(i,j,k) + + idetg = one / ( gxx(i,j,k)*tmp1 + gxy(i,j,k)*tmp2 + gxz(i,j,k)*tmp3 ) + + g3xx = tmp1 * idetg + g3xy = tmp2 * idetg + g3xz = tmp3 * idetg + g3yy = ( gxx(i,j,k)*gzz(i,j,k) - gxz(i,j,k)**2 ) * idetg + g3yz = ( gxy(i,j,k)*gxz(i,j,k) - gxx(i,j,k)*gyz(i,j,k) ) * idetg + g3zz = ( gxx(i,j,k)*gyy(i,j,k) - gxy(i,j,k)**2 ) * idetg + + tmp1 = betax(i,j,k) * dfx(i,j,k,l) + & + betay(i,j,k) * dfy(i,j,k,l) + & + betaz(i,j,k) * dfz(i,j,k,l) + + tmp2 = g3xx * dfx(i,j,k,l)**2 + & + g3yy * dfy(i,j,k,l)**2 + & + g3zz * dfz(i,j,k,l)**2 + & + two * ( g3xy * dfx(i,j,k,l) * dfy(i,j,k,l) + & + g3xz * dfx(i,j,k,l) * dfz(i,j,k,l) + & + g3yz * dfy(i,j,k,l) * dfz(i,j,k,l) ) + if ( tmp2 .ge. zero ) then + sf(i,j,k,l) = tmp1 - ssign * sqrt ( alp2 * tmp2 ) + else + call CCTK_WARN ( 0, "3-metric not positive definite: Stopping" ) + end if else - call CCTK_WARN ( 0, "3-metric not positive definite: Stopping" ) + sf(i,j,k,l) = zero end if - else - sf(i,j,k) = zero - end if + end do end do - end do - end do + end do + end do end if if ( CCTK_EQUALS ( metric_type, 'static conformal' ) ) then - do k = kzl, kzr - do j = jyl, jyr - do i = ixl, ixr - if ( eh_mask(i,j,k) .ge. 0 ) then - alp2 = alp(i,j,k)**2 + do l = 1, eh_number_level_sets + do k = kzl, kzr + do j = jyl, jyr + do i = ixl, ixr + if ( eh_mask(i,j,k,l) .ge. 0 ) then + alp2 = alp(i,j,k)**2 - psito4 = psi(i,j,k)**4 - gxxc = gxx(i,j,k) * psito4 - gxyc = gxy(i,j,k) * psito4 - gxzc = gxz(i,j,k) * psito4 - gyyc = gyy(i,j,k) * psito4 - gyzc = gyz(i,j,k) * psito4 - gzzc = gzz(i,j,k) * psito4 + psito4 = psi(i,j,k)**4 + gxxc = gxx(i,j,k) * psito4 + gxyc = gxy(i,j,k) * psito4 + gxzc = gxz(i,j,k) * psito4 + gyyc = gyy(i,j,k) * psito4 + gyzc = gyz(i,j,k) * psito4 + gzzc = gzz(i,j,k) * psito4 - tmp1 = gyyc*gzzc - gyzc**2 - tmp2 = gxzc*gyzc - gxyc*gzzc - tmp3 = gxyc*gyzc - gxzc*gyyc + tmp1 = gyyc*gzzc - gyzc**2 + tmp2 = gxzc*gyzc - gxyc*gzzc + tmp3 = gxyc*gyzc - gxzc*gyyc - idetg = one / ( gxxc*tmp1 + gxyc*tmp2 + gxzc*tmp3 ) - - g3xx = tmp1 * idetg - g3xy = tmp2 * idetg - g3xz = tmp3 * idetg - g3yy = ( gxxc*gzzc - gxzc**2 ) * idetg - g3yz = ( gxyc*gxzc - gxxc*gyzc ) * idetg - g3zz = ( gxxc*gyyc - gxyc**2 ) * idetg - - tmp1 = betax(i,j,k) * dfx(i,j,k) + & - betay(i,j,k) * dfy(i,j,k) + & - betaz(i,j,k) * dfz(i,j,k) - - tmp2 = g3xx * dfx(i,j,k)**2 + & - g3yy * dfy(i,j,k)**2 + & - g3zz * dfz(i,j,k)**2 + & - two * ( g3xy * dfx(i,j,k) * dfy(i,j,k) + & - g3xz * dfx(i,j,k) * dfz(i,j,k) + & - g3yz * dfy(i,j,k) * dfz(i,j,k) ) - if ( tmp2 .ge. zero ) then - sf(i,j,k) = tmp1 - ssign * sqrt ( alp2 * tmp2 ) + idetg = one / ( gxxc*tmp1 + gxyc*tmp2 + gxzc*tmp3 ) + + g3xx = tmp1 * idetg + g3xy = tmp2 * idetg + g3xz = tmp3 * idetg + g3yy = ( gxxc*gzzc - gxzc**2 ) * idetg + g3yz = ( gxyc*gxzc - gxxc*gyzc ) * idetg + g3zz = ( gxxc*gyyc - gxyc**2 ) * idetg + + tmp1 = betax(i,j,k) * dfx(i,j,k,l) + & + betay(i,j,k) * dfy(i,j,k,l) + & + betaz(i,j,k) * dfz(i,j,k,l) + + tmp2 = g3xx * dfx(i,j,k,l)**2 + & + g3yy * dfy(i,j,k,l)**2 + & + g3zz * dfz(i,j,k,l)**2 + & + two * ( g3xy * dfx(i,j,k,l) * dfy(i,j,k,l) + & + g3xz * dfx(i,j,k,l) * dfz(i,j,k,l) + & + g3yz * dfy(i,j,k,l) * dfz(i,j,k,l) ) + if ( tmp2 .ge. zero ) then + sf(i,j,k,l) = tmp1 - ssign * sqrt ( alp2 * tmp2 ) + else + call CCTK_WARN ( 0, "3-metric not positive definite: Stopping" ) + end if else - call CCTK_WARN ( 0, "3-metric not positive definite: Stopping" ) + sf(i,j,k,l) = zero end if - else - sf(i,j,k) = zero - end if + end do end do - end do - end do + end do + end do end if return diff --git a/src/EHFinder_mod.F90 b/src/EHFinder_mod.F90 index 6d6b319..e4d3e27 100644 --- a/src/EHFinder_mod.F90 +++ b/src/EHFinder_mod.F90 @@ -33,5 +33,6 @@ module EHFinder_mod logical :: mask_first = .true. logical :: symx, symy, symz logical :: s_symx, s_symy, s_symz - logical :: reparam_undone + logical, dimension(:), allocatable :: reparam_this_level_set, & + reparam_undone end module EHFinder_mod diff --git a/src/MoL_Init.F90 b/src/MoL_Init.F90 index 6a43ddd..4683c09 100644 --- a/src/MoL_Init.F90 +++ b/src/MoL_Init.F90 @@ -11,34 +11,37 @@ subroutine EHFinder_MoLRegister(CCTK_ARGUMENTS) CCTK_INT :: ierr, ierr_cum, varindex, rhsindex CCTK_INT :: MoLRegisterEvolved, MoLRegisterEvolvedGroup + CCTK_INT :: l + + character(len=15) :: vname DECLARE_CCTK_PARAMETERS DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS ierr_cum = 0 - call CCTK_VarIndex(varindex, 'ehfinder::f') - call CCTK_VarIndex(rhsindex, 'ehfinder::sf') -! call MoL_RegisterVar(ierr, varindex, rhsindex) -! ierr_cum = ierr_cum + ierr - ierr_cum = ierr_cum + MoLRegisterEvolved(varindex, rhsindex) + do l = 0, eh_number_level_sets - 1 + write(vname,'(a12,i1,a1)') 'ehfinder::f[', l, ']' + call CCTK_VarIndex(varindex, vname ) + write(vname,'(a13,i1,a1)') 'ehfinder::sf[', l, ']' + call CCTK_VarIndex(rhsindex, vname) + ierr_cum = ierr_cum + MoLRegisterEvolved(varindex, rhsindex) + end do if ( evolve_generators .gt. 0 ) then - call CCTK_GroupIndex (varindex, 'ehfinder::generators') - call CCTK_GroupIndex(rhsindex, 'ehfinder::generators_rhs') - -!!$ call CCTK_VarIndex(varindex, 'ehfinder::xg') -!!$ call CCTK_VarIndex(rhsindex, 'ehfinder::dxg') -!!$ ierr_cum = ierr_cum + MoLRegisterEvolved(varindex, rhsindex) -!!$ -!!$ call CCTK_VarIndex(varindex, 'ehfinder::yg') -!!$ call CCTK_VarIndex(rhsindex, 'ehfinder::dyg') -!!$ ierr_cum = ierr_cum + MoLRegisterEvolved(varindex, rhsindex) -!!$ -!!$ call CCTK_VarIndex(varindex, 'ehfinder::zg') -!!$ call CCTK_VarIndex(rhsindex, 'ehfinder::dzg') -!!$ ierr_cum = ierr_cum + MoLRegisterEvolved(varindex, rhsindex) + call CCTK_GroupIndex (varindex, 'ehfinder::xg') + call CCTK_GroupIndex(rhsindex, 'ehfinder::dxg') + + ierr_cum = ierr_cum + MoLRegisterEvolvedGroup(varindex, rhsindex) + + call CCTK_GroupIndex (varindex, 'ehfinder::yg') + call CCTK_GroupIndex(rhsindex, 'ehfinder::dyg') + + ierr_cum = ierr_cum + MoLRegisterEvolvedGroup(varindex, rhsindex) + + call CCTK_GroupIndex (varindex, 'ehfinder::zg') + call CCTK_GroupIndex(rhsindex, 'ehfinder::dzg') ierr_cum = ierr_cum + MoLRegisterEvolvedGroup(varindex, rhsindex) -- cgit v1.2.3