Major cleanup and fixes. The current version is able to read in example Meudon data and converges in both hamiltonian and momentum constraint on the Cactus grid until it hits the spectral residuum

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/Meudon_Bin_NS/trunk@9 f2ea251b-07d6-4a45-8f2c-a162b3fa7596

8 files changed, 137 insertions, 178 deletions

diff --git a/README b/README
index 8191c71..e62409b 100644
--- a/README
+++ b/README
@@ -1,6 +1,7 @@
-Cactus Code Thorn ID_Bin_NS
-Author(s)    : Erik Schnetter <schnetter@cct.lsu.edu>
-Maintainer(s): Erik Schnetter <schnetter@cct.lsu.edu>
+Cactus Code Thorn Meudon_Bin_NS
+Author(s)    : Erik Schnetter
+               Frank Löffler
+Maintainer(s): Einstein Toolkit Maintainers
 Licence      : GPL
 --------------------------------------------------------------------------
 
@@ -8,5 +9,5 @@ Licence      : GPL
 
 Import LORENE Bin_NS binary neutron star initial data.
 
-This thorn is one of three closely related thorns ID_Bin_BH,
-ID_Bin_NS, and ID_Mag_NS which all import LORENE initial data.
+This thorn is one of three closely related thorns Meudon_Bin_BH,
+Meudon_Bin_NS, and Meudon_Mag_NS which all import LORENE initial data.
diff --git a/configuration.ccl b/configuration.ccl
index ba8ae64..98a01db 100644
--- a/configuration.ccl
+++ b/configuration.ccl
@@ -1,3 +1,3 @@
-# Configuration definition for thorn ID_Bin_NS
+# Configuration definition for thorn Meudon_Bin_NS
 
 REQUIRES LORENE
diff --git a/interface.ccl b/interface.ccl
index e4ac277..3b75a5a 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -1,6 +1,6 @@
-# Interface definition for thorn ID_Bin_NS
+# Interface definition for thorn Meudon_Bin_NS
 
-IMPLEMENTS: ID_Bin_NS
+IMPLEMENTS: Meudon_Bin_NS
 
 INHERITS: grid SummationByParts ADMBase HydroBase
 
diff --git a/param.ccl b/param.ccl
index d3c1588..d7ac849 100644
--- a/param.ccl
+++ b/param.ccl
@@ -1,30 +1,30 @@
-# Parameter definitions for thorn ID_Bin_NS
+# Parameter definitions for thorn Meudon_Bin_NS
 
 SHARES: ADMBase
 
 EXTENDS KEYWORD initial_data
 {
-  "ID_Bin_NS" :: ""
+  "Meudon_Bin_NS" :: ""
 }
 
 EXTENDS KEYWORD initial_lapse
 {
-  "ID_Bin_NS" :: ""
+  "Meudon_Bin_NS" :: ""
 }
 
 EXTENDS KEYWORD initial_shift
 {
-  "ID_Bin_NS" :: ""
+  "Meudon_Bin_NS" :: ""
 }
 
 EXTENDS KEYWORD initial_dtlapse
 {
-  "ID_Bin_NS" :: ""
+  "Meudon_Bin_NS" :: ""
 }
 
 EXTENDS KEYWORD initial_dtshift
 {
-  "ID_Bin_NS" :: ""
+  "Meudon_Bin_NS" :: ""
 }
 
 
@@ -33,7 +33,7 @@ SHARES: HydroBase
 
 EXTENDS KEYWORD initial_hydro
 {
-  "ID_Bin_NS" :: ""
+  "Meudon_Bin_NS" :: ""
 }
 
 
diff --git a/schedule.ccl b/schedule.ccl
index a174c93..64164f8 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -1,18 +1,13 @@
-# Schedule definitions for thorn ID_Bin_NS
+# Schedule definitions for thorn Meudon_Bin_NS
 
