User: rhaas Date: 2010/05/28 07:07 AM

Modified: /trunk/ param.ccl /trunk/src/ external.inc

Log: TOVSolver: fix computation of Lorentz factor

Use the Valencia three velocity to compute the Lorentz factor. Add references to the definitions to code and param.ccl. Makes no difference for initial data that is at rest, so the testsuites are unchanged.

File Changes:

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

File [modified]: external.inc Delta lines: +22 -15 =================================================================== --- trunk/src/external.inc 2010-05-03 03:53:11 UTC (rev 110) +++ trunk/src/external.inc 2010-05-28 12:07:30 UTC (rev 111) @@ -91,9 +91,10 @@ rho = pow(press/TOV_K[star], 1.0/TOV_Gamma[star]); eps = (rho > 0.0) ? press/(TOV_Gamma[star] - 1.0) / rho : 0.0; gxx = r / (r_to_star + 1.0e-30) * r / (r_to_star + 1.0e-30); - w_lorentz_2 = 1. - gxx*TOV_Velocity_x[star]*TOV_Velocity_x[star] - - gxx*TOV_Velocity_y[star]*TOV_Velocity_y[star] - - gxx*TOV_Velocity_z[star]*TOV_Velocity_z[star]; + /* velocity components as in gr-qc/9811015 */ + w_lorentz_2 = 1. / (1. - gxx*TOV_Velocity_x[star]*TOV_Velocity_x[star] + - gxx*TOV_Velocity_y[star]*TOV_Velocity_y[star] + - gxx*TOV_Velocity_z[star]*TOV_Velocity_z[star]); if (rho > 0.0) /* source[i]=gxx*gxx*rho*(1.0+eps); */ source[i]=gxx*gxx* (w_lorentz_2*rho + @@ -160,9 +161,10 @@ rho = pow(press/TOV_K[star], 1.0/TOV_Gamma[star]); gxx = r / (r_to_star + 1.0e-30) * r / (r_to_star + 1.0e-30); psip = pow(gxx, TOV_Momentum_Psi_Power/4.); - w_lorentz_2 = 1. - gxx*TOV_Velocity_x[star]*TOV_Velocity_x[star] - - gxx*TOV_Velocity_y[star]*TOV_Velocity_y[star] - - gxx*TOV_Velocity_z[star]*TOV_Velocity_z[star]; + /* velocity components as in gr-qc/9811015 */ + w_lorentz_2 = 1./(1. - gxx*TOV_Velocity_x[star]*TOV_Velocity_x[star] + - gxx*TOV_Velocity_y[star]*TOV_Velocity_y[star] + - gxx*TOV_Velocity_z[star]*TOV_Velocity_z[star]); rho_eps = w_lorentz_2 * rho + w_lorentz_2*TOV_Gamma[star]/(TOV_Gamma[star]-1.) * press; switch(dir) @@ -266,7 +268,8 @@ *vz = mom_source_z/my_psi4/ (source/my_psi4/my_psi4+K * pow(*rhonew, Gamma)); *v_2 = my_psi4*(*vx**vx+*vy**vy+*vz**vz); - *w_lorentz_2 = 1.- *v_2; + /* velocity components as in gr-qc/9811015 */ + *w_lorentz_2 = 1./(1.- *v_2); f = my_psi4*my_psi4*( *w_lorentz_2 * *rhonew + (*w_lorentz_2*Gamma/(Gamma-1.)-1.) * @@ -316,15 +319,17 @@ my_psi4, x, y, z, rho_orig); return; } - *w_lorentz_2 = 1.- *v_2; + /* velocity components as in gr-qc/9811015 */ + *w_lorentz_2 = 1./(1.- *v_2); f = my_psi4*my_psi4*( *w_lorentz_2 * *rhonew + (*w_lorentz_2*Gamma/(Gamma-1.)-1.) * K * pow(*rhonew, Gamma)) - source; + /* d_w_lorentz_2/drhonew */ d_w_lorentz_2 = - 2. * (*w_lorentz_2 - 1.) / - (source/my_psi4/my_psi4+K*pow(*rhonew, Gamma)) * - K * Gamma * pow(*rhonew, Gamma-1.); + -2 * (*w_lorentz_2)*(*w_lorentz_2) * + (*v_2) * K*Gamma*pow(*rhonew, Gamma-1) / + (source/my_psi4/my_psi4 + K*pow(*rhonew, Gamma)); df = my_psi4*my_psi4*( d_w_lorentz_2**rhonew + *w_lorentz_2 + d_w_lorentz_2*Gamma/(Gamma-1.)* @@ -475,7 +480,8 @@ " orig_rho:%g rho:%g, v^2=%g, rel_diff=%g", \ x[i3D], y[i3D], z[i3D], rho[i3D], rhonew, v_2, fabs(rhonew-rhoold)/fabs(rhonew)); - w_lorentz_2 = 1.- v_2; + /* velocity components as in gr-qc/9811015 */ + w_lorentz_2 = 1./(1.- v_2); /* f = TOV_K[0] / (TOV_Gamma[0]-1.0) * pow(rhonew, TOV_Gamma[0]) + rhonew - source[i3D]; @@ -491,10 +497,11 @@ (w_lorentz_2*TOV_Gamma[0]/(TOV_Gamma[0]-1.)-1.) * TOV_K[0] * TOV_Gamma[0] * pow(rhonew, TOV_Gamma[0]-1.); */ + /* d_w_lorentz_2/drhonew */ d_w_lorentz_2 = - 2. * (w_lorentz_2 - 1.) / - (source[i3D]/my_psi4/my_psi4+TOV_K[0]*pow(rhonew, TOV_Gamma[0])) * - TOV_K[0] * TOV_Gamma[0] * pow(rhonew, TOV_Gamma[0]-1.); + -2 * w_lorentz_2*w_lorentz_2 * + v_2 * TOV_K[0]*TOV_Gamma[0]*pow(rhonew, TOV_Gamma[0]-1) / + (source[i3D]/my_psi4/my_psi4 + TOV_K[0]*pow(rhonew, TOV_Gamma[0])); df = my_psi4*my_psi4*( d_w_lorentz_2*rhonew + w_lorentz_2 + d_w_lorentz_2*TOV_Gamma[0]/(TOV_Gamma[0]-1.)*

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

File [modified]: param.ccl Delta lines: +3 -0 =================================================================== --- trunk/param.ccl 2010-05-03 03:53:11 UTC (rev 110) +++ trunk/param.ccl 2010-05-28 12:07:30 UTC (rev 111) @@ -63,6 +63,9 @@ *:* :: "real" } 0.0

+# contravariant fluid three velocity as measured by the Eulerian observers: v^i = u^i / (alpha u^t) + beta^i / alpha +# as used by HydroBase. Follows the Valencia formulation eg. eqs. 26 and 27 of Font et al's paper (gr-qc/9811015) or +# below Equ. 31 in http://relativity.livingreviews.org/Articles/lrr-2008-7/articlesu1.html#x6-3... CCTK_REAL TOV_Velocity_x[10] "(fixed) Velocity of neutron star, x coordinate (caution!)" STEERABLE=always { *:* :: "real"