From e74d0febe4c8a26c5a1decd92e9c4c637990aaea Mon Sep 17 00:00:00 2001 From: diener Date: Fri, 13 Feb 2004 13:40:01 +0000 Subject: First attempt at using 2d grid arrays to hold the generator information. Something fishy is going on. Don't use yet. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/EHFinder/trunk@165 2a26948c-0e4f-0410-aee8-f1d3e353619c --- interface.ccl | 32 ++++ param.ccl | 11 ++ schedule.ccl | 71 +++++--- src/EHFinder_Generator_Sources.F90 | 345 ++++++++++++++++++++++++++++++++++++- src/EHFinder_Init.F90 | 108 ++++++++---- src/MoL_Init.F90 | 39 ++++- 6 files changed, 535 insertions(+), 71 deletions(-) diff --git a/interface.ccl b/interface.ccl index f52afdc..9d887be 100644 --- a/interface.ccl +++ b/interface.ccl @@ -186,3 +186,35 @@ CCTK_REAL generator_gf TYPE=GF TIMELEVELS=1 { xgf, ygf, zgf } "Temporary grid function used in calculating the right hand side of the generator evolution equation" + +# The following is for a 2-d distribution of generators on the surface. + +# The X-position of the generators of the event horizons. +CCTK_REAL xg2[eh_number_level_sets] TYPE=ARRAY DIM=2 TIMELEVELS=2 SIZE=number_of_generators_theta,number_of_generators_phi GHOSTSIZE=0,0 DISTRIB=DEFAULT + +# The Y-position of the generators of the event horizons. +CCTK_REAL yg2[eh_number_level_sets] TYPE=ARRAY DIM=2 TIMELEVELS=2 SIZE=number_of_generators_theta,number_of_generators_phi GHOSTSIZE=0,0 DISTRIB=DEFAULT + +# The Z-position of the generators of the event horizons. +CCTK_REAL zg2[eh_number_level_sets] TYPE=ARRAY DIM=2 TIMELEVELS=2 SIZE=number_of_generators_theta,number_of_generators_phi GHOSTSIZE=0,0 DISTRIB=DEFAULT + +# The right hand side of the X-position of the generators. +CCTK_REAL dxg2[eh_number_level_sets] TYPE=ARRAY DIM=2 TIMELEVELS=1 SIZE=number_of_generators_theta,number_of_generators_phi GHOSTSIZE=0,0 DISTRIB=DEFAULT + +# The right hand side of the Y-position of the generators. +CCTK_REAL dyg2[eh_number_level_sets] TYPE=ARRAY DIM=2 TIMELEVELS=1 SIZE=number_of_generators_theta,number_of_generators_phi GHOSTSIZE=0,0 DISTRIB=DEFAULT + +# The right hand side of the Z-position of the generators. +CCTK_REAL dzg2[eh_number_level_sets] TYPE=ARRAY DIM=2 TIMELEVELS=1 SIZE=number_of_generators_theta,number_of_generators_phi GHOSTSIZE=0,0 DISTRIB=DEFAULT + + +CCTK_REAL generator_arrays2 TYPE=ARRAY DIM=2 TIMELEVELS=1 SIZE=number_of_generators_theta,number_of_generators_phi GHOSTSIZE=0,0 DISTRIB=DEFAULT +{ + alpg2, betaxg2, betayg2, betazg2, gxxg2, gxyg2, gxzg2, gyyg2, gyzg2, gzzg2, dfxg2, dfyg2, dfzg2, psig2 +} "Arrays to hold the interpolated metric, gauge and level set data" + +CCTK_REAL generator_gf2 TYPE=GF TIMELEVELS=1 +{ + xgf2, ygf2, zgf2 +} "Temporary grid function used in calculating the right hand side of the generator evolution equation" + diff --git a/param.ccl b/param.ccl index b08888f..0007865 100644 --- a/param.ccl +++ b/param.ccl @@ -215,9 +215,20 @@ CCTK_INT number_of_generators "How many generators should be evolved" 1:* :: "Postive please" } 1 +CCTK_INT number_of_generators_theta "How many generators in the theta direction" +{ + 1:* :: "Positive please" +} 1 + +CCTK_INT number_of_generators_phi "How many generators in the phi direction" +{ + 1:* :: "Positive please" +} 1 + KEYWORD generator_distribution "What initial distribution should be used" { "line" :: "Put the generators on a line in the xz-plane" + "2D array" :: "Put the generators on a surface with spherical topology" } "line" STRING surface_interpolator "What interpolator should be used to locate the surface" diff --git a/schedule.ccl b/schedule.ccl index 3e9f6ab..f522f5c 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -19,8 +19,16 @@ else STORAGE: center_arrays if ( evolve_generators ) { - STORAGE: xg[2], yg[2], zg[2] - STORAGE: dxg, dyg, dzg + if (CCTK_Equals(generator_distribution,"line")) + { + STORAGE: xg[2], yg[2], zg[2] + STORAGE: dxg, dyg, dzg + } + if (CCTK_Equals(generator_distribution,"2D array")) + { + STORAGE: xg2[2], yg2[2], zg2[2] + STORAGE: dxg2, dyg2, dzg2 + } } } @@ -341,11 +349,22 @@ if (CCTK_Equals(mode,"normal")) { if (CCTK_Equals(generator_tracking_method,"interpolate_before")) { - schedule EHFinder_Generator_Sources in MoL_CalcRHS + if (CCTK_Equals(generator_distribution,"line")) { - LANG: Fortran - STORAGE: generator_arrays - } "Calculate the source terms for the generator evolution" + schedule EHFinder_Generator_Sources in MoL_CalcRHS + { + LANG: Fortran + STORAGE: generator_arrays + } "Calculate the source terms for the generator evolution" + } + if (CCTK_Equals(generator_distribution,"2D array")) + { + schedule EHFinder_Generator_Sources_2D in MoL_CalcRHS + { + LANG: Fortran + STORAGE: generator_arrays2 + } "Calculate the source terms for the 2D generator evolution" + } } if (CCTK_Equals(generator_tracking_method,"interpolate_after")) { @@ -353,7 +372,7 @@ if (CCTK_Equals(mode,"normal")) { LANG: Fortran STORAGE: generator_arrays - } "Calculate the source terms for the generator evolution" + } "Calculate the source terms2 for the generator evolution" } } @@ -533,23 +552,23 @@ if (CCTK_Equals(mode,"normal")) LANG: Fortran } "Apply symmetry boundaries" - if ( evolve_generators) - { - if (CCTK_Equals(generator_tracking_method,"interpolate_before")) - { - schedule EHFinder_Generator_Sources as EGS in MoL_CalcRHS - { - LANG: Fortran - STORAGE: generator_arrays - } "Calculate the source terms for the generator evolution" - } - if (CCTK_Equals(generator_tracking_method,"interpolate_after")) - { - schedule EHFinder_Generator_Sources2 as EGS2 in MoL_CalcRHS - { - LANG: Fortran - STORAGE: generator_arrays - } "Calculate the source terms for the generator evolution" - } - } +# if ( evolve_generators) +# { +# if (CCTK_Equals(generator_tracking_method,"interpolate_before")) +# { +# schedule EHFinder_Generator_Sources as EGS in MoL_CalcRHS +# { +# LANG: Fortran +# STORAGE: generator_arrays +# } "Calculate the source terms for the generator evolution" +# } +# if (CCTK_Equals(generator_tracking_method,"interpolate_after")) +# { +# schedule EHFinder_Generator_Sources2 as EGS2 in MoL_CalcRHS +# { +# LANG: Fortran +# STORAGE: generator_arrays +# } "Calculate the source terms for the generator evolution" +# } +# } } diff --git a/src/EHFinder_Generator_Sources.F90 b/src/EHFinder_Generator_Sources.F90 index 984d80f..3dee327 100644 --- a/src/EHFinder_Generator_Sources.F90 +++ b/src/EHFinder_Generator_Sources.F90 @@ -53,9 +53,9 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS) call CCTK_InterpHandle ( interp_handle, gen_interp ) if ( interp_handle .lt. 0 ) then - warn_message = 'Cannot get handle for interpolation. ' - warn_message = trim(warn_message)//'Forgot to activate an implementation ' - warn_message = trim(warn_message)//'providing interpolation operators?' + warn_message = 'Cannot get handle for interpolation.' + warn_message = trim(warn_message)//' Forgot to activate an implementation' + warn_message = trim(warn_message)//' providing interpolation operators?' call CCTK_WARN( 0, trim(warn_message) ) end if @@ -73,9 +73,9 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS) ! Get the 3D coordinate system handle. 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 = trim(warn_message)//'Forgot to activate an implementation ' - warn_message = trim(warn_message)//'providing coordinates?' + warn_message = 'Cannot get handle for cart3d coordinate system.' + warn_message = trim(warn_message)//' Forgot to activate an implementation' + warn_message = trim(warn_message)//' providing coordinates?' call CCTK_WARN( 0, trim(warn_message) ) endif @@ -194,6 +194,18 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS) dxg(i,l) = - betaxg(i) + ssign * factor * dfux dyg(i,l) = - betayg(i) + ssign * factor * dfuy dzg(i,l) = - betazg(i) + ssign * factor * dfuz + if ( i == 5 ) then + print*,'Debug : ', xg(i,l), yg(i,l), zg(i,l), & + dxg(i,l), dyg(i,l), dzg(i,l) + print*,'lapse : ', alp2 + print*,'shift : ', betaxg(i), betayg(i), betazg(i) + print*,'metric : ', gxxg(i), gxyg(i), gxzg(i), & + gyyg(i), gyzg(i), gzzg(i) + print*,'imetric : ', guxx, guxy, guxz, guyy, guyz, guzz + print*,'df : ', dfxg(i), dfyg(i), dfzg(i) + print*,'dfu : ', dfux, dfuy, dfuz + print*,'factor : ', factor + end if end do else if ( CCTK_EQUALS ( metric_type, 'static conformal' ) ) then call Util_TableSetIntArray ( status, table_handle, 14, & @@ -260,3 +272,324 @@ subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS) return end subroutine EHFinder_Generator_Sources + + +subroutine EHFinder_Generator_Sources_2D(CCTK_ARGUMENTS) + + use EHFinder_mod + + implicit none + + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + + CCTK_INT :: i, j, l + 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 + + CCTK_INT, dimension(2) :: lsh + CCTK_POINTER, dimension(3) :: interp_coords + CCTK_POINTER, dimension(14) :: out_arrays + CCTK_INT, dimension(12) :: in_arrays + CCTK_INT, dimension(14), parameter :: op_indices = (/ 0, 1, 2, 3, 4, & + 5, 6, 7, 8, 9, & + 10, 10, 10, 11 /), & + op_codes = (/ 0, 0, 0, 0, 0, & + 0, 0, 0, 0, 0, & + 1, 2, 3, 0 /) + CCTK_INT, dimension(14) :: out_types + CCTK_REAL :: alp2, psi4, dfux, dfuy, dfuz, factor, ssign + CCTK_REAL :: idetg, guxx, guxy, guxz, guyy, guyz, guzz + CCTK_REAL, dimension(:,:), allocatable :: tmp, postmp + CCTK_INT, dimension(1) :: shape + +! 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. + call CCTK_FortranString ( gen_interp_len, generator_interpolator, & + gen_interp ) + +! Get the corresponding interpolator handle. + call CCTK_InterpHandle ( interp_handle, gen_interp ) + + if ( interp_handle .lt. 0 ) then + warn_message = 'Cannot get handle for interpolation.' + warn_message = trim(warn_message)//' Forgot to activate an implementation' + warn_message = trim(warn_message)//' providing interpolation operators?' + call CCTK_WARN( 0, trim(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. + 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 + +! Get the 3D coordinate system handle. + 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 = trim(warn_message)//' Forgot to activate an implementation' + warn_message = trim(warn_message)//' providing coordinates?' + call CCTK_WARN( 0, trim(warn_message) ) + endif + +! Find out how many interpolation points are located on this processor. + call CCTK_GrouplshGN ( status, cctkGH, 2, lsh, 'ehfinder::xg2' ) + if ( status .lt. 0 ) then + call CCTK_WARN ( 0, 'cannot get local size for surface arrays' ) + end if + + shape(1) = lsh(1)*lsh(2) + allocate ( tmp(lsh(1)*lsh(2),14), postmp(lsh(1)*lsh(2),3) ) + + do l = 1, 14 + out_arrays(l) = CCTK_PointerTo(tmp(:,l)) + end do + out_arrays(1) = CCTK_PointerTo(alpg2) + out_arrays(2) = CCTK_PointerTo(betaxg2) + out_arrays(3) = CCTK_PointerTo(betayg2) + out_arrays(4) = CCTK_PointerTo(betazg2) + out_arrays(5) = CCTK_PointerTo(gxxg2) + out_arrays(6) = CCTK_PointerTo(gxyg2) + out_arrays(7) = CCTK_PointerTo(gxzg2) + out_arrays(8) = CCTK_PointerTo(gyyg2) + out_arrays(9) = CCTK_PointerTo(gyzg2) + out_arrays(10) = CCTK_PointerTo(gzzg2) + out_arrays(11) = CCTK_PointerTo(dfxg2) + out_arrays(12) = CCTK_PointerTo(dfyg2) + out_arrays(13) = CCTK_PointerTo(dfzg2) + out_arrays(14) = CCTK_PointerTo(psig2) + +!! Set the pointers to the output arrays. + +! Loop over the level sets + do l = 1, eh_number_level_sets + +! Set the pointers to the points to be interpolated to. +! postmp(:,1) = reshape(xg2(:,:,l),shape) +! postmp(:,2) = reshape(yg2(:,:,l),shape) +! postmp(:,3) = reshape(zg2(:,:,l),shape) +! interp_coords(1) = CCTK_PointerTo(postmp(:,1)) +! interp_coords(2) = CCTK_PointerTo(postmp(:,2)) +! interp_coords(3) = CCTK_PointerTo(postmp(:,3)) + + interp_coords(1) = CCTK_PointerTo(xg2(1,1,l)) + interp_coords(2) = CCTK_PointerTo(yg2(1,1,l)) + interp_coords(3) = CCTK_PointerTo(zg2(1,1,l)) + +! Set the indices to the input grid functions. + call CCTK_VarIndex ( in_arrays(1), 'admbase::alp' ) + call CCTK_VarIndex ( in_arrays(2), 'admbase::betax' ) + call CCTK_VarIndex ( in_arrays(3), 'admbase::betay' ) + call CCTK_VarIndex ( in_arrays(4), 'admbase::betaz' ) + call CCTK_VarIndex ( in_arrays(5), 'admbase::gxx' ) + call CCTK_VarIndex ( in_arrays(6), 'admbase::gxy' ) + call CCTK_VarIndex ( in_arrays(7), 'admbase::gxz' ) + call CCTK_VarIndex ( in_arrays(8), 'admbase::gyy' ) + call CCTK_VarIndex ( in_arrays(9), 'admbase::gyz' ) + call CCTK_VarIndex ( in_arrays(10), 'admbase::gzz' ) + write(vname,'(a12,i1,a1)') 'ehfinder::f[', l-1,']' + call CCTK_VarIndex ( in_arrays(11), vname ) + call CCTK_VarIndex ( in_arrays(12), 'staticconformal::psi' ) + +! Check the metric type. At present physical and static_conformal are +! supported. + if ( CCTK_EQUALS ( metric_type, 'physical' ) ) then + +! Set the operand indices table entry, corresponding +! to interpolation of admbase::alp (1), admbase::shift (3), +! admbase::metric(6) and the first derivatives of ehfinder::f (3) for +! a total of 13 output arrays. + call Util_TableSetIntArray ( status, table_handle, 13, & + op_indices(1:13), 'operand_indices' ) + if ( status .lt. 0 ) then + warn_message = 'Cannot set operand indices array in parameter table' + call CCTK_WARN ( 0, trim(warn_message) ) + 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 + warn_message = 'Cannot set operation codes array in parameter table' + call CCTK_WARN ( 0, trim(warn_message) ) + endif + +! Call the interpolator. + call CCTK_InterpGridArrays ( status, cctkGH, 3, interp_handle, & + table_handle, coord_system_handle, & + lsh(1)*lsh(2), CCTK_VARIABLE_REAL, & + interp_coords, 11, in_arrays(1:11), & + 13, out_types(1:13), out_arrays(1:13) ) + + if ( status .lt. 0 ) then + call CCTK_INFO ( 'Interpolation failed.' ) + end if + +! alpg2 = reshape(tmp(:,1),lsh) +! betaxg2 = reshape(tmp(:,2),lsh) +! betayg2 = reshape(tmp(:,3),lsh) +! betazg2 = reshape(tmp(:,4),lsh) +! gxxg2 = reshape(tmp(:,5),lsh) +! gxyg2 = reshape(tmp(:,6),lsh) +! gxzg2 = reshape(tmp(:,7),lsh) +! gyyg2 = reshape(tmp(:,8),lsh) +! gyzg2 = reshape(tmp(:,9),lsh) +! gzzg2 = reshape(tmp(:,10),lsh) +! dfxg2 = reshape(tmp(:,11),lsh) +! dfyg2 = reshape(tmp(:,12),lsh) +! dfzg2 = reshape(tmp(:,13),lsh) + +! For each point on this processor calculate the right hand side of the +! characteristic evolution equation. + do j = 1, lsh(2) + do i = 1, lsh(1) + +! calculate the square of the lapse. + alp2 = alpg2(i,j)**2 + +! Calculate the inverse of the 3-metric. + guxx = gyyg2(i,j) * gzzg2(i,j) - gyzg2(i,j)**2 + guxy = gxzg2(i,j) * gyzg2(i,j) - gxyg2(i,j) * gzzg2(i,j) + guxz = gxyg2(i,j) * gyzg2(i,j) - gxzg2(i,j) * gyyg2(i,j) + + idetg = one / ( gxxg2(i,j) * guxx + gxyg2(i,j) * guxy & + + gxzg2(i,j) * guxz ) + + guxx = idetg * guxx + guxy = idetg * guxy + guxz = idetg * guxz + + guyy = ( gxxg2(i,j) * gzzg2(i,j) - gxzg2(i,j)**2 ) * idetg + guyz = ( gxyg2(i,j) * gxzg2(i,j) - gxxg2(i,j) * gyzg2(i,j) ) * idetg + guzz = ( gxxg2(i,j) * gyyg2(i,j) - gxyg2(i,j)**2 ) * idetg + +! Raise the index of the partial derivatives of f. + dfux = guxx * dfxg2(i,j) + guxy * dfyg2(i,j) + guxz * dfzg2(i,j) + dfuy = guxy * dfxg2(i,j) + guyy * dfyg2(i,j) + guyz * dfzg2(i,j) + dfuz = guxz * dfxg2(i,j) + guyz * dfyg2(i,j) + guzz * dfzg2(i,j) + +! Calculate the overall multiplication factor. + factor = alp2 / sqrt ( alp2 * ( dfux * dfxg2(i,j) + & + dfuy * dfyg2(i,j) + & + dfuz * dfzg2(i,j) ) ) + +! Finally obtain dx^i/dt. + dxg2(i,j,l) = - betaxg2(i,j) + ssign * factor * dfux + dyg2(i,j,l) = - betayg2(i,j) + ssign * factor * dfuy + dzg2(i,j,l) = - betazg2(i,j) + ssign * factor * dfuz + if ( j == 1 .and. i == 5 ) then + + print*,'debug2 : ', xg2(i,j,l), yg2(i,j,l), zg2(i,j,l), & + dxg2(i,j,l), dyg2(i,j,l), dzg2(i,j,l) + print*,'lapse : ', alp2 + print*,'shift : ', betaxg2(i,j), betayg2(i,j), betazg2(i,j) + print*,'metric : ', gxxg2(i,j), gxyg2(i,j), gxzg2(i,j), & + gyyg2(i,j), gyzg2(i,j), gzzg2(i,j) + print*,'imetric : ', guxx, guxy, guxz, guyy, guyz, guzz + print*,'df : ', dfxg2(i,j), dfyg2(i,j), dfzg2(i,j) + print*,'dfu : ', dfux, dfuy, dfuz + print*,'factor : ', factor + end if + end do + 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 + warn_message = 'Cannot set operand indices array in parameter table' + call CCTK_WARN ( 0, trim(warn_message) ) + endif + + call Util_TableSetIntArray ( status, table_handle, 14, & + op_codes, 'operation_codes' ) + if ( status .lt. 0 ) then + warn_message = 'Cannot set operation codes array in parameter table' + call CCTK_WARN ( 0, trim(warn_message) ) + endif + + call CCTK_InterpGridArrays ( status, cctkGH, 3, interp_handle, & + table_handle, coord_system_handle, & + lsh(1)*lsh(2), 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 + +! alpg2 = reshape(tmp(:,1),lsh) +! betaxg2 = reshape(tmp(:,2),lsh) +! betayg2 = reshape(tmp(:,3),lsh) +! betazg2 = reshape(tmp(:,4),lsh) +! gxxg2 = reshape(tmp(:,5),lsh) +! gxyg2 = reshape(tmp(:,6),lsh) +! gxzg2 = reshape(tmp(:,7),lsh) +! gyyg2 = reshape(tmp(:,8),lsh) +! gyzg2 = reshape(tmp(:,9),lsh) +! gzzg2 = reshape(tmp(:,10),lsh) +! dfxg2 = reshape(tmp(:,11),lsh) +! dfyg2 = reshape(tmp(:,12),lsh) +! dfzg2 = reshape(tmp(:,13),lsh) +! psig2 = reshape(tmp(:,14),lsh) + + do j = 1, lsh(2) + do i = 1, lsh(1) + alp2 = alpg2(i,j)**2 + +! The inverse of psi^4 + psi4 = one / psig2(i,j)**4 + + guxx = gyyg2(i,j) * gzzg2(i,j) - gyzg2(i,j)**2 + guxy = gxzg2(i,j) * gyzg2(i,j) - gxyg2(i,j) * gzzg2(i,j) + guxz = gxyg2(i,j) * gyzg2(i,j) - gxzg2(i,j) * gyyg2(i,j) + +! The inverse of the determinant divided by psi^4. + idetg = psi4 / ( gxxg2(i,j) * guxx + & + gxyg2(i,j) * guxy + & + gxzg2(i,j) * 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 = ( gxxg2(i,j) * gzzg2(i,j) - gxzg2(i,j)**2 ) * idetg + guyz = ( gxyg2(i,j) * gxzg2(i,j) - gxxg2(i,j) * gyzg2(i,j) ) * idetg + guzz = ( gxxg2(i,j) * gyyg2(i,j) - gxyg2(i,j)**2 ) * idetg + + dfux = guxx * dfxg2(i,j) + guxy * dfyg2(i,j) + guxz * dfzg2(i,j) + dfuy = guxy * dfxg2(i,j) + guyy * dfyg2(i,j) + guyz * dfzg2(i,j) + dfuz = guxz * dfxg2(i,j) + guyz * dfyg2(i,j) + guzz * dfzg2(i,j) + factor = alp2 / sqrt ( alp2 * ( dfux * dfxg2(i,j) + & + dfuy * dfyg2(i,j) + & + dfuz * dfzg2(i,j) ) ) + dxg2(i,j,l) = - betaxg2(i,j) + ssign * factor * dfux + dyg2(i,j,l) = - betayg2(i,j) + ssign * factor * dfuy + dzg2(i,j,l) = - betazg2(i,j) + ssign * factor * dfuz + + end do + end do + end if + end do + + deallocate ( tmp, postmp ) + return +end subroutine EHFinder_Generator_Sources_2D diff --git a/src/EHFinder_Init.F90 b/src/EHFinder_Init.F90 index b7f2a25..c9dd60c 100644 --- a/src/EHFinder_Init.F90 +++ b/src/EHFinder_Init.F90 @@ -22,7 +22,9 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS) CCTK_REAL :: cosa, sina, cosb, sinb, cosc, sinc CCTK_REAL :: last_time CCTK_REAL :: theta, dtheta, thetamin, thetamax, r_el + CCTK_REAL :: phi, dphi, phimin, phimax CCTK_INT, dimension(1) :: lsh, lbnd + CCTK_INT, dimension(2) :: lsh2, lbnd2 ! Get the size of the local grid. nx = cctk_lsh(1) @@ -52,41 +54,76 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS) if ( evolve_generators .gt. 0 ) then - call CCTK_GrouplbndGN ( status, cctkGH, 1, lbnd, 'ehfinder::xg' ) - if ( status .lt. 0 ) then - call CCTK_WARN ( 0, 'cannot get lower bounds for generator arrays' ) - end if - call CCTK_GrouplshGN ( status, cctkGH, 1, lsh, 'ehfinder::xg' ) - if ( status .lt. 0 ) then - call CCTK_WARN ( 0, 'cannot get local size for generator arrays' ) - end if - - if ( CCTK_EQUALS( generator_distribution, 'line' ) ) then - - if ( CCTK_EQUALS( domain, 'full' ) ) then + if ( CCTK_EQUALS( domain, 'full' ) ) then + thetamin = zero; thetamax = pi + phimin = zero; phimax = 2*pi + else if ( CCTK_EQUALS( domain, 'bitant') ) then + if ( CCTK_EQUALS( bitant_plane, 'xy' ) ) then + thetamin = zero; thetamax = half * pi + phimin = zero; phimax = 2*pi + else if ( CCTK_EQUALS( bitant_plane, 'xz' ) ) then thetamin = zero; thetamax = pi - else if ( CCTK_EQUALS( domain, 'bitant') ) then - if ( CCTK_EQUALS( bitant_plane, 'xy' ) ) then - thetamin = zero; thetamax = half * pi - else - thetamin = zero; thetamax = pi - end if - else if ( CCTK_EQUALS( domain, 'quadrant' ) ) then - if ( CCTK_EQUALS( quadrant_direction, 'x' ) .or. & - CCTK_EQUALS( quadrant_direction, 'y' ) ) then - thetamin = zero; thetamax = half * pi - else - thetamin = zero; thetamax = pi - end if - else if ( CCTK_EQUALS( domain, 'octant' ) ) then + phimin = zero; phimax = pi + else + thetamin = zero; thetamax = pi + phimin = -half * pi; phimax = half * pi + end if + else if ( CCTK_EQUALS( domain, 'quadrant' ) ) then + if ( CCTK_EQUALS( quadrant_direction, 'x' ) ) then thetamin = zero; thetamax = half * pi + phimin = zero; phimax = pi + else if ( CCTK_EQUALS( quadrant_direction, 'y' ) ) then + thetamin = zero; thetamax = half * pi + phimin = -half * pi; phimax = half * pi + else + thetamin = zero; thetamax = pi + phimin = zero; phimax = half * pi end if + else if ( CCTK_EQUALS( domain, 'octant' ) ) then + thetamin = zero; thetamax = half * pi + phimin = zero; phimax = half * pi + end if + if ( CCTK_EQUALS( generator_distribution, 'line' ) ) then + + call CCTK_GrouplbndGN ( status, cctkGH, 1, lbnd, 'ehfinder::xg' ) + if ( status .lt. 0 ) then + call CCTK_WARN ( 0, 'cannot get lower bounds for generator arrays' ) + end if + call CCTK_GrouplshGN ( status, cctkGH, 1, lsh, 'ehfinder::xg' ) + if ( status .lt. 0 ) then + call CCTK_WARN ( 0, 'cannot get local size for generator arrays' ) + end if + if ( number_of_generators .eq. 1 ) then theta = half * ( thetamax - thetamin ) + thetamin else dtheta = ( thetamax - thetamin ) / ( number_of_generators - 1 ) end if + + else if ( CCTK_EQUALS( generator_distribution, '2D array' ) ) then + + call CCTK_GrouplbndGN ( status, cctkGH, 2, lbnd2, 'ehfinder::xg2' ) + if ( status .lt. 0 ) then + call CCTK_WARN ( 0, 'cannot get lower bounds for generator arrays' ) + end if + call CCTK_GrouplshGN ( status, cctkGH, 2, lsh2, 'ehfinder::xg2' ) + if ( status .lt. 0 ) then + call CCTK_WARN ( 0, 'cannot get local size for generator arrays' ) + end if + + if ( number_of_generators_theta .eq. 1 ) then + theta = half * ( thetamax - thetamin ) + thetamin + else + dtheta = ( thetamax - thetamin ) / ( number_of_generators_theta - 1 ) + end if + + if ( number_of_generators_phi .eq. 1 ) then + phi = half * ( phimax - phimin ) + phimin + else + dphi = ( phimax - phimin ) / ( number_of_generators_phi - 1 ) + end if + end if end if @@ -103,10 +140,21 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS) if ( evolve_generators .gt. 0 ) then if ( CCTK_EQUALS( generator_distribution, 'line' ) ) then do i = 1, lsh(1) - theta = thetamin + dtheta * ( i + lbnd(1) - 1 ) - xg(i,l) = initial_rad(l) * sin(theta) + translate_x(l) - yg(i,l) = translate_y(l) - zg(i,l) = initial_rad(l) * cos(theta) + translate_z(l) + theta = thetamin + dtheta * ( i + lbnd(1) - 1 ) + xg(i,l) = initial_rad(l) * sin(theta) + translate_x(l) + yg(i,l) = translate_y(l) + zg(i,l) = initial_rad(l) * cos(theta) + translate_z(l) + end do + else if ( CCTK_EQUALS( generator_distribution, '2D array' ) ) then + do j = 1, lsh2(2) + phi = phimin + dphi * ( j + lbnd2(2) - 1 ) + do i = 1, lsh2(1) + theta = thetamin + dtheta * ( i + lbnd2(1) - 1 ) + xg2(i,j,l) = initial_rad(l) * sin(theta) * cos(phi) + translate_x(l) + yg2(i,j,l) = initial_rad(l) * sin(theta) * sin(phi) + translate_x(l) + zg2(i,j,l) = initial_rad(l) * cos(theta) + translate_z(l) + print*,xg2(i,j,l),yg2(i,j,l),zg2(i,j,l) + end do end do end if end if diff --git a/src/MoL_Init.F90 b/src/MoL_Init.F90 index 61a84e6..d01eec1 100644 --- a/src/MoL_Init.F90 +++ b/src/MoL_Init.F90 @@ -26,20 +26,41 @@ subroutine EHFinder_MoLRegister(CCTK_ARGUMENTS) if ( evolve_generators .gt. 0 ) then - call CCTK_GroupIndex (varindex, 'ehfinder::xg') - call CCTK_GroupIndex(rhsindex, 'ehfinder::dxg') + if ( CCTK_EQUALS( generator_distribution, 'line' ) ) then - ierr_cum = ierr_cum + MoLRegisterEvolvedGroup(varindex, rhsindex) + call CCTK_GroupIndex (varindex, 'ehfinder::xg') + call CCTK_GroupIndex(rhsindex, 'ehfinder::dxg') - call CCTK_GroupIndex (varindex, 'ehfinder::yg') - call CCTK_GroupIndex(rhsindex, 'ehfinder::dyg') + ierr_cum = ierr_cum + MoLRegisterEvolvedGroup(varindex, rhsindex) - ierr_cum = ierr_cum + MoLRegisterEvolvedGroup(varindex, rhsindex) + call CCTK_GroupIndex (varindex, 'ehfinder::yg') + call CCTK_GroupIndex(rhsindex, 'ehfinder::dyg') - call CCTK_GroupIndex (varindex, 'ehfinder::zg') - call CCTK_GroupIndex(rhsindex, 'ehfinder::dzg') + ierr_cum = ierr_cum + MoLRegisterEvolvedGroup(varindex, rhsindex) - ierr_cum = ierr_cum + MoLRegisterEvolvedGroup(varindex, rhsindex) + call CCTK_GroupIndex (varindex, 'ehfinder::zg') + call CCTK_GroupIndex(rhsindex, 'ehfinder::dzg') + + ierr_cum = ierr_cum + MoLRegisterEvolvedGroup(varindex, rhsindex) + + else if ( CCTK_EQUALS( generator_distribution, '2D array' ) ) then + + call CCTK_GroupIndex (varindex, 'ehfinder::xg2') + call CCTK_GroupIndex(rhsindex, 'ehfinder::dxg2') + + ierr_cum = ierr_cum + MoLRegisterEvolvedGroup(varindex, rhsindex) + + call CCTK_GroupIndex (varindex, 'ehfinder::yg2') + call CCTK_GroupIndex(rhsindex, 'ehfinder::dyg2') + + ierr_cum = ierr_cum + MoLRegisterEvolvedGroup(varindex, rhsindex) + + call CCTK_GroupIndex (varindex, 'ehfinder::zg2') + call CCTK_GroupIndex(rhsindex, 'ehfinder::dzg2') + + ierr_cum = ierr_cum + MoLRegisterEvolvedGroup(varindex, rhsindex) + + end if end if -- cgit v1.2.3