From 8425a5848ba1716273567a6fbc177ec1c6fcf763 Mon Sep 17 00:00:00 2001 From: diener Date: Tue, 3 Jun 2003 11:42:41 +0000 Subject: Added support for a conformal metric. Added support for the new interpolator parameters. Cleaned up the code a bit. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/EHFinder/trunk@113 2a26948c-0e4f-0410-aee8-f1d3e353619c --- src/EHFinder_Generator_Sources.F90 | 218 +++++++++++++++++++++++++------------ 1 file changed, 150 insertions(+), 68 deletions(-) diff --git a/src/EHFinder_Generator_Sources.F90 b/src/EHFinder_Generator_Sources.F90 index 924588f..95cb22e 100644 --- a/src/EHFinder_Generator_Sources.F90 +++ b/src/EHFinder_Generator_Sources.F90 @@ -17,6 +17,11 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS) CCTK_INT :: i CCTK_INT :: interp_handle, table_handle, status, coord_system_handle + + character(len=200) :: gen_interp + CCTK_INT :: gen_interp_len + character(len=7) :: gen_order + CCTK_INT, dimension(1) :: lsh CCTK_POINTER, dimension(3) :: interp_coords CCTK_POINTER, dimension(14) :: out_arrays @@ -28,78 +33,102 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS) 0, 0, 0, 0, 0, & 1, 2, 3, 0 /) CCTK_INT, dimension(14) :: out_types - CCTK_REAL :: alp2, dfux, dfuy, dfuz, factor + CCTK_REAL :: alp2, psi4, dfux, dfuy, dfuz, factor CCTK_REAL :: idetg, guxx, guxy, guxz, guyy, guyz, guzz out_types = CCTK_VARIABLE_REAL - if ( CCTK_EQUALS ( metric_type, 'physical' ) ) then - call CCTK_GrouplshGN ( status, cctkGH, 1, lsh, "ehfinder::generators" ) - if ( status .lt. 0 ) then - call CCTK_WARN ( 0, "cannot get local size for surface arrays" ) - end if -! print*,lsh + ! Convert the generator_interpolator string parameter to a Fortran string. + call CCTK_FortranString ( gen_interp_len, generator_interpolator, & + gen_interp ) - call CCTK_InterpHandle ( interp_handle, & - "Lagrange polynomial interpolation" ) - if ( interp_handle .lt. 0 ) then - call CCTK_WARN( 0, "Cannot get handle for interpolation. Forgot to activate an implementation providing interpolation operators??" ) - end if - ! For now order=2 is hard wired. - call Util_TableCreateFromString ( table_handle, "order=3" ) - if ( table_handle .lt. 0 ) then - call CCTK_WARN( 0, "Cannot create parameter table for interpolator" ) - end if + ! Get the corresponding interpolator handle. + call CCTK_InterpHandle ( interp_handle, gen_interp ) - ! 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 ??" ) - endif + if ( interp_handle .lt. 0 ) then + call CCTK_WARN( 0, "Cannot get handle for interpolation. Forgot to activate an implementation providing interpolation operators??" ) + 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. + write(gen_order,'(a6,i1)') 'order=',generator_interpolation_order + + ! 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" ) + end if - interp_coords(1) = CCTK_PointerTo(xg) - interp_coords(2) = CCTK_PointerTo(yg) - interp_coords(3) = CCTK_PointerTo(zg) - - out_arrays(1) = CCTK_PointerTo(alpg) - out_arrays(2) = CCTK_PointerTo(betaxg) - out_arrays(3) = CCTK_PointerTo(betayg) - out_arrays(4) = CCTK_PointerTo(betazg) - out_arrays(5) = CCTK_PointerTo(gxxg) - out_arrays(6) = CCTK_PointerTo(gxyg) - out_arrays(7) = CCTK_PointerTo(gxzg) - out_arrays(8) = CCTK_PointerTo(gyyg) - out_arrays(9) = CCTK_PointerTo(gyzg) - out_arrays(10) = CCTK_PointerTo(gzzg) - out_arrays(11) = CCTK_PointerTo(dfxg) - out_arrays(12) = CCTK_PointerTo(dfyg) - out_arrays(13) = CCTK_PointerTo(dfzg) - - call CCTK_VarIndex ( in_arrays(1), "admbase::alp" ) - call CCTK_VarIndex ( in_arrays(2), "admbase::betax" ) - call CCTK_VarIndex ( in_arrays(3), "admbase::betay" ) - call CCTK_VarIndex ( in_arrays(4), "admbase::betaz" ) - call CCTK_VarIndex ( in_arrays(5), "admbase::gxx" ) - call CCTK_VarIndex ( in_arrays(6), "admbase::gxy" ) - call CCTK_VarIndex ( in_arrays(7), "admbase::gxz" ) - call CCTK_VarIndex ( in_arrays(8), "admbase::gyy" ) - call CCTK_VarIndex ( in_arrays(9), "admbase::gyz" ) - call CCTK_VarIndex ( in_arrays(10), "admbase::gzz" ) - call CCTK_VarIndex ( in_arrays(11), "ehfinder::f" ) + ! 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 ??" ) + endif + + ! Find out how many interpolation points are located on this processor. + call CCTK_GrouplshGN ( status, cctkGH, 1, lsh, "ehfinder::generators" ) + if ( status .lt. 0 ) then + call CCTK_WARN ( 0, "cannot get local size for surface arrays" ) + end if + + ! Set the pointers to the points to be interpolated to. + interp_coords(1) = CCTK_PointerTo(xg) + interp_coords(2) = CCTK_PointerTo(yg) + interp_coords(3) = CCTK_PointerTo(zg) + + ! Set the pointers to the output arrays. + out_arrays(1) = CCTK_PointerTo(alpg) + out_arrays(2) = CCTK_PointerTo(betaxg) + out_arrays(3) = CCTK_PointerTo(betayg) + out_arrays(4) = CCTK_PointerTo(betazg) + out_arrays(5) = CCTK_PointerTo(gxxg) + out_arrays(6) = CCTK_PointerTo(gxyg) + out_arrays(7) = CCTK_PointerTo(gxzg) + out_arrays(8) = CCTK_PointerTo(gyyg) + out_arrays(9) = CCTK_PointerTo(gyzg) + out_arrays(10) = CCTK_PointerTo(gzzg) + out_arrays(11) = CCTK_PointerTo(dfxg) + out_arrays(12) = CCTK_PointerTo(dfyg) + out_arrays(13) = CCTK_PointerTo(dfzg) + out_arrays(14) = CCTK_PointerTo(psig) + + ! Set the indices to the input grid functions. + call CCTK_VarIndex ( in_arrays(1), "admbase::alp" ) + call CCTK_VarIndex ( in_arrays(2), "admbase::betax" ) + call CCTK_VarIndex ( in_arrays(3), "admbase::betay" ) + call CCTK_VarIndex ( in_arrays(4), "admbase::betaz" ) + call CCTK_VarIndex ( in_arrays(5), "admbase::gxx" ) + call CCTK_VarIndex ( in_arrays(6), "admbase::gxy" ) + call CCTK_VarIndex ( in_arrays(7), "admbase::gxz" ) + call CCTK_VarIndex ( in_arrays(8), "admbase::gyy" ) + call CCTK_VarIndex ( in_arrays(9), "admbase::gyz" ) + call CCTK_VarIndex ( in_arrays(10), "admbase::gzz" ) + call CCTK_VarIndex ( in_arrays(11), "ehfinder::f" ) + call CCTK_VarIndex ( in_arrays(12), "staticconformal::psi" ) + + ! Check the metric type. At present physical and static_conformal are + ! supported. + if ( CCTK_EQUALS ( metric_type, 'physical' ) ) then + ! Set the operand indices table entry, corresponding + ! to interpolation of admbase::alp (1), admbase::shift (3), + ! admbase::metric(6) and the first derivatives of ehfinder::f (3) for + ! a total of 13 output arrays. call Util_TableSetIntArray ( status, table_handle, 13, & op_indices(1:13), "operand_indices" ) if ( status .lt. 0 ) then call CCTK_WARN ( 0, "Cannot set operand indices array in parameter table" ) endif + ! Set the corresponding table entry for the operation codes. call Util_TableSetIntArray ( status, table_handle, 13, & op_codes(1:13), "operation_codes" ) if ( status .lt. 0 ) then call CCTK_WARN ( 0, "Cannot set operation codes array in parameter table" ) endif - + ! Call the interpolator. call CCTK_InterpGridArrays ( status, cctkGH, 3, interp_handle, & table_handle, coord_system_handle, & lsh(1), CCTK_VARIABLE_REAL, & @@ -109,9 +138,15 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS) if ( status .lt. 0 ) then call CCTK_INFO ( 'Interpolation failed.' ) end if + + ! For each point on this processor calculate the right hand side of the + ! characteristic evolution equation. do i = 1, lsh(1) + + ! calculate the square of the lapse. alp2 = alpg(i)**2 + ! Calculate the inverse of the 3-metric. guxx = gyyg(i) * gzzg(i) - gyzg(i)**2 guxy = gxzg(i) * gyzg(i) - gxyg(i) * gzzg(i) guxz = gxyg(i) * gyzg(i) - gxzg(i) * gyyg(i) @@ -126,6 +161,69 @@ 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. + dfux = guxx * dfxg(i) + guxy * dfyg(i) + guxz * dfzg(i) + dfuy = guxy * dfxg(i) + guyy * dfyg(i) + guyz * dfzg(i) + dfuz = guxz * dfxg(i) + guyz * dfyg(i) + guzz * dfzg(i) + + ! Calculate the overall multiplication factor. + factor = alp2 / sqrt ( alp2 * ( dfux * dfxg(i) + & + dfuy * dfyg(i) + & + dfuz * dfzg(i) ) ) + + ! Finally obtain dx^i/dt. + dxg(i) = - betaxg(i) + factor * dfux + dyg(i) = - betayg(i) + factor * dfuy + dzg(i) = - betazg(i) + factor * dfuz + end do + else if ( CCTK_EQUALS ( metric_type, 'static conformal' ) ) then + call Util_TableSetIntArray ( status, table_handle, 14, & + op_indices, "operand_indices" ) + if ( status .lt. 0 ) then + call CCTK_WARN ( 0, "Cannot set operand indices array in parameter table" ) + endif + + call Util_TableSetIntArray ( status, table_handle, 14, & + op_codes, "operation_codes" ) + if ( status .lt. 0 ) then + call CCTK_WARN ( 0, "Cannot set operation codes array in parameter table" ) + endif + + call CCTK_InterpGridArrays ( status, cctkGH, 3, interp_handle, & + table_handle, coord_system_handle, & + lsh(1), CCTK_VARIABLE_REAL, & + interp_coords, 12, in_arrays, & + 14, out_types, out_arrays ) + + if ( status .lt. 0 ) then + call CCTK_INFO ( 'Interpolation failed.' ) + end if + + do i = 1, lsh(1) + alp2 = alpg(i)**2 + +! The inverse of psi^4 + psi4 = one / psig(i)**4 + + guxx = gyyg(i) * gzzg(i) - gyzg(i)**2 + guxy = gxzg(i) * gyzg(i) - gxyg(i) * gzzg(i) + guxz = gxyg(i) * gyzg(i) - gxzg(i) * gyyg(i) + +! The determinant divided by psi^4. + idetg = psi4 / ( gxxg(i) * guxx + & + gxyg(i) * guxy + & + gxzg(i) * guxz ) + +! The inverse metric. Since the determinant is already divided +! by psi^4, this gives the inverse of the physical metric. + guxx = idetg * guxx + guxy = idetg * guxy + guxz = idetg * guxz + + guyy = ( gxxg(i) * gzzg(i) - gxzg(i)**2 ) * idetg + guyz = ( gxyg(i) * gxzg(i) - gxxg(i) * gyzg(i) ) * idetg + guzz = ( gxxg(i) * gyyg(i) - gxyg(i)**2 ) * idetg + dfux = guxx * dfxg(i) + guxy * dfyg(i) + guxz * dfzg(i) dfuy = guxy * dfxg(i) + guyy * dfyg(i) + guyz * dfzg(i) dfuz = guxz * dfxg(i) + guyz * dfyg(i) + guzz * dfzg(i) @@ -135,22 +233,6 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS) dxg(i) = - betaxg(i) + factor * dfux dyg(i) = - betayg(i) + factor * dfuy dzg(i) = - betazg(i) + factor * dfuz -! print*, alpg -! print* -! print*, betaxg, betayg, betazg -! print* -! print*, gxxg, gxyg, gxzg, gyyg, gyzg, gzzg -! print* -! print*, dfxg, dfyg, dfzg -! print* -! print*, alp2 -! print* -! print*, dfux, dfuy, dfuz -! print* -! print*, factor -! print* -! print*, dxg(i), dyg(i), dzg(i) -! stop end do end if return -- cgit v1.2.3