User: rhaas Date: 2013/03/27 08:48 PM

Modified: /trunk/ param.ccl /trunk/src/ GRHydro_BondiM.c

Log: GRHydro_InitData: set Bvec according to UIUC formula

From: Roland Haas rhaas@tapir.caltech.edu

File Changes:

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

File [modified]: GRHydro_BondiM.c Delta lines: +46 -26 =================================================================== --- trunk/src/GRHydro_BondiM.c 2013-03-28 01:48:53 UTC (rev 197) +++ trunk/src/GRHydro_BondiM.c 2013-03-28 01:48:55 UTC (rev 198) @@ -901,6 +901,8 @@ -- driver routine for the Bondi solution;

***********************************************************************************/ +static void GRHydro_BondiM_Internal(CCTK_ARGUMENTS, CCTK_REAL range_min, CCTK_REAL range_max); + void GRHydro_BondiM(CCTK_ARGUMENTS) { GRHydro_BondiM_Internal(CCTK_PASS_CTOC, 1e100, -1e100); @@ -914,17 +916,18 @@ }

-void GRHydro_BondiM_Internal(CCTK_ARGUMENTS, CCTK_REAL range_min, CCTK_REAL range_max) +static void GRHydro_BondiM_Internal(CCTK_ARGUMENTS, CCTK_REAL range_min, CCTK_REAL range_max) { DECLARE_CCTK_ARGUMENTS; DECLARE_CCTK_PARAMETERS; - CCTK_REAL det; + CCTK_REAL det, inv_alp; CCTK_REAL rhotmp, utmp, vtmp, rspher, xpos[NDIM], x_spher[NDIM], ucon[NDIM], bcon[NDIM]; int retval; CCTK_REAL *r_bondi, *logr_bondi, *rho_bondi, *u_bondi, *v_bondi; CCTK_REAL dlogr,logrmin; CCTK_REAL rmin_bondi,rmax_bondi; -// CCTK_INT GRHydro_reflevel=0; + int bvec_method; + enum {BVEC_DIRECT, BVEC_TRANSFORM};

int nbondi;

@@ -946,14 +949,26 @@

if( CCTK_EQUALS(bondi_coordinates,"Boyer-Lindquist") ) { coord_type = COORD_BOYERLINDQUIST ; - } - if( CCTK_EQUALS(bondi_coordinates,"Kerr-Schild") ) { + } else if( CCTK_EQUALS(bondi_coordinates,"Kerr-Schild") ) { coord_type = COORD_KERRSCHILD ; - } - if( CCTK_EQUALS(bondi_coordinates,"Isotropic") ) { + } else if( CCTK_EQUALS(bondi_coordinates,"Isotropic") ) { coord_type = COORD_ISOTROPIC ; + } else { + CCTK_VWarn(CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING, + "Unknown coordinate type '%s'", bondi_coordinates); }

+ if( CCTK_EQUALS(bondi_Bvec_method, "direct")) { + bvec_method = BVEC_DIRECT; + } else if( CCTK_EQUALS(bondi_Bvec_method, "transform")) { + bvec_method = BVEC_TRANSFORM; + } else { + CCTK_VWarn(CCTK_WARN_ABORT, __LINE__, __FILE__, CCTK_THORNSTRING, + "Unknown bvec setup method '%s'", bondi_Bvec_method); + } + + + /* xyz location of the black hole : */ i = 0; x_cen = bh_bondi_pos_x[i] ; @@ -1215,11 +1230,11 @@ } }

- det = 1./alp[i]; /* temp var */ + inv_alp = 1./alp[i];

- velx(i) = (ucon[XX]/ucon[TT] + betax[i]) * det; - vely(i) = (ucon[YY]/ucon[TT] + betay[i]) * det; - velz(i) = (ucon[ZZ]/ucon[TT] + betaz[i]) * det; + velx(i) = (ucon[XX]/ucon[TT] + betax[i]) * inv_alp; + vely(i) = (ucon[YY]/ucon[TT] + betay[i]) * inv_alp; + velz(i) = (ucon[ZZ]/ucon[TT] + betaz[i]) * inv_alp;