-if (CCTK_EQUALS (initial_data,    "ID_Bin_NS") ||
-    CCTK_EQUALS (initial_lapse,   "ID_Bin_NS") ||
-    CCTK_EQUALS (initial_shift,   "ID_Bin_NS") ||
-    CCTK_EQUALS (initial_dtlapse, "ID_Bin_NS") ||
-    CCTK_EQUALS (initial_dtshift, "ID_Bin_NS") ||
-    CCTK_EQUALS (initial_hydro,   "ID_Bin_NS"))
+if (CCTK_EQUALS (initial_data,    "Meudon_Bin_NS") ||
+    CCTK_EQUALS (initial_lapse,   "Meudon_Bin_NS") ||
+    CCTK_EQUALS (initial_shift,   "Meudon_Bin_NS") ||
+    CCTK_EQUALS (initial_dtlapse, "Meudon_Bin_NS") ||
+    CCTK_EQUALS (initial_dtshift, "Meudon_Bin_NS") ||
+    CCTK_EQUALS (initial_hydro,   "Meudon_Bin_NS"))
 {
-  SCHEDULE ID_Bin_NS_check_parameters AT paramcheck
-  {
-    LANG: C
-  } "Check parameters"
-  
-  SCHEDULE ID_Bin_NS_initialise IN HydroBase_Initial
+  SCHEDULE Meudon_Bin_NS_initialise IN HydroBase_Initial
   {
     LANG: C
   } "Set up binary neutron star initial data"
diff --git a/src/Bin_NS.cc b/src/Bin_NS.cc
index 66dd9df..64fb2d9 100644
--- a/src/Bin_NS.cc
+++ b/src/Bin_NS.cc
@@ -1,3 +1,6 @@
+/* (c) 2009 Erik Schnetter
+ * (c) 2010 Frank Loeffler */
+
 #include <cstdio>
 #include <cassert>
 #include <vector>
@@ -7,23 +10,23 @@
 #include <cctk_Parameters.h>
 
 #include <bin_ns.h>
+#include <unites.h>
 
 using namespace std;
 
-
 static void set_dt_from_domega (CCTK_ARGUMENTS,
                                 CCTK_REAL const* const var,
                                 CCTK_REAL      * const dtvar,
                                 CCTK_REAL const& omega)
 {
   DECLARE_CCTK_ARGUMENTS;
-  
+
   int const npoints = cctk_lsh[0] * cctk_lsh[1] * cctk_lsh[2];
   vector<CCTK_REAL> dxvar(npoints), dyvar(npoints);
-  
+
   Diff_gv (cctkGH, 0, var, &dxvar[0], -1);
   Diff_gv (cctkGH, 1, var, &dyvar[0], -1);
-  
+
 #pragma omp parallel for
   for (int i=0; i<npoints; ++i) {
     CCTK_REAL const ephix = +y[i];
@@ -33,59 +36,52 @@ static void set_dt_from_domega (CCTK_ARGUMENTS,
   }
 }
 
-
-
 extern "C"
-void ID_Bin_NS_initialise (CCTK_ARGUMENTS)
+void Meudon_Bin_NS_initialise (CCTK_ARGUMENTS)
 {
   DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
-  
+
   CCTK_INFO ("Setting up LORENE Bin_NS initial data");
-  
+
   // Meudon data are distributed in SI units (MKSA).  Here are some
   // conversion factors.
-  
-  // Defined constants
-  CCTK_REAL const c_light = 299792458.0; // speed of light [m/s]
-  
-  // Constants of nature (IAU, CODATA):
-  CCTK_REAL const G_grav = 6.67428e-11; // gravitational constant [m^3/kg/s^2]
-  CCTK_REAL const M_sun  = 1.98892e+30; // solar mass [kg]
-  
+  // Be aware: these are the constants Lorene uses. They do differ from other
+  // conventions, but they gave the best results in some tests.
+
+  CCTK_REAL const c_light  = Unites::c_si;      // speed of light [m/s]
+  CCTK_REAL const nuc_dens = Unites::rhonuc_si; // Nuclear density as used in Lorene units [kg/m^3]
+  CCTK_REAL const G_grav   = Unites::g_si;      // gravitational constant [m^3/kg/s^2]
+  CCTK_REAL const M_sun    = Unites::msol_si;   // solar mass [kg]
+
   // Cactus units in terms of SI units:
   // (These are derived from M = M_sun, c = G = 1, and using 1/M_sun
   // for the magnetic field)
   CCTK_REAL const cactusM = M_sun;
   CCTK_REAL const cactusL = cactusM * G_grav / pow(c_light,2);
   CCTK_REAL const cactusT = cactusL / c_light;
-  
+
   // Other quantities in terms of Cactus units
-  CCTK_REAL const coord_unit = cactusL / 1.0e+3;         // from km
+  CCTK_REAL const coord_unit = cactusL / 1.0e+3;         // from km (~1.477)
   CCTK_REAL const rho_unit   = cactusM / pow(cactusL,3); // from kg/m^3
-  CCTK_REAL const ener_unit  = pow(cactusL,2);           // from c^2
-  CCTK_REAL const vel_unit   = cactusL / cactusT / c_light; // from c
-  
-  
-  
+
+
   CCTK_INFO ("Setting up coordinates");
-  
+
   int const npoints = cctk_lsh[0] * cctk_lsh[1] * cctk_lsh[2];
   vector<double> xx(npoints), yy(npoints), zz(npoints);
-  
+
 #pragma omp parallel for
   for (int i=0; i<npoints; ++i) {
     xx[i] = x[i] * coord_unit;
     yy[i] = y[i] * coord_unit;
     zz[i] = z[i] * coord_unit;
   }
-  
-  
-  
+
   CCTK_VInfo (CCTK_THORNSTRING, "Reading from file \"%s\"", filename);
-  
+
   Bin_NS bin_ns (npoints, &xx[0], &yy[0], &zz[0], filename);
-  
+
   CCTK_VInfo (CCTK_THORNSTRING, "omega [rad/s]:       %g", bin_ns.omega);
   CCTK_VInfo (CCTK_THORNSTRING, "dist [km]:           %g", bin_ns.dist);
   CCTK_VInfo (CCTK_THORNSTRING, "dist_mass [km]:      %g", bin_ns.dist_mass);
@@ -101,113 +97,104 @@ void ID_Bin_NS_initialise (CCTK_ARGUMENTS)
   CCTK_VInfo (CCTK_THORNSTRING, "rad2_y [km]:         %g", bin_ns.rad2_y);
   CCTK_VInfo (CCTK_THORNSTRING, "rad2_z [km]:         %g", bin_ns.rad2_z);
   CCTK_VInfo (CCTK_THORNSTRING, "rad2_x_opp [km]:     %g", bin_ns.rad2_x_opp);
+  double K = bin_ns.kappa_poly1 * pow(c_light, 6.0) /
+             ( pow(G_grav, 3.0) * M_sun * M_sun *
+               pow(nuc_dens, bin_ns.gamma_poly1-1.) );
+  CCTK_VInfo (CCTK_THORNSTRING, "K [ET unit]:         %.15g", K);
+
   assert (bin_ns.np == npoints);
-  
-  
-  
+
   CCTK_INFO ("Filling in Cactus grid points");
-  
+
 #pragma omp parallel for
   for (int i=0; i<npoints; ++i) {
-    
-    alp[i] = bin_ns.nnn[i];
-    
-    betax[i] = -bin_ns.beta_x[i];
-    betay[i] = -bin_ns.beta_y[i];
-    betaz[i] = -bin_ns.beta_z[i];
-    
-    CCTK_REAL g[3][3];
-    g[0][0] = bin_ns.g_xx[i];
-    g[0][1] = bin_ns.g_xy[i];
-    g[0][2] = bin_ns.g_xz[i];
-    g[1][1] = bin_ns.g_yy[i];
-    g[1][2] = bin_ns.g_yz[i];
-    g[2][2] = bin_ns.g_zz[i];
-    g[1][0] = g[0][1];
-    g[2][0] = g[0][2];
-    g[2][1] = g[1][2];
-    
-    CCTK_REAL ku[3][3];
-    ku[0][0] = bin_ns.k_xx[i];
-    ku[0][1] = bin_ns.k_xy[i];
-    ku[0][2] = bin_ns.k_xz[i];
-    ku[1][1] = bin_ns.k_yy[i];
-    ku[1][2] = bin_ns.k_yz[i];
-    ku[2][2] = bin_ns.k_zz[i];
-    ku[1][0] = ku[0][1];
-    ku[2][0] = ku[0][2];
-    ku[2][1] = ku[1][2];
-    
-    CCTK_REAL k[3][3];
-    for (int a=0; a<3; ++a) {
-      for (int b=0; b<3; ++b) {
-        k[a][b] = 0.0;
-        for (int c=0; c<3; ++c) {
-          for (int d=0; d<3; ++d) {
-            k[a][b] += g[a][c] * g[b][d] * ku[c][d];
-          }
-        }
-      }
+
+    if (CCTK_EQUALS(initial_lapse, "Meudon_Bin_NS")) { 
+      alp[i] = bin_ns.nnn[i];
     }
-    
-    gxx[i] = g[0][0];
-    gxy[i] = g[0][1];
-    gxz[i] = g[0][2];
-    gyy[i] = g[1][1];
-    gyz[i] = g[1][2];
-    gzz[i] = g[2][2];
-    
-    kxx[i] = ku[0][0] * coord_unit;
-    kxy[i] = ku[0][1] * coord_unit;
-    kxz[i] = ku[0][2] * coord_unit;
-    kyy[i] = ku[1][1] * coord_unit;
-    kyz[i] = ku[1][2] * coord_unit;
-    kzz[i] = ku[2][2] * coord_unit;
-    
-    rho[i] = bin_ns.nbar[i] / rho_unit;
-
-    eps[i] = rho[i] * bin_ns.ener_spec[i] / ener_unit;
-    
-    vel[i          ] = bin_ns.u_euler_x[i] / vel_unit;
-    vel[i+  npoints] = bin_ns.u_euler_y[i] / vel_unit;
-    vel[i+2*npoints] = bin_ns.u_euler_z[i] / vel_unit;
-
-    if (rho[i] < 1.e-20) {
-      rho[i          ] = 1.e-20;
-      vel[i          ] = 0.0;
-      vel[i+  npoints] = 0.0;
-      vel[i+2*npoints] = 0.0;
+
+    if (CCTK_EQUALS(initial_shift, "Meudon_Bin_NS")) { 
+      betax[i] = -bin_ns.beta_x[i];
+      betay[i] = -bin_ns.beta_y[i];
+      betaz[i] = -bin_ns.beta_z[i];
+    }
+
+    if (CCTK_EQUALS(initial_data, "Meudon_Bin_NS")) {
+      gxx[i] = bin_ns.g_xx[i];
+      gxy[i] = bin_ns.g_xy[i];
+      gxz[i] = bin_ns.g_xz[i];
+      gyy[i] = bin_ns.g_yy[i];
+      gyz[i] = bin_ns.g_yz[i];
+      gzz[i] = bin_ns.g_zz[i];
+
+      kxx[i] = bin_ns.k_xx[i] * coord_unit;
+      kxy[i] = bin_ns.k_xy[i] * coord_unit;
+      kxz[i] = bin_ns.k_xz[i] * coord_unit;
+      kyy[i] = bin_ns.k_yy[i] * coord_unit;
+      kyz[i] = bin_ns.k_yz[i] * coord_unit;
+      kzz[i] = bin_ns.k_zz[i] * coord_unit;
+    }
+
+    if (CCTK_EQUALS(initial_data, "Meudon_Bin_NS")) {
+      rho[i] = bin_ns.nbar[i] / rho_unit;
+      eps[i] = bin_ns.ener_spec[i];
+      // The following would recompute eps using the polytopic EOS, but
+      // setting eps directly seems to work as well
+      // eps[i] = K * pow(rho[i], bin_ns.gamma_poly1-1.) / (bin_ns.gamma_poly1-1.);
+
+      // TODO: we should really use some EOS calls for the pressure, but for the
+      //       moment this works with polytropes at least
+      press[i] = K * pow(rho[i], bin_ns.gamma_poly1);
+
+      vel[i          ] = bin_ns.u_euler_x[i];
+      vel[i+  npoints] = bin_ns.u_euler_y[i];
+      vel[i+2*npoints] = bin_ns.u_euler_z[i];
+
+      // Especially the velocity is set to strange values outside of the
+      // matter region, so take care of this in the following way
+      if (rho[i] < 1.e-20) {
+        rho[i          ] = 1.e-20;
+        vel[i          ] = 0.0;
+        vel[i+  npoints] = 0.0;
+        vel[i+2*npoints] = 0.0;
+        eps[i          ] = K * pow(rho[i], bin_ns.gamma_poly1-1.) / (bin_ns.gamma_poly1-1.);
+        press[i        ] = K * pow(rho[i], bin_ns.gamma_poly1);
+      }
     }
 
   } // for i
-  
-  CCTK_INFO ("Calculating time derivatives of lapse and shift");
+
   {
     // Angular velocity
     CCTK_REAL const omega = bin_ns.omega * cactusT;
-    
+
     // These initial data assume a helical Killing vector field
-    
-    if (CCTK_EQUALS (initial_dtlapse, "ID_Bin_NS")) {
-      set_dt_from_domega (CCTK_PASS_CTOC, alp, dtalp, omega);
-    } else if (CCTK_EQUALS (initial_dtlapse, "none")) {
-      // do nothing
-    } else {
-      CCTK_WARN (CCTK_WARN_ABORT, "internal error");
+
+    if (CCTK_EQUALS(initial_lapse, "Meudon_Bin_NS")) { 
+      if (CCTK_EQUALS (initial_dtlapse, "Meudon_Bin_NS")) {
+        CCTK_INFO ("Calculating time derivatives of lapse");
+        set_dt_from_domega (CCTK_PASS_CTOC, alp, dtalp, omega);
+      } else if (CCTK_EQUALS (initial_dtlapse, "none")) {
+        // do nothing
+      } else {
+        CCTK_WARN (CCTK_WARN_ABORT, "internal error");
+      }
     }
-    
-    if (CCTK_EQUALS (initial_dtshift, "ID_Bin_NS")) {
-      set_dt_from_domega (CCTK_PASS_CTOC, betax, dtbetax, omega);
-      set_dt_from_domega (CCTK_PASS_CTOC, betay, dtbetay, omega);
-      set_dt_from_domega (CCTK_PASS_CTOC, betaz, dtbetaz, omega);
-    } else if (CCTK_EQUALS (initial_dtshift, "none")) {
-      // do nothing
-    } else {
-      CCTK_WARN (CCTK_WARN_ABORT, "internal error");
+
+    if (CCTK_EQUALS(initial_shift, "Meudon_Bin_NS")) { 
+      if (CCTK_EQUALS (initial_dtshift, "Meudon_Bin_NS")) {
+        CCTK_INFO ("Calculating time derivatives of shift");
+        set_dt_from_domega (CCTK_PASS_CTOC, betax, dtbetax, omega);
+        set_dt_from_domega (CCTK_PASS_CTOC, betay, dtbetay, omega);
+        set_dt_from_domega (CCTK_PASS_CTOC, betaz, dtbetaz, omega);
+      } else if (CCTK_EQUALS (initial_dtshift, "none")) {
+        // do nothing
+      } else {
+        CCTK_WARN (CCTK_WARN_ABORT, "internal error");
+      }
     }
   }
-  
-  
-  
+
   CCTK_INFO ("Done.");
 }
+
diff --git a/src/check_parameters.cc b/src/check_parameters.cc
deleted file mode 100644
index d96379d..0000000
--- a/src/check_parameters.cc
+++ /dev/null
@@ -1,24 +0,0 @@
-#include <cctk.h>
-#include <cctk_Arguments.h>
-#include <cctk_Parameters.h>
-
-
-
-extern "C"
-void ID_Bin_NS_check_parameters (CCTK_ARGUMENTS)
-{
-  DECLARE_CCTK_ARGUMENTS;
-  DECLARE_CCTK_PARAMETERS;
-  
-  if (not CCTK_EQUALS (initial_data,    "ID_Bin_NS") or
-      not CCTK_EQUALS (initial_lapse,   "ID_Bin_NS") or
-      not CCTK_EQUALS (initial_shift,   "ID_Bin_NS") or
-      not (CCTK_EQUALS (initial_dtlapse, "ID_Bin_NS") or
-           CCTK_EQUALS (initial_dtlapse, "none")) or
-      not (CCTK_EQUALS (initial_dtshift, "ID_Bin_NS") or
-           CCTK_EQUALS (initial_dtshift, "none")) or
-      not CCTK_EQUALS (initial_hydro,   "ID_Bin_NS"))
-  {
-    CCTK_PARAMWARN ("The parameters ADMBase::initial_data, ADMBase::initial_lapse, ADMBase::initial_shift, ADMBase::initial_dtlapse, ADMBase::initial_dtshift, and HydroBase::initial_hydro must all be set to the value \"ID_Bin_NS\" or \"none\"");
-  }
-}
diff --git a/src/make.code.defn b/src/make.code.defn
index bc238db..8656099 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -1,7 +1,7 @@
-# Main make.code.defn file for thorn ID_Bin_NS
+# Main make.code.defn file for thorn Meudon_Bin_NS
 
 # Source files in this directory
-SRCS = Bin_NS.cc check_parameters.cc
+SRCS = Bin_NS.cc
 
 # Subdirectories containing source files
 SUBDIRS =