26 May
2010
26 May
'10
4:49 a.m.
Hello all,
After looking at the equations carefully yesterday, I would strongly suggest that you apply your patch to all definitions of w_lorentz_2 in external.inc. The equations there do follow the Valencia formulation as, for example, it appears in eqs. 26 and 27 of Font et al's paper (gr-qc/9811015).
Attached please find an updated version of my original patch that hopefully catches all occurrences of the Lorentz factor. Bruno, Frank: does this look ok to you?
You could add comments to the code, either to the definitions of these variables or to the overall README. Pointing to literature in the code will help future readers a lot.
I added the references and definition to param.ccl and some notes to the code.
Ok to commit?
Yours, Roland
--
My email is as private as my paper mail. I therefore support encrypting
and signing email messages. Get my PGP key from http://keys.gnupg.net.
Index: param.ccl
===================================================================
--- param.ccl (revision 110)
+++ param.ccl (working copy)
@@ -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"
Index: src/external.inc
===================================================================
--- src/external.inc (revision 110)
+++ src/external.inc (working copy)
@@ -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.);
+ (*w_lorentz_2)*(*w_lorentz_2) *
+ (-2 * (*v_2) / (source/my_psi4/my_psi4 + K*pow(*rhonew, Gamma)))*
+ K*Gamma*pow(*rhonew, Gamma-1);
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.);
+ w_lorentz_2*w_lorentz_2 *
+ (-2 * v_2 / (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);
df = my_psi4*my_psi4*(
d_w_lorentz_2*rhonew + w_lorentz_2 +
d_w_lorentz_2*TOV_Gamma[0]/(TOV_Gamma[0]-1.)*
@@ -649,3 +656,4 @@
*exported_TOV_r_1d = TOV_r_1d;
}
+