From 9296db060960cf6e6aa48a182efe492d9e409823 Mon Sep 17 00:00:00 2001 From: Yohai Meiron Date: Tue, 17 Mar 2020 12:24:53 -0400 Subject: [PATCH] Now config file handles most parameters (but not external potential) --- Makefile | 2 +- config.cpp | 126 +++++++++++++++++++++++++++++++++++++++++ phigrape.h => config.h | 14 ++++- io.cpp | 123 ---------------------------------------- phigrape.conf | 2 +- phigrape.cpp | 86 +++++++++++++++------------- 6 files changed, 188 insertions(+), 165 deletions(-) create mode 100644 config.cpp rename phigrape.h => config.h (51%) delete mode 100644 io.cpp diff --git a/Makefile b/Makefile index 1739015..e109393 100644 --- a/Makefile +++ b/Makefile @@ -21,7 +21,7 @@ MPICXX ?= mpic++ EXECUTABLE ?= phigrape default: - $(MPICXX) $(CPPFLAGS) $(CXXFLAGS) -DETICS_DTSCF=$(ETICS_DTSCF) $(INC) phigrape.cpp -o $(EXECUTABLE) $(LIB) + $(MPICXX) $(CPPFLAGS) $(CXXFLAGS) -DETICS_DTSCF=$(ETICS_DTSCF) $(INC) config.cpp phigrape.cpp -o $(EXECUTABLE) $(LIB) yebisu: CPPFLAGS := $(filter-out -DETICS, $(CPPFLAGS)) yebisu: default diff --git a/config.cpp b/config.cpp new file mode 100644 index 0000000..c4c229e --- /dev/null +++ b/config.cpp @@ -0,0 +1,126 @@ +#include "config.h" +#include +#include +#include + +using Dictionary = std::unordered_map; + +std::string Config::strip(const std::string str) +{ + std::string str_new = str; + auto pos = str_new.find_first_not_of(" \t"); + if (pos != std::string::npos) str_new = str_new.substr(pos, str_new.size()); + pos = str_new.find_last_not_of(" \t"); + if (pos != std::string::npos) str_new = str_new.substr(0, pos+1); + return str_new; +} + +Dictionary Config::read_config_file(const std::string file_name) +{ + std::unordered_map dictionary; + std::ifstream file(file_name); + if (!file.good()) throw std::runtime_error("File not found."); + std::string str; + int line_number = 0; + while (std::getline(file, str)) { + line_number++; + auto pos = str.find('#'); + if (pos != std::string::npos) str = str.substr(0, pos); + str = strip(str); + if (str.size() == 0) continue; + pos = str.find_first_of("="); + if (pos == std::string::npos) throw std::runtime_error("Error: expected a key-value pair in line " + std::to_string(line_number) + " of file " + file_name); + std::string key = strip(str.substr(0, pos)); + pos = str.find_first_not_of(" \t", pos+1); + std::string val = strip(str.substr(pos, str.size())); + dictionary[key] = val; + } + return dictionary; +} + +template<> std::string Config::string_cast(const std::string str) +{ + return str; +} + +template<> double Config::string_cast(const std::string str) +{ + size_t idx; + auto value = std::stod(str, &idx); + if (idx == str.size()) return value; + else throw std::runtime_error("Cannot convert \"" + str + "\" into a double"); +} + +template<> int Config::string_cast(const std::string str) +{ + size_t idx; + auto value = std::stoi(str, &idx); + if (idx == str.size()) return value; + else throw std::runtime_error("Cannot convert \"" + str + "\" into an int"); +} + +template<> bool Config::string_cast(const std::string str) +{ + if ((str=="true") || (str=="True") || (str=="yes") || (str=="Yes") || (str=="1")) return true; + else if ((str=="false") || (str=="False") || (str=="no") || (str=="No") || (str=="0")) return false; + throw std::runtime_error("Cannot convert \"" + str + "\" into a bool"); +} + + +// For mandatory parameters +template +T Config::get_parameter(Dictionary dictionary, std::string name) +{ + auto item = dictionary.find(name); + if (item==dictionary.end()) throw std::runtime_error("Mandatory parameter " + name + " must be defined"); + else return string_cast((*item).second); +} + +// For optional parameters +template +T Config::get_parameter(Dictionary dictionary, std::string name, T default_value) +{ + auto item = dictionary.find(name); + if (item==dictionary.end()) return default_value; + else return string_cast((*item).second); +} + +void Config::error_checking() +{ + if ((live_smbh_count < 0) || (live_smbh_count > 2)) + 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 (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)"); + if (live_smbh_neighbor_output && (live_smbh_count == 0)) + throw std::runtime_error("Black hole neighbour output (live_smbh_neighbor_output=true) requires at least one live black hole (live_smbh_count)"); + if (binary_smbh_influence_sphere_output && (live_smbh_count != 2)) + throw std::runtime_error("Binary black hole influence sphere output (binary_smbh_influence_sphere_output=true) requires exactly two live black holes (live_smbh_count=2)"); +} + +Config::Config(std::string file_name) +{ + auto dictionary = read_config_file(file_name); + + // TODO check if dt_disk and dt_contr are powers of two + eps = get_parameter(dictionary, "eps"); + t_end = get_parameter(dictionary, "t_end"); + dt_disk = get_parameter(dictionary, "dt_disk"); + dt_contr = get_parameter(dictionary, "dt_contr"); + dt_bh = get_parameter(dictionary, "dt_bh", dt_contr); + 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.); + + error_checking(); +} diff --git a/phigrape.h b/config.h similarity index 51% rename from phigrape.h rename to config.h index 2c68183..5c06a7d 100644 --- a/phigrape.h +++ b/config.h @@ -1,7 +1,11 @@ #pragma once #include +#include + +class Config { +public: + Config(std::string file_name); -struct Parameters { double eps; double t_end; double dt_disk; @@ -18,4 +22,12 @@ struct Parameters { bool binary_smbh_pn; bool binary_smbh_influence_sphere_output; double binary_smbh_influence_radius_factor; +private: + using Dictionary = std::unordered_map; + Dictionary read_config_file(const std::string file_name); + std::string strip(const std::string str); + template T string_cast(const std::string str); + template T get_parameter(Dictionary dictionary, std::string name); + template T get_parameter(Dictionary dictionary, std::string name, T default_value); + void error_checking(); }; diff --git a/io.cpp b/io.cpp deleted file mode 100644 index 29967ec..0000000 --- a/io.cpp +++ /dev/null @@ -1,123 +0,0 @@ -#include "phigrape.h" -#include -#include -#include -#include -#include - -using Dictionary = std::unordered_map; - -std::string strip(const std::string str) -{ - std::string str_new = str; - auto pos = str_new.find_first_not_of(" \t"); - if (pos != std::string::npos) str_new = str_new.substr(pos, str_new.size()); - pos = str_new.find_last_not_of(" \t"); - if (pos != std::string::npos) str_new = str_new.substr(0, pos+1); - return str_new; -} - -Dictionary read_config_file(const std::string file_name) -{ - std::unordered_map dictionary; - std::ifstream file(file_name); - if (!file.good()) throw std::runtime_error("File not found."); - std::string str; - int line_number = 0; - while (std::getline(file, str)) { - line_number++; - auto pos = str.find('#'); - if (pos != std::string::npos) str = str.substr(0, pos); - str = strip(str); - if (str.size() == 0) continue; - pos = str.find_first_of("="); - if (pos == std::string::npos) throw std::runtime_error("Error: expected a key-value pair in line " + std::to_string(line_number) + " of file " + file_name); - std::string key = strip(str.substr(0, pos)); - pos = str.find_first_not_of(" \t", pos+1); - std::string val = strip(str.substr(pos, str.size())); - dictionary[key] = val; - } - return dictionary; -} - -template T string_cast(const std::string str); - -template<> std::string string_cast(const std::string str) -{ - return str; -} - -template<> double string_cast(const std::string str) -{ - size_t idx; - auto value = std::stod(str, &idx); - if (idx == str.size()) return value; - else throw std::runtime_error("Cannot convert \"" + str + "\" into a double"); -} - -template<> int string_cast(const std::string str) -{ - size_t idx; - auto value = std::stoi(str, &idx); - if (idx == str.size()) return value; - else throw std::runtime_error("Cannot convert \"" + str + "\" into an int"); -} - -template<> bool string_cast(const std::string str) -{ - if ((str=="true") || (str=="True") || (str=="yes") || (str=="Yes") || (str=="1")) return true; - else if ((str=="false") || (str=="False") || (str=="no") || (str=="No") || (str=="0")) return false; - throw std::runtime_error("Cannot convert \"" + str + "\" into a bool"); -} - - -// For mandatory parameters -template -T get_parameter(Dictionary dictionary, std::string name) -{ - auto item = dictionary.find(name); - if (item==dictionary.end()) throw std::runtime_error("Mandatory parameter " + name + " must be defined"); - else return string_cast((*item).second); -} - -// For optional parameters -template -T get_parameter(Dictionary dictionary, std::string name, T default_value) -{ - auto item = dictionary.find(name); - if (item==dictionary.end()) return default_value; - else return string_cast((*item).second); -} - -int main() -{ - auto dictionary = read_config_file("phigrape.conf"); - - for (auto key_value : dictionary) { - auto key = key_value.first; - auto value = key_value.second; - printf("dictionary[\"%s\"] = \"%s\"\n", key.c_str(), value.c_str()); - } - - - Parameters parameters; - - // TODO check if dt_disk and dt_contr are powers of two - parameters.dt_bh = get_parameter(dictionary, "eps"); - parameters.t_end = get_parameter(dictionary, "t_end"); - parameters.dt_disk = get_parameter(dictionary, "dt_disk"); - parameters.dt_contr = get_parameter(dictionary, "dt_contr"); - parameters.dt_bh = get_parameter(dictionary, "dt_bh", parameters.dt_contr); - parameters.eta = get_parameter(dictionary, "eta"); - parameters.input_file_name = get_parameter(dictionary, "dt_bh", "data.con"); - parameters.dt_min_warning = get_parameter(dictionary, "dt_min_warning", false); - parameters.live_smbh_count = get_parameter(dictionary, "live_smbh_count", 0); - parameters.live_smbh_custom_eps = get_parameter(dictionary, "live_smbh_custom_eps", -1); - parameters.live_smbh_output = get_parameter(dictionary, "live_smbh_output", false); - parameters.live_smbh_neighbor_output = get_parameter(dictionary, "live_smbh_neighbor_output", false); - parameters.live_smbh_neighbor_number = get_parameter(dictionary, "live_smbh_neighbor_number", 10); - parameters.binary_smbh_pn = get_parameter(dictionary, "binary_smbh_pn", false); - parameters.binary_smbh_influence_sphere_output = get_parameter(dictionary, "binary_smbh_influence_sphere_output", false); - parameters.binary_smbh_influence_radius_factor = get_parameter(dictionary, "binary_smbh_influence_radius_factor", 10.); - -} diff --git a/phigrape.conf b/phigrape.conf index dd01e95..732518f 100644 --- a/phigrape.conf +++ b/phigrape.conf @@ -6,7 +6,7 @@ eps = 1E-4 # End time of the calculation -t_end = 4.0 +t_end = 0.25 # Interval of snapshot files output (xxxxxx.dat) dt_disk = 1.0 diff --git a/phigrape.cpp b/phigrape.cpp index 962d9b3..b194e0b 100644 --- a/phigrape.cpp +++ b/phigrape.cpp @@ -53,8 +53,8 @@ Coded by : Peter Berczik Version number : 19.04 Last redaction : 2019.04.16 12:55 *****************************************************************************/ -#include "phigrape.h" -Parameters parameters; +#include "config.h" +Config *config; // struct Parameters { // double eps = 1E-4; @@ -665,7 +665,7 @@ fclose(out); void write_bh_data() { - if (parameters.live_smbh_count == 2) { + if (config->live_smbh_count == 2) { out = fopen("bh.dat", "a"); @@ -693,12 +693,12 @@ void write_bh_data() tmp_a = sqrt( SQR(a[i][0]) + SQR(a[i][1]) + SQR(a[i][2]) ); tmp_adot = sqrt( SQR(adot[i][0]) + SQR(adot[i][1]) + SQR(adot[i][2]) ); - if (parameters.live_smbh_custom_eps >= 0) { + if (config->live_smbh_custom_eps >= 0) { tmp_a_bh = sqrt( SQR(a_bh[0]) + SQR(a_bh[1]) + SQR(a_bh[2]) ); tmp_adot_bh = sqrt( SQR(adot_bh[0]) + SQR(adot_bh[1]) + SQR(adot_bh[2]) ); - if (parameters.binary_smbh_spin) { + if (config->binary_smbh_pn) { tmp_a_bh_pn0 = sqrt( SQR(a_pn[0][0]) + SQR(a_pn[0][1]) + SQR(a_pn[0][2]) ); tmp_a_bh_pn1 = sqrt( SQR(a_pn[1][0]) + SQR(a_pn[1][1]) + SQR(a_pn[1][2]) ); tmp_a_bh_pn2 = sqrt( SQR(a_pn[2][0]) + SQR(a_pn[2][1]) + SQR(a_pn[2][2]) ); @@ -766,7 +766,7 @@ void write_bh_data() } fprintf(out,"\n"); fclose(out); - } else if (parameters.live_smbh_count == 1) { + } else if (config->live_smbh_count == 1) { out = fopen("bh.dat", "a"); tmp_r = sqrt( SQR(x[0][0]) + SQR(x[0][1]) + SQR(x[0][2]) ); tmp_v = sqrt( SQR(v[0][0]) + SQR(v[0][1]) + SQR(v[0][2]) ); @@ -787,7 +787,7 @@ void write_bh_data() void write_bh_nb_data() { - int i_bh, nb = parameters.live_smbh_neighbor_number; + int i_bh, nb = config->live_smbh_neighbor_number; double tmp, tmp_r, tmp_v; out = fopen("bh_neighbors.dat", "a"); @@ -796,7 +796,7 @@ void write_bh_nb_data() i_bh = 0; - for (int i_bh=0; i_bh < parameters.live_smbh_count; i_bh++) { + for (int i_bh=0; i_bh < config->live_smbh_count; i_bh++) { for (i=0; ieps; + t_end = config->t_end; + dt_disk = config->dt_disk; + dt_contr = config->dt_contr; + dt_bh = config->dt_bh; + eta = config->eta; + strcpy(inp_fname, config->input_file_name.c_str()); + + if (config->binary_smbh_influence_sphere_output) for (i=0; i 0)) { + if (config->live_smbh_output && (config->live_smbh_count > 0)) { out = fopen("bh.dat", "w"); fclose(out); } - if ((parameters.live_smbh_neighbor_output) && (parameters.live_smbh_count > 0)) { + if ((config->live_smbh_neighbor_output) && (config->live_smbh_count > 0)) { out = fopen("bh_neighbors.dat", "w"); fclose(out); } - if (parameters.binary_smbh_influence_sphere_output) { + if (config->binary_smbh_influence_sphere_output) { out = fopen("bbh_inf.dat", "w"); fclose(out); } @@ -1872,11 +1880,11 @@ int main(int argc, char *argv[]) /* Wait to all processors to finish his works... */ MPI_Barrier(MPI_COMM_WORLD); - if (parameters.live_smbh_count == 2) { + if (config->live_smbh_count == 2) { i_bh1 = 0; i_bh2 = 1; - if (parameters.live_smbh_custom_eps >= 0) { + if (config->live_smbh_custom_eps >= 0) { m_bh1 = m[i_bh1]; m_bh2 = m[i_bh2]; @@ -1912,7 +1920,7 @@ int main(int argc, char *argv[]) tmp_i = calc_force_n_BH(m_bh1, x_bh1, v_bh1, m_bh2, x_bh2, v_bh2, - parameters.live_smbh_custom_eps, + config->live_smbh_custom_eps, &pot_bh1, a_bh1, adot_bh1, &pot_bh2, a_bh2, adot_bh2); @@ -1928,7 +1936,7 @@ int main(int argc, char *argv[]) } } - if (parameters.binary_smbh_pn) { + if (config->binary_smbh_pn) { // calculate and "plus" the new BH <-> BH : PN1, PN2, PN2.5, PN3, PN3.5 : acc, jerk @@ -2028,7 +2036,7 @@ int main(int argc, char *argv[]) dt[i] = dt_tmp; - if (parameters.dt_min_warning && (myRank == 0)) { + if (config->dt_min_warning && (myRank == 0)) { if (dt[i] == dt_min) { printf("!!! Warning0: dt = dt_min = %.6E \t ind = %07d \n", dt[i], ind[i]); fflush(stdout); @@ -2037,9 +2045,9 @@ int main(int argc, char *argv[]) } /* i */ - if (parameters.live_smbh_count > 0) { + if (config->live_smbh_count > 0) { double min_dt = *std::min_element(dt, dt+N); - for (int i=0; ilive_smbh_count; i++) dt[i] = min_dt; } /* Wait to all processors to finish his works... */ @@ -2067,10 +2075,10 @@ int main(int argc, char *argv[]) if (myRank == rootRank) { /* Write BH data... */ - if (parameters.live_smbh_output) write_bh_data(); + if (config->live_smbh_output) write_bh_data(); /* Write BH NB data... */ - if (parameters.live_smbh_neighbor_output) write_bh_nb_data(); + if (config->live_smbh_neighbor_output) write_bh_nb_data(); } /* if (myRank == rootRank) */ @@ -2207,7 +2215,7 @@ int main(int argc, char *argv[]) exit(1); } #else - if (parameters.live_smbh_count > 0) { + if (config->live_smbh_count > 0) { i_bh1 = 0; i_bh2 = 1; } @@ -2288,8 +2296,8 @@ int main(int argc, char *argv[]) DT_ACT_REDUCE += (CPU_tmp_user - CPU_tmp_user0); #endif - if (parameters.live_smbh_count == 2) { - if (parameters.live_smbh_custom_eps >= 0) { + if (config->live_smbh_count == 2) { + if (config->live_smbh_custom_eps >= 0) { m_bh1 = m_act[i_bh1]; m_bh2 = m_act[i_bh2]; @@ -2320,7 +2328,7 @@ int main(int argc, char *argv[]) tmp_i = calc_force_n_BH(m_bh1, x_bh1, v_bh1, m_bh2, x_bh2, v_bh2, - parameters.live_smbh_custom_eps, + config->live_smbh_custom_eps, &pot_bh1, a_bh1, adot_bh1, &pot_bh2, a_bh2, adot_bh2); @@ -2334,7 +2342,7 @@ int main(int argc, char *argv[]) adot_act_new[i_bh2] += adot_bh2; } - if (parameters.binary_smbh_pn) { + if (config->binary_smbh_pn) { // calculate and "plus" the new BH <-> BH : PN1, PN2, PN2.5, PN3, PN3.5 : acc, jerk dt_bh_tmp = dt[0]; @@ -2400,7 +2408,7 @@ int main(int argc, char *argv[]) a2dot1abs = a2dot1.norm(); a3dot1abs = a3.norm(); - if ((parameters.live_smbh_count > 0) && (ind_act[i] < parameters.live_smbh_count)) + if ((config->live_smbh_count > 0) && (ind_act[i] < config->live_smbh_count)) dt_new = sqrt(eta_bh*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs)); else dt_new = sqrt(eta*(a1abs*a2dot1abs+adot1abs*adot1abs)/(adot1abs*a3dot1abs+a2dot1abs*a2dot1abs)); @@ -2428,7 +2436,7 @@ int main(int argc, char *argv[]) a_act[i] = a_act_new[i]; adot_act[i] = adot_act_new[i]; - if (parameters.dt_min_warning && (myRank == 0)) { + if (config->dt_min_warning && (myRank == 0)) { if (dt_act[i] == dt_min) { printf("!!! Warning1: dt_act = dt_min = %.6E \t ind_act = %07d \n", dt[i], ind_act[i]); fflush(stdout); @@ -2439,13 +2447,13 @@ int main(int argc, char *argv[]) /* define the min. dt over all the act. part. and set it also for the BH... */ - if (parameters.live_smbh_count > 0) { + if (config->live_smbh_count > 0) { double min_dt = *std::min_element(dt_act, dt_act+n_act); - if (parameters.live_smbh_count>=1) dt_act[i_bh1] = min_dt; - if (parameters.live_smbh_count==2) dt_act[i_bh2] = min_dt; + if (config->live_smbh_count>=1) dt_act[i_bh1] = min_dt; + if (config->live_smbh_count==2) dt_act[i_bh2] = min_dt; } - if (parameters.binary_smbh_influence_sphere_output) { + if (config->binary_smbh_influence_sphere_output) { if (myRank == rootRank) { out = fopen("bbh_inf.dat", "a"); @@ -2476,7 +2484,7 @@ int main(int argc, char *argv[]) iii = ind_act[i]; - if (tmp_r2 < SEMI_a2*SQR(parameters.binary_smbh_influence_radius_factor)) { + if (tmp_r2 < SEMI_a2*SQR(config->binary_smbh_influence_radius_factor)) { if (inf_event[iii] == 0) { @@ -2586,10 +2594,10 @@ int main(int argc, char *argv[]) if (time_cur >= t_bh) { if (myRank == rootRank) { /* Write BH data... */ - if (parameters.live_smbh_output) write_bh_data(); + if (config->live_smbh_output) write_bh_data(); /* Write BH NB data... */ - if (parameters.live_smbh_neighbor_output) write_bh_nb_data(); + if (config->live_smbh_neighbor_output) write_bh_nb_data(); } /* if (myRank == rootRank) */