User: knarf Date: 2010/08/04 10:34 PM
Removed: /trunk/src/ check_parameters.cc
Modified: /trunk/ README, configuration.ccl, interface.ccl, param.ccl, schedule.ccl /trunk/src/ Bin_NS.cc, make.code.defn
Log: 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
File Changes:
Directory: /trunk/src/ ======================
File [modified]: Bin_NS.cc Delta lines: +109 -122 =================================================================== --- trunk/src/Bin_NS.cc 2010-07-30 03:19:48 UTC (rev 8) +++ trunk/src/Bin_NS.cc 2010-08-05 03:34:36 UTC (rev 9) @@ -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 @@ } }
- - 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 @@ 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 (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 (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_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."); } +
File [removed]: check_parameters.cc Delta lines: +0 -24 =================================================================== --- trunk/src/check_parameters.cc 2010-07-30 03:19:48 UTC (rev 8) +++ trunk/src/check_parameters.cc 2010-08-05 03:34:36 UTC (rev 9) @@ -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""); - } -}
File [modified]: make.code.defn Delta lines: +2 -2 =================================================================== --- trunk/src/make.code.defn 2010-07-30 03:19:48 UTC (rev 8) +++ trunk/src/make.code.defn 2010-08-05 03:34:36 UTC (rev 9) @@ -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 =
Directory: /trunk/ ==================
File [modified]: README Delta lines: +6 -5 =================================================================== --- trunk/README 2010-07-30 03:19:48 UTC (rev 8) +++ trunk/README 2010-08-05 03:34:36 UTC (rev 9) @@ -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 @@
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.
File [modified]: configuration.ccl Delta lines: +1 -1 =================================================================== --- trunk/configuration.ccl 2010-07-30 03:19:48 UTC (rev 8) +++ trunk/configuration.ccl 2010-08-05 03:34:36 UTC (rev 9) @@ -1,3 +1,3 @@ -# Configuration definition for thorn ID_Bin_NS +# Configuration definition for thorn Meudon_Bin_NS
REQUIRES LORENE
File [modified]: interface.ccl Delta lines: +2 -2 =================================================================== --- trunk/interface.ccl 2010-07-30 03:19:48 UTC (rev 8) +++ trunk/interface.ccl 2010-08-05 03:34:36 UTC (rev 9) @@ -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
File [modified]: param.ccl Delta lines: +7 -7 =================================================================== --- trunk/param.ccl 2010-07-30 03:19:48 UTC (rev 8) +++ trunk/param.ccl 2010-08-05 03:34:36 UTC (rev 9) @@ -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 @@
EXTENDS KEYWORD initial_hydro { - "ID_Bin_NS" :: "" + "Meudon_Bin_NS" :: "" }
File [modified]: schedule.ccl Delta lines: +8 -13 =================================================================== --- trunk/schedule.ccl 2010-07-30 03:19:48 UTC (rev 8) +++ trunk/schedule.ccl 2010-08-05 03:34:36 UTC (rev 9) @@ -1,19 +1,14 @@ -# 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 + SCHEDULE Meudon_Bin_NS_initialise IN HydroBase_Initial { LANG: C - } "Check parameters" - - SCHEDULE ID_Bin_NS_initialise IN HydroBase_Initial - { - LANG: C } "Set up binary neutron star initial data" }