aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authordiener <diener@2a26948c-0e4f-0410-aee8-f1d3e353619c>2003-08-13 17:11:22 +0000
committerdiener <diener@2a26948c-0e4f-0410-aee8-f1d3e353619c>2003-08-13 17:11:22 +0000
commitade581e1acd584ea6e54ff7a2aad11bc23477f9d (patch)
treedf39165329f5ff9aeb33d4d24bb9194b9dec9cde /src
parentd032c16e368d85d834388d81ab5171df134ad238 (diff)
Major modification of almost all source files to use vector groups for various
variables. Not tested in all detail, but the standard features seem to work. The version before all these changes was tagged with PRE_MULTI. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/EHFinder/trunk@126 2a26948c-0e4f-0410-aee8-f1d3e353619c
Diffstat (limited to 'src')
-rw-r--r--src/EHFinder_Check.F90140
-rw-r--r--src/EHFinder_FindSurface.F90161
-rw-r--r--src/EHFinder_Generator_Sources.F90299
-rw-r--r--src/EHFinder_Generator_Sources2.F90217
-rw-r--r--src/EHFinder_Init.F90165
-rw-r--r--src/EHFinder_Integrate2.F90140
-rw-r--r--src/EHFinder_ReParametrize.F90906
-rw-r--r--src/EHFinder_SetMask.F901452
-rw-r--r--src/EHFinder_SetSym.F90100
-rw-r--r--src/EHFinder_Sources.F90192
-rw-r--r--src/EHFinder_mod.F903
-rw-r--r--src/MoL_Init.F9041
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)
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-19 19:11:35 +0000