diff options
author | schnetter <schnetter@9a8dfbc3-8997-403f-b745-98a780c092cc> | 2009-11-24 21:52:24 +0000 |
---|---|---|
committer | schnetter <schnetter@9a8dfbc3-8997-403f-b745-98a780c092cc> | 2009-11-24 21:52:24 +0000 |
commit | 610bb2b5de5eeccc0e78cdcaf99dbe0be92ad80c (patch) | |
tree | 66e85da925b9849cffe5142e99fea843f923f5a9 /src | |
parent | 9e08ae4b7c22f48c50786326beb45020f7aa8452 (diff) |
Add three thorns ID_Bin_BH, ID_Bin_NS, and ID_Mag_NS that import
LORENE initial data.
git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/Meudon_Mag_NS/trunk@2 9a8dfbc3-8997-403f-b745-98a780c092cc
Diffstat (limited to 'src')
-rw-r--r-- | src/Mag_NS.cc | 166 | ||||
-rw-r--r-- | src/check_parameters.cc | 21 | ||||
-rw-r--r-- | src/make.code.defn | 7 |
3 files changed, 194 insertions, 0 deletions
diff --git a/src/Mag_NS.cc b/src/Mag_NS.cc new file mode 100644 index 0000000..af67c4c --- /dev/null +++ b/src/Mag_NS.cc @@ -0,0 +1,166 @@ +#include <cassert> +#include <cstdio> +#include <vector> + +#include <cctk.h> +#include <cctk_Arguments.h> +#include <cctk_Parameters.h> + +#include <mag_ns.h> + +using namespace std; + + + +extern "C" +void ID_Mag_NS_initialise (CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + CCTK_INFO ("Setting up LORENE Mag_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] + CCTK_REAL const mu0 = 4*M_PI * 1.0e-7; // vacuum permeability [N/A^2] + CCTK_REAL const eps0 = 1 / (mu0 * pow(c_light,2)); + + // 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] + + // 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; + CCTK_REAL const cactusB = + sqrt(4*M_PI) / cactusL / sqrt(4*M_PI * eps0 * G_grav / pow(c_light,2)); + + // Other quantities in terms of Cactus units + CCTK_REAL const coord_unit = cactusL / 1.0e+3; // from km + 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_REAL const B_unit = cactusB / 1.0e+9; // from 10^9 T + + + + 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); + + Mag_NS mag_ns (npoints, &xx[0], &yy[0], &zz[0], filename); + + CCTK_VInfo (CCTK_THORNSTRING, "omega [rad/s]: %g", mag_ns.omega); + CCTK_VInfo (CCTK_THORNSTRING, "rho_c [kg/m^3] %g", mag_ns.rho_c); + CCTK_VInfo (CCTK_THORNSTRING, "eps_c [c^2]: %g", mag_ns.eps_c); + CCTK_VInfo (CCTK_THORNSTRING, "mass_b [M_sun]: %g", mag_ns.mass_b); + CCTK_VInfo (CCTK_THORNSTRING, "mass_g [M_sun]: %g", mag_ns.mass_g); + CCTK_VInfo (CCTK_THORNSTRING, "r_eq [km]: %g", mag_ns.r_eq); + CCTK_VInfo (CCTK_THORNSTRING, "r_p [km]: %g", mag_ns.r_p); + CCTK_VInfo (CCTK_THORNSTRING, "L [G M_sum^2/c]: %g", mag_ns.angu_mom); + CCTK_VInfo (CCTK_THORNSTRING, "T/W: %g", mag_ns.T_over_W); + CCTK_VInfo (CCTK_THORNSTRING, "mu [A m^2]: %g", mag_ns.magn_mom); + CCTK_VInfo (CCTK_THORNSTRING, "B_z_pole [10^9 T]: %g", mag_ns.b_z_pole); + CCTK_VInfo (CCTK_THORNSTRING, "B_z_eq [10^9 T]: %g", mag_ns.b_z_eq); + assert (mag_ns.np == npoints); + + + + CCTK_INFO ("Filling in Cactus grid points"); + +#pragma omp parallel for + for (int i=0; i<npoints; ++i) { + + alp[i] = mag_ns.nnn[i]; + + betax[i] = mag_ns.beta_x[i]; + betay[i] = mag_ns.beta_y[i]; + betaz[i] = mag_ns.beta_z[i]; + + CCTK_REAL g[3][3]; + g[0][0] = mag_ns.g_xx[i]; + g[0][1] = mag_ns.g_xy[i]; + g[0][2] = mag_ns.g_xz[i]; + g[1][1] = mag_ns.g_yy[i]; + g[1][2] = mag_ns.g_yz[i]; + g[2][2] = mag_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] = mag_ns.k_xx[i]; + ku[0][1] = mag_ns.k_xy[i]; + ku[0][2] = mag_ns.k_xz[i]; + ku[1][1] = mag_ns.k_yy[i]; + ku[1][2] = mag_ns.k_yz[i]; + ku[2][2] = mag_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]; + } + } + } + } + + 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] = k[0][0]; + kxy[i] = k[0][1]; + kxz[i] = k[0][2]; + kyy[i] = k[1][1]; + kyz[i] = k[1][2]; + kzz[i] = k[2][2]; + + rho[i] = mag_ns.nbar[i] / rho_unit; + + eps[i] = rho[i] * mag_ns.ener_spec[i] / ener_unit; + + vel[i ] = mag_ns.u_euler_x[i] / vel_unit; + vel[i+ npoints] = mag_ns.u_euler_y[i] / vel_unit; + vel[i+2*npoints] = mag_ns.u_euler_z[i] / vel_unit; + + Bmag[i ] = mag_ns.bb_x[i] / B_unit; + Bmag[i+ npoints] = mag_ns.bb_y[i] / B_unit; + Bmag[i+2*npoints] = mag_ns.bb_z[i] / B_unit; + + } // for i + + + + CCTK_INFO ("Done."); +} diff --git a/src/check_parameters.cc b/src/check_parameters.cc new file mode 100644 index 0000000..fac34a1 --- /dev/null +++ b/src/check_parameters.cc @@ -0,0 +1,21 @@ +#include <cctk.h> +#include <cctk_Arguments.h> +#include <cctk_Parameters.h> + + + +extern "C" +void ID_Mag_NS_check_parameters (CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + if (not CCTK_EQUALS (initial_data , "ID_Mag_NS") or + not CCTK_EQUALS (initial_lapse, "ID_Mag_NS") or + not CCTK_EQUALS (initial_shift, "ID_Mag_NS") or + not CCTK_EQUALS (initial_hydro, "ID_Mag_NS") or + not CCTK_EQUALS (initial_Bmag, "ID_Mag_NS")) + { + CCTK_PARAMWARN ("The parameters ADMBase::initial_data, ADMBase::initial_lapse, ADMBase::initial_shift, HydroBase::initial_hydro, and HydroBase::initial_Bmag must all be set to the value \"ID_Mag_NS\""); + } +} diff --git a/src/make.code.defn b/src/make.code.defn new file mode 100644 index 0000000..329f9f1 --- /dev/null +++ b/src/make.code.defn @@ -0,0 +1,7 @@ +# Main make.code.defn file for thorn ID_Mag_NS + +# Source files in this directory +SRCS = Mag_NS.cc check_parameters.cc + +# Subdirectories containing source files +SUBDIRS = |