From be12e621fe75e9b7fabfe531c6f03dabf98882ad Mon Sep 17 00:00:00 2001 From: diener Date: Tue, 30 Sep 2003 09:03:44 +0000 Subject: Preparing for release. Added comments, deleted obsolete pieces of code, fixed a few bugs. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/EHFinder/trunk@143 2a26948c-0e4f-0410-aee8-f1d3e353619c --- src/EHFinder_Check.F90 | 21 +- src/EHFinder_Constants.F90 | 3 - src/EHFinder_FindSurface.F90 | 259 ++++++------- src/EHFinder_Generator_Sources.F90 | 118 +++--- src/EHFinder_Generator_Sources2.F90 | 113 +++--- src/EHFinder_Init.F90 | 93 ++--- src/EHFinder_Integrate.F90 | 704 ++++++++++++++++++++++++++++++++++++ src/EHFinder_ParamCheck.F90 | 29 +- src/EHFinder_ReadData.F90 | 138 +++---- src/EHFinder_SetMask.F90 | 86 ++--- src/EHFinder_SetSym.F90 | 151 +++----- src/EHFinder_Sources.F90 | 197 +++++----- src/EHFinder_mod.F90 | 20 +- src/MoL_Init.F90 | 47 +-- 14 files changed, 1269 insertions(+), 710 deletions(-) create mode 100644 src/EHFinder_Integrate.F90 diff --git a/src/EHFinder_Check.F90 b/src/EHFinder_Check.F90 index 22b3703..1fd07d4 100644 --- a/src/EHFinder_Check.F90 +++ b/src/EHFinder_Check.F90 @@ -17,7 +17,7 @@ subroutine EHFinder_ReInitialize_Check(CCTK_ARGUMENTS) DECLARE_CCTK_FUNCTIONS CCTK_INT :: i, j, k, l - CCTK_REAL :: fmin_extremum, fmin_extremum_loc, fmax_loc, fmin_bound_loc + CCTK_REAL :: fmin_extremum, fmin_extremum_loc, fmax_loc character(len=128) :: info_message logical :: check_re_init CCTK_INT, dimension(3) :: ex_loc @@ -45,15 +45,20 @@ subroutine EHFinder_ReInitialize_Check(CCTK_ARGUMENTS) ! values in the extremum. fmax_loc = zero +! On the local processor find the maximum value of f and the maximum value +! of f in any local extrema where f is negative. 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 + 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 @@ -61,6 +66,7 @@ subroutine EHFinder_ReInitialize_Check(CCTK_ARGUMENTS) end do end do +! Maximum reduction over all processors. call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, max_handle, & fmin_extremum_loc, fmin_extremum, CCTK_VARIABLE_REAL ) @@ -68,15 +74,18 @@ subroutine EHFinder_ReInitialize_Check(CCTK_ARGUMENTS) call CCTK_WARN(0,'Reduction of fmax_extremum failed') end if +! Write an info message about the local extremum. 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) +! If the local extremum indicate and imminent pinch off undo the +! re-initialization if this is requested by the user. if ( re_init_undo .gt. 0 ) then if ( fmin_extremum .gt. -two * minval(cctk_delta_space) ) then - write(info_message,'(a45,i3,a6,f12.5)') & + write(info_message,'(a45,i3,a33)') & 'Re-initialization of level set ', l, & - 'undone due to imminent pinch-off' + ' undone due to imminent pinch-off' call CCTK_INFO(info_message) f(:,:,:,l) = fbak(:,:,:,l) eh_mask(:,:,:,l) = eh_mask_bak(:,:,:,l) diff --git a/src/EHFinder_Constants.F90 b/src/EHFinder_Constants.F90 index fdd5b6f..2162631 100644 --- a/src/EHFinder_Constants.F90 +++ b/src/EHFinder_Constants.F90 @@ -10,13 +10,10 @@ module EHFinder_Constants CCTK_REAL, parameter :: three = 3.0d0 CCTK_REAL, parameter :: four = 4.0d0 CCTK_REAL, parameter :: eight = 8.0d0 - CCTK_REAL, parameter :: ten = 10.0d0 CCTK_REAL, parameter :: half = 0.5d0 CCTK_REAL, parameter :: onethird = one / three CCTK_REAL, parameter :: twothirds = two / three CCTK_REAL, parameter :: fourthirds = four / three CCTK_REAL, parameter :: quarter = 0.25d0 - CCTK_REAL, parameter :: tenth = 0.1d0 - CCTK_REAL, parameter :: huge = 1d23 CCTK_REAL, parameter :: pi = 3.1415926535897932385d0 end module EHFinder_Constants diff --git a/src/EHFinder_FindSurface.F90 b/src/EHFinder_FindSurface.F90 index 2537408..d53db6e 100644 --- a/src/EHFinder_FindSurface.F90 +++ b/src/EHFinder_FindSurface.F90 @@ -21,13 +21,14 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) CCTK_INT :: interp_handle2 character(len=80) :: info_message + character(len=128) :: warn_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 + CCTK_INT, dimension(2) :: gsh, lsh, lbnd, ubnd, nghost CCTK_REAL, dimension(3) :: center_loc CCTK_INT :: nc_loc, nc @@ -55,52 +56,52 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) ! set to -1. find_surface_status = 0 sn = surface_counter - integrate_counter + 1 + +! Set the index to the current level set function. l = levelset_counter + +! Write an info message about what is going on. 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 - ! Grid arrays. - call CCTK_GroupbboxGN ( status, cctkGH, 4, bbox, "ehfinder::surface_arrays" ) +! Find information about the local and global properties of the 2D +! Grid arrays. + call CCTK_GroupbboxGN ( status, cctkGH, 4, bbox, 'ehfinder::surface_arrays' ) if ( status .lt. 0 ) then - call CCTK_WARN ( 0, "cannot get bounding box for surface arrays" ) + call CCTK_WARN ( 0, 'cannot get bounding box for surface arrays' ) end if -! print*,bbox - call CCTK_GroupgshGN ( status, cctkGH, 2, gsh, "ehfinder::surface_arrays" ) + call CCTK_GroupgshGN ( status, cctkGH, 2, gsh, 'ehfinder::surface_arrays' ) if ( status .lt. 0 ) then - call CCTK_WARN ( 0, "cannot get global size for surface arrays" ) + call CCTK_WARN ( 0, 'cannot get global size for surface arrays' ) end if -! print*,gsh - call CCTK_GrouplbndGN ( status, cctkGH, 2, lbnd, "ehfinder::surface_arrays" ) + call CCTK_GrouplbndGN ( status, cctkGH, 2, lbnd, 'ehfinder::surface_arrays' ) if ( status .lt. 0 ) then - call CCTK_WARN ( 0, "cannot get lower bounds for surface arrays" ) + call CCTK_WARN ( 0, 'cannot get lower bounds for surface arrays' ) end if -! print*,lbnd - call CCTK_GroupubndGN ( status, cctkGH, 2, ubnd, "ehfinder::surface_arrays" ) + call CCTK_GroupubndGN ( status, cctkGH, 2, ubnd, 'ehfinder::surface_arrays' ) if ( status .lt. 0 ) then - call CCTK_WARN ( 0, "cannot get upper bounds for surface arrays" ) + call CCTK_WARN ( 0, 'cannot get upper bounds for surface arrays' ) end if -! print*,ubnd - call CCTK_GrouplshGN ( status, cctkGH, 2, lsh, "ehfinder::surface_arrays" ) + call CCTK_GrouplshGN ( status, cctkGH, 2, lsh, 'ehfinder::surface_arrays' ) if ( status .lt. 0 ) then - call CCTK_WARN ( 0, "cannot get local size for surface arrays" ) + call CCTK_WARN ( 0, 'cannot get local size for surface arrays' ) end if -! print*,lsh - call CCTK_GroupnghostzonesGN ( status, cctkGH, 2, nghost, "ehfinder::surface_arrays" ) + call CCTK_GroupnghostzonesGN ( status, cctkGH, 2, nghost, & + 'ehfinder::surface_arrays' ) if ( status .lt. 0 ) then - call CCTK_WARN ( 0, "cannot get local size for surface arrays" ) + call CCTK_WARN ( 0, 'cannot get local size for surface arrays' ) end if - ! If the user defined center point should be used... +! If the user defined center point should be used... if ( use_user_center .gt. 0 ) then center(1) = center_x center(2) = center_y center(3) = center_z else - ! IF not find an approximation for it by averaging the coodinates that - ! are within one grid cell distance from the surface in question. +! If not find an approximation for it by averaging the coodinates that +! are within one grid cell distance from the surface in question. nc_loc = 0; center_loc = zero min_delta = minval ( cctk_delta_space ) do k = kzl, kzr @@ -130,7 +131,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) end do end do - ! Reduce over the processors. +! Reduce over the processors. call CCTK_ReduceLocScalar ( status, cctkGH, -1, sum_handle, & nc_loc, nc, CCTK_VARIABLE_INT ) call CCTK_ReduceLocArrayToArray1D ( status, cctkGH, -1, sum_handle, & @@ -138,14 +139,10 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) CCTK_VARIABLE_REAL ) ncinv = one / nc center = ncinv * center -! print*,'EHFinder_FindSurfaces : ', sn, center -! print* end if -! print*,'Center : ',center -! print*,'use_user_center : ',use_user_center - ! Try to figure out the symmetries. First find the range in coordinates - ! for the points inside the current surface. +! Try to figure out the symmetries. First find the range in coordinates +! for the points inside the current surface. crange_min_loc(1) = minval ( x, mask = sc .eq. sn ) crange_min_loc(2) = minval ( y, mask = sc .eq. sn ) crange_min_loc(3) = minval ( z, mask = sc .eq. sn ) @@ -153,18 +150,15 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) crange_max_loc(2) = maxval ( y, mask = sc .eq. sn ) crange_max_loc(3) = maxval ( z, mask = sc .eq. sn ) +! Reduce over all processors. call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, min_handle, & crange_min_loc, crange_min_glob, & 3, CCTK_VARIABLE_REAL ) call CCTK_ReduceLocArrayToArray1D ( ierr, cctkGH, -1, max_handle, & crange_max_loc, crange_max_glob, & 3, CCTK_VARIABLE_REAL ) -! print*,crange_min_loc,crange_max_loc -! print*,crange_min_glob,crange_max_glob -! print* -! print*,'Symmetries : ', symx, symy, symz - ! Full mode. Modify if necessary. +! Full mode. Modify if necessary. ltheta = zero; utheta = pi lphi = zero; uphi = two * pi sym_factor = one @@ -172,9 +166,12 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) if ( CCTK_EQUALS ( domain, 'bitant' ) ) then - ! Bitant with symmetry in the z-direction +! Bitant with symmetry in the z-direction if ( CCTK_EQUALS ( bitant_plane, 'xy' ) ) then +! If the minimum z-coordinate of the points inside the surface is less +! than zero, then the surface must share the grid symmetry across the +! xy-plane. if ( crange_min_glob(3) .lt. zero ) then ltheta = zero; utheta = half * pi lphi = zero; uphi = two * pi @@ -184,7 +181,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) end if end if - ! Bitant with symmetry in the y-direction +! Bitant with symmetry in the y-direction if ( CCTK_EQUALS ( bitant_plane, 'xz' ) ) then if ( crange_min_glob(2) .lt. zero ) then @@ -196,7 +193,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) end if end if - ! Bitant with symmetry in the x-direction +! Bitant with symmetry in the x-direction if ( CCTK_EQUALS ( bitant_plane, 'yz' ) ) then if ( crange_min_glob(1) .lt. zero ) then @@ -211,7 +208,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) if ( CCTK_EQUALS ( domain, 'quadrant' ) ) then - ! Quadrant in x-direction +! Quadrant in x-direction if ( CCTK_EQUALS ( quadrant_direction, 'x' ) ) then if ( ( crange_min_glob(3) .lt. zero ) .and. & @@ -236,7 +233,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) end if end if - ! Quadrant in y-direction +! Quadrant in y-direction if ( CCTK_EQUALS ( quadrant_direction, 'y' ) ) then if ( ( crange_min_glob(3) .lt. zero ) .and. & @@ -261,7 +258,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) end if end if - ! Quadrant in z-direction +! Quadrant in z-direction if ( CCTK_EQUALS ( quadrant_direction,'z' ) ) then if ( ( crange_min_glob(2) .lt. zero ) .and. & @@ -287,7 +284,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) end if end if - ! Octant +! Octant if ( CCTK_EQUALS ( domain, 'octant' ) ) then if ( ( crange_min_glob(3) .lt. zero ) .and. & @@ -340,24 +337,30 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) end if end if +! Set the theta index to be used for the equatorial circumference integration +! based on the symmetries. if ( .not. s_symz ) then githeta = ( ntheta - 1 ) / 2 + 1 else githeta = ntheta end if +! Set the phi index to be used for the polar circumference integration based +! on the symmetries. if ( ( .not. s_symy ) .and. ( s_symx ) ) then gjphi = ( nphi - 1 ) / 2 + 1 else gjphi = 1 end if +! Set up the symmetry muliplication factors based on the range in the +! angular coordinates. theta_sym_factor = two * pi / ( utheta - ltheta ) phi_sym_factor = two * pi / ( uphi - lphi ) - ! 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. +! 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. dtheta = ( utheta - ltheta ) / ( ntheta - 1 ) dphi = ( uphi - lphi ) / ( nphi - 1 ) @@ -370,64 +373,71 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) cphi(:,j) = lphi + dphi * ( j +lbnd(2) - 1 ) end do -! print*,'Delta : ',dtheta,dphi - ! Calculate sines and cosines and store them in grid arrays since they - ! are expensive. +! Calculate sines and cosines and store them in grid arrays since they +! are computationally expensive and memory for 2D grid arrays is cheap. sintheta = sin(ctheta) costheta = cos(ctheta) sinphi = sin(cphi) cosphi = cos(cphi) - ! 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 - ! and is initialised to 2*Delta. The grid array n_since_last_reduction is - ! a help in keeping track on how long ago drsurf for a given point has been - ! reduced and is initialised to -1 to indicate that no reduction has taken - ! place. +! 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 +! and is initialised to 2*Delta. The grid array n_since_last_reduction is +! a help in keeping track on how long ago drsurf for a given point has been +! reduced and is initialised to -1 to indicate that no reduction has taken +! place. drsurf = two * min_delta rsurf = zero n_since_last_reduction = -1 - ! Get and interpolation handle for required polynomial interpolation. - ! If Lagrange polynomial interpolation does not work try Hermite - ! polynomial interpolation, which avoids problems with non-continuous - ! derivatives in the Newton iteration later on. +! Get and interpolation handle for required polynomial interpolation. +! If Lagrange polynomial interpolation does not work try Hermite +! polynomial interpolation, which avoids problems with non-continuous +! derivatives in the Newton iteration later on. - ! First convert the string parameter to a Fortran string. +! First convert the string parameter to a Fortran string. call CCTK_FortranString ( surface_interp_len, surface_interpolator, & surface_interp ) - ! Then get the hande ... +! Then get the handle ... call CCTK_InterpHandle ( interp_handle, surface_interp(1:surface_interp_len) ) +! The following is done because some C-preprocessors seem to have problems +! handling character strings that are continued on several lines. if ( interp_handle .lt. 0 ) then - call CCTK_WARN( 0, "Cannot get handle for interpolation. Forgot to activate an implementation providing interpolation operators??" ) + warn_message = 'Cannot get handle for interpolation. ' + warn_message = warn_message//'Forgot to activate an implementation ' + warn_message = warn_message//'providing interpolation operators??' + call CCTK_WARN( 0, warn_message ) end if - ! Write the required interpolation order parameter into a fortran string. +! Write the required interpolation order parameter into a fortran string. write(surface_order,'(a6,i1)') 'order=',surface_interpolation_order - ! Then create a table from the string. +! Then create a table from the string. call Util_TableCreateFromString ( table_handle, surface_order ) if ( table_handle .lt. 0 ) then - call CCTK_WARN( 0, "Cannot create parameter table for interpolator" ) + call CCTK_WARN( 0, 'Cannot create parameter table for interpolator' ) end if - ! Get the 3D coordinate system handle. - call CCTK_CoordSystemHandle ( coord_system_handle, "cart3d" ) +! Get the 3D coordinate system handle. + call CCTK_CoordSystemHandle ( coord_system_handle, 'cart3d' ) if ( coord_system_handle .lt. 0) then - call CCTK_WARN( 0, "Cannot get handle for cart3d coordinate system. Forgot to activate an implementation providing coordinates ??" ) + warn_message = 'Cannot get handle for cart3d coordinate system. ' + warn_message = warn_message//'Forgot to activate an implementation ' + warn_message = warn_message//'providing coordinates??' + call CCTK_WARN( 0, warn_message ) endif - ! Get pointers to the grid arrays containing the interpolation points. +! Get pointers to the grid arrays containing the interpolation points. interp_coords(1) = CCTK_PointerTo(interp_x) interp_coords(2) = CCTK_PointerTo(interp_y) interp_coords(3) = CCTK_PointerTo(interp_z) - ! Get pointer to the interpolation output array +! Get pointer to the interpolation output array out_array(1) = CCTK_PointerTo(f_interp) - ! Get the variable index for f +! Get the variable index for f write(vname,'(a12,i1,a1)') 'ehfinder::f[', l-1, ']' call CCTK_VarIndex ( in_array, vname ) @@ -442,22 +452,22 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) do j = 1, lsh(2) do i = 1, lsh(1) - ! The radius that are tested in this direction. +! The radius that are tested in this direction. rloc = rsurf(i,j) + drsurf(i,j) - ! Calculate the interpolation points. +! Calculate the interpolation points. interp_x(i,j) = center(1) + rloc * sintheta(i,j) * cosphi(i,j) interp_y(i,j) = center(2) + rloc * sintheta(i,j) * sinphi(i,j) interp_z(i,j) = center(3) + rloc * costheta(i,j) - ! If the point has been reduced once, increment n_since_last_reduction. - ! +! If the point has been reduced once, increment n_since_last_reduction. if ( n_since_last_reduction(i,j) .ge. 0 ) then n_since_last_reduction(i,j) = n_since_last_reduction(i,j) + 1 end if end do end do +! Call the interpolator. call CCTK_InterpGridArrays ( status, cctkGH, 3, interp_handle, & table_handle, coord_system_handle, & lsh(1) * lsh(2), CCTK_VARIABLE_REAL, & @@ -470,18 +480,22 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) return end if +! Find the maximal change in rsurf. maxdr_loc = maxval ( drsurf ) call CCTK_ReduceLocScalar ( status, cctkGH, -1, max_handle, & maxdr_loc, maxdr, CCTK_VARIABLE_REAL ) +! Find the maximal value of interpolated f values. maxf_loc = maxval ( abs ( f_interp ) ) call CCTK_ReduceLocScalar ( status, cctkGH, -1, max_handle, & maxf_loc, maxf, CCTK_VARIABLE_REAL ) - maxl = maxloc ( abs ( f_interp ) ) - +! If the maximal change in rsurf is small enough, we are satisfied that +! the Newton iteration will be able to find the root. if ( maxdr .lt. min_delta ) exit find_starting_points +! Check if we have crossed the zero point and if we have reduce the +! change in rsurf we are going to try next. do j = 1, lsh(2) do i = 1, lsh(1) if ( f_interp(i,j) .gt. zero ) then @@ -497,28 +511,28 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) end do end do find_starting_points -! print*,'Max values: ',maxval(rsurf),maxval(f_interp) -! print*,'Min values: ',minval(rsurf),minval(f_interp) -! print*,'found starting points' + +! Now for the Newton iteration we also need the interpolates derivatives. out_array(2) = CCTK_PointerTo(dfdx_interp) out_array(3) = CCTK_PointerTo(dfdy_interp) out_array(4) = CCTK_PointerTo(dfdz_interp) call Util_TableSetIntArray ( status, table_handle, 4, & - op_indices, "operand_indices" ) + op_indices, 'operand_indices' ) if ( status .lt. 0 ) then - call CCTK_WARN ( 0, "Cannot set operand indices array in parameter table" ) + call CCTK_WARN ( 0, 'Cannot set operand indices array in parameter table' ) endif call Util_TableSetIntArray ( status, table_handle, 4, & - op_codes, "operation_codes" ) + op_codes, 'operation_codes' ) if ( status .lt. 0 ) then - call CCTK_WARN ( 0, "Cannot set operation codes array in parameter table" ) + call CCTK_WARN ( 0, 'Cannot set operation codes array in parameter table' ) endif ! find the surface points in all directions simultaneously. find_surface_points: do k = 1, 15 +! Calculate the interpolation points. 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) @@ -527,6 +541,7 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) end do end do +! Call the interpolator. call CCTK_InterpGridArrays ( status, cctkGH, 3, interp_handle, & table_handle, coord_system_handle, & lsh(1) * lsh(2), CCTK_VARIABLE_REAL, & @@ -539,12 +554,13 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) return end if +! Find the maximum value of interpolated f values. maxf_loc = maxval ( abs ( f_interp ) ) call CCTK_ReduceLocScalar ( status, cctkGH, -1, max_handle, & maxf_loc, maxf, CCTK_VARIABLE_REAL ) - maxl = maxloc ( abs ( f_interp ) ) - +! If the maximum interpolated value of f is to big, perform a Newton +! iteration if ( maxf .gt. eps ) then do j = 1, lsh(2) do i = 1, lsh(1) @@ -555,16 +571,22 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) end do end do else +! Else we are done exit find_surface_points end if -! call CCTK_SyncGroup ( status, cctkGH, "ehfinder::surface_arrays" ) end do find_surface_points - if ( k .gt. 15 ) call CCTK_WARN(1,'Newton iteration did not converge! Area may be inaccurate') -! print*,'found_surface_points' + +! If we used to many iterations, the Newton iteration did not converge, so +! print a warning. + if ( k .gt. 15 ) then + warn_message = 'Newton iteration did not converge! Area may be inaccurate' + call CCTK_WARN(1, warn_message ) + end if drdtheta = zero drdphi = zero + ! Calculate the derivatives of r. ! First for interior points @@ -661,13 +683,11 @@ subroutine EHFinder_FindSurface(CCTK_ARGUMENTS) rsurf(im,jm-2) ) end if - call CCTK_SyncGroup ( status, cctkGH, "ehfinder::surface_arrays" ) - +! I do not remember why this is done. I will have to investigate. integrate_counter = integrate_counter - 1 -! print*,'Calculated derivatives' -! print*,rsurf + +! Get rid of the interpolation table. call Util_TableDestroy ( status, table_handle ) -! print*,'gjphi = ',gjphi end subroutine EHFinder_FindSurface @@ -718,7 +738,6 @@ subroutine EHFinder_CountSurfacesInit(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS -! print*,'EHFinder_CountSurfacesInit called' more_surfaces = 1 surface_counter = 0 sc = 0 @@ -744,27 +763,27 @@ subroutine EHFinder_CountSurfaces(CCTK_ARGUMENTS) ! the 3D grid functions. #include "include/physical_part.h" - ! Find the location of the minimum among active points not already found - ! to be inside a surface. +! 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,levelset_counter), & ( sc(ixl:ixr,jyl:jyr,kzl:kzr) .eq. 0 ) .and. & ( eh_mask(ixl:ixr,jyl:jyr,kzl:kzr,levelset_counter) .ge. 0 ) ) - ! Find the value of f at that location. +! Find the value of f at that location. minf_loc = f(minl(1)+ixl-1,minl(2)+jyl-1,minl(3)+kzl-1,levelset_counter) - ! Find the minimum over all processors. +! Find the minimum over all processors. call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, min_handle, & minf_loc, minf_glob, CCTK_VARIABLE_REAL ) if ( ierr .ne. 0 ) then call CCTK_WARN(0,'Reduction of minf_glob failed') end if - ! To figure out on which processor the minimum is located I set procpos_loc - ! to the local processor number if the local minimum is equal to the global - ! minimum otherwise I set it to the total number of processors. When I then - ! find the minimum of procpos_loc over all processors I get the processor - ! containing the minimum. +! To figure out on which processor the minimum is located I set procpos_loc +! to the local processor number if the local minimum is equal to the global +! minimum otherwise I set it to the total number of processors. When I then +! find the minimum of procpos_loc over all processors I get the processor +! containing the minimum. if ( minf_loc .eq. minf_glob ) then procpos_loc = CCTK_MyProc ( cctkGH ) else @@ -777,8 +796,8 @@ subroutine EHFinder_CountSurfaces(CCTK_ARGUMENTS) call CCTK_WARN(0,'Reduction of minf_glob failed') end if - ! If minf_glob<0 there must exist another surface so lets set up the variables - ! to start marking and counting these points. +! If minf_glob<0 there must exist another surface so lets set up the variables +! to start marking and counting these points. if ( minf_glob .lt. zero ) then surface_counter = surface_counter + 1 points_counter = 0 @@ -813,13 +832,13 @@ subroutine EHFinder_MarkSurfaces(CCTK_ARGUMENTS) ! the 3D grid functions. #include "include/physical_part.h" - ! Initialise the local point counter. +! Initialise the local point counter. n_loc = 0 - ! Check if any points have a marked point as a neighbour. Only look at - ! active points with f<0 that are not already marked. Take care at the - ! boundaries (physical and excision) not to look at inactive or not - ! exisiting points. +! Check if any points have a marked point as a neighbour. Only look at +! active points with f<0 that are not already marked. Take care at the +! boundaries (physical and excision) not to look at inactive or not +! exisiting points. do k = kzl, kzr do j = jyl, jyr do i = ixl, ixr @@ -828,7 +847,7 @@ subroutine EHFinder_MarkSurfaces(CCTK_ARGUMENTS) ( 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 +! Find out if we are on the boundary if ( i + cctk_lbnd(1) .eq. 1 ) then i1 = 0 else @@ -881,31 +900,25 @@ subroutine EHFinder_MarkSurfaces(CCTK_ARGUMENTS) end do end do - ! Reduce the number of newly marked points over all processors +! Reduce the number of newly marked points over all processors call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, sum_handle, & n_loc, n_glob, CCTK_VARIABLE_INT ) if ( ierr .ne. 0 ) then call CCTK_WARN(0,'Reduction of n_glob failed') end if - ! If the total number of newly marked points is 0, set more_points to 0 - ! to indicate that we are finished. Otherwise continue and add the number - ! of points to the points_counter. +! If the total number of newly marked points is 0, set more_points to 0 +! to indicate that we are finished. Otherwise continue and add the number +! of points to the points_counter. if ( n_glob .eq. 0 ) then more_points = 0 end if points_counter = points_counter + n_glob - ! This should not be necessary, since I have a SYNC in the schedule.ccl, - ! but for some reason that does not work. - call CCTK_SyncGroup ( ierr, cctkGH, "ehfinder::surface_index" ) - - call CartSymVN(ierr,cctkGH,'ehfinder::sc') - end subroutine EHFinder_MarkSurfaces -!This routine prints out information about the found surfaces. +! This routine prints out information about the found surfaces. subroutine EHFinder_InfoSurfaces(CCTK_ARGUMENTS) use EHFinder_mod @@ -918,7 +931,7 @@ subroutine EHFinder_InfoSurfaces(CCTK_ARGUMENTS) character(len=80) :: info_message - ! Write out how many surfaces were found. +! Write out how many surfaces were found. if ( surface_counter .eq. 1 ) then write(info_message,'(a9,i3,a22,i3)') 'There is ',surface_counter, & ' surface in level set ',levelset_counter diff --git a/src/EHFinder_Generator_Sources.F90 b/src/EHFinder_Generator_Sources.F90 index bb8e350..28ce33a 100644 --- a/src/EHFinder_Generator_Sources.F90 +++ b/src/EHFinder_Generator_Sources.F90 @@ -19,6 +19,7 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS) CCTK_INT :: interp_handle, table_handle, status, coord_system_handle character(len=200) :: gen_interp + character(len=128) :: warn_message CCTK_INT :: gen_interp_len character(len=7) :: gen_order character(len=15) :: vname @@ -37,46 +38,53 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS) CCTK_REAL :: alp2, psi4, dfux, dfuy, dfuz, factor, ssign CCTK_REAL :: idetg, guxx, guxy, guxz, guyy, guyz, guzz +! Set the sign correctly depending on the surface direction. if ( CCTK_EQUALS ( surface_direction, 'outward' ) ) ssign = one if ( CCTK_EQUALS ( surface_direction, 'inward' ) ) ssign = -one out_types = CCTK_VARIABLE_REAL - ! Convert the generator_interpolator string parameter to a Fortran string. +! Convert the generator_interpolator string parameter to a Fortran string. call CCTK_FortranString ( gen_interp_len, generator_interpolator, & gen_interp ) - ! Get the corresponding interpolator handle. +! Get the corresponding interpolator handle. call CCTK_InterpHandle ( interp_handle, gen_interp ) if ( interp_handle .lt. 0 ) then - call CCTK_WARN( 0, "Cannot get handle for interpolation. Forgot to activate an implementation providing interpolation operators??" ) + warn_message = 'Cannot get handle for interpolation. ' + warn_message = warn_message//'Forgot to activate an implementation ' + warn_message = warn_message//'providing interpolation operators??' + call CCTK_WARN( 0, warn_message ) end if - ! Convert the interpolation order parameter to a Fortran string to be placed - ! in the interpolator table. Note that the order is assumed to contain only - ! 1 digit. +! Convert the interpolation order parameter to a Fortran string to be placed +! in the interpolator table. Note that the order is assumed to contain only +! 1 digit. write(gen_order,'(a6,i1)') 'order=',generator_interpolation_order - ! Create the table directly from the string. +! Create the table directly from the string. call Util_TableCreateFromString ( table_handle, gen_order ) if ( table_handle .lt. 0 ) then - call CCTK_WARN( 0, "Cannot create parameter table for interpolator" ) + call CCTK_WARN( 0, 'Cannot create parameter table for interpolator' ) end if - ! Get the 3D coordinate system handle. - call CCTK_CoordSystemHandle ( coord_system_handle, "cart3d" ) +! Get the 3D coordinate system handle. + call CCTK_CoordSystemHandle ( coord_system_handle, 'cart3d' ) if ( coord_system_handle .lt. 0) then - call CCTK_WARN( 0, "Cannot get handle for cart3d coordinate system. Forgot to activate an implementation providing coordinates ??" ) + warn_message = 'Cannot get handle for cart3d coordinate system. ' + warn_message = warn_message//'Forgot to activate an implementation ' + warn_message = warn_message//'providing coordinates??' + call CCTK_WARN( 0, warn_message ) endif - ! Find out how many interpolation points are located on this processor. - call CCTK_GrouplshGN ( status, cctkGH, 1, lsh, "ehfinder::xg" ) +! Find out how many interpolation points are located on this processor. + 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" ) + call CCTK_WARN ( 0, 'cannot get local size for surface arrays' ) end if - ! Set the pointers to the output arrays. +! Set the pointers to the output arrays. out_arrays(1) = CCTK_PointerTo(alpg) out_arrays(2) = CCTK_PointerTo(betaxg) out_arrays(3) = CCTK_PointerTo(betayg) @@ -92,51 +100,53 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS) out_arrays(13) = CCTK_PointerTo(dfzg) out_arrays(14) = CCTK_PointerTo(psig) - ! Loop over the level sets +! Loop over the level sets do l = 1, eh_number_level_sets - ! Set the pointers to the points to be interpolated to. +! 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" ) +! 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" ) + call CCTK_VarIndex ( in_arrays(12), 'staticconformal::psi' ) - ! Check the metric type. At present physical and static_conformal are - ! supported. +! 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. +! 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" ) + op_indices(1:13), 'operand_indices' ) if ( status .lt. 0 ) then - call CCTK_WARN ( 0, "Cannot set operand indices array in parameter table" ) + warn_message = 'Cannot set operand indices array in parameter table' + call CCTK_WARN ( 0, warn_message ) endif - ! Set the corresponding table entry for the operation codes. +! Set the corresponding table entry for the operation codes. call Util_TableSetIntArray ( status, table_handle, 13, & - op_codes(1:13), "operation_codes" ) + op_codes(1:13), 'operation_codes' ) if ( status .lt. 0 ) then - call CCTK_WARN ( 0, "Cannot set operation codes array in parameter table" ) + warn_message = 'Cannot set operation codes array in parameter table' + call CCTK_WARN ( 0, warn_message ) endif - ! Call the interpolator. +! Call the interpolator. call CCTK_InterpGridArrays ( status, cctkGH, 3, interp_handle, & table_handle, coord_system_handle, & lsh(1), CCTK_VARIABLE_REAL, & @@ -147,14 +157,14 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS) call CCTK_INFO ( 'Interpolation failed.' ) end if - ! For each point on this processor calculate the right hand side of the - ! characteristic evolution equation. +! 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. +! calculate the square of the lapse. alp2 = alpg(i)**2 - ! Calculate the inverse of the 3-metric. +! 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) @@ -169,32 +179,34 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS) 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. +! 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. +! Calculate the overall multiplication factor. factor = alp2 / sqrt ( alp2 * ( dfux * dfxg(i) + & dfuy * dfyg(i) + & dfuz * dfzg(i) ) ) - ! Finally obtain dx^i/dt. +! 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" ) + op_indices, 'operand_indices' ) if ( status .lt. 0 ) then - call CCTK_WARN ( 0, "Cannot set operand indices array in parameter table" ) + warn_message = 'Cannot set operand indices array in parameter table' + call CCTK_WARN ( 0, warn_message ) endif call Util_TableSetIntArray ( status, table_handle, 14, & - op_codes, "operation_codes" ) + op_codes, 'operation_codes' ) if ( status .lt. 0 ) then - call CCTK_WARN ( 0, "Cannot set operation codes array in parameter table" ) + warn_message = 'Cannot set operation codes array in parameter table' + call CCTK_WARN ( 0, warn_message ) endif call CCTK_InterpGridArrays ( status, cctkGH, 3, interp_handle, & @@ -217,7 +229,7 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS) guxy = gxzg(i) * gyzg(i) - gxyg(i) * gzzg(i) guxz = gxyg(i) * gyzg(i) - gxzg(i) * gyyg(i) -! The determinant divided by psi^4. +! The inverse of the determinant divided by psi^4. idetg = psi4 / ( gxxg(i) * guxx + & gxyg(i) * guxy + & gxzg(i) * guxz ) diff --git a/src/EHFinder_Generator_Sources2.F90 b/src/EHFinder_Generator_Sources2.F90 index b557381..27a9f90 100644 --- a/src/EHFinder_Generator_Sources2.F90 +++ b/src/EHFinder_Generator_Sources2.F90 @@ -19,6 +19,7 @@ subroutine EHFinder_Generator_Sources2(CCTK_ARGUMENTS) CCTK_INT :: interp_handle, table_handle, status, coord_system_handle character(len=200) :: gen_interp + character(len=128) :: warn_message CCTK_INT :: gen_interp_len character(len=7) :: gen_order @@ -32,82 +33,89 @@ subroutine EHFinder_Generator_Sources2(CCTK_ARGUMENTS) CCTK_REAL :: alp2, psi4, dfux, dfuy, dfuz, factor, ssign CCTK_REAL :: idetg, guxx, guxy, guxz, guyy, guyz, guzz +! Set the sign depending on the surface direction. if ( CCTK_EQUALS ( surface_direction, 'outward' ) ) ssign = one if ( CCTK_EQUALS ( surface_direction, 'inward' ) ) ssign = -one out_types = CCTK_VARIABLE_REAL - ! Convert the generator_interpolator string parameter to a Fortran string. +! Convert the generator_interpolator string parameter to a Fortran string. call CCTK_FortranString ( gen_interp_len, generator_interpolator, & gen_interp ) - ! Get the corresponding interpolator handle. +! Get the corresponding interpolator handle. call CCTK_InterpHandle ( interp_handle, gen_interp ) if ( interp_handle .lt. 0 ) then - call CCTK_WARN( 0, "Cannot get handle for interpolation. Forgot to activate an implementation providing interpolation operators??" ) + warn_message = 'Cannot get handle for interpolation. ' + warn_message = warn_message//'Forgot to activate an implementation ' + warn_message = warn_message//'providing interpolation operators??' + call CCTK_WARN( 0, warn_message ) end if - ! Convert the interpolation order parameter to a Fortran string to be placed - ! in the interpolator table. Note that the order is assumed to contain only - ! 1 digit. +! Convert the interpolation order parameter to a Fortran string to be placed +! in the interpolator table. Note that the order is assumed to contain only +! 1 digit. write(gen_order,'(a6,i1)') 'order=',generator_interpolation_order - ! Create the table directly from the string. +! Create the table directly from the string. call Util_TableCreateFromString ( table_handle, gen_order ) if ( table_handle .lt. 0 ) then - call CCTK_WARN( 0, "Cannot create parameter table for interpolator" ) + call CCTK_WARN( 0, 'Cannot create parameter table for interpolator' ) end if - ! Get the 3D coordinate system handle. - call CCTK_CoordSystemHandle ( coord_system_handle, "cart3d" ) +! Get the 3D coordinate system handle. + call CCTK_CoordSystemHandle ( coord_system_handle, 'cart3d' ) if ( coord_system_handle .lt. 0) then - call CCTK_WARN( 0, "Cannot get handle for cart3d coordinate system. Forgot to activate an implementation providing coordinates ??" ) + warn_message = 'Cannot get handle for cart3d coordinate system. ' + warn_message = warn_message//'Forgot to activate an implementation ' + warn_message = warn_message//'providing coordinates??' + call CCTK_WARN( 0, warn_message ) endif #include "include/physical_part.h" - ! Find out how many interpolation points are located on this processor. - call CCTK_GrouplshGN ( status, cctkGH, 1, lsh, "ehfinder::xg" ) +! Find out how many interpolation points are located on this processor. + 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" ) + call CCTK_WARN ( 0, 'cannot get local size for surface arrays' ) end if - ! Set the indices to the input grid functions. - call CCTK_VarIndex ( in_arrays(1), "ehfinder::xgf" ) - call CCTK_VarIndex ( in_arrays(2), "ehfinder::ygf" ) - call CCTK_VarIndex ( in_arrays(3), "ehfinder::zgf" ) +! Set the indices to the input grid functions. + call CCTK_VarIndex ( in_arrays(1), 'ehfinder::xgf' ) + call CCTK_VarIndex ( in_arrays(2), 'ehfinder::ygf' ) + call CCTK_VarIndex ( in_arrays(3), 'ehfinder::zgf' ) - ! Set the operand indices table entry, corresponding - ! to interpolation of ehfinder::generator_gf (3) +! Set the operand indices table entry, corresponding +! to interpolation of ehfinder::generator_gf (3) call Util_TableSetIntArray ( status, table_handle, 3, & - op_indices, "operand_indices" ) + op_indices, 'operand_indices' ) if ( status .lt. 0 ) then - call CCTK_WARN ( 0, "Cannot set operand indices array in parameter table" ) + call CCTK_WARN ( 0, 'Cannot set operand indices array in parameter table' ) endif - ! Set the corresponding table entry for the operation codes. +! Set the corresponding table entry for the operation codes. call Util_TableSetIntArray ( status, table_handle, 3, & - op_codes, "operation_codes" ) + op_codes, 'operation_codes' ) if ( status .lt. 0 ) then - call CCTK_WARN ( 0, "Cannot set operation codes array in parameter table" ) + call CCTK_WARN ( 0, 'Cannot set operation codes array in parameter table' ) endif - ! Loop over the level sets +! Loop over the level sets do l = 1, eh_number_level_sets - ! Set the pointers to the points to be interpolated to. +! 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. +! 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. +! Check the metric type. At present physical and static_conformal are +! supported. if ( CCTK_EQUALS ( metric_type, 'physical' ) ) then do k = kzl, kzr @@ -115,10 +123,10 @@ subroutine EHFinder_Generator_Sources2(CCTK_ARGUMENTS) do i = ixl, ixr if ( eh_mask(i,j,k,l) .ge. 0 ) then - ! calculate the square of the lapse. +! calculate the square of the lapse. alp2 = alp(i,j,k)**2 - ! Calculate the inverse of the 3-metric. +! 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) @@ -132,20 +140,25 @@ subroutine EHFinder_Generator_Sources2(CCTK_ARGUMENTS) 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 + 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. +! 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. +! 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 @@ -166,14 +179,14 @@ subroutine EHFinder_Generator_Sources2(CCTK_ARGUMENTS) if ( eh_mask(i,j,k,l) .ge. 0 ) then alp2 = alp(i,j,k)**2 - ! The inverse of psi^4 +! The inverse of psi^4 psi4 = one / psi(i,j,k)**4 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. +! The inverse of the determinant divided by psi^4. idetg = psi4 / ( gxx(i,j,k) * guxx + & gxy(i,j,k) * guxy + & gxz(i,j,k) * guxz ) @@ -185,12 +198,16 @@ subroutine EHFinder_Generator_Sources2(CCTK_ARGUMENTS) 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 + 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 * 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) + 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) ) ) @@ -207,7 +224,7 @@ subroutine EHFinder_Generator_Sources2(CCTK_ARGUMENTS) end do end if - ! Call the interpolator. +! Call the interpolator. call CCTK_InterpGridArrays ( status, cctkGH, 3, interp_handle, & table_handle, coord_system_handle, & lsh(1), CCTK_VARIABLE_REAL, & diff --git a/src/EHFinder_Init.F90 b/src/EHFinder_Init.F90 index c291376..32c50a4 100644 --- a/src/EHFinder_Init.F90 +++ b/src/EHFinder_Init.F90 @@ -23,25 +23,22 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS) CCTK_REAL :: theta, dtheta, thetamin, thetamax, r_el CCTK_INT, dimension(1) :: lsh, lbnd - ! Initialize the current_iteration variable. - current_iteration = last_iteration_number - - ! Initialize the last_time varible. Note the parameters should be - ! chosen to be consistent with the run producing the numerical data. - ! F.ex. if dt was 0.1 but the data was only stored every 4 iterations, the - ! parameters should be chosen so that dt now is 0.4. +! Initialize the last_time varible. Note the parameters should be +! chosen to be consistent with the run producing the numerical data. +! F.ex. if dt was 0.1 but the data was only stored every 4 iterations, the +! parameters should be chosen so that dt now is 0.4. last_time = abs(cctk_delta_time) * last_iteration_number / & saved_iteration_every cctk_time = last_time - ! Allocate the logical array containing the flag determining if the - ! corresponding level set should be re-initialized. +! 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 ) + if ( allocated(re_init_this_level_set) ) then + deallocate ( re_init_this_level_set ) end if - allocate ( reparam_this_level_set(eh_number_level_sets) ) + allocate ( re_init_this_level_set(eh_number_level_sets) ) if ( allocated(re_initialize_undone) ) then deallocate ( re_initialize_undone ) end if @@ -49,13 +46,13 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS) if ( evolve_generators .gt. 0 ) then - call CCTK_GrouplbndGN ( status, cctkGH, 1, lbnd, "ehfinder::xg" ) + 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" ) + call CCTK_WARN ( 0, 'cannot get lower bounds for generator arrays' ) end if - call CCTK_GrouplshGN ( status, cctkGH, 1, lsh, "ehfinder::xg" ) + 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" ) + call CCTK_WARN ( 0, 'cannot get local size for generator arrays' ) end if if ( CCTK_EQUALS( generator_distribution, 'line' ) ) then @@ -89,11 +86,11 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS) do l = 1, eh_number_level_sets - ! If a sphere is requested... +! 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. +! 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(l) )**2 + & ( y - translate_y(l) )**2 + & ( z - translate_z(l) )**2 ) - initial_rad(l) @@ -109,10 +106,10 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS) end if end if - ! If an ellipsoid is requested... +! If an ellipsoid is requested... if ( CCTK_EQUALS( initial_f(l), 'ellipsoid' ) ) then - ! Calculate sines and cosines of the rotation parameters. +! Calculate sines and cosines of the rotation parameters. cosa = cos(rotation_alpha(l)) sina = sin(rotation_alpha(l)) cosb = cos(rotation_beta(l)) @@ -120,8 +117,8 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS) cosc = cos(rotation_gamma(l)) sinc = sin(rotation_gamma(l)) - ! 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. +! 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 @@ -131,11 +128,10 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS) 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. +! 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 @@ -168,14 +164,14 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS) end if end if - ! if an ovaloid of Cassini is requested... +! 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(l)**4 - & 2*cas_a(l)**2*(x**2 - (y**2+z**2)) - cas_b(l)**4 end if end do - ! Initialise the internal mask. +! Initialise the internal mask. eh_mask = 0 return @@ -191,26 +187,13 @@ subroutine EHFinder_Init(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS - ! Set up the value used in interiour inactive cells. +! Set up the value used in interiour inactive cells. ex_value = - ( one + shell_width ) * maxval(cctk_delta_space) - ! Initialize nx, ny and nz based on the local gris size. They are defined - ! in the module EHFinder_mod and will therefore be known from now on by all - ! routine that uses this module. - nx = cctk_lsh(1) - ny = cctk_lsh(2) - nz = cctk_lsh(3) - - ! Find the maximal grid spacing. +! Find the maximal grid spacing. delta = maxval ( cctk_delta_space ) - ! Initialise the ghostzone information. These variables are defined in - ! the module EHFinder_mod. - ngx = cctk_nghostzones(1) - ngy = cctk_nghostzones(2) - ngz = cctk_nghostzones(3) - - ! Get handles for various reduction operations. +! Get handles for various reduction operations. call CCTK_ReductionArrayHandle ( max_handle, 'maximum' ) if ( max_handle .lt. 0 ) then call CCTK_WARN(0,'Could not obtain a handle for maximum reduction') @@ -224,24 +207,4 @@ subroutine EHFinder_Init(CCTK_ARGUMENTS) call CCTK_WARN(0,'Could not obtain a handle for sum reduction') end if -! ! Register a coodinatsystem for the 2D surface grid arrays. This has no -! ! effect yet, but should do something useful, when TAGS can be used to -! ! assign coordinates to grid arrays. -! call CCTK_CoordRegisterSystem ( ierr, 2, 'cart2d' ) -! if ( ierr .lt. 0 ) then -! call CCTK_WARN(1,'Could not register a 2D coordinate system as "cart2d"') -! end if -! call CCTK_CoordRegisterData ( ierr, 1, 'ehfinder::ctheta', & -! 'x', 'cart2d' ) -! if ( ierr .lt. 0 ) then -! call CCTK_WARN(1,'Could not register ctheta as the 1st coordinate for "cart2d"') -! end if -! call CCTK_CoordRegisterData ( ierr, 2, 'ehfinder::cphi', & -! 'y', 'cart2d' ) -! if ( ierr .lt. 0 ) then -! call CCTK_WARN(1,'Could not register cphi as the 2nd coordinate for "cart2d"') -! end if -! -! call CCTK_INFO('2d coordinate system registered') - end subroutine EHFinder_Init diff --git a/src/EHFinder_Integrate.F90 b/src/EHFinder_Integrate.F90 new file mode 100644 index 0000000..db2a315 --- /dev/null +++ b/src/EHFinder_Integrate.F90 @@ -0,0 +1,704 @@ +! Integrating over the surface(s). Right now only in full mode. +! $Header$ + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" + +subroutine EHFinder_FindSurfaceElement(CCTK_ARGUMENTS) + + use EHFinder_mod + + implicit none + + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + + CCTK_INT :: i, j, k, l, im, jm + CCTK_INT :: interp_handle, table_handle, coord_system_handle + + character(len=200) :: area_interp + character(len=128) :: warn_message + CCTK_INT :: area_interp_len + character(len=7) :: area_order + + CCTK_INT, dimension(4) :: bbox + CCTK_INT, dimension(2) :: gsh, lsh, lbnd, ubnd, nghost + + CCTK_REAL :: dtheta, dphi, dthetainv, dphiinv + CCTK_REAL :: dxdth, dxdph, dydth, dydph, dzdth, dzdph + CCTK_POINTER, dimension(3) :: interp_coords + CCTK_POINTER, dimension(7) :: out_array + CCTK_INT, dimension(7) :: in_array + CCTK_INT, dimension(7), parameter :: out_types = (/ CCTK_VARIABLE_REAL, & + CCTK_VARIABLE_REAL, & + CCTK_VARIABLE_REAL, & + CCTK_VARIABLE_REAL, & + CCTK_VARIABLE_REAL, & + CCTK_VARIABLE_REAL, & + CCTK_VARIABLE_REAL /) + +! If finding of surface failed do not try to integrate but exit. + if ( find_surface_status .lt. 0 ) then + return + endif + +! Store the levelset_counter in a shorter named variable for convenience. + l = levelset_counter + +! Write an Info message about what we are doing. + call CCTK_INFO ( 'Finding surface element' ) + +! Obtain the dimensions and characteristics of the grid arrays. + 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 + 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 + 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 + 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 + 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 + 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' ) + end if + +! Setup the theta and phi spacings and their inverse. + dtheta = ctheta(2,1) - ctheta(1,1) + dphi = cphi(1,2) - cphi(1,1) + dthetainv = one / dtheta + dphiinv = one / dphi + +! Convert the interpolator parameter into a fortran string. + call CCTK_FortranString ( area_interp_len, area_interpolator, & + area_interp ) + +! Get an interpolation handle. + call CCTK_InterpHandle ( interp_handle, area_interp(1:area_interp_len) ) + + if ( interp_handle .lt. 0 ) then + warn_message = 'Cannot get handle for interpolation. ' + warn_message = warn_message//'Forgot to activate an implementation ' + warn_message = warn_message//'providing interpolation operators??' + call CCTK_WARN( 0, warn_message ) + end if + +! Write the interpolation order parameter into a string + write(area_order,'(a6,i1)') 'order=',area_interpolation_order + +! Create a table from the interpolation order string. + call Util_TableCreateFromString ( table_handle, area_order ) + if ( table_handle .lt. 0 ) then + call CCTK_WARN( 0, 'Cannot create parameter table for interpolator' ) + end if + +! Get a handle for the coordinate system. + call CCTK_CoordSystemHandle ( coord_system_handle, 'cart3d' ) + if ( coord_system_handle .lt. 0) then + warn_message = 'Cannot get handle for cart3d coordinate system. ' + warn_message = warn_message//'Forgot to activate an implementation ' + warn_message = warn_message//'providing coordinates??' + call CCTK_WARN( 0, warn_message ) + endif + +! Get the pointers to the interpolation points. + interp_coords(1) = CCTK_PointerTo(interp_x) + interp_coords(2) = CCTK_PointerTo(interp_y) + interp_coords(3) = CCTK_PointerTo(interp_z) + +! Get the pointers to the interpolation return arrays. + out_array(1) = CCTK_PointerTo(gxxi) + out_array(2) = CCTK_PointerTo(gxyi) + out_array(3) = CCTK_PointerTo(gxzi) + out_array(4) = CCTK_PointerTo(gyyi) + out_array(5) = CCTK_PointerTo(gyzi) + out_array(6) = CCTK_PointerTo(gzzi) + out_array(7) = CCTK_PointerTo(psii) + +! find the cartesian coordinates for the interpolation points + 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) + interp_y(i,j) = center(2) + rsurf(i,j) * sintheta(i,j) * sinphi(i,j) + interp_z(i,j) = center(3) + rsurf(i,j) * costheta(i,j) + end do + end do + +! Get the indices for the grid functions to be interpolated. + call CCTK_VarIndex ( in_array(1), 'admbase::gxx' ) + call CCTK_VarIndex ( in_array(2), 'admbase::gxy' ) + call CCTK_VarIndex ( in_array(3), 'admbase::gxz' ) + call CCTK_VarIndex ( in_array(4), 'admbase::gyy' ) + call CCTK_VarIndex ( in_array(5), 'admbase::gyz' ) + call CCTK_VarIndex ( in_array(6), 'admbase::gzz' ) + +! If the static conformal factor is used, we also need to interpolate it. + if ( CCTK_EQUALS ( metric_type, 'static conformal' ) ) then + call CCTK_VarIndex ( in_array(7), 'staticconformal::psi' ) + +! Call the interpolator for the metric and the conformal factor. + call CCTK_InterpGridArrays ( ierr, cctkGH, 3, interp_handle, & + table_handle, coord_system_handle, & + lsh(1) * lsh(2), CCTK_VARIABLE_REAL, & + interp_coords, 7, in_array, & + 7, out_types, out_array ) + +! Get the physical metric. + gxxi = psii**4 * gxxi + gxyi = psii**4 * gxyi + gxzi = psii**4 * gxzi + gyyi = psii**4 * gyyi + gyzi = psii**4 * gyzi + gzzi = psii**4 * gzzi + else +! Else call the interpolater only for the metric. + call CCTK_InterpGridArrays ( ierr, cctkGH, 3, interp_handle, & + table_handle, coord_system_handle, & + lsh(1) * lsh(2), CCTK_VARIABLE_REAL, & + interp_coords, 6, in_array(1:6), & + 6, out_types(1:6), out_array(1:6) ) + end if + +! For all points on the surface... + do j = 1, lsh(2) + do i = 1, lsh(1) + +! Calculate the coordinate transformation between the 3-d cartesian +! coordinates and the 2-d (theta,phi) coordinates on the surface. + dxdth = ( drdtheta(i,j) * sintheta(i,j) + & + 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) + dzdth = drdtheta(i,j) * costheta(i,j) - rsurf(i,j) * sintheta(i,j) + dxdph = ( drdphi(i,j) * cosphi(i,j) - & + rsurf(i,j) * sinphi(i,j) ) * sintheta(i,j) + dydph = ( drdphi(i,j) * sinphi(i,j) + & + rsurf(i,j) * cosphi(i,j) ) * sintheta(i,j) + dzdph = drdphi(i,j) * costheta(i,j) + +! Using this coordinate transformation, calculate the induced 2-metric +! on the surface. + gtt(i,j) = dxdth**2 * gxxi(i,j) + dydth**2 * gyyi(i,j) + & + dzdth**2 * gzzi(i,j) + & + two * ( dxdth * dydth * gxyi(i,j) + & + dxdth * dzdth * gxzi(i,j) + & + dydth * dzdth * gyzi(i,j) ) + gtp(i,j) = dxdth * dxdph * gxxi(i,j) + & + ( dxdth * dydph + dydth * dxdph ) * gxyi(i,j) + & + ( dxdth * dzdph + dzdth * dxdph ) * gxzi(i,j) + & + dydth * dydph * gyyi(i,j) + & + ( dydth * dzdph + dzdth * dydph ) * gyzi(i,j) + & + dzdth * dzdph * gzzi(i,j) + gpp(i,j) = dxdph**2 * gxxi(i,j) + dydph**2 * gyyi(i,j) + & + dzdph**2 * gzzi(i,j) + & + two * ( dxdph * dydph * gxyi(i,j) + & + dxdph * dzdph * gxzi(i,j) + & + dydph * dzdph * gyzi(i,j) ) + +! Calculate the area element as the square-root of the determinant +! of the two metric multiplied by dtheta*dphi. + da(i,j) = dtheta * dphi * sqrt ( gtt(i,j) * gpp(i,j) - gtp(i,j)**2 ) + +! Calculate the arc length element for the polar circumference integration. + dltheta(i,j) = dtheta * sqrt ( gtt(i,j) ) + +! Calculate the arc length element for the equatorial circumference +! integration. + dlphi(i,j) = dphi * sqrt ( gpp(i,j) ) + + end do + end do + +end subroutine EHFinder_FindSurfaceElement + + +subroutine EHFinder_IntegrateArea(CCTK_ARGUMENTS) + + use EHFinder_mod + + implicit none + + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + + CCTK_INT :: int_var, sn, l + CCTK_REAL :: eh_area_tmp + character(len=30) :: info_message + CCTK_INT, dimension(1) :: lbnd, ubnd, lsh + +! If finding of surface failed do not try to integrate but exit. + if ( find_surface_status .lt. 0 ) then + return + endif + +! Store the levelset_counter in a shorter named variable for convenience. + l = levelset_counter + +! Write an info message, to indicate what we are doing. + call CCTK_INFO ( 'Integrating area' ) + + sn = surface_counter - integrate_counter + +! Taking into account the sym_factor and the weights, the area is simply +! The sum of the elements of this array. + int_tmp = sym_factor * weights * da + +! Get the index to the area integration array. + call CCTK_VarIndex ( int_var, 'ehfinder::int_tmp' ) + if ( int_var .lt. 0 ) then + call CCTK_WARN ( 0, 'Could not get index to grid array int_tmp' ) + end if + +! Do a sum reduce over all processors. The result is the area. + call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, & + eh_area_tmp, 1, int_var ) + +! Write an info message reporting the area found. + write(info_message,'(a15,f10.5)') 'Horizon area = ',eh_area_tmp + + call CCTK_INFO(info_message) + +! Figure out the bounds of the area storage grid array. + call CCTK_GrouplbndGN ( ierr, cctkGH, 1, lbnd, 'ehfinder::eh_area2' ) + if ( ierr .lt. 0 ) then + call CCTK_WARN ( 0, 'cannot get lower bounds for area array' ) + end if + call CCTK_GroupubndGN ( ierr, cctkGH, 1, ubnd, 'ehfinder::eh_area2' ) + if ( ierr .lt. 0 ) then + call CCTK_WARN ( 0, 'cannot get upper bounds for area array' ) + end if + call CCTK_GrouplshGN ( ierr, cctkGH, 1, lsh, 'ehfinder::eh_area2' ) + if ( ierr .lt. 0 ) then + call CCTK_WARN ( 0, 'cannot get local size for area array' ) + end if + +! If this processor contains the right part of the area storage grid array +! save the area, otherwise do nothing. + if ( ( sn .ge. lbnd(1) + 1 ) .and. ( sn .le. ubnd(1) + 1 ) ) then + eh_area2(sn-lbnd(1),l) = eh_area_tmp + end if + +end subroutine EHFinder_IntegrateArea + + +subroutine EHFinder_IntegrateCentroid(CCTK_ARGUMENTS) + + use EHFinder_mod + + implicit none + + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + + 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 + +! If finding of surface failed do not try to integrate but exit. + if ( find_surface_status .lt. 0 ) then + return + endif + +! Store the levelset_counter in a shorter named variable for convenience. + l = levelset_counter + +! Write an info message, to indicate what we are doing. + call CCTK_INFO ( 'Integrating centroid' ) + + sn = surface_counter - integrate_counter + +! Calculate the (x,y,z) coordinates of the points on the surface. + interp_x = center(1) + rsurf * sintheta * cosphi + interp_y = center(2) + rsurf * sintheta * sinphi + interp_z = center(3) + rsurf * costheta + +! Get the index to the temporary integration array. + call CCTK_VarIndex ( int_var, 'ehfinder::int_tmp' ) + if ( int_var .lt. 0 ) then + call CCTK_WARN ( 0, 'Could not get index to grid array int_tmp' ) + end if + +! Take into account the symmetry, weights, x-coordinate and area element +! to calculate the integral of x*da over the surface. This gives +! centroid_x * Area, so later we divide by the area. + int_tmp = sym_factor * interp_x * weights * da + +! Sum up the all elements of int_tmp. + call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, & + centroid_x, 1, int_var ) + +! Take into account the symmetry, weights, y-coordinate and area element +! to calculate the integral of y*da over the surface. This gives +! centroid_y * Area, so later we divide by the area. + int_tmp = sym_factor * interp_y * weights * da + +! Sum up the all elements of int_tmp. + call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, & + centroid_y, 1, int_var ) + +! Take into account the symmetry, weights, z-coordinate and area element +! to calculate the integral of z*da over the surface. This gives +! centroid_z * Area, so later we divide by the area. + int_tmp = sym_factor * interp_z * weights * da + +! Sum up the all elements of int_tmp. + call CCTK_Reduce ( ierr, cctkGH, -1, sum_handle, 1, CCTK_VARIABLE_REAL, & + centroid_z, 1, int_var ) +! If this surface has symmetries set the corresponding centroids to zero. + if ( s_symx ) centroid_x = zero + if ( s_symy ) centroid_y = zero + if ( s_symz ) centroid_z = zero + +! Get the bounds for the centroid arrays. It should be enough to do +! it just for the x-array, since they are all defined in the same way. + 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_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_x' ) + if ( ierr .lt. 0 ) then + call CCTK_WARN ( 0, 'cannot get local size for centroid array' ) + end if + +! If we are on the right processor, store the centroid, otherwise do nothing. + if ( ( sn .ge. lbnd(1) + 1 ) .and. ( sn .le. ubnd(1) + 1 ) ) then + k = sn - lbnd(1) + 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 an info message with the calculated centroids. + write(info_message,'(a19,3(f10.5))') 'Horizon centroid = ', & + eh_centroid2_x(sn,l), & + eh_centroid2_y(sn,l), & + eh_centroid2_z(sn,l) + call CCTK_INFO(info_message) + +end subroutine EHFinder_IntegrateCentroid + + +subroutine EHFinder_IntegrateCircumference(CCTK_ARGUMENTS) + + use EHFinder_mod + + implicit none + + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + + 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 + CCTK_INT, dimension(4) :: bbox + CCTK_INT :: ithl, ithr, jphl, jphr, itha, jpha + CCTK_REAL :: eq_circ_loc, eq_circ, pol_circ_loc, pol_circ + +! If finding of surface failed do not try to integrate but exit. + if ( find_surface_status .lt. 0 ) then + return + endif + +! Store the levelset_counter in a shorter named variable for convenience. + l = levelset_counter + + sn = surface_counter - integrate_counter + + ! Find out the lower bounds of the distributed integration grid arrays. + 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 + + ! Find out the local size of the distributed integration grid arrays + 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 + + ! Find out the bounding box of the distributed integration grid arrays + 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 + + ! Find out the ghost size of the distributed integration grid arrays + call CCTK_GroupnghostzonesGN ( ierr, cctkGH, 2, nghost, & + 'ehfinder::surface_arrays' ) + if ( ierr .lt. 0 ) then + call CCTK_WARN ( 0, 'cannot get ghost zones for surface arrays' ) + end if + +! Initialize the lower and upper local index for the theta and phi arrays. + ithl = 1; ithr = lsh(1) + jphl = 1; jphr = lsh(2) + +! If the boundaries are internal, adjust the lower and upper local index +! to take into account the ghost cells. + if ( bbox(1) .eq. 0 ) ithl = ithl + nghost(1) + if ( bbox(2) .eq. 0 ) ithr = ithr - nghost(1) + if ( bbox(3) .eq. 0 ) jphl = jphl + nghost(2) + if ( bbox(4) .eq. 0 ) jphr = jphr - nghost(2) + +! Write and info message about what we are doing. + call CCTK_INFO ( 'Integrating equatorial circumference' ) + +! Setup the integration array for the equatorial circumference +! calculation. Note this is still a 2d-array, so we have to +! do the sum only for the correct 1d-slice of the array. + int_tmp = phi_sym_factor * phiweights * dlphi + +! If we are on a processor that contains part of the correct 1-d slice +! sum up the appropriate elemets, otherwise this processor will +! contribute zero. + if ( ( ithl+lbnd(1) .le. githeta ) .and. ( githeta .le. ithr+lbnd(1) ) ) then + itha = githeta - lbnd(1) + eq_circ_loc = sum ( int_tmp(itha,jphl:jphr) ) + else + eq_circ_loc = zero + end if + +! Then reduce the results over all processors and send the result to +! all processors. + call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, sum_handle, & + eq_circ_loc, eq_circ, CCTK_VARIABLE_REAL ) + +! Get the upper and lower bounds for the circumference arrays. + 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_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_circ_eq2' ) + if ( ierr .lt. 0 ) then + call CCTK_WARN ( 0, 'cannot get local size for area array' ) + end if + +! If we are on the correct processor store the result. + if ( ( sn .ge. lbnd1(1) + 1 ) .and. ( sn .le. ubnd1(1) + 1 ) ) then + k = sn - lbnd1(1) + eh_circ_eq2(k,l) = eq_circ + end if + +! Write an info message with the calculated circunference. + write(info_message,'(a35,f10.5)') 'Horizon equatorial circumference = ', & + eq_circ + call CCTK_INFO(info_message) + +! Now do the same for the polar circomference. + call CCTK_INFO ( 'Integrating polar circumference' ) + +! Setup the integration array for the polar circumference +! calculation. Note this is still a 2d-array, so we have to +! do the sum only for the correct 1d-slice of the array. + int_tmp = theta_sym_factor * thetaweights * dltheta + +! If we are on a processor that contains part of the correct 1-d slice +! sum up the appropriate elemets, otherwise this processor will +! contribute zero. + 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) ) + else + pol_circ_loc = zero + end if + +! Then reduce the results over all processors and send the result to +! all processors. + call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, sum_handle, & + pol_circ_loc, pol_circ, CCTK_VARIABLE_REAL ) + + +! If we are on the correct processor store the result. + if ( ( sn .ge. lbnd1(1) + 1 ) .and. ( sn .le. ubnd1(1) + 1 ) ) then + k = sn - lbnd1(1) + eh_circ_pol2(k,l) = pol_circ + end if + +! Write an info message with the calculated circumference. + write(info_message,'(a30,f10.5)') 'Horizon polar circumference = ',pol_circ + call CCTK_INFO(info_message) + +end subroutine EHFinder_IntegrateCircumference + + +! This routine sets up the weights for the Simpsons rule integration +! over the surface. +subroutine EHFinder_InitWeights(CCTK_ARGUMENTS) + + use EHFinder_mod + + implicit none + + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + + CCTK_INT :: i, j, k, im, jm + + CCTK_INT, dimension(4) :: bbox + CCTK_INT, dimension(2) :: gsh, lsh, lbnd, ubnd, nghost + +! Find out the lower bounds of the distributed integration grid arrays. + 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 + +! Find out the local size of the distributed integration grid arrays + 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 + + weights = one + +! Initialise the weight grid array for the 2D Simpsons rule integration. +! To do this I need to figure out the global location of the given point. +! There are 3 cases in the one dimensional case. If the point is on the +! boundary the weight is 1/3. If it is at an even position the weight +! is 4/3 and if it is at an odd position the weight is 2/3. + + do j = 1, lsh(2) + +! This is first done in the phi direction. Meaning that all points with +! the same theta coordinate are set to the same weight. + if ( ( lbnd(2)+j .eq. 1 ) .or. ( lbnd(2)+j .eq. nphi ) ) then + weights(:,j) = onethird + phiweights(:,j) = onethird + else if ( mod(lbnd(2)+j,2) .eq. 0 ) then + weights(:,j) = fourthirds + phiweights(:,j) = fourthirds + else + weights(:,j) = twothirds + phiweights(:,j) = twothirds + end if + +! Then it is done in the theta direction with the one-dimensional +! weights beeing multiplied. + do i = 1, lsh(1) + if ( ( lbnd(1)+i .eq. 1 ) .or. ( lbnd(1)+i .eq. ntheta ) ) then + weights(i,j) = onethird * weights(i,j) + thetaweights(i,j) = onethird + else if ( mod(lbnd(1)+i,2) .eq. 0 ) then + weights(i,j) = fourthirds * weights(i,j) + thetaweights(i,j) = fourthirds + else + weights(i,j) = twothirds * weights(i,j) + thetaweights(i,j) = twothirds + end if + end do + end do + +! The end result is a 2D array with the weights in the following pattern. + +! 1/9 4/9 2/9 4/9 2/9 4/9 1/9 +! 4/9 16/9 8/9 16/9 8/9 16/9 4/9 +! 2/9 8/9 4/9 8/9 4/9 8/9 2/9 +! 4/9 16/9 8/9 16/9 8/9 16/9 4/9 +! 2/9 8/9 4/9 8/9 4/9 8/9 2/9 +! 4/9 16/9 8/9 16/9 8/9 16/9 4/9 +! 1/9 4/9 2/9 4/9 2/9 4/9 1/9 +end subroutine EHFinder_InitWeights + + +! This is the routine where the stored areas (in eh_area2) are copied +! into the trigger variable eh_area. +subroutine EHFinder_CopyArea(CCTK_ARGUMENTS) + + use EHFinder_mod + + implicit none + + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + +! If finding of surface failed set area to zero. + if ( find_surface_status .lt. 0 ) then + eh_area = zero + return + else + call CCTK_INFO('Copying area data') + eh_area = eh_area2 + endif +end subroutine EHFinder_CopyArea + + +! This is the routine where the stored centroids (in eh_centroid2_[xyz]) +! are copied into the trigger variable eh_centroid_[xyz]. +subroutine EHFinder_CopyCentroid(CCTK_ARGUMENTS) + + use EHFinder_mod + + implicit none + + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + +! If finding of surface failed set centroids to zero. + if ( find_surface_status .lt. 0 ) then + eh_centroid_x = zero + eh_centroid_y = zero + eh_centroid_z = zero + return + else + call CCTK_INFO('Copying centroid data') + eh_centroid_x = eh_centroid2_x + eh_centroid_y = eh_centroid2_y + eh_centroid_z = eh_centroid2_z + end if + +end subroutine EHFinder_CopyCentroid + + +! This is the routine where the stored circumferences (in eh_circ_[eq,pol]2) +! are copied into the trigger variable eh_circ_[eq,pol]. +subroutine EHFinder_CopyCircumference(CCTK_ARGUMENTS) + + use EHFinder_mod + + implicit none + + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + + if ( find_surface_status .lt. 0 ) then + eh_circ_eq = zero + eh_circ_pol = zero + return + else + call CCTK_INFO('Copying circumference data') + eh_circ_eq = eh_circ_eq2 + eh_circ_pol = eh_circ_pol2 + end if + +end subroutine EHFinder_CopyCircumference diff --git a/src/EHFinder_ParamCheck.F90 b/src/EHFinder_ParamCheck.F90 index ae9ec1e..9640d45 100644 --- a/src/EHFinder_ParamCheck.F90 +++ b/src/EHFinder_ParamCheck.F90 @@ -12,40 +12,21 @@ #include "cctk_Arguments.h" #include "cctk_Parameters.h" -! @@ -! @routine EHFinder_ParamCheck -! @date Tue May 21 21:26:45 2002 -! @author Peter Diener -! @desc -! Scheduled routine to detect invalid parameter settings. -! @enddesc -! @calls -! @calledby -! @history -! -! @endhistory -! -! @@ - subroutine EHFinder_ParamCheck(CCTK_ARGUMENTS) DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS + character(len=128) :: param_warn + ! Check if the metric is one of the known types: physical or ! static conformal if ( ( .not. CCTK_EQUALS(metric_type,'physical') ) .and. & ( .not. CCTK_EQUALS(metric_type,'static conformal') ) ) then - call CCTK_PARAMWARN('Unknown ADMBase::metric_type - known types are "physical" and "static conformal"') + param_warn = 'Unknown ADMBase::metric_type - known types are "physical" ' + param_warn = param_warn//'and "static conformal"' + call CCTK_PARAMWARN ( param_warn ) end if - ! Check if ntheta and nphi are odd as required for the Simpsons rulei - ! integration. -! if ( mod(ntheta,2) .eq. 0 ) then -! call CCTK_PARAMWARN('EHFinder::ntheta has to be an odd number') -! end if -! if ( mod(nphi,2) .eq. 0 ) then -! call CCTK_PARAMWARN('EHFinder::nphi has to be an odd number') -! end if end subroutine EHFinder_ParamCheck diff --git a/src/EHFinder_ReadData.F90 b/src/EHFinder_ReadData.F90 index 08e3ce8..0185f5b 100644 --- a/src/EHFinder_ReadData.F90 +++ b/src/EHFinder_ReadData.F90 @@ -21,25 +21,24 @@ subroutine EHFinder_Read_Metric(CCTK_ARGUMENTS) character(len=10) :: iteration_string CCTK_INT :: i, nc, ntot, res - ! Figure out which iteration number to read, based on the parameters - ! last_iteration_number (the last iteration in the numerical evolution - ! producing the metric) and saved_iteration_every (how often was the metric - ! saved) and the current iteration and save it in a string variable. - +! Figure out which iteration number to read, based on the parameters +! last_iteration_number (the last iteration in the numerical evolution +! producing the metric) and saved_iteration_every (how often was the metric +! saved) and the current iteration and save it in a string variable. write(iteration_string,'(i10)') last_iteration_number - & saved_iteration_every * cctk_iteration - ! Trim the string variable. +! Trim the string variable. iteration_string = adjustl(iteration_string) nc = len_trim(iteration_string) - ! Generate a string with the filenames. Note at present the requirement - ! therefore is that all iterations of one variable is written in the - ! same file. +! Generate a string with the filenames. Note at present the requirement +! therefore is that all iterations of one variable is written in the +! same file. in_files = 'gxx gxy gxz gyy gyz gzz' - ! Fill in the string array, used to requesting the metric components at - ! the specified iteration number. +! Fill in the string array, used to requesting the metric components at +! the specified iteration number. var_names(1) = 'admbase::gxx[cctk_iteration='//iteration_string(1:nc)//']' var_names(2) = 'admbase::gxy[cctk_iteration='//iteration_string(1:nc)//']' var_names(3) = 'admbase::gxz[cctk_iteration='//iteration_string(1:nc)//']' @@ -47,7 +46,7 @@ subroutine EHFinder_Read_Metric(CCTK_ARGUMENTS) var_names(5) = 'admbase::gyz[cctk_iteration='//iteration_string(1:nc)//']' var_names(6) = 'admbase::gzz[cctk_iteration='//iteration_string(1:nc)//']' - ! merge all the varible names into a single string. +! merge all the variable names into a single string. in_vars = ' ' ntot = 0 do i = 1, 6 @@ -56,8 +55,8 @@ subroutine EHFinder_Read_Metric(CCTK_ARGUMENTS) ntot = ntot + nc + 1 end do - ! Call the routine that actualle does the file access. Note failures are - ! just silently ignored. +! Call the routine that actualle does the file access. Note failures are +! just silently ignored. call IOUtil_RecoverVarsFromDatafiles ( res, cctkGH, in_files, in_vars ) end subroutine EHFinder_Read_Metric @@ -78,14 +77,28 @@ subroutine EHFinder_Read_Lapse(CCTK_ARGUMENTS) character(len=10) :: iteration_string CCTK_INT :: nc, res +! Figure out which iteration number to read, based on the parameters +! last_iteration_number (the last iteration in the numerical evolution +! producing the metric) and saved_iteration_every (how often was the lapse +! saved) and the current iteration and save it in a string variable. write(iteration_string,'(i10)') last_iteration_number - & saved_iteration_every * cctk_iteration + +! Trim the string variable. iteration_string = adjustl(iteration_string) nc = len_trim(iteration_string) +! Generate a string with the filename. Note at present the requirement +! therefore is that all iterations of one variable is written in the +! same file. in_files = 'alp' + +! Fill in the string used to requesting the lapse at the specified +! iteration number. in_vars = 'admbase::alp[cctk_iteration='//iteration_string(1:nc)//']' +! Call the routine that actualle does the file access. Note failures are +! just silently ignored. call IOUtil_RecoverVarsFromDatafiles ( res, cctkGH, in_files, in_vars ) end subroutine EHFinder_Read_Lapse @@ -107,15 +120,28 @@ subroutine EHFinder_Read_Shift(CCTK_ARGUMENTS) character(len=10) :: iteration_string CCTK_INT :: i, nc, ntot, res +! Figure out which iteration number to read, based on the parameters +! last_iteration_number (the last iteration in the numerical evolution +! producing the metric) and saved_iteration_every (how often was the metric +! saved) and the current iteration and save it in a string variable. write(iteration_string,'(i10)') last_iteration_number - & saved_iteration_every * cctk_iteration +! Trim the string variable. iteration_string = adjustl(iteration_string) nc = len_trim(iteration_string) +! Generate a string with the filenames. Note at present the requirement +! therefore is that all iterations of one variable is written in the +! same file. in_files = 'betax betay betaz' + +! Fill in the string array, used to requesting the shift components at +! the specified iteration number. var_names(1) = 'admbase::betax[cctk_iteration='//iteration_string(1:nc)//']' var_names(2) = 'admbase::betay[cctk_iteration='//iteration_string(1:nc)//']' var_names(3) = 'admbase::betaz[cctk_iteration='//iteration_string(1:nc)//']' + +! merge all the variable names into a single string. in_vars = ' ' ntot = 0 do i = 1, 3 @@ -124,6 +150,8 @@ subroutine EHFinder_Read_Shift(CCTK_ARGUMENTS) ntot = ntot + nc + 1 end do +! Call the routine that actualle does the file access. Note failures are +! just silently ignored. call IOUtil_RecoverVarsFromDatafiles ( res, cctkGH, in_files, in_vars ) end subroutine EHFinder_Read_Shift @@ -144,18 +172,34 @@ subroutine EHFinder_Read_Conformal(CCTK_ARGUMENTS) character(len=10) :: iteration_string CCTK_INT :: nc, res +! Figure out which iteration number to read, based on the parameters +! last_iteration_number (the last iteration in the numerical evolution +! producing the metric) and saved_iteration_every (how often was the metric +! saved) and the current iteration and save it in a string variable. Note +! that for this variable the default option is to only to read it in once at +! the first timestep. if ( read_conformal_factor_once .gt. 0 ) then write(iteration_string,'(i10)') 0 else write(iteration_string,'(i10)') last_iteration_number - & saved_iteration_every * cctk_iteration end if + +! Trim the string variable. iteration_string = adjustl(iteration_string) nc = len_trim(iteration_string) +! Generate a string with the filenames. Note at present the requirement +! therefore is that all iterations of one variable is written in the +! same file. in_files = 'psi' + +! Fill in the string array, used to requesting the conformal factor at +! the specified iteration number. in_vars = 'staticconformal::psi[cctk_iteration='//iteration_string(1:nc)//']' +! Call the routine that actualle does the file access. Note failures are +! just silently ignored. call IOUtil_RecoverVarsFromDatafiles ( res, cctkGH, in_files, in_vars ) end subroutine EHFinder_Read_Conformal @@ -176,66 +220,28 @@ subroutine EHFinder_Read_Mask(CCTK_ARGUMENTS) character(len=10) :: iteration_string CCTK_INT :: nc, res +! Figure out which iteration number to read, based on the parameters +! last_iteration_number (the last iteration in the numerical evolution +! producing the metric) and saved_iteration_every (how often was the metric +! saved) and the current iteration and save it in a string variable. write(iteration_string,'(i10)') last_iteration_number - & saved_iteration_every * cctk_iteration + +! Trim the string variable. iteration_string = adjustl(iteration_string) nc = len_trim(iteration_string) +! Generate a string with the filename. Note at present the requirement +! therefore is that all iterations of one variable is written in the +! same file. in_files = 'emask' + +! Fill in the string array, used to requesting the old style mask at +! the specified iteration number. in_vars = 'spacemask::emask[cctk_iteration='//iteration_string(1:nc)//']' +! Call the routine that actualle does the file access. Note failures are +! just silently ignored. call IOUtil_RecoverVarsFromDatafiles ( res, cctkGH, in_files, in_vars ) end subroutine EHFinder_Read_Mask - - -! This routine reads in everything. -subroutine EHFinder_ReadData(CCTK_ARGUMENTS) - - use EHFinder_mod - - implicit none - - DECLARE_CCTK_PARAMETERS - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_FUNCTIONS - - character(len=1024) :: in_files, in_vars - character(len=64), dimension(11) :: var_names - character(len=10) :: iteration_string - CCTK_INT :: i, nc, ntot, res - - write(iteration_string,'(i10)') last_iteration_number - & - saved_iteration_every * cctk_iteration - iteration_string = adjustl(iteration_string) - nc = len_trim(iteration_string) - - in_files = 'alp betax betay betaz gxx gxy gxz gyy gyz gzz psi' - var_names(1) = 'admbase::alp[cctk_iteration='//iteration_string(1:nc)//']' - var_names(2) = 'admbase::betax[cctk_iteration='//iteration_string(1:nc)//']' - var_names(3) = 'admbase::betay[cctk_iteration='//iteration_string(1:nc)//']' - var_names(4) = 'admbase::betaz[cctk_iteration='//iteration_string(1:nc)//']' - var_names(5) = 'admbase::gxx[cctk_iteration='//iteration_string(1:nc)//']' - var_names(6) = 'admbase::gxy[cctk_iteration='//iteration_string(1:nc)//']' - var_names(7) = 'admbase::gxz[cctk_iteration='//iteration_string(1:nc)//']' - var_names(8) = 'admbase::gyy[cctk_iteration='//iteration_string(1:nc)//']' - var_names(9) = 'admbase::gyz[cctk_iteration='//iteration_string(1:nc)//']' - var_names(10) = 'admbase::gzz[cctk_iteration='//iteration_string(1:nc)//']' - var_names(11) = 'staticconformal::psi[cctk_iteration='//iteration_string(1:nc)//']' - - in_vars = ' ' - ntot = 0 - do i = 1, 11 - nc = len_trim(var_names(i)) - in_vars(ntot+1:ntot+1+nc+1) = var_names(i)(1:nc+1) - ntot = ntot + nc + 1 - end do - - print*,in_vars(1:len_trim(in_vars)),len_trim(in_vars) - - call IOUtil_RecoverVarsFromDatafiles ( res, cctkGH, in_files, in_vars ) - print*,res - current_iteration = current_iteration - 1 - - return -end subroutine EHFinder_ReadData diff --git a/src/EHFinder_SetMask.F90 b/src/EHFinder_SetMask.F90 index 566f654..eb13858 100644 --- a/src/EHFinder_SetMask.F90 +++ b/src/EHFinder_SetMask.F90 @@ -65,9 +65,8 @@ end subroutine EHFinder_MaskInit ! This routine will excise points and re-activate excised points -! after need. EHFinder_Setmask2 will then find the boundary of the +! as needed. EHFinder_Setmask2 will then find the boundary of the ! excsion region and make sure that the mask is correct there. - subroutine EHFinder_SetMask1(CCTK_ARGUMENTS) use EHFinder_mod @@ -218,9 +217,6 @@ subroutine EHFinder_SetMask1(CCTK_ARGUMENTS) call CCTK_WARN ( 0, 'Reduction of array fimax_loc failed' ) end if -! print*, imin_glob(1), imax_glob(1), fimin_glob(1), fimax_glob(1) -! print*, imin_glob(2), imax_glob(2), fimin_glob(2), fimax_glob(2) -! print*, imin_glob(3), imax_glob(3), fimin_glob(3), fimax_glob(3) end if ! Now check and see if any interior excised points need to be activated. @@ -307,11 +303,6 @@ subroutine EHFinder_SetMask1(CCTK_ARGUMENTS) end if end do -! print*,'New range' -! print*,imin_n(1),imax_n(1) -! print*,imin_n(2),imax_n(2) -! print*,imin_n(3),imax_n(3) - ! Use the new indices to actually activate points, taking care to ! activate points that are actually on the local grid. First do the ! faces of the rectangular box. @@ -322,11 +313,9 @@ subroutine EHFinder_SetMask1(CCTK_ARGUMENTS) j2 = min ( min ( imax_n(2), imax_glob(2) ) - cctk_lbnd(2), jyr ) k1 = max ( max ( imin_n(3), imin_glob(3) ) - cctk_lbnd(3), kzl ) k2 = min ( min ( imax_n(3), imax_glob(3) ) - cctk_lbnd(3), kzr ) -! print*,'x1 ',i1,j1,j2,k1,k2 if ( ( j1 .le. j2 ) .and. ( k1 .le. k2 ) ) then tm_mask(i1,j1:j2,k1:k2,l) = 0 ftmp(i1,j1:j2,k1:k2,l) = ftmp(i1+1,j1:j2,k1:k2,l) + delta -! print*,'done' end if end if end if @@ -337,11 +326,9 @@ subroutine EHFinder_SetMask1(CCTK_ARGUMENTS) j2 = min ( min ( imax_n(2), imax_glob(2) ) - cctk_lbnd(2), jyr ) k1 = max ( max ( imin_n(3), imin_glob(3) ) - cctk_lbnd(3), kzl ) k2 = min ( min ( imax_n(3), imax_glob(3) ) - cctk_lbnd(3), kzr ) -! print*,'x2 ',i2,j1,j2,k1,k2 if ( ( j1 .le. j2 ) .and. ( k1 .le. k2 ) ) then tm_mask(i2,j1:j2,k1:k2,l) = 0 ftmp(i2,j1:j2,k1:k2,l) = ftmp(i2-1,j1:j2,k1:k2,l) + delta -! print*,'done' end if end if end if @@ -352,11 +339,9 @@ subroutine EHFinder_SetMask1(CCTK_ARGUMENTS) i2 = min ( min ( imax_n(1), imax_glob(1) ) - cctk_lbnd(1), ixr ) k1 = max ( max ( imin_n(3), imin_glob(3) ) - cctk_lbnd(3), kzl ) k2 = min ( min ( imax_n(3), imax_glob(3) ) - cctk_lbnd(3), kzr ) -! print*,'y1 ',j1,i1,i2,k1,k2 if ( ( i1 .le. i2 ) .and. ( k1 .le. k2 ) ) then tm_mask(i1:i2,j1,k1:k2,l) = 0 ftmp(i1:i2,j1,k1:k2,l) = ftmp(i1:i2,j1+1,k1:k2,l) + delta -! print*,'done' end if end if end if @@ -367,11 +352,9 @@ subroutine EHFinder_SetMask1(CCTK_ARGUMENTS) i2 = min ( min ( imax_n(1), imax_glob(1) ) - cctk_lbnd(1), ixr ) k1 = max ( max ( imin_n(3), imin_glob(3) ) - cctk_lbnd(3), kzl ) k2 = min ( min ( imax_n(3), imax_glob(3) ) - cctk_lbnd(3), kzr ) -! print*,'y2 ',j2,i1,i2,k1,k2 if ( ( i1 .le. i2 ) .and. ( k1 .le. k2 ) ) then tm_mask(i1:i2,j2,k1:k2,l) = 0 ftmp(i1:i2,j2,k1:k2,l) = f(i1:i2,j2-1,k1:k2,l) + delta -! print*,'done' end if end if end if @@ -382,11 +365,9 @@ subroutine EHFinder_SetMask1(CCTK_ARGUMENTS) i2 = min ( min ( imax_n(1), imax_glob(1) ) - cctk_lbnd(1), ixr ) j1 = max ( max ( imin_n(2), imin_glob(2) ) - cctk_lbnd(2), jyl ) j2 = min ( min ( imax_n(2), imax_glob(2) ) - cctk_lbnd(2), jyr ) -! print*,'z1 ',k1,i1,i2,j1,j2 if ( ( i1 .le. i2 ) .and. ( j1 .le. j2 ) ) then tm_mask(i1:i2,j1:j2,k1,l) = 0 ftmp(i1:i2,j1:j2,k1,l) = ftmp(i1:i2,j1:j2,k1+1,l) + delta -! print*,'done' end if end if end if @@ -397,11 +378,9 @@ subroutine EHFinder_SetMask1(CCTK_ARGUMENTS) i2 = min ( min ( imax_n(1), imax_glob(1) ) - cctk_lbnd(1), ixr ) j1 = max ( max ( imin_n(2), imin_glob(2) ) - cctk_lbnd(2), jyl ) j2 = min ( min ( imax_n(2), imax_glob(2) ) - cctk_lbnd(2), jyr ) -! print*,'z2 ',k2,i1,i2,j1,j2 if ( ( i1 .le. i2 ) .and. ( j1 .le. j2 ) ) then tm_mask(i1:i2,j1:j2,k2,l) = 0 ftmp(i1:i2,j1:j2,k2,l) = ftmp(i1:i2,j1:j2,k2-1,l) + delta -! print*,'done' end if end if end if @@ -674,7 +653,6 @@ subroutine EHFinder_SetMask1(CCTK_ARGUMENTS) i2 = imax_n(1) - cctk_lbnd(1) j2 = imax_n(2) - cctk_lbnd(2) k2 = imax_n(3) - cctk_lbnd(3) -! print*,'Debug 3 : ',i2,j2,k2 if ( ( ixl .le. i2 ) .and. ( i2 .le. min(ixr,nx-1) ) .and. & ( jyl .le. j2 ) .and. ( j2 .le. min(jyr,ny-1) ) .and. & ( kzl .le. k2 ) .and. ( k2 .le. min(kzr,nz-1) ) ) then @@ -829,7 +807,6 @@ subroutine EHFinder_SetMask2(CCTK_ARGUMENTS) tm_mask(ixl:ixr,jyl:jyr,kzl:kzr,l) end where end if -! call CCTK_INFO('Mask modified') end do ! Indicate that it is not anymore the first time the mask is set. @@ -839,6 +816,10 @@ subroutine EHFinder_SetMask2(CCTK_ARGUMENTS) end subroutine EHFinder_SetMask2 +! This routine checks if any of the excision regions are to close together +! so that some points do not have enough active neigbouring points to be +! able to calculate derivatives. If this is the case these points are +! excised as well. subroutine EHFinder_SetMask3(CCTK_ARGUMENTS) use EHFinder_mod @@ -884,81 +865,76 @@ subroutine EHFinder_SetMask3(CCTK_ARGUMENTS) do j = jyl, jyr do i = ixl, ixr if ( eh_mask(i,j,k,l) .ge. 0 ) then - ! If we are on an excision boundary in the -x-direction... + +! If we are on an excision boundary in the -x-direction... if ( btest(eh_mask(i,j,k,l),0) ) then - ! If any of the two closest neighbours in the +x-direction is - ! excised we have to excise this point + +! If any of the two closest neighbours in the +x-direction is +! excised we have to excise this point if ( ( eh_mask(i+1,j,k,l) .eq. -1 ) .or. & ( eh_mask(i+2,j,k,l) .eq. -1 ) ) then tm_mask(i,j,k,l) = -1 mask_modified = .true. - print*, i, j, k, 'x-' - print*, eh_mask(i:i+2,j,k,l) end if end if - ! If we are on an excision boundary in the +x-direction... +! If we are on an excision boundary in the +x-direction... if ( btest(eh_mask(i,j,k,l),1) ) then - ! If any of the two closest neighbours in the -x-direction is - ! excised we have to excise this point + +! If any of the two closest neighbours in the -x-direction is +! excised we have to excise this point if ( ( eh_mask(i-1,j,k,l) .eq. -1 ) .or. & ( eh_mask(i-2,j,k,l) .eq. -1 ) ) then tm_mask(i,j,k,l) = -1 mask_modified = .true. - print*, i, j, k, 'x+' - print*, eh_mask(i:i-2,j,k,l) end if end if - ! If we are on an excision boundary in the -y-direction... +! If we are on an excision boundary in the -y-direction... if ( btest(eh_mask(i,j,k,l),2) ) then - ! If any of the two closest neighbours in the +y-direction is - ! excised we have to excise this point + +! If any of the two closest neighbours in the +y-direction is +! excised we have to excise this point if ( ( eh_mask(i,j+1,k,l) .eq. -1 ) .or. & ( eh_mask(i,j+2,k,l) .eq. -1 ) ) then tm_mask(i,j,k,l) = -1 mask_modified = .true. - print*, i, j, k, 'y-' - print*, eh_mask(i,j:j+2,k,l) end if end if - ! If we are on an excision boundary in the +y-direction... +! If we are on an excision boundary in the +y-direction... if ( btest(eh_mask(i,j,k,l),3) ) then - ! If any of the two closest neighbours in the -y-direction is - ! excised we have to excise this point + +! If any of the two closest neighbours in the -y-direction is +! excised we have to excise this point if ( ( eh_mask(i,j-1,k,l) .eq. -1 ) .or. & ( eh_mask(i,j-2,k,l) .eq. -1 ) ) then tm_mask(i,j,k,l) = -1 mask_modified = .true. - print*, i, j, k, 'y+' - print*, eh_mask(i,j:j-2,k,l) end if end if - ! If we are on an excision boundary in the -z-direction... +! If we are on an excision boundary in the -z-direction... if ( btest(eh_mask(i,j,k,l),4) ) then - ! If any of the two closest neighbours in the +z-direction is - ! excised we have to excise this point + +! If any of the two closest neighbours in the +z-direction is +! excised we have to excise this point if ( ( eh_mask(i,j,k+1,l) .eq. -1 ) .or. & ( eh_mask(i,j,k+2,l) .eq. -1 ) ) then tm_mask(i,j,k,l) = -1 mask_modified = .true. - print*, i, j, k, 'z-' - print*, eh_mask(i,j,k:k+2,l) end if end if - ! If we are on an excision boundary in the +z-direction... +! If we are on an excision boundary in the +z-direction... if ( btest(eh_mask(i,j,k,l),5) ) then - ! If any of the two closest neighbours in the -z-direction is - ! excised we have to excise this point + +! If any of the two closest neighbours in the -z-direction is +! excised we have to excise this point if ( ( eh_mask(i,j,k-1,l) .eq. -1 ) .or. & ( eh_mask(i,j,k-2,l) .eq. -1 ) ) then tm_mask(i,j,k,l) = -1 mask_modified = .true. - print*, i, j, k, 'z+' - print*, eh_mask(i,j,k:k-2,l) end if end if end if @@ -968,7 +944,7 @@ subroutine EHFinder_SetMask3(CCTK_ARGUMENTS) end do if ( mask_modified ) then - call CCTK_WARN ( 1, "Mask modified because it was not convex" ) + call CCTK_WARN ( 1, 'Mask modified because it was not convex' ) end if end if end do diff --git a/src/EHFinder_SetSym.F90 b/src/EHFinder_SetSym.F90 index b04672d..7811406 100644 --- a/src/EHFinder_SetSym.F90 +++ b/src/EHFinder_SetSym.F90 @@ -1,4 +1,4 @@ -! Initialisation of the level set function +! Registration of symmetries for the necessary grid functions. ! $Header$ #include "cctk.h" @@ -15,33 +15,28 @@ subroutine EHFinder_SetSym(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS CCTK_INT, dimension(3) :: sym - CCTK_INT :: l - character(len=14) :: fname +! All grid functions have even symmetries in all directions. sym = 1 - 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 +! Set up symmetries for the level set function. + call SetCartSymGN ( 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 + +! Set up symmetries for the mask function. + call SetCartSymGN( ierr, cctkGH, sym, 'ehfinder::eh_mask' ) + if ( ierr .gt. 0 ) then + call CCTK_WARN(1,'Failed to register symmetry for eh_mask') + end if +! Set up symmetries for the surface index function. 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 -! call SetCartSymVN(ierr,cctkGH,sym,'ehfinder::rep_mask') -! if ( ierr .gt. 0 ) then -! call CCTK_WARN(1,'Failed to register symmetry for the level reparametrization mask') -! end if - return end subroutine EHFinder_SetSym @@ -56,28 +51,22 @@ subroutine EHFinder_ApplySymAll(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS - 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 + character(len=80) :: warn_message -# 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,:,:,:) - 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,:,:) +! Apply symmetry for the level set function. + call CartSymGN ( ierr, cctkGH, 'ehfinder::f' ) + if ( ierr .gt. 0 ) then + warn_message = 'Failed to perform symmetry operation on the level ' + warn_message = warn_message//'set function' + call CCTK_WARN( 1, warn_message ) 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,:) + +! Apply symmetry for the mask function. + call CartSymGN ( ierr, cctkGH, 'ehfinder::eh_mask' ) + if ( ierr .gt. 0 ) then + call CCTK_WARN( 1, 'Failed to perform symmetry operation on the mask') end if + return end subroutine EHFinder_ApplySymAll @@ -92,19 +81,21 @@ subroutine EHFinder_ApplySymF(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS - CCTK_INT :: l - character(len=14) :: fname + character(len=80) :: warn_message - do l = 0, eh_number_level_sets - 1 - write(fname,'(a12,i1,a1)') 'ehfinder::f[', l, ']' - call CartSymVN(ierr,cctkGH,fname) - end do +! Apply symmetry for the level set function. + call CartSymGN ( ierr, cctkGH, 'ehfinder::f' ) + if ( ierr .gt. 0 ) then + warn_message = 'Failed to perform symmetry operation on the level ' + warn_message = warn_message//'set function' + call CCTK_WARN( 1, warn_message ) + end if return end subroutine EHFinder_ApplySymF -subroutine EHFinder_ApplySymRep(CCTK_ARGUMENTS) +subroutine EHFinder_ApplySymMask(CCTK_ARGUMENTS) use EHFinder_mod @@ -114,55 +105,17 @@ subroutine EHFinder_ApplySymRep(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS -# 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 - return -end subroutine EHFinder_ApplySymRep - - -subroutine EHFinder_ApplySymFRep(CCTK_ARGUMENTS) - - use EHFinder_mod - - implicit none - - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_PARAMETERS - DECLARE_CCTK_FUNCTIONS +! Apply symmetry for the mask function. + call CartSymGN ( ierr, cctkGH, 'ehfinder::eh_mask' ) + if ( ierr .gt. 0 ) then + call CCTK_WARN( 1, 'Failed to perform symmetry operation on the mask') + end if - CCTK_INT :: l - character(len=14) :: fname - -# include "include/physical_part.h" - - 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 return -end subroutine EHFinder_ApplySymFRep +end subroutine EHFinder_ApplySymMask -subroutine EHFinder_ApplySymMask(CCTK_ARGUMENTS) +subroutine EHFinder_ApplySymSC(CCTK_ARGUMENTS) use EHFinder_mod @@ -172,16 +125,14 @@ subroutine EHFinder_ApplySymMask(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS -# include "include/physical_part.h" + character(len=80) :: warn_message - if ( symx ) then - 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,:,:) - end if - if ( symz ) then - eh_mask(:,:,ngz:1:-1,:) = eh_mask(:,:,ngz+1:ngz+ngz,:) +! Apply the symmetry for the surface index function. + call CartSymGN ( ierr, cctkGH, 'ehfinder::surface_index' ) + if ( ierr .gt. 0 ) then + warn_message = 'Failed to perform symmetry operation on the surface index' + call CCTK_WARN( 1, warn_message ) end if + return -end subroutine EHFinder_ApplySymMask +end subroutine EHFinder_ApplySymSC diff --git a/src/EHFinder_Sources.F90 b/src/EHFinder_Sources.F90 index 65e0488..69a5464 100644 --- a/src/EHFinder_Sources.F90 +++ b/src/EHFinder_Sources.F90 @@ -27,37 +27,44 @@ subroutine EHFinder_Sources(CCTK_ARGUMENTS) CCTK_REAL :: alminus, alplus, blminus, blplus, clminus, clplus CCTK_REAL :: arminus, arplus, brminus, brplus, crminus, crplus -! calculate 1/(2*delta) in each direction - #include "include/physical_part.h" +! calculate 1/(2*delta) in each direction idx = half / cctk_delta_space(1) idy = half / cctk_delta_space(2) idz = half / cctk_delta_space(3) mdelta = maxval ( cctk_delta_space ) +! Set the sign depending on the surface direction. if ( CCTK_EQUALS ( surface_direction, 'outward' ) ) ssign = one if ( CCTK_EQUALS ( surface_direction, 'inward' ) ) ssign = -one - if ( CCTK_EQUALS ( upwind_type, 'intrinsic' ) ) then - do l = 1, eh_number_level_sets + do l = 1, eh_number_level_sets + + do k = kzl, kzr + do j = jyl, jyr + do i = ixl, ixr +! Calculate the inverse of the 3-metric +# include "include/metric.h" + end do + end do + end do + +! Calculate the derivatives of the level set function using the intrinsic +! scheme. Note, this should never be used and may disappear in later versions. + if ( CCTK_EQUALS ( upwind_type, 'intrinsic' ) ) then 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 #include "include/upwind_second2.h" - end if end do end do end do - end do - end if + end if - if ( CCTK_EQUALS ( upwind_type, 'shift' ) ) then - do l = 1, eh_number_level_sets +! Calculate the derivatives of the level set function using shift upwinding. + if ( CCTK_EQUALS ( upwind_type, 'shift' ) ) then do k = kzl, kzr do j = jyl, jyr do i = ixl, ixr @@ -65,119 +72,89 @@ subroutine EHFinder_Sources(CCTK_ARGUMENTS) end do end do end do - end do - end if + end if - if ( CCTK_EQUALS ( upwind_type, 'characteristic' ) ) then - do l = 1, eh_number_level_sets +! If the three metric is the static conformal metric we convert the inverse +! three metric to the physical inverse three metric by multiplying with +! psi^(-4). + if ( CCTK_EQUALS ( metric_type, 'static conformal' ) ) then do k = kzl, kzr - do j = jyl, jyr + do j = jyl,jyr do i = ixl, ixr -#include "include/metric.h" -#include "include/upwind_characteristic_second2.h" + if ( eh_mask(i,j,k,l) .ge. 0 ) then + psito4 = psi(i,j,k)**(-4) + g3xx(i,j,k) = g3xx(i,j,k) * psito4 + g3xy(i,j,k) = g3xy(i,j,k) * psito4 + g3xz(i,j,k) = g3xz(i,j,k) * psito4 + g3yy(i,j,k) = g3yy(i,j,k) * psito4 + g3yz(i,j,k) = g3yz(i,j,k) * psito4 + g3zz(i,j,k) = g3zz(i,j,k) * psito4 + end if end do end do end do - end do - end if + end if - if ( CCTK_EQUALS ( metric_type, 'physical' ) ) then - do l = 1, eh_number_level_sets +! Calculate the derivatives of the level set using characteristic upwinding. + if ( CCTK_EQUALS ( upwind_type, 'characteristic' ) ) then 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 - - 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 - sf(i,j,k,l) = zero - end if + +! We use centered derivatives to figure out which direction +! to upwind in. +#include "include/upwind_characteristic_second2.h" end do end do - end do - end do - end if + end do + end if - if ( CCTK_EQUALS ( metric_type, 'static conformal' ) ) then - 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 - - 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,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 + do k = kzl, kzr + do j = jyl, jyr + do i = ixl, ixr + +! If the current point is active ... + if ( eh_mask(i,j,k,l) .ge. 0 ) then + +! Square the lapse. + alp2 = alp(i,j,k)**2 + +! Calculate beta^i df_i. + 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) + +! Calculate gamma^ij df_i df_j. + tmp2 = g3xx(i,j,k) * dfx(i,j,k,l)**2 + & + g3yy(i,j,k) * dfy(i,j,k,l)**2 + & + g3zz(i,j,k) * dfz(i,j,k,l)**2 + & + two * ( g3xy(i,j,k) * dfx(i,j,k,l) * dfy(i,j,k,l) + & + g3xz(i,j,k) * dfx(i,j,k,l) * dfz(i,j,k,l) + & + g3yz(i,j,k) * dfy(i,j,k,l) * dfz(i,j,k,l) ) + +! If the metric is positive definite ... + if ( tmp2 .ge. zero ) then + +! Calculate the right hand side. + sf(i,j,k,l) = tmp1 - ssign * sqrt ( alp2 * tmp2 ) + +! If the lapse is negative we change the sign of the right hand +! side function. This is done to be able to handle for example +! Schwarzschild in isotropic coordinates with the isotropic +! lapse. + sf(i,j,k,l) = sf(i,j,k,l) * sign ( one, alp(i,j,k) ) else - sf(i,j,k,l) = zero + +! Otherwise print a level 0 warning. + call CCTK_WARN ( 0, '3-metric not positive definite: Stopping' ) end if - end do + else + sf(i,j,k,l) = zero + end if end do - end do - end do - end if - + end do + end do + end do + return end subroutine EHFinder_Sources diff --git a/src/EHFinder_mod.F90 b/src/EHFinder_mod.F90 index 9c731a0..3fe10a5 100644 --- a/src/EHFinder_mod.F90 +++ b/src/EHFinder_mod.F90 @@ -6,33 +6,21 @@ module EHFinder_mod use EHFinder_Constants - CCTK_INT, dimension(0:5), parameter :: ix = (/ -1, 1, 0, 0, 0, 0 /) - CCTK_INT, dimension(0:5), parameter :: jy = (/ 0, 0, -1, 1, 0, 0 /) - CCTK_INT, dimension(0:5), parameter :: kz = (/ 0, 0, 0, 0, -1, 1 /) CCTK_INT, dimension(0:5), parameter :: ll = (/ 1, 2, 4, 8, 16, 32 /) CCTK_INT :: ixl, ixr, jyl, jyr, kzl, kzr -! CCTK_INT, parameter :: ixl = 1 -! CCTK_INT, parameter :: ixr = 2 -! CCTK_INT, parameter :: iyl = 4 -! CCTK_INT, parameter :: iyr = 8 -! CCTK_INT, parameter :: izl = 16 -! CCTK_INT, parameter :: izr = 32 CCTK_REAL :: hfac - CCTK_INT :: nx, ny, nz, ngx, ngy, ngz, ierr + CCTK_INT :: nx, ny, nz, ngx, ngy, ngz CCTK_INT :: max_handle, min_handle, sum_handle - CCTK_INT :: ncalls, rep_count, rep_total + CCTK_INT :: ncalls, ierr CCTK_REAL :: delta - CCTK_REAL :: h, hlocal + CCTK_REAL :: h CCTK_REAL :: ex_value - CCTK_REAL :: rdel = 0.8d0 - CCTK_REAL :: fmin_bound CCTK_INT, dimension(3) :: imin_glob, imax_glob CCTK_REAL, dimension(3) :: fimin_glob, fimax_glob - CCTK_INT :: current_iteration CCTK_INT :: githeta, gjphi logical :: mask_first = .true. logical :: symx, symy, symz logical :: s_symx, s_symy, s_symz - logical, dimension(:), allocatable :: reparam_this_level_set, & + logical, dimension(:), allocatable :: re_init_this_level_set, & re_initialize_undone end module EHFinder_mod diff --git a/src/MoL_Init.F90 b/src/MoL_Init.F90 index 4683c09..58a5611 100644 --- a/src/MoL_Init.F90 +++ b/src/MoL_Init.F90 @@ -10,23 +10,19 @@ subroutine EHFinder_MoLRegister(CCTK_ARGUMENTS) implicit none CCTK_INT :: ierr, ierr_cum, varindex, rhsindex - CCTK_INT :: MoLRegisterEvolved, MoLRegisterEvolvedGroup + CCTK_INT :: MoLRegisterEvolvedGroup!, MoLRegisterEvolved CCTK_INT :: l - character(len=15) :: vname - DECLARE_CCTK_PARAMETERS DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS ierr_cum = 0 - 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 + + call CCTK_GroupIndex ( varindex, 'ehfinder::f' ) + call CCTK_GroupIndex ( rhsindex, 'ehfinder::sf' ) + + ierr_cum = ierr_cum + MolRegisterEvolvedGroup ( varindex, rhsindex ) if ( evolve_generators .gt. 0 ) then @@ -47,37 +43,6 @@ subroutine EHFinder_MoLRegister(CCTK_ARGUMENTS) end if -! call CCTK_VarIndex(varindex, 'admbase::gxx') -! call MoL_RegisterPrimitive(ierr, varindex) -! ierr_cum = ierr_cum + ierr -! call CCTK_VarIndex(varindex, 'admbase::gxy') -! call MoL_RegisterPrimitive(ierr, varindex) -! ierr_cum = ierr_cum + ierr -! call CCTK_VarIndex(varindex, 'admbase::gxz') -! call MoL_RegisterPrimitive(ierr, varindex) -! ierr_cum = ierr_cum + ierr -! call CCTK_VarIndex(varindex, 'admbase::gyy') -! call MoL_RegisterPrimitive(ierr, varindex) -! ierr_cum = ierr_cum + ierr -! call CCTK_VarIndex(varindex, 'admbase::gyz') -! call MoL_RegisterPrimitive(ierr, varindex) -! ierr_cum = ierr_cum + ierr -! call CCTK_VarIndex(varindex, 'admbase::gzz') -! call MoL_RegisterPrimitive(ierr, varindex) -! ierr_cum = ierr_cum + ierr -! call CCTK_VarIndex(varindex, 'admbase::alp') -! call MoL_RegisterPrimitive(ierr, varindex) -! ierr_cum = ierr_cum + ierr -! call CCTK_VarIndex(varindex, 'admbase::betax') -! call MoL_RegisterPrimitive(ierr, varindex) -! ierr_cum = ierr_cum + ierr -! call CCTK_VarIndex(varindex, 'admbase::betay') -! call MoL_RegisterPrimitive(ierr, varindex) -! ierr_cum = ierr_cum + ierr -! call CCTK_VarIndex(varindex, 'admbase::betaz') -! call MoL_RegisterPrimitive(ierr, varindex) -! ierr_cum = ierr_cum + ierr - print*,'MoLRegister ', ierr_cum if ( ierr_cum .gt. 0 ) then call CCTK_WARN(0,'Problems registering variables with MoL') end if -- cgit v1.2.3