Now config file handles most parameters (but not external potential)

This commit is contained in:
Yohai Meiron 2020-03-17 12:24:53 -04:00
parent bb9a343763
commit 9296db0609
6 changed files with 188 additions and 165 deletions

View file

@ -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

126
config.cpp Normal file
View file

@ -0,0 +1,126 @@
#include "config.h"
#include <stdexcept>
#include <cstdio>
#include <fstream>
using Dictionary = std::unordered_map<std::string,std::string>;
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<std::string,std::string> 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<typename T>
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<T>((*item).second);
}
// For optional parameters
template<typename T>
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<T>((*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<double>(dictionary, "eps");
t_end = get_parameter<double>(dictionary, "t_end");
dt_disk = get_parameter<double>(dictionary, "dt_disk");
dt_contr = get_parameter<double>(dictionary, "dt_contr");
dt_bh = get_parameter<double>(dictionary, "dt_bh", dt_contr);
eta = get_parameter<double>(dictionary, "eta");
input_file_name = get_parameter<std::string>(dictionary, "input_file_name", "data.con");
dt_min_warning = get_parameter<bool>(dictionary, "dt_min_warning", false);
live_smbh_count = get_parameter<int>(dictionary, "live_smbh_count", 0);
live_smbh_custom_eps = get_parameter<double>(dictionary, "live_smbh_custom_eps", -1);
live_smbh_output = get_parameter<bool>(dictionary, "live_smbh_output", false);
live_smbh_neighbor_output = get_parameter<bool>(dictionary, "live_smbh_neighbor_output", false);
live_smbh_neighbor_number = get_parameter<int>(dictionary, "live_smbh_neighbor_number", 10);
binary_smbh_pn = get_parameter<bool>(dictionary, "binary_smbh_pn", false);
binary_smbh_influence_sphere_output = get_parameter<bool>(dictionary, "binary_smbh_influence_sphere_output", false);
binary_smbh_influence_radius_factor = get_parameter<double>(dictionary, "binary_smbh_influence_radius_factor", 10.);
error_checking();
}

View file

@ -1,7 +1,11 @@
#pragma once
#include <string>
#include <unordered_map>
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<std::string,std::string>;
Dictionary read_config_file(const std::string file_name);
std::string strip(const std::string str);
template<typename T> T string_cast(const std::string str);
template<typename T> T get_parameter(Dictionary dictionary, std::string name);
template<typename T> T get_parameter(Dictionary dictionary, std::string name, T default_value);
void error_checking();
};

123
io.cpp
View file

@ -1,123 +0,0 @@
#include "phigrape.h"
#include <string>
#include <unordered_map>
#include <stdexcept>
#include <cstdio>
#include <fstream>
using Dictionary = std::unordered_map<std::string,std::string>;
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<std::string,std::string> 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<typename T> 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<typename T>
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<T>((*item).second);
}
// For optional parameters
template<typename T>
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<T>((*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<double>(dictionary, "eps");
parameters.t_end = get_parameter<double>(dictionary, "t_end");
parameters.dt_disk = get_parameter<double>(dictionary, "dt_disk");
parameters.dt_contr = get_parameter<double>(dictionary, "dt_contr");
parameters.dt_bh = get_parameter<double>(dictionary, "dt_bh", parameters.dt_contr);
parameters.eta = get_parameter<double>(dictionary, "eta");
parameters.input_file_name = get_parameter<std::string>(dictionary, "dt_bh", "data.con");
parameters.dt_min_warning = get_parameter<bool>(dictionary, "dt_min_warning", false);
parameters.live_smbh_count = get_parameter<int>(dictionary, "live_smbh_count", 0);
parameters.live_smbh_custom_eps = get_parameter<double>(dictionary, "live_smbh_custom_eps", -1);
parameters.live_smbh_output = get_parameter<bool>(dictionary, "live_smbh_output", false);
parameters.live_smbh_neighbor_output = get_parameter<bool>(dictionary, "live_smbh_neighbor_output", false);
parameters.live_smbh_neighbor_number = get_parameter<int>(dictionary, "live_smbh_neighbor_number", 10);
parameters.binary_smbh_pn = get_parameter<bool>(dictionary, "binary_smbh_pn", false);
parameters.binary_smbh_influence_sphere_output = get_parameter<bool>(dictionary, "binary_smbh_influence_sphere_output", false);
parameters.binary_smbh_influence_radius_factor = get_parameter<double>(dictionary, "binary_smbh_influence_radius_factor", 10.);
}

View file

@ -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

View file

@ -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; i<N; i++) var_sort[i] = (x[i]-x[i_bh]).norm();
std::iota(ind_sort, ind_sort+N, 0);
std::partial_sort(ind_sort, ind_sort + nb, ind_sort + N, [&](int i, int j) {return var_sort[i] < var_sort[j];});
@ -1511,12 +1511,20 @@ int main(int argc, char *argv[])
/* read the input parameters to the rootRank */
if (myRank == rootRank) {
inp = fopen("phi-GRAPE.cfg","r");
fscanf(inp,"%lE %lE %lE %lE %lE %lE %s", &eps, &t_end, &dt_disk, &dt_contr, &dt_bh, &eta, inp_fname);
fclose(inp);
config = new Config("phigrape.conf");
if (parameters.binary_smbh_influence_sphere_output) for (i=0; i<N; i++) inf_event[i] = 0;
if (myRank == rootRank) {
//TODO move it out of (myRank == rootRank) so you don't need to communicate them.
eps = config->eps;
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<N; i++) inf_event[i] = 0; // WARNING N wasn't set yet!
/*
eps : Plummer softening parameter (can be even 0)
@ -1634,15 +1642,15 @@ int main(int argc, char *argv[])
fclose(out);
#endif
if (parameters.live_smbh_output && (parameters.live_smbh_count > 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; i<parameters.live_smbh_count; i++) dt[i] = min_dt;
for (int i=0; i<config->live_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) */