[Commits] [svn:einsteintoolkit] NullNews/trunk/ (Rev. 15)
reisswig at tapir.caltech.edu
reisswig at tapir.caltech.edu
Wed Jun 5 02:58:20 CDT 2013
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
More information about the Commits
mailing list