From 2c16c36fdc02812e28464559197babca599c8cf7 Mon Sep 17 00:00:00 2001 From: Jacob Date: Thu, 28 Oct 2021 13:43:59 +0200 Subject: [PATCH] Minor fixes should run --- C++/src/CMakeLists.txt | 6 +++ C++/src/Calculation.cpp | 88 ++++++++++++++++++++-------------- C++/src/Calculation.h | 35 ++++++++------ C++/src/Compute.cpp | 101 ++++++++++++++++++++++++++++------------ C++/src/Compute.h | 35 ++++++++------ C++/src/Rod2d.cpp | 2 +- C++/src/main.cpp | 95 ++++++++++++++++++++----------------- 7 files changed, 228 insertions(+), 134 deletions(-) diff --git a/C++/src/CMakeLists.txt b/C++/src/CMakeLists.txt index 9249c05..9275e10 100644 --- a/C++/src/CMakeLists.txt +++ b/C++/src/CMakeLists.txt @@ -2,6 +2,8 @@ find_package(fmt) find_package(spdlog) find_package(Eigen3) + + # Generic test that uses conan libs add_library(integratoren SHARED LiveAgg.cpp Rod2d.cpp Simulation.cpp Simulation.h Integratoren2d.cpp @@ -16,6 +18,10 @@ target_link_libraries( Eigen3::Eigen3) add_executable(main main.cpp) +find_package(OpenMP) +if(OpenMP_CXX_FOUND) + target_link_libraries(main PUBLIC OpenMP::OpenMP_CXX) +endif() target_link_libraries( main PRIVATE project_options diff --git a/C++/src/Calculation.cpp b/C++/src/Calculation.cpp index cdc2c5d..b858535 100644 --- a/C++/src/Calculation.cpp +++ b/C++/src/Calculation.cpp @@ -9,7 +9,7 @@ void Calculation::run(size_t steps) { for (size_t step = 0; step < steps; ++step) { m_integrator(rod, sim); - for (auto &comp: computes) { + for (auto &comp : computes) { comp.eval(rod); } } @@ -17,56 +17,74 @@ void Calculation::run(size_t steps) { const Rod2d &Calculation::getRod() const { return rod; } -Calculation::Calculation( - std::function t_integrator, - std::initializer_list> t_computes, - double deltaT, size_t seed, double length) - : rod(length), start_rod(rod), sim(deltaT, seed), m_integrator(std::move(t_integrator)) { - for (const auto &pair: t_computes) { - computes.emplace_back(Compute(rod, pair.first, pair.second, sim)); - } -} - -Calculation::Calculation(inte_force_type t_integrator, - std::initializer_list> t_computes, - double deltaT, size_t seed, force_type force, torque_type torque, double length) - : rod(length), start_rod(rod), sim(deltaT, seed) { - m_integrator = [&](Rod2d &t_rod, Simulation &t_sim) { - t_integrator(t_rod, t_sim, force, torque); - }; - for (const auto &pair: t_computes) { - computes.emplace_back(Compute(rod, pair.first, pair.second, sim)); - } -} - const std::vector &Calculation::getComputes() const { return computes; } const Simulation &Calculation::getSim() const { return sim; } -Calculation::Calculation(std::function t_integrator, const std::vector>& t_computes, - double deltaT, size_t seed, double length) : rod(length), start_rod(rod), sim(deltaT, seed), - m_integrator(std::move(t_integrator)){ - for (const auto &pair: t_computes) { - computes.emplace_back(Compute(rod, pair.first, pair.second, sim)); +Calculation::Calculation( + std::function t_integrator, + std::initializer_list> t_computes, + double deltaT, size_t seed, double length) + : rod(length), + start_rod(rod), + sim(deltaT, seed), + m_integrator(std::move(t_integrator)) { + for (const auto &pair : t_computes) { + computes.emplace_back( + Compute(rod, pair.first, Compute::Force::zero_F, pair.second, sim)); } } -Calculation::Calculation(Calculation::inte_force_type t_integrator, std::vector> t_computes, double deltaT, - size_t seed, Calculation::force_type force, Calculation::torque_type torque, double length) - : rod(length), start_rod(rod), sim(deltaT, seed) { +Calculation::Calculation( + inte_force_type t_integrator, + std::initializer_list> t_computes, + double deltaT, size_t seed, force_type force, Compute::Force type_force, + torque_type torque, double length) + : rod(length), start_rod(rod), sim(deltaT, seed) { m_integrator = [&](Rod2d &t_rod, Simulation &t_sim) { t_integrator(t_rod, t_sim, force, torque); }; - for (const auto &pair: t_computes) { - computes.emplace_back(Compute(rod, pair.first, pair.second, sim)); + for (const auto &pair : t_computes) { + computes.emplace_back( + Compute(rod, pair.first, type_force, pair.second, sim)); + } +} + + +Calculation::Calculation( + std::function t_integrator, + const std::vector> &t_computes, + double deltaT, size_t seed, double length) + : rod(length), + start_rod(rod), + sim(deltaT, seed), + m_integrator(std::move(t_integrator)) { + for (const auto &pair : t_computes) { + computes.emplace_back( + Compute(rod, pair.first, Compute::Force::zero_F, pair.second, sim)); + } +} + +Calculation::Calculation( + Calculation::inte_force_type t_integrator, + const std::vector>& t_computes, double deltaT, + size_t seed, Calculation::force_type force, Compute::Force type_force, + Calculation::torque_type torque, double length) + : rod(length), start_rod(rod), sim(deltaT, seed) { + m_integrator = [&](Rod2d &t_rod, Simulation &t_sim) { + t_integrator(t_rod, t_sim, force, torque); + }; + for (const auto &pair : t_computes) { + computes.emplace_back( + Compute(rod, pair.first, type_force, pair.second, sim)); } } void Calculation::reset() { - rod =start_rod; - for(auto & c:computes){ + rod = start_rod; + for (auto &c : computes) { c.reset(start_rod); } } diff --git a/C++/src/Calculation.h b/C++/src/Calculation.h index 973a3ea..a745b6d 100644 --- a/C++/src/Calculation.h +++ b/C++/src/Calculation.h @@ -10,35 +10,44 @@ #include "Simulation.h" class Calculation { -private: + private: Rod2d rod; const Rod2d start_rod; Simulation sim; std::function m_integrator; std::vector computes; -public: - + public: void reset(); [[nodiscard]] const Simulation &getSim() const; - using force_type = std::function; + using force_type = + std::function; using torque_type = std::function; - using inte_force_type = std::function; + using inte_force_type = + std::function; - Calculation(inte_force_type t_integrator, std::initializer_list> t_computes, - double deltaT, size_t seed, force_type force, torque_type torque, double length = 1.0); + Calculation( + inte_force_type t_integrator, + std::initializer_list> t_computes, + double deltaT, size_t seed, force_type force, Compute::Force force_type, + torque_type torque, double length = 1.0); + + Calculation( + std::function t_integrator, + std::initializer_list> t_computes, + double deltaT, size_t seed, double length = 1.0); Calculation(std::function t_integrator, - std::initializer_list> t_computes, + const std::vector> &t_computes, double deltaT, size_t seed, double length = 1.0); - Calculation(std::function t_integrator, - const std::vector>& t_computes, double deltaT, size_t seed, double length = 1.0); - - Calculation(inte_force_type t_integrator, std::vector> t_computes, - double deltaT, size_t seed, force_type force, torque_type torque, double length = 1.0); + Calculation(inte_force_type t_integrator, + const std::vector>& t_computes, + double deltaT, size_t seed, force_type force, + Compute::Force force_type, torque_type torque, + double length = 1.0); [[nodiscard]] const std::vector &getComputes() const; diff --git a/C++/src/Compute.cpp b/C++/src/Compute.cpp index 5cf59b4..68eab69 100644 --- a/C++/src/Compute.cpp +++ b/C++/src/Compute.cpp @@ -19,8 +19,13 @@ void Compute::evalMSD(const Rod2d &rod2D) { void Compute::evalX(const Rod2d &rod2D) { const auto &new_pos = rod2D.getPos(); - auto old_pos = start_rod.getPos(); - auto msd = (new_pos - old_pos).norm(); + auto msd = (new_pos).sum(); + agg.feed(msd); +} + +void Compute::evalX_squared(const Rod2d &rod2D) { + const auto &new_pos = rod2D.getPos(); + auto msd = (new_pos).squaredNorm(); agg.feed(msd); } @@ -35,7 +40,7 @@ void Compute::eval_empXX(const Rod2d &rod2D) { const auto &new_pos = rod2D.getPos(); auto old_pos = start_rod.getPos(); auto old_e = start_rod.getE(); - double empxx = pow((new_pos - old_pos).dot(old_e), 2.0) / 2.0 / time; + double empxx = pow((new_pos - old_pos).dot(old_e), 2) / 2 / time; agg.feed(empxx); } @@ -44,7 +49,7 @@ void Compute::eval_empYY(const Rod2d &rod2D) { auto old_pos = start_rod.getPos(); auto old_e = start_rod.getE(); Eigen::Vector2d old_e_ortho = {-old_e[1], old_e[0]}; - double empxx = pow((new_pos - old_pos).dot(old_e_ortho), 2.0) / 2.0 / time; + double empxx = pow((new_pos - old_pos).dot(old_e_ortho), 2) / 2 / time; agg.feed(empxx); } @@ -64,11 +69,6 @@ void Compute::eval(const Rod2d &rod2D) { case Type::empyy: eval_empYY(rod2D); break; - case Type::msd_harmonic_k_1: - evalMSD(rod2D); - break; - case Type::msd_const: - break; case Type::x: if (time_step == every) { evalX(rod2D); @@ -76,7 +76,7 @@ void Compute::eval(const Rod2d &rod2D) { break; case Type::x_squared: if (time_step == every) { - evalMSD(rod2D); + evalX_squared(rod2D); } break; } @@ -87,7 +87,8 @@ void Compute::eval(const Rod2d &rod2D) { } } -Compute::Compute(Rod2d rod, Type t_type, size_t t_every, Simulation &sim) +Compute::Compute(Rod2d rod, Type t_type, Force force, size_t t_every, + Simulation &sim) : start_rod(std::move(rod)), every(t_every), time_step(0), @@ -97,9 +98,20 @@ Compute::Compute(Rod2d rod, Type t_type, size_t t_every, Simulation &sim) case Type::msd: { resetting = true; type_str = "msd"; - target = 2 * (rod.getDiff().trace()) * time; + switch (force) { + case Force::zero_F: + target = 2 * (rod.getDiff().trace()) * time; + break; + case Force::const_F: + target = static_cast(NAN); + break; + case Force::harmonic_F: + target = 2 * (1.0 - std::exp(-2 * time)); + break; + } break; } + case Type::oaf: { resetting = true; type_str = "oaf"; @@ -112,8 +124,17 @@ Compute::Compute(Rod2d rod, Type t_type, size_t t_every, Simulation &sim) const double Dmean = 0.5 * (rod.getDiff().trace()); const double u = 4.0; const double deltaD = rod.getDiff()(1, 1) - rod.getDiff()(0, 0); - target = Dmean - deltaD / 2 * (1 - exp(-u * rod.getDRot() * time)) / - u / rod.getDRot() / time; + switch (force) { + case Force::zero_F: + target = Dmean - deltaD / 2 * + (1 - exp(-u * rod.getDRot() * time)) / + u / rod.getDRot() / time; + break; + case Force::const_F: + case Force::harmonic_F: + target = static_cast(NAN); + break; + } } break; case Type::empyy: { resetting = true; @@ -121,29 +142,51 @@ Compute::Compute(Rod2d rod, Type t_type, size_t t_every, Simulation &sim) const double Dmean = 0.5 * (rod.getDiff().trace()); const double u = 4.0; const double deltaD = rod.getDiff()(1, 1) - rod.getDiff()(0, 0); - target = Dmean + deltaD / 2 * (1 - exp(-u * rod.getDRot() * time)) / - u / rod.getDRot() / time; - } break; - case Type::msd_harmonic_k_1: { - resetting = true; - type_str = "msd_harmonic"; - target = 2 * (1.0 - std::exp(-2 * time)); - } break; - case Type::msd_const: { - resetting = true; - type_str = "msd_const"; - target = -1.0; + switch (force) { + case Force::zero_F: + target = Dmean + deltaD / 2 * (1 - exp(-u * rod.getDRot() * time)) / + u / rod.getDRot() / time; + break; + case Force::const_F: + case Force::harmonic_F: + target = static_cast(NAN); + break; + } } break; case Type::x: { resetting = false; type_str = ""; - target = rod.getPos().norm()*exp(-1.0*time); + switch (force) { + case Force::zero_F: + target = rod.getPos().sum(); + break; + case Force::const_F: + target = static_cast(NAN); + break; + case Force::harmonic_F: + target = rod.getPos().sum() * exp(-1 * time); + break; + } + + + } break; case Type::x_squared: { resetting = false; type_str = ""; - target = rod.getPos().squaredNorm()*exp(-2.0*time)+ - 2.0*(1-exp(-2.0*time)); + + switch (force) { + case Force::zero_F: + target = static_cast(NAN); + break; + case Force::const_F: + target = static_cast(NAN); + break; + case Force::harmonic_F: + target = rod.getPos().squaredNorm() * exp(-2 * time) + + 2 * (1 - exp(-2 * time)); + break; + } } break; } } diff --git a/C++/src/Compute.h b/C++/src/Compute.h index 309b336..a92b55b 100644 --- a/C++/src/Compute.h +++ b/C++/src/Compute.h @@ -12,22 +12,26 @@ class Compute { public: - enum class Type { msd, oaf, empxx, empyy, msd_harmonic_k_1,msd_const, x,x_squared }; + enum class Type { + msd, + oaf, + empxx, + empyy, + x, + x_squared + }; + enum class Force { + zero_F, + const_F, + harmonic_F + }; - explicit Compute(Rod2d rod, Type t_type, size_t t_every, Simulation &sim); + explicit Compute(Rod2d rod, Type t_type, Force force,size_t t_every, Simulation &sim); void eval(const Rod2d &rod2D); - void evalMSD(const Rod2d &rod2D); - - void evalOAF(const Rod2d &rod2D); - - void eval_empXX(const Rod2d &rod2D); - [[nodiscard]] double getTarget() const; - void eval_empYY(const Rod2d &rod2D); - [[nodiscard]] const LiveAgg &getAgg() const; [[nodiscard]] Type getType() const; @@ -40,7 +44,7 @@ class Compute { void reset(const Rod2d &rod); -private: + private: Rod2d start_rod; LiveAgg agg; size_t every; @@ -49,11 +53,14 @@ private: double target; double time; std::string type_str; - + bool resetting; void evalX(const Rod2d &rod2D); - - bool resetting; + void evalMSD(const Rod2d &rod2D); + void evalOAF(const Rod2d &rod2D); + void eval_empXX(const Rod2d &rod2D); + void eval_empYY(const Rod2d &rod2D); + void evalX_squared(const Rod2d &rod2D); }; #endif // MYPROJECT_COMPUTE_H diff --git a/C++/src/Rod2d.cpp b/C++/src/Rod2d.cpp index 62e0709..f306612 100644 --- a/C++/src/Rod2d.cpp +++ b/C++/src/Rod2d.cpp @@ -10,7 +10,7 @@ void sqrt(Eigen::Matrix2d &mat) { mat.data()[i] = sqrt(mat.data()[i]); } } -Rod2d::Rod2d(double L) : m_pos({0, 0}), m_e({1, 0}) { +Rod2d::Rod2d(double L) : m_pos({1, 1}), m_e({1, 0}) { assert(L >= 1.0); if (L == 1.0) { m_D_rot = 3.0; diff --git a/C++/src/main.cpp b/C++/src/main.cpp index 837d22c..f4c7032 100644 --- a/C++/src/main.cpp +++ b/C++/src/main.cpp @@ -11,12 +11,14 @@ #include "force_lambdas.h" constexpr size_t SEED = 1234; constexpr double stepSize = 0.001; -constexpr int n_computes = 100; -constexpr size_t num_runs = 5; -constexpr int numStep = 1000000; -constexpr int delta_compute = 4; -void run_no_force(size_t integrator_index, size_t force_index, - size_t length_index) { +constexpr size_t n_computes = 100; +constexpr size_t num_integratoren = 5; +constexpr size_t num_forces = 3; +constexpr size_t num_length = 3; +constexpr size_t delta_compute = 4; +constexpr size_t numStep = 10000; +inline const std::string header = {"time,val,target,std\n"}; +void run(size_t integrator_index, size_t force_index, size_t length_index) { using type = std::pair; std::vector vec({type(Integratoren2d::Set1_Euler_f, "euler"), type(Integratoren2d::Set2_Heun_f, "heun"), @@ -30,90 +32,100 @@ void run_no_force(size_t integrator_index, size_t force_index, {zero_Force, const_Force, harmonic_Force})[force_index]; auto force_name = std::vector( {"zero_force_", "const_force_", "harmonic_force_"})[force_index]; + auto force_type = std::vector( + {Compute::Force::zero_F, Compute::Force::const_F, + Compute::Force::harmonic_F})[force_index]; namespace fs = std::filesystem; - std::string folder = "out/" + force_name + name + "_L"+std::to_string(length_index)+"_"; + std::string folder = + "out/" + force_name + name + "_L" + std::to_string(length_index) + "_"; fs::create_directories("out"); double length = std::vector({1.0, 1.5, 2.0})[length_index]; { - Compute::Type MSD_VARIANT; - if (force_index == 0) MSD_VARIANT = Compute::Type::msd; - if (force_index == 1) MSD_VARIANT = Compute::Type::msd_const; - if (force_index == 2) MSD_VARIANT = Compute::Type::msd_harmonic_k_1; - std::vector> computes; - for (int i = 1; i < n_computes; ++i) { - computes.emplace_back(MSD_VARIANT, i * delta_compute); + for (size_t i = 1; i < n_computes; ++i) { + computes.emplace_back(Compute::Type::msd, i * delta_compute); computes.emplace_back(Compute::Type::oaf, i * delta_compute); computes.emplace_back(Compute::Type::empxx, i * delta_compute); computes.emplace_back(Compute::Type::empyy, i * delta_compute); } Calculation euler(integrator, computes, stepSize, SEED, force, - zero_Torque, length); + force_type, zero_Torque, length); - euler.run(numStep); + for (size_t i = 0; i < numStep; ++i) { + euler.run(n_computes * delta_compute); + } std::ofstream msdFile(folder + "msd.dat"); std::ofstream oafFile(folder + "oaf.dat"); std::ofstream empxxFile(folder + "empxx.dat"); std::ofstream empyyFile(folder + "empyy.dat"); - msdFile << "time,val,target\n"; - oafFile << "time,val,target\n"; - empxxFile << "time,val,target\n"; - empyyFile << "time,val,target\n"; + msdFile << header; + oafFile << header; + empxxFile << header; + empyyFile << header; for (const auto &com : euler.getComputes()) { - if (com.getType() == MSD_VARIANT) { + if (com.getType() == Compute::Type::msd) { msdFile << com.getTime() << ", " << com.getAgg().getMean() - << ", " << com.getTarget() << std::endl; + << ", " << com.getTarget() << ", " + << com.getAgg().getSEM() << std::endl; } if (com.getType() == Compute::Type::oaf) { oafFile << com.getTime() << ", " << com.getAgg().getMean() - << ", " << com.getTarget() << std::endl; + << ", " << com.getTarget() << ", " + << com.getAgg().getSEM() << std::endl; } if (com.getType() == Compute::Type::empxx) { empxxFile << com.getTime() << ", " << com.getAgg().getMean() - << ", " << com.getTarget() << std::endl; + << ", " << com.getTarget() << ", " + << com.getAgg().getSEM() << std::endl; } if (com.getType() == Compute::Type::empyy) { empyyFile << com.getTime() << ", " << com.getAgg().getMean() - << ", " << com.getTarget() << std::endl; + << ", " << com.getTarget() << ", " + << com.getAgg().getSEM() << std::endl; } if (com.getType() == Compute::Type::msd) { msdFile << com.getTime() << ", " << com.getAgg().getMean() - << ", " << com.getTarget() << std::endl; + << ", " << com.getTarget() << ", " + << com.getAgg().getSEM() << std::endl; } } } { std::vector> repeating_computes; - for (int i = 1; i < n_computes; ++i) { - repeating_computes.emplace_back(Compute::Type::x, i * delta_compute); - repeating_computes.emplace_back(Compute::Type::x_squared, i * delta_compute); + for (size_t i = 1; i < n_computes; ++i) { + repeating_computes.emplace_back(Compute::Type::x, + i * delta_compute); + repeating_computes.emplace_back(Compute::Type::x_squared, + i * delta_compute); } Calculation calc_repeat(integrator, repeating_computes, stepSize, SEED, - force, zero_Torque, 1.0); - for (int i = 0; i < numStep / n_computes; ++i) { - calc_repeat.run(n_computes); + force, force_type, zero_Torque, 1.0); + for (size_t i = 0; i < numStep; ++i) { + calc_repeat.run(n_computes * delta_compute); calc_repeat.reset(); } std::ofstream xFile(folder + "x.dat"); std::ofstream xsqFile(folder + "x_squared.dat"); - xFile << "time,val,target\n"; - xsqFile << "time,val,target\n"; + xFile << header; + xsqFile << header; for (const auto &com : calc_repeat.getComputes()) { if (com.getType() == Compute::Type::x) { xFile << com.getTime() << ", " << com.getAgg().getMean() << ", " - << com.getTarget() << std::endl; + << com.getTarget() << ", " << com.getAgg().getSEM() + << std::endl; } if (com.getType() == Compute::Type::x_squared) { xsqFile << com.getTime() << ", " << com.getAgg().getMean() - << ", " << com.getTarget() << std::endl; + << ", " << com.getTarget() << ", " + << com.getAgg().getSEM() << std::endl; } } } @@ -122,11 +134,10 @@ void run_no_force(size_t integrator_index, size_t force_index, } int main() { - for (size_t j = 0; j < 3; ++j) { - for (size_t i = 0; i < num_runs; ++i) { - for (size_t k = 0; k < 1; ++k) { - run_no_force(i, j, k); - } - } +#pragma omp parallel for default(none) + for (size_t j = 0; j < num_forces * num_length * num_integratoren; ++j) { + run(j % (num_integratoren), + (j % (num_integratoren * num_forces)) / num_integratoren, + j / (num_integratoren * num_forces)); } } \ No newline at end of file