From 842957bbba8d535aa92bb421f2c4d7b73a3d01c0 Mon Sep 17 00:00:00 2001 From: Jacob Date: Tue, 26 Oct 2021 15:00:28 +0200 Subject: [PATCH] 2d working without MBD --- src/Calculation.cpp | 50 +++++++++++++++--------- src/Calculation.h | 30 ++++++++------- src/Compute.cpp | 33 ++++++++++------ src/Compute.h | 2 +- src/Integratoren2d_force.cpp | 9 ++--- src/main.cpp | 74 +++++++++++++++++++++--------------- 6 files changed, 119 insertions(+), 79 deletions(-) diff --git a/src/Calculation.cpp b/src/Calculation.cpp index 92cdfc5..18a0f17 100644 --- a/src/Calculation.cpp +++ b/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); } } @@ -18,30 +18,23 @@ 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), sim(deltaT, seed),m_integrator(std::move(t_integrator)) { - for (const auto &pair : t_computes) { + std::function t_integrator, + std::initializer_list> t_computes, + double deltaT, size_t seed, double length) + : rod(length), 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< - void(Rod2d &, Simulation &, - std::function, - std::function)> - t_integrator, - std::initializer_list> t_computes, - double deltaT, size_t seed, - std::function force, - std::function torque, double length) - :rod(length), sim(deltaT, seed) { +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), 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) { + for (const auto &pair: t_computes) { computes.emplace_back(Compute(rod, pair.first, pair.second, sim)); } } @@ -51,3 +44,24 @@ const std::vector &Calculation::getComputes() const { } 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), 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(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), 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)); + } + rod.setPos({1000,1000}); + +} diff --git a/src/Calculation.h b/src/Calculation.h index b59ac88..3e4a366 100644 --- a/src/Calculation.h +++ b/src/Calculation.h @@ -10,29 +10,31 @@ #include "Simulation.h" class Calculation { - private: +private: Rod2d rod; Simulation sim; std::function m_integrator; std::vector computes; - public: +public: [[nodiscard]] const Simulation &getSim() const; + using force_type = std::function; using torque_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); + using inte_force_type = std::function; - explicit Calculation( - std::function t_integrator, - std::initializer_list> t_computes, - double deltaT, size_t seed, double length = 1.0); + 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(std::function t_integrator, + std::initializer_list> 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); [[nodiscard]] const std::vector &getComputes() const; diff --git a/src/Compute.cpp b/src/Compute.cpp index 1e2cafd..433f595 100644 --- a/src/Compute.cpp +++ b/src/Compute.cpp @@ -28,7 +28,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.0) / 2.0 / time; agg.feed(empxx); } @@ -37,7 +37,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.0) / 2.0 / time; agg.feed(empxx); } @@ -57,19 +57,23 @@ void Compute::eval(const Rod2d &rod2D) { case Type::empyy: eval_empYY(rod2D); break; + case Type::msd_harmonic_k_10: + evalMSD(rod2D); + break; } start_rod = rod2D; } } Compute::Compute(Rod2d rod, Type t_type, size_t t_every, Simulation &sim) - : start_rod(std::move(rod)), every(t_every), time_step(0), type(t_type) { + : start_rod(std::move(rod)), every(t_every), time_step(0), type(t_type) { time = sim.getMDeltaT() * static_cast(every); switch (type) { - case Type::msd:{ + case Type::msd: { type_str = "msd"; target = 4.0 * 0.5 * (rod.getDiff().trace()) * time; - break;} + break; + } case Type::oaf: { type_str = "oaf"; target = std::exp(-rod.getDRot() * time); @@ -81,18 +85,25 @@ Compute::Compute(Rod2d rod, Type t_type, size_t t_every, Simulation &sim) const double u = 4.0; const double deltaD = rod.getDiff()(1, 1) - rod.getDiff()(0, 0); target = Dmean - deltaD / 2.0 * - (1 - exp(-u * rod.getDRot() * time)) / u / - rod.getDRot() / time; - } break; + (1 - exp(-u * rod.getDRot() * time)) / u / + rod.getDRot() / time; + } + break; case Type::empyy: { type_str = "empyy"; 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.0 * - (1 - exp(-u * rod.getDRot() * time)) / u / - rod.getDRot() / time ; - } break; + (1 - exp(-u * rod.getDRot() * time)) / u / + rod.getDRot() / time; + } + break; + case Type::msd_harmonic_k_10: { + type_str = "msd_harmonic"; + target = 2.0/10.0*(1.0 - std::exp(-20 * time)); + } + break; } } diff --git a/src/Compute.h b/src/Compute.h index cc4cfa0..ae3213e 100644 --- a/src/Compute.h +++ b/src/Compute.h @@ -12,7 +12,7 @@ class Compute { public: - enum class Type { msd, oaf, empxx, empyy }; + enum class Type { msd, oaf, empxx, empyy, msd_harmonic_k_10 }; explicit Compute(Rod2d rod, Type t_type, size_t t_every, Simulation &sim); diff --git a/src/Integratoren2d_force.cpp b/src/Integratoren2d_force.cpp index a85e4aa..07d358b 100644 --- a/src/Integratoren2d_force.cpp +++ b/src/Integratoren2d_force.cpp @@ -32,7 +32,7 @@ void Integratoren2d_force::Set1_Euler( // Force Vector F_lab = force(rod2D.getPos(), e); Vector F_part = rotation_Matrix(e).inverse() * F_lab; - Vector F_trans = rod2D.getDiff_Sqrt() * F_part; + Vector F_trans = rotation_Matrix(e)*rod2D.getDiff_Sqrt() * F_part; F_trans *= sim.getMDeltaT(); // rotation @@ -59,7 +59,7 @@ Vector Integratoren2d_force::Heun_predictor_pos( Vector F_pred_lab = force(rod2D.getPos(), rod2D.getE()); Vector F_pred_part = rotation_Matrix(rod2D.getE()).inverse() * F_pred_lab; - Vector F_pred_trans = rod2D.getDiff_Sqrt() * F_pred_part; + Vector F_pred_trans = rotation_Matrix(rod2D.getE())*rod2D.getDiff_Sqrt() * F_pred_part; F_pred_trans *= sim.getMDeltaT() ; return F_pred_trans + trans_pred_lab; @@ -96,7 +96,6 @@ void Integratoren2d_force::Set2_Heun( Vector delta_e_future = Heun_predictor_rot(pred_rod, sim, torque); - Vector e_integrated = rod2D.getE() + 0.5 * (delta_e_predictor + delta_e_future); // apply @@ -117,7 +116,7 @@ void Integratoren2d_force::Set3_Exact( Vector F_lab = force(rod2D.getPos(), rod2D.getE()); Vector F_part = rotation_Matrix(rod2D.getE()).inverse() * F_lab; - Vector F_trans = rod2D.getDiff_Sqrt() * F_part; + Vector F_trans = rotation_Matrix(rod2D.getE())*rod2D.getDiff_Sqrt() * F_part; F_trans *= sim.getMDeltaT(); Vector trans_lab = rotation_Matrix(rod2D.getE()) * trans_part; @@ -149,7 +148,7 @@ void Integratoren2d_force::Set4_BDAS( // Force Vector F_lab = force(rod2D.getPos(), rod2D.getE()); Vector F_part = rotation_Matrix(rod2D.getE()).inverse() * F_lab; - Vector F_trans = rod2D.getDiff_Sqrt() * F_part; + Vector F_trans = rotation_Matrix(rod2D.getE())*rod2D.getDiff_Sqrt() * F_part; F_trans *= sim.getMDeltaT(); // rotation diff --git a/src/main.cpp b/src/main.cpp index 06bed95..847a6c3 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -3,6 +3,7 @@ // #include +#include #include "Calculation.h" #include "Integratoren2d_force.h" @@ -10,42 +11,55 @@ int main() { constexpr int numStep = 10000000; - constexpr double k = 0.01; + constexpr double k = 10; [[maybe_unused]] auto harmonic_Force = - [](const Eigen::Vector2d &pos, - const Eigen::Vector2d & /*torque*/) -> Eigen::Vector2d { - return -k * pos; - }; + [](const Eigen::Vector2d &pos, + const Eigen::Vector2d & /*torque*/) -> Eigen::Vector2d { + return -k * pos; + }; [[maybe_unused]] auto zero_Force = - [](const Eigen::Vector2d & /*pos*/, - const Eigen::Vector2d & /*torque*/) -> Eigen::Vector2d { - return {0.0, 0.0}; - }; + [](const Eigen::Vector2d & /*pos*/, + const Eigen::Vector2d & /*torque*/) -> Eigen::Vector2d { + return {0.0, 0.0}; + }; [[maybe_unused]] auto zero_Torque = [](const Eigen::Vector2d & /*pos*/, - const Eigen::Vector2d & /*torque*/) -> double { + const Eigen::Vector2d & /*torque*/) -> double { return 0.0; }; + std::vector> computes; + for (int i = 1; i < 100; ++i) { + computes.emplace_back(Compute::Type::msd_harmonic_k_10, i * 1); +//computes.emplace_back(Compute::Type::oaf, i * 1); + //computes.emplace_back(Compute::Type::empxx, i * 1); + //computes.emplace_back(Compute::Type::empyy, i * 1); + } + + Calculation euler(Integratoren2d_force::Set4_BDAS, + computes, 0.01, 12345, harmonic_Force, zero_Torque, 1.0); - constexpr Compute::Type compute = Compute::Type::oaf; - Calculation euler(Integratoren2d_force::Set2_Heun, - {{compute, 1}, - {compute, 2}, - {compute, 4}, - {compute, 8}, - {compute, 16}, - {compute, 32}, - {compute, 64}, - {compute, 128}, - {compute, 256}, - {compute, 512}, - {compute, 1024}}, - 0.01, 12345, zero_Force,zero_Torque,1.0); euler.run(numStep); - for (const auto &com : euler.getComputes()) { - if (com.getType() != compute) { continue; -} - std::cout - << com.getAgg().getMean() << " , " << com.getTarget() - << std::endl; + + std::ofstream msdFile("msd.dat"); + std::ofstream msd_harmFile("msd_harm.dat"); + std::ofstream oafFile("oaf.dat"); + std::ofstream empxxFile("empxx.dat"); + std::ofstream empyyFile("empyy.dat"); + + for (const auto &com: euler.getComputes()) { + if (com.getType() == Compute::Type::msd) { + msdFile << com.getTime() << ", " << com.getAgg().getMean() << ", " << com.getTarget() << std::endl; + } + if (com.getType() == Compute::Type::oaf) { + oafFile << com.getTime() << ", " << com.getAgg().getMean() << ", " << com.getTarget() << std::endl; + } + if (com.getType() == Compute::Type::empxx) { + empxxFile << com.getTime() << ", " << com.getAgg().getMean() << ", " << com.getTarget() << std::endl; + } + if (com.getType() == Compute::Type::empyy) { + empyyFile << com.getTime() << ", " << com.getAgg().getMean() << ", " << com.getTarget() << std::endl; + } + if (com.getType() == Compute::Type::msd_harmonic_k_10) { + msd_harmFile << com.getTime() << ", " << com.getAgg().getMean() << ", " << com.getTarget() << std::endl; + } } } \ No newline at end of file