From fd5a4312226842434f81a445910e3284b7a77b43 Mon Sep 17 00:00:00 2001 From: schnetter Date: Wed, 7 Apr 2010 17:00:34 +0000 Subject: Correct nans git-svn-id: https://svn.cct.lsu.edu/repos/numrel/LSUThorns/QuasiLocalMeasures/trunk@5 ef6f4158-a8ee-47d1-ba14-cb73256398e0 --- interface.ccl | 26 +- schedule.ccl | 12 +- src/qlm_3determinant.F90 | 2 +- src/qlm_broadcast.c | 3 - src/qlm_calculate.F90 | 3 + src/qlm_coordinates.F90 | 10 +- src/qlm_import_surface.F90 | 32 --- src/qlm_interpolate.F90 | 688 ++++++++++++++++++++++++++------------------- src/qlm_multipoles.F90 | 161 ++++++----- src/qlm_tetrad.F90 | 22 +- src/qlm_variables.F90 | 47 +++- src/qlm_weyl_scalars.F90 | 57 ++++ 12 files changed, 618 insertions(+), 445 deletions(-) diff --git a/interface.ccl b/interface.ccl index fffdbf8..0c01900 100644 --- a/interface.ccl +++ b/interface.ccl @@ -2,7 +2,7 @@ IMPLEMENTS: QuasiLocalMeasures -INHERITS: ADMBase SphericalSurface +INHERITS: ADMBase SphericalSurface TmunuBase @@ -86,24 +86,6 @@ COMPLEX qlm_tetrad_m[num_surfaces] TYPE=array DIM=2 SIZE=SphericalSurface::maxnt qlm_m0 qlm_m1 qlm_m2 qlm_m3 # null normal within surface } "Tetrad vector m^mu" -REAL qlm_tetrad_l_p[num_surfaces] TYPE=array DIM=2 SIZE=SphericalSurface::maxntheta,SphericalSurface::maxnphi DISTRIB=constant TAGS='convergence_power=1' -{ - qlm_l0_p qlm_l1_p qlm_l2_p qlm_l3_p - qlm_l0_p_p qlm_l1_p_p qlm_l2_p_p qlm_l3_p_p -} "Previous tetrad vectors l^mu" - -REAL qlm_tetrad_n_p[num_surfaces] TYPE=array DIM=2 SIZE=SphericalSurface::maxntheta,SphericalSurface::maxnphi DISTRIB=constant TAGS='convergence_power=1' -{ - qlm_n0_p qlm_n1_p qlm_n2_p qlm_n3_p - qlm_n0_p_p qlm_n1_p_p qlm_n2_p_p qlm_n3_p_p -} "Previous tetrad vectors n^mu" - -COMPLEX qlm_tetrad_m_p[num_surfaces] TYPE=array DIM=2 SIZE=SphericalSurface::maxntheta,SphericalSurface::maxnphi DISTRIB=constant TAGS='convergence_power=1' -{ - qlm_m0_p qlm_m1_p qlm_m2_p qlm_m3_p - qlm_m0_p_p qlm_m1_p_p qlm_m2_p_p qlm_m3_p_p -} "Previous tetrad vectors m^mu" - COMPLEX qlm_newman_penrose[num_surfaces] TYPE=array DIM=2 SIZE=SphericalSurface::maxntheta,SphericalSurface::maxnphi DISTRIB=constant TAGS='Checkpoint="no"' TAGS='convergence_power=1' { qlm_npkappa qlm_nptau qlm_npsigma qlm_nprho @@ -156,12 +138,6 @@ REAL qlm_killing_vector[num_surfaces] TYPE=array DIM=2 SIZE=SphericalSurface::ma qlm_chi } "Killing vector field" -REAL qlm_killing_vector_p[num_surfaces] TYPE=array DIM=2 SIZE=SphericalSurface::maxntheta,SphericalSurface::maxnphi DISTRIB=constant TAGS='convergence_power=1' -{ - qlm_xi_t_p qlm_xi_p_p - qlm_xi_t_p_p qlm_xi_p_p_p -} "Previous Killing vector field" - REAL qlm_killed_twometric[num_surfaces] TYPE=array DIM=2 SIZE=SphericalSurface::maxntheta,SphericalSurface::maxnphi DISTRIB=constant TAGS='Checkpoint="no"' TAGS='convergence_power=1' { qlm_lqtt qlm_lqtp qlm_lqpp diff --git a/schedule.ccl b/schedule.ccl index 92ca28f..630aa04 100644 --- a/schedule.ccl +++ b/schedule.ccl @@ -4,8 +4,8 @@ STORAGE: qlm_state qlm_scalars STORAGE: qlm_state_p qlm_scalars_p STORAGE: qlm_grid_int qlm_grid_real qlm_grid_real_p STORAGE: qlm_shapes qlm_tetrad_l qlm_tetrad_n qlm_tetrad_m -STORAGE: qlm_shapes_p qlm_tetrad_l_p qlm_tetrad_n_p qlm_tetrad_m_p -STORAGE: qlm_killing_vector qlm_killing_vector_p +STORAGE: qlm_shapes_p +STORAGE: qlm_killing_vector # Ensure that there is enough storage for the ADM variables STORAGE: ADMBase::metric[3] ADMBase::curv[3] ADMBase::lapse[3] ADMBase::shift[3] @@ -23,15 +23,15 @@ SCHEDULE qlm_paramcheck AT paramcheck { OPTIONS: global LANG: Fortran -} "Check Isolated and Dynamical Horizon parameter settings" +} "Check quasi-local parameter settings" SCHEDULE qlm_init AT initial { OPTIONS: global LANG: Fortran -} "Initialise Isolated and Dynamical Horizon calculations" +} "Initialise quasi-local calculations" -SCHEDULE qlm_calculate AT analysis AFTER SphericalSurface_HasBeenSet +SCHEDULE qlm_calculate AT analysis AFTER (SphericalSurface_HasBeenSet SetTmunu) { OPTIONS: global LANG: Fortran @@ -47,4 +47,4 @@ SCHEDULE qlm_calculate AT analysis AFTER SphericalSurface_HasBeenSet TRIGGERS: qlm_ricci_scalars qlm_twometric qlm_killing_vector qlm_killed_twometric TRIGGERS: qlm_invariant_coordinates qlm_multipole_moments qlm_3determinant TRIGGERS: qlm_scalars -} "Calculate Isolated and Dynamical Horizon quantities" +} "Calculate quasi-local quantities" diff --git a/src/qlm_3determinant.F90 b/src/qlm_3determinant.F90 index d263009..1b575c1 100644 --- a/src/qlm_3determinant.F90 +++ b/src/qlm_3determinant.F90 @@ -37,7 +37,7 @@ SUBROUTINE qlm_calc_3determinant (CCTK_ARGUMENTS, hn) CCTK_REAL :: theta, phi IF (veryverbose/=0) THEN - CALL CCTK_INFO ("Computing 3-Volume Element of Horizon") + CALL CCTK_INFO ("Computing 3-Volume Element of Surface") END IF t0 = qlm_time(hn) diff --git a/src/qlm_broadcast.c b/src/qlm_broadcast.c index a491fb5..3266177 100644 --- a/src/qlm_broadcast.c +++ b/src/qlm_broadcast.c @@ -147,9 +147,6 @@ CCTK_FNAME(qlm_broadcast) (CCTK_POINTER_TO_CONST * restrict const cctkGH_) bcast (cctkGH, "QuasiLocalMeasures::qlm_tetrad_l" , hn, root); bcast (cctkGH, "QuasiLocalMeasures::qlm_tetrad_n" , hn, root); bcast (cctkGH, "QuasiLocalMeasures::qlm_tetrad_m" , hn, root); - bcast (cctkGH, "QuasiLocalMeasures::qlm_tetrad_l_p" , hn, root); - bcast (cctkGH, "QuasiLocalMeasures::qlm_tetrad_n_p" , hn, root); - bcast (cctkGH, "QuasiLocalMeasures::qlm_tetrad_m_p" , hn, root); bcast (cctkGH, "QuasiLocalMeasures::qlm_newman_penrose" , hn, root); bcast (cctkGH, "QuasiLocalMeasures::qlm_weyl_scalars" , hn, root); bcast (cctkGH, "QuasiLocalMeasures::qlm_ricci_scalars" , hn, root); diff --git a/src/qlm_calculate.F90 b/src/qlm_calculate.F90 index eb4c4ba..19b0ffb 100644 --- a/src/qlm_calculate.F90 +++ b/src/qlm_calculate.F90 @@ -106,6 +106,9 @@ subroutine qlm_calculate (CCTK_ARGUMENTS) end if call qlm_calc_3determinant (CCTK_PASS_FTOF, hn) call qlm_analyse (CCTK_PASS_FTOF, hn) + if (qlm_have_killing_vector(hn) /= 0) then + call qlm_multipoles_normalise (CCTK_PASS_FTOF, hn) + end if 9999 continue diff --git a/src/qlm_coordinates.F90 b/src/qlm_coordinates.F90 index 6807bcd..17b7b0c 100644 --- a/src/qlm_coordinates.F90 +++ b/src/qlm_coordinates.F90 @@ -9,6 +9,7 @@ subroutine qlm_calc_coordinates (CCTK_ARGUMENTS, hn) use cctk use constants use qlm_boundary + use tensor2 implicit none DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_FUNCTIONS @@ -43,7 +44,8 @@ subroutine qlm_calc_coordinates (CCTK_ARGUMENTS, hn) qq(1,2) = qlm_qtp(i,j,hn) qq(2,2) = qlm_qpp(i,j,hn) qq(2,1) = qq(1,2) - dtq = qq(1,1) * qq(2,2) - qq(1,2) * qq(2,1) + + call calc_2det (qq, dtq) area = area + sqrt(dtq) * qlm_delta_theta(hn) * qlm_delta_phi(hn) @@ -98,7 +100,8 @@ subroutine qlm_calc_coordinates (CCTK_ARGUMENTS, hn) qq(1,2) = qlm_qtp(i,j,hn) qq(2,2) = qlm_qpp(i,j,hn) qq(2,1) = qq(1,2) - dtq = qq(1,1) * qq(2,2) - qq(1,2) * qq(2,1) + + call calc_2det (qq, dtq) integral_z = integral_z + qlm_inv_z(i,j,hn) * sqrt(dtq) * qlm_delta_theta(hn) * qlm_delta_phi(hn) @@ -172,7 +175,8 @@ contains qq(1,2) = qlm_qtp(i,j,hn) qq(2,2) = qlm_qpp(i,j,hn) qq(2,1) = qq(1,2) - dtq = qq(1,1) * qq(2,2) - qq(1,2) * qq(2,1) + + call calc_2det (qq, dtq) zdot = vv(1) * (- (1/radius**2) * qlm_xi_p(i,j,hn) * sqrt(dtq)) & + vv(2) * (+ (1/radius**2) * qlm_xi_t(i,j,hn) * sqrt(dtq)) diff --git a/src/qlm_import_surface.F90 b/src/qlm_import_surface.F90 index 8eec453..d10b736 100644 --- a/src/qlm_import_surface.F90 +++ b/src/qlm_import_surface.F90 @@ -57,38 +57,6 @@ subroutine qlm_import_surface (CCTK_ARGUMENTS, hn) qlm_shape_p_p(:,:,hn) = qlm_shape_p(:,:,hn) qlm_shape_p (:,:,hn) = qlm_shape (:,:,hn) - - qlm_l0_p_p(:,:,hn) = qlm_l0_p(:,:,hn) - qlm_l0_p (:,:,hn) = qlm_l0 (:,:,hn) - qlm_l1_p_p(:,:,hn) = qlm_l1_p(:,:,hn) - qlm_l1_p (:,:,hn) = qlm_l1 (:,:,hn) - qlm_l2_p_p(:,:,hn) = qlm_l2_p(:,:,hn) - qlm_l2_p (:,:,hn) = qlm_l2 (:,:,hn) - qlm_l3_p_p(:,:,hn) = qlm_l3_p(:,:,hn) - qlm_l3_p (:,:,hn) = qlm_l3 (:,:,hn) - - qlm_n0_p_p(:,:,hn) = qlm_n0_p(:,:,hn) - qlm_n0_p (:,:,hn) = qlm_n0 (:,:,hn) - qlm_n1_p_p(:,:,hn) = qlm_n1_p(:,:,hn) - qlm_n1_p (:,:,hn) = qlm_n1 (:,:,hn) - qlm_n2_p_p(:,:,hn) = qlm_n2_p(:,:,hn) - qlm_n2_p (:,:,hn) = qlm_n2 (:,:,hn) - qlm_n3_p_p(:,:,hn) = qlm_n3_p(:,:,hn) - qlm_n3_p (:,:,hn) = qlm_n3 (:,:,hn) - - qlm_m0_p_p(:,:,hn) = qlm_m0_p(:,:,hn) - qlm_m0_p (:,:,hn) = qlm_m0 (:,:,hn) - qlm_m1_p_p(:,:,hn) = qlm_m1_p(:,:,hn) - qlm_m1_p (:,:,hn) = qlm_m1 (:,:,hn) - qlm_m2_p_p(:,:,hn) = qlm_m2_p(:,:,hn) - qlm_m2_p (:,:,hn) = qlm_m2 (:,:,hn) - qlm_m3_p_p(:,:,hn) = qlm_m3_p(:,:,hn) - qlm_m3_p (:,:,hn) = qlm_m3 (:,:,hn) - - qlm_xi_t_p_p(:,:,hn) = qlm_xi_t_p(:,:,hn) - qlm_xi_t_p (:,:,hn) = qlm_xi_t (:,:,hn) - qlm_xi_p_p_p(:,:,hn) = qlm_xi_p_p(:,:,hn) - qlm_xi_p_p (:,:,hn) = qlm_xi_p (:,:,hn) end if diff --git a/src/qlm_interpolate.F90 b/src/qlm_interpolate.F90 index 82ec27e..8317de7 100644 --- a/src/qlm_interpolate.F90 +++ b/src/qlm_interpolate.F90 @@ -29,6 +29,8 @@ subroutine qlm_interpolate (CCTK_ARGUMENTS, hn) integer, parameter :: sk = kind(interpolator) CCTK_REAL, parameter :: one = 1 + CCTK_REAL, parameter :: poison_value = -42 + integer :: len_coordsystem integer :: len_interpolator integer :: len_interpolator_options @@ -54,14 +56,17 @@ subroutine qlm_interpolate (CCTK_ARGUMENTS, hn) integer :: ind_kxx, ind_kxy, ind_kxz, ind_kyy, ind_kyz, ind_kzz integer :: ind_alpha integer :: ind_betax, ind_betay, ind_betaz + integer :: ind_ttt + integer :: ind_ttx, ind_tty, ind_ttz + integer :: ind_txx, ind_txy, ind_txz, ind_tyy, ind_tyz, ind_tzz integer :: coord_type CCTK_POINTER :: coords(3) - CCTK_INT :: inputs(16) - CCTK_INT :: output_types(88) - CCTK_POINTER :: outputs(88) - CCTK_INT :: operand_indices(88) - CCTK_INT :: operation_codes(88) + CCTK_INT :: inputs(26) + CCTK_INT :: output_types(98) + CCTK_POINTER :: outputs(98) + CCTK_INT :: operand_indices(98) + CCTK_INT :: operation_codes(98) integer :: npoints character :: msg*1000 @@ -82,6 +87,10 @@ subroutine qlm_interpolate (CCTK_ARGUMENTS, hn) call CCTK_WARN (0, "The shift must have storage") end if +!!$ if (stress_energy_state==0) then +!!$ call CCTK_WARN (0, "The stress-energy tensor must have storage") +!!$ end if + ! Get coordinate system @@ -137,22 +146,45 @@ subroutine qlm_interpolate (CCTK_ARGUMENTS, hn) ! TODO: check the excision mask ! Get variable indices - call CCTK_VarIndex (ind_gxx , "ADMBase::gxx" ) - call CCTK_VarIndex (ind_gxy , "ADMBase::gxy" ) - call CCTK_VarIndex (ind_gxz , "ADMBase::gxz" ) - call CCTK_VarIndex (ind_gyy , "ADMBase::gyy" ) - call CCTK_VarIndex (ind_gyz , "ADMBase::gyz" ) - call CCTK_VarIndex (ind_gzz , "ADMBase::gzz" ) - call CCTK_VarIndex (ind_kxx , "ADMBase::kxx" ) - call CCTK_VarIndex (ind_kxy , "ADMBase::kxy" ) - call CCTK_VarIndex (ind_kxz , "ADMBase::kxz" ) - call CCTK_VarIndex (ind_kyy , "ADMBase::kyy" ) - call CCTK_VarIndex (ind_kyz , "ADMBase::kyz" ) - call CCTK_VarIndex (ind_kzz , "ADMBase::kzz" ) - call CCTK_VarIndex (ind_alpha, "ADMBase::alp" ) - call CCTK_VarIndex (ind_betax, "ADMBase::betax") - call CCTK_VarIndex (ind_betay, "ADMBase::betay") - call CCTK_VarIndex (ind_betaz, "ADMBase::betaz") + call CCTK_VarIndex (ind_gxx , "ADMBase::gxx" ) + call CCTK_VarIndex (ind_gxy , "ADMBase::gxy" ) + call CCTK_VarIndex (ind_gxz , "ADMBase::gxz" ) + call CCTK_VarIndex (ind_gyy , "ADMBase::gyy" ) + call CCTK_VarIndex (ind_gyz , "ADMBase::gyz" ) + call CCTK_VarIndex (ind_gzz , "ADMBase::gzz" ) + call CCTK_VarIndex (ind_kxx , "ADMBase::kxx" ) + call CCTK_VarIndex (ind_kxy , "ADMBase::kxy" ) + call CCTK_VarIndex (ind_kxz , "ADMBase::kxz" ) + call CCTK_VarIndex (ind_kyy , "ADMBase::kyy" ) + call CCTK_VarIndex (ind_kyz , "ADMBase::kyz" ) + call CCTK_VarIndex (ind_kzz , "ADMBase::kzz" ) + call CCTK_VarIndex (ind_alpha, "ADMBase::alp" ) + call CCTK_VarIndex (ind_betax, "ADMBase::betax" ) + call CCTK_VarIndex (ind_betay, "ADMBase::betay" ) + call CCTK_VarIndex (ind_betaz, "ADMBase::betaz" ) + if (stress_energy_state /= 0) then + call CCTK_VarIndex (ind_ttt , "TmunuBase::eTtt") + call CCTK_VarIndex (ind_ttx , "TmunuBase::eTtx") + call CCTK_VarIndex (ind_tty , "TmunuBase::eTty") + call CCTK_VarIndex (ind_ttz , "TmunuBase::eTtz") + call CCTK_VarIndex (ind_txx , "TmunuBase::eTxx") + call CCTK_VarIndex (ind_txy , "TmunuBase::eTxy") + call CCTK_VarIndex (ind_txz , "TmunuBase::eTxz") + call CCTK_VarIndex (ind_tyy , "TmunuBase::eTyy") + call CCTK_VarIndex (ind_tyz , "TmunuBase::eTyz") + call CCTK_VarIndex (ind_tzz , "TmunuBase::eTzz") + else + ind_ttt = -1 + ind_ttx = -1 + ind_tty = -1 + ind_ttz = -1 + ind_txx = -1 + ind_txy = -1 + ind_txz = -1 + ind_tyy = -1 + ind_tyz = -1 + ind_tzz = -1 + end if @@ -170,7 +202,11 @@ subroutine qlm_interpolate (CCTK_ARGUMENTS, hn) ind_gxx, ind_gxy, ind_gxz, ind_gyy, ind_gyz, ind_gzz, & ind_kxx, ind_kxy, ind_kxz, ind_kyy, ind_kyz, ind_kzz, & ind_alpha, & - ind_betax, ind_betay, ind_betaz /) + ind_betax, ind_betay, ind_betaz, & + ind_ttt, & + ind_ttx, ind_tty, ind_ttz, & + ind_txx, ind_txy, ind_txz, ind_tyy, ind_tyz, ind_tzz /) + call CCTK_NumVars (nvars) if (nvars < 0) call CCTK_WARN (0, "internal error") if (any(inputs /= -1 .and. (inputs < 0 .or. inputs >= nvars))) then @@ -193,7 +229,10 @@ subroutine qlm_interpolate (CCTK_ARGUMENTS, hn) 06, 07, 08, 09, 10, 11, & 06, 07, 08, 09, 10, 11, & 12, & ! alp - 13, 14, 15 /) ! beta^i + 13, 14, 15, & ! beta^i + 16, & ! T_tt + 17, 18, 19, & ! T_ti + 20, 21, 22, 23, 24, 25 /) ! T_ij operation_codes = (/ & 0, 0, 0, 0, 0, 0, & ! g_ij @@ -211,7 +250,10 @@ subroutine qlm_interpolate (CCTK_ARGUMENTS, hn) 2, 2, 2, 2, 2, 2, & 3, 3, 3, 3, 3, 3, & 0, & ! alp - 0, 0, 0 /) ! beta^i + 0, 0, 0, & ! beta^i + 0, & ! T_tt + 0, 0, 0, & ! T_ti + 0, 0, 0, 0, 0, 0 /) ! T_ij output_types(:) = CCTK_VARIABLE_REAL if (hn > 0) then @@ -231,7 +273,10 @@ subroutine qlm_interpolate (CCTK_ARGUMENTS, hn) P(qlm_dkxxy), P(qlm_dkxyy), P(qlm_dkxzy), P(qlm_dkyyy), P(qlm_dkyzy), P(qlm_dkzzy), & P(qlm_dkxxz), P(qlm_dkxyz), P(qlm_dkxzz), P(qlm_dkyyz), P(qlm_dkyzz), P(qlm_dkzzz), & P(qlm_alpha), & - P(qlm_betax), P(qlm_betay), P(qlm_betaz) /) + P(qlm_betax), P(qlm_betay), P(qlm_betaz), & + P(qlm_ttt), & + P(qlm_ttx), P(qlm_tty), P(qlm_ttz), & + P(qlm_txx), P(qlm_txy), P(qlm_txz), P(qlm_tyy), P(qlm_tyz), P(qlm_tzz) /) else outputs(:) = CCTK_NullPointer() end if @@ -245,95 +290,104 @@ subroutine qlm_interpolate (CCTK_ARGUMENTS, hn) #if 0 ! Poison the output variables -#define poison -42 - qlm_gxx = poison - qlm_gxy = poison - qlm_gxz = poison - qlm_gyy = poison - qlm_gyz = poison - qlm_gzz = poison - qlm_dgxxx = poison - qlm_dgxyx = poison - qlm_dgxzx = poison - qlm_dgyyx = poison - qlm_dgyzx = poison - qlm_dgzzx = poison - qlm_dgxxy = poison - qlm_dgxyy = poison - qlm_dgxzy = poison - qlm_dgyyy = poison - qlm_dgyzy = poison - qlm_dgzzy = poison - qlm_dgxxz = poison - qlm_dgxyz = poison - qlm_dgxzz = poison - qlm_dgyyz = poison - qlm_dgyzz = poison - qlm_dgzzz = poison - qlm_ddgxxxx = poison - qlm_ddgxyxx = poison - qlm_ddgxzxx = poison - qlm_ddgyyxx = poison - qlm_ddgyzxx = poison - qlm_ddgzzxx = poison - qlm_ddgxxxy = poison - qlm_ddgxyxy = poison - qlm_ddgxzxy = poison - qlm_ddgyyxy = poison - qlm_ddgyzxy = poison - qlm_ddgzzxy = poison - qlm_ddgxxxz = poison - qlm_ddgxyxz = poison - qlm_ddgxzxz = poison - qlm_ddgyyxz = poison - qlm_ddgyzxz = poison - qlm_ddgzzxz = poison - qlm_ddgxxyy = poison - qlm_ddgxyyy = poison - qlm_ddgxzyy = poison - qlm_ddgyyyy = poison - qlm_ddgyzyy = poison - qlm_ddgzzyy = poison - qlm_ddgxxyz = poison - qlm_ddgxyyz = poison - qlm_ddgxzyz = poison - qlm_ddgyyyz = poison - qlm_ddgyzyz = poison - qlm_ddgzzyz = poison - qlm_ddgxxzz = poison - qlm_ddgxyzz = poison - qlm_ddgxzzz = poison - qlm_ddgyyzz = poison - qlm_ddgyzzz = poison - qlm_ddgzzzz = poison - qlm_kxx = poison - qlm_kxy = poison - qlm_kxz = poison - qlm_kyy = poison - qlm_kyz = poison - qlm_kzz = poison - qlm_dkxxx = poison - qlm_dkxyx = poison - qlm_dkxzx = poison - qlm_dkyyx = poison - qlm_dkyzx = poison - qlm_dkzzx = poison - qlm_dkxxy = poison - qlm_dkxyy = poison - qlm_dkxzy = poison - qlm_dkyyy = poison - qlm_dkyzy = poison - qlm_dkzzy = poison - qlm_dkxxz = poison - qlm_dkxyz = poison - qlm_dkxzz = poison - qlm_dkyyz = poison - qlm_dkyzz = poison - qlm_dkzzz = poison - qlm_alpha = poison - qlm_betax = poison - qlm_betay = poison - qlm_betaz = poison + call poison (qlm_gxx ) + call poison (qlm_gxy ) + call poison (qlm_gxz ) + call poison (qlm_gyy ) + call poison (qlm_gyz ) + call poison (qlm_gzz ) + call poison (qlm_dgxxx ) + call poison (qlm_dgxyx ) + call poison (qlm_dgxzx ) + call poison (qlm_dgyyx ) + call poison (qlm_dgyzx ) + call poison (qlm_dgzzx ) + call poison (qlm_dgxxy ) + call poison (qlm_dgxyy ) + call poison (qlm_dgxzy ) + call poison (qlm_dgyyy ) + call poison (qlm_dgyzy ) + call poison (qlm_dgzzy ) + call poison (qlm_dgxxz ) + call poison (qlm_dgxyz ) + call poison (qlm_dgxzz ) + call poison (qlm_dgyyz ) + call poison (qlm_dgyzz ) + call poison (qlm_dgzzz ) + call poison (qlm_ddgxxxx) + call poison (qlm_ddgxyxx) + call poison (qlm_ddgxzxx) + call poison (qlm_ddgyyxx) + call poison (qlm_ddgyzxx) + call poison (qlm_ddgzzxx) + call poison (qlm_ddgxxxy) + call poison (qlm_ddgxyxy) + call poison (qlm_ddgxzxy) + call poison (qlm_ddgyyxy) + call poison (qlm_ddgyzxy) + call poison (qlm_ddgzzxy) + call poison (qlm_ddgxxxz) + call poison (qlm_ddgxyxz) + call poison (qlm_ddgxzxz) + call poison (qlm_ddgyyxz) + call poison (qlm_ddgyzxz) + call poison (qlm_ddgzzxz) + call poison (qlm_ddgxxyy) + call poison (qlm_ddgxyyy) + call poison (qlm_ddgxzyy) + call poison (qlm_ddgyyyy) + call poison (qlm_ddgyzyy) + call poison (qlm_ddgzzyy) + call poison (qlm_ddgxxyz) + call poison (qlm_ddgxyyz) + call poison (qlm_ddgxzyz) + call poison (qlm_ddgyyyz) + call poison (qlm_ddgyzyz) + call poison (qlm_ddgzzyz) + call poison (qlm_ddgxxzz) + call poison (qlm_ddgxyzz) + call poison (qlm_ddgxzzz) + call poison (qlm_ddgyyzz) + call poison (qlm_ddgyzzz) + call poison (qlm_ddgzzzz) + call poison (qlm_kxx ) + call poison (qlm_kxy ) + call poison (qlm_kxz ) + call poison (qlm_kyy ) + call poison (qlm_kyz ) + call poison (qlm_kzz ) + call poison (qlm_dkxxx ) + call poison (qlm_dkxyx ) + call poison (qlm_dkxzx ) + call poison (qlm_dkyyx ) + call poison (qlm_dkyzx ) + call poison (qlm_dkzzx ) + call poison (qlm_dkxxy ) + call poison (qlm_dkxyy ) + call poison (qlm_dkxzy ) + call poison (qlm_dkyyy ) + call poison (qlm_dkyzy ) + call poison (qlm_dkzzy ) + call poison (qlm_dkxxz ) + call poison (qlm_dkxyz ) + call poison (qlm_dkxzz ) + call poison (qlm_dkyyz ) + call poison (qlm_dkyzz ) + call poison (qlm_dkzzz ) + call poison (qlm_alpha ) + call poison (qlm_betax ) + call poison (qlm_betay ) + call poison (qlm_betaz ) + call poison (qlm_ttt ) + call poison (qlm_ttx ) + call poison (qlm_tty ) + call poison (qlm_ttz ) + call poison (qlm_txx ) + call poison (qlm_txy ) + call poison (qlm_txz ) + call poison (qlm_tyy ) + call poison (qlm_tyz ) + call poison (qlm_tzz ) #endif @@ -368,187 +422,220 @@ subroutine qlm_interpolate (CCTK_ARGUMENTS, hn) ! Unpack the variables if (hn > 0) then - call unpack (qlm_gxx , ni, nj) - call unpack (qlm_gxy , ni, nj) - call unpack (qlm_gxz , ni, nj) - call unpack (qlm_gyy , ni, nj) - call unpack (qlm_gyz , ni, nj) - call unpack (qlm_gzz , ni, nj) - call unpack (qlm_dgxxx , ni, nj) - call unpack (qlm_dgxyx , ni, nj) - call unpack (qlm_dgxzx , ni, nj) - call unpack (qlm_dgyyx , ni, nj) - call unpack (qlm_dgyzx , ni, nj) - call unpack (qlm_dgzzx , ni, nj) - call unpack (qlm_dgxxy , ni, nj) - call unpack (qlm_dgxyy , ni, nj) - call unpack (qlm_dgxzy , ni, nj) - call unpack (qlm_dgyyy , ni, nj) - call unpack (qlm_dgyzy , ni, nj) - call unpack (qlm_dgzzy , ni, nj) - call unpack (qlm_dgxxz , ni, nj) - call unpack (qlm_dgxyz , ni, nj) - call unpack (qlm_dgxzz , ni, nj) - call unpack (qlm_dgyyz , ni, nj) - call unpack (qlm_dgyzz , ni, nj) - call unpack (qlm_dgzzz , ni, nj) - call unpack (qlm_ddgxxxx , ni, nj) - call unpack (qlm_ddgxyxx , ni, nj) - call unpack (qlm_ddgxzxx , ni, nj) - call unpack (qlm_ddgyyxx , ni, nj) - call unpack (qlm_ddgyzxx , ni, nj) - call unpack (qlm_ddgzzxx , ni, nj) - call unpack (qlm_ddgxxxy , ni, nj) - call unpack (qlm_ddgxyxy , ni, nj) - call unpack (qlm_ddgxzxy , ni, nj) - call unpack (qlm_ddgyyxy , ni, nj) - call unpack (qlm_ddgyzxy , ni, nj) - call unpack (qlm_ddgzzxy , ni, nj) - call unpack (qlm_ddgxxxz , ni, nj) - call unpack (qlm_ddgxyxz , ni, nj) - call unpack (qlm_ddgxzxz , ni, nj) - call unpack (qlm_ddgyyxz , ni, nj) - call unpack (qlm_ddgyzxz , ni, nj) - call unpack (qlm_ddgzzxz , ni, nj) - call unpack (qlm_ddgxxyy , ni, nj) - call unpack (qlm_ddgxyyy , ni, nj) - call unpack (qlm_ddgxzyy , ni, nj) - call unpack (qlm_ddgyyyy , ni, nj) - call unpack (qlm_ddgyzyy , ni, nj) - call unpack (qlm_ddgzzyy , ni, nj) - call unpack (qlm_ddgxxyz , ni, nj) - call unpack (qlm_ddgxyyz , ni, nj) - call unpack (qlm_ddgxzyz , ni, nj) - call unpack (qlm_ddgyyyz , ni, nj) - call unpack (qlm_ddgyzyz , ni, nj) - call unpack (qlm_ddgzzyz , ni, nj) - call unpack (qlm_ddgxxzz , ni, nj) - call unpack (qlm_ddgxyzz , ni, nj) - call unpack (qlm_ddgxzzz , ni, nj) - call unpack (qlm_ddgyyzz , ni, nj) - call unpack (qlm_ddgyzzz , ni, nj) - call unpack (qlm_ddgzzzz , ni, nj) - call unpack (qlm_kxx , ni, nj) - call unpack (qlm_kxy , ni, nj) - call unpack (qlm_kxz , ni, nj) - call unpack (qlm_kyy , ni, nj) - call unpack (qlm_kyz , ni, nj) - call unpack (qlm_kzz , ni, nj) - call unpack (qlm_dkxxx , ni, nj) - call unpack (qlm_dkxyx , ni, nj) - call unpack (qlm_dkxzx , ni, nj) - call unpack (qlm_dkyyx , ni, nj) - call unpack (qlm_dkyzx , ni, nj) - call unpack (qlm_dkzzx , ni, nj) - call unpack (qlm_dkxxy , ni, nj) - call unpack (qlm_dkxyy , ni, nj) - call unpack (qlm_dkxzy , ni, nj) - call unpack (qlm_dkyyy , ni, nj) - call unpack (qlm_dkyzy , ni, nj) - call unpack (qlm_dkzzy , ni, nj) - call unpack (qlm_dkxxz , ni, nj) - call unpack (qlm_dkxyz , ni, nj) - call unpack (qlm_dkxzz , ni, nj) - call unpack (qlm_dkyyz , ni, nj) - call unpack (qlm_dkyzz , ni, nj) - call unpack (qlm_dkzzz , ni, nj) - call unpack (qlm_alpha , ni, nj) - call unpack (qlm_betax , ni, nj) - call unpack (qlm_betay , ni, nj) - call unpack (qlm_betaz , ni, nj) + call unpack (qlm_gxx , ni, nj) + call unpack (qlm_gxy , ni, nj) + call unpack (qlm_gxz , ni, nj) + call unpack (qlm_gyy , ni, nj) + call unpack (qlm_gyz , ni, nj) + call unpack (qlm_gzz , ni, nj) + call unpack (qlm_dgxxx , ni, nj) + call unpack (qlm_dgxyx , ni, nj) + call unpack (qlm_dgxzx , ni, nj) + call unpack (qlm_dgyyx , ni, nj) + call unpack (qlm_dgyzx , ni, nj) + call unpack (qlm_dgzzx , ni, nj) + call unpack (qlm_dgxxy , ni, nj) + call unpack (qlm_dgxyy , ni, nj) + call unpack (qlm_dgxzy , ni, nj) + call unpack (qlm_dgyyy , ni, nj) + call unpack (qlm_dgyzy , ni, nj) + call unpack (qlm_dgzzy , ni, nj) + call unpack (qlm_dgxxz , ni, nj) + call unpack (qlm_dgxyz , ni, nj) + call unpack (qlm_dgxzz , ni, nj) + call unpack (qlm_dgyyz , ni, nj) + call unpack (qlm_dgyzz , ni, nj) + call unpack (qlm_dgzzz , ni, nj) + call unpack (qlm_ddgxxxx, ni, nj) + call unpack (qlm_ddgxyxx, ni, nj) + call unpack (qlm_ddgxzxx, ni, nj) + call unpack (qlm_ddgyyxx, ni, nj) + call unpack (qlm_ddgyzxx, ni, nj) + call unpack (qlm_ddgzzxx, ni, nj) + call unpack (qlm_ddgxxxy, ni, nj) + call unpack (qlm_ddgxyxy, ni, nj) + call unpack (qlm_ddgxzxy, ni, nj) + call unpack (qlm_ddgyyxy, ni, nj) + call unpack (qlm_ddgyzxy, ni, nj) + call unpack (qlm_ddgzzxy, ni, nj) + call unpack (qlm_ddgxxxz, ni, nj) + call unpack (qlm_ddgxyxz, ni, nj) + call unpack (qlm_ddgxzxz, ni, nj) + call unpack (qlm_ddgyyxz, ni, nj) + call unpack (qlm_ddgyzxz, ni, nj) + call unpack (qlm_ddgzzxz, ni, nj) + call unpack (qlm_ddgxxyy, ni, nj) + call unpack (qlm_ddgxyyy, ni, nj) + call unpack (qlm_ddgxzyy, ni, nj) + call unpack (qlm_ddgyyyy, ni, nj) + call unpack (qlm_ddgyzyy, ni, nj) + call unpack (qlm_ddgzzyy, ni, nj) + call unpack (qlm_ddgxxyz, ni, nj) + call unpack (qlm_ddgxyyz, ni, nj) + call unpack (qlm_ddgxzyz, ni, nj) + call unpack (qlm_ddgyyyz, ni, nj) + call unpack (qlm_ddgyzyz, ni, nj) + call unpack (qlm_ddgzzyz, ni, nj) + call unpack (qlm_ddgxxzz, ni, nj) + call unpack (qlm_ddgxyzz, ni, nj) + call unpack (qlm_ddgxzzz, ni, nj) + call unpack (qlm_ddgyyzz, ni, nj) + call unpack (qlm_ddgyzzz, ni, nj) + call unpack (qlm_ddgzzzz, ni, nj) + call unpack (qlm_kxx , ni, nj) + call unpack (qlm_kxy , ni, nj) + call unpack (qlm_kxz , ni, nj) + call unpack (qlm_kyy , ni, nj) + call unpack (qlm_kyz , ni, nj) + call unpack (qlm_kzz , ni, nj) + call unpack (qlm_dkxxx , ni, nj) + call unpack (qlm_dkxyx , ni, nj) + call unpack (qlm_dkxzx , ni, nj) + call unpack (qlm_dkyyx , ni, nj) + call unpack (qlm_dkyzx , ni, nj) + call unpack (qlm_dkzzx , ni, nj) + call unpack (qlm_dkxxy , ni, nj) + call unpack (qlm_dkxyy , ni, nj) + call unpack (qlm_dkxzy , ni, nj) + call unpack (qlm_dkyyy , ni, nj) + call unpack (qlm_dkyzy , ni, nj) + call unpack (qlm_dkzzy , ni, nj) + call unpack (qlm_dkxxz , ni, nj) + call unpack (qlm_dkxyz , ni, nj) + call unpack (qlm_dkxzz , ni, nj) + call unpack (qlm_dkyyz , ni, nj) + call unpack (qlm_dkyzz , ni, nj) + call unpack (qlm_dkzzz , ni, nj) + call unpack (qlm_alpha , ni, nj) + call unpack (qlm_betax , ni, nj) + call unpack (qlm_betay , ni, nj) + call unpack (qlm_betaz , ni, nj) + if (stress_energy_state /= 0) then + call unpack (qlm_ttt , ni, nj) + call unpack (qlm_ttx , ni, nj) + call unpack (qlm_tty , ni, nj) + call unpack (qlm_ttz , ni, nj) + call unpack (qlm_txx , ni, nj) + call unpack (qlm_txy , ni, nj) + call unpack (qlm_txz , ni, nj) + call unpack (qlm_tyy , ni, nj) + call unpack (qlm_tyz , ni, nj) + call unpack (qlm_tzz , ni, nj) + else + qlm_ttt = 0 + qlm_ttx = 0 + qlm_tty = 0 + qlm_ttz = 0 + qlm_txx = 0 + qlm_txy = 0 + qlm_txz = 0 + qlm_tyy = 0 + qlm_tyz = 0 + qlm_tzz = 0 + end if #if 0 ! Check for poison - if (any(qlm_gxx == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_gxy == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_gxz == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_gyy == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_gyz == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_gzz == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_dgxxx == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_dgxyx == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_dgxzx == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_dgyyx == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_dgyzx == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_dgzzx == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_dgxxy == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_dgxyy == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_dgxzy == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_dgyyy == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_dgyzy == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_dgzzy == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_dgxxz == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_dgxyz == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_dgxzz == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_dgyyz == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_dgyzz == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_dgzzz == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_ddgxxxx == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_ddgxyxx == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_ddgxzxx == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_ddgyyxx == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_ddgyzxx == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_ddgzzxx == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_ddgxxxy == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_ddgxyxy == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_ddgxzxy == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_ddgyyxy == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_ddgyzxy == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_ddgzzxy == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_ddgxxxz == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_ddgxyxz == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_ddgxzxz == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_ddgyyxz == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_ddgyzxz == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_ddgzzxz == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_ddgxxyy == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_ddgxyyy == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_ddgxzyy == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_ddgyyyy == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_ddgyzyy == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_ddgzzyy == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_ddgxxyz == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_ddgxyyz == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_ddgxzyz == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_ddgyyyz == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_ddgyzyz == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_ddgzzyz == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_ddgxxzz == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_ddgxyzz == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_ddgxzzz == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_ddgyyzz == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_ddgyzzz == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_ddgzzzz == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_kxx == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_kxy == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_kxz == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_kyy == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_kyz == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_kzz == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_dkxxx == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_dkxyx == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_dkxzx == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_dkyyx == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_dkyzx == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_dkzzx == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_dkxxy == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_dkxyy == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_dkxzy == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_dkyyy == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_dkyzy == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_dkzzy == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_dkxxz == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_dkxyz == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_dkxzz == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_dkyyz == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_dkyzz == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_dkzzz == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_alpha == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_betax == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_betay == poison)) call CCTK_WARN (0, "poison found") - if (any(qlm_betaz == poison)) call CCTK_WARN (0, "poison found") + call poison_check (qlm_gxx , "qlm_gxx ") + call poison_check (qlm_gxy , "qlm_gxy ") + call poison_check (qlm_gxz , "qlm_gxz ") + call poison_check (qlm_gyy , "qlm_gyy ") + call poison_check (qlm_gyz , "qlm_gyz ") + call poison_check (qlm_gzz , "qlm_gzz ") + call poison_check (qlm_dgxxx , "qlm_dgxxx ") + call poison_check (qlm_dgxyx , "qlm_dgxyx ") + call poison_check (qlm_dgxzx , "qlm_dgxzx ") + call poison_check (qlm_dgyyx , "qlm_dgyyx ") + call poison_check (qlm_dgyzx , "qlm_dgyzx ") + call poison_check (qlm_dgzzx , "qlm_dgzzx ") + call poison_check (qlm_dgxxy , "qlm_dgxxy ") + call poison_check (qlm_dgxyy , "qlm_dgxyy ") + call poison_check (qlm_dgxzy , "qlm_dgxzy ") + call poison_check (qlm_dgyyy , "qlm_dgyyy ") + call poison_check (qlm_dgyzy , "qlm_dgyzy ") + call poison_check (qlm_dgzzy , "qlm_dgzzy ") + call poison_check (qlm_dgxxz , "qlm_dgxxz ") + call poison_check (qlm_dgxyz , "qlm_dgxyz ") + call poison_check (qlm_dgxzz , "qlm_dgxzz ") + call poison_check (qlm_dgyyz , "qlm_dgyyz ") + call poison_check (qlm_dgyzz , "qlm_dgyzz ") + call poison_check (qlm_dgzzz , "qlm_dgzzz ") + call poison_check (qlm_ddgxxxx, "qlm_ddgxxxx") + call poison_check (qlm_ddgxyxx, "qlm_ddgxyxx") + call poison_check (qlm_ddgxzxx, "qlm_ddgxzxx") + call poison_check (qlm_ddgyyxx, "qlm_ddgyyxx") + call poison_check (qlm_ddgyzxx, "qlm_ddgyzxx") + call poison_check (qlm_ddgzzxx, "qlm_ddgzzxx") + call poison_check (qlm_ddgxxxy, "qlm_ddgxxxy") + call poison_check (qlm_ddgxyxy, "qlm_ddgxyxy") + call poison_check (qlm_ddgxzxy, "qlm_ddgxzxy") + call poison_check (qlm_ddgyyxy, "qlm_ddgyyxy") + call poison_check (qlm_ddgyzxy, "qlm_ddgyzxy") + call poison_check (qlm_ddgzzxy, "qlm_ddgzzxy") + call poison_check (qlm_ddgxxxz, "qlm_ddgxxxz") + call poison_check (qlm_ddgxyxz, "qlm_ddgxyxz") + call poison_check (qlm_ddgxzxz, "qlm_ddgxzxz") + call poison_check (qlm_ddgyyxz, "qlm_ddgyyxz") + call poison_check (qlm_ddgyzxz, "qlm_ddgyzxz") + call poison_check (qlm_ddgzzxz, "qlm_ddgzzxz") + call poison_check (qlm_ddgxxyy, "qlm_ddgxxyy") + call poison_check (qlm_ddgxyyy, "qlm_ddgxyyy") + call poison_check (qlm_ddgxzyy, "qlm_ddgxzyy") + call poison_check (qlm_ddgyyyy, "qlm_ddgyyyy") + call poison_check (qlm_ddgyzyy, "qlm_ddgyzyy") + call poison_check (qlm_ddgzzyy, "qlm_ddgzzyy") + call poison_check (qlm_ddgxxyz, "qlm_ddgxxyz") + call poison_check (qlm_ddgxyyz, "qlm_ddgxyyz") + call poison_check (qlm_ddgxzyz, "qlm_ddgxzyz") + call poison_check (qlm_ddgyyyz, "qlm_ddgyyyz") + call poison_check (qlm_ddgyzyz, "qlm_ddgyzyz") + call poison_check (qlm_ddgzzyz, "qlm_ddgzzyz") + call poison_check (qlm_ddgxxzz, "qlm_ddgxxzz") + call poison_check (qlm_ddgxyzz, "qlm_ddgxyzz") + call poison_check (qlm_ddgxzzz, "qlm_ddgxzzz") + call poison_check (qlm_ddgyyzz, "qlm_ddgyyzz") + call poison_check (qlm_ddgyzzz, "qlm_ddgyzzz") + call poison_check (qlm_ddgzzzz, "qlm_ddgzzzz") + call poison_check (qlm_kxx , "qlm_kxx ") + call poison_check (qlm_kxy , "qlm_kxy ") + call poison_check (qlm_kxz , "qlm_kxz ") + call poison_check (qlm_kyy , "qlm_kyy ") + call poison_check (qlm_kyz , "qlm_kyz ") + call poison_check (qlm_kzz , "qlm_kzz ") + call poison_check (qlm_dkxxx , "qlm_dkxxx ") + call poison_check (qlm_dkxyx , "qlm_dkxyx ") + call poison_check (qlm_dkxzx , "qlm_dkxzx ") + call poison_check (qlm_dkyyx , "qlm_dkyyx ") + call poison_check (qlm_dkyzx , "qlm_dkyzx ") + call poison_check (qlm_dkzzx , "qlm_dkzzx ") + call poison_check (qlm_dkxxy , "qlm_dkxxy ") + call poison_check (qlm_dkxyy , "qlm_dkxyy ") + call poison_check (qlm_dkxzy , "qlm_dkxzy ") + call poison_check (qlm_dkyyy , "qlm_dkyyy ") + call poison_check (qlm_dkyzy , "qlm_dkyzy ") + call poison_check (qlm_dkzzy , "qlm_dkzzy ") + call poison_check (qlm_dkxxz , "qlm_dkxxz ") + call poison_check (qlm_dkxyz , "qlm_dkxyz ") + call poison_check (qlm_dkxzz , "qlm_dkxzz ") + call poison_check (qlm_dkyyz , "qlm_dkyyz ") + call poison_check (qlm_dkyzz , "qlm_dkyzz ") + call poison_check (qlm_dkzzz , "qlm_dkzzz ") + call poison_check (qlm_alpha , "qlm_alpha ") + call poison_check (qlm_betax , "qlm_betax ") + call poison_check (qlm_betay , "qlm_betay ") + call poison_check (qlm_betaz , "qlm_betaz ") + call poison_check (qlm_ttt , "qlm_ttt ") + call poison_check (qlm_ttx , "qlm_ttx ") + call poison_check (qlm_tty , "qlm_tty ") + call poison_check (qlm_ttz , "qlm_ttz ") + call poison_check (qlm_txx , "qlm_txx ") + call poison_check (qlm_txy , "qlm_txy ") + call poison_check (qlm_txz , "qlm_txz ") + call poison_check (qlm_tyy , "qlm_tyy ") + call poison_check (qlm_tyz , "qlm_tyz ") + call poison_check (qlm_tzz , "qlm_tzz ") #endif end if @@ -599,4 +686,25 @@ contains a = b end subroutine copy + subroutine poison (arr) + CCTK_REAL, intent(out) :: arr(:,:) + arr = poison_value + end subroutine poison + + subroutine poison_check (arr, name) + CCTK_REAL, intent(in) :: arr(:,:) + character(*), intent(in) :: name + character*1000 :: msg +!!$ integer :: i, j + if (any(arr==poison_value)) then + write (msg, '("Poison found in ",a)') trim(name) + call CCTK_WARN (CCTK_WARN_ALERT, msg) +!!$ do j=1,size(arr,2) +!!$ do i=1,size(arr,1) +!!$ print '(2i6)', i,j +!!$ end do +!!$ end do + end if + end subroutine poison_check + end subroutine qlm_interpolate diff --git a/src/qlm_multipoles.F90 b/src/qlm_multipoles.F90 index 79b75d5..3183d90 100644 --- a/src/qlm_multipoles.F90 +++ b/src/qlm_multipoles.F90 @@ -25,9 +25,9 @@ subroutine qlm_multipoles (CCTK_ARGUMENTS, hn) CCTK_REAL, parameter :: o128=one/128 CCTK_REAL :: qq(2,2), dtq, rsc - CCTK_REAL :: ee(3,2) - CCTK_REAL :: kk(3,3) - CCTK_REAL :: ll(0:3), nn(0:3), ss(0:3) +!!$ CCTK_REAL :: ee(3,2) +!!$ CCTK_REAL :: kk(3,3) +!!$ CCTK_REAL :: ll(0:3), nn(0:3), ss(0:3) CCTK_COMPLEX :: psi2 CCTK_REAL :: zz, dzz(2) CCTK_REAL :: area, mass, spin @@ -74,31 +74,31 @@ subroutine qlm_multipoles (CCTK_ARGUMENTS, hn) call calc_2det (qq, dtq) - ee(1,1:2) = deriv (qlm_x(:,:,hn), i, j, delta_space) - ee(2,1:2) = deriv (qlm_y(:,:,hn), i, j, delta_space) - ee(3,1:2) = deriv (qlm_z(:,:,hn), i, j, delta_space) - - kk(1,1) = qlm_kxx(i,j) - kk(1,2) = qlm_kxy(i,j) - kk(1,3) = qlm_kxz(i,j) - kk(2,2) = qlm_kyy(i,j) - kk(2,3) = qlm_kyz(i,j) - kk(3,3) = qlm_kzz(i,j) - kk(2,1) = kk(1,2) - kk(3,1) = kk(1,3) - kk(3,2) = kk(2,3) - - ll(0) = qlm_l0(i,j,hn) - ll(1) = qlm_l1(i,j,hn) - ll(2) = qlm_l2(i,j,hn) - ll(3) = qlm_l3(i,j,hn) - - nn(0) = qlm_n0(i,j,hn) - nn(1) = qlm_n1(i,j,hn) - nn(2) = qlm_n2(i,j,hn) - nn(3) = qlm_n3(i,j,hn) - - ss = (ll - nn) / sqrt(two) +!!$ ee(1,1:2) = deriv (qlm_x(:,:,hn), i, j, delta_space) +!!$ ee(2,1:2) = deriv (qlm_y(:,:,hn), i, j, delta_space) +!!$ ee(3,1:2) = deriv (qlm_z(:,:,hn), i, j, delta_space) +!!$ +!!$ kk(1,1) = qlm_kxx(i,j) +!!$ kk(1,2) = qlm_kxy(i,j) +!!$ kk(1,3) = qlm_kxz(i,j) +!!$ kk(2,2) = qlm_kyy(i,j) +!!$ kk(2,3) = qlm_kyz(i,j) +!!$ kk(3,3) = qlm_kzz(i,j) +!!$ kk(2,1) = kk(1,2) +!!$ kk(3,1) = kk(1,3) +!!$ kk(3,2) = kk(2,3) +!!$ +!!$ ll(0) = qlm_l0(i,j,hn) +!!$ ll(1) = qlm_l1(i,j,hn) +!!$ ll(2) = qlm_l2(i,j,hn) +!!$ ll(3) = qlm_l3(i,j,hn) +!!$ +!!$ nn(0) = qlm_n0(i,j,hn) +!!$ nn(1) = qlm_n1(i,j,hn) +!!$ nn(2) = qlm_n2(i,j,hn) +!!$ nn(3) = qlm_n3(i,j,hn) +!!$ +!!$ ss = (ll - nn) / sqrt(two) rsc = qlm_rsc(i,j,hn) psi2 = qlm_psi2(i,j,hn) @@ -157,50 +157,6 @@ subroutine qlm_multipoles (CCTK_ARGUMENTS, hn) end do end do - ! Normalise - -!!$ ! This is the normalisation for I_n and L_n -!!$ qlm_mp_m0(hn) = qlm_mp_m0(hn) / sqrt(4*pi/ 1) -!!$ qlm_mp_m1(hn) = qlm_mp_m1(hn) / sqrt(4*pi/ 3) -!!$ qlm_mp_m2(hn) = qlm_mp_m2(hn) / sqrt(4*pi/ 5) -!!$ qlm_mp_m3(hn) = qlm_mp_m3(hn) / sqrt(4*pi/ 7) -!!$ qlm_mp_m4(hn) = qlm_mp_m4(hn) / sqrt(4*pi/ 9) -!!$ qlm_mp_m5(hn) = qlm_mp_m5(hn) / sqrt(4*pi/11) -!!$ qlm_mp_m6(hn) = qlm_mp_m6(hn) / sqrt(4*pi/13) -!!$ qlm_mp_m7(hn) = qlm_mp_m7(hn) / sqrt(4*pi/15) -!!$ qlm_mp_m8(hn) = qlm_mp_m8(hn) / sqrt(4*pi/17) -!!$ -!!$ qlm_mp_j0(hn) = qlm_mp_j0(hn) / sqrt(4*pi/ 1) -!!$ qlm_mp_j1(hn) = qlm_mp_j1(hn) / sqrt(4*pi/ 3) -!!$ qlm_mp_j2(hn) = qlm_mp_j2(hn) / sqrt(4*pi/ 5) -!!$ qlm_mp_j3(hn) = qlm_mp_j3(hn) / sqrt(4*pi/ 7) -!!$ qlm_mp_j4(hn) = qlm_mp_j4(hn) / sqrt(4*pi/ 9) -!!$ qlm_mp_j5(hn) = qlm_mp_j5(hn) / sqrt(4*pi/11) -!!$ qlm_mp_j6(hn) = qlm_mp_j6(hn) / sqrt(4*pi/13) -!!$ qlm_mp_j7(hn) = qlm_mp_j7(hn) / sqrt(4*pi/15) -!!$ qlm_mp_j8(hn) = qlm_mp_j8(hn) / sqrt(4*pi/17) - - ! This is the normalisation for M_n and J_n - qlm_mp_m0(hn) = qlm_mp_m0(hn) * qlm_mass(hn) * qlm_radius(hn)**0 / (2*pi) - qlm_mp_m1(hn) = qlm_mp_m1(hn) * qlm_mass(hn) * qlm_radius(hn)**1 / (2*pi) - qlm_mp_m2(hn) = qlm_mp_m2(hn) * qlm_mass(hn) * qlm_radius(hn)**2 / (2*pi) - qlm_mp_m3(hn) = qlm_mp_m3(hn) * qlm_mass(hn) * qlm_radius(hn)**3 / (2*pi) - qlm_mp_m4(hn) = qlm_mp_m4(hn) * qlm_mass(hn) * qlm_radius(hn)**4 / (2*pi) - qlm_mp_m5(hn) = qlm_mp_m5(hn) * qlm_mass(hn) * qlm_radius(hn)**5 / (2*pi) - qlm_mp_m6(hn) = qlm_mp_m6(hn) * qlm_mass(hn) * qlm_radius(hn)**6 / (2*pi) - qlm_mp_m7(hn) = qlm_mp_m7(hn) * qlm_mass(hn) * qlm_radius(hn)**7 / (2*pi) - qlm_mp_m8(hn) = qlm_mp_m8(hn) * qlm_mass(hn) * qlm_radius(hn)**8 / (2*pi) - - qlm_mp_j0(hn) = qlm_mp_j0(hn) * qlm_radius(hn)**1 / (4*pi) - qlm_mp_j1(hn) = qlm_mp_j1(hn) * qlm_radius(hn)**2 / (4*pi) - qlm_mp_j2(hn) = qlm_mp_j2(hn) * qlm_radius(hn)**3 / (4*pi) - qlm_mp_j3(hn) = qlm_mp_j3(hn) * qlm_radius(hn)**4 / (4*pi) - qlm_mp_j4(hn) = qlm_mp_j4(hn) * qlm_radius(hn)**5 / (4*pi) - qlm_mp_j5(hn) = qlm_mp_j5(hn) * qlm_radius(hn)**6 / (4*pi) - qlm_mp_j6(hn) = qlm_mp_j6(hn) * qlm_radius(hn)**7 / (4*pi) - qlm_mp_j7(hn) = qlm_mp_j7(hn) * qlm_radius(hn)**8 / (4*pi) - qlm_mp_j8(hn) = qlm_mp_j8(hn) * qlm_radius(hn)**9 / (4*pi) - contains ! Legendre polynomials @@ -321,3 +277,64 @@ contains end function dp8 end subroutine qlm_multipoles + + + +subroutine qlm_multipoles_normalise (CCTK_ARGUMENTS, hn) + use cctk + use constants + implicit none + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_FUNCTIONS + DECLARE_CCTK_PARAMETERS + integer :: hn + + if (veryverbose/=0) then + call CCTK_INFO ("Normalising multipole moments") + end if + + ! Normalise + +!!$ ! This is the normalisation for I_n and L_n +!!$ qlm_mp_m0(hn) = qlm_mp_m0(hn) / sqrt(4*pi/ 1) +!!$ qlm_mp_m1(hn) = qlm_mp_m1(hn) / sqrt(4*pi/ 3) +!!$ qlm_mp_m2(hn) = qlm_mp_m2(hn) / sqrt(4*pi/ 5) +!!$ qlm_mp_m3(hn) = qlm_mp_m3(hn) / sqrt(4*pi/ 7) +!!$ qlm_mp_m4(hn) = qlm_mp_m4(hn) / sqrt(4*pi/ 9) +!!$ qlm_mp_m5(hn) = qlm_mp_m5(hn) / sqrt(4*pi/11) +!!$ qlm_mp_m6(hn) = qlm_mp_m6(hn) / sqrt(4*pi/13) +!!$ qlm_mp_m7(hn) = qlm_mp_m7(hn) / sqrt(4*pi/15) +!!$ qlm_mp_m8(hn) = qlm_mp_m8(hn) / sqrt(4*pi/17) +!!$ +!!$ qlm_mp_j0(hn) = qlm_mp_j0(hn) / sqrt(4*pi/ 1) +!!$ qlm_mp_j1(hn) = qlm_mp_j1(hn) / sqrt(4*pi/ 3) +!!$ qlm_mp_j2(hn) = qlm_mp_j2(hn) / sqrt(4*pi/ 5) +!!$ qlm_mp_j3(hn) = qlm_mp_j3(hn) / sqrt(4*pi/ 7) +!!$ qlm_mp_j4(hn) = qlm_mp_j4(hn) / sqrt(4*pi/ 9) +!!$ qlm_mp_j5(hn) = qlm_mp_j5(hn) / sqrt(4*pi/11) +!!$ qlm_mp_j6(hn) = qlm_mp_j6(hn) / sqrt(4*pi/13) +!!$ qlm_mp_j7(hn) = qlm_mp_j7(hn) / sqrt(4*pi/15) +!!$ qlm_mp_j8(hn) = qlm_mp_j8(hn) / sqrt(4*pi/17) + + ! This is the normalisation for M_n and J_n + qlm_mp_m0(hn) = qlm_mp_m0(hn) * qlm_mass(hn) * qlm_radius(hn)**0 / (2*pi) + qlm_mp_m1(hn) = qlm_mp_m1(hn) * qlm_mass(hn) * qlm_radius(hn)**1 / (2*pi) + qlm_mp_m2(hn) = qlm_mp_m2(hn) * qlm_mass(hn) * qlm_radius(hn)**2 / (2*pi) + qlm_mp_m3(hn) = qlm_mp_m3(hn) * qlm_mass(hn) * qlm_radius(hn)**3 / (2*pi) + qlm_mp_m4(hn) = qlm_mp_m4(hn) * qlm_mass(hn) * qlm_radius(hn)**4 / (2*pi) + qlm_mp_m5(hn) = qlm_mp_m5(hn) * qlm_mass(hn) * qlm_radius(hn)**5 / (2*pi) + qlm_mp_m6(hn) = qlm_mp_m6(hn) * qlm_mass(hn) * qlm_radius(hn)**6 / (2*pi) + qlm_mp_m7(hn) = qlm_mp_m7(hn) * qlm_mass(hn) * qlm_radius(hn)**7 / (2*pi) + qlm_mp_m8(hn) = qlm_mp_m8(hn) * qlm_mass(hn) * qlm_radius(hn)**8 / (2*pi) + + qlm_mp_j0(hn) = qlm_mp_j0(hn) * qlm_radius(hn)**1 / (4*pi) + qlm_mp_j1(hn) = qlm_mp_j1(hn) * qlm_radius(hn)**2 / (4*pi) + qlm_mp_j2(hn) = qlm_mp_j2(hn) * qlm_radius(hn)**3 / (4*pi) + qlm_mp_j3(hn) = qlm_mp_j3(hn) * qlm_radius(hn)**4 / (4*pi) + qlm_mp_j4(hn) = qlm_mp_j4(hn) * qlm_radius(hn)**5 / (4*pi) + qlm_mp_j5(hn) = qlm_mp_j5(hn) * qlm_radius(hn)**6 / (4*pi) + qlm_mp_j6(hn) = qlm_mp_j6(hn) * qlm_radius(hn)**7 / (4*pi) + qlm_mp_j7(hn) = qlm_mp_j7(hn) * qlm_radius(hn)**8 / (4*pi) + qlm_mp_j8(hn) = qlm_mp_j8(hn) * qlm_radius(hn)**9 / (4*pi) + +end subroutine qlm_multipoles_normalise diff --git a/src/qlm_tetrad.F90 b/src/qlm_tetrad.F90 index 85837dc..d309906 100644 --- a/src/qlm_tetrad.F90 +++ b/src/qlm_tetrad.F90 @@ -46,8 +46,8 @@ subroutine qlm_calc_tetrad (CCTK_ARGUMENTS, hn) CCTK_REAL :: nabla_ll(0:3,0:3), nabla_nn(0:3,0:3) CCTK_COMPLEX :: nabla_mm(0:3,0:3) - CCTK_REAL :: t0, t1, t2 - logical :: ce0, ce1, ce2 + !CCTK_REAL :: t0, t1, t2 + !logical :: ce0, ce1, ce2 CCTK_REAL :: delta_space(2) CCTK_REAL :: count, accuracy @@ -153,16 +153,17 @@ subroutine qlm_calc_tetrad (CCTK_ARGUMENTS, hn) ee(1,2) = qlm_y(i,j,hn) - qlm_origin_y(hn) ee(1,3) = qlm_z(i,j,hn) - qlm_origin_z(hn) - ee_p(1,1) = qlm_x_p(i,j,hn) - qlm_origin_x_p(hn) - ee_p(1,2) = qlm_y_p(i,j,hn) - qlm_origin_y_p(hn) - ee_p(1,3) = qlm_z_p(i,j,hn) - qlm_origin_z_p(hn) + !ee_p(1,1) = qlm_x_p(i,j,hn) - qlm_origin_x_p(hn) + !ee_p(1,2) = qlm_y_p(i,j,hn) - qlm_origin_y_p(hn) + !ee_p(1,3) = qlm_z_p(i,j,hn) - qlm_origin_z_p(hn) - ee_p_p(1,1) = qlm_x_p_p(i,j,hn) - qlm_origin_x_p_p(hn) - ee_p_p(1,2) = qlm_y_p_p(i,j,hn) - qlm_origin_y_p_p(hn) - ee_p_p(1,3) = qlm_z_p_p(i,j,hn) - qlm_origin_z_p_p(hn) + !ee_p_p(1,1) = qlm_x_p_p(i,j,hn) - qlm_origin_x_p_p(hn) + !ee_p_p(1,2) = qlm_y_p_p(i,j,hn) - qlm_origin_y_p_p(hn) + !ee_p_p(1,3) = qlm_z_p_p(i,j,hn) - qlm_origin_z_p_p(hn) dee(1,0,:) = 0 - dee(1,1:3,0) = timederiv (ee(1,1:3), ee_p(1,1:3), ee_p_p(1,1:3), t0,t1,t2, ce0,ce1,ce2) + !dee(1,1:3,0) = timederiv (ee(1,1:3), ee_p(1,1:3), ee_p_p(1,1:3), t0,t1,t2, ce0,ce1,ce2) + dee(1,1:3,0) = 0 dee_spher(1,:,1) = 0 ! this is a choice dee_spher(1,1,2:3) = deriv (qlm_x(:,:,hn), i,j, delta_space) dee_spher(1,2,2:3) = deriv (qlm_y(:,:,hn), i,j, delta_space) @@ -184,7 +185,8 @@ subroutine qlm_calc_tetrad (CCTK_ARGUMENTS, hn) ee_p_p(2:3,3) = deriv (qlm_z_p_p(:,:,hn), i,j, delta_space) dee(2:3,0,:) = 0 - dee(2:3,1:3,0) = timederiv (ee(2:3,1:3), ee_p(2:3,1:3), ee_p_p(2:3,1:3), t0,t1,t2, ce0,ce1,ce2) + !dee(2:3,1:3,0) = timederiv (ee(2:3,1:3), ee_p(2:3,1:3), ee_p_p(2:3,1:3), t0,t1,t2, ce0,ce1,ce2) + dee(2:3,1:3,0) = 0 dee_spher(2:3,:,1) = 0 ! this is a choice dee_spher(2:3,1,2:3) = deriv2 (qlm_x(:,:,hn), i,j, delta_space) dee_spher(2:3,2,2:3) = deriv2 (qlm_y(:,:,hn), i,j, delta_space) diff --git a/src/qlm_variables.F90 b/src/qlm_variables.F90 index d21842c..9343666 100644 --- a/src/qlm_variables.F90 +++ b/src/qlm_variables.F90 @@ -21,7 +21,10 @@ module qlm_variables qlm_dkxxy, qlm_dkxyy, qlm_dkxzy, qlm_dkyyy, qlm_dkyzy, qlm_dkzzy, & qlm_dkxxz, qlm_dkxyz, qlm_dkxzz, qlm_dkyyz, qlm_dkyzz, qlm_dkzzz, & qlm_alpha, & - qlm_betax, qlm_betay, qlm_betaz + qlm_betax, qlm_betay, qlm_betaz, & + qlm_ttt, & + qlm_ttx, qlm_tty, qlm_ttz, & + qlm_txx, qlm_txy, qlm_txz, qlm_tyy, qlm_tyz, qlm_tzz type tetrad_derivs ! nabla_ll(a,b) = D_b l_a @@ -34,8 +37,11 @@ module qlm_variables contains subroutine allocate_variables(ntheta, nphi) - integer, intent(in) :: ntheta, nphi - + integer, intent(in) :: ntheta, nphi + CCTK_REAL, parameter :: zero = 0 + integer, parameter :: rk = kind(zero) + type(tetrad_derivs) :: tetrad_derivs_nan + allocate(qlm_gxx(ntheta,nphi)) allocate(qlm_gxy(ntheta,nphi)) allocate(qlm_gxz(ntheta,nphi)) @@ -124,6 +130,16 @@ contains allocate(qlm_betax(ntheta,nphi)) allocate(qlm_betay(ntheta,nphi)) allocate(qlm_betaz(ntheta,nphi)) + allocate(qlm_ttt(ntheta,nphi)) + allocate(qlm_ttx(ntheta,nphi)) + allocate(qlm_tty(ntheta,nphi)) + allocate(qlm_ttz(ntheta,nphi)) + allocate(qlm_txx(ntheta,nphi)) + allocate(qlm_txy(ntheta,nphi)) + allocate(qlm_txz(ntheta,nphi)) + allocate(qlm_tyy(ntheta,nphi)) + allocate(qlm_tyz(ntheta,nphi)) + allocate(qlm_tzz(ntheta,nphi)) allocate(qlm_tetrad_derivs(ntheta,nphi)) @@ -215,6 +231,21 @@ contains qlm_betax = TAT_nan() qlm_betay = TAT_nan() qlm_betaz = TAT_nan() + qlm_ttt = TAT_nan() + qlm_ttx = TAT_nan() + qlm_tty = TAT_nan() + qlm_ttz = TAT_nan() + qlm_txx = TAT_nan() + qlm_txy = TAT_nan() + qlm_txz = TAT_nan() + qlm_tyy = TAT_nan() + qlm_tyz = TAT_nan() + qlm_tzz = TAT_nan() + + tetrad_derivs_nan%nabla_ll = TAT_nan() + tetrad_derivs_nan%nabla_nn = TAT_nan() + tetrad_derivs_nan%nabla_mm = cmplx(TAT_nan(),TAT_nan(),rk) + qlm_tetrad_derivs = tetrad_derivs_nan end subroutine allocate_variables subroutine deallocate_variables @@ -306,6 +337,16 @@ contains deallocate(qlm_betax) deallocate(qlm_betay) deallocate(qlm_betaz) + deallocate(qlm_ttt) + deallocate(qlm_ttx) + deallocate(qlm_tty) + deallocate(qlm_ttz) + deallocate(qlm_txx) + deallocate(qlm_txy) + deallocate(qlm_txz) + deallocate(qlm_tyy) + deallocate(qlm_tyz) + deallocate(qlm_tzz) deallocate(qlm_tetrad_derivs) end subroutine deallocate_variables diff --git a/src/qlm_weyl_scalars.F90 b/src/qlm_weyl_scalars.F90 index d79751c..ade2e31 100644 --- a/src/qlm_weyl_scalars.F90 +++ b/src/qlm_weyl_scalars.F90 @@ -24,6 +24,9 @@ subroutine qlm_calc_weyl_scalars (CCTK_ARGUMENTS, hn) CCTK_REAL, parameter :: two=2, four=4 CCTK_REAL :: gg(3,3), dgg(3,3,3), ddgg(3,3,3,3), gg_dot(3,3), gg_dot2(3,3), dgg_dot(3,3,3) + CCTK_REAL :: kk(3,3), dkk(3,3,3), kk_dot(3,3) + CCTK_REAL :: tt(3,3) + CCTK_REAL :: dtg, gu(3,3), dgu(3,3,3), gamma(3,3,3), dgamma(3,3,3,3), ri(3,3) CCTK_REAL :: g4(0:3,0:3), dg4(0:3,0:3,0:3), ddg4(0:3,0:3,0:3,0:3) CCTK_REAL :: gu4(0:3,0:3), dgu4(0:3,0:3,0:3) CCTK_REAL :: gamma4(0:3,0:3,0:3), dgamma4(0:3,0:3,0:3,0:3) @@ -122,6 +125,48 @@ subroutine qlm_calc_weyl_scalars (CCTK_ARGUMENTS, hn) ddgg(:,:,3,1) = ddgg(:,:,1,3) ddgg(:,:,3,2) = ddgg(:,:,2,3) + kk(1,1) = qlm_kxx(i,j) + kk(1,2) = qlm_kxy(i,j) + kk(1,3) = qlm_kxz(i,j) + kk(2,2) = qlm_kyy(i,j) + kk(2,3) = qlm_kyz(i,j) + kk(3,3) = qlm_kzz(i,j) + kk(2,1) = kk(1,2) + kk(3,1) = kk(1,3) + kk(3,2) = kk(2,3) + + dkk(1,1,1) = qlm_dkxxx(i,j) + dkk(1,2,1) = qlm_dkxyx(i,j) + dkk(1,3,1) = qlm_dkxzx(i,j) + dkk(2,2,1) = qlm_dkyyx(i,j) + dkk(2,3,1) = qlm_dkyzx(i,j) + dkk(3,3,1) = qlm_dkzzx(i,j) + dkk(1,1,2) = qlm_dkxxy(i,j) + dkk(1,2,2) = qlm_dkxyy(i,j) + dkk(1,3,2) = qlm_dkxzy(i,j) + dkk(2,2,2) = qlm_dkyyy(i,j) + dkk(2,3,2) = qlm_dkyzy(i,j) + dkk(3,3,2) = qlm_dkzzy(i,j) + dkk(1,1,3) = qlm_dkxxz(i,j) + dkk(1,2,3) = qlm_dkxyz(i,j) + dkk(1,3,3) = qlm_dkxzz(i,j) + dkk(2,2,3) = qlm_dkyyz(i,j) + dkk(2,3,3) = qlm_dkyzz(i,j) + dkk(3,3,3) = qlm_dkzzz(i,j) + dkk(2,1,:) = dkk(1,2,:) + dkk(3,1,:) = dkk(1,3,:) + dkk(3,2,:) = dkk(2,3,:) + + tt(1,1) = qlm_txx(i,j) + tt(1,2) = qlm_txy(i,j) + tt(1,3) = qlm_txz(i,j) + tt(2,2) = qlm_tyy(i,j) + tt(2,3) = qlm_tyz(i,j) + tt(3,3) = qlm_tzz(i,j) + tt(2,1) = tt(1,2) + tt(3,1) = tt(1,3) + tt(3,2) = tt(2,3) + ll(0) = qlm_l0(i,j,hn) ll(1) = qlm_l1(i,j,hn) ll(2) = qlm_l2(i,j,hn) @@ -140,6 +185,18 @@ subroutine qlm_calc_weyl_scalars (CCTK_ARGUMENTS, hn) ! Calculate 4-metric + call calc_det (gg, dtg) + call calc_inv (gg, dtg, gu) + call calc_invderiv (gu, dgg, dgu) + call calc_connections (gu, dgg, gamma) + call calc_connectionderivs (gu, dgg, dgu, ddgg, dgamma) + call calc_ricci (gamma, dgamma, ri) + + call calc_3metricdot_simple (kk, gg_dot) + call calc_3metricderivdot_simple (dkk, dgg_dot) + call calc_extcurvdot_simple (gg,gu,ri, kk, tt, kk_dot) + call calc_3metricdot2_simple (kk_dot, gg_dot2) + call calc_4metricderivs2_simple (gg, dgg, & ddgg, gg_dot, gg_dot2, dgg_dot, g4,dg4,ddg4) call calc_4inv (g4, gu4) -- cgit v1.2.3