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.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/TOVSolver/trunk@111 1bdb13ef-5d69-4035-bb54-08abeb3aa7f1

2 files changed, 25 insertions, 15 deletions

diff --git a/param.ccl b/param.ccl
index 3c05d85..2fb9026 100644
--- a/param.ccl
+++ b/param.ccl
@@ -63,6 +63,9 @@ CCTK_REAL TOV_Position_z[10] "Position of neutron star, z coordinate" STEERABLE=
   *:* :: "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-30002.1
 CCTK_REAL TOV_Velocity_x[10] "(fixed) Velocity of neutron star, x coordinate (caution!)" STEERABLE=always
 {
   *:* :: "real"
diff --git a/src/external.inc b/src/external.inc
index 08aab6b..a1318f4 100644
--- a/src/external.inc
+++ b/src/external.inc
@@ -91,9 +91,10 @@ CCTK_INT TOV_Set_Rho_ADM(CCTK_POINTER_TO_CONST cctkGH,
     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 @@ CCTK_INT TOV_Set_Momentum_Source(
     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 @@ void bisection     (CCTK_REAL *vx, CCTK_REAL *vy, CCTK_REAL *vz,
         *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 @@ void newton_raphson(CCTK_REAL *vx, CCTK_REAL *vy, CCTK_REAL *vz,
                     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 @@ CCTK_INT TOV_Rescale_Sources(
                      " 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 @@ CCTK_INT TOV_Rescale_Sources(
             (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.)*