User: reisswig Date: 2013/06/05 02:58 AM

Modified: /trunk/ interface.ccl, param.ccl, schedule.ccl /trunk/src/ NullNews_GetNews.F90, NullNews_SphHarmDecomp.F90, NullNews_interp_to_constant_uBondi.cc, set_pointers.h

Log: Implementation of lin. strain computation.

File Changes:

Directory: /trunk/src/ ======================

File [modified]: NullNews_GetNews.F90 Delta lines: +18 -2 =================================================================== --- trunk/src/NullNews_GetNews.F90 2013-05-27 22:33:26 UTC (rev 14) +++ trunk/src/NullNews_GetNews.F90 2013-06-05 07:58:20 UTC (rev 15) @@ -19,10 +19,10 @@

CCTK_COMPLEX, dimension (:,:,:), allocatable, save::& J, J_u, J_l, J_l_u, cB, U, comega, eth_J,& - beth_J, eth_U, beth_U, U_u, U_l, U_l_l + beth_J, eth_U, beth_U, U_u, U_l, U_l_l, etheth_u0

CCTK_REAL, dimension (:,:,:), allocatable, save::& - beta, omega, beta_u + beta, omega, beta_u, u0

CCTK_REAL, dimension(3) :: BUF logical, save :: FirstTime = .True. @@ -62,6 +62,12 @@ eth_J=0; beth_J=0; eth_U=0; beth_U=0; beta=0; omega=0 beta_u = 0; U_u=0; U_l=0; U_l_l=0

+ if (compute_lin_strain .ne. 0) then + allocate(u0(null_lsh(1), null_lsh(2), 2), etheth_u0(null_lsh(1), null_lsh(2), 2)) + u0=0 + etheth_u0=0 + endif + if (minval(abs(zeta - dcmplx(1,1))) < 1.0d-10) then ThisProc = .true. Tarr = minloc(abs(stereo_q(:,1)-1.)) @@ -196,6 +202,16 @@ end do end if

+ if (compute_lin_strain .ne. 0) then + u0 = uBondi + call NullInterp_d2 (etheth_u0(:,:,1), dcmplx(u0(:,:,1)), 0, 1, 1) + call NullInterp_d2 (etheth_u0(:,:,2), dcmplx(u0(:,:,2)), 0, 1, 1) + call NullInterp_cinterp(cctkGH, tmp_cgfn, tmp_cgfs, etheth_u0(:,:,1), etheth_u0(:,:,2), 1) + call NullNews_ResetInactive(null_lsh, etheth_u0) + linStrain(:,:,1) = J_l(:,:,1) + tmp_cgfn + linStrain(:,:,2) = J_l(:,:,2) + tmp_cgfs + endif + !! at this point its counterproductive to change spin-basis ! ! La definicion es z=\tan(\theta/2) e^{i\phi}, asi que aca le ! ! sacamos el z/zb que le pusimos en Spheroid_RJ.

File [modified]: NullNews_SphHarmDecomp.F90 Delta lines: +15 -0 =================================================================== --- trunk/src/NullNews_SphHarmDecomp.F90 2013-05-27 22:33:26 UTC (rev 14) +++ trunk/src/NullNews_SphHarmDecomp.F90 2013-06-05 07:58:20 UTC (rev 15) @@ -31,6 +31,7 @@ CCTK_COMPLEX, dimension(0:lmax, -lmax:lmax) :: betanSphCoeff CCTK_COMPLEX, dimension(0:lmax, -lmax:lmax) :: WnSphCoeff CCTK_COMPLEX, dimension(0:lmax, -lmax:lmax) :: uBondiSphCoeff + CCTK_COMPLEX, dimension(2:lmax, -lmax:lmax) :: linStrainSphCoeff CCTK_REAL :: NewsNorm CCTK_INT :: spin, lmode, nchar character(len=256) :: io_out_dir, filestring_re, filestring_im, filestring, tmp @@ -83,6 +84,11 @@ zeta, NewsB, NewsBSphCoeff) call SpinDecompCoefs(cctkGH, null_lsh(1), null_lsh(2), NewsSpinWeight, & zeta, Psi4, Psi4SphCoeff) + + if (compute_lin_strain .ne. 0) then + call SpinDecompCoefs(cctkGH, null_lsh(1), null_lsh(2), NewsSpinWeight, & + zeta, linStrain, linStrainSphCoeff) + endif

