User: zetienne Date: 2013/10/15 03:35 PM

Modified: /trunk/src/ Lorenz_psi6phi_rhs__add_gauge_terms_to_A_i_rhs.C, compute_tau_rhs_extrinsic_curvature_terms_and_TUPmunu.C, driver_conserv_to_prims.C, driver_evaluate_MHD_rhs.C, driver_prims_to_conserv.C, evaluate_MHD_rhs_headers.h, harm_primitives_headers.h, harm_primitives_lowlevel.C, harm_primitives_recompute_conservs_tau_floor.C

Log: (Ignore prev commit.) Fixed variable names (gij -> gtij, Aij,trK -> kij) for compatibility with ADMBase. Next step: fix gupij & write ADM->BSSN translation routine

File Changes:

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

File [modified]: Lorenz_psi6phi_rhs__add_gauge_terms_to_A_i_rhs.C Delta lines: +1 -0 =================================================================== --- trunk/src/Lorenz_psi6phi_rhs__add_gauge_terms_to_A_i_rhs.C 2013-10-15 20:31:57 UTC (rev 11) +++ trunk/src/Lorenz_psi6phi_rhs__add_gauge_terms_to_A_i_rhs.C 2013-10-15 20:35:06 UTC (rev 12) @@ -151,6 +151,7 @@ shiftz_iphjphkph[index]=shiftz_iphjphkphL; }

+ // This loop requires two additional ghostzones in every direction. Hence the following loop definition: #pragma omp parallel for for(int k=3;k<ext[2]-3;k++) for(int j=3;j<ext[1]-3;j++) for(int i=3;i<ext[0]-3;i++) { int index = CCTK_GFINDEX3D(cctkGH,i,j,k);

