From 1997edeae7f810fdfa3a0ae0a70716b11191acb2 Mon Sep 17 00:00:00 2001 From: Jacob Date: Wed, 27 Oct 2021 16:11:55 +0200 Subject: [PATCH] A lot of refactoring --- .clang-tidy | 2 +- benchmark/bench.cpp | 2 +- benchmark/bench_euler.cpp | 2 +- benchmark/bench_mbd.cpp | 2 +- src/CMakeLists.txt | 2 +- src/Calculation.cpp | 17 +- src/Calculation.h | 4 + src/Compute.cpp | 65 +++- src/Compute.h | 15 +- src/Integratoren2d.cpp | 205 ++++++++++ src/Integratoren2d.h | 44 +++ src/Rod2d.cpp | 5 - src/Rod2d.hpp | 1 - src/force_lambdas.h | 25 ++ src/{ => legacy}/Integratoren2d_force.cpp | 11 +- src/{ => legacy}/Integratoren2d_force.h | 7 +- src/{ => legacy}/Integratoren2d_forceless.cpp | 0 src/{ => legacy}/Integratoren2d_forceless.h | 0 src/main.cpp | 152 +++++--- src/vec_trafos.h | 2 +- test/.clang-tidy | 7 + test/CMakeLists.txt | 2 +- test/test_Calculation.cpp | 64 ---- test/test_Compute.cpp | 1 - test/test_Integratoren.cpp | 359 ++---------------- 25 files changed, 510 insertions(+), 486 deletions(-) create mode 100644 src/Integratoren2d.cpp create mode 100644 src/Integratoren2d.h create mode 100644 src/force_lambdas.h rename src/{ => legacy}/Integratoren2d_force.cpp (94%) rename src/{ => legacy}/Integratoren2d_force.h (90%) rename src/{ => legacy}/Integratoren2d_forceless.cpp (100%) rename src/{ => legacy}/Integratoren2d_forceless.h (100%) create mode 100644 test/.clang-tidy delete mode 100644 test/test_Calculation.cpp diff --git a/.clang-tidy b/.clang-tidy index c5b7916..1b968e1 100644 --- a/.clang-tidy +++ b/.clang-tidy @@ -1,5 +1,5 @@ --- -Checks: '*,-fuchsia-*,-google-*,-zircon-*,-abseil-*,-modernize-use-trailing-return-type,-llvm*,-readability-magic-numbers,-cppcoreguidelines-avoid-magic-numbers' +Checks: '*,-fuchsia-*,-google-*,-zircon-*,-abseil-*,-modernize-use-trailing-return-type,-llvm*' WarningsAsErrors: '' HeaderFilterRegex: '' FormatStyle: none diff --git a/benchmark/bench.cpp b/benchmark/bench.cpp index f351722..85a1400 100644 --- a/benchmark/bench.cpp +++ b/benchmark/bench.cpp @@ -6,7 +6,7 @@ #define CATCH_CONFIG_ENABLE_BENCHMARKING #include -#include "Integratoren2d_forceless.h" +#include "legacy/Integratoren2d_forceless.h" #include "Rod2d.hpp" TEST_CASE("Euler - Baseline", "[benchmark]") { diff --git a/benchmark/bench_euler.cpp b/benchmark/bench_euler.cpp index 5034487..0f2cff6 100644 --- a/benchmark/bench_euler.cpp +++ b/benchmark/bench_euler.cpp @@ -3,7 +3,7 @@ // #include "Rod2d.hpp" #include "Simulation.h" -#include "Integratoren2d_forceless.h" +#include "legacy/Integratoren2d_forceless.h" int main(){ Rod2d rod(1.0); Simulation sim(0.01, 1); diff --git a/benchmark/bench_mbd.cpp b/benchmark/bench_mbd.cpp index d45e6e9..362c744 100644 --- a/benchmark/bench_mbd.cpp +++ b/benchmark/bench_mbd.cpp @@ -4,7 +4,7 @@ #include "Rod2d.hpp" #include "Simulation.h" -#include "Integratoren2d_forceless.h" +#include "legacy/Integratoren2d_forceless.h" int main(){ Rod2d rod(1.0); Simulation sim(0.01, 1); diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 3c3038a..9249c05 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -4,7 +4,7 @@ find_package(Eigen3) # Generic test that uses conan libs add_library(integratoren SHARED LiveAgg.cpp Rod2d.cpp Simulation.cpp - Simulation.h Integratoren2d_forceless.cpp Integratoren2d_force.cpp + Simulation.h Integratoren2d.cpp Calculation.cpp Compute.cpp) target_include_directories(integratoren PUBLIC .) target_link_libraries( diff --git a/src/Calculation.cpp b/src/Calculation.cpp index 18a0f17..cdc2c5d 100644 --- a/src/Calculation.cpp +++ b/src/Calculation.cpp @@ -21,7 +21,7 @@ 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)) { + : 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)); } @@ -30,7 +30,7 @@ Calculation::Calculation( 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) { + : rod(length), start_rod(rod), sim(deltaT, seed) { m_integrator = [&](Rod2d &t_rod, Simulation &t_sim) { t_integrator(t_rod, t_sim, force, torque); }; @@ -46,7 +46,7 @@ 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), + 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)); @@ -55,13 +55,18 @@ Calculation::Calculation(std::function t_integrator 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) { + : 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)); } - rod.setPos({1000,1000}); - +} + +void Calculation::reset() { + rod =start_rod; + for(auto & c:computes){ + c.reset(start_rod); + } } diff --git a/src/Calculation.h b/src/Calculation.h index 3e4a366..973a3ea 100644 --- a/src/Calculation.h +++ b/src/Calculation.h @@ -12,11 +12,15 @@ class Calculation { private: Rod2d rod; + const Rod2d start_rod; Simulation sim; std::function m_integrator; std::vector computes; public: + + void reset(); + [[nodiscard]] const Simulation &getSim() const; using force_type = std::function; diff --git a/src/Compute.cpp b/src/Compute.cpp index 433f595..f2b17d1 100644 --- a/src/Compute.cpp +++ b/src/Compute.cpp @@ -17,6 +17,13 @@ void Compute::evalMSD(const Rod2d &rod2D) { agg.feed(msd); } +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(); + agg.feed(msd); +} + void Compute::evalOAF(const Rod2d &rod2D) { const auto &new_e = rod2D.getE(); auto old_e = start_rod.getE(); @@ -57,51 +64,86 @@ void Compute::eval(const Rod2d &rod2D) { case Type::empyy: eval_empYY(rod2D); break; - case Type::msd_harmonic_k_10: + case Type::msd_harmonic_k_1: + evalMSD(rod2D); + break; + case Type::msd_const: + break; + case Type::x: + evalX(rod2D); + break; + case Type::x_squared: evalMSD(rod2D); break; } - start_rod = rod2D; + if (resetting) { + start_rod = rod2D; + } + return; } } 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) { - time = sim.getMDeltaT() * static_cast(every); + : 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: { + resetting = true; type_str = "msd"; - target = 4.0 * 0.5 * (rod.getDiff().trace()) * time; + target = 2 * (rod.getDiff().trace()) * time; break; } case Type::oaf: { + resetting = true; type_str = "oaf"; target = std::exp(-rod.getDRot() * time); break; } case Type::empxx: { + resetting = true; 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); - target = Dmean - deltaD / 2.0 * + target = Dmean - deltaD / 2 * (1 - exp(-u * rod.getDRot() * time)) / u / rod.getDRot() / time; } break; case Type::empyy: { + resetting = true; 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 * + target = Dmean + deltaD / 2 * (1 - exp(-u * rod.getDRot() * time)) / u / rod.getDRot() / time; } break; - case Type::msd_harmonic_k_10: { + case Type::msd_harmonic_k_1: { + resetting = true; type_str = "msd_harmonic"; - target = 2.0/10.0*(1.0 - std::exp(-20 * time)); + target = 2 * (1.0 - std::exp(-2 * time)); + } + break; + case Type::msd_const: { + resetting = true; + type_str = "msd_const"; + target = -1.0; + } + break; + case Type::x: { + resetting = false; + type_str = ""; + target = -1.0; + } + break; + case Type::x_squared: { + resetting = false; + type_str = ""; + target = -1.0; } break; } @@ -121,3 +163,8 @@ std::ostream &operator<<(std::ostream &os, const Compute &compute) { os << "agg: " << compute.agg << " type: " << compute.type_str; return os; } + +void Compute::reset(const Rod2d &rod) { + time_step = 0; + start_rod = rod; +} diff --git a/src/Compute.h b/src/Compute.h index ae3213e..309b336 100644 --- a/src/Compute.h +++ b/src/Compute.h @@ -12,7 +12,7 @@ class Compute { public: - enum class Type { msd, oaf, empxx, empyy, msd_harmonic_k_10 }; + enum class Type { msd, oaf, empxx, empyy, msd_harmonic_k_1,msd_const, x,x_squared }; explicit Compute(Rod2d rod, Type t_type, size_t t_every, Simulation &sim); @@ -24,7 +24,7 @@ class Compute { void eval_empXX(const Rod2d &rod2D); - double getTarget() const; + [[nodiscard]] double getTarget() const; void eval_empYY(const Rod2d &rod2D); @@ -32,12 +32,14 @@ class Compute { [[nodiscard]] Type getType() const; - double getTime() const; + [[nodiscard]] double getTime() const; - double getDifference() const; + [[nodiscard]] double getDifference() const; friend std::ostream &operator<<(std::ostream &os, const Compute &compute); + void reset(const Rod2d &rod); + private: Rod2d start_rod; LiveAgg agg; @@ -47,6 +49,11 @@ private: double target; double time; std::string type_str; + + + void evalX(const Rod2d &rod2D); + + bool resetting; }; #endif // MYPROJECT_COMPUTE_H diff --git a/src/Integratoren2d.cpp b/src/Integratoren2d.cpp new file mode 100644 index 0000000..d0fdd3c --- /dev/null +++ b/src/Integratoren2d.cpp @@ -0,0 +1,205 @@ +// +// Created by jholder on 21.10.21. +// + +#include "Integratoren2d.h" +#include "vec_trafos.h" + +using Vector = Eigen::Vector2d; +using Matrix = Eigen::Matrix2d; + +Vector +force_euler(const Rod2d &rod2D, const Simulation &sim, const std::function &force) { + if (not force) { + return {0.0, 0.0}; + } + const auto &e = rod2D.getE(); + Vector F_lab = force(rod2D.getPos(), e); + Vector F_part = rotation_Matrix(e).inverse() * F_lab; + Vector F_trans = rotation_Matrix(e) * rod2D.getDiff_Sqrt() * F_part; + F_trans *= sim.getMDeltaT(); + return F_trans; +} + +Vector torque_euler(const Rod2d &rod2D, const Simulation &sim, const std::function &torque) { + if (not torque) { + return {0.0, 0.0}; + } + const auto &e = rod2D.getE(); + return torque(rod2D.getPos(), e) * rod2D.getDRot() * sim.getMDeltaT() * ortho(rod2D.getE()); +} + +Vector +rand_trans_euler(const Rod2d &rod2D, Simulation &sim) { + const auto std = sim.getSTD(); + const auto &e = rod2D.getE(); + Vector rand = {sim.getNorm(std), sim.getNorm(std)}; // gaussian noise + Vector trans_part = rod2D.getDiff_Sqrt() * rand; + return rotation_Matrix(e) * trans_part; +} + +Vector rand_rot_euler(const Rod2d &rod2D, Simulation &sim) { + const auto std = sim.getSTD(); + return sim.getNorm(std) * rod2D.getDRot_Sqrt() * ortho(rod2D.getE()); +} + +void Integratoren2d::Set1_Euler_f( + Rod2d &rod2D, Simulation &sim, + const std::function &force, + const std::function &torque) { + + // Gaussian - translation + auto trans = rand_trans_euler(rod2D, sim); + + // Force + auto F_trans = force_euler(rod2D, sim, force); + + // Gaussian rotation + auto rot = rand_rot_euler(rod2D, sim); + + // torque + auto M_rot = torque_euler(rod2D, sim, torque); + + Vector integrated_e = rod2D.getE() + rot + M_rot; + integrated_e.normalize(); + Vector integrated_pos = rod2D.getPos() + trans + F_trans; + + // apply + rod2D.setPos(integrated_pos); + rod2D.setE(integrated_e); +} + +void Integratoren2d::Set1_Euler( + Rod2d &rod2D, Simulation &sim) { + Set1_Euler_f(rod2D, sim, {}, {}); +} + + +void Integratoren2d::Set2_Heun_f( + Rod2d &rod2D, Simulation &sim, + const std::function &force, + const std::function &torque) { + //position + Vector pos_predictor = rod2D.getPos() + rand_trans_euler(rod2D, sim) + force_euler(rod2D, sim, force); + //rotation + Vector delta_e_predictor = rand_rot_euler(rod2D, sim) + torque_euler(rod2D, sim, torque); + Vector e_predictor = rod2D.getE() + delta_e_predictor; + e_predictor.normalize(); + + //make predicted Particle + Rod2d pred_rod = rod2D; + pred_rod.setPos(pos_predictor); + pred_rod.setE(e_predictor); + + //integrate euler for future + Vector delta_e_future = rand_rot_euler(rod2D, sim) + torque_euler(rod2D, sim, torque); + + //Heun integration + Vector e_integrated = + rod2D.getE() + 0.5 * (delta_e_predictor + delta_e_future); + // apply + rod2D.setPos(pos_predictor); + rod2D.setE(e_integrated.normalized()); +} + +void Integratoren2d::Set2_Heun( + Rod2d &rod2D, Simulation &sim) { + Set2_Heun_f(rod2D, sim, {}, {}); +} + +void Integratoren2d::Set3_Exact_f( + Rod2d &rod2D, Simulation &sim, + const std::function &force, + const std::function &torque) { + // Gaussian - translation + auto trans = rand_trans_euler(rod2D, sim); + + // Force + auto F_trans = force_euler(rod2D, sim, force); + + Vector integrated_pos = rod2D.getPos() + trans + F_trans; + + // Gaussian rotation + auto rot = rand_rot_euler(rod2D, sim); + + // torque + auto M_rot = torque_euler(rod2D, sim, torque); + + auto correction = + -0.5 * pow(rod2D.getDRot(), 2) * sim.getMDeltaT() * rod2D.getE(); + + Vector integrated_e = rod2D.getE() + rot + M_rot + correction; + integrated_e.normalize(); + // apply + rod2D.setPos(integrated_pos); + rod2D.setE(integrated_e); +} + +void Integratoren2d::Set3_Exact(Rod2d &rod2D, Simulation &sim) { + Set3_Exact_f(rod2D, sim, {}, {}); +} + +void Integratoren2d::Set4_BDAS_f( + Rod2d &rod2D, Simulation &sim, + const std::function &force, + const std::function &torque) { + // Gaussian - translation + auto trans = rand_trans_euler(rod2D, sim); + // Force + auto F_trans = force_euler(rod2D, sim, force); + Vector integrated_pos = rod2D.getPos() + trans + F_trans; + + // rotation + double rot = sim.getNorm(sim.getSTD()) * rod2D.getDRot_Sqrt(); + + double M_rot = torque ? torque(rod2D.getPos(), rod2D.getE()) * rod2D.getDRot() * sim.getMDeltaT() : 0.0; + + Vector integrated_e = Eigen::Rotation2D(rot + M_rot) * rod2D.getE(); + integrated_e.normalize(); + // apply + rod2D.setPos(integrated_pos); + rod2D.setE(integrated_e); +} + +void Integratoren2d::Set4_BDAS(Rod2d &rod2D, Simulation &sim) { + Set4_BDAS_f(rod2D, sim, {}, {}); +} + + +void Integratoren2d::Set5_MBD_f(Rod2d &rod2D, Simulation &sim, + const std::function &force, + const std::function &torque) { + + //position + Vector pos_predictor = rod2D.getPos() + rand_trans_euler(rod2D, sim) + force_euler(rod2D, sim, force); + //rotation + Vector delta_e_predictor = rand_rot_euler(rod2D, sim) + torque_euler(rod2D, sim, torque); + Vector e_predictor = rod2D.getE() + delta_e_predictor; + e_predictor.normalize(); + + //make predicted Particle + Rod2d pred_rod = rod2D; + pred_rod.setPos(pos_predictor); + pred_rod.setE(e_predictor); + + auto std = sim.getSTD(); + Vector e = rod2D.getE(); + + auto rot_Mat = [](const Matrix& mat,const Matrix& rotMat) { return rotMat * mat * rotMat.transpose(); }; + Eigen::Matrix2d Diff_combined = 0.5 * (rot_Mat(rod2D.getDiff(), rotation_Matrix(e)) + + rot_Mat(rod2D.getDiff(), rotation_Matrix(e_predictor))); + Eigen::Matrix2d D_combined_sqrt = Diff_combined.llt().matrixL(); + auto rand = Eigen::Vector2d(sim.getNorm(std), sim.getNorm(std)); + double rot = sim.getNorm(std) * rod2D.getDRot_Sqrt(); + Eigen::Vector2d e_integrated = e + 0.5 * (rot * ortho(e) + rot * ortho(e_predictor)); + + Eigen::Vector2d pos_integrated = rod2D.getPos() + D_combined_sqrt * rand + 0.5*(force_euler(rod2D,sim,force)+force_euler(pred_rod,sim,force)); + + //apply + rod2D.setPos(pos_integrated); + rod2D.setE(e_integrated.normalized()); +} + +void Integratoren2d::Set5_MBD(Rod2d &rod2D, Simulation &sim) { + Set5_MBD_f(rod2D, sim, {}, {}); +} \ No newline at end of file diff --git a/src/Integratoren2d.h b/src/Integratoren2d.h new file mode 100644 index 0000000..4142807 --- /dev/null +++ b/src/Integratoren2d.h @@ -0,0 +1,44 @@ +// +// Created by jholder on 21.10.21. +// + +#pragma once + +#include "Rod2d.hpp" +#include "Simulation.h" + +class Integratoren2d { + using Vector = Eigen::Vector2d; +public: + + static void Set1_Euler_f(Rod2d &rod2D, Simulation &sim, const std::function &force, + const std::function &torque); + + static void Set1_Euler(Rod2d &rod2D, Simulation &sim); + + static void Set2_Heun_f(Rod2d &rod2D, Simulation &sim, const std::function &force, + const std::function &torque); + + static void Set2_Heun(Rod2d &rod2D, Simulation &sim); + + static void + Set3_Exact_f(Rod2d &rod2D, Simulation &sim, + const std::function &force, + const std::function &torque); + + static void Set3_Exact(Rod2d &rod2D, Simulation &sim); + + static void + Set4_BDAS_f(Rod2d &rod2D, Simulation &sim, + const std::function &force, + const std::function &torque); + + static void Set4_BDAS(Rod2d &rod2D, Simulation &sim); + + static void + Set5_MBD_f(Rod2d &rod2D, Simulation &sim, + const std::function &force, + const std::function &torque); + + static void Set5_MBD(Rod2d &rod2D, Simulation &sim); +}; diff --git a/src/Rod2d.cpp b/src/Rod2d.cpp index cf29906..b803c66 100644 --- a/src/Rod2d.cpp +++ b/src/Rod2d.cpp @@ -31,11 +31,6 @@ Rod2d::Rod2d(double L) : m_pos({0, 0}), m_e({1, 0}) { sqrt(m_Diff_sqrt); } -void Rod2d::reset() { - m_pos = {0, 0}; - m_e = {1, 0}; -} - void Rod2d::setPos(const Eigen::Vector2d &Pos) { m_pos = Pos; } double Rod2d::getDRot() const { return m_D_rot; } diff --git a/src/Rod2d.hpp b/src/Rod2d.hpp index dfb8079..10bce4e 100644 --- a/src/Rod2d.hpp +++ b/src/Rod2d.hpp @@ -4,7 +4,6 @@ class Rod2d { public: explicit Rod2d(double L); - void reset(); [[nodiscard]] const Eigen::Matrix2d &getDiff() const; diff --git a/src/force_lambdas.h b/src/force_lambdas.h new file mode 100644 index 0000000..d805705 --- /dev/null +++ b/src/force_lambdas.h @@ -0,0 +1,25 @@ +// +// Created by jholder on 27.10.21. +// + +#pragma once +[[maybe_unused]] const static auto harmonic_Force = + [](const Eigen::Vector2d &pos, + const Eigen::Vector2d & /*torque*/) -> Eigen::Vector2d { + return -1.0 * pos; + }; +[[maybe_unused]] const static auto const_Force = + [](const Eigen::Vector2d & /*pos*/, + const Eigen::Vector2d & /*torque*/) -> Eigen::Vector2d { + return {1,1}; + }; +[[maybe_unused]] const static auto zero_Force = + [](const Eigen::Vector2d & /*pos*/, + const Eigen::Vector2d & /*torque*/) -> Eigen::Vector2d { + return {0.0, 0.0}; + }; +[[maybe_unused]] const static auto zero_Torque = [](const Eigen::Vector2d & /*pos*/, + const Eigen::Vector2d & /*torque*/) -> double { + return 0.0; +}; + diff --git a/src/Integratoren2d_force.cpp b/src/legacy/Integratoren2d_force.cpp similarity index 94% rename from src/Integratoren2d_force.cpp rename to src/legacy/Integratoren2d_force.cpp index 07d358b..142af06 100644 --- a/src/Integratoren2d_force.cpp +++ b/src/legacy/Integratoren2d_force.cpp @@ -8,14 +8,6 @@ using Vector = Eigen::Vector2d; using Matrix = Eigen::Matrix2d; -void Integratoren2d_force::Set0_deterministic(Rod2d &rod2D, - Simulation & /*sim*/) { - auto trans_lab = Vector({0, 0.01}); - rod2D.setPos(rod2D.getPos() + trans_lab); - Eigen::Rotation2D rot(0.1); - rod2D.setE(rot * rod2D.getE()); -} - const Vector Integratoren2d_force::unitVec = {1, 0}; void Integratoren2d_force::Set1_Euler( @@ -36,8 +28,7 @@ void Integratoren2d_force::Set1_Euler( F_trans *= sim.getMDeltaT(); // rotation - double rot = - sim.getNorm(std) * rod2D.getDRot_Sqrt(); // rotationsmatrix verwenden. + double rot = sim.getNorm(std) * rod2D.getDRot_Sqrt(); double M_rot = torque(rod2D.getPos(), e) * rod2D.getDRot() * sim.getMDeltaT(); Vector new_e = e + (rot + M_rot) * ortho(e); new_e.normalize(); diff --git a/src/Integratoren2d_force.h b/src/legacy/Integratoren2d_force.h similarity index 90% rename from src/Integratoren2d_force.h rename to src/legacy/Integratoren2d_force.h index a3bc9c2..5816954 100644 --- a/src/Integratoren2d_force.h +++ b/src/legacy/Integratoren2d_force.h @@ -37,9 +37,7 @@ class Integratoren2d_force { const std::function& force, const std::function& torque); - static void Set0_deterministic(Rod2d &rod2D, Simulation &sim); - - private: +private: static Eigen::Vector2d Heun_predictor_pos( const Rod2d &rod2D, Simulation &sim, const std::function @@ -48,4 +46,7 @@ class Integratoren2d_force { static Eigen::Vector2d Heun_predictor_rot( const Rod2d &rod2D, Simulation &sim, const std::function &torque); + + void Set1_Euler(Rod2d &rod2D, Simulation &sim, const std::function &force, + const std::function &torque); }; diff --git a/src/Integratoren2d_forceless.cpp b/src/legacy/Integratoren2d_forceless.cpp similarity index 100% rename from src/Integratoren2d_forceless.cpp rename to src/legacy/Integratoren2d_forceless.cpp diff --git a/src/Integratoren2d_forceless.h b/src/legacy/Integratoren2d_forceless.h similarity index 100% rename from src/Integratoren2d_forceless.h rename to src/legacy/Integratoren2d_forceless.h diff --git a/src/main.cpp b/src/main.cpp index 847a6c3..7dafd08 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -6,60 +6,114 @@ #include #include "Calculation.h" -#include "Integratoren2d_force.h" -#include "Integratoren2d_forceless.h" +#include "Integratoren2d.h" +#include "force_lambdas.h" -int main() { - constexpr int numStep = 10000000; - constexpr double k = 10; - [[maybe_unused]] auto harmonic_Force = - [](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}; - }; - [[maybe_unused]] auto zero_Torque = [](const Eigen::Vector2d & /*pos*/, - 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); +constexpr size_t SEED = 1234; +constexpr double stepSize = 0.01; +constexpr int n_computes = 100; +constexpr size_t num_runs = 5; +constexpr int numStep = 1000000; +void run_no_force(size_t integrator_index, size_t force_index) { + + + using type = std::pair; + std::vector vec({ + type(Integratoren2d::Set1_Euler_f, "euler"), + type(Integratoren2d::Set2_Heun_f, "heun"), + type(Integratoren2d::Set3_Exact_f, "exact"), + type(Integratoren2d::Set4_BDAS_f, "bdas"), + type(Integratoren2d::Set5_MBD_f, "mbd")}); + auto[integrator, name] = vec[integrator_index]; + + auto force = std::vector>({zero_Force, const_Force, + harmonic_Force})[force_index]; + auto force_name = std::vector({"zero_force_", "const_force_", "harmonic_force_"})[force_index]; + std::string folder = force_name + name + "_"; + { + 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 * 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(integrator, + computes, stepSize, SEED, force, zero_Torque, 1.0); + + + + euler.run(numStep); + std::ofstream msdFile(folder + "msd.dat"); + std::ofstream msd(folder + "msd.dat"); + std::ofstream oafFile(folder + "oaf.dat"); + std::ofstream empxxFile(folder + "empxx.dat"); + std::ofstream empyyFile(folder + "empyy.dat"); + + for (const auto &com: euler.getComputes()) { + if (com.getType() == MSD_VARIANT) { + 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) { + msdFile<< com.getTime() << ", " << com.getAgg().getMean() << ", " << com.getTarget() << std::endl; + } + } + } + { + std::vector> repeating_computes; + for (int i = 1; i < n_computes; ++i) { + repeating_computes.emplace_back(Compute::Type::x, i * 1); + repeating_computes.emplace_back(Compute::Type::x_squared, i * 1); + + } + Calculation calc_repeat(integrator, + repeating_computes, stepSize, SEED, force, zero_Torque, 1.0); + for (int i = 0; i - -#include "Integratoren2d_forceless.h" -#include "LiveAgg.hpp" - -TEST_CASE("Basic run of Calculation") { - Calculation euler(Integratoren2d_forceless::Set1_Euler, - {{Compute::Type::msd, 10}, {Compute::Type::oaf, 10}}, 1, - 1); - Calculation heun(Integratoren2d_forceless::Set2_Heun, - {{Compute::Type::msd, 10}, {Compute::Type::oaf, 10}}, 1, - 1); - Calculation exact(Integratoren2d_forceless::Set3_Exact, - {{Compute::Type::msd, 10}, {Compute::Type::oaf, 10}}, 1, - 1); - Calculation bdas(Integratoren2d_forceless::Set4_BDAS, - {{Compute::Type::msd, 10}, {Compute::Type::oaf, 10}}, 1, - 1); - SECTION("Euler") { - euler.run(100); - CHECK(euler.getRod().getE().norm() == Catch::Detail::Approx(1.0)); - CHECK(euler.getComputes()[0].getType() == Compute::Type::msd); - CHECK(euler.getComputes()[1].getType() == Compute::Type::oaf); - } - SECTION("Heun") { - heun.run(100); - CHECK(heun.getRod().getE().norm() == Catch::Detail::Approx(1.0)); - } - SECTION("Exact") { - exact.run(100); - CHECK(exact.getRod().getE().norm() == Catch::Detail::Approx(1.0)); - } - SECTION("Euler") { - bdas.run(100); - CHECK(bdas.getRod().getE().norm() == Catch::Detail::Approx(1.0)); - } - - SECTION("MSD") { - euler.run(100); - CHECK(euler.getComputes()[0].getAgg().getNumPoints() == 10); - CHECK(euler.getComputes()[1].getAgg().getNumPoints() == 10); - } - SECTION("Deterministic") { - Calculation determ(Integratoren2d_forceless::Set0_deterministic, - {{Compute::Type::msd, 1}, {Compute::Type::oaf, 1}}, - 1, 1); - determ.run(10); - auto time = 1; - - auto targetMSD = pow(0.01, 2) * time; - auto targetOAF = cos(0.1); - - CHECK(determ.getComputes()[0].getAgg().getMean() == - Catch::Detail::Approx(targetMSD)); - CHECK(determ.getComputes()[1].getAgg().getMean() == - Catch::Detail::Approx(targetOAF)); - } -} diff --git a/test/test_Compute.cpp b/test/test_Compute.cpp index 6e06456..5437766 100644 --- a/test/test_Compute.cpp +++ b/test/test_Compute.cpp @@ -5,7 +5,6 @@ #include "Compute.h" #include -#include "Integratoren2d_forceless.h" TEST_CASE("Compute") { Rod2d rod(1.0); diff --git a/test/test_Integratoren.cpp b/test/test_Integratoren.cpp index e1bc7c3..ef59e70 100644 --- a/test/test_Integratoren.cpp +++ b/test/test_Integratoren.cpp @@ -1,357 +1,62 @@ + // // Created by jholder on 24.10.21. // #include "Calculation.h" #include "Compute.h" -#include "Integratoren2d_force.h" -#include "Integratoren2d_forceless.h" +#include "Integratoren2d.h" +#include "force_lambdas.h" #include +#include -TEST_CASE("Forceless Integratoren") { +TEST_CASE("Integrator Check") { const size_t SEED = Catch::rngSeed(); - const double length = GENERATE(1.0,1.4,2.0); + const double length = GENERATE(1.0, 1.4, 2.0); constexpr int steps = 10000; constexpr double delta = 0.1; + using type = std::tuple, Calculation::inte_force_type, std::string>; + auto[integrator, integrator_force, name] = GENERATE( + type(Integratoren2d::Set1_Euler, Integratoren2d::Set1_Euler_f, "Euler"), + type(Integratoren2d::Set2_Heun, Integratoren2d::Set2_Heun_f, "Heun"), + type(Integratoren2d::Set3_Exact, Integratoren2d::Set3_Exact_f, "Exact"), + type(Integratoren2d::Set4_BDAS, Integratoren2d::Set4_BDAS_f, "BDAS"), + type(Integratoren2d::Set5_MBD, Integratoren2d::Set5_MBD_f, "MBD")); - 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, + Calculation calc(integrator, {{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") { + 0.01, SEED, length); + calc.run(steps); + SECTION("ForceLess") { + CAPTURE(name); CAPTURE(length); size_t i = 0; - for (auto &c: euler_zero.getComputes()) { + for (const auto &c: calc.getComputes()) { CAPTURE(c); - CHECK(c.getDifference() <= delta*c.getAgg().getMean()); + 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") { + Calculation calc_f(integrator_force, + {{Compute::Type::msd, 1}, + {Compute::Type::oaf, 1}, + {Compute::Type::empxx, 1}, + {Compute::Type::empyy, 1}}, + 0.01, SEED, zero_Force, zero_Torque, length); + calc_f.run(steps); + + SECTION("ForceLess") { CAPTURE(length); - size_t i = 0; - for (auto &c: exact_zero.getComputes()) { + CAPTURE(name); + for (const auto &c: calc_f.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; + CHECK(c.getDifference() <= delta * c.getAgg().getMean()); } } } - -/* - - - - - - - SECTION("Euler") { - - { - - } - - Calculation euler_force( - Integratoren2d_force::Set1_Euler, - {{Compute::Type::msd, 1}, {Compute::Type::oaf, 1}}, 0.01, SEED, - force, torque); - { - euler_force.run(10000); - for (auto &c : euler_force.getComputes()) { - CHECK(c.getDifference() <= 0.005); - } - } - CHECK(euler.getComputes()[0].getAgg().getMean() == - euler_force.getComputes()[0].getAgg().getMean()); - CHECK(euler.getComputes()[1].getAgg().getMean() == - euler_force.getComputes()[1].getAgg().getMean()); - } - SECTION("Heun") { - 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); - { - heun.run(10000); - for (auto &c : heun.getComputes()) { - CHECK(c.getDifference() <= 0.005); - } - } - - Calculation heun_force(Integratoren2d_force::Set2_Heun, - {{Compute::Type::msd, 1}, - {Compute::Type::oaf, 1}, - {Compute::Type::empxx, 1}, - {Compute::Type::empyy, 1}}, - 0.01, SEED, force, torque); - { - heun_force.run(10000); - for (auto &c : heun_force.getComputes()) { - CHECK(c.getDifference() <= 0.005); - } - } - CHECK(heun.getComputes()[0].getAgg().getMean() == - heun_force.getComputes()[0].getAgg().getMean()); - CHECK(heun.getComputes()[1].getAgg().getMean() == - heun_force.getComputes()[1].getAgg().getMean()); - } - SECTION("Exact") { - 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); - { - exact.run(10000); - for (auto &c : exact.getComputes()) { - CHECK(c.getDifference() <= 0.005); - } - } - - Calculation exact_force(Integratoren2d_force::Set3_Exact, - {{Compute::Type::msd, 1}, - {Compute::Type::oaf, 1}, - {Compute::Type::empxx, 1}, - {Compute::Type::empyy, 1}}, - 0.01, SEED, force, torque); - { - exact_force.run(10000); - for (auto &c : exact_force.getComputes()) { - CHECK(c.getDifference() <= 0.005); - } - } - CHECK(exact.getComputes()[0].getAgg().getMean() == - Approx(exact_force.getComputes()[0].getAgg().getMean()) - .epsilon(1e-10)); - CHECK(exact.getComputes()[1].getAgg().getMean() == - exact_force.getComputes()[1].getAgg().getMean()); - } - SECTION("BDAS") { - 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); - { - bdas.run(10000); - for (auto &c : bdas.getComputes()) { - CHECK(c.getDifference() <= 0.005); - } - } - - Calculation bdas_force(Integratoren2d_force::Set4_BDAS, - {{Compute::Type::msd, 1}, - {Compute::Type::oaf, 1}, - {Compute::Type::empxx, 1}, - {Compute::Type::empyy, 1}}, - 0.01, SEED, force, torque); - { - bdas_force.run(10000); - for (auto &c : bdas_force.getComputes()) { - CHECK(c.getDifference() <= 0.005); - } - } - CHECK(bdas.getComputes()[0].getAgg().getMean() == - Approx(bdas_force.getComputes()[0].getAgg().getMean()) - .epsilon(10e-10)); - CHECK(bdas.getComputes()[1].getAgg().getMean() == - 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