else time = constant_uBondi @@ -93,6 +99,11 @@ zeta, NewsB_uBondi, NewsBSphCoeff) call SpinDecompCoefs(cctkGH, null_lsh(1), null_lsh(2), NewsSpinWeight, & zeta, Psi4_uBondi, Psi4SphCoeff) + + if (compute_lin_strain .ne. 0) then + call SpinDecompCoefs(cctkGH, null_lsh(1), null_lsh(2), NewsSpinWeight, & + zeta, linStrain_uBondi, linStrainSphCoeff) + endif endif

if (debug_output == 1) then @@ -138,6 +149,10 @@ call NullDecomp_WriteCoefFile('NewsB_scri', truncate, 2, lmax, NewsBSphCoeff, time) call NullDecomp_WriteCoefFile('Psi4_scri', truncate, 2, lmax, Psi4SphCoeff, time)

+ if (compute_lin_strain .ne. 0) then + call NullDecomp_WriteCoefFile('linStrain_scri', truncate, 2, lmax, linStrainSphCoeff, time) + endif + if (debug_output == 1) then

call NullDecomp_WriteCoefFile('J_scri', truncate, 2, lmax, JnSphCoeff, time)

File [modified]: NullNews_interp_to_constant_uBondi.cc Delta lines: +15 -0 =================================================================== --- trunk/src/NullNews_interp_to_constant_uBondi.cc 2013-05-27 22:33:26 UTC (rev 14) +++ trunk/src/NullNews_interp_to_constant_uBondi.cc 2013-06-05 07:58:20 UTC (rev 15) @@ -159,6 +159,16 @@ NewsBP[tl-2][ij].imag(), NewsBP[tl-1][ij].imag(), NewsBP[tl][ij].imag(), NewsBP[tl+1][ij].imag(), NewsBP[tl+2][ij].imag(), uBondiP[tl-2][ij], uBondiP[tl-1][ij], uBondiP[tl][ij], uBondiP[tl+1][ij], uBondiP[tl+2][ij]); NewsB_uBondi[ij] = CCTK_Cmplx(real, imag); + + if (compute_lin_strain) { + real = interpolant(*constant_uBondi, + linStrainP[tl-2][ij].real(), linStrainP[tl-1][ij].real(), linStrainP[tl][ij].real(), linStrainP[tl+1][ij].real(), linStrainP[tl+2][ij].real(), + uBondiP[tl-2][ij], uBondiP[tl-1][ij], uBondiP[tl][ij], uBondiP[tl+1][ij], uBondiP[tl+2][ij]); + imag = interpolant(*constant_uBondi, + linStrainP[tl-2][ij].imag(), linStrainP[tl-1][ij].imag(), linStrainP[tl][ij].imag(), linStrainP[tl+1][ij].imag(), linStrainP[tl+2][ij].imag(), + uBondiP[tl-2][ij], uBondiP[tl-1][ij], uBondiP[tl][ij], uBondiP[tl+1][ij], uBondiP[tl+2][ij]); + linStrain_uBondi[ij] = CCTK_Cmplx(real, imag); + } }

} @@ -187,6 +197,9 @@ NewsP[tl][ij] = CCTK_Cmplx(0,0);

NewsBP[tl][ij] = CCTK_Cmplx(0,0); + + if (compute_lin_strain) + linStrainP[tl][ij] = CCTK_Cmplx(0,0); }

*constant_uBondi = uBondiP[2][0]; @@ -215,6 +228,8 @@ Psi4P[tl+1][ij] = Psi4P[tl][ij]; NewsP[tl+1][ij] = NewsP[tl][ij]; NewsBP[tl+1][ij] = NewsBP[tl][ij]; + if (compute_lin_strain) + linStrainP[tl+1][ij] = linStrainP[tl][ij]; } }

