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);
