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

1 files changed, 22 insertions, 15 deletions

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.)*