From a454628dc04978583a1189f51b62493418b9f48b Mon Sep 17 00:00:00 2001 From: diener Date: Mon, 26 May 2003 07:48:15 +0000 Subject: Changed to use the newest version of MoL. Used its support for evolving grid arrays to introduce tracking of the generators of the horizon along with the horizon itself. Still needs some work on setting up the initial generators (now it only does one) but first it needs to be tested. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/EHFinder/trunk@108 2a26948c-0e4f-0410-aee8-f1d3e353619c --- interface.ccl | 21 ++++- param.ccl | 40 +++++++++- schedule.ccl | 14 ++++ src/EHFinder_Generator_Sources.F90 | 157 +++++++++++++++++++++++++++++++++++++ src/EHFinder_Init.F90 | 5 ++ src/EHFinder_Integrate2.F90 | 8 +- src/MoL_Init.F90 | 25 +++++- src/make.code.defn | 1 + src/make.code.deps | 29 +++---- 9 files changed, 276 insertions(+), 24 deletions(-) create mode 100644 src/EHFinder_Generator_Sources.F90 diff --git a/interface.ccl b/interface.ccl index b061a46..ebb03d0 100644 --- a/interface.ccl +++ b/interface.ccl @@ -5,7 +5,11 @@ implements: ehfinder inherits: grid admbase coordgauge staticconformal spacemask boundary USES INCLUDE: Boundary.h -USES INCLUDE: MoL.h +#USES INCLUDE: MoL.h + +CCTK_INT FUNCTION MoLRegisterEvolved(CCTK_INT IN EvolvedIndex, \ + CCTK_INT IN RHSIndex) +USES FUNCTION MoLRegisterEvolved public: @@ -163,3 +167,18 @@ CCTK_REAL eh_circumference2 TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=maximum_surface_n CCTK_REAL eh_group_test[maximum_surface_number] TYPE=SCALAR CCTK_REAL eh_group_array_test[maximum_surface_number] TYPE=ARRAY DIM=1 SIZE=100 GHOSTSIZE=0 DISTRIB=DEFAULT + +CCTK_REAL generators TYPE=ARRAY DIM=1 TIMELEVELS=2 SIZE=number_of_generators GHOSTSIZE=0 DISTRIB=DEFAULT +{ + xg, yg, zg +} "The position of the generators of the event horizon" + +CCTK_REAL generators_rhs TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=number_of_generators GHOSTSIZE=0 DISTRIB=DEFAULT +{ + dxg, dyg, dzg +} + +CCTK_REAL generators_arrays TYPE=ARRAY DIM=1 TIMELEVELS=1 SIZE=number_of_generators GHOSTSIZE=0 DISTRIB=DEFAULT +{ + alpg, betaxg, betayg, betazg, gxxg, gxyg, gxzg, gyyg, gyzg, gzzg, dfxg, dfyg, dfzg, psig +} diff --git a/param.ccl b/param.ccl index 0eefed4..f54d5e8 100644 --- a/param.ccl +++ b/param.ccl @@ -198,19 +198,28 @@ BOOLEAN use_user_center "Should the user prescribed center be used" CCTK_REAL center_x "The x-coordinate of the center" { -*:* :: "Anything" + *:* :: "Anything" } 0.0 CCTK_REAL center_y "The y-coordinate of the center" { -*:* :: "Anything" + *:* :: "Anything" } 0.0 CCTK_REAL center_z "The z-coordinate of the center" { -*:* :: "Anything" + *:* :: "Anything" } 0.0 +BOOLEAN evolve_generators "Should the generators be evolved" +{ +} "no" + +CCTK_INT number_of_generators "How many generators should be evolved" +{ + 1:* :: "Postive please" +} 1 + shares: grid USES KEYWORD domain @@ -239,3 +248,28 @@ EXTENDS KEYWORD initial_shift shares: spacemask USES KEYWORD use_mask + +shares: MethodOfLines + +USES CCTK_INT MoL_Num_Evolved_Vars +USES CCTK_INT MoL_Max_Evolved_Array_Size +USES CCTK_INT MoL_Num_ArrayEvolved_Vars + +restricted: + +CCTK_INT EHFinder_MaxNumEvolvedVars "The maximum number of evolved variables used by EHFinder" ACCUMULATOR-BASE=MethodOfLines::MoL_Num_Evolved_Vars +{ + 1:1 :: "Only evolve the level set function" +} 1 + +CCTK_INT EHFinder_Max_Evolved_Array_Size "The maximum size of evolved grid arrays used by EHFinder" ACCUMULATOR-BASE=MethodOfLines::MoL_Max_Evolved_Array_Size +{ + 1:* :: "The size of the generator grid arrays" +} 1 + +CCTK_INT EHFinder_Num_ArrayEvolved_Vars "The maximum number of evolved grid arrays used by EHFinder" ACCUMULATOR-BASE=MethodOfLines::MoL_Num_ArrayEvolved_Vars +{ + 0:3 :: "Should be exactly zero or three" +} 3 + + diff --git a/schedule.ccl b/schedule.ccl index 3e5f227..c234ae5 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -21,6 +21,12 @@ else STORAGE: eh_area2, eh_centroid2, eh_circumference2 STORAGE: find_surface_status # STORAGE: surface_index + if ( evolve_generators ) + { + STORAGE: generators[2] + STORAGE: generators_rhs + STORAGE: generators_arrays + } } # Check for metric_state @@ -297,6 +303,14 @@ if (CCTK_Equals(mode,"normal")) LANG: Fortran } "Calculate the source terms" + if ( evolve_generators) + { + schedule EHFinder_Generator_Sources in MoL_CalcRHS + { + LANG: Fortran + } "Calculate the source terms for the generator evolution" + } + schedule GROUP EHFinder_PostStep in MoL_PostStep { } "Schedule group for symmetry boundaries and syncing" diff --git a/src/EHFinder_Generator_Sources.F90 b/src/EHFinder_Generator_Sources.F90 new file mode 100644 index 0000000..924588f --- /dev/null +++ b/src/EHFinder_Generator_Sources.F90 @@ -0,0 +1,157 @@ +! Calculation of the sources for the level set function. +! $Header$ + +#include "cctk.h" +#include "cctk_Parameters.h" +#include "cctk_Arguments.h" + +subroutine EHFinder_Generator_Sources(CCTK_ARGUMENTS) + + use EHFinder_mod + + implicit none + + DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + + CCTK_INT :: i + CCTK_INT :: interp_handle, table_handle, status, coord_system_handle + CCTK_INT, dimension(1) :: 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, 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 + + 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 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 + + 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" ) + + 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 + + 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 CCTK_InterpGridArrays ( status, cctkGH, 3, interp_handle, & + table_handle, coord_system_handle, & + lsh(1), CCTK_VARIABLE_REAL, & + interp_coords, 11, in_arrays(1:11), & + 13, out_types(1:13), out_arrays(1:13) ) + + if ( status .lt. 0 ) then + call CCTK_INFO ( 'Interpolation failed.' ) + end if + do i = 1, lsh(1) + alp2 = alpg(i)**2 + + guxx = gyyg(i) * gzzg(i) - gyzg(i)**2 + guxy = gxzg(i) * gyzg(i) - gxyg(i) * gzzg(i) + guxz = gxyg(i) * gyzg(i) - gxzg(i) * gyyg(i) + + idetg = one / ( gxxg(i) * guxx + gxyg(i) * guxy + gxzg(i) * guxz ) + + guxx = idetg * guxx + guxy = idetg * guxy + guxz = idetg * guxz + + guyy = ( gxxg(i) * gzzg(i) - gxzg(i)**2 ) * idetg + guyz = ( gxyg(i) * gxzg(i) - gxxg(i) * gyzg(i) ) * idetg + guzz = ( gxxg(i) * gyyg(i) - gxyg(i)**2 ) * idetg + + dfux = guxx * dfxg(i) + guxy * dfyg(i) + guxz * dfzg(i) + dfuy = guxy * dfxg(i) + guyy * dfyg(i) + guyz * dfzg(i) + dfuz = guxz * dfxg(i) + guyz * dfyg(i) + guzz * dfzg(i) + factor = alp2 / sqrt ( alp2 * ( dfux * dfxg(i) + & + dfuy * dfyg(i) + & + dfuz * dfzg(i) ) ) + dxg(i) = - betaxg(i) + 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 +end subroutine EHFinder_Generator_Sources diff --git a/src/EHFinder_Init.F90 b/src/EHFinder_Init.F90 index cf700f3..ebb7cb7 100644 --- a/src/EHFinder_Init.F90 +++ b/src/EHFinder_Init.F90 @@ -41,6 +41,11 @@ subroutine EHFinder_Init_F(CCTK_ARGUMENTS) f = sqrt( ( x - translate_x )**2 + & ( y - translate_y )**2 + & ( z - translate_z )**2 ) - initial_rad + if ( evolve_generators .gt. 0 ) then + xg(1) = initial_rad + translate_x + yg(1) = translate_y + zg(1) = translate_z + end if end if ! If an ellipsoid is requested... diff --git a/src/EHFinder_Integrate2.F90 b/src/EHFinder_Integrate2.F90 index 81745ea..48bd609 100644 --- a/src/EHFinder_Integrate2.F90 +++ b/src/EHFinder_Integrate2.F90 @@ -35,7 +35,7 @@ subroutine EHFinder_FindSurfaceElement(CCTK_ARGUMENTS) CCTK_VARIABLE_REAL, & CCTK_VARIABLE_REAL /) -! If finding of surface failed don't try to integrate but exit. +! If finding of surface failed do not try to integrate but exit. if ( find_surface_status .lt. 0 ) then return endif @@ -212,7 +212,7 @@ subroutine EHFinder_IntegrateArea(CCTK_ARGUMENTS) character(len=30) :: info_message CCTK_INT, dimension(1) :: lbnd, ubnd, lsh -! If finding of surface failed don't try to integrate but exit. +! If finding of surface failed do not try to integrate but exit. if ( find_surface_status .lt. 0 ) then return endif @@ -274,7 +274,7 @@ subroutine EHFinder_IntegrateCentroid(CCTK_ARGUMENTS) CCTK_INT, dimension(1) :: lbnd, ubnd, lsh CCTK_REAL :: centroid_x, centroid_y, centroid_z -! If finding of surface failed don't try to integrate but exit. +! If finding of surface failed do not try to integrate but exit. if ( find_surface_status .lt. 0 ) then return endif @@ -382,7 +382,7 @@ subroutine EHFinder_IntegrateEquatorial(CCTK_ARGUMENTS) 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 don't try to integrate but exit. +! If finding of surface failed do not try to integrate but exit. if ( find_surface_status .lt. 0 ) then return endif diff --git a/src/MoL_Init.F90 b/src/MoL_Init.F90 index 679c8fb..bccc870 100644 --- a/src/MoL_Init.F90 +++ b/src/MoL_Init.F90 @@ -10,6 +10,7 @@ subroutine EHFinder_MoLRegister(CCTK_ARGUMENTS) implicit none CCTK_INT :: ierr, ierr_cum, varindex, rhsindex + CCTK_INT :: MoLRegisterEvolved DECLARE_CCTK_PARAMETERS DECLARE_CCTK_ARGUMENTS @@ -18,8 +19,28 @@ subroutine EHFinder_MoLRegister(CCTK_ARGUMENTS) ierr_cum = 0 call CCTK_VarIndex(varindex, 'ehfinder::f') call CCTK_VarIndex(rhsindex, 'ehfinder::sf') - call MoL_RegisterVar(ierr, varindex, rhsindex) - ierr_cum = ierr_cum + ierr +! call MoL_RegisterVar(ierr, varindex, rhsindex) +! ierr_cum = ierr_cum + ierr + ierr_cum = ierr_cum + MoLRegisterEvolved(varindex, rhsindex) + + if ( evolve_generators .gt. 0 ) then + +! call CCTK_GroupIndex (varindex, 'ehfinder::generators') +! call CCTK_GroupIndex(rhsindex, 'ehfinder::generators_rhs') + + call CCTK_VarIndex(varindex, 'ehfinder::xg') + call CCTK_VarIndex(rhsindex, 'ehfinder::dxg') + ierr_cum = ierr_cum + MoLRegisterEvolved(varindex, rhsindex) + + call CCTK_VarIndex(varindex, 'ehfinder::yg') + call CCTK_VarIndex(rhsindex, 'ehfinder::dyg') + ierr_cum = ierr_cum + MoLRegisterEvolved(varindex, rhsindex) + + call CCTK_VarIndex(varindex, 'ehfinder::zg') + call CCTK_VarIndex(rhsindex, 'ehfinder::dzg') + ierr_cum = ierr_cum + MoLRegisterEvolved(varindex, rhsindex) + end if + ! call CCTK_VarIndex(varindex, 'admbase::gxx') ! call MoL_RegisterPrimitive(ierr, varindex) ! ierr_cum = ierr_cum + ierr diff --git a/src/make.code.defn b/src/make.code.defn index c54b696..c8e8ce4 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -7,6 +7,7 @@ SRCS = EHFinder_mod.F90 \ EHFinder_Init.F90 \ MoL_Init.F90 \ EHFinder_Sources.F90 \ + EHFinder_Generator_Sources.F90 \ EHFinder_ReParametrize.F90 \ EHFinder_Check.F90 \ EHFinder_Sources2.F90 \ diff --git a/src/make.code.deps b/src/make.code.deps index 1c09018..d8e8918 100644 --- a/src/make.code.deps +++ b/src/make.code.deps @@ -1,17 +1,18 @@ # Dependencies file for thorn EHFinder. # $Header$ -EHFinder_mod.F90.o: EHFinder_Constants.F90.o -EHFinder_Init.F90.o: EHFinder_mod.F90.o -EHFinder_Sources.F90.o: EHFinder_mod.F90.o -EHFinder_ReParametrize.F90.o: EHFinder_mod.F90.o -EHFinder_Check.F90.o: EHFinder_mod.F90.o -EHFinder_Sources2.F90.o: EHFinder_mod.F90.o -EHFinder_Sources3.F90.o: EHFinder_mod.F90.o -EHFinder_SetMask.F90.o: EHFinder_mod.F90.o -EHFinder_SetSym.F90.o: EHFinder_mod.F90.o -EHFinder_ReadData.F90.o: EHFinder_mod.F90.o -EHFinder_Check.F90.o: EHFinder_mod.F90.o -EHFinder_Integrate.F90.o: EHFinder_mod.F90.o -EHFinder_Integrate2.F90.o: EHFinder_mod.F90.o -EHFinder_FindSurface.F90.o: EHFinder_mod.F90.o +EHFinder_mod.F90.o: EHFinder_Constants.F90.o +EHFinder_Init.F90.o: EHFinder_mod.F90.o +EHFinder_Sources.F90.o: EHFinder_mod.F90.o +EHFinder_Generator_Sources.F90.o: EHFinder_mod.F90.o +EHFinder_ReParametrize.F90.o: EHFinder_mod.F90.o +EHFinder_Check.F90.o: EHFinder_mod.F90.o +EHFinder_Sources2.F90.o: EHFinder_mod.F90.o +EHFinder_Sources3.F90.o: EHFinder_mod.F90.o +EHFinder_SetMask.F90.o: EHFinder_mod.F90.o +EHFinder_SetSym.F90.o: EHFinder_mod.F90.o +EHFinder_ReadData.F90.o: EHFinder_mod.F90.o +EHFinder_Check.F90.o: EHFinder_mod.F90.o +EHFinder_Integrate.F90.o: EHFinder_mod.F90.o +EHFinder_Integrate2.F90.o: EHFinder_mod.F90.o +EHFinder_FindSurface.F90.o: EHFinder_mod.F90.o -- cgit v1.2.3