Make my old interface from the TOVSolver to TwoPunctures public

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

4 files changed, 711 insertions, 0 deletions

diff --git a/interface.ccl b/interface.ccl
index ce1f981..ccd7415 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -12,3 +12,59 @@ void FUNCTION SpatialDet(CCTK_REAL IN gxx, CCTK_REAL IN gxy, \
                          CCTK_REAL IN gyz, CCTK_REAL IN gzz, \
                          CCTK_REAL OUT det)
 USES FUNCTION SpatialDet
+
+CCTK_INT FUNCTION Set_Rho_ADM(              \
+    CCTK_POINTER_TO_CONST IN  cctkGH,       \
+    CCTK_INT              IN  size,         \
+    CCTK_REAL ARRAY       OUT source,       \
+    CCTK_REAL ARRAY       IN  x,            \
+    CCTK_REAL ARRAY       IN  y,            \
+    CCTK_REAL ARRAY       IN  z             \
+    )
+PROVIDES FUNCTION Set_Rho_ADM  \
+         WITH TOV_Set_Rho_ADM  \
+         LANGUAGE C
+
+CCTK_INT FUNCTION Set_Momentum_Source(      \
+    CCTK_POINTER_TO_CONST IN  cctkGH,       \
+    CCTK_INT              IN  dir,          \
+    CCTK_INT              IN  size,         \
+    CCTK_REAL ARRAY       OUT source,       \
+    CCTK_REAL ARRAY       IN  x,            \
+    CCTK_REAL ARRAY       IN  y,            \
+    CCTK_REAL ARRAY       IN  z             \
+    )
+PROVIDES FUNCTION Set_Momentum_Source  \
+         WITH TOV_Set_Momentum_Source  \
+         LANGUAGE C
+
+CCTK_INT FUNCTION Set_Initial_Guess_for_u(  \
+    CCTK_POINTER_TO_CONST IN  cctkGH,       \
+    CCTK_INT              IN  size,         \
+    CCTK_REAL ARRAY       OUT u,            \
+    CCTK_REAL ARRAY       IN  x,            \
+    CCTK_REAL ARRAY       IN  y,            \
+    CCTK_REAL ARRAY       IN  z             \
+    )
+PROVIDES FUNCTION Set_Initial_Guess_for_u \
+         WITH TOV_Set_Initial_Guess_for_u \
+         LANGUAGE C
+
+CCTK_INT FUNCTION Rescale_Sources(  \
+    CCTK_POINTER_TO_CONST IN  cctkGH,       \
+    CCTK_INT              IN  size,         \
+    CCTK_REAL ARRAY       IN  x,            \
+    CCTK_REAL ARRAY       IN  y,            \
+    CCTK_REAL ARRAY       IN  z,            \
+    CCTK_REAL ARRAY       IN  psi,          \
+    CCTK_REAL ARRAY       IN  gxx,          \
+    CCTK_REAL ARRAY       IN  gyy,          \
+    CCTK_REAL ARRAY       IN  gzz,          \
+    CCTK_REAL ARRAY       IN  gxy,          \
+    CCTK_REAL ARRAY       IN  gxz,          \
+    CCTK_REAL ARRAY       IN  gyz           \
+    )
+PROVIDES FUNCTION Rescale_Sources \
+         WITH TOV_Rescale_Sources \
+         LANGUAGE C
+
diff --git a/param.ccl b/param.ccl
index 55602b3..77dc77f 100644
--- a/param.ccl
+++ b/param.ccl
@@ -119,6 +119,10 @@ CCTK_INT TOV_fake_evolution "Fake evolution by setting ID at every step" STEERAB
   *:* :: "anything, 0 as off (default), everything else as on"
 } 0
 
+STRING TOV_save_to_datafile "Only save data to file and exit"
+{
+  ".*" :: "Any filename, not used if empty"
+} ""
 
 shares:HydroBase
 
diff --git a/schedule.ccl b/schedule.ccl
index f4c2446..3de374f 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -42,6 +42,15 @@ schedule TOV_C_Integrate_RHS IN TOV_Initial_Data
   OPTIONS: GLOBAL
 } "Integrate the 1d equations for the TOV star"
 
