diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/EHFinder_Check.F90 | 140 | ||||
| -rw-r--r-- | src/EHFinder_FindSurface.F90 | 161 | ||||
| -rw-r--r-- | src/EHFinder_Generator_Sources.F90 | 299 | ||||
| -rw-r--r-- | src/EHFinder_Generator_Sources2.F90 | 217 | ||||
| -rw-r--r-- | src/EHFinder_Init.F90 | 165 | ||||
| -rw-r--r-- | src/EHFinder_Integrate2.F90 | 140 | ||||
| -rw-r--r-- | src/EHFinder_ReParametrize.F90 | 906 | ||||
| -rw-r--r-- | src/EHFinder_SetMask.F90 | 1452 | ||||
| -rw-r--r-- | src/EHFinder_SetSym.F90 | 100 | ||||
| -rw-r--r-- | src/EHFinder_Sources.F90 | 192 | ||||
| -rw-r--r-- | src/EHFinder_mod.F90 | 3 | ||||
| -rw-r--r-- | src/MoL_Init.F90 | 41 |
12 files changed, 1906 insertions, 1910 deletions
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 and find their minimum -! and maximum global cell index. - do k = kzl, kzr - do j = jyl, jyr - do i = ixl, ixr - if ( eh_mask(i,j,k) .ge. 0 ) then - if ( f(i,j,k) .lt. shell_width * delta ) then - if ( i+cctk_lbnd(1) .le. imin_loc(1) ) then - imin_loc(1) = i+cctk_lbnd(1) - end if - if ( i+cctk_lbnd(1) .ge. imax_loc(1) ) then - imax_loc(1) = i+cctk_lbnd(1) - end if - if ( j+cctk_lbnd(2) .le. imin_loc(2) ) then - imin_loc(2) = j+cctk_lbnd(2) - end if - if ( j+cctk_lbnd(2) .ge. imax_loc(2) ) then - imax_loc(2) = j+cctk_lbnd(2) + loop_over_l: do l = 1, eh_number_level_sets + +! Store the current mask and level set function + tm_mask(ixl:ixr,jyl:jyr,kzl:kzr,l) = eh_mask(ixl:ixr,jyl:jyr,kzl:kzr,l) + ftmp(ixl:ixr,jyl:jyr,kzl:kzr,l) = f(ixl:ixr,jyl:jyr,kzl:kzr,l) + + not_undone: if ( .not. reparam_undone(l) ) then + +! 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 and find their minimum +! and maximum global cell index. + do k = kzl, kzr + do j = jyl, jyr + do i = ixl, ixr + if ( eh_mask(i,j,k,l) .ge. 0 ) then + if ( f(i,j,k,l) .lt. shell_width * delta ) then + if ( i+cctk_lbnd(1) .le. imin_loc(1) ) then + imin_loc(1) = i+cctk_lbnd(1) + end if + if ( i+cctk_lbnd(1) .ge. imax_loc(1) ) then + imax_loc(1) = i+cctk_lbnd(1) + end if + if ( j+cctk_lbnd(2) .le. imin_loc(2) ) then + imin_loc(2) = j+cctk_lbnd(2) + end if + if ( j+cctk_lbnd(2) .ge. imax_loc(2) ) then + imax_loc(2) = j+cctk_lbnd(2) + end if + if ( k+cctk_lbnd(3) .le. imin_loc(3) ) then + imin_loc(3) = k+cctk_lbnd(3) + end if + if ( k+cctk_lbnd(3) .ge. imax_loc(3) ) then + imax_loc(3) = k+cctk_lbnd(3) + end if + end if end if - if ( k+cctk_lbnd(3) .le. imin_loc(3) ) then - imin_loc(3) = k+cctk_lbnd(3) - end if - if ( k+cctk_lbnd(3) .ge. imax_loc(3) ) then - imax_loc(3) = k+cctk_lbnd(3) - end if - end if - end if + end do + end do end do - end do - end do -! Reduce over all processors to get the global indeces of the -! rectangular box containing all cells with f<shell_width*delta - call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, min_handle, & +! Reduce over all processors to get the global indeces of the +! rectangular box containing all cells with f<shell_width*delta + call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, min_handle, & imin_loc, imin_glob, 3, CCTK_VARIABLE_INT ) - if ( ierr .ne. zero ) then - call CCTK_WARN ( 0, 'Reduction of array imin_loc failed' ) - end if + if ( ierr .ne. zero ) then + call CCTK_WARN ( 0, 'Reduction of array imin_loc failed' ) + end if - call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, max_handle, & + call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, max_handle, & imax_loc, imax_glob, 3, CCTK_VARIABLE_INT ) - if ( ierr .ne. zero ) then - call CCTK_WARN ( 0, 'Reduction of array imax_loc failed' ) - end if + if ( ierr .ne. zero ) then + call CCTK_WARN ( 0, 'Reduction of array imax_loc failed' ) + end if -! Convert into local indeces. Note that these might be less than 1 or -! larger than the local grid size if the box does not overlap with the -! local grid. - i1 = imin_glob(1)-cctk_lbnd(1) - i2 = imax_glob(1)-cctk_lbnd(1) - j1 = imin_glob(2)-cctk_lbnd(2) - j2 = imax_glob(2)-cctk_lbnd(2) - k1 = imin_glob(3)-cctk_lbnd(3) - k2 = imax_glob(3)-cctk_lbnd(3) - -! Find the minimum value of f on the various faces of the rectangular box -! if part of the face is present on the current grid. - if ( ( 1 .le. i1 ) .and. ( i1 .le. nx ) ) then - fimin_loc(1) = minval ( f(i1,jyl:jyr,kzl:kzr) ) - end if - if ( ( 1 .le. i2 ) .and. ( i2 .le. nx ) ) then - fimax_loc(1) = minval ( f(i2,jyl:jyr,kzl:kzr) ) - end if - if ( ( 1 .le. j1 ) .and. ( j1 .le. ny ) ) then - fimin_loc(2) = minval ( f(ixl:ixr,j1,kzl:kzr) ) - end if - if ( ( 1 .le. j2 ) .and. ( j2 .le. ny ) ) then - fimax_loc(2) = minval ( f(ixl:ixr,j2,kzl:kzr) ) - end if - if ( ( 1 .le. k1 ) .and. ( k1 .le. nz ) ) then - fimin_loc(3) = minval ( f(ixl:ixr,jyl:jyr,k1) ) - end if - if ( ( 1 .le. k2 ) .and. ( k2 .le. nz ) ) then - fimax_loc(3) = minval ( f(ixl:ixr,jyl:jyr,k2) ) - end if +! Convert into local indeces. Note that these might be less than 1 or +! larger than the local grid size if the box does not overlap with the +! local grid. + i1 = imin_glob(1)-cctk_lbnd(1) + i2 = imax_glob(1)-cctk_lbnd(1) + j1 = imin_glob(2)-cctk_lbnd(2) + j2 = imax_glob(2)-cctk_lbnd(2) + k1 = imin_glob(3)-cctk_lbnd(3) + k2 = imax_glob(3)-cctk_lbnd(3) + +! Find the minimum value of f on the various faces of the rectangular +! box if part of the face is present on the current grid. + if ( ( 1 .le. i1 ) .and. ( i1 .le. nx ) ) then + fimin_loc(1) = minval ( f(i1,jyl:jyr,kzl:kzr,l) ) + end if + if ( ( 1 .le. i2 ) .and. ( i2 .le. nx ) ) then + fimax_loc(1) = minval ( f(i2,jyl:jyr,kzl:kzr,l) ) + end if + if ( ( 1 .le. j1 ) .and. ( j1 .le. ny ) ) then + fimin_loc(2) = minval ( f(ixl:ixr,j1,kzl:kzr,l) ) + end if + if ( ( 1 .le. j2 ) .and. ( j2 .le. ny ) ) then + fimax_loc(2) = minval ( f(ixl:ixr,j2,kzl:kzr,l) ) + end if + if ( ( 1 .le. k1 ) .and. ( k1 .le. nz ) ) then + fimin_loc(3) = minval ( f(ixl:ixr,jyl:jyr,k1,l) ) + end if + if ( ( 1 .le. k2 ) .and. ( k2 .le. nz ) ) then + fimax_loc(3) = minval ( f(ixl:ixr,jyl:jyr,k2,l) ) + end if -! Reduce to get the minimum values of f on the faces of the box - call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, min_handle, & +! Reduce to get the minimum values of f on the faces of the box + call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, min_handle, & fimin_loc, fimin_glob, 3, CCTK_VARIABLE_REAL ) - if ( ierr .ne. zero ) then - call CCTK_WARN ( 0, 'Reduction of array fimin_loc failed' ) - end if + if ( ierr .ne. zero ) then + call CCTK_WARN ( 0, 'Reduction of array fimin_loc failed' ) + end if - call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, min_handle, & + call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, min_handle, & fimax_loc, fimax_glob, 3, CCTK_VARIABLE_REAL ) - if ( ierr .ne. zero ) then - call CCTK_WARN ( 0, 'Reduction of array fimax_loc failed' ) - end if + if ( ierr .ne. zero ) then + call CCTK_WARN ( 0, 'Reduction of array fimax_loc failed' ) + end if - print*, imin_glob(1), imax_glob(1), fimin_glob(1), fimax_glob(1) - print*, imin_glob(2), imax_glob(2), fimin_glob(2), fimax_glob(2) - print*, imin_glob(3), imax_glob(3), fimin_glob(3), fimax_glob(3) - end if +! print*, imin_glob(1), imax_glob(1), fimin_glob(1), fimax_glob(1) +! print*, imin_glob(2), imax_glob(2), fimin_glob(2), fimax_glob(2) +! print*, imin_glob(3), imax_glob(3), fimin_glob(3), fimax_glob(3) + end if -! Now check and see if any interior excised points need to be activated. - 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) .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) |