diff --git a/Makefile b/Makefile index 933616e..792b77e 100644 --- a/Makefile +++ b/Makefile @@ -20,8 +20,12 @@ LIB = $(GRAPELIB) $(CUDALIB) -lm MPICXX ?= mpic++ EXECUTABLE ?= phigrape +# HDF5 +CPPFLAGS += -DHAS_HDF5 +LIB += -lhdf5 -lz -ldl + default: - $(MPICXX) $(CPPFLAGS) $(CXXFLAGS) -DETICS_DTSCF=$(ETICS_DTSCF) $(INC) config.cpp phigrape.cpp -o $(EXECUTABLE) $(LIB) + $(MPICXX) $(CPPFLAGS) $(CXXFLAGS) -DETICS_DTSCF=$(ETICS_DTSCF) $(INC) io.cpp config.cpp phigrape.cpp -o $(EXECUTABLE) $(LIB) yebisu: CPPFLAGS := $(filter-out -DETICS, $(CPPFLAGS)) yebisu: default diff --git a/config.cpp b/config.cpp index 25f79cb..f66a014 100644 --- a/config.cpp +++ b/config.cpp @@ -106,6 +106,10 @@ T Config::get_parameter(Dictionary dictionary, std::string name, T default_value void Config::error_checking() { +#ifndef HAS_HDF5 + if (output_hdf5) + throw std::runtime_error("HDF5 output format (output_hdf5=true) requires the code to be compiled with HAS_HDF5"); +#endif 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)) @@ -148,6 +152,7 @@ Config::Config(std::string file_name) 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"); + output_hdf5 = get_parameter(dictionary, "output_hdf5", false); dt_min_warning = get_parameter(dictionary, "dt_min_warning", false); live_smbh_count = get_parameter(dictionary, "live_smbh_count", 0); diff --git a/config.h b/config.h index 5275bf7..d4b68b6 100644 --- a/config.h +++ b/config.h @@ -14,6 +14,7 @@ public: double dt_bh; double eta; std::string input_file_name; + bool output_hdf5; bool dt_min_warning; int live_smbh_count; double live_smbh_custom_eps; diff --git a/double3.h b/double3.h new file mode 100644 index 0000000..ce0f3ef --- /dev/null +++ b/double3.h @@ -0,0 +1,82 @@ +#pragma once +#include + +struct double3 { + double data[3]; + double3() {} + double3(const double x, const double y, const double z) + { + data[0] = x; + data[1] = y; + data[2] = z; + } + double& operator[](int i) {return data[i];} + const double operator[](int i) const {return data[i];} + double3& operator=(const double3& a) + { + data[0] = a.data[0]; + data[1] = a.data[1]; + data[2] = a.data[2]; + return *this; + } + double3& operator+=(const double3& a) + { + data[0] += a.data[0]; + data[1] += a.data[1]; + data[2] += a.data[2]; + return *this; + } + double3& operator-=(const double3& a) + { + data[0] -= a.data[0]; + data[1] -= a.data[1]; + data[2] -= a.data[2]; + return *this; + } + double3& operator/=(const double& c) + { + data[0] /= c; + data[1] /= c; + data[2] /= c; + return *this; + } + 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]); + } + operator double*() {return data;} +}; + +inline double3 operator*(const double& c, const double3& a) +{ + return double3(a.data[0]*c, a.data[1]*c, a.data[2]*c); +} + +inline double3 operator*(const double3& a, const double& c) +{ + return double3(a.data[0]*c, a.data[1]*c, a.data[2]*c); +} + +inline 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]; +} + +inline double3 operator/(const double3& a, const double& c) +{ + return double3(a.data[0]/c, a.data[1]/c, a.data[2]/c); +} + +inline double3 operator+(const double3& a, const double3& b) +{ + return double3(a.data[0]+b.data[0], a.data[1]+b.data[1], a.data[2]+b.data[2]); +} + +inline double3 operator-(const double3& a, const double3& b) +{ + return double3(a.data[0]-b.data[0], a.data[1]-b.data[1], a.data[2]-b.data[2]); +} diff --git a/io.cpp b/io.cpp index dce1d55..122415e 100644 --- a/io.cpp +++ b/io.cpp @@ -1,71 +1,79 @@ +#ifdef HAS_HDF5 #include "hdf5.h" +#endif + +#include "double3.h" +#include +#include #include -#include #include #include +#include #include -struct double3 {double x,y,z; -}; + + +bool is_hdf5(std::string file_name) +{ + std::ifstream file(file_name, std::ifstream::binary); + const char hdf5_magic[] = "\x89HDF\x0d\x0a\x1a\x0a"; + char buffer[8]; + file.read(buffer, 8); + if (!file) throw std::runtime_error("Unable to read file " + file_name); + bool result = (memcmp(buffer, hdf5_magic, 8)==0); + file.close(); + return result; +} + +void ascii_read(const std::string file_name, int *step_num, int *N, double *t, double *m, double3 *x, double3 *v) +{ + char rest[512]; + int result; + + std::ifstream file(file_name); + if (!file.good()) throw std::runtime_error("File " + file_name + " not found."); + std::string str; + + std::getline(file, str); + result = sscanf(str.c_str(), "%d%s", step_num, rest); + if (result!=1) throw std::runtime_error("Error parsing line 1: expected one integer"); + + std::getline(file, str); + result = sscanf(str.c_str(), "%d%s", N, rest); + if (result!=1) throw std::runtime_error("Error parsing line 2: expected one integer"); + + std::getline(file, str); + result = sscanf(str.c_str(), "%lf%s", t, rest); + if (result!=1) throw std::runtime_error("Error parsing line 3: expected one real number"); + + int i = -1; + while (std::getline(file, str)) { + if (++i > *N) throw std::runtime_error("Error parsing line " + std::to_string(i+4) + ": particle out of range"); + result = sscanf(str.c_str(), "%*s %lf %lf %lf %lf %lf %lf %lf%s", &m[i], &x[i][0], &x[i][1], &x[i][2], &v[i][0], &v[i][1], &v[i][2], rest); + } + file.close(); +} void ascii_write(const std::string file_name, const int step_num, const int N, const double t, const double *m, const double3 *x, const double3 *v, int precision=10) { + auto file = std::ofstream(file_name); + if (!file.is_open()) throw std::runtime_error("Cannot open file for output"); int id_width = (int)log10(N-1) + 1; - char string_template[256]; + char string_template[256], output_string[256]; + file << step_num << '\n'; + file << N << '\n'; + file.precision(16); + file << std::scientific << t << '\n'; sprintf(string_template, "%%0%dd%%%d.%dE%%%d.%dE%%%d.%dE%%%d.%dE%%%d.%dE%%%d.%dE%%%d.%dE\n", id_width, precision+7, precision, precision+8, precision, precision+8, precision, precision+8, precision, precision+8, precision, precision+8, precision, precision+8, precision); - printf("%d\n", step_num); - printf("%d\n", N); - printf("%.16E\n", t); for (int i=0; i>1)%2; - bool write_jrk = (write_mode>>2)%2; - - hid_t file_id, group_id, attribute_id, dataset_id, dataspace_id; - hsize_t dims[2] = {(hsize_t)N, 3}; - - hid_t h5_float_type; - if (use_double_precision) h5_float_type = H5T_IEEE_F64LE; - else h5_float_type = H5T_IEEE_F32LE; - - file_id = H5Fcreate(file_name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); - char group_name[32]; - sprintf(group_name, "/Step#%d", step_num); - group_id = H5Gcreate2(file_id, group_name, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); - dataspace_id = H5Screate(H5S_SCALAR); - attribute_id = H5Acreate2 (group_id, "Time", H5T_IEEE_F64LE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT); - H5Awrite(attribute_id, H5T_NATIVE_DOUBLE, &t); - H5Sclose(dataspace_id); - - auto write_dataset = [&](const char dataset_name[], int ndims, double *data) { - hid_t dataspace_id = H5Screate_simple(ndims, dims, NULL); - char dataset_path[32]; - sprintf(dataset_path, "%s/%s", group_name, dataset_name); - hid_t dataset_id = H5Dcreate2(file_id, dataset_path, h5_float_type, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); - H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, data); - H5Dclose(dataset_id); - H5Sclose(dataspace_id); - }; - - write_dataset("MASS", 1, (double*)m); - write_dataset("X", 2, (double*)x); - write_dataset("V", 2, (double*)v); - if (write_pot) write_dataset("POT", 1, (double*)pot); - if (write_acc) write_dataset("ACC", 2, (double*)acc); - if (write_jrk) write_dataset("JRK", 2, (double*)jrk); - - H5Gclose(group_id); - H5Fclose(file_id); -} - - void h5_read(const std::string file_name, int *step_num, int *N, double *t, double m[], double3 x[], double3 v[]) { +#ifdef HAS_HDF5 // Open file and root group; count number of top level objects hid_t file_id = H5Fopen(file_name.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); hid_t group_id = H5Gopen2(file_id, "/", H5P_DEFAULT); @@ -91,6 +99,7 @@ void h5_read(const std::string file_name, int *step_num, int *N, double *t, doub sprintf(path, "/Step#%d", *step_num); attr_id = H5Aopen_by_name(file_id, path, "Time", H5P_DEFAULT, H5P_DEFAULT); H5Aread(attr_id, H5T_NATIVE_DOUBLE, t); + H5Aclose(attr_id); sprintf(path, "/Step#%d/MASS", *step_num); dataset_id = H5Dopen2(file_id, path, H5P_DEFAULT); @@ -100,6 +109,8 @@ void h5_read(const std::string file_name, int *step_num, int *N, double *t, doub throw std::runtime_error("Dataset MASS in Step#" + std::to_string(*step_num) + " of file " + file_name + " is " + std::to_string(ndims) + "-dimensional (expected 1-dimensional)"); H5Sget_simple_extent_dims(dataspace_id, dims, NULL); H5Dread(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, m); + H5Sclose(dataspace_id); + H5Dclose(dataset_id); *N = dims[0]; sprintf(path, "/Step#%d/X", *step_num); @@ -110,6 +121,8 @@ void h5_read(const std::string file_name, int *step_num, int *N, double *t, doub H5Sget_simple_extent_dims(dataspace_id, dims, NULL); if ((dims[0] != *N) || (dims[1] != 3)) throw std::runtime_error("Bad dimensionality"); H5Dread(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, (double*)x); + H5Sclose(dataspace_id); + H5Dclose(dataset_id); sprintf(path, "/Step#%d/V", *step_num); dataset_id = H5Dopen2(file_id, path, H5P_DEFAULT); @@ -119,36 +132,65 @@ void h5_read(const std::string file_name, int *step_num, int *N, double *t, doub H5Sget_simple_extent_dims(dataspace_id, dims, NULL); if ((dims[0] != *N) || (dims[1] != 3)) throw std::runtime_error("Bad dimensionality"); H5Dread(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, (double*)v); + H5Sclose(dataspace_id); + H5Dclose(dataset_id); + + H5Gclose(group_id); + H5Fclose(file_id); +#else + throw std::runtime_error("h5_read was called but compiled without HDF5 support"); +#endif } - -int main() +void h5_write(const std::string file_name, const int step_num, const int N, const double t, const double *m, const double3 *x, const double3 *v, const double *pot, const double3 *acc, const double3 *jrk, const int write_mode=0, const bool use_double_precision=true) { -// double3 aaa; -// double pot; - std::string file_name = "new_file.h5"; - int N = 25; - double3 x[N]; - for (int i=0; i> 1) % 2; + bool write_jrk = (write_mode >> 2) % 2; + if (write_pot) write_dataset("POT", 1, (double*)pot); + if (write_acc) write_dataset("ACC", 2, (double*)acc); + if (write_jrk) write_dataset("JRK", 2, (double*)jrk); + + H5Gclose(group_id); + H5Fclose(file_id); + H5close(); // If we don't do that (HDF5 1.10.5) then the file isn't really closed... +#else + throw std::runtime_error("h5_write was called but compiled without HDF5 support"); +#endif } diff --git a/io.h b/io.h new file mode 100644 index 0000000..1baa152 --- /dev/null +++ b/io.h @@ -0,0 +1,16 @@ +#pragma once +#include +#include "double3.h" + +bool is_hdf5(std::string file_name); +// This function is implemented independently of the HDF5 library + +void ascii_read(const std::string file_name, int *step_num, int *N, double *t, double *m, double3 *x, double3 *v); + +void ascii_write(const std::string file_name, const int step_num, const int N, const double t, const double *m, const double3 *x, const double3 *v, int precision=10); + +void h5_read(const std::string file_name, int *step_num, int *N, double *t, double m[], double3 x[], double3 v[]); +// In case the code is compiled without HDF5 support, the implementation of this function just throws an error + +void h5_write(const std::string file_name, const int step_num, const int N, const double t, const double *m, const double3 *x, const double3 *v, const double *pot, const double3 *acc, const double3 *jrk, const int write_mode=0, const bool use_double_precision=true); +// In case the code is compiled without HDF5 support, the implementation of this function just throws an error diff --git a/phigrape.conf b/phigrape.conf index c9394b4..0ceaa73 100644 --- a/phigrape.conf +++ b/phigrape.conf @@ -21,11 +21,11 @@ dt_bh = 0.125 eta = 0.01 # Name of the input file; use "data.con" in most cases -inp_data = data.con +input_file_name = data.con +output_hdf5 = true ##### NOT IMPLEMENTED ####################### -output_format = HDF5 dt_min_warning = false ############################################# @@ -33,9 +33,9 @@ dt_min_warning = false # EXTERNAL GRAVITY # #################### -ext_units_physical = true -unit_mass = 1e5 -unit_length = 10 +ext_units_physical = false +#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! @@ -47,9 +47,9 @@ ext_a_disk = 0 ext_b_disk = 0 ext_m_halo_plummer = 0 ext_b_halo_plummer = 0 -ext_log_halo_v = 30 +ext_log_halo_v = 0 # here the length scale for the logarithmic halo is in kpc, not pc like in the old phigrape. -ext_log_halo_r = 0.04 +ext_log_halo_r = 0 diff --git a/phigrape.cpp b/phigrape.cpp index 04059b8..81eda25 100644 --- a/phigrape.cpp +++ b/phigrape.cpp @@ -53,7 +53,10 @@ Coded by : Peter Berczik Version number : 19.04 Last redaction : 2019.04.16 12:55 *****************************************************************************/ +#include "double3.h" #include "config.h" +#include "io.h" + Config *config; #define NORM // Physical normalization @@ -177,84 +180,6 @@ Config *config; #define N_MAX (6*MB) #define N_MAX_loc (2*MB) -struct double3 { - double data[3]; - double3() {} - double3(const double x, const double y, const double z) - { - data[0] = x; - data[1] = y; - data[2] = z; - } - double& operator[](int i) {return data[i];} - double3& operator=(const double3& a) - { - data[0] = a.data[0]; - data[1] = a.data[1]; - data[2] = a.data[2]; - return *this; - } - double3& operator+=(const double3& a) - { - data[0] += a.data[0]; - data[1] += a.data[1]; - data[2] += a.data[2]; - return *this; - } - double3& operator-=(const double3& a) - { - data[0] -= a.data[0]; - data[1] -= a.data[1]; - data[2] -= a.data[2]; - return *this; - } - double3& operator/=(const double& c) - { - data[0] /= c; - data[1] /= c; - data[2] /= c; - return *this; - } - 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]); - } - operator double*() {return data;} -}; - -double3 operator*(const double& c, const double3& a) -{ - return double3(a.data[0]*c, a.data[1]*c, a.data[2]*c); -} - -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); -} - -double3 operator+(const double3& a, const double3& b) -{ - return double3(a.data[0]+b.data[0], a.data[1]+b.data[1], a.data[2]+b.data[2]); -} - -double3 operator-(const double3& a, const double3& b) -{ - return double3(a.data[0]-b.data[0], a.data[1]-b.data[1], a.data[2]-b.data[2]); -} @@ -366,16 +291,6 @@ void get_CPU_time(double *time_real, double *time_user, double *time_syst) *time_user = *time_real; } -void read_data(char inp_fname[], int *diskstep, int *N, double *time_cur, int ind[], double m[], double3 x[], double3 v[]) -{ - auto inp = fopen(inp_fname, "r"); - fscanf(inp, "%d \n", diskstep); - fscanf(inp, "%d \n", N); - fscanf(inp, "%lE \n", time_cur); - for (int i=0; i<*N; i++) fscanf(inp,"%d %lE %lE %lE %lE %lE %lE %lE \n", &ind[i], &m[i], &x[i][0], &x[i][1], &x[i][2], &v[i][0], &v[i][1], &v[i][2]); - fclose(inp); -} - void write_snap_data(char out_fname[], int diskstep, int N, double time_cur, int ind[], double m[], double3 x[], double3 v[]) { auto out = fopen(out_fname, "w"); @@ -962,6 +877,17 @@ int main(int argc, char *argv[]) C_NB = config->pn_c; std::copy(config->pn_usage.begin(), config->pn_usage.end(), usedOrNot); + if (is_hdf5(config->input_file_name)) { +#ifndef HAS_HDF5 + fprintf(stderr, "ERROR: input file is in HDF5 format, but the code was compiled without HDF5 support\n") + return -1; +#endif + h5_read(config->input_file_name, &diskstep, &N, &time_cur, m, x, v); + } + else + ascii_read(config->input_file_name, &diskstep, &N, &time_cur, m, x, v); + std::iota(ind, ind+N, 0); + if (myRank == rootRank) { //TODO move it out of (myRank == rootRank) so you don't need to communicate them. @@ -995,12 +921,6 @@ int main(int argc, char *argv[]) fscanf(inp,"%lE %lE %lE", &m_ext, &r_ext, &g_ext); #endif - /* read the global data for particles to the rootRank */ - - read_data(inp_fname, &diskstep, &N, &time_cur, ind, m, x, v); - - /* possible coordinate & velocity limits for ALL particles !!! */ - printf("\n"); printf("Begin the calculation of phi-GRAPE program on %03d processors\n", n_proc); printf("\n"); @@ -1994,7 +1914,8 @@ int main(int argc, char *argv[]) /* write cont data */ - write_snap_data((char*)"data.con", diskstep, N, time_cur, ind, m, x, v); + if (config->output_hdf5) h5_write("data.con", diskstep, N, time_cur, m, x, v, pot, a, adot, 0, true); + else ascii_write("data.con", diskstep, N, time_cur, m, x, v, 16); /* possible OUT for timing !!! */ @@ -2049,8 +1970,10 @@ int main(int argc, char *argv[]) if (myRank == rootRank) { diskstep++; char out_fname[256]; - sprintf(out_fname, "%06d.dat", diskstep); - write_snap_data(out_fname, diskstep, N, time_cur, ind, m, x, v); + sprintf(out_fname, "%06d", diskstep); + if (config->output_hdf5) h5_write(std::string(out_fname) + ".h5", diskstep, N, time_cur, m, x, v, pot, a, adot, 0, true); + else ascii_write(std::string(out_fname) + ".dat", diskstep, N, time_cur, m, x, v, 10); + // TODO custom precision } /* if (myRank == rootRank) */ #ifdef ETICS_DUMP