diff --git a/.gitignore b/.gitignore index 773bfd7..9a8c7cc 100644 --- a/.gitignore +++ b/.gitignore @@ -8,3 +8,4 @@ *.mask grapite-dev-exec-threshold phigrape +*.h5 diff --git a/config.cpp b/config.cpp index f66a014..ce229d9 100644 --- a/config.cpp +++ b/config.cpp @@ -151,8 +151,14 @@ Config::Config(std::string file_name) 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"); + output_hdf5 = get_parameter(dictionary, "output_hdf5", false); + output_hdf5_double_precision = get_parameter(dictionary, "output_hdf5_double_precision", true); + output_ascii_precision = get_parameter(dictionary, "output_ascii_precision", 10); + output_extra_mode = get_parameter(dictionary, "output_extra_mode", 10); + 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 d4b68b6..832d517 100644 --- a/config.h +++ b/config.h @@ -15,6 +15,9 @@ public: double eta; std::string input_file_name; bool output_hdf5; + bool output_hdf5_double_precision; + int output_ascii_precision; + int output_extra_mode; bool dt_min_warning; int live_smbh_count; double live_smbh_custom_eps; diff --git a/io.cpp b/io.cpp index 122415e..a2c0d3d 100644 --- a/io.cpp +++ b/io.cpp @@ -142,7 +142,7 @@ void h5_read(const std::string file_name, int *step_num, int *N, double *t, doub #endif } -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) +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 extra_mode=0, const bool use_double_precision=true) { #ifdef HAS_HDF5 hid_t file_id, group_id, attribute_id, dataset_id, dataspace_id; @@ -151,13 +151,7 @@ void h5_write(const std::string file_name, const int step_num, const int N, cons hid_t h5_float_type; if (use_double_precision) h5_float_type = H5T_IEEE_F64LE; else h5_float_type = H5T_IEEE_F32LE; - -// static int zzz = 5; file_id = H5Fcreate(file_name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); - -// file_id = H5Fcreate(std::to_string(zzz).c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); - -// zzz++; char group_name[32]; sprintf(group_name, "/Step#%d", step_num); group_id = H5Gcreate2(file_id, group_name, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); @@ -180,9 +174,9 @@ void h5_write(const std::string file_name, const int step_num, const int N, cons write_dataset("X", 2, (double*)x); write_dataset("V", 2, (double*)v); - bool write_pot = (write_mode ) % 2; - bool write_acc = (write_mode >> 1) % 2; - bool write_jrk = (write_mode >> 2) % 2; + bool write_pot = (extra_mode ) % 2; + bool write_acc = (extra_mode >> 1) % 2; + bool write_jrk = (extra_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); diff --git a/phigrape.conf b/phigrape.conf index 0ceaa73..acb3a82 100644 --- a/phigrape.conf +++ b/phigrape.conf @@ -8,8 +8,8 @@ eps = 1E-4 # End time of the calculation t_end = 0.25 -# Interval of snapshot files output (xxxxxx.dat) -dt_disk = 1.0 +# Interval of snapshot files output +dt_disk = 0.125 # Interval for the energy control output (contr.dat) dt_contr = 0.125 @@ -20,100 +20,132 @@ dt_bh = 0.125 # Parameter for timestep determination eta = 0.01 -# Name of the input file; use "data.con" in most cases +# Name of the input file; use "data.con" in most cases [default: data.con] input_file_name = data.con -output_hdf5 = true -##### NOT IMPLEMENTED ####################### -dt_min_warning = false -############################################# +########## +# OUTPUT # +########## + +# Whether to use HDF5 format for snapshot and restart; regular ASCII snapshorts are saved if false [default: false] +#output_hdf5 = true + +# If using HDF5 output, use double precision or not [default: true] +# Consider setting to false to save disk space. Restart file is always saved in double precision. +#output_hdf5_double_precision = true + +# If using ASCII output, the number of digits after the decimal point [default: 10] +# Restart file is saved with 16 digits after the decimal point. +#output_ascii_precision = 6 + +# Extra output: optionally save potential, acceleration and jerk in snapshot files [default: 0] +# This is a number between 0 and 7 that encodes the output options in the following way: +# [value] = [save jerk]*4 + [save acceleration]*2 + [save potential] +# Example: choose 5 if it is needed to save the jerk and the potential, but not the acceleration for some reason. +# Currently implemented in HDF5 output only. +#output_extra_mode = 7 + +# Whether to output a warning on the screen when the minimum time step is encountered. [default: false] +#dt_min_warning = false #################### # EXTERNAL GRAVITY # #################### -ext_units_physical = false -#unit_mass = 1e5 -#unit_length = 10 -# TODO add the option to normalize using other units +# Remember that external gravity models are applied at the same coordinate system as the particles. If the idea is to simulate a globular cluster orbiting in an external field, be sure to set the initial conditions appropriately (applying a shift to the coordinates and velocities). -# Notice that if physical units are used, the "a" and "b" parameters are in kiloparsec, not parsec! +# Whether the parameters for the external gravitational field given below are in physical units or Hénon units. If true, the system used is {kiloparsec, solar mass, kilometre per second} [default: false] +#ext_units_physical = true -ext_m_bulge = 0 -ext_b_bulge = 0 -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 = 0 -# here the length scale for the logarithmic halo is in kpc, not pc like in the old phigrape. -ext_log_halo_r = 0 +# If Physical units were selected, specify the simulation's unit mass (is solar masses) and unit lenght (in parsec; not kiloparsec) +# TODO: add the option to normalize using other units. +#unit_mass = 4E5 # MSun +#unit_length = 15 # pc + + +# The bulge is a Plummer potential with the following total mass and radius. +#ext_m_bulge = 5E9 # MSun +#ext_b_bulge = 1.9 # kpc + + +# The disk is a Miyamoto-Nagai potential with the following total mass, scale length, and scale height +#ext_m_disk = 6.8E10 # MSun +#ext_a_disk = 3.00 # kpc +#ext_b_disk = 0.28 # kpc + + +# This halo option is yet another Plummer potential with the following total mass and radius. +#ext_m_halo_plummer = 8E11 # MSun +#ext_b_halo_plummer = 245 # kpc + + +# This halo option is a logarithmic potential with the following velocity and radius parameters. +#ext_log_halo_v = 240 # km/s +#ext_log_halo_r = 1 # kpc -#################################### -## 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 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 = 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 + +# Output additional diagnostics about live SMBHs. [default: false] #live_smbh_output = true -# -## Output additional diagnostics about the SMBH's (or SMBHs') nearest neighbours (number could be set as shown below). [default: false] + +# 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] + +# 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 SUPERMASSIVE BLACK HOLE # +################################## + +# The following parameters can be set when live_smbh_count is 2. + +# Output additional diagnostics about the 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] + +# 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] +#binary_smbh_influence_radius_factor = 10 + +# 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} + +# 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} [default: {1,1,1,1,1,1,1}] #pn_usage = {1, 1, 1, 1, 0, 0, 0} -# -## The speed of light in N-body units [default: 500] + +# 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 # -##################################### + +#################################### +# 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 81eda25..91143f9 100644 --- a/phigrape.cpp +++ b/phigrape.cpp @@ -291,22 +291,6 @@ void get_CPU_time(double *time_real, double *time_user, double *time_syst) *time_user = *time_real; } -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"); - fprintf(out,"%06d \n", diskstep); - fprintf(out,"%07d \n", N); - fprintf(out,"%.10E \n", time_cur); - for (int i=0; ilive_smbh_count == 2) { @@ -1971,9 +1955,8 @@ int main(int argc, char *argv[]) diskstep++; char out_fname[256]; 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 (config->output_hdf5) h5_write(std::string(out_fname) + ".h5", diskstep, N, time_cur, m, x, v, pot, a, adot, config->output_extra_mode, config->output_hdf5_double_precision); + else ascii_write(std::string(out_fname) + ".dat", diskstep, N, time_cur, m, x, v, config->output_ascii_precision); } /* if (myRank == rootRank) */ #ifdef ETICS_DUMP