diff options
| -rw-r--r-- | interface.ccl | 21 | ||||
| -rw-r--r-- | param.ccl | 40 | ||||
| -rw-r--r-- | schedule.ccl | 14 | ||||
| -rw-r--r-- | src/EHFinder_Generator_Sources.F90 | 157 | ||||
| -rw-r--r-- | src/EHFinder_Init.F90 | 5 | ||||
| -rw-r--r-- | src/EHFinder_Integrate2.F90 | 8 | ||||
| -rw-r--r-- | src/MoL_Init.F90 | 25 | ||||
| -rw-r--r-- | src/make.code.defn | 1 | ||||
| -rw-r--r-- | src/make.code.deps | 29 |
9 files changed, 276 insertions, 24 deletions
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 +} @@ -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 |