w_lorentz[i] = 1./sqrt(1-(gxx[i]*SQR(velx(i))+gyy[i]*SQR(vely(i))+gzz[i]*SQR(velz(i))+ 2*(gxy[i]*velx(i)*vely(i)+gxz[i]*velx(i)*velz(i)+gyz[i]*vely(i)*velz(i))) ); @@ -1230,24 +1245,29 @@ 2*(gxy[i]*velx(i)*vely(i)+gxz[i]*velx(i)*velz(i)+gyz[i]*vely(i)*velz(i)))); }

- - Bvecx(i) = w_lorentz[i]*bcon[XX] - alp[i]*bcon[TT]*ucon[XX]; - Bvecy(i) = w_lorentz[i]*bcon[YY] - alp[i]*bcon[TT]*ucon[YY]; - Bvecz(i) = w_lorentz[i]*bcon[ZZ] - alp[i]*bcon[TT]*ucon[ZZ]; - SpatialDet(gxx[i],gxy[i],gxz[i],gyy[i],gyz[i],gzz[i],&det);

- Prim2ConGenM(*GRHydro_eos_handle,gxx[i],gxy[i], - gxz[i],gyy[i],gyz[i],gzz[i], - det, &dens[i],&sx(i),&sy(i),&sz(i), - &tau[i], - &Bconsx(i),&Bconsy(i),&Bconsz(i), - rho[i], - velx(i),vely(i),velz(i), - eps[i],&press[i], - Bvecx(i),Bvecy(i),Bvecz(i), - &w_lorentz[i]); + if(bvec_method == BVEC_TRANSFORM) { + Bvecx(i) = w_lorentz[i]*bcon[XX] - alp[i]*bcon[TT]*ucon[XX]; + Bvecy(i) = w_lorentz[i]*bcon[YY] - alp[i]*bcon[TT]*ucon[YY]; + Bvecz(i) = w_lorentz[i]*bcon[ZZ] - alp[i]*bcon[TT]*ucon[ZZ]; + } else { + Bvecx(i) = bondi_bmag*SQR(M)*x[i]/sqrt(det)/CUBE(r[i]); + Bvecy(i) = bondi_bmag*SQR(M)*y[i]/sqrt(det)/CUBE(r[i]); + Bvecz(i) = bondi_bmag*SQR(M)*z[i]/sqrt(det)/CUBE(r[i]); + }

+ Prim2ConGenM(*GRHydro_eos_handle,gxx[i],gxy[i], + gxz[i],gyy[i],gyz[i],gzz[i], + det, &dens[i],&sx(i),&sy(i),&sz(i), + &tau[i], + &Bconsx(i),&Bconsy(i),&Bconsz(i), + rho[i], + velx(i),vely(i),velz(i), + eps[i],&press[i], + Bvecx(i),Bvecy(i),Bvecz(i), + &w_lorentz[i]); + }

free(r_bondi); free(logr_bondi); free(rho_bondi); free(u_bondi); free(v_bondi);

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

File [modified]: param.ccl Delta lines: +6 -0 =================================================================== --- trunk/param.ccl 2013-03-28 01:48:53 UTC (rev 197) +++ trunk/param.ccl 2013-03-28 01:48:55 UTC (rev 198) @@ -395,6 +395,12 @@ (0:* :: "positive" } 1.0

+CCTK_KEYWORD bondi_Bvec_method "how to compute the magnetic field vector" +{ + "direct" :: "directly from Cartesian metric" + "transform" :: "transform Schwarzschild solution to Kerr Schild" +} "direct" + CCTK_BOOLEAN bondi_evolve_only_annulus "reset to initial data outside of bondi_freeze_inner_radius and bondi_freeze_outer_radius" { } "no"