diff --git a/config.cpp b/config.cpp index fb121be..aac05c2 100644 --- a/config.cpp +++ b/config.cpp @@ -110,6 +110,8 @@ void Config::error_checking() throw std::runtime_error("The number of live black holes (live_smbh_count) can be 0, 1, or 2"); if (binary_smbh_pn && (live_smbh_count!=2)) throw std::runtime_error("Post-Newtonian gravity (binary_smbh_pn=true) requires live_smbh_count=2"); + if (binary_smbh_pn && (pn_c <= 0)) + throw std::runtime_error("Post-Newtonian gravity (binary_smbh_pn=true) requires pn_c > 0"); if (live_smbh_custom_eps == eps) live_smbh_custom_eps = -1; if (live_smbh_output && (live_smbh_count == 0)) throw std::runtime_error("Black hole output (live_smbh_output=true) requires at least one live black hole (live_smbh_count)"); @@ -122,6 +124,14 @@ void Config::error_checking() for (int i; i<7; i++) if (!((pn_usage[i] == 0) || (pn_usage[i] == 1))) throw std::runtime_error("PN usage array (pn_usage) must have ones and zeros only"); + if (ext_units_physical && ((norm_mass == 0) || (norm_length == 0))) + throw std::runtime_error("Physical units for external gravity (ext_units_physical) requires ext_norm_mass and ext_norm_length to be positive numbers"); + if ((ext_m_bulge > 0) && (ext_b_bulge < 0)) + throw std::runtime_error("To use external bulge gravity, please specify positive ext_m_bulge and ext_b_bulge"); + if ((ext_m_halo_plummer > 0) && (ext_b_halo_plummer < 0)) + throw std::runtime_error("To use external Plummer halo gravity, please specify positive ext_m_halo_plummer and ext_b_halo_plummer"); + if ((ext_m_disk > 0) && ((ext_a_disk < 0) || (ext_b_disk < 0))) + throw std::runtime_error("To use external disk gravity, please specify positive ext_m_disk, ext_a_disk and ext_b_disk"); } Config::Config(std::string file_name) @@ -137,16 +147,30 @@ Config::Config(std::string file_name) eta = get_parameter(dictionary, "eta"); input_file_name = get_parameter(dictionary, "input_file_name", "data.con"); dt_min_warning = get_parameter(dictionary, "dt_min_warning", false); + live_smbh_count = get_parameter(dictionary, "live_smbh_count", 0); live_smbh_custom_eps = get_parameter(dictionary, "live_smbh_custom_eps", -1); live_smbh_output = get_parameter(dictionary, "live_smbh_output", false); live_smbh_neighbor_output = get_parameter(dictionary, "live_smbh_neighbor_output", false); live_smbh_neighbor_number = get_parameter(dictionary, "live_smbh_neighbor_number", 10); - binary_smbh_pn = get_parameter(dictionary, "binary_smbh_pn", false); + binary_smbh_influence_sphere_output = get_parameter(dictionary, "binary_smbh_influence_sphere_output", false); binary_smbh_influence_radius_factor = get_parameter(dictionary, "binary_smbh_influence_radius_factor", 10.); + + binary_smbh_pn = get_parameter(dictionary, "binary_smbh_pn", false); pn_usage = get_parameter>(dictionary, "pn_usage", std::vector({1,1,1,1,1,1,1})); - pn_c = get_parameter(dictionary, "pn_c", 500); + pn_c = get_parameter(dictionary, "pn_c", 0); + + ext_units_physical = get_parameter(dictionary, "ext_units_physical", false); + norm_mass = get_parameter(dictionary, "norm_mass", !ext_units_physical); + norm_length = get_parameter(dictionary, "norm_length", !ext_units_physical); + ext_m_bulge = get_parameter(dictionary, "ext_m_bulge", 0); + ext_b_bulge = get_parameter(dictionary, "ext_b_bulge", -1); + ext_m_disk = get_parameter(dictionary, "ext_m_disk", 0); + ext_a_disk = get_parameter(dictionary, "ext_a_disk", -1); + ext_b_disk = get_parameter(dictionary, "ext_b_disk", -1); + ext_m_halo_plummer = get_parameter(dictionary, "ext_m_halo_plummer", 0); + ext_b_halo_plummer = get_parameter(dictionary, "ext_b_halo_plummer", -1); error_checking(); } diff --git a/config.h b/config.h index 6ce5d60..f92d21b 100644 --- a/config.h +++ b/config.h @@ -25,6 +25,17 @@ public: double binary_smbh_influence_radius_factor; std::vector pn_usage; double pn_c; + bool ext_units_physical; + double norm_mass; + double norm_length; + double ext_m_bulge; + double ext_b_bulge; + double ext_m_disk; + double ext_a_disk; + double ext_b_disk; + double ext_m_halo_plummer; + double ext_b_halo_plummer; + private: using Dictionary = std::unordered_map; Dictionary read_config_file(const std::string file_name); diff --git a/phigrape.conf b/phigrape.conf index 245a644..337b154 100644 --- a/phigrape.conf +++ b/phigrape.conf @@ -29,67 +29,87 @@ output_format = HDF5 dt_min_warning = false ############################################# +#################### +# EXTERNAL GRAVITY # +#################### + +ext_units_physical = true +norm_mass = 1e5 +norm_length = 10 +# TODO add the option to normalize using other units + +# Notice that if physical units are used, the "a" and "b" parameters are in kiloparsec, not parsec! + +ext_m_bulge = 1e5 +ext_b_bulge = 20 +ext_m_disk = 1e6 +ext_a_disk = 80 +ext_b_disk = 70 +ext_m_halo_plummer = 1e6 +ext_b_halo_plummer = 200. + + +#################################### +## LIVE SUPERMASSIVE BLACK HOLE(S) # +#################################### +# +## There is special treatment for particles representing supermassive black holes (SMBHs): they are integrated at every time step, they can have custom softening in SMBH-SMBH interactions, and post Newtonian terms can be added to the gravity. +# +## The number of SMBH particles. Can be 0 (no SMBH), 1, or 2. [default: 0] +#live_smbh_count = 0 +# +## Custom softening length for SMBH-SMBH interactions (can also be zero). If non-negative, the custom softening is applied. [default: -1] +#live_smbh_custom_eps = 0 +# +## Output additional diagnostics about live SMBHs. [default: false] +##TODO# dt_bh +#live_smbh_output = true +# +## Output additional diagnostics about the SMBH's (or SMBHs') nearest neighbours (number could be set as shown below). [default: false] +#live_smbh_neighbor_output = true +# +## Number of nearest neighbours to the SMBH (or SMBHs) to include in output. [default: 10] +#live_smbh_neighbor_number = 10 +# ################################### -# LIVE SUPERMASSIVE BLACK HOLE(S) # +## BINARY SUPERMASSIVE BLACK HOLE # ################################### - -# There is special treatment for particles representing supermassive black holes (SMBHs): they are integrated at every time step, they can have custom softening in SMBH-SMBH interactions, and post Newtonian terms can be added to the gravity. - -# The number of SMBH particles. Can be 0 (no SMBH), 1, or 2. [default: 0] -live_smbh_count = 2 - -# Custom softening length for SMBH-SMBH interactions (can also be zero). If non-negative, the custom softening is applied. [default: -1] -live_smbh_custom_eps = 0 - -# Output additional diagnostics about live SMBHs. [default: false] -#TODO# dt_bh -live_smbh_output = true - -# Output additional diagnostics about the SMBH's (or SMBHs') nearest neighbours (number could be set as shown below). [default: false] -live_smbh_neighbor_output = true - -# Number of nearest neighbours to the SMBH (or SMBHs) to include in output. [default: 10] -live_smbh_neighbor_number = 10 - -################################## -# BINARY SUPERMASSIVE BLACK HOLE # -################################## - -# The following parameters can be set when `live_smbh_count` is 2. - -# Output additional diagnostics about the SMBH's sphere of influence (size could be set as shown below). [default: false] -binary_smbh_influence_sphere_output = true - -# The influence sphere is centred at the binary SMBH's centre of mass, and its radius is the semi-major axis of the binary times the factor below. [default: 10.0] -binary_smbh_influence_radius_factor = 3.162277660168379497918067e+03 - -# Add post Newtonian terms to SMBH-SMBH gravity. [default: false] -binary_smbh_pn = true - -# A mask array (zeros and ones) determining whether or not to use specific post-Newtonian terms. -# The elements represent {0, 1, 2, 2.5, 3, 3.5, spin} -pn_usage = {1, 1, 1, 1, 0, 0, 0} - -# The speed of light in N-body units [default: 500] -pn_c = 477.12 - -#################################### -# Negative powers of two # -#################################### -# -1 1/2 0.5 # -# -2 1/4 0.25 # -# -3 1/8 0.125 # -# -4 1/16 0.0625 # -# -5 1/32 0.03125 # -# -6 1/64 0.015625 # -# -7 1/128 0.0078125 # -# -8 1/256 0.00390625 # -# -9 1/512 0.001953125 # -# -10 1/1024 0.0009765625 # -# -11 1/2048 0.00048828125 # -# -12 1/4096 0.000244140625 # -# -13 1/8192 0.0001220703125 # -# -14 1/16384 0.00006103515625 # -# -15 1/32768 0.000030517578125 # -# -16 1/65536 0.0000152587890625 # -#################################### +# +## The following parameters can be set when `live_smbh_count` is 2. +# +## Output additional diagnostics about the SMBH's sphere of influence (size could be set as shown below). [default: false] +#binary_smbh_influence_sphere_output = true +# +## The influence sphere is centred at the binary SMBH's centre of mass, and its radius is the semi-major axis of the binary times the factor below. [default: 10.0] +#binary_smbh_influence_radius_factor = 3.162277660168379497918067e+03 +# +## Add post Newtonian terms to SMBH-SMBH gravity. [default: false] +#binary_smbh_pn = true +# +## A mask array (zeros and ones) determining whether or not to use specific post-Newtonian terms. +## The elements represent {0, 1, 2, 2.5, 3, 3.5, spin} +#pn_usage = {1, 1, 1, 1, 0, 0, 0} +# +## The speed of light in N-body units [default: 500] +#pn_c = 477.12 +# +##################################### +## Negative powers of two # +##################################### +## -1 1/2 0.5 # +## -2 1/4 0.25 # +## -3 1/8 0.125 # +## -4 1/16 0.0625 # +## -5 1/32 0.03125 # +## -6 1/64 0.015625 # +## -7 1/128 0.0078125 # +## -8 1/256 0.00390625 # +## -9 1/512 0.001953125 # +## -10 1/1024 0.0009765625 # +## -11 1/2048 0.00048828125 # +## -12 1/4096 0.000244140625 # +## -13 1/8192 0.0001220703125 # +## -14 1/16384 0.00006103515625 # +## -15 1/32768 0.000030517578125 # +## -16 1/65536 0.0000152587890625 # +##################################### diff --git a/phigrape.cpp b/phigrape.cpp index dbf9514..f46dc36 100644 --- a/phigrape.cpp +++ b/phigrape.cpp @@ -56,12 +56,11 @@ Last redaction : 2019.04.16 12:55 #include "config.h" Config *config; -//#define NORM // Physical normalization +#define NORM // Physical normalization -//#define EXTPOT // external potential (BH? or galactic?) +#define EXTPOT // external potential (BH? or galactic?) //#define EXTPOT_BH // BH - usually NB units -//#define EXTPOT_GAL // Galactic B+D+H PK - usually physical units //#define EXTPOT_GAL_DEH // Dehnen Galactic - usually physical units //#define EXTPOT_GAL_LOG // Log Galactic - usually physical units @@ -112,6 +111,7 @@ Config *config; #include #include +#include #define G6_NPIPE 2048 #include "grape6.h" @@ -217,7 +217,11 @@ struct double3 { data[2] /= c; return *this; } - double norm() + double norm2() const + { + return data[0]*data[0]+data[1]*data[1]+data[2]*data[2]; + } + double norm() const { return sqrt(data[0]*data[0]+data[1]*data[1]+data[2]*data[2]); } @@ -234,6 +238,11 @@ double3 operator*(const double3& a, const double& c) return double3(a.data[0]*c, a.data[1]*c, a.data[2]*c); } +double operator*(const double3& a, const double3& b) +{ + return a.data[0]*b.data[0]+a.data[1]*b.data[1]+a.data[2]*b.data[2]; +} + double3 operator/(const double3& a, const double& c) { return double3(a.data[0]/c, a.data[1]/c, a.data[2]/c); @@ -292,9 +301,7 @@ double3 a2by18, a1by6, aby2; /* normalization... */ -#ifdef NORM double m_norm, r_norm, v_norm, t_norm; -#endif double eps_BH=0.0; @@ -302,14 +309,6 @@ double eps_BH=0.0; #ifdef EXTPOT -#ifdef EXTPOT_GAL -double m_bulge, a_bulge, b_bulge, - m_disk, a_disk, b_disk, - m_halo, a_halo, b_halo, // NOTE there are other "halo" variables in EXTPOT_GAL_LOG - x2_ij, y2_ij, z2_ij, // NOTE this variable exist also in EXTPOT_GAL_LOG - r_tmp, r2_tmp, z_tmp, z2_tmp; -#endif - #ifdef EXTPOT_GAL_DEH double m_ext, r_ext, g_ext, tmp_r2, tmp_r3, dum, dum2, dum3, dum_g, tmp_g; @@ -624,254 +623,72 @@ void write_bh_nb_data(double time_cur, int N, double m[], double3 x[], double3 v fclose(out); } -void calc_ext_grav_zero() -{ - -#ifdef EXTPOT - -/* Define the external potential for all particles on all nodes */ - -ni = n_act; - -for (i=0; i 0) + external_gravity.calc_plummer(ni, x_act_new,v_act_new, pot_act_ext, a_act_new, adot_act_new, 1); // Bulge + if (external_gravity.m_disk > 0) + external_gravity.calc_disk(ni, x_act_new,v_act_new, pot_act_ext, a_act_new, adot_act_new); // Disk + if (external_gravity.m_halo > 0) + external_gravity.calc_plummer(ni, x_act_new,v_act_new, pot_act_ext, a_act_new, adot_act_new, 2); // Halo + + + +for (int i=0; iext_m_bulge/config->norm_mass; + external_gravity.b_bulge = config->ext_b_bulge/config->norm_length*1000; + external_gravity.m_disk = config->ext_m_disk/config->norm_mass; + external_gravity.a_disk = config->ext_a_disk/config->norm_length*1000; + external_gravity.b_disk = config->ext_b_disk/config->norm_length*1000; + external_gravity.m_halo = config->ext_m_halo_plummer/config->norm_mass; + external_gravity.b_halo = config->ext_b_halo_plummer/config->norm_length*1000; + + if (external_gravity.m_bulge > 0) + printf("m_bulge = %.4E b_bulge = %.4E\n", external_gravity.m_bulge, external_gravity.b_bulge); + if (external_gravity.m_disk > 0) + printf("m_disk = %.4E a_disk = %.4E b_disk = %.4E\n", external_gravity.m_disk, external_gravity.a_disk, external_gravity.b_disk); + if (external_gravity.m_halo > 0) + printf("m_halo = %.4E b_halo = %.4E\n", external_gravity.m_halo, external_gravity.b_halo); + + eta_s = eta/ETA_S_CORR; eta_bh = eta/ETA_BH_CORR; @@ -1636,7 +1323,7 @@ int main(int argc, char *argv[]) /* init the local GRAPE's */ #ifdef MPI_OVERRIDE - numGPU = 1; + numGPU = 1; // TODO get this from config file clusterid = myRank % numGPU; #else MPI_Comm shmcomm; @@ -1816,7 +1503,8 @@ int main(int argc, char *argv[]) } #ifdef EXTPOT - calc_ext_grav_zero(); +// /*/*/*/*/*/*calc_ext_grav_zero();*/*/*/*/*/*/ + calc_ext_grav(n_act, x_act, v_act, a, adot); #endif /* Wait to all processors to finish his works... */ @@ -2223,7 +1911,7 @@ int main(int argc, char *argv[]) } #ifdef EXTPOT - calc_ext_grav(); + calc_ext_grav(n_act, x_act_new, v_act_new, a_act_new, adot_act_new); #endif /* correct the active particles positions etc... on all the nodes */