diff options
Diffstat (limited to 'src/EHFinder_Generator_Sources.F90')
| -rw-r--r-- | src/EHFinder_Generator_Sources.F90 | 118 |
1 files changed, 65 insertions, 53 deletions
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 ) |