+if (!CCTK_Equals(TOV_save_to_datafile,""))
+{
+  schedule TOV_write_1D_datafile IN TOV_Initial_Data AFTER TOV_C_Integrate_RHS BEFORE TOV_C_Exact
+  {
+    LANG: C
+    OPTIONS: GLOBAL
+  } "Save data to file and exit"
+}
+
 if (CCTK_Equals(initial_data, "tov") || CCTK_Equals(initial_hydro, "tov") ||
     (TOV_Use_Old_Initial_Data > 0) || TOV_Enforce_Interpolation)
 {
diff --git a/src/external.inc b/src/external.inc
index 156c758..cc49146 100644
--- a/src/external.inc
+++ b/src/external.inc
@@ -1,6 +1,648 @@
 
+/* for debug output the location of the points */
+void TOV_debug_input_points(CCTK_INT size,
+                       const CCTK_REAL *x,
+                       const CCTK_REAL *y,
+                       const CCTK_REAL *z)
+{
+    FILE *f;
+    int i;
+    f=fopen("test.out", "w");
+    for (i=0; i<size; i++)
+      fprintf(f, "%.8g %.8g %.8g\n", x[i], y[i], z[i]);
+    fclose(f);
+}
+
+inline int get_nearest_star(const CCTK_REAL *x, const CCTK_REAL *y,
+                            const CCTK_REAL *z, const int i)
+{
+    DECLARE_CCTK_PARAMETERS
+    int star = 0, star_loop;
+    CCTK_REAL r_tmp, r_to_star=1.e40;  /* something very large */
+    for (star_loop=0; star_loop<TOV_Num_TOVs; star_loop++)
+    {
+        r_tmp = sqrt( (x[i]-TOV_Position_x[star_loop]) *
+                      (x[i]-TOV_Position_x[star_loop]) +
+                      (y[i]-TOV_Position_y[star_loop]) *
+                      (y[i]-TOV_Position_y[star_loop]) +
+                      (z[i]-TOV_Position_z[star_loop]) *
+                      (z[i]-TOV_Position_z[star_loop]) );
+        if (r_tmp < r_to_star)
+        {
+            r_to_star = r_tmp;
+            star = star_loop;
+        }
+    }
+    return star;
+}
+/* This function provides an interface for other thorns which can input a list
+ * of points (p_x, p_y, p_z) with a size of 'size' and get back the rho_adm
+ * at these points (source) */
+/* NOTE: when used with two or more NS, there might be problems if they are
+ * very close to each other, because this function does not use the fancy
+ * algorithms above to determine from which NS the data for a perticular grid
+ * point should be taken, but just uses the closest one */
+CCTK_INT TOV_Set_Rho_ADM(CCTK_POINTER_TO_CONST cctkGH,
+                         CCTK_INT  size,
+                         CCTK_REAL *source,
+                         const CCTK_REAL *x,
+                         const CCTK_REAL *y,
+                         const CCTK_REAL *z )
+{
+  DECLARE_CCTK_PARAMETERS
+  CCTK_REAL r_to_star, press, phi, r;
+  CCTK_REAL rho, eps, gxx, w_lorentz_2;
+  CCTK_INT  star_i=0, star=0;
+  CCTK_INT  i;
+
+  assert(TOV_Surface!=0);
+  assert(TOV_R_Surface!=0);
+  assert(TOV_Atmosphere!=0);
+
+  assert(TOV_r_1d!=0);
+  assert(TOV_rbar_1d!=0);
+  assert(TOV_press_1d!=0);
+  assert(TOV_phi_1d!=0);
+  assert(TOV_m_1d!=0);
+
+  /* CCTK_INFO("Request for Sources from TOVSolverC"); */
+  /* TOV_debug_input_points(size,x,y,z); */
+  /* loop over all points we got */
+
+  for (i=0; i<size; i++)
+  {
+    /* calculate the distance to the star */
+    star=get_nearest_star(x,y,z,i);
+    star_i = star * TOV_Num_Radial;
+    r_to_star = sqrt( (x[i]-TOV_Position_x[star]) *
+                      (x[i]-TOV_Position_x[star]) +
+                      (y[i]-TOV_Position_y[star]) *
+                      (y[i]-TOV_Position_y[star]) +
+                      (z[i]-TOV_Position_z[star]) *
+                      (z[i]-TOV_Position_z[star]) );
+    /* call the interpolator */
+    TOV_C_interp_tov_isotropic(star,
+                               &(TOV_press_1d[star_i]), &(TOV_phi_1d[star_i]),
+                               &(TOV_rbar_1d[star_i]),  &(TOV_r_1d[star_i]),
+                               &r_to_star, TOV_Surface[star],
+                               &press, &phi, &r,
+                               0, 1, NULL);
+    press = (press > 0.0) ? press : 0.0;
+    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];
+    if (rho > 0.0)
+      /* source[i]=gxx*gxx*rho*(1.0+eps); */
+      source[i]=gxx*gxx* (w_lorentz_2*rho +
+                          (w_lorentz_2*TOV_Gamma[star]/(TOV_Gamma[star]-1.)-1.)*
+                          press);
+    else
+      source[i]=0.0;
+  }
+  return 1;
+}
+
+/* This function provides an interface for other thorns which can input a list
+ * of points (p_x, p_y, p_z) with a size of 'size' and get back the momentum
+ * sources at these points (source) */
+/* NOTE: when used with two or more NS, there might be problems if they are
+ * very close to each other, because this function does not use the fancy
+ * algorithms above to determine from which NS the data for a perticular grid
+ * point should be taken, but just uses the closest one */
+CCTK_INT TOV_Set_Momentum_Source(
+                         CCTK_POINTER_TO_CONST cctkGH,
+                         CCTK_INT  dir,
+                         CCTK_INT  size,
+                         CCTK_REAL *source,
+                         const CCTK_REAL *x,
+                         const CCTK_REAL *y,
+                         const CCTK_REAL *z )
+{
+  DECLARE_CCTK_PARAMETERS
+  CCTK_REAL r_to_star, press, phi, r, gxx, rho, rho_eps, psip, w_lorentz_2;
+  CCTK_INT  star_i=0, star=0;
+  CCTK_INT  i;
+
+  assert(TOV_Surface!=0);
+  assert(TOV_R_Surface!=0);
+  assert(TOV_Atmosphere!=0);
+
+  assert(TOV_r_1d!=0);
+  assert(TOV_rbar_1d!=0);
+  assert(TOV_press_1d!=0);
+  assert(TOV_phi_1d!=0);
+  assert(TOV_m_1d!=0);
+
+  /* CCTK_INFO("Request for Momentum-Sources from TOVSolverC"); */
+  /* loop over all points we got */
+
+  for (i=0; i<size; i++)
+  {
+    /* calculate the distance to the star */
+    star=get_nearest_star(x,y,z,i);
+    star_i = star * TOV_Num_Radial;
+    r_to_star = sqrt( (x[i]-TOV_Position_x[star]) *
+                      (x[i]-TOV_Position_x[star]) +
+                      (y[i]-TOV_Position_y[star]) *
+                      (y[i]-TOV_Position_y[star]) +
+                      (z[i]-TOV_Position_z[star]) *
+                      (z[i]-TOV_Position_z[star]) );
+    /* call the interpolator */
+    TOV_C_interp_tov_isotropic(star,
+                               &(TOV_press_1d[star_i]), &(TOV_phi_1d[star_i]),
+                               &(TOV_rbar_1d[star_i]),  &(TOV_r_1d[star_i]),
+                               &r_to_star, TOV_Surface[star],
+                               &press, &phi, &r,
+                               0, 1, NULL);
+    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];
+    rho_eps = w_lorentz_2 * rho +
+              w_lorentz_2*TOV_Gamma[star]/(TOV_Gamma[star]-1.) * press;
+    switch(dir)
+    {
+      case 0: source[i]=psip*rho_eps*gxx*TOV_Velocity_x[star]; break;
+      case 1: source[i]=psip*rho_eps*gxx*TOV_Velocity_y[star]; break;
+      case 2: source[i]=psip*rho_eps*gxx*TOV_Velocity_z[star]; break;
+      default: CCTK_WARN(0,
+                         "dir has to be in [0:2] in TOV_Set_Momentum_Source");
+    }
+    /* Only set velocity where matter is */
+    if (press<1.e-15)
+      source[i]=0.0;
+  }
+  return 1;
+}
+
+/* NOTE: when used with two or more NS, there might be problems if they are
+ * very close to each other, because this function does not use the fancy
+ * algorithms above to determine from which NS the data for a perticular grid
+ * point should be taken, but just uses the closest one */
+CCTK_INT TOV_Set_Initial_Guess_for_u(
+                         CCTK_POINTER_TO_CONST cctkGH,
+                         CCTK_INT  size,
+                         CCTK_REAL *u,
+                         const CCTK_REAL *x,
+                         const CCTK_REAL *y,
+                         const CCTK_REAL *z )
+{
+  DECLARE_CCTK_PARAMETERS
+  CCTK_REAL r_to_star, press, phi, r;
+  CCTK_INT  star_i=0, star=0;
+  CCTK_INT  i;
+
+  assert(TOV_Surface!=0);
+  assert(TOV_R_Surface!=0);
+  assert(TOV_Atmosphere!=0);
+
+  assert(TOV_r_1d!=0);
+  assert(TOV_rbar_1d!=0);
+  assert(TOV_press_1d!=0);
+  assert(TOV_phi_1d!=0);
+  assert(TOV_m_1d!=0);
+
+  CCTK_INFO("Using initial guess from TOVSolverC");
+  /* TOV_debug_input_points(size,x,y,z); */
+
+  /* loop over all points we got */
+  for (i=0; i<size; i++)
+  {
+    CCTK_REAL rho, eps, gxx, r_tmp;
+    int star_loop;
+    /* calculate the distance to the stars and choose the closest */
+    star=get_nearest_star(x,y,z,i);
+    star_i = star * TOV_Num_Radial;
+    r_to_star = sqrt( (x[i]-TOV_Position_x[star]) *
+                      (x[i]-TOV_Position_x[star]) +
+                      (y[i]-TOV_Position_y[star]) *
+                      (y[i]-TOV_Position_y[star]) +
+                      (z[i]-TOV_Position_z[star]) *
+                      (z[i]-TOV_Position_z[star]) );
+    /* call the interpolator */
+    TOV_C_interp_tov_isotropic(star,
+                               &(TOV_press_1d[star_i]), &(TOV_phi_1d[star_i]),
+                               &(TOV_rbar_1d[star_i]),  &(TOV_r_1d[star_i]),
+                               &r_to_star, TOV_Surface[star],
+                               &press, &phi, &r,
+                               0, 1, NULL);
+    rho = pow(press/TOV_K[star], 1.0/TOV_Gamma[star]);
+    if (rho>0.0)
+      eps = press/(TOV_Gamma[star] - 1.0) / rho;
+    else
+      eps = 0.0;
+    gxx = r / (r_to_star + 1.0e-30) * r / (r_to_star + 1.0e-30);
+    u[i]=pow(gxx, 0.25) - 1.0;
+  }
+  return 1;
+}
+
+
+void bisection     (CCTK_REAL *vx, CCTK_REAL *vy, CCTK_REAL *vz,
+                    CCTK_REAL *rhoold, CCTK_REAL *rhonew,
+                    CCTK_REAL *w_lorentz_2, CCTK_REAL *v_2,
+                    CCTK_REAL Gamma, CCTK_REAL K,
+                    CCTK_REAL source,
+                    CCTK_REAL mom_source_x,
+                    CCTK_REAL mom_source_y,
+                    CCTK_REAL mom_source_z,
+                    CCTK_REAL my_psi4,
+                    CCTK_REAL x, CCTK_REAL y, CCTK_REAL z, CCTK_REAL rho_orig)
+{
+    CCTK_REAL lb, ub, f;
+    lb=0.0;
+    ub=2*rho_orig;
+    do
+    {
+        *rhonew = (ub-lb)/2.;
+        *vx = mom_source_x/my_psi4/
+              (source/my_psi4/my_psi4+K * pow(*rhonew, Gamma));
+        *vy = mom_source_y/my_psi4/
+              (source/my_psi4/my_psi4+K * pow(*rhonew, Gamma));
+        *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;
+        f = my_psi4*my_psi4*(
+            *w_lorentz_2 * *rhonew +
+            (*w_lorentz_2*Gamma/(Gamma-1.)-1.) *
+            K * pow(*rhonew, Gamma)) - source;
+        if (f < 0.0)
+          ub = *rhonew;
+        else
+          lb = *rhonew;
+        *rhoold = *rhonew;
+    } while (ub-lb < 1.e-10);
+}
+
+void newton_raphson(CCTK_REAL *vx, CCTK_REAL *vy, CCTK_REAL *vz,
+                    CCTK_REAL *rhoold, CCTK_REAL *rhonew,
+                    CCTK_REAL *w_lorentz_2, CCTK_REAL *v_2,
+                    CCTK_REAL Gamma, CCTK_REAL K,
+                    CCTK_REAL source,
+                    CCTK_REAL mom_source_x,
+                    CCTK_REAL mom_source_y,
+                    CCTK_REAL mom_source_z,
+                    CCTK_REAL my_psi4,
+                    CCTK_REAL x, CCTK_REAL y, CCTK_REAL z, CCTK_REAL rho_orig)
+{
+    int count=0;
+    CCTK_REAL d_w_lorentz_2, f, df;
+      do
+      {
+        count++;
+        *vx = mom_source_x/my_psi4/
+              (source/my_psi4/my_psi4+K * pow(*rhonew, Gamma));
+        *vy = mom_source_y/my_psi4/
+              (source/my_psi4/my_psi4+K * pow(*rhonew, Gamma));
+        *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);
+        if (count > 100)
+        {
+          CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING, \
+                     "Failed to converge IVP rescaling\nat (%.3g %.3g %.3g)\n"\
+                     " orig_rho:%g rho:%g, v^2=%g, rel_diff=%g", \
+                     x, y, z, rho_orig, *rhonew, *v_2,
+                     fabs(*rhonew-*rhoold)/fabs(*rhonew));
+          CCTK_WARN(1, "trying slow bisection method");
+          bisection(vx,vy,vz,rhoold,rhonew,w_lorentz_2,v_2,
+                    Gamma, K,
+                    source, mom_source_x, mom_source_y, mom_source_z,
+                    my_psi4, x, y, z, rho_orig);
+          return;
+        }
+        *w_lorentz_2 = 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 =
+             2. * (*w_lorentz_2 - 1.) /
+             (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.)*
+              K * pow(*rhonew, Gamma) +
+             (*w_lorentz_2*Gamma/(Gamma-1.)-1.) *
+               K * Gamma * pow(*rhonew, Gamma-1.));
+        *rhoold = *rhonew;
+        *rhonew = *rhoold - f/df;
+      }
+      while ( (fabs(*rhonew - *rhoold)/fabs(*rhonew) > 1.e-12) &&
+              (fabs(*rhonew - *rhoold) > 1.e-12)
+            );
+}
+
+CCTK_INT TOV_Rescale_Sources(
+    CCTK_POINTER_TO_CONST cctkGH,
+    CCTK_INT              size,
+    const CCTK_REAL       *x,
+    const CCTK_REAL       *y,
+    const CCTK_REAL       *z,
+    const CCTK_REAL       *psi,
+    const CCTK_REAL       *gxx,
+    const CCTK_REAL       *gyy,
+    const CCTK_REAL       *gzz,
+    const CCTK_REAL       *gxy,
+    const CCTK_REAL       *gxz,
+    const CCTK_REAL       *gyz)
+{
+  DECLARE_CCTK_PARAMETERS
+  CCTK_REAL *rho, *press, *eps, *dens, *tau, *w_lorentz, *vel0, *vel1, *vel2,
+            *scon0, *scon1, *scon2;
+  CCTK_REAL *rho_p, *press_p, *eps_p, *dens_p, *tau_p, *w_lorentz_p, *vel0_p,
+            *vel1_p, *vel2_p, *scon0_p, *scon1_p, *scon2_p;
+  CCTK_REAL *rho_p_p, *press_p_p, *eps_p_p, *dens_p_p, *tau_p_p, *w_lorentz_p_p,
+            *vel0_p_p, *vel1_p_p, *vel2_p_p, *scon0_p_p, *scon1_p_p, *scon2_p_p;
+  CCTK_REAL *u, *source, *mom_source[3];
+  CCTK_REAL psip, my_psi4, sqrt_det, vx, vy, vz, v_2,
+            w_lorentz_2, d_w_lorentz_2, D_h_w;
+  CCTK_REAL vlowx, vlowy, vlowz;
+  CCTK_INT  i, type;
+  CCTK_REAL rhoold, rhonew, f, df;
+  int i3D;
+  static int debug = 1;
+  FILE *debugfile = NULL;
+
+  CCTK_INFO("Rescaling Sources");
+
+  /* get Whisky variables */
+  rho      =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 0, "HydroBase::rho");
+  press    =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 0, "HydroBase::press");
+  eps      =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 0, "HydroBase::eps");
+  dens     =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 0, "Whisky::dens");
+  tau      =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 0, "Whisky::tau");
+  w_lorentz=(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 0, "Whisky::w_lorentz");
+  vel0     =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 0, "HydroBase::vel[0]");
+  vel1     =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 0, "HydroBase::vel[1]");
+  vel2     =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 0, "HydroBase::vel[2]");
+  scon0       =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 0, "Whisky::scon[0]");
+  scon1       =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 0, "Whisky::scon[1]");
+  scon2       =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 0, "Whisky::scon[2]");
+
+  rho_p    =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 1, "HydroBase::rho");
+  press_p  =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 1, "HydroBase::press");
+  eps_p    =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 1, "HydroBase::eps");
+  dens_p   =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 1, "Whisky::dens");
+  tau_p    =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 1, "Whisky::tau");
+  w_lorentz_p=(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 1, "Whisky::w_lorentz");
+  vel0_p   =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 1, "HydroBase::vel[0]");
+  vel1_p   =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 1, "HydroBase::vel[1]");
+  vel2_p   =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 1, "HydroBase::vel[2]");
+  scon0_p     =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 1, "Whisky::scon[0]");
+  scon1_p     =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 1, "Whisky::scon[1]");
+  scon2_p     =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 1, "Whisky::scon[2]");
+
+  rho_p_p  =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 2, "HydroBase::rho");
+  press_p_p=(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 2, "HydroBase::press");
+  eps_p_p  =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 2, "HydroBase::eps");
+  dens_p_p =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 2, "Whisky::dens");
+  tau_p_p  =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 2, "Whisky::tau");
+  w_lorentz_p_p=(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 2, "Whisky::w_lorentz");
+  vel0_p_p =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 2, "HydroBase::vel[0]");
+  vel1_p_p =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 2, "HydroBase::vel[1]");
+  vel2_p_p =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 2, "HydroBase::vel[2]");
+  scon0_p_p   =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 2, "Whisky::scon[0]");
+  scon1_p_p   =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 2, "Whisky::scon[1]");
+  scon2_p_p   =(CCTK_REAL*)CCTK_VarDataPtr(cctkGH, 2, "Whisky::scon[2]");
+
+  /* get u */
+  u=calloc(size, sizeof(CCTK_REAL));
+  TOV_Set_Initial_Guess_for_u(cctkGH, size, u, x, y, z);
+
+  /* get the old source */
+  source=calloc(size, sizeof(CCTK_REAL));
+  TOV_Set_Rho_ADM(cctkGH, size, source, x, y, z);
+  mom_source[0]=calloc(size, sizeof(CCTK_REAL));
+  mom_source[1]=calloc(size, sizeof(CCTK_REAL));
+  mom_source[2]=calloc(size, sizeof(CCTK_REAL));
+  TOV_Set_Momentum_Source(cctkGH, 0, size, mom_source[0], x, y, z);
+  TOV_Set_Momentum_Source(cctkGH, 1, size, mom_source[1], x, y, z);
+  TOV_Set_Momentum_Source(cctkGH, 2, size, mom_source[2], x, y, z);
+
+  if (debug)
+  {
+    debugfile=fopen("tovdebug.dat", "w");
+  }
+  /* loop over all points we got */
+  for (i3D=0; i3D<size; i3D++)
+  {
+    if (source[i3D]>0.0)
+    {
+      int count=0;
+      if (debug && (fabs(z[i3D])<1.e-15))
+      {
+        fprintf(debugfile,
+                "%.15g %.15g %.15g %.15g %.15g %.15g ",
+                x[i3D], y[i3D], source[i3D],
+                mom_source[0][i3D], mom_source[1][i3D], mom_source[1][i3D]);
+      }
+      if (psi != NULL)
+        my_psi4 = pow(psi[i3D], 4.0)*gxx[i3D];
+      else
+        my_psi4 = gxx[i3D];
+      /* we do not need u, we need gxx here, which will be stored in u! */
+      u[i3D] = pow(u[i3D]+1.0, 4.0);
+      /* We also store the original J^i in mom_source[i] */
+      psip = pow(u[i3D], TOV_Momentum_Psi_Power/4.);
+      mom_source[0][i3D] /= psip;
+      mom_source[1][i3D] /= psip;
+      mom_source[2][i3D] /= psip;
+
+      if (rho[i3D] < 1.e-10)
+          rho[i3D] = 1.e-10;
+      rhoold = rhonew = rho[i3D];
+#ifdef OLD_NR
+      do
+      {
+        count++;
+        vx = mom_source[0][i3D]/my_psi4/
+             (source[i3D]/my_psi4/my_psi4+TOV_K[0] * pow(rhonew, TOV_Gamma[0]));
+        vy = mom_source[1][i3D]/my_psi4/
+             (source[i3D]/my_psi4/my_psi4+TOV_K[0] * pow(rhonew, TOV_Gamma[0]));
+        vz = mom_source[2][i3D]/my_psi4/
+             (source[i3D]/my_psi4/my_psi4+TOV_K[0] * pow(rhonew, TOV_Gamma[0]));
+        v_2 = my_psi4*(vx*vx+vy*vy+vz*vz);
+        if (count > 100)
+          CCTK_VWarn(count<110, __LINE__, __FILE__, CCTK_THORNSTRING, \
+                     "Failed to converge IVP rescaling\nat (%.3g %.3g %.3g)\n"\
+                     " 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;
+        /*
+        f = TOV_K[0] / (TOV_Gamma[0]-1.0) * pow(rhonew, TOV_Gamma[0]) +
+            rhonew - source[i3D];
+        df = TOV_K[0] * TOV_Gamma[0] / (TOV_Gamma[0]-1.0) *
+             pow(rhonew, TOV_Gamma[0]-1.0) + 1.0;
+        */
+        f = my_psi4*my_psi4*(
+            w_lorentz_2 * rhonew +
+            (w_lorentz_2*TOV_Gamma[0]/(TOV_Gamma[0]-1.)-1.) *
+            TOV_K[0] * pow(rhonew, TOV_Gamma[0])) - source[i3D];
+        /*
+        df= w_lorentz_2 +
+            (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 =
+             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.);
+        df = my_psi4*my_psi4*(
+             d_w_lorentz_2*rhonew + w_lorentz_2 +
+             d_w_lorentz_2*TOV_Gamma[0]/(TOV_Gamma[0]-1.)*
+              TOV_K[0] * pow(rhonew, TOV_Gamma[0]) +
+             (w_lorentz_2*TOV_Gamma[0]/(TOV_Gamma[0]-1.)-1.) *
+               TOV_K[0] * TOV_Gamma[0] * pow(rhonew, TOV_Gamma[0]-1.));
+        /*
+        df= 1. + my_psi4*v_2/(rhonew*rhonew) +
+            ((TOV_Gamma[0]+(2.-TOV_Gamma[0])*my_psi4*v_2/(rhonew*rhonew))
+             / (TOV_Gamma[0]-1.) - 1.)
+            * TOV_K[0] * TOV_Gamma[0] * pow(rhonew, TOV_Gamma[0]-1.);
+        */
+
+        rhoold = rhonew;
+        rhonew = rhoold - f/df;
+      }
+      while ( (fabs(rhonew-rhoold)/fabs(rhonew) > 1.e-12) &&
+              (fabs(rhonew - rhoold) > 1.e-12)
+            );
+      if (count==1)
+          printf("Fast rescale! %g %g %g %g\n", u[i3D], my_psi4,
+                 mom_source[0][i3D]/u[i3D]/
+                  (source[i3D]/u[i3D]/u[i3D]+TOV_K[0]*
+                  pow(rhonew, TOV_Gamma[0])),
+                 mom_source[0][i3D]/my_psi4/
+                  (source[i3D]/u[i3D]/my_psi4+TOV_K[0]*
+                  pow(rhonew, TOV_Gamma[0])));
+#else
+      newton_raphson(&vx,&vy,&vz,&rhoold,&rhonew,&w_lorentz_2,&v_2,
+                     TOV_Gamma[0], TOV_K[0],
+                     source[i3D],
+                     mom_source[0][i3D],mom_source[1][i3D],mom_source[2][i3D],
+                     my_psi4, x[i3D], y[i3D], z[i3D], rho[i3D]);
+#endif
+      rho[i3D]   = rhonew;
+      press[i3D] = TOV_K[0] * pow(rhonew, TOV_Gamma[0]);
+      eps[i3D]   = TOV_K[0] * pow(rhonew, TOV_Gamma[0] - 1.0) /
+                              (TOV_Gamma[0] - 1.0);
+      sqrt_det   = pow(my_psi4, 1.5);
+      w_lorentz[i3D] = sqrt(w_lorentz_2);
+      dens[i3D] = sqrt_det * w_lorentz[i3D] * rho[i3D];
+      tau[i3D]  = sqrt_det * (rho[i3D]  *w_lorentz[i3D]*w_lorentz[i3D] *
+                                         (1.0 + eps[i3D]) +
+                              press[i3D]*(w_lorentz[i3D]*w_lorentz[i3D]-1.0) ) -
+                  dens[i3D];
+      vel0[i3D] = vx;
+      vel1[i3D] = vy;
+      vel2[i3D] = vz;
+      D_h_w = dens[i3D] * (1 + eps[i3D] + press[i3D]/rho[i3D]) * w_lorentz[i3D];
+
+      vlowx = gxx[i3D]*vel0[i3D] + gxy[i3D]*vel1[i3D] + gxz[i3D]*vel2[i3D];
+      vlowy = gxy[i3D]*vel0[i3D] + gyy[i3D]*vel1[i3D] + gyz[i3D]*vel2[i3D];
+      vlowz = gxz[i3D]*vel0[i3D] + gyz[i3D]*vel1[i3D] + gzz[i3D]*vel2[i3D];
+
+      scon0[i3D] = D_h_w * vlowx;
+      scon1[i3D] = D_h_w * vlowy;
+      scon2[i3D] = D_h_w * vlowz;
+
+      if (debug && (fabs(z[i3D])<1.e-15))
+      {
+        fprintf(debugfile,
+                "%.15g %.15g %.15g %.15g\n", /*7*/
+                gxx[i3D]*gxx[i3D]*( /* ham - source*/
+                  w_lorentz_2*rho[i3D]+
+                  (w_lorentz_2*TOV_Gamma[0]/(TOV_Gamma[0]-1.)-1.)*press[i3D]),
+                gxx[i3D]*gxx[i3D]*( /* ham - source */
+                  w_lorentz_2*(rho[i3D]*(1.+eps[i3D])+press[i3D])-press[i3D]),
+                (w_lorentz_2 * rho[i3D] + /* momx - source */
+                 w_lorentz_2*TOV_Gamma[0]/(TOV_Gamma[0]-1.) * press[i3D])*
+                 gxx[i3D]*vel0[i3D],
+                w_lorentz_2 * (rho[i3D]* /* momx - source */
+                 (1.+eps[i3D])+ press[i3D])*
+                 gxx[i3D]*vel0[i3D]
+                );
+      }
+
+      /* FIXME: check for Atmosphere */
+    }
+  }
+  if (debug)
+    fclose(debugfile);
+
+  switch(*(const CCTK_INT*)CCTK_ParameterGet("TOV_Populate_Timelevels",
+                                     "Whisky_TOVSolverC", &type))
+  {
+    case 3:
+        TOV_Copy(size, rho_p_p,  rho);
+        TOV_Copy(size, eps_p_p,  eps);
+        TOV_Copy(size, vel0_p_p, vel0);
+        TOV_Copy(size, vel1_p_p, vel1);
+        TOV_Copy(size, vel2_p_p, vel2);
+        TOV_Copy(size, dens_p_p, dens);
+        TOV_Copy(size, tau_p_p,  tau);
+        TOV_Copy(size, scon0_p_p,   scon0);
+        TOV_Copy(size, scon1_p_p,   scon1);
+        TOV_Copy(size, scon2_p_p,   scon2);
+        TOV_Copy(size, w_lorentz_p_p, w_lorentz);
+    case 2:
+        TOV_Copy(size, rho_p,  rho);
+        TOV_Copy(size, eps_p,  eps);
+        TOV_Copy(size, vel0_p, vel0);
+        TOV_Copy(size, vel1_p, vel1);
+        TOV_Copy(size, vel2_p, vel2);
+        TOV_Copy(size, dens_p, dens);
+        TOV_Copy(size, tau_p,  tau);
+        TOV_Copy(size, scon0_p,   scon0);
+        TOV_Copy(size, scon1_p,   scon1);
+        TOV_Copy(size, scon2_p,   scon2);
+        TOV_Copy(size, w_lorentz_p, w_lorentz);
+  }
+  free(source);
+  free(u);
+  return 1;
+}
+
+void TOV_write_1D_datafile(CCTK_ARGUMENTS)
+{
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_PARAMETERS
+  int i;
+
+  if (TOV_Num_TOVs > 1)
+    CCTK_WARN(0, "Writing a data file for multiple stars is not (yet) "
+                 "supported");
+  FILE *file;
+  file = fopen(TOV_save_to_datafile, "w");
+  fprintf(file, "TOVSolver data file\n");
+  fprintf(file, "version 1.0\n");
+  fprintf(file, "TOV_Num_Radial %d\n", TOV_Num_Radial);
+  fprintf(file, "\n");
+  for (i=0; i < TOV_Num_Radial; i++)
+  {
+      fprintf(file, "%g %g %g %g %g\n",
+              TOV_rbar_1d[i],
+              TOV_r_1d[i],
+              TOV_m_1d[i],
+              TOV_phi_1d[i],
+              TOV_press_1d[i]);
+  }
+  fclose(file);
+  CCTK_WARN(0, "Simulation stopped as requested after writing data to file");
+}
+
 void TOVSolverC_export_local_variables(CCTK_REAL **exported_TOV_press_1d, CCTK_REAL **exported_TOV_m_1d, CCTK_REAL **exported_TOV_phi_1d, CCTK_REAL **exported_TOV_rbar_1d, CCTK_REAL **exported_TOV_r_1d)
 {
+  int i;
+
   *exported_TOV_press_1d = TOV_press_1d;
   *exported_TOV_m_1d     = TOV_m_1d;
   *exported_TOV_phi_1d   = TOV_phi_1d;