diff --git a/.gitignore b/.gitignore index d0c6438..32f95dc 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,6 @@ .* -main -libmain.so +scattering_experiment +libreplay.so data __pycache__ *.png diff --git a/Makefile b/Makefile index 94d4882..59ad505 100644 --- a/Makefile +++ b/Makefile @@ -1,9 +1,9 @@ OPTIMIZATION ?= 3 -CXXFLAGS += -O$(OPTIMIZATION) -pthread -INC += -I/home/meiron/boost_1_72_0 -LIB += -lgsl +CXXFLAGS += --std=c++17 -O$(OPTIMIZATION) +INC += +LIB += -lgsl -pthread -EXECUTABLE ?= main +EXECUTABLE ?= scattering_experiment defualt: loadtxt.o replay.o scattering_experiment.o $(CXX) $(CPPFLAGS) $(CXXFLAGS) $(INC) $? -o $(EXECUTABLE) $(LIB) diff --git a/scattering_experiment.cpp b/scattering_experiment.cpp index 7d43715..f390daa 100644 --- a/scattering_experiment.cpp +++ b/scattering_experiment.cpp @@ -1,6 +1,8 @@ #include #include +#include #include +#include #include #include #include @@ -9,6 +11,8 @@ #include // We just want to print the timestamp #include +using double3 = Eigen::Vector3d; + extern "C" void integrate(const double y0[], const double t_max, const double stride_size, double y[]); @@ -36,32 +40,82 @@ std::string iso_time(std::chrono::system_clock::time_point time) return std::string(time_string); } -class Initial_conditions { +class In_box_cyl { public: - Initial_conditions(unsigned long long seed, double x_max, double v_max) + In_box_cyl(const double pos_tol, const double vel_tol, const double3 target_pos, const double3 target_vel) + : pos_tol(pos_tol), vel_tol(vel_tol) + { + target_d = sqrt(target_pos[0]*target_pos[0] + target_pos[1]*target_pos[1]); + target_z = target_pos[2]; + target_vd = (target_pos[0]*target_vel[0] + target_pos[1]*target_vel[1])/target_d; + target_vt = (target_pos[0]*target_vel[1] - target_pos[1]*target_vel[0])/target_d; + target_vz = target_vel[2]; + } + bool operator()(const double w[6]) { + const double d = sqrt(w[0]*w[0] + w[1]*w[1]); + const double& z = w[2]; + const double vd = (w[0]*w[3] + w[1]*w[4])/d; + const double vt = (w[0]*w[4] - w[1]*w[3])/d; + const double& vz = w[5]; + return (abs(d - target_d) < pos_tol) + && (abs(z - target_z) < pos_tol) + && (abs(vd - target_vd) < vel_tol) + && (abs(vt - target_vt) < vel_tol) + && (abs(vz - target_vz) < vel_tol); + } + friend std::ostream& operator<<(std::ostream& stream, const In_box_cyl& obj) + { + stream.setf(std::ios::scientific); + stream << "# Acceptor type: box in cylindrical coordinates\n"; + stream << "# Acceptor parameters:\n"; + stream << "# pos_tol = " << std::setw(13) << obj.pos_tol << '\n'; + stream << "# vel_tol = " << std::setw(13) << obj.vel_tol << '\n'; + stream << "# target_d = " << std::setw(13) << obj.target_d << '\n'; + stream << "# target_z = " << std::setw(13) << obj.target_z << '\n'; + stream << "# target_vd = " << std::setw(13) << obj.target_vd << '\n'; + stream << "# target_vt = " << std::setw(13) << obj.target_vt << '\n'; + stream << "# target_vz = " << std::setw(13) << obj.target_vz << '\n'; + } +private: + double pos_tol, vel_tol; + double target_d, target_z, target_vd, target_vt, target_vz; +}; + +class Box_initial_conditions { +public: + Box_initial_conditions(double x_max, double v_max) : x_max(x_max), v_max(v_max) { - generator = std::default_random_engine(seed); - position_distribution = std::uniform_real_distribution(-x_max, x_max); - velocity_distribution = std::uniform_real_distribution(-v_max, v_max); + position_distribution = std::uniform_real_distribution(-x_max, x_max); + velocity_distribution = std::uniform_real_distribution(-v_max, v_max); } - std::array generate() + template + std::array operator()(Generator& generator) { std::array result; for (int i=0; i<3; i++) result[i] = position_distribution(generator); for (int i=3; i<6; i++) result[i] = velocity_distribution(generator); return result; } + friend std::ostream& operator<<(std::ostream& stream, const Box_initial_conditions& obj) + { + stream.setf(std::ios::scientific); + stream << "# Initial conditions type: uniform box\n"; + stream << "# Initial conditions parameters:\n"; + stream << "# x_max = " << std::setw(13) << obj.x_max << '\n'; + stream << "# v_max = " << std::setw(13) << obj.v_max << '\n'; + } private: double x_max, v_max; - std::default_random_engine generator; - std::uniform_real_distribution position_distribution, velocity_distribution; + std::uniform_real_distribution<> position_distribution, velocity_distribution; }; + +template class Scattering_experiment { public: - Scattering_experiment(std::string file_name) - : bunch_size(16384), counter_target(10000000L), report_interval_seconds(3) + Scattering_experiment(const std::string file_name, const double t_max, const unsigned long long n_experiments, Initial_conditions& initial_conditions, Acceptor& accept) + : bunch_size(256), n_experiments(n_experiments), report_interval_seconds(300), t_max(t_max), accept(accept), initial_conditions(initial_conditions) { auto time = std::chrono::system_clock::now(); file = std::ofstream(file_name/*, std::ofstream::app*/); @@ -79,36 +133,40 @@ public: const unsigned long big_prime = 840580612873L; unique_id = time_count + big_prime*pid; // TODO hash it file << "# Unique identifier: " << unique_id << '\n'; - n_threads = 4; + n_threads = std::thread::hardware_concurrency(); file << "# Concurrency: " << n_threads << '\n'; - // TODO print information about target and galaxy model + file << initial_conditions; + file << accept; file.flush(); + file.setf(std::ios::scientific | std::ios::right); + file.precision(16); + // TODO information about galaxy model } ~Scattering_experiment() { file.close(); } - void output_data(int aaa) - { - std::lock_guard lock(mtx); - file << aaa << '\n'; - file.flush(); - } void thread_task(int tid) { auto thread_begin_time = std::chrono::system_clock::now(); auto report_interval = std::chrono::seconds(report_interval_seconds); auto next_report_time = thread_begin_time + report_interval; - // std::default_random_engine generator(unique_id + tid); - // std::uniform_real_distribution position_distribution(0, 1000); - // auto random_position_component = [&]() { return position_distribution(generator); }; - Initial_conditions initial_conditions(unique_id + tid, 100, 100); + std::default_random_engine generator(unique_id + tid); - while (counter < counter_target) { + while (counter < n_experiments) { for (int i=0; i lock(mtx); + char buffer[512]; + file << std::setw(3) << tid; + for (int j=0; j<12; j++) file << std::setw(24) << y[j]; + file << std::endl; + } } counter += bunch_size; if (tid == 0) { @@ -137,21 +195,30 @@ public: int n_threads; std::vector threads; std::atomic counter; - unsigned long long counter_target; + unsigned long long n_experiments; unsigned long unique_id; int bunch_size; int report_interval_seconds; - - int i; + double t_max; + Acceptor accept; + Initial_conditions initial_conditions; }; - - int main() { std::cout << "Hi\n"; - Scattering_experiment scattering_experiment("aaa.out"); + const double t_max = 11.65551390; // Gyr + double3 target_pos = {-9.1817914423447837E+03, -2.3161972143758221E+03, 4.2321813890923086E+03}; // pc + double3 target_vel = { 1.6995805137819784E+02, 1.3476658021349482E+02, -1.3342192079702110E+02}; // km/s + // Unit adjustment + target_pos *= 0.001; // Now in kpc + target_vel *= 1.0226911647958985; // Now in kpc/Gyr + + In_box_cyl accept(2, 20, target_pos, target_vel); + Box_initial_conditions box_initial_conditions(75, 250); + + Scattering_experiment scattering_experiment("results.dat", t_max, 1000000000L, box_initial_conditions, accept); scattering_experiment.launch(); std::cout << "Bye\n";