From 131642aa715a3a72c4f367489a699d97968dfe20 Mon Sep 17 00:00:00 2001 From: Yohai Meiron Date: Wed, 25 Mar 2020 21:42:29 -0400 Subject: [PATCH] Added logarithmic external potential --- Makefile | 2 +- config.cpp | 13 +++++-- config.h | 6 ++- phigrape.conf | 14 ++++--- phigrape.cpp | 100 ++++++++++++++------------------------------------ 5 files changed, 51 insertions(+), 84 deletions(-) diff --git a/Makefile b/Makefile index e109393..933616e 100644 --- a/Makefile +++ b/Makefile @@ -1,5 +1,5 @@ CUDAHOME ?= /usr/local/cuda -CPPFLAGS += -DYEBISU -DETICS +CPPFLAGS += -DETICS OPTIMIZATION ?= 3 ETICS_DTSCF ?= 0.015625 diff --git a/config.cpp b/config.cpp index 301e259..25f79cb 100644 --- a/config.cpp +++ b/config.cpp @@ -124,14 +124,16 @@ void Config::error_checking() for (int i=0; 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_units_physical && ((unit_mass == 0) || (unit_length == 0))) + throw std::runtime_error("Physical units for external gravity (ext_units_physical) requires ext_unit_mass and ext_unit_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"); + if (((ext_log_halo_r > 0) && (ext_log_halo_v <= 0)) || ((ext_log_halo_r <= 0) && (ext_log_halo_v > 0))) + throw std::runtime_error("To use external logarithmic halo gravity, please specify positive ext_log_halo_r and ext_log_halo_v"); } Config::Config(std::string file_name) @@ -162,8 +164,8 @@ Config::Config(std::string file_name) 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); + unit_mass = get_parameter(dictionary, "unit_mass", !ext_units_physical); + unit_length = get_parameter(dictionary, "unit_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); @@ -171,6 +173,9 @@ Config::Config(std::string file_name) 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); + ext_log_halo_v = get_parameter(dictionary, "ext_log_halo_v", 0); + ext_log_halo_r = get_parameter(dictionary, "ext_log_halo_r", 0); + error_checking(); } diff --git a/config.h b/config.h index f92d21b..5275bf7 100644 --- a/config.h +++ b/config.h @@ -26,8 +26,8 @@ public: std::vector pn_usage; double pn_c; bool ext_units_physical; - double norm_mass; - double norm_length; + double unit_mass; + double unit_length; double ext_m_bulge; double ext_b_bulge; double ext_m_disk; @@ -35,6 +35,8 @@ public: double ext_b_disk; double ext_m_halo_plummer; double ext_b_halo_plummer; + double ext_log_halo_r; + double ext_log_halo_v; private: using Dictionary = std::unordered_map; diff --git a/phigrape.conf b/phigrape.conf index 2414d1e..c9394b4 100644 --- a/phigrape.conf +++ b/phigrape.conf @@ -34,19 +34,23 @@ dt_min_warning = false #################### ext_units_physical = true -norm_mass = 1e5 -norm_length = 10 +unit_mass = 1e5 +unit_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 = 0 ext_b_bulge = 0 -ext_m_disk = 1e9 -ext_a_disk = 0.8 -ext_b_disk = 0.2 +ext_m_disk = 0 +ext_a_disk = 0 +ext_b_disk = 0 ext_m_halo_plummer = 0 ext_b_halo_plummer = 0 +ext_log_halo_v = 30 +# here the length scale for the logarithmic halo is in kpc, not pc like in the old phigrape. +ext_log_halo_r = 0.04 + #################################### diff --git a/phigrape.cpp b/phigrape.cpp index 99d88e7..04059b8 100644 --- a/phigrape.cpp +++ b/phigrape.cpp @@ -62,7 +62,6 @@ Config *config; //#define EXTPOT_GAL_DEH // Dehnen Galactic - usually physical units -//#define EXTPOT_GAL_LOG // Log Galactic - usually physical units #ifdef ETICS #include "grapite.h" @@ -312,12 +311,6 @@ double m_ext, r_ext, g_ext, tmp_r2, tmp_r3, dum, dum2, dum3, dum_g, tmp_g; #endif -#ifdef EXTPOT_GAL_LOG -double v_halo, r_halo, - v2_halo, r2_halo, r2_r2_halo, - x2_ij, y2_ij, z2_ij; -#endif - #ifdef EXTPOT_BH double r2, rv_ij, m_bh, b_bh, eps_bh; // NOTE different from eps_BH @@ -591,6 +584,24 @@ private: double m, a, b; }; +class Logarithmic_halo : public External_gravity { +public: + Logarithmic_halo(double v_halo, double r_halo) : v2_halo(v_halo*v_halo), r2_halo(r_halo*r_halo) {is_active=(r_halo>0);} +private: + void calc_gravity() + { + auto r2 = x.norm2(); + auto rv_ij = 2.0*(v*x); + auto r2_r2_halo = (r2 + r2_halo); + potential = - 0.5*v2_halo * log(1.0 + r2/r2_halo); + auto tmp = v2_halo/r2_r2_halo; + acceleration = - tmp * x; + tmp /= (r2 + r2_halo); + jerk = tmp * (rv_ij * x - r2_r2_halo * v); + } + double v2_halo, r2_halo; +}; + void calc_self_grav(double t, double eps2, double &g6_calls, int n_loc, int n_act, int ind_act[], double3 x_act_new[], double3 v_act_new[], @@ -688,44 +699,6 @@ for (i=0; iext_units_physical) { - normalization_mass = 1/config->norm_mass; - normalization_length = 1000/config->norm_length; + normalization_mass = 1/config->unit_mass; + normalization_length = 1000/config->unit_length; + normalization_velocity = 1.52484071426404437233e+01*sqrt(config->unit_length/config->unit_mass); } Plummer ext_bulge(config->ext_m_bulge*normalization_mass, config->ext_b_bulge*normalization_length); Miyamoto_Nagai ext_disk(config->ext_m_disk*normalization_mass, config->ext_a_disk*normalization_length, config->ext_b_disk*normalization_length); Plummer ext_halo_plummer(config->ext_m_halo_plummer*normalization_mass, config->ext_b_halo_plummer*normalization_length); + Logarithmic_halo ext_log_halo(config->ext_log_halo_v*normalization_velocity, config->ext_log_halo_r*normalization_length); std::vector external_gravity_components; external_gravity_components.push_back(&ext_bulge); external_gravity_components.push_back(&ext_disk); external_gravity_components.push_back(&ext_halo_plummer); + external_gravity_components.push_back(&ext_log_halo); - if (config->ext_m_bulge > 0) printf("m_bulge = %.4E b_bulge = %.4E\n", config->ext_m_bulge*normalization_mass, config->ext_b_bulge*normalization_length); - if (config->ext_m_disk > 0) printf("m_disk = %.4E a_disk = %.4E b_disk = %.4E\n", config->ext_m_disk*normalization_mass, config->ext_a_disk*normalization_length, config->ext_b_disk*normalization_length); - if (config->ext_m_halo_plummer > 0) printf("m_halo = %.4E b_halo = %.4E\n", config->ext_m_halo_plummer*normalization_mass, config->ext_b_halo_plummer*normalization_length); + if (ext_bulge.is_active) printf("m_bulge = %.4E b_bulge = %.4E\n", config->ext_m_bulge*normalization_mass, config->ext_b_bulge*normalization_length); + if (ext_disk.is_active) printf("m_disk = %.4E a_disk = %.4E b_disk = %.4E\n", config->ext_m_disk*normalization_mass, config->ext_a_disk*normalization_length, config->ext_b_disk*normalization_length); + if (ext_halo_plummer.is_active) printf("m_halo = %.4E b_halo = %.4E\n", config->ext_m_halo_plummer*normalization_mass, config->ext_b_halo_plummer*normalization_length); + if (ext_log_halo.is_active) printf("v_halo = %.6E r_halo = %.6E \n", config->ext_log_halo_v*normalization_velocity, config->ext_log_halo_r*normalization_length); printf("\n"); fflush(stdout);