File [modified]: set_pointers.h Delta lines: +15 -0 =================================================================== --- trunk/src/set_pointers.h 2013-05-27 22:33:26 UTC (rev 14) +++ trunk/src/set_pointers.h 2013-06-05 07:58:20 UTC (rev 15) @@ -4,6 +4,7 @@ vector<CCTK_COMPLEX*> Psi4P(max_timelevels, (CCTK_COMPLEX*)NULL); vector<CCTK_COMPLEX*> NewsP(max_timelevels, (CCTK_COMPLEX*)NULL); vector<CCTK_COMPLEX*> NewsBP(max_timelevels, (CCTK_COMPLEX*)NULL); + vector<CCTK_COMPLEX*> linStrainP(max_timelevels, (CCTK_COMPLEX*)NULL);

// Get current timelevel... // We assume here that the two vars uBondi[0] amd uBondi[1] are right next to each @@ -12,15 +13,22 @@ const int varindex_Psi4 = CCTK_VarIndex("NullNews::Psi4[0]"); const int varindex_News = CCTK_VarIndex("NullNews::News[0]"); const int varindex_NewsB = CCTK_VarIndex("NullNews::NewsB[0]"); + const int varindex_linStrain = CCTK_VarIndex("NullNews::linStrain[0]");

if (varindex_uBondi < 0 || varindex_Psi4 < 0 || varindex_News < 0 || varindex_NewsB < 0) CCTK_WARN(0, "Error getting variable index!");

+ if (compute_lin_strain && varindex_linStrain < 0) + CCTK_WARN(0, "Error getting variable index!"); + uBondiP[0] = (CCTK_REAL*) CCTK_VarDataPtrI(cctkGH, 0, varindex_uBondi); Psi4P[0] = (CCTK_COMPLEX*) CCTK_VarDataPtrI(cctkGH, 0, varindex_Psi4); NewsP[0] = (CCTK_COMPLEX*) CCTK_VarDataPtrI(cctkGH, 0, varindex_News); NewsBP[0] = (CCTK_COMPLEX*) CCTK_VarDataPtrI(cctkGH, 0, varindex_NewsB);

+ if (compute_lin_strain) + linStrainP[0] = (CCTK_COMPLEX*) CCTK_VarDataPtrI(cctkGH, 0, varindex_linStrain); + // Get past timelevels... // We assume here that the two vars uBondi[0] amd uBondi[1] are right next to each // other in memory, i.e. the pointer to uBondi[0] can be used to access uBondi[1] as well! @@ -28,10 +36,14 @@ const int varindex_Psi4_past = CCTK_VarIndex("NullNews::Psi4_past[0]"); const int varindex_News_past = CCTK_VarIndex("NullNews::News_past[0]"); const int varindex_NewsB_past = CCTK_VarIndex("NullNews::NewsB_past[0]"); + const int varindex_linStrain_past = CCTK_VarIndex("NullNews::linStrain_past[0]");

if (varindex_uBondi_past < 0 || varindex_Psi4_past < 0 || varindex_News_past < 0 || varindex_NewsB_past < 0) CCTK_WARN(0, "Error getting variable index!");

+ if (compute_lin_strain && varindex_linStrain_past < 0) + CCTK_WARN(0, "Error getting variable index!"); + const int offset = 2*null_lsh[0]*null_lsh[1];

for (int tl=1; tl < max_timelevels; ++tl) @@ -40,4 +52,7 @@ Psi4P[tl] = &((CCTK_COMPLEX*) CCTK_VarDataPtrI(cctkGH, 0, varindex_Psi4_past))[offset*(tl-1)]; NewsP[tl] = &((CCTK_COMPLEX*) CCTK_VarDataPtrI(cctkGH, 0, varindex_News_past))[offset*(tl-1)]; NewsBP[tl] = &((CCTK_COMPLEX*) CCTK_VarDataPtrI(cctkGH, 0, varindex_NewsB_past))[offset*(tl-1)]; + + if (compute_lin_strain) + linStrainP[tl] = &((CCTK_COMPLEX*) CCTK_VarDataPtrI(cctkGH, 0, varindex_linStrain_past))[offset*(tl-1)]; }