File [modified]: compute_tau_rhs_extrinsic_curvature_terms_and_TUPmunu.C Delta lines: +16 -25 =================================================================== --- trunk/src/compute_tau_rhs_extrinsic_curvature_terms_and_TUPmunu.C 2013-10-15 20:31:57 UTC (rev 11) +++ trunk/src/compute_tau_rhs_extrinsic_curvature_terms_and_TUPmunu.C 2013-10-15 20:35:06 UTC (rev 12) @@ -50,7 +50,7 @@ void compute_tau_rhs_extrinsic_curvature_terms_and_TUPmunu (const cGH *cctkGH,int *ext,double *dX,double **metric,gf_and_gz_struct *prims,double **TUPmunu,double *u0,eos_struct &eos, double *gupxy,double *gupxz,double *gupyz, - double *Axx,double *Axy,double *Axz,double *Ayy,double *Ayz,double *Azz,double *trK, + double *kxx,double *kxy,double *kxz,double *kyy,double *kyz,double *kzz, double *tau_rhs) {

// These loop extents must be consistent with add_fluxes_and_source_terms_to_hydro_rhss(), since we use TUPmunu there as well. @@ -69,13 +69,12 @@ // forces us to read in only once, storing in local variables // U{p,m}{1,2,3}, METRIC{p,m}{1,2}, U, and METRIC.

- double trKL = trK[index]; - double AxxL = Axx[index]; - double AxyL = Axy[index]; - double AxzL = Axz[index]; - double AyyL = Ayy[index]; - double AyzL = Ayz[index]; - double AzzL = Azz[index]; + double KxxL = kxx[index]; + double KxyL = kxy[index]; + double KxzL = kxz[index]; + double kyyL = kyy[index]; + double kyzL = kyz[index]; + double kzzL = kzz[index];

//----------------------------------------------------------------------------- // Compute T^{\mu \nu} @@ -140,26 +139,18 @@ // First, compute extrinsic curvature terms for tau_rhs: // \alpha \sqrt{\gamma} (T^{00} \beta^i \beta ^j + 2T^{0i} \beta^j + T^{ij})*K_{ij} double alpha_sqrtgamma = 2.0*half_alpha_sqrtgamma; -#define ONETHIRD 0.33333333333333333333333333333333333 - double Psi4 = 1.0/Psim4; - double Kxx = Psi4*(AxxL+ONETHIRD*METRIC[GXX]*trKL); - double Kxy = Psi4*(AxyL+ONETHIRD*METRIC[GXY]*trKL); - double Kxz = Psi4*(AxzL+ONETHIRD*METRIC[GXZ]*trKL); - double Kyy = Psi4*(AyyL+ONETHIRD*METRIC[GYY]*trKL); - double Kyz = Psi4*(AyzL+ONETHIRD*METRIC[GYZ]*trKL); - double Kzz = Psi4*(AzzL+ONETHIRD*METRIC[GZZ]*trKL);

// \alpha \sqrt{\gamma} (T^{00} \beta^i \beta ^j + 2T^{0i} \beta^j + T^{ij})*K_{ij} double tau_rhs_extrinsic_curvature_terms = alpha_sqrtgamma* ( - TUP[0][0] * (SQR(METRIC[SHIFTX])*Kxx + SQR(METRIC[SHIFTY])*Kyy + SQR(METRIC[SHIFTZ])*Kzz + - 2.0*(METRIC[SHIFTX]*METRIC[SHIFTY]*Kxy + METRIC[SHIFTX]*METRIC[SHIFTZ]*Kxz + METRIC[SHIFTY]*METRIC[SHIFTZ]*Kyz) ) + - 2.0*(TUP[0][1]*METRIC[SHIFTX]*Kxx + TUP[0][2]*METRIC[SHIFTY]*Kyy + TUP[0][3]*METRIC[SHIFTZ]*Kzz + - (TUP[0][1]*METRIC[SHIFTY] + TUP[0][2]*METRIC[SHIFTX])*Kxy + - (TUP[0][1]*METRIC[SHIFTZ] + TUP[0][3]*METRIC[SHIFTX])*Kxz + - (TUP[0][2]*METRIC[SHIFTZ] + TUP[0][3]*METRIC[SHIFTY])*Kyz ) + - TUP[1][1]*Kxx + TUP[2][2]*Kyy + TUP[3][3]*Kzz + - 2.0*(TUP[1][2]*Kxy + TUP[1][3]*Kxz + TUP[2][3]*Kyz) ); + TUP[0][0] * (SQR(METRIC[SHIFTX])*KxxL + SQR(METRIC[SHIFTY])*KyyL + SQR(METRIC[SHIFTZ])*KzzL + + 2.0*(METRIC[SHIFTX]*METRIC[SHIFTY]*KxyL + METRIC[SHIFTX]*METRIC[SHIFTZ]*KxzL + METRIC[SHIFTY]*METRIC[SHIFTZ]*KyzL) ) + + 2.0*(TUP[0][1]*METRIC[SHIFTX]*KxxL + TUP[0][2]*METRIC[SHIFTY]*KyyL + TUP[0][3]*METRIC[SHIFTZ]*KzzL + + (TUP[0][1]*METRIC[SHIFTY] + TUP[0][2]*METRIC[SHIFTX])*KxyL + + (TUP[0][1]*METRIC[SHIFTZ] + TUP[0][3]*METRIC[SHIFTX])*KxzL + + (TUP[0][2]*METRIC[SHIFTZ] + TUP[0][3]*METRIC[SHIFTY])*KyzL ) + + TUP[1][1]*KxxL + TUP[2][2]*KyyL + TUP[3][3]*KzzL + + 2.0*(TUP[1][2]*KxyL + TUP[1][3]*KxzL + TUP[2][3]*KyzL) );

tau_rhs[index] = tau_rhs_extrinsic_curvature_terms;

@@ -201,7 +192,7 @@ smallb[SMALLBT] = alpha_sqrt_4pi_bt * ONE_OVER_LAPSE_SQRT_4PI;

// b^2 = g_{\mu \nu} b^{\mu} b^{\nu} - // = gtt bt^2 + gxx bx^2 + gyy by^2 + gzz bz^2 + 2 (gtx bt bx + gty bt by + gtz bt bz + gxy bx by + gxz bx bz + gyz by bz) + // = gtt bt^2 + gtxx bx^2 + gtyy by^2 + gtzz bz^2 + 2 (gtx bt bx + gty bt by + gtz bt bz + gtxy bx by + gtxz bx bz + gtyz by bz) // = (-al^2 + gamma_{ij} betai betaj) bt^2 + b^i b^j gamma_{ij} + 2 g_{t i} b^t b^i // = - (alpha b^t)^2 + (b^t)^2 gamma_{ij} beta^i beta^j + b^i b^j gamma_{ij} + 2 b^t g_{t i} b^i // = - (alpha b^t)^2 + (b^t)^2 gamma_{ij} beta^i beta^j + b^i b^j gamma_{ij} + 2 b^t (gamma_{ij} beta^j) b^i

File [modified]: driver_conserv_to_prims.C Delta lines: +16 -16 =================================================================== --- trunk/src/driver_conserv_to_prims.C 2013-10-15 20:31:57 UTC (rev 11) +++ trunk/src/driver_conserv_to_prims.C 2013-10-15 20:35:06 UTC (rev 12) @@ -41,7 +41,7 @@ #define SQR(x) ((x) * (x))

extern "C" void BSSN_compute_gupij(const cGH *cctkGH,int *cctk_lsh, - double *gxx,double *gxy,double *gxz,double *gyy,double *gyz,double *gzz, + double *gtxx,double *gtxy,double *gtxz,double *gtyy,double *gtyz,double *gtzz, double *gupxx,double *gupxy,double *gupxz,double *gupyy,double *gupyz,double *gupzz);

@@ -54,7 +54,7 @@

// gupij is not evolved, and so is not defined anywhere that the grid has moved. // Here we recompute gupij from the newly updated gij's: - BSSN_compute_gupij(cctkGH,cctk_lsh, gxx,gxy,gxz,gyy,gyz,gzz, gupxx,gupxy,gupxz,gupyy,gupyz,gupzz); + BSSN_compute_gupij(cctkGH,cctk_lsh, gtxx,gtxy,gtxz,gtyy,gtyz,gtzz, gupxx,gupxy,gupxz,gupyy,gupyz,gupzz);

if(Symmetry==1) { // SET SYMMETRY GHOSTZONES ON ALL CONSERVATIVE VARIABLES! @@ -169,12 +169,12 @@ double By_f1o4pa = By[index]*f1o4pa; double Bz_f1o4pa = Bz[index]*f1o4pa;

- double gxx_phys = gxx[index]*Psi4; - double gxy_phys = gxy[index]*Psi4; - double gxz_phys = gxz[index]*Psi4; - double gyy_phys = gyy[index]*Psi4; - double gyz_phys = gyz[index]*Psi4; - double gzz_phys = gzz[index]*Psi4; + double gxx_phys = gtxx[index]*Psi4; + double gxy_phys = gtxy[index]*Psi4; + double gxz_phys = gtxz[index]*Psi4; + double gyy_phys = gtyy[index]*Psi4; + double gyz_phys = gtyz[index]*Psi4; + double gzz_phys = gtzz[index]*Psi4;

double B_xl = (gxx_phys *Bx[index] + gxy_phys * By[index] + gxz_phys * Bz[index] ); @@ -214,7 +214,7 @@ apply_tau_floor(index,tau_atm,rho_b_atm, Bx,By,Bz, tau,rho_star,mhd_st_x,mhd_st_y,mhd_st_z, - phi,gxx,gxy,gxz,gyy,gyz,gzz, + phi,gtxx,gtxy,gtxz,gtyy,gtyz,gtzz, gupxx,gupxy,gupxz,gupyy,gupyz,gupzz, new_primitives, lapm1,shiftx,shifty,shiftz, @@ -236,7 +236,7 @@ stats.font_fixed=0; for(int ii=0;ii<3;ii++) { check = harm_primitives_gammalaw_lowlevel(index,i,j,k,x,y,z, - phi,gxx,gxy,gxz,gyy,gyz,gzz, + phi,gtxx,gtxy,gtxz,gtyy,gtyz,gtzz, gupxx,gupxy,gupxz,gupyy,gupyz,gupzz, lapm1,shiftx,shifty,shiftz, Bx,By,Bz, @@ -343,12 +343,12 @@ myfile.write((char*)y, (fullsize)*sizeof(double)); myfile.write((char*)z, (fullsize)*sizeof(double)); myfile.write((char*)phi, (fullsize)*sizeof(double)); - myfile.write((char*)gxx, (fullsize)*sizeof(double)); - myfile.write((char*)gxy, (fullsize)*sizeof(double)); - myfile.write((char*)gxz, (fullsize)*sizeof(double)); - myfile.write((char*)gyy, (fullsize)*sizeof(double)); - myfile.write((char*)gyz, (fullsize)*sizeof(double)); - myfile.write((char*)gzz, (fullsize)*sizeof(double)); + myfile.write((char*)gtxx, (fullsize)*sizeof(double)); + myfile.write((char*)gtxy, (fullsize)*sizeof(double)); + myfile.write((char*)gtxz, (fullsize)*sizeof(double)); + myfile.write((char*)gtyy, (fullsize)*sizeof(double)); + myfile.write((char*)gtyz, (fullsize)*sizeof(double)); + myfile.write((char*)gtzz, (fullsize)*sizeof(double));

myfile.write((char*)gupxx, (fullsize)*sizeof(double)); myfile.write((char*)gupxy, (fullsize)*sizeof(double));

File [modified]: driver_evaluate_MHD_rhs.C Delta lines: +18 -12 =================================================================== --- trunk/src/driver_evaluate_MHD_rhs.C 2013-10-15 20:31:57 UTC (rev 11) +++ trunk/src/driver_evaluate_MHD_rhs.C 2013-10-15 20:35:06 UTC (rev 12) @@ -102,12 +102,12 @@ ww=0; metric[ww]=phi; ww++; metric[ww]=psi; ww++; - metric[ww]=gxx; ww++; - metric[ww]=gxy; ww++; - metric[ww]=gxz; ww++; - metric[ww]=gyy; ww++; - metric[ww]=gyz; ww++; - metric[ww]=gzz; ww++; + metric[ww]=gtxx; ww++; + metric[ww]=gtxy; ww++; + metric[ww]=gtxz; ww++; + metric[ww]=gtyy; ww++; + metric[ww]=gtyz; ww++; + metric[ww]=gtzz; ww++; metric[ww]=lapm1; ww++; metric[ww]=shiftx;ww++; metric[ww]=shifty;ww++; @@ -149,11 +149,11 @@ // Here, we: // 1) Compute tau_rhs extrinsic curvature terms, and // 2) Compute TUPmunu. - // This function is housed in the file: "initialize_hydro_rhss__compute_tau_rhs_extrinsic_curvature_terms_and_TUPmunu.C" - initialize_hydro_rhss__compute_tau_rhs_extrinsic_curvature_terms_and_TUPmunu(cctkGH,cctk_lsh,dX,metric,in_prims,TUPmunu,u0,eos, - gupxy,gupxz,gupyz, - Axx,Axy,Axz,Ayy,Ayz,Azz,trK, - rho_star_rhs,tau_rhs,st_x_rhs,st_y_rhs,st_z_rhs); + // This function is housed in the file: "compute_tau_rhs_extrinsic_curvature_terms_and_TUPmunu.C" + compute_tau_rhs_extrinsic_curvature_terms_and_TUPmunu(cctkGH,cctk_lsh,dX,metric,in_prims,TUPmunu,u0,eos, + gupxy,gupxz,gupyz, + kxx,kxy,kxz,kyy,kyz,kzz, + tau_rhs}; int flux_dirn; flux_dirn=1; // First compute ftilde, which is used for flattening left and right face values @@ -189,6 +189,8 @@ // This function is housed in the file: "reconstruct_set_of_prims_PPM.C" reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct, eos,in_prims,out_prims_r,out_prims_l,ftilde_gf,temporary); + //Right and left face values of BI_CENTER are used in mhdflux computation (first to compute b^a). + // Instead of reconstructing, we simply set B^x face values to be consistent with BX_STAGGER. #pragma omp parallel for for(int k=0;k<cctk_lsh[2];k++) for(int j=0;j<cctk_lsh[1];j++) for(int i=0;i<cctk_lsh[0];i++) { int index=CCTK_GFINDEX3D(cctkGH,i,j,k), indexim1=CCTK_GFINDEX3D(cctkGH,i-1+(i==0),j,k); @@ -263,6 +265,8 @@ // This function is housed in the file: "reconstruct_set_of_prims_PPM.C" reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct, eos,in_prims,out_prims_r,out_prims_l,ftilde_gf,temporary); + //Right and left face values of BI_CENTER are used in mhdflux computation (first to compute b^a). + // Instead of reconstructing, we simply set B^y face values to be consistent with BY_STAGGER. #pragma omp parallel for for(int k=0;k<cctk_lsh[2];k++) for(int j=0;j<cctk_lsh[1];j++) for(int i=0;i<cctk_lsh[0];i++) { int index=CCTK_GFINDEX3D(cctkGH,i,j,k), indexjm1=CCTK_GFINDEX3D(cctkGH,i,j-1+(j==0),k); @@ -358,6 +362,8 @@ // This function is housed in the file: "reconstruct_set_of_prims_PPM.C" reconstruct_set_of_prims_PPM(cctkGH,cctk_lsh,flux_dirn,num_prims_to_reconstruct,which_prims_to_reconstruct, eos,in_prims,out_prims_r,out_prims_l,ftilde_gf,temporary); + //Right and left face values of BI_CENTER are used in mhdflux computation (first to compute b^a). + // Instead of reconstructing, we simply set B^z face values to be consistent with BZ_STAGGER. #pragma omp parallel for for(int k=0;k<cctk_lsh[2];k++) for(int j=0;j<cctk_lsh[1];j++) for(int i=0;i<cctk_lsh[0];i++) { int index=CCTK_GFINDEX3D(cctkGH,i,j,k), indexkm1=CCTK_GFINDEX3D(cctkGH,i,j,k-1+(k==0)); @@ -503,7 +509,7 @@ // convenient to edit a 600 line file than an 1800 line file, so I'd prefer to leave this unconventional structure // alone. #include "reconstruct_set_of_prims_PPM.C" -#include "initialize_hydro_rhss__compute_tau_rhs_extrinsic_curvature_terms_and_TUPmunu.C" +#include "compute_tau_rhs_extrinsic_curvature_terms_and_TUPmunu.C" #include "add_fluxes_and_source_terms_to_hydro_rhss.C" #include "mhdflux.C" #include "A_i_rhs_no_gauge_terms.C"

File [modified]: driver_prims_to_conserv.C Delta lines: +6 -6 =================================================================== --- trunk/src/driver_prims_to_conserv.C 2013-10-15 20:31:57 UTC (rev 11) +++ trunk/src/driver_prims_to_conserv.C 2013-10-15 20:35:06 UTC (rev 12) @@ -42,12 +42,12 @@ double Psim4 = 1.0/Psi4;

double METRIC[NUMVARS_FOR_METRIC]; - METRIC[GXX]=gxx[index]; - METRIC[GXY]=gxy[index]; - METRIC[GXZ]=gxz[index]; - METRIC[GYY]=gyy[index]; - METRIC[GYZ]=gyz[index]; - METRIC[GZZ]=gzz[index]; + METRIC[GXX]=gtxx[index]; + METRIC[GXY]=gtxy[index]; + METRIC[GXZ]=gtxz[index]; + METRIC[GYY]=gtyy[index]; + METRIC[GYZ]=gtyz[index]; + METRIC[GZZ]=gtzz[index]; METRIC[GUPXX]=gupxx[index]; METRIC[GUPXY]=gupxy[index]; METRIC[GUPXZ]=gupxz[index];

File [modified]: evaluate_MHD_rhs_headers.h Delta lines: +3 -3 =================================================================== --- trunk/src/evaluate_MHD_rhs_headers.h 2013-10-15 20:31:57 UTC (rev 11) +++ trunk/src/evaluate_MHD_rhs_headers.h 2013-10-15 20:35:06 UTC (rev 12) @@ -60,11 +60,11 @@ eos_struct &eosi,gf_and_gz_struct *in_prims,gf_and_gz_struct *out_prims_r,gf_and_gz_struct *out_prims_l, double *ftilde_gf,double *temporary);

-void initialize_hydro_rhss__compute_tau_rhs_extrinsic_curvature_terms_and_TUPmunu +void compute_tau_rhs_extrinsic_curvature_terms_and_TUPmunu (const cGH *cctkGH,int *ext,double *dX,double **metric,gf_and_gz_struct *prims,double **TUPmunu,double *u0,eos_struct &eos, double *gupxy,double *gupxz,double *gupyz, - double *Axx,double *Axy,double *Axz,double *Ayy,double *Ayz,double *Azz,double *trK, - double *rho_star_rhs,double *tau_rhs,double *st_x_rhs,double *st_y_rhs,double *st_z_rhs); + double *kxx,double *kxy,double *kxz,double *kyy,double *kyz,double *kzz, + double *tau_rhs); void A_i_rhs_no_gauge_terms(int &A_dirn, const cGH *cctkGH,int *ext,gf_and_gz_struct *out_prims_r,gf_and_gz_struct *out_prims_l, double *phi_interped,double *cmax_1,double *cmin_1,double *cmax_2,double *cmin_2, double *A3_rhs);

File [modified]: harm_primitives_headers.h Delta lines: +2 -2 =================================================================== --- trunk/src/harm_primitives_headers.h 2013-10-15 20:31:57 UTC (rev 11) +++ trunk/src/harm_primitives_headers.h 2013-10-15 20:35:06 UTC (rev 12) @@ -30,7 +30,7 @@ int apply_tau_floor(int &index,const double &tau_atm,const double &rhobatm, double *Bx,double *By,double *Bz, double *tau,double *rho_star,double *mhd_st_x,double *mhd_st_y,double *mhd_st_z, - double *phi,double *gxx,double *gxy,double *gxz,double *gyy,double *gyz,double *gzz, + double *phi,double *gtxx,double *gtxy,double *gtxz,double *gtyy,double *gtyz,double *gtzz, double *gupxx,double *gupxy,double *gupxz,double *gupyy,double *gupyz,double *gupzz, struct output_primitives &new_primitives, double *lapm1,double *shiftx,double *shifty,double *shiftz, @@ -38,7 +38,7 @@ const double &Psi6threshold);

int harm_primitives_gammalaw_lowlevel(int &index,int &i,int &j,int &k,double *X,double *Y,double *Z, - double *phi,double *gxx,double *gxy,double *gxz,double *gyy,double *gyz,double *gzz, + double *phi,double *gtxx,double *gtxy,double *gtxz,double *gtyy,double *gtyz,double *gtzz, double *gupxx,double *gupxy,double *gupxz,double *gupyy,double *gupyz,double *gupzz, double *lapm1,double *shiftx,double *shifty,double *shiftz, double *Bx,double *By,double *Bz,

File [modified]: harm_primitives_lowlevel.C Delta lines: +8 -8 =================================================================== --- trunk/src/harm_primitives_lowlevel.C 2013-10-15 20:31:57 UTC (rev 11) +++ trunk/src/harm_primitives_lowlevel.C 2013-10-15 20:35:06 UTC (rev 12) @@ -2,7 +2,7 @@ #include "harm_primitives_recompute_conservs_tau_floor.C"

int harm_primitives_gammalaw_lowlevel(int &index,int &i,int &j,int &k,double *X,double *Y,double *Z, - double *phi,double *gxx,double *gxy,double *gxz,double *gyy,double *gyz,double *gzz, + double *phi,double *gtxx,double *gtxy,double *gtxz,double *gtyy,double *gtyz,double *gtzz, double *gupxx,double *gupxy,double *gupxz,double *gupyy,double *gupyz,double *gupzz, double *lapm1,double *shiftx,double *shifty,double *shiftz, double *Bx,double *By,double *Bz, @@ -43,12 +43,12 @@ double shiftzL = shiftz[index];

//Define physical metric & its inverse: - double gxx_physL = gxx[index]*Psi4; - double gxy_physL = gxy[index]*Psi4; - double gxz_physL = gxz[index]*Psi4; - double gyy_physL = gyy[index]*Psi4; - double gyz_physL = gyz[index]*Psi4; - double gzz_physL = gzz[index]*Psi4; + double gxx_physL = gtxx[index]*Psi4; + double gxy_physL = gtxy[index]*Psi4; + double gxz_physL = gtxz[index]*Psi4; + double gyy_physL = gtyy[index]*Psi4; + double gyz_physL = gtyz[index]*Psi4; + double gzz_physL = gtzz[index]*Psi4; double gupxx_physL = gupxx[index]*Psim4; double gupxy_physL = gupxy[index]*Psim4; double gupxz_physL = gupxz[index]*Psim4; @@ -58,7 +58,7 @@

// Check to see if the metric is positive-defitive double lam1,lam2,lam3; - double M11 = gxx[index], M12=gxy[index], M13=gxz[index], M22=gyy[index], M23=gyz[index], M33=gzz[index]; + double M11 = gtxx[index], M12=gtxy[index], M13=gtxz[index], M22=gtyy[index], M23=gtyz[index], M33=gtzz[index]; eigenvalues_3by3_real_sym_matrix(lam1, lam2, lam3,M11, M12, M13, M22, M23, M33); if (lam1 < 0.0 || lam2 < 0.0 || lam3 < 0) { // Metric is not positive-defitive, reset the physical metric to be conformally-flat.

File [modified]: harm_primitives_recompute_conservs_tau_floor.C Delta lines: +30 -30 =================================================================== --- trunk/src/harm_primitives_recompute_conservs_tau_floor.C 2013-10-15 20:31:57 UTC (rev 11) +++ trunk/src/harm_primitives_recompute_conservs_tau_floor.C 2013-10-15 20:35:06 UTC (rev 12) @@ -6,7 +6,7 @@ int apply_tau_floor(int &index,const double &tau_atm,const double &rhobatm, double *Bx,double *By,double *Bz, double *tau,double *rho_star,double *mhd_st_x,double *mhd_st_y,double *mhd_st_z, - double *phi,double *gxx,double *gxy,double *gxz,double *gyy,double *gyz,double *gzz, + double *phi,double *gtxx,double *gtxy,double *gtxz,double *gtyy,double *gtyz,double *gtzz, double *gupxx,double *gupxy,double *gupxz,double *gupyy,double *gupyz,double *gupzz, struct output_primitives &new_primitives, double *lapm1,double *shiftx,double *shifty,double *shiftz, @@ -28,32 +28,32 @@ //rho_star = alpha u0 Psi6 rho_b, alpha u0 > 1, so if rho_star < Psi6 rhobatm, then we are GUARANTEED that we can reset to atmosphere. //if(rho_star[index] < 1e4*Psi6*rhobatm) { //if(rho_star[index] < 2*Psi6*rhobatm) { - double gxxL = gxx[index]; - double gxyL = gxy[index]; - double gxzL = gxz[index]; - double gyyL = gyy[index]; - double gyzL = gyz[index]; - double gzzL = gzz[index]; + double gtxxL = gtxx[index]; + double gtxyL = gtxy[index]; + double gtxzL = gtxz[index]; + double gtyyL = gtyy[index]; + double gtyzL = gtyz[index]; + double gtzzL = gtzz[index];

double lam1,lam2,lam3; - eigenvalues_3by3_real_sym_matrix(lam1, lam2, lam3,gxxL, gxyL, gxzL, gyyL, gyzL, gzzL); + eigenvalues_3by3_real_sym_matrix(lam1, lam2, lam3,gtxxL, gtxyL, gtxzL, gtyyL, gtyzL, gtzzL);

- double gxxi=gupxx[index]*Psim4,gyyi=gupyy[index]*Psim4,gzzi=gupzz[index]*Psim4,gxyi=gupxy[index]*Psim4,gxzi=gupxz[index]*Psim4,gyzi=gupyz[index]*Psim4; + double gtxxi=gupxx[index]*Psim4,gtyyi=gupyy[index]*Psim4,gtzzi=gupzz[index]*Psim4,gtxyi=gupxy[index]*Psim4,gtxzi=gupxz[index]*Psim4,gtyzi=gupyz[index]*Psim4;

if (lam1 < 0.0 || lam2 < 0.0 || lam3 < 0.0) { // Metric is not positive-defitive, reset the metric to be conformally-flat. - gxxL = 1.0; - gxyL = 0.0; - gxzL = 0.0; - gyyL = 1.0; - gyzL = 0.0; - gzzL = 1.0; - gxxi = Psim4; - gxyi = 0.0; - gxzi = 0.0; - gyyi = Psim4; - gyzi = 0.0; - gzzi = Psim4; + gtxxL = 1.0; + gtxyL = 0.0; + gtxzL = 0.0; + gtyyL = 1.0; + gtyzL = 0.0; + gtzzL = 1.0; + gtxxi = Psim4; + gtxyi = 0.0; + gtxzi = 0.0; + gtyyi = Psim4; + gtyzi = 0.0; + gtzzi = Psim4; }

/* DO NOT UNCOMMENT THESE LINES; THIS IS A TERRIBLE PRESCRIPTION. @@ -77,12 +77,12 @@ double ByL = By[index]; double BzL = Bz[index];

- double B2 = Psi4*(gxxL*SQR(BxL) + - 2.0*gxyL*BxL*ByL + - 2.0*gxzL*BxL*BzL + - gyyL*SQR(ByL) + - 2.0*gyzL*ByL*BzL + - gzzL*SQR(BzL) ) ; + double B2 = Psi4*(gtxxL*SQR(BxL) + + 2.0*gtxyL*BxL*ByL + + 2.0*gtxzL*BxL*BzL + + gtyyL*SQR(ByL) + + 2.0*gtyzL*ByL*BzL + + gtzzL*SQR(BzL) ) ;

new_primitives.u0_new = 1.0/(lapm1[index]+1.0); new_primitives.vx_new = -shiftx[index]; @@ -109,7 +109,7 @@ double f1os4p = 1.0/sqrt(4.0*M_PI); double Bxbar = Bx[index]*f1os4p,Bybar = By[index]*f1os4p,Bzbar = Bz[index]*f1os4p;

- double gxx_physL=gxxL*Psi4,gyy_physL=gyyL*Psi4,gzz_physL=gzzL*Psi4,gxy_physL=gxyL*Psi4,gxz_physL=gxzL*Psi4,gyz_physL=gyzL*Psi4; + double gxx_physL=gtxxL*Psi4,gyy_physL=gtyyL*Psi4,gzz_physL=gtzzL*Psi4,gxy_physL=gtxyL*Psi4,gxz_physL=gtxzL*Psi4,gyz_physL=gtyzL*Psi4; double Bbar_x = gxx_physL*Bxbar + gxy_physL*Bybar + gxz_physL*Bzbar; double Bbar_y = gxy_physL*Bxbar + gyy_physL*Bybar + gyz_physL*Bzbar; double Bbar_z = gxz_physL*Bxbar + gyz_physL*Bybar + gzz_physL*Bzbar; @@ -153,8 +153,8 @@ // mhd_st_x[index] = stxi; mhd_st_y[index] = styi; mhd_st_z[index] = stzi; //}

- double sdots= gxxi*SQR(stxi)+gyyi*SQR(styi)+gzzi*SQR(stzi)+2.0* - (gxyi*stxi*styi+gxzi*stxi*stzi+gyzi*styi*stzi); + double sdots= gtxxi*SQR(stxi)+gtyyi*SQR(styi)+gtzzi*SQR(stzi)+2.0* + (gtxyi*stxi*styi+gtxzi*stxi*stzi+gtyzi*styi*stzi);

double Wm = sqrt(SQR(hatBbardotS)+ SQR(rho_s))/Psi6; double Sm2 = (SQR(Wm)*sdots + SQR(BbardotS)*(Bbar2+2.0*Wm))/SQR(Wm+Bbar2);