From 8af7eee26783617b510747de3762a92c011d1a54 Mon Sep 17 00:00:00 2001 From: Jacob Date: Mon, 25 Oct 2021 18:45:54 +0200 Subject: [PATCH] Multiple cahnges :D --- benchmark/bench.cpp | 3 + src/Calculation.cpp | 8 +- src/Calculation.h | 19 ++- src/Compute.cpp | 26 +++- src/Compute.h | 7 +- src/Integratoren2d_force.cpp | 137 +++++++++-------- src/Integratoren2d_force.h | 16 +- src/Integratoren2d_forceless.cpp | 139 +++++------------- src/Integratoren2d_forceless.h | 5 - src/LiveAgg.cpp | 5 + src/LiveAgg.hpp | 7 +- src/Rod2d.cpp | 5 - src/Rod2d.hpp | 2 - src/legacy/Integratoren.cpp | 130 ----------------- src/legacy/Integratoren.hpp | 19 --- src/legacy/Simulation2d.cpp | 6 - src/legacy/Simulation2d.hpp | 157 -------------------- src/legacy/main.cpp | 72 --------- src/main.cpp | 66 +++------ src/vec_trafos.h | 13 ++ test/test_Integratoren.cpp | 243 +++++++++++++++++++++++++++++-- test/test_Rod2d.cpp | 40 +---- 22 files changed, 440 insertions(+), 685 deletions(-) delete mode 100644 src/legacy/Integratoren.cpp delete mode 100644 src/legacy/Integratoren.hpp delete mode 100644 src/legacy/Simulation2d.cpp delete mode 100644 src/legacy/Simulation2d.hpp delete mode 100644 src/legacy/main.cpp create mode 100644 src/vec_trafos.h diff --git a/benchmark/bench.cpp b/benchmark/bench.cpp index 528d1ea..f351722 100644 --- a/benchmark/bench.cpp +++ b/benchmark/bench.cpp @@ -24,4 +24,7 @@ TEST_CASE("Euler - Baseline", "[benchmark]") { BENCHMARK("BDAS without force") { Integratoren2d_forceless::Set4_BDAS(rod, sim); }; + BENCHMARK("MBD without force") { + Integratoren2d_forceless::Set5_MBD(rod, sim); + }; } \ No newline at end of file diff --git a/src/Calculation.cpp b/src/Calculation.cpp index 9eaf12f..92cdfc5 100644 --- a/src/Calculation.cpp +++ b/src/Calculation.cpp @@ -20,8 +20,8 @@ const Rod2d &Calculation::getRod() const { return rod; } Calculation::Calculation( std::function t_integrator, std::initializer_list> t_computes, - double deltaT, size_t seed) - : sim(deltaT, seed), m_integrator(std::move(t_integrator)) { + 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)); } @@ -36,8 +36,8 @@ Calculation::Calculation( std::initializer_list> t_computes, double deltaT, size_t seed, std::function force, - std::function torque) - : sim(deltaT, seed) { + std::function torque, double length) + :rod(length), sim(deltaT, seed) { m_integrator = [&](Rod2d &t_rod, Simulation &t_sim) { t_integrator(t_rod, t_sim, force, torque); }; diff --git a/src/Calculation.h b/src/Calculation.h index 0078f5b..b59ac88 100644 --- a/src/Calculation.h +++ b/src/Calculation.h @@ -11,29 +11,28 @@ class Calculation { private: - Rod2d rod = Rod2d(1.0); + Rod2d rod; Simulation sim; std::function m_integrator; std::vector computes; public: - const Simulation &getSim() const; - + [[nodiscard]] const Simulation &getSim() const; + using force_type = std::function; + using torque_type = std::function; + using inte_force_type = std::function; Calculation( - std::function, - std::function)> + inte_force_type t_integrator, std::initializer_list> t_computes, double deltaT, size_t seed, - std::function force, - std::function torque); + force_type force, + torque_type torque, double length = 1.0); explicit Calculation( std::function t_integrator, std::initializer_list> t_computes, - double deltaT, size_t seed); + double deltaT, size_t seed, double length = 1.0); [[nodiscard]] const std::vector &getComputes() const; diff --git a/src/Compute.cpp b/src/Compute.cpp index 8db6a05..1e2cafd 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); + 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); + double empxx = pow((new_pos - old_pos).dot(old_e_ortho), 2.0)/2.0 / time; agg.feed(empxx); } @@ -64,15 +64,19 @@ void Compute::eval(const Rod2d &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) { - auto time = sim.getMDeltaT() * static_cast(every); + 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; - case Type::oaf: + break;} + case Type::oaf: { + type_str = "oaf"; target = std::exp(-rod.getDRot() * time); break; + } case Type::empxx: { + type_str = "empxx"; const double Dmean = 0.5 * (rod.getDiff().trace()); const double u = 4.0; const double deltaD = rod.getDiff()(1, 1) - rod.getDiff()(0, 0); @@ -81,12 +85,13 @@ Compute::Compute(Rod2d rod, Type t_type, size_t t_every, Simulation &sim) 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; + rod.getDRot() / time ; } break; } } @@ -99,4 +104,9 @@ double Compute::getTime() const { return static_cast(every); } double Compute::getTarget() const { return target; } -double Compute::getDifference() const { return abs(agg.getMean() - target); } \ No newline at end of file +double Compute::getDifference() const { return abs(agg.getMean() - target); } + +std::ostream &operator<<(std::ostream &os, const Compute &compute) { + os << "agg: " << compute.agg << " type: " << compute.type_str; + return os; +} diff --git a/src/Compute.h b/src/Compute.h index e121171..cc4cfa0 100644 --- a/src/Compute.h +++ b/src/Compute.h @@ -5,6 +5,7 @@ #ifndef MYPROJECT_COMPUTE_H #define MYPROJECT_COMPUTE_H +#include #include "LiveAgg.hpp" #include "Rod2d.hpp" #include "Simulation.h" @@ -35,13 +36,17 @@ class Compute { double getDifference() const; - private: + friend std::ostream &operator<<(std::ostream &os, const Compute &compute); + +private: Rod2d start_rod; LiveAgg agg; size_t every; size_t time_step; Type type; double target; + double time; + std::string type_str; }; #endif // MYPROJECT_COMPUTE_H diff --git a/src/Integratoren2d_force.cpp b/src/Integratoren2d_force.cpp index 85e5687..9bcd80e 100644 --- a/src/Integratoren2d_force.cpp +++ b/src/Integratoren2d_force.cpp @@ -3,17 +3,10 @@ // #include "Integratoren2d_force.h" +#include "vec_trafos.h" -constexpr double kBT = 1.0; using Vector = Eigen::Vector2d; - -Vector Integratoren2d_force::ortho(Vector e) { return Vector{-e[1], e[0]}; } - -Eigen::Matrix2d Integratoren2d_force::MatTraf(Vector e) { - Eigen::Matrix2d mat; - mat << e[0], -e[1], e[1], e[0]; - return mat; -} +using Matrix = Eigen::Matrix2d; void Integratoren2d_force::Set0_deterministic(Rod2d &rod2D, Simulation & /*sim*/) { @@ -25,31 +18,27 @@ void Integratoren2d_force::Set0_deterministic(Rod2d &rod2D, const Vector Integratoren2d_force::unitVec = {1, 0}; -void Integratoren2d_force::Set5_MBD(Rod2d & /*rod2D*/, Simulation & /*sim*/) {} - void Integratoren2d_force::Set1_Euler( - Rod2d &rod2D, Simulation &sim, - const std::function &force, - const std::function &torque) { + Rod2d &rod2D, Simulation &sim, + const std::function &force, + const std::function &torque) { auto std = sim.getSTD(); // sqrt(2*delta_T) // translation Vector rand = {sim.getNorm(std), sim.getNorm(std)}; // gaussian noise - Vector trans_part = - rod2D.getDiff_Sqrt() * rand; // translation vector in Particle System - Vector trans_lab = rod2D.getE_Base_matrix() * trans_part; + Vector trans_part = rod2D.getDiff_Sqrt() * rand; + Vector trans_lab = rotation_Matrix(rod2D.getE()) * trans_part; // Force Vector F_lab = force(rod2D.getPos(), rod2D.getE()); - Vector F_part = rod2D.getE_Base_matrix().inverse() * F_lab; + Vector F_part = rotation_Matrix(rod2D.getE()).inverse() * F_lab; Vector F_trans = rod2D.getDiff_Sqrt() * F_part; - F_trans *= sim.getMDeltaT() / kBT; + F_trans *= sim.getMDeltaT(); // rotation double rot = sim.getNorm(std) * rod2D.getDRot_Sqrt(); // rotationsmatrix verwenden. Vector e = rod2D.getE(); - double M_rot = torque(rod2D.getPos(), rod2D.getE()) * rod2D.getDRot() / - kBT * sim.getMDeltaT(); + double M_rot = torque(rod2D.getPos(), rod2D.getE()) * rod2D.getDRot() * sim.getMDeltaT(); Vector new_e = e + (rot + M_rot) * ortho(e); new_e.normalize(); // apply @@ -66,12 +55,12 @@ Vector Integratoren2d_force::Heun_predictor_pos( Vector trans_pred_part = rod2D.getDiff_Sqrt() * rand_pred; // translation vector in Particle System - Vector trans_pred_lab = rod2D.getE_Base_matrix() * trans_pred_part; + Vector trans_pred_lab = rotation_Matrix(rod2D.getE()) * trans_pred_part; Vector F_pred_lab = force(rod2D.getPos(), rod2D.getE()); - Vector F_pred_part = rod2D.getE_Base_matrix().inverse() * F_pred_lab; + Vector F_pred_part = rotation_Matrix(rod2D.getE()).inverse() * F_pred_lab; Vector F_pred_trans = rod2D.getDiff_Sqrt() * F_pred_part; - F_pred_trans *= sim.getMDeltaT() / kBT; + F_pred_trans *= sim.getMDeltaT() ; return F_pred_trans + trans_pred_lab; } @@ -82,9 +71,9 @@ Vector Integratoren2d_force::Heun_predictor_rot( auto std = sim.getSTD(); double rot_predict = sim.getNorm(std) * rod2D.getDRot_Sqrt(); // rotationsmatrix verwenden. - Vector e = rod2D.getE(); + const Vector& e = rod2D.getE(); double M_predict_rot = torque(rod2D.getPos(), rod2D.getE()) * - rod2D.getDRot() / kBT * sim.getMDeltaT(); + rod2D.getDRot() * sim.getMDeltaT(); Vector e_change_predict = (rot_predict + M_predict_rot) * ortho(e); return e_change_predict; @@ -98,30 +87,27 @@ void Integratoren2d_force::Set2_Heun( Vector pos_predictor = rod2D.getPos() + delta_pos_predictor; Vector delta_e_predictor = Heun_predictor_rot(rod2D, sim, torque); - Vector e_predict = rod2D.getE() + delta_e_predictor; - e_predict.normalize(); + Vector e_predictor = rod2D.getE() + delta_e_predictor; + e_predictor.normalize(); Rod2d pred_rod = rod2D; pred_rod.setPos(pos_predictor); - pred_rod.setE(e_predict); + pred_rod.setE(e_predictor); - Vector delta_pos_future = Heun_predictor_pos(pred_rod, sim, force); Vector delta_e_future = Heun_predictor_rot(pred_rod, sim, torque); - // integration - Vector pos_integrated = - 0.5 * rod2D.getPos() + 0.5 * pos_predictor + 0.5 * delta_pos_future; + Vector e_integrated = rod2D.getE() + 0.5 * (delta_e_predictor + delta_e_future); // apply - rod2D.setPos(pos_integrated); + rod2D.setPos(pos_predictor); rod2D.setE(e_integrated.normalized()); } void Integratoren2d_force::Set3_Exact( Rod2d &rod2D, Simulation &sim, - std::function force, - std::function torque) { + const std::function& force, + const std::function& torque) { auto std = sim.getSTD(); // sqrt(2*delta_T) // translation Vector rand = {sim.getNorm(std), sim.getNorm(std)}; // gaussian noise @@ -129,20 +115,19 @@ void Integratoren2d_force::Set3_Exact( rod2D.getDiff_Sqrt() * rand; // translation vector in Particle System Vector F_lab = force(rod2D.getPos(), rod2D.getE()); - Vector F_part = rod2D.getE_Base_matrix().inverse() * F_lab; + Vector F_part = rotation_Matrix(rod2D.getE()).inverse() * F_lab; Vector F_trans = rod2D.getDiff_Sqrt() * F_part; - F_trans *= sim.getMDeltaT() / kBT; - Vector trans_lab = rod2D.getE_Base_matrix() * trans_part; + F_trans *= sim.getMDeltaT(); + Vector trans_lab = rotation_Matrix(rod2D.getE()) * trans_part; // rotation double rot = sim.getNorm(std) * rod2D.getDRot_Sqrt(); // rotationsmatrix verwenden. Vector e = rod2D.getE(); - double M_rot = torque(rod2D.getPos(), rod2D.getE()) * rod2D.getDRot() / - kBT * sim.getMDeltaT(); + double M_rot = torque(rod2D.getPos(), rod2D.getE()) * rod2D.getDRot() * sim.getMDeltaT(); auto correction = - -0.5 * pow(rod2D.getDRot(), 2) / pow(kBT, 2) * sim.getMDeltaT(); + -0.5 * pow(rod2D.getDRot(), 2) * sim.getMDeltaT(); Vector new_e = e + (rot + M_rot) * ortho(e) + correction * e; new_e.normalize(); // apply @@ -152,30 +137,58 @@ void Integratoren2d_force::Set3_Exact( void Integratoren2d_force::Set4_BDAS( Rod2d &rod2D, Simulation &sim, - std::function force, - std::function /*torque*/) { - auto std = sim.getSTD(); + const std::function& force, + const std::function& torque) { + auto std = sim.getSTD(); // sqrt(2*delta_T) + Vector e = rod2D.getE(); // translation - auto rand = Vector(sim.getNorm(std), sim.getNorm(std)); - auto trans_part = - rod2D.getDiff_Sqrt() * rand; // translation vector in Particle System - Eigen::Rotation2D rotMat(acos(unitVec.dot(rod2D.getE()))); - auto trans_lab = rotMat * trans_part; + Vector rand = {sim.getNorm(std), sim.getNorm(std)}; // gaussian noise + Vector trans_part = rod2D.getDiff_Sqrt() * rand; + Vector trans_lab = rotation_Matrix(rod2D.getE()) * trans_part; + // Force Vector F_lab = force(rod2D.getPos(), rod2D.getE()); - Vector F_part = rotMat.inverse() * F_lab; - + Vector F_part = rotation_Matrix(rod2D.getE()).inverse() * F_lab; Vector F_trans = rod2D.getDiff_Sqrt() * F_part; - F_trans *= sim.getMDeltaT() / kBT; + F_trans *= sim.getMDeltaT(); - auto rot = - Eigen::Rotation2D(sim.getNorm(std) * rod2D.getDRot_Sqrt()); - Eigen::Rotation2Dd rotation(rot); - auto e_new = (rotation.toRotationMatrix() * rod2D.getE()).normalized(); + // rotation + double rot = + sim.getNorm(std) * rod2D.getDRot_Sqrt(); // rotationsmatrix verwenden. - // Normalisation should not be necessary if a proper angular representation - // is used. But with vector e it is done in case of numerical errors - - rod2D.setPos(rod2D.getPos() + trans_lab); - rod2D.setE(e_new); + double M_rot = torque(rod2D.getPos(), rod2D.getE()) * rod2D.getDRot() * sim.getMDeltaT(); + Vector new_e = Eigen::Rotation2D(rot + M_rot) * e; + new_e.normalize(); + // apply + rod2D.setPos(rod2D.getPos() + trans_lab + F_trans); + rod2D.setE(new_e); +} + +void Integratoren2d_force::Set5_MBD(Rod2d &rod2D, Simulation &sim, + const std::function& force, + const std::function& torque) { + + Vector delta_pos_predictor = Heun_predictor_pos(rod2D, sim, force); + [[maybe_unused]] Vector pos_predictor = rod2D.getPos() + delta_pos_predictor; + + Vector delta_e_predictor = Heun_predictor_rot(rod2D, sim, torque); + Vector e_predictor = rod2D.getE() + delta_e_predictor; + e_predictor.normalize(); + + + auto std = sim.getSTD(); + Eigen::Vector2d e = rod2D.getE(); + + + Eigen::Matrix2d Diff_combined = 0.5*(rod2D.getDiff() + rotation_Matrix(e).inverse() * rotation_Matrix(e_predictor)*rod2D.getDiff()); + Eigen::Matrix2d D_combined_sqrt = Diff_combined.llt().matrixL(); + + auto rand = Eigen::Vector2d(sim.getNorm(std), sim.getNorm(std)); + Eigen::Vector2d pos_integrated =rod2D.getPos() + rotation_Matrix(e) * D_combined_sqrt * rand; + + double rot = sim.getNorm(std) * rod2D.getDRot_Sqrt(); + Eigen::Vector2d e_integrated = e + 0.5 * (rot * ortho(e) + rot * ortho(e_predictor)); + //apply + rod2D.setPos(pos_integrated); + rod2D.setE(e_integrated.normalized()); } diff --git a/src/Integratoren2d_force.h b/src/Integratoren2d_force.h index d1c56fb..a3bc9c2 100644 --- a/src/Integratoren2d_force.h +++ b/src/Integratoren2d_force.h @@ -23,25 +23,23 @@ class Integratoren2d_force { static void Set3_Exact( Rod2d &rod2D, Simulation &sim, - std::function force, - std::function torque); + const std::function& force, + const std::function& torque); static void Set4_BDAS( Rod2d &rod2D, Simulation &sim, - std::function force, - std::function torque); + const std::function& force, + const std::function& torque); static const Eigen::Vector2d unitVec; - static void Set5_MBD(Rod2d &rod2D, Simulation &sim); + static void Set5_MBD(Rod2d &rod2D, Simulation &sim, + const std::function& force, + const std::function& torque); static void Set0_deterministic(Rod2d &rod2D, Simulation &sim); private: - static Eigen::Vector2d ortho(Eigen::Vector2d e); - - static Eigen::Matrix2d MatTraf(Eigen::Matrix e); - static Eigen::Vector2d Heun_predictor_pos( const Rod2d &rod2D, Simulation &sim, const std::function diff --git a/src/Integratoren2d_forceless.cpp b/src/Integratoren2d_forceless.cpp index cbc2c99..836220a 100644 --- a/src/Integratoren2d_forceless.cpp +++ b/src/Integratoren2d_forceless.cpp @@ -3,10 +3,7 @@ // #include "Integratoren2d_forceless.h" - -Eigen::Vector2d Integratoren2d_forceless::ortho(Eigen::Vector2d e) { - return Eigen::Vector2d{-e[1], e[0]}; -} +#include "vec_trafos.h" void Integratoren2d_forceless::Set0_deterministic(Rod2d &rod2D, Simulation & /*sim*/) { @@ -18,165 +15,109 @@ void Integratoren2d_forceless::Set0_deterministic(Rod2d &rod2D, void Integratoren2d_forceless::Set1_Euler(Rod2d &rod2D, Simulation &sim) { auto std = sim.getSTD(); // sqrt(2*delta_T) + Eigen::Vector2d e = rod2D.getE(); // translation - Eigen::Vector2d rand = {sim.getNorm(std), - sim.getNorm(std)}; // gaussian noise - Eigen::Vector2d trans_part = - rod2D.getDiff_Sqrt() * rand; // translation vector in Particle System - Eigen::Vector2d trans_lab = rod2D.getE_Base_matrix() * trans_part; - /* TODO - - Statt Zeile 22 und 23 kannst du trans_lab = paralleler_matrix_eintrag * e - plus senkrechter_matrix_eintrag * ortho(e); - - */ + Eigen::Vector2d rand = {sim.getNorm(std), sim.getNorm(std)}; + Eigen::Vector2d trans_part = rod2D.getDiff_Sqrt() * rand; + Eigen::Vector2d trans_lab = rotation_Matrix(e) * trans_part; // rotation - double rot = - sim.getNorm(std) * rod2D.getDRot_Sqrt(); // rotationsmatrix verwenden. - Eigen::Vector2d e = rod2D.getE(); + double rot = sim.getNorm(std) * rod2D.getDRot_Sqrt(); Eigen::Vector2d new_e = e + rot * ortho(e); - new_e.normalize(); + // apply rod2D.setPos(rod2D.getPos() + trans_lab); - rod2D.setE(new_e); + rod2D.setE(new_e.normalized()); } void Integratoren2d_forceless::Set2_Heun(Rod2d &rod2D, Simulation &sim) { // Predict with Euler auto std = sim.getSTD(); + Eigen::Vector2d e = rod2D.getE(); Eigen::Vector2d rand = Eigen::Vector2d(sim.getNorm(std), sim.getNorm(std)); Eigen::Vector2d trans_part = - rod2D.getDiff_Sqrt() * rand; // translation vector in Particle System - Eigen::Vector2d trans_lab = rod2D.getE_Base_matrix() * trans_part; + rod2D.getDiff_Sqrt() * rand; + Eigen::Vector2d trans_lab = rotation_Matrix(e) * trans_part; Eigen::Vector2d pos_predictor = rod2D.getPos() + trans_lab; - /* - - TODO - - Siehe Set_1: Hier auch über e und orth(e) - - */ - // rotation - double rot_predict = - sim.getNorm(std) * rod2D.getDRot_Sqrt(); // rotationsmatrix verwenden. - Eigen::Vector2d e = rod2D.getE(); + double rot_predict = sim.getNorm(std) * rod2D.getDRot_Sqrt(); Eigen::Vector2d delta_e_predict = rot_predict * ortho(e); Eigen::Vector2d e_predict = (e + delta_e_predict).normalized(); - double rot = sim.getNorm(std) * rod2D.getDRot_Sqrt(); Eigen::Vector2d e_integrated = e + 0.5 * (rot * ortho(e) + rot * ortho(e_predict)); - // TODO pos integration with future system ??? - // Eigen::Vector2d pos_integrated = rod2D.getPos() + 0.5 * (rot * ortho(e) - // + rot * ortho(e_predict)); // apply - rod2D.setPos(pos_predictor); + rod2D.setPos(pos_predictor); // pos with Euler rod2D.setE(e_integrated.normalized()); } -constexpr double kBT = 1.0; - void Integratoren2d_forceless::Set3_Exact(Rod2d &rod2D, Simulation &sim) { auto std = sim.getSTD(); + Eigen::Vector2d e = rod2D.getE(); // translation - Eigen::Vector2d rand = {sim.getNorm(std), - sim.getNorm(std)}; // gaussian noise - Eigen::Vector2d trans_part = - rod2D.getDiff_Sqrt() * rand; // translation vector in Particle System - Eigen::Vector2d trans_lab = rod2D.getE_Base_matrix() * trans_part; - - /* - - TODO: Siehe Set1 - - */ + Eigen::Vector2d rand = {sim.getNorm(std), sim.getNorm(std)}; + Eigen::Vector2d trans_part = rod2D.getDiff_Sqrt() * rand; + Eigen::Vector2d trans_lab = rotation_Matrix(e) * trans_part; // rotation double rot = sim.getNorm(std) * rod2D.getDRot_Sqrt(); - Eigen::Vector2d e = rod2D.getE(); auto correction = - -0.5 * pow(rod2D.getDRot(), 2) / pow(kBT, 2) * sim.getMDeltaT(); - + -0.5 * pow(rod2D.getDRot(), 2) * sim.getMDeltaT(); Eigen::Vector2d delta_e = correction * e + rot * ortho(e); - /* - TODO Warum hier kBT definiert? - */ - // apply rod2D.setPos(rod2D.getPos() + trans_lab); rod2D.setE((e + delta_e).normalized()); } -const Eigen::Vector2d Integratoren2d_forceless::unitVec = {1, 0}; - // this is a slow implementation to keep the same data structure for performance // analysis this must be altered +// Normalisation should not be necessary if a proper angular representation +// is used. But with vector e it is done in case of numerical errors void Integratoren2d_forceless::Set4_BDAS(Rod2d &rod2D, Simulation &sim) { auto std = sim.getSTD(); + Eigen::Vector2d e = rod2D.getE(); // translation auto rand = Eigen::Vector2d(sim.getNorm(std), sim.getNorm(std)); auto trans_part = rod2D.getDiff_Sqrt() * rand; // translation vector in Particle System - Eigen::Rotation2D rotMat(acos(unitVec.dot(rod2D.getE()))); - auto trans_lab = rotMat * trans_part; - - /* - TODO: Warum hier nicht die Matrix von oben? - - */ + Eigen::Vector2d trans_lab = rotation_Matrix(e) * trans_part; auto rot = Eigen::Rotation2D(sim.getNorm(std) * rod2D.getDRot_Sqrt()); Eigen::Rotation2Dd rotation(rot); - auto e_new = (rotation.toRotationMatrix() * rod2D.getE()).normalized(); - - /* - TODO Warum Normierung? Ist der Fehler so gros? Wie siehts mit einer - Winkelvariable aus, dann musst du den Winkel nicht jedesmal berechnen? - */ - - // Normalisation should not be necessary if a proper angular representation - // is used. But with vector e it is done in case of numerical errors + auto e_new = (rotation.toRotationMatrix() * e).normalized(); rod2D.setPos(rod2D.getPos() + trans_lab); rod2D.setE(e_new); } -Eigen::Matrix2d Integratoren2d_forceless::e_2_matrix(Eigen::Vector2d m_e) { - Eigen::Matrix2d mat; - mat << m_e[0], -m_e[1], m_e[1], m_e[0]; - return mat; -} +void Integratoren2d_forceless::Set5_MBD(Rod2d & rod2D, + Simulation & sim) { -void Integratoren2d_forceless::Set5_MBD(Rod2d & /*rod2D*/, - Simulation & /*sim*/) { - /* auto std = sim.getSTD(); - - auto rand = Eigen::Vector2d(sim.getNorm(std), sim.getNorm(std)); - auto trans_part = rod2D.getDiff_Sqrt() * rand; //translation vector in - Particle System auto trans_lab = rod2D.getE_Base_matrix() * trans_part; auto - pos_predictor = rod2D.getPos() + trans_lab; + Eigen::Vector2d e = rod2D.getE(); //rotation - auto rot_predict = Eigen::Rotation2D(sim.getNorm(std) * - rod2D.getDRot_Sqrt());// rotationsmatrix verwenden. auto e = rod2D.getE(); - auto delta_e_predict = rot_predict * ortho(e); auto e_predict = (e + - delta_e_predict).normalized(); + double rot_predict = sim.getNorm(std) * + rod2D.getDRot_Sqrt(); + auto delta_e_predict = rot_predict * ortho(e); + auto e_predict = (e + delta_e_predict).normalized(); - auto std_combined = rod2D.getDiff() + rod2D.getE_Base_matrix() * - rod2D.getDiff(); //old + new diff - auto rot = Eigen::Rotation2D(sim.getNorm(std) * rod2D.getDRot_Sqrt()); - auto e_integrated = e + 0.5 * (rot * ortho(e) + rot * ortho(e_predict)); + Eigen::Matrix2d Diff_combined = 0.5*(rod2D.getDiff() + rotation_Matrix(e).inverse() * rotation_Matrix(e_predict)*rod2D.getDiff()); + Eigen::Matrix2d D_combined_sqrt = Diff_combined.llt().matrixL(); + + auto rand = Eigen::Vector2d(sim.getNorm(std), sim.getNorm(std)); + Eigen::Vector2d pos_integrated =rod2D.getPos() + rotation_Matrix(e) * D_combined_sqrt * rand; + + double rot = sim.getNorm(std) * rod2D.getDRot_Sqrt(); + Eigen::Vector2d e_integrated = e + 0.5 * (rot * ortho(e) + rot * ortho(e_predict)); //apply - rod2D.setPos(pos_predictor); + rod2D.setPos(pos_integrated); rod2D.setE(e_integrated.normalized()); - */ + } diff --git a/src/Integratoren2d_forceless.h b/src/Integratoren2d_forceless.h index 9ba8344..efeb594 100644 --- a/src/Integratoren2d_forceless.h +++ b/src/Integratoren2d_forceless.h @@ -19,12 +19,7 @@ class Integratoren2d_forceless { static const Eigen::Vector2d unitVec; - static Eigen::Matrix2d e_2_matrix(Eigen::Vector2d m_e); - static void Set5_MBD(Rod2d &rod2D, Simulation &sim); static void Set0_deterministic(Rod2d &rod2D, Simulation &sim); - - private: - static Eigen::Vector2d ortho(Eigen::Vector2d e); }; diff --git a/src/LiveAgg.cpp b/src/LiveAgg.cpp index eda1b71..8b4fc34 100644 --- a/src/LiveAgg.cpp +++ b/src/LiveAgg.cpp @@ -19,3 +19,8 @@ void LiveAgg::feed(double value) noexcept { S += delta * delta2; } int LiveAgg::getNumPoints() const noexcept { return num_of_data; } + +std::ostream &operator<<(std::ostream &os, const LiveAgg &agg) { + os << "num_of_data: " << agg.num_of_data << " mean: " << agg.mean << " S: " << agg.S; + return os; +} diff --git a/src/LiveAgg.hpp b/src/LiveAgg.hpp index aa15eaf..7c591d1 100644 --- a/src/LiveAgg.hpp +++ b/src/LiveAgg.hpp @@ -1,4 +1,7 @@ #pragma once + +#include + class LiveAgg { public: LiveAgg(); @@ -8,7 +11,9 @@ class LiveAgg { [[nodiscard]] int getNumPoints() const noexcept; void feed(double value) noexcept; - private: + friend std::ostream &operator<<(std::ostream &os, const LiveAgg &agg); + +private: int num_of_data; double mean; double S; diff --git a/src/Rod2d.cpp b/src/Rod2d.cpp index 05e53e9..cf29906 100644 --- a/src/Rod2d.cpp +++ b/src/Rod2d.cpp @@ -52,8 +52,3 @@ const Eigen::Vector2d &Rod2d::getE() const { return m_e; } void Rod2d::setE(const Eigen::Vector2d &mE) { m_e = mE; } -Eigen::Matrix2d Rod2d::getE_Base_matrix() const { - Eigen::Matrix2d mat; - mat << m_e[0], -m_e[1], m_e[1], m_e[0]; - return mat; -} diff --git a/src/Rod2d.hpp b/src/Rod2d.hpp index f594b4c..dfb8079 100644 --- a/src/Rod2d.hpp +++ b/src/Rod2d.hpp @@ -22,8 +22,6 @@ class Rod2d { [[nodiscard]] double getDRot_Sqrt() const; - Eigen::Matrix2d getE_Base_matrix() const; - private: Eigen::Matrix2d m_Diff; Eigen::Matrix2d m_Diff_sqrt; diff --git a/src/legacy/Integratoren.cpp b/src/legacy/Integratoren.cpp deleted file mode 100644 index 3234c6e..0000000 --- a/src/legacy/Integratoren.cpp +++ /dev/null @@ -1,130 +0,0 @@ -#include "Integratoren.hpp" -constexpr double kBT = 1.0; - -void Integratoren::integrateITOEuler_1(Rod2d &rod) -{ - //DLOG_SCOPE_F(1,"Integrateing Rod2d:"); - double movePara = std::sqrt(2.0 * rod.Dpara * rod.deltaT) * rod.dis(rod.generator); - double moveOrtho = std::sqrt(2.0 * rod.Dortho * rod.deltaT) * rod.dis(rod.generator); - - double tmpx1 = rod.x1 + movePara * rod.e1 - moveOrtho * rod.e2; - double tmpx2 = rod.x2 + movePara * rod.e2 + moveOrtho * rod.e1; - - double rodOrtho = std::sqrt(2.0 * rod.Drot * rod.deltaT) * rod.dis(rod.generator); - double tmpe1 = rod.e1 - rodOrtho * rod.e2; - double tmpe2 = rod.e2 + rodOrtho * rod.e1; - - rod.x1 = tmpx1; - rod.x2 = tmpx2; - rod.e1 = tmpe1; - rod.e2 = tmpe2; - rod.normalize(); -} - -void Integratoren::integrateITOHeun_2(Rod2d &rod) -{ - //DLOG_SCOPE_F(1,"Integrateing Rod2d:"); - - double movePara = std::sqrt(2.0 * rod.Dpara * rod.deltaT) * rod.dis(rod.generator); - double moveOrtho = std::sqrt(2.0 * rod.Dortho * rod.deltaT) * rod.dis(rod.generator); - - double tmpx1 = rod.x1 + movePara * rod.e1 - moveOrtho * rod.e2; - double tmpx2 = rod.x2 + movePara * rod.e2 + moveOrtho * rod.e1; - - - double rodOrtho = std::sqrt(2.0 * rod.Drot * rod.deltaT) * rod.dis(rod.generator); - double e1 = rod.e1 - rodOrtho * rod.e1; - double e2 = rod.e2 + rodOrtho * rod.e2; - double norm = std::sqrt(e1 * e1 + e2 * e2); - e1 /= norm; - e2 /= norm; - - double rodOrtho1 = std::sqrt(2.0 * rod.Drot * rod.deltaT); - double rodOrthoRand = rod.dis(rod.generator); - double tmpe1 = rod.e1 - 0.5 * (rodOrtho1 * rod.e2 + rodOrtho1 * e2) * rodOrthoRand;//TODO Check if 0.5 faktor - double tmpe2 = rod.e2 + 0.5 * (rodOrtho1 * rod.e1 + rodOrtho1 * e1) * rodOrthoRand; - - rod.x1 = tmpx1; - rod.x2 = tmpx2; - rod.e1 = tmpe1; - rod.e2 = tmpe2; - rod.normalize(); -} - - -void Integratoren::integrateExact_3(Rod2d &rod) -{ - //DLOG_SCOPE_F(1,"Integrateing Rod2d:"); - - double movePara = std::sqrt(2.0 * rod.Dpara * rod.deltaT) * rod.dis(rod.generator); - double moveOrtho = std::sqrt(2.0 * rod.Dortho * rod.deltaT) * rod.dis(rod.generator); - - double tmpx1 = rod.x1 + movePara * rod.e1 - moveOrtho * rod.e2; - double tmpx2 = rod.x2 + movePara * rod.e2 + moveOrtho * rod.e1; - - double rodOrtho = std::sqrt(2.0 * rod.Drot * rod.deltaT) * rod.dis(rod.generator); - double korr = -0.5 * rod.Drot * rod.Drot / kBT / kBT * rod.deltaT; - double tmpe1 = rod.e1 - rodOrtho * rod.e2 + korr * rod.e1; - double tmpe2 = rod.e2 + rodOrtho * rod.e1 + korr * rod.e2; - - rod.x1 = tmpx1; - rod.x2 = tmpx2; - rod.e1 = tmpe1; - rod.e2 = tmpe2; - rod.normalize(); -} - -void Integratoren::integrateBDAS_4(Rod2d &rod) -{ - double deltaRTpara = std::sqrt(2.0 * rod.Dpara * rod.deltaT) * rod.dis(rod.generator); - double deltaRTortho = std::sqrt(2.0 * rod.Dortho * rod.deltaT) * rod.dis(rod.generator); - double tmpx1 = rod.x1 + deltaRTpara * rod.e1 - deltaRTortho * rod.e2; - double tmpx2 = rod.x2 + deltaRTpara * rod.e2 + deltaRTortho * rod.e1; - - double deltaPhi = std::sqrt(2.0 * rod.Drot * rod.deltaT) * rod.dis(rod.generator); - double tmpe1 = rod.e1 * std::cos(deltaPhi) - rod.e2 * std::sin(deltaPhi); - double tmpe2 = rod.e1 * std::sin(deltaPhi) + rod.e2 * std::cos(deltaPhi); - - rod.x1 = tmpx1; - rod.x2 = tmpx2; - rod.e1 = tmpe1; - rod.e2 = tmpe2; - rod.normalize(); -} - -void Integratoren::integrateMBD_5(Rod2d &rod) -{ - /* - double movePara = std::sqrt(2.0*rod.Dpara*rod.deltaT)*rod.dis(rod.generator); // Taken from euler - double moveOrtho = std::sqrt(2.0*rod.Dortho*rod.deltaT)*rod.dis(rod.generator); // taken from euler - double x1 = rod.x1 + movePara*rod.e1- moveOrtho*rod.e2; - double x2 = rod.x2 + movePara*rod.e2+ moveOrtho*rod.e1; - */ - // unused without force - double rodOrtho = std::sqrt(2.0 * rod.Drot * rod.deltaT) * rod.dis(rod.generator); - double e1 = rod.e1 - rodOrtho * rod.e2; - double e2 = rod.e2 + rodOrtho * rod.e1; - double norm = std::sqrt(e1 * e1 + e2 * e2); - e1 /= norm; - e2 /= norm; - - double movePara1 = std::sqrt(1.0 * rod.Dpara * rod.deltaT) * rod.dis(rod.generator);// If difussion would be time dependend, average of the predicted D and now D - double moveOrtho1 = std::sqrt(1.0 * rod.Dortho * rod.deltaT) * rod.dis(rod.generator); - double movePara2 = std::sqrt(1.0 * rod.Dpara * rod.deltaT) * rod.dis(rod.generator);// If difussion would be time dependend, average of the predicted D and now D - double moveOrtho2 = std::sqrt(1.0 * rod.Dortho * rod.deltaT) * rod.dis(rod.generator); - double tmpx1 = rod.x1 + (movePara1 * rod.e1 - moveOrtho1 * rod.e2) + (movePara2 * e1 - moveOrtho2 * e2); - double tmpx2 = rod.x2 + (movePara1 * rod.e2 + moveOrtho1 * rod.e1) + (movePara2 * e2 + moveOrtho2 * e1); - - - double rodOrtho1 = std::sqrt(2.0 * rod.Drot * rod.deltaT); - double rodOrtho2 = rod.dis(rod.generator); - - double tmpe1 = rod.e1 - 0.5 * (rodOrtho1 * rod.e2 + rodOrtho1 * e2) * rodOrtho2; - double tmpe2 = rod.e2 + 0.5 * (rodOrtho1 * rod.e1 + rodOrtho1 * e1) * rodOrtho2; - - rod.x1 = tmpx1; - rod.x2 = tmpx2; - rod.e1 = tmpe1; - rod.e2 = tmpe2; - rod.normalize(); -} diff --git a/src/legacy/Integratoren.hpp b/src/legacy/Integratoren.hpp deleted file mode 100644 index bcb0168..0000000 --- a/src/legacy/Integratoren.hpp +++ /dev/null @@ -1,19 +0,0 @@ -#ifndef INTEGRATOREN_HPP -#define INTEGRATOREN_HPP -#include "Rod2d.hpp" - - -class Integratoren -{ -public: - static void integrateITOEuler_1(Rod2d &rod); - - static void integrateITOHeun_2(Rod2d &rod); - - static void integrateExact_3(Rod2d &rod); - - static void integrateBDAS_4(Rod2d &rod); - - static void integrateMBD_5(Rod2d &rod); -}; -#endif// INTEGRATOREN_HPP diff --git a/src/legacy/Simulation2d.cpp b/src/legacy/Simulation2d.cpp deleted file mode 100644 index 3c9b5ed..0000000 --- a/src/legacy/Simulation2d.cpp +++ /dev/null @@ -1,6 +0,0 @@ -#include "Simulation2d.hpp" -#include -#include -#include "LiveAgg.hpp" -template -std::vector simulate(double deltaT, int seed, int meanSteps); diff --git a/src/legacy/Simulation2d.hpp b/src/legacy/Simulation2d.hpp deleted file mode 100644 index 2c66b4b..0000000 --- a/src/legacy/Simulation2d.hpp +++ /dev/null @@ -1,157 +0,0 @@ -#ifndef SIMULATION2D_HPP -#define SIMULATION2D_HPP -#include "Rod2d.hpp" -#include "LiveAgg.hpp" -#include -#include -/* This is used to calculate the convergence order of the particles -*/ - -inline void printProgress(double progress) -{ - const int barWidth = 70; - std::cout << "["; - const int pos = static_cast(barWidth * progress); - for (int i = 0; i < barWidth; ++i) { - if (i < pos) { - std::cout << "="; - } else if (i == pos) { - std::cout << ">"; - } else { - std::cout << " "; - } - } - std::cout << "] " << int(progress * 100.0) << " % \r"; - std::cout.flush(); -} - -template -void simulate(double L, int seed, const std::string &filename) -{ - // we now simulate the msd after 10 steps - - - Rod2d rod(L, 1e-5, seed); - const int sampleWidth = 1000; - std::vector msd(sampleWidth, LiveAgg()); - std::vector oaf(sampleWidth, LiveAgg()); - std::vector empxx(sampleWidth, LiveAgg()); - std::vector empyy(sampleWidth, LiveAgg()); - std::vector xPos(sampleWidth, 0.0); - std::vector yPos(sampleWidth, 0.0); - std::vector e1(sampleWidth, 0.0); - std::vector e2(sampleWidth, 0.0); - std::vector targetMSD(sampleWidth, 0.0); - std::vector targetOAF(sampleWidth, 0.0); - std::vector targetEMPXX(sampleWidth, 0.0); - std::vector targetEMPYY(sampleWidth, 0.0); - std::vector samplePoints(sampleWidth, 1); - size_t temp = 1; - const double Dmean = 0.5 * (rod.Dpara + rod.Dortho); - const double deltaD = rod.Dortho - rod.Dpara; - for (size_t i = 0; i < samplePoints.size(); ++i) { - samplePoints[i] = temp; - temp += 100; - const double time = static_cast(temp) * rod.deltaT; - targetMSD[i] = 4.0 * 0.5 * (rod.Dpara + rod.Dortho) * time; - targetOAF[i] = std::exp(-rod.Drot * time); - const double u = 4.0; - targetEMPXX[i] = Dmean - deltaD / 2.0 * (1 - exp(-u * rod.Drot * time)) / u / rod.Drot / time; - targetEMPYY[i] = Dmean + deltaD / 2.0 * (1 - exp(-u * rod.Drot * time)) / u / rod.Drot / time; - } - unsigned long long int iteration = 0; - //const unsigned long long int maxItr = ((long int)10000)*1000000; - const unsigned long long int maxItr = 10000000; - const unsigned long long int updateInter = 10000; - while (iteration < maxItr) { - T(rod); - iteration++; - for (size_t i = 0; i < samplePoints.size(); ++i) { - const size_t sample = samplePoints[i]; - if (iteration % sample == 0) { - msd[i].feed((rod.x1 - xPos[i]) * (rod.x1 - xPos[i]) + (rod.x2 - yPos[i]) * (rod.x2 - yPos[i])); - oaf[i].feed(rod.e1 * e1[i] + rod.e2 * e2[i]); - empxx[i].feed(std::pow((rod.x1 - xPos[i]) * e1[i] + (rod.x2 - yPos[i]) * e2[i], 2)); - empyy[i].feed(std::pow((rod.x1 - xPos[i]) * e2[i] - (rod.x2 - yPos[i]) * e1[i], 2)); - xPos[i] = rod.x1; - yPos[i] = rod.x2; - e1[i] = rod.e1; - e2[i] = rod.e2; - } - } - if (iteration % updateInter == 0) { - printProgress(static_cast(iteration)/ maxItr); - } - } - - std::cout << "\n"; - - std::ofstream fileMSD(filename + "_" + std::to_string(L) + "_MSD.out"); - std::ofstream fileOAF(filename + "_" + std::to_string(L) + "_OAF.out"); - std::ofstream fileXX(filename + "_" + std::to_string(L) + "_XX.out"); - std::ofstream fileYY(filename + "_" + std::to_string(L) + "_YY.out"); - - for (size_t i = 0; i < samplePoints.size(); ++i) { - fileMSD << static_cast(samplePoints[i]) * rod.deltaT << "\t" << targetMSD[i] << "\t" << msd[i].getMean() << "\t" << msd[i].getSEM() << "\t" << msd[i].getSD() << "\t" << msd[i].getNumPoints() << "\n"; - fileOAF << static_cast(samplePoints[i]) * rod.deltaT << "\t" << targetOAF[i] << "\t" << oaf[i].getMean() << "\t" << oaf[i].getSEM() << "\t" << oaf[i].getSD() << "\t" << oaf[i].getNumPoints() << "\n"; - fileXX << static_cast(samplePoints[i]) * rod.deltaT << "\t" << targetEMPXX[i] << "\t" << empxx[i].getMean() << "\t" << empxx[i].getSEM() << "\t" << empxx[i].getSD() << "\t" << empxx[i].getNumPoints() << "\n"; - fileYY << static_cast(samplePoints[i]) * rod.deltaT << "\t" << targetEMPYY[i] << "\t" << empyy[i].getMean() << "\t" << empyy[i].getSEM() << "\t" << empyy[i].getSD() << "\t" << empyy[i].getNumPoints() << "\n"; - } - fileMSD.close(); - fileOAF.close(); - fileXX.close(); - fileYY.close(); -} - - -template -void konvergenz(double L, int seed, double deltaT, const std::string &filename) -{ - // we now simulate the msd after 10 steps - Rod2d rod(L, deltaT, seed); - const double time = 2; - const int steps = static_cast(time / deltaT); - const int minSteps = 1000; - - //const double targetMSD = 4.0*0.5*(rod.Dpara+rod.Dortho)*time; - const double targetOAF = std::exp(-rod.Drot * time); - /*const double u = 4.0; - const double Dmean = 0.5*(rod.Dpara+rod.Dortho); - const double deltaD = rod.Dortho - rod.Dpara; - const double targetEMPXX = Dmean-deltaD/2.0*(1 - exp(-u*rod.Drot*time))/u/rod.Drot / time; - const double targetEMPYY = Dmean+deltaD/2.0*(1 - exp(-u*rod.Drot*time))/u/rod.Drot / time;*/ - //LiveAgg msd; - LiveAgg oaf; - //LiveAgg empxx; - //LiveAgg empyy; - int meanCount = 0; - while (meanCount < minSteps or std::abs(targetOAF - oaf.getMean()) < oaf.getSEM() * 10) { - //reset the position - rod.reset(); - for (int i = 0; i < steps; ++i) { - T(rod); - } - //msd.feed(rod.x1*rod.x1+rod.x2*rod.x2); - oaf.feed(rod.e1); - //.feed(rod.x1*rod.x1/2/time); - //empyy.feed(rod.x2*rod.x2/2/time); - if (meanCount % 100000 == 0) { - //std::cout< -#include -#include -#include "legacy/Integratoren.hpp" -#include "legacy/Simulation2d.hpp" -const int seed = 12345; -[[maybe_unused]] void overTime(int sim) -{ - const std::vector Length = { 1.0, 1.3974, 20.0 }; - const int numCases = 5; - int i = 0; - for (auto L : Length) { - if (sim == numCases * i) { - simulate(L, seed, "1"); - } - if (sim == numCases * i + 1) { - simulate(L, seed, "2"); - } - if (sim == numCases * i + 2) { - simulate(L, seed, "3"); - } - if (sim == numCases * i + 3) { - simulate(L, seed, "4"); - } - if (sim == numCases * i + 4) { - simulate(L, seed, "5"); - } - i++; - } -} - -[[maybe_unused]]void konvergenz(int sim) -{ - const std::vector Length = { 1.0, 1.3974, 20.0 }; - const std::vector DELTAT = { 2.0, 1.0, 0.5, 0.25, 0.2, 0.1, 0.05, 0.025, 0.02, 0.01, 0.005, 0.0025, 0.200, 0.001, 0.0005, 0.00025 }; - int i = 0; - for (auto dt : DELTAT) { - - for (auto L : Length) { - if (sim == i++) { - konvergenz(L, seed, dt, "1"); - } - if (sim == i++) { - konvergenz(L, seed, dt, "2"); - } - if (sim == i++) { - konvergenz(L, seed, dt, "3"); - } - if (sim == i++) { - konvergenz(L, seed, dt, "4"); - } - if (sim == i++) { - konvergenz(L, seed, dt, "5"); - } - } - } -} - - -int main(int argc, char *argv[]) -{ - int sim; - auto args = std::span(argv, size_t(argc)); - if (argc > 1) { - sim = std::stoi(args[1]) - 1; - } else { - sim = 0; - } - overTime(sim); - //konvergenz(sim); - return 0; -} diff --git a/src/main.cpp b/src/main.cpp index df9b018..06bed95 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -9,7 +9,7 @@ #include "Integratoren2d_forceless.h" int main() { - constexpr int numStep = 1000000; + constexpr int numStep = 10000000; constexpr double k = 0.01; [[maybe_unused]] auto harmonic_Force = [](const Eigen::Vector2d &pos, @@ -21,55 +21,31 @@ int main() { const Eigen::Vector2d & /*torque*/) -> Eigen::Vector2d { return {0.0, 0.0}; }; - auto zero_Torque = [](const Eigen::Vector2d & /*pos*/, + [[maybe_unused]] auto zero_Torque = [](const Eigen::Vector2d & /*pos*/, const Eigen::Vector2d & /*torque*/) -> double { return 0.0; }; - /* auto euler_Zero = [&](Rod2d &rod, Simulation &sim) { - return Integratoren2d_force::Set1_Euler(rod, sim, zero_Force, - zero_Torque); - };*/ - auto heun_Zero = [&](Rod2d &rod, Simulation &sim) { - return Integratoren2d_force::Set1_Euler(rod, sim, harmonic_Force, - zero_Torque); - }; - /*auto exact_Zero = [&](Rod2d &rod, Simulation &sim) { - return Integratoren2d_force::Set3_Exact(rod, sim, zero_Force, - zero_Torque); - };*/ - /* - auto bdas_Zero = [&](Rod2d &rod, Simulation &sim) { - return Integratoren2d_force::Set4_BDAS(rod, sim, zero_Force, zero_Torque); - };*/ - Calculation euler(heun_Zero, - {{Compute::Type::msd, 1}, - {Compute::Type::msd, 2}, - {Compute::Type::msd, 4}, - {Compute::Type::msd, 8}, - {Compute::Type::msd, 16}, - {Compute::Type::msd, 32}, - {Compute::Type::msd, 64}, - {Compute::Type::msd, 128}, - {Compute::Type::msd, 256}, - {Compute::Type::msd, 512}, - {Compute::Type::msd, 1024}, - {Compute::Type::oaf, 1}, - {Compute::Type::oaf, 2}, - {Compute::Type::oaf, 4}, - {Compute::Type::oaf, 8}, - {Compute::Type::oaf, 16}}, - 0.01, 12345); + + 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::Type::msd) continue; - std::cout << "MSD: " << com.getDifference() << " " - << com.getAgg().getMean() << " <-> " << com.getTarget() - << std::endl; - } - for (const auto &com : euler.getComputes()) { - if (com.getType() != Compute::Type::oaf) continue; - std::cout << "OAF: " << com.getDifference() << " " - << com.getAgg().getMean() << " <-> " << com.getTarget() + if (com.getType() != compute) { continue; +} + std::cout + << com.getAgg().getMean() << " , " << com.getTarget() << std::endl; } } \ No newline at end of file diff --git a/src/vec_trafos.h b/src/vec_trafos.h new file mode 100644 index 0000000..963b812 --- /dev/null +++ b/src/vec_trafos.h @@ -0,0 +1,13 @@ +// +// Created by jholder on 10/25/21. +// +#pragma once +#include +inline static Eigen::Vector2d ortho(Eigen::Vector2d e) { + return Eigen::Vector2d{-e[1], e[0]}; +} +inline static Eigen::Matrix2d rotation_Matrix(Eigen::Vector2d e) { + Eigen::Matrix2d mat; + mat << e[0], -e[1], e[1], e[0]; + return mat; +} \ No newline at end of file diff --git a/test/test_Integratoren.cpp b/test/test_Integratoren.cpp index a618bd8..e1bc7c3 100644 --- a/test/test_Integratoren.cpp +++ b/test/test_Integratoren.cpp @@ -10,26 +10,202 @@ TEST_CASE("Forceless Integratoren") { const size_t SEED = Catch::rngSeed(); - auto force = [](const Eigen::Vector2d & /*pos*/, - const Eigen::Vector2d & /*torque*/) -> Eigen::Vector2d { - return {0.0, 0.0}; - }; - auto torque = [](const Eigen::Vector2d & /*pos*/, - const Eigen::Vector2d & /*torque*/) -> double { + const double length = GENERATE(1.0,1.4,2.0); + constexpr int steps = 10000; + constexpr double delta = 0.1; + + + Calculation euler(Integratoren2d_forceless::Set1_Euler, + {{Compute::Type::msd, 1}, + {Compute::Type::oaf, 1}, + {Compute::Type::empxx, 1}, + {Compute::Type::empyy, 1}}, + 0.01, SEED,length); + euler.run(steps); + + Calculation heun(Integratoren2d_forceless::Set2_Heun, + {{Compute::Type::msd, 1}, + {Compute::Type::oaf, 1}, + {Compute::Type::empxx, 1}, + {Compute::Type::empyy, 1}}, + 0.01, SEED,length); + heun.run(steps); + + Calculation exact(Integratoren2d_forceless::Set3_Exact, + {{Compute::Type::msd, 1}, + {Compute::Type::oaf, 1}, + {Compute::Type::empxx, 1}, + {Compute::Type::empyy, 1}}, + 0.01, SEED,length); + exact.run(steps); + Calculation bdas(Integratoren2d_forceless::Set4_BDAS, + {{Compute::Type::msd, 1}, + {Compute::Type::oaf, 1}, + {Compute::Type::empxx, 1}, + {Compute::Type::empyy, 1}}, + 0.01, SEED,length); + bdas.run(steps); + Calculation mbd(Integratoren2d_forceless::Set5_MBD, + {{Compute::Type::msd, 1}, + {Compute::Type::oaf, 1}, + {Compute::Type::empxx, 1}, + {Compute::Type::empyy, 1}}, + 0.01, SEED,length); + mbd.run(steps); + SECTION("ForceLess Euler") { + CAPTURE(length); + size_t i = 0; + for (auto &c: euler.getComputes()) { + CAPTURE(c); + CHECK(c.getDifference() <= delta*c.getAgg().getMean()); + ++i; + } + } + SECTION("ForceLess Heun") { + CAPTURE(length); + size_t i = 0; + for (auto &c: heun.getComputes()) { + CAPTURE(c); + CHECK(c.getDifference() <= delta*c.getAgg().getMean()); + ++i; + } + } + SECTION("ForceLess Exact") { + CAPTURE(length); + size_t i = 0; + for (auto &c: exact.getComputes()) { + CAPTURE(c); + CHECK(c.getDifference() <= delta*c.getAgg().getMean()); + ++i; + } + } + SECTION("ForceLess BDAS") { + CAPTURE(length); + size_t i = 0; + for (auto &c: bdas.getComputes()) { + CAPTURE(c); + CHECK(c.getDifference() <= delta*c.getAgg().getMean()); + ++i; + } + } + SECTION("ForceLess MBD") { + CAPTURE(length); + size_t i = 0; + for (auto &c: mbd.getComputes()) { + CAPTURE(c); + CHECK(c.getDifference() <= delta*c.getAgg().getMean()); + ++i; + } + } + + [[maybe_unused]] auto zero_Force = + [](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 { return 0.0; }; + Calculation euler_zero(Integratoren2d_force::Set1_Euler, + {{Compute::Type::msd, 1}, + {Compute::Type::oaf, 1}, + {Compute::Type::empxx, 1}, + {Compute::Type::empyy, 1}}, + 0.01, SEED,zero_Force,zero_Torque,length); + euler_zero.run(steps); + + Calculation heun_zero(Integratoren2d_force::Set2_Heun, + {{Compute::Type::msd, 1}, + {Compute::Type::oaf, 1}, + {Compute::Type::empxx, 1}, + {Compute::Type::empyy, 1}}, + 0.01, SEED,zero_Force,zero_Torque,length); + heun_zero.run(steps); + + Calculation exact_zero(Integratoren2d_force::Set3_Exact, + {{Compute::Type::msd, 1}, + {Compute::Type::oaf, 1}, + {Compute::Type::empxx, 1}, + {Compute::Type::empyy, 1}}, + 0.01, SEED,zero_Force,zero_Torque,length); + exact_zero.run(steps); + + Calculation bdas_zero(Integratoren2d_force::Set4_BDAS, + {{Compute::Type::msd, 1}, + {Compute::Type::oaf, 1}, + {Compute::Type::empxx, 1}, + {Compute::Type::empyy, 1}}, + 0.01, SEED,zero_Force,zero_Torque,length); + bdas_zero.run(steps); + + Calculation mbd_zero(Integratoren2d_force::Set5_MBD, + {{Compute::Type::msd, 1}, + {Compute::Type::oaf, 1}, + {Compute::Type::empxx, 1}, + {Compute::Type::empyy, 1}}, + 0.01, SEED,zero_Force,zero_Torque,length); + mbd_zero.run(steps); + SECTION("Force Euler") { + CAPTURE(length); + size_t i = 0; + for (auto &c: euler_zero.getComputes()) { + CAPTURE(c); + CHECK(c.getDifference() <= delta*c.getAgg().getMean()); + ++i; + } + } + SECTION("Force Heun") { + CAPTURE(length); + + size_t i = 0; + for (auto &c: heun_zero.getComputes()) { + CAPTURE(c); + CHECK(c.getDifference() <= delta*c.getAgg().getMean()); + ++i; + } + } + SECTION("Force Exact") { + CAPTURE(length); + size_t i = 0; + for (auto &c: exact_zero.getComputes()) { + CAPTURE(c); + CHECK(c.getDifference() <= delta*c.getAgg().getMean()); + ++i; + } + } + SECTION("Force BDAS") { + CAPTURE(length); + size_t i = 0; + for (auto &c: bdas_zero.getComputes()) { + CAPTURE(c); + CHECK(c.getDifference() <= delta*c.getAgg().getMean()); + ++i; + } + } + SECTION("Force MBD") { + CAPTURE(length); + size_t i = 0; + for (auto &c: mbd_zero.getComputes()) { + CAPTURE(c); + CHECK(c.getDifference() <= delta*c.getAgg().getMean()); + ++i; + } + } + +} + +/* + + + + + + SECTION("Euler") { - Calculation euler(Integratoren2d_forceless::Set1_Euler, - {{Compute::Type::msd, 1}, - {Compute::Type::oaf, 1}, - {Compute::Type::empxx, 1}, - {Compute::Type::empyy, 1}}, - 0.01, SEED); + { - euler.run(10000); - for (auto &c : euler.getComputes()) { - CHECK(c.getDifference() <= 0.005); - } + } Calculation euler_force( @@ -143,4 +319,39 @@ TEST_CASE("Forceless Integratoren") { Approx(bdas_force.getComputes()[1].getAgg().getMean()) .epsilon(10e-10)); } + + SECTION("MBD") { + Calculation mbd(Integratoren2d_forceless::Set5_MBD, + {{Compute::Type::msd, 1}, + {Compute::Type::oaf, 1}, + {Compute::Type::empxx, 1}, + {Compute::Type::empyy, 1}}, + 0.01, SEED); + { + mbd.run(10000); + for (auto &c : mbd.getComputes()) { + CHECK(c.getDifference() <= 0.005); + } + } + + Calculation mbd_force(Integratoren2d_force::Set5_MBD, + {{Compute::Type::msd, 1}, + {Compute::Type::oaf, 1}, + {Compute::Type::empxx, 1}, + {Compute::Type::empyy, 1}}, + 0.01, SEED, force, torque); + { + mbd_force.run(10000); + for (auto &c : mbd_force.getComputes()) { + CHECK(c.getDifference() <= 0.005); + } + } + CHECK(mbd.getComputes()[0].getAgg().getMean() == + Approx(mbd_force.getComputes()[0].getAgg().getMean()) + .epsilon(10e-10)); + CHECK(mbd.getComputes()[1].getAgg().getMean() == + Approx(mbd_force.getComputes()[1].getAgg().getMean()) + .epsilon(10e-10)); + } } +*/ \ No newline at end of file diff --git a/test/test_Rod2d.cpp b/test/test_Rod2d.cpp index 3cbdffd..1a77b1a 100644 --- a/test/test_Rod2d.cpp +++ b/test/test_Rod2d.cpp @@ -6,7 +6,7 @@ #include #include "Rod2d.hpp" - +#include "vec_trafos.h" TEST_CASE("Rods") { Rod2d sphere(1); Rod2d rod(2); @@ -21,41 +21,13 @@ TEST_CASE("Rods") { REQUIRE(sphere.getDiff_Sqrt()(1, 1) == 1.0); } SECTION("Checking rotation") { - SECTION("Forwards == e 90") { - // - // ===== - // - rod.setE({0, 1}); // Lab system is 90 degree rotated + double f1 = GENERATE(take(10, random(-100, 100))); + double f2 = GENERATE(take(10, random(-100, 100))); + rod.setE(Eigen::Vector2d({f1, f2}).normalized()); auto test = Eigen::Vector2d({1, 0}); - auto erg = (rod.getE_Base_matrix() * test).normalized(); auto e = rod.getE(); + auto erg = (rotation_Matrix(e) * test).normalized(); REQUIRE(erg.isApprox(e)); - } - SECTION("Forwards == e 45") { - // o - // o - // o - rod.setE(Eigen::Vector2d({1, 1}) - .normalized()); // Lab system is 90 degree rotated - auto test = Eigen::Vector2d({1, 0}); - auto erg = (rod.getE_Base_matrix() * test).normalized(); - auto e = rod.getE(); - REQUIRE(erg.isApprox(e)); - } - } - SECTION("Angle Stuff") { - const Eigen::Vector2d unitVec = {1, 0}; - double f1 = GENERATE(take(10, random(-100, 100))); - double f2 = GENERATE(take(10, random(0, 100))); - sphere.setE(Eigen::Vector2d({f1, f2}).normalized()); - Eigen::Rotation2D rotMat(acos(unitVec.dot(sphere.getE()))); - CHECK(rotMat.toRotationMatrix()(0, 0) == - Approx(sphere.getE_Base_matrix()(0, 0)).epsilon(10e-10)); - CHECK(rotMat.toRotationMatrix()(1, 0) == - Approx(sphere.getE_Base_matrix()(1, 0)).epsilon(10e-10)); - CHECK(rotMat.toRotationMatrix()(0, 1) == - Approx(sphere.getE_Base_matrix()(0, 1)).epsilon(10e-10)); - CHECK(rotMat.toRotationMatrix()(1, 1) == - Approx(sphere.getE_Base_matrix()(1, 1)).epsilon(10e-10)); + } } \ No newline at end of file