[Commits] [svn:einsteintoolkit] IllinoisGRMHD/trunk/src/ (Rev. 13)
zachetie at gmail.com
zachetie at gmail.com
Tue Oct 15 20:41:19 CDT 2013
User: zetienne
Date: 2013/10/15 08:41 PM
Added:
/trunk/src/
convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.C
Modified:
/trunk/src/
driver_conserv_to_prims.C, driver_evaluate_MHD_rhs.C, driver_prims_to_conserv.C, evaluate_MHD_rhs_headers.h, harm_primitives_recompute_conservs_tau_floor.C, make.code.defn, mhdflux.C
Log:
Install (hopefully) all necessary changes to allow for ADMBase variables as input.
File Changes:
Directory: /trunk/src/
======================
File [added]: convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.C
Delta lines: +73 -0
===================================================================
--- trunk/src/convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.C (rev 0)
+++ trunk/src/convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.C 2013-10-16 01:41:18 UTC (rev 13)
@@ -0,0 +1,73 @@
+void convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij(double *gxx,double *gxy,double *gxz,double *gyy,double *gyz,double *gzz,double *alp,
+ double *gtxx,double *gtxy,double *gtxz,double *gtyy,double *gtyz,double *gtzz,
+ double *gtupxx,double *gtupxy,double *gtupxz,double *gtupyy,double *gtupyz,double *gtupzz,
+ double *phi,double *psi,double *lapm1) {
+#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);
+ double gxx_physL=gxx[index];
+ double gxy_physL=gxy[index];
+ double gxz_physL=gxz[index];
+ double gyy_physL=gyy[index];
+ double gyz_physL=gyz[index];
+ double gzz_physL=gzz[index];
+
+ /**********************************************************************
+ * Compute \tilde{\gamma_{ij}}, phi, and psi (BSSN) from g_{ij} (ADM) *
+ **********************************************************************/
+ double gijdet = gxx_physL * gyy_physL * gzz_physL + gxy_physL * gyz_physL * gxz_physL + gxz_physL * gxy_physL * gyz_physL -
+ gxz_physL * gyy_physL * gxz_physL - gxy_physL * gxy_physL * gzz_physL - gxx_physL * gyz_physL * gyz_physL;
+
+ static const double F1o12 = 0.08333333333333333333;
+ double phiL = F1o12 * log(gijdet);
+ double psiL = exp(phiL);
+
+ double Psim4 = 1.0/(psiL*psiL*psiL*psiL);
+ double gtxxL = gxx_physL*Psim4;
+ double gtxyL = gxy_physL*Psim4;
+ double gtxzL = gxz_physL*Psim4;
+ double gtyyL = gyy_physL*Psim4;
+ double gtyzL = gyz_physL*Psim4;
+ double gtzzL = gzz_physL*Psim4;
+
+ /*********************************
+ * Apply det gtij = 1 constraint *
+ *********************************/
+ double gtijdet = gtxxL * gtyyL * gtzzL + gtxyL * gtyzL * gtxzL + gtxzL * gtxyL * gtyzL -
+ gtxzL * gtyyL * gtxzL - gtxyL * gtxyL * gtzzL - gtxxL * gtyzL * gtyzL;
+
+ static const double Fm1o3 = -0.3333333333333333333333333333333;
+ double gtijdet_Fm1o3 = pow(fabs(gtijdet),Fm1o3);
+
+ gtxxL = gtxxL * gtijdet_Fm1o3;
+ gtxyL = gtxyL * gtijdet_Fm1o3;
+ gtxzL = gtxzL * gtijdet_Fm1o3;
+ gtyyL = gtyyL * gtijdet_Fm1o3;
+ gtyzL = gtyzL * gtijdet_Fm1o3;
+ gtzzL = gtzzL * gtijdet_Fm1o3;
+
+ /*****************************************
+ * Set all the needed BSSN gridfunctions *
+ *****************************************/
+ phi[index] = phiL;
+ psi[index] = psiL;
+
+ lapm1[index] = alp[index] - 1.0;
+
+ gtxx[index] = gtxxL;
+ gtxy[index] = gtxyL;
+ gtxz[index] = gtxzL;
+ gtyy[index] = gtyyL;
+ gtyz[index] = gtyzL;
+ gtzz[index] = gtzzL;
+
+ gtupxx[index] = ( gtyyL * gtzzL - gtyzL * gtyzL );
+ gtupxy[index] = - ( gtxyL * gtzzL - gtyzL * gtxzL );
+ gtupxz[index] = ( gtxyL * gtyzL - gtyyL * gtxzL );
+ gtupyy[index] = ( gtxxL * gtzzL - gtxzL * gtxzL );
+ gtupyz[index] = - ( gtxxL * gtyzL - gtxyL * gtxzL );
+ gtupzz[index] = ( gtxxL * gtyyL - gtxyL * gtxyL );
+ }
+}
+
+
File [modified]: driver_conserv_to_prims.C
Delta lines: +19 -21
===================================================================
--- trunk/src/driver_conserv_to_prims.C 2013-10-15 20:35:06 UTC (rev 12)
+++ trunk/src/driver_conserv_to_prims.C 2013-10-16 01:41:18 UTC (rev 13)
@@ -40,11 +40,6 @@
#define SQR(x) ((x) * (x))
-extern "C" void BSSN_compute_gupij(const cGH *cctkGH,int *cctk_lsh,
- double *gtxx,double *gtxy,double *gtxz,double *gtyy,double *gtyz,double *gtzz,
- double *gupxx,double *gupxy,double *gupxz,double *gupyy,double *gupyz,double *gupzz);
-
-
extern "C" void conserv_to_prims(CCTK_ARGUMENTS) {
DECLARE_CCTK_ARGUMENTS;
DECLARE_CCTK_PARAMETERS;
@@ -52,9 +47,12 @@
// We use proper C++ here, for file I/O later.
using namespace std;
- // 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, gtxx,gtxy,gtxz,gtyy,gtyz,gtzz, gupxx,gupxy,gupxz,gupyy,gupyz,gupzz);
+ // These BSSN-based variables are not evolved, and so are not defined anywhere that the grid has moved.
+ // Here we convert ADM variables (from ADMBase) to the BSSN-based variables expected by this routine.
+ convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij(gxx,gxy,gxz,gyy,gyz,gzz,alp,
+ gtxx,gtxy,gtxz,gtyy,gtyz,gtzz,
+ gtupxx,gtupxy,gtupxz,gtupyy,gtupyz,gtupzz,
+ phi,psi,lapm1);
if(Symmetry==1) {
// SET SYMMETRY GHOSTZONES ON ALL CONSERVATIVE VARIABLES!
@@ -215,9 +213,9 @@
Bx,By,Bz,
tau,rho_star,mhd_st_x,mhd_st_y,mhd_st_z,
phi,gtxx,gtxy,gtxz,gtyy,gtyz,gtzz,
- gupxx,gupxy,gupxz,gupyy,gupyz,gupzz,
+ gtupxx,gtupxy,gtupxz,gtupyy,gtupyz,gtupzz,
new_primitives,
- lapm1,shiftx,shifty,shiftz,
+ lapm1,betax,betay,betaz,
neos,rho_tab,P_tab,eps_tab,k_tab,gamma_tab,
Psi6threshold);
@@ -237,8 +235,8 @@
for(int ii=0;ii<3;ii++) {
check = harm_primitives_gammalaw_lowlevel(index,i,j,k,x,y,z,
phi,gtxx,gtxy,gtxz,gtyy,gtyz,gtzz,
- gupxx,gupxy,gupxz,gupyy,gupyz,gupzz,
- lapm1,shiftx,shifty,shiftz,
+ gtupxx,gtupxy,gtupxz,gtupyy,gtupyz,gtupzz,
+ lapm1,betax,betay,betaz,
Bx,By,Bz,
mhd_st_x,mhd_st_y,mhd_st_z,tau,rho_star,
vx,vy,vz,P,rho_b,u0,
@@ -350,16 +348,16 @@
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));
- myfile.write((char*)gupxz, (fullsize)*sizeof(double));
- myfile.write((char*)gupyy, (fullsize)*sizeof(double));
- myfile.write((char*)gupyz, (fullsize)*sizeof(double));
- myfile.write((char*)gupzz, (fullsize)*sizeof(double));
+ myfile.write((char*)gtupxx, (fullsize)*sizeof(double));
+ myfile.write((char*)gtupxy, (fullsize)*sizeof(double));
+ myfile.write((char*)gtupxz, (fullsize)*sizeof(double));
+ myfile.write((char*)gtupyy, (fullsize)*sizeof(double));
+ myfile.write((char*)gtupyz, (fullsize)*sizeof(double));
+ myfile.write((char*)gtupzz, (fullsize)*sizeof(double));
- myfile.write((char*)shiftx, (fullsize)*sizeof(double));
- myfile.write((char*)shifty, (fullsize)*sizeof(double));
- myfile.write((char*)shiftz, (fullsize)*sizeof(double));
+ myfile.write((char*)betax, (fullsize)*sizeof(double));
+ myfile.write((char*)betay, (fullsize)*sizeof(double));
+ myfile.write((char*)betaz, (fullsize)*sizeof(double));
myfile.write((char*)lapm1, (fullsize)*sizeof(double));
File [modified]: driver_evaluate_MHD_rhs.C
Delta lines: +28 -28
===================================================================
--- trunk/src/driver_evaluate_MHD_rhs.C 2013-10-15 20:35:06 UTC (rev 12)
+++ trunk/src/driver_evaluate_MHD_rhs.C 2013-10-16 01:41:18 UTC (rev 13)
@@ -92,29 +92,29 @@
for(int i=0;i<MAXNUMVARS;i++) for(int j=1;j<=3;j++) { out_prims_r[i].gz_lo[j]=0; out_prims_r[i].gz_hi[j]=0; }
for(int i=0;i<MAXNUMVARS;i++) for(int j=1;j<=3;j++) { out_prims_l[i].gz_lo[j]=0; out_prims_l[i].gz_hi[j]=0; }
-#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);
- psi[index] = exp(phi[index]);
- }
+ // Convert ADM variables (from ADMBase) to the BSSN-based variables expected by this routine.
+ convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij(gxx,gxy,gxz,gyy,gyz,gzz,alp,
+ gtxx,gtxy,gtxz,gtyy,gtyz,gtzz,
+ gtupxx,gtupxy,gtupxz,gtupyy,gtupyz,gtupzz,
+ phi,psi,lapm1);
double *metric[NUMVARS_FOR_METRIC_FACEVALS]; // "metric" here is array of pointers to the actual gridfunctions.
ww=0;
metric[ww]=phi; ww++;
metric[ww]=psi; 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]=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++;
- metric[ww]=shiftz;ww++;
- metric[ww]=gupxx; ww++;
- metric[ww]=gupyy; ww++;
- metric[ww]=gupzz; ww++;
+ metric[ww]=betax; ww++;
+ metric[ww]=betay; ww++;
+ metric[ww]=betaz; ww++;
+ metric[ww]=gtupxx; ww++;
+ metric[ww]=gtupyy; ww++;
+ metric[ww]=gtupzz; ww++;
double *TUPmunu[10];// "TUPmunu" here is array of pointers to the actual gridfunctions.
ww=0;
@@ -151,7 +151,7 @@
// 2) Compute TUPmunu.
// 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,
+ gtupxy,gtupxz,gtupyz,
kxx,kxy,kxz,kyy,kyz,kzz,
tau_rhs};
int flux_dirn;
@@ -461,18 +461,18 @@
// Next compute psi6phi_rhs, and add gauge terms to A_i_rhs terms!
// Note that in the following function, we don't bother with reconstruction, instead interpolating.
- // We need A^i, but only have A_i. So we add gupij to the list of input variables.
+ // We need A^i, but only have A_i. So we add gtupij to the list of input variables.
double *interp_vars[MAXNUMINTERP];
int w=0;
- interp_vars[w]=shiftx; w++;
- interp_vars[w]=shifty; w++;
- interp_vars[w]=shiftz; w++;
- interp_vars[w]=gupxx; w++;
- interp_vars[w]=gupxy; w++;
- interp_vars[w]=gupxz; w++;
- interp_vars[w]=gupyy; w++;
- interp_vars[w]=gupyz; w++;
- interp_vars[w]=gupzz; w++;
+ interp_vars[w]=betax; w++;
+ interp_vars[w]=betay; w++;
+ interp_vars[w]=betaz; w++;
+ interp_vars[w]=gtupxx; w++;
+ interp_vars[w]=gtupxy; w++;
+ interp_vars[w]=gtupxz; w++;
+ interp_vars[w]=gtupyy; w++;
+ interp_vars[w]=gtupyz; w++;
+ interp_vars[w]=gtupzz; w++;
interp_vars[w]=psi; w++;
interp_vars[w]=lapm1; w++;
interp_vars[w]=Ax; w++;
File [modified]: driver_prims_to_conserv.C
Delta lines: +12 -6
===================================================================
--- trunk/src/driver_prims_to_conserv.C 2013-10-15 20:35:06 UTC (rev 12)
+++ trunk/src/driver_prims_to_conserv.C 2013-10-16 01:41:18 UTC (rev 13)
@@ -28,6 +28,12 @@
exit(1);
}
+ // Convert ADM variables (from ADMBase) to the BSSN-based variables expected by this routine.
+ convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij(gxx,gxy,gxz,gyy,gyz,gzz,alp,
+ gtxx,gtxy,gtxz,gtyy,gtyz,gtzz,
+ gtupxx,gtupxy,gtupxz,gtupyy,gtupyz,gtupzz,
+ phi,psi,lapm1);
+
// for(int index=0;index<cctk_lsh[0]*cctk_lsh[1]*cctk_lsh[2];index++) {
#pragma omp parallel for
for(int k=0;k<cctk_lsh[2];k++)
@@ -48,12 +54,12 @@
METRIC[GYY]=gtyy[index];
METRIC[GYZ]=gtyz[index];
METRIC[GZZ]=gtzz[index];
- METRIC[GUPXX]=gupxx[index];
- METRIC[GUPXY]=gupxy[index];
- METRIC[GUPXZ]=gupxz[index];
- METRIC[GUPYY]=gupyy[index];
- METRIC[GUPYZ]=gupyz[index];
- METRIC[GUPZZ]=gupzz[index];
+ METRIC[GUPXX]=gtupxx[index];
+ METRIC[GUPXY]=gtupxy[index];
+ METRIC[GUPXZ]=gtupxz[index];
+ METRIC[GUPYY]=gtupyy[index];
+ METRIC[GUPYZ]=gtupyz[index];
+ METRIC[GUPZZ]=gtupzz[index];
METRIC[SHIFTX]=shiftx[index];
METRIC[SHIFTY]=shifty[index];
File [modified]: evaluate_MHD_rhs_headers.h
Delta lines: +4 -0
===================================================================
--- trunk/src/evaluate_MHD_rhs_headers.h 2013-10-15 20:35:06 UTC (rev 12)
+++ trunk/src/evaluate_MHD_rhs_headers.h 2013-10-16 01:41:18 UTC (rev 13)
@@ -84,3 +84,7 @@
void impose_speed_limit_output_u0(double *METRIC,double *U,double &psi4,double &ONE_OVER_LAPSE,double &u0_out);
+void convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij(double *gxx,double *gxy,double *gxz,double *gyy,double *gyz,double *gzz,double *alp,
+ double *gtxx,double *gtxy,double *gtxz,double *gtyy,double *gtyz,double *gtzz,
+ double *gtupxx,double *gtupxy,double *gtupxz,double *gtupyy,double *gtupyz,double *gtupzz,
+ double *phi,double *psi,double *lapm1);
File [modified]: harm_primitives_recompute_conservs_tau_floor.C
Delta lines: +6 -6
===================================================================
--- trunk/src/harm_primitives_recompute_conservs_tau_floor.C 2013-10-15 20:35:06 UTC (rev 12)
+++ trunk/src/harm_primitives_recompute_conservs_tau_floor.C 2013-10-16 01:41:18 UTC (rev 13)
@@ -7,9 +7,9 @@
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 *gtxx,double *gtxy,double *gtxz,double *gtyy,double *gtyz,double *gtzz,
- double *gupxx,double *gupxy,double *gupxz,double *gupyy,double *gupyz,double *gupzz,
+ double *gtupxx,double *gtupxy,double *gtupxz,double *gtupyy,double *gtupyz,double *gtupzz,
struct output_primitives &new_primitives,
- double *lapm1,double *shiftx,double *shifty,double *shiftz,
+ double *lapm1,double *betax,double *betay,double *betaz,
const int &neos,double *rho_tab,double *P_tab,double *eps_tab,double *k_tab,double *gamma_tab,
const double &Psi6threshold) {
@@ -38,7 +38,7 @@
double lam1,lam2,lam3;
eigenvalues_3by3_real_sym_matrix(lam1, lam2, lam3,gtxxL, gtxyL, gtxzL, gtyyL, gtyzL, gtzzL);
- 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;
+ double gtxxi=gtupxx[index]*Psim4,gtyyi=gtupyy[index]*Psim4,gtzzi=gtupzz[index]*Psim4,gtxyi=gtupxy[index]*Psim4,gtxzi=gtupxz[index]*Psim4,gtyzi=gtupyz[index]*Psim4;
if (lam1 < 0.0 || lam2 < 0.0 || lam3 < 0.0) {
// Metric is not positive-defitive, reset the metric to be conformally-flat.
@@ -85,9 +85,9 @@
gtzzL*SQR(BzL) ) ;
new_primitives.u0_new = 1.0/(lapm1[index]+1.0);
- new_primitives.vx_new = -shiftx[index];
- new_primitives.vy_new = -shifty[index];
- new_primitives.vz_new = -shiftz[index];
+ new_primitives.vx_new = -betax[index];
+ new_primitives.vy_new = -betay[index];
+ new_primitives.vz_new = -betaz[index];
double f1o4p = 1.0/(4.0*M_PI);
double sb2 = B2*f1o4p;
File [modified]: make.code.defn
Delta lines: +2 -0
===================================================================
--- trunk/src/make.code.defn 2013-10-15 20:35:06 UTC (rev 12)
+++ trunk/src/make.code.defn 2013-10-16 01:41:18 UTC (rev 13)
@@ -6,6 +6,8 @@
postinitialdata.C \
driver_prims_to_conserv.C \
\
+ convert_ADM_to_BSSN__enforce_detgtij_eq_1__and_compute_gtupij.C \
+ \
driver_evaluate_MHD_rhs.C \
compute_B_and_Bstagger_from_A.C \
driver_conserv_to_prims.C \
File [modified]: mhdflux.C
Delta lines: +1 -1
===================================================================
--- trunk/src/mhdflux.C 2013-10-15 20:35:06 UTC (rev 12)
+++ trunk/src/mhdflux.C 2013-10-16 01:41:18 UTC (rev 13)
@@ -233,7 +233,7 @@
// X^2 ( (1-v02) u0^2 - v02 g^{00}) + X (2 ux u0 (1-v02) - 2 v02 g^{x0}) + (1-v02) ux^2 - v02 g^{xx}
// a = (1-v02) u0^2 - v02 g^{00} = (1-v02) u0^2 + v02/lapse^2 <-- VERIFIED
// b = 2 ux u0 (1-v02) - 2 v02 shiftx/lapse^2 <-- VERIFIED, X->-X, because X = -w/k_1, and we are solving for -X.
- // c = (1-v02) ux^2 - v02 (gupxx*psim4 - (shiftx/lapse)^2) <-- VERIFIED
+ // c = (1-v02) ux^2 - v02 (gtupxx*psim4 - (betax/lapse)^2) <-- VERIFIED
// v02 = v_A^2 + c_s^2 (1 - v_A^2)
double u0_SQUARED=SQR(u0);
More information about the Commits
mailing list