[Commits] [svn:einsteintoolkit] IllinoisGRMHD/trunk/src/ (Rev. 13)

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