Directory: /trunk/ ==================

File [modified]: interface.ccl Delta lines: +19 -0 =================================================================== --- trunk/interface.ccl 2013-05-27 22:33:26 UTC (rev 14) +++ trunk/interface.ccl 2013-06-05 07:58:20 UTC (rev 15) @@ -49,7 +49,20 @@ "Psi4 at constant Bondi time"

+cctk_complex linStrain[2] Type=ARRAY dim=2 size=\ +(NullGrid::N_ang_pts_inside_eq+2*(NullGrid::N_ang_ev_outside_eq+NullGrid::N_ang_stencil_size)),\ +(NullGrid::N_ang_pts_inside_eq+2*(NullGrid::N_ang_ev_outside_eq+NullGrid::N_ang_stencil_size)) \ + GHOSTSIZE=NullGrid::N_ang_ghost_pts,NullGrid::N_ang_ghost_pts \ +"linearized strain h"

+cctk_complex linStrain_uBondi[2] Type=ARRAY dim=2 size=\ +(NullGrid::N_ang_pts_inside_eq+2*(NullGrid::N_ang_ev_outside_eq+NullGrid::N_ang_stencil_size)),\ +(NullGrid::N_ang_pts_inside_eq+2*(NullGrid::N_ang_ev_outside_eq+NullGrid::N_ang_stencil_size)) \ + GHOSTSIZE=NullGrid::N_ang_ghost_pts,NullGrid::N_ang_ghost_pts \ +"linearized strain h at constant Bondi time" + + + cctk_real uBondi[2] timelevels=1 Type=ARRAY dim=2 size=\ (NullGrid::N_ang_pts_inside_eq+2*(NullGrid::N_ang_ev_outside_eq+NullGrid::N_ang_stencil_size)),\ (NullGrid::N_ang_pts_inside_eq+2*(NullGrid::N_ang_ev_outside_eq+NullGrid::N_ang_stencil_size)) \ @@ -81,8 +94,14 @@ GHOSTSIZE=NullGrid::N_ang_ghost_pts,NullGrid::N_ang_ghost_pts \ "past levels of Bondi time"

+cctk_complex linStrain_past[2*max_timelevels] Type=ARRAY timelevels=1 dim=2 size=\ +(NullGrid::N_ang_pts_inside_eq+2*(NullGrid::N_ang_ev_outside_eq+NullGrid::N_ang_stencil_size)),\ +(NullGrid::N_ang_pts_inside_eq+2*(NullGrid::N_ang_ev_outside_eq+NullGrid::N_ang_stencil_size)) \ + GHOSTSIZE=NullGrid::N_ang_ghost_pts,NullGrid::N_ang_ghost_pts \ +"past levels of lin. strain"

+ cctk_real constant_uBondi TYPE=SCALAR "constant Bondi-time to which we interpolate"

cctk_real time_of_news TYPE=SCALAR "time of news"

File [modified]: param.ccl Delta lines: +5 -0 =================================================================== --- trunk/param.ccl 2013-05-27 22:33:26 UTC (rev 14) +++ trunk/param.ccl 2013-06-05 07:58:20 UTC (rev 15) @@ -105,3 +105,8 @@ *:* :: "must be between -l and +l" } 0

+BOOLEAN compute_lin_strain "Compute linearized strain?" +{ +} no + +

File [modified]: schedule.ccl Delta lines: +11 -0 =================================================================== --- trunk/schedule.ccl 2013-05-27 22:33:26 UTC (rev 14) +++ trunk/schedule.ccl 2013-06-05 07:58:20 UTC (rev 15) @@ -14,11 +14,22 @@ storage: constant_uBondi storage: uBondi_past News_past NewsB_past Psi4_past

+ if (compute_lin_strain) + { + storage: linStrain_uBondi + storage: linStrain_past + } }

storage: uBondi[1] NewsB[1] News[1] Psi4[1]

+if (compute_lin_strain) +{ + storage: linStrain[1] +}

+ + schedule NullNews_Init at PostInitial after NullEvol_Initial { LANG: Fortran