From 29779754e9faeeea4b2a2f6e6d957cb07a3da2be Mon Sep 17 00:00:00 2001 From: Jacob Date: Sun, 24 Oct 2021 17:41:05 +0200 Subject: [PATCH] Wokring Without Heunwith force and MBD --- README.md | 32 ++++++ src/CMakeLists.txt | 4 +- src/Calculation.cpp | 2 +- src/Compute.cpp | 30 +++++- src/Compute.h | 10 +- src/Integratoren2d_force.cpp | 171 +++++++++++++++++++++++++++++++ src/Integratoren2d_force.h | 48 +++++++++ src/Integratoren2d_forceless.cpp | 12 +-- src/Integratoren2d_forceless.h | 6 +- src/main.cpp | 57 ++++++----- test/CMakeLists.txt | 2 +- test/test_Calculation.cpp | 8 +- test/test_Compute.cpp | 5 +- test/test_Integratoren.cpp | 76 ++++++++++++++ 14 files changed, 417 insertions(+), 46 deletions(-) create mode 100644 README.md create mode 100644 src/Integratoren2d_force.cpp create mode 100644 src/Integratoren2d_force.h create mode 100644 test/test_Integratoren.cpp diff --git a/README.md b/README.md new file mode 100644 index 0000000..7576fb4 --- /dev/null +++ b/README.md @@ -0,0 +1,32 @@ +Brownian Motion Integrator Implementation +=============== + +This is a sample Implementation for different BD Algorithm + +1. Classic Euler algorithm +2. Heun algorithm +3. Exact solution +4. BDAS - RotationMatrix +5. MBD + +The current implementation allows only the use for a 2D case + +build +----- + +To build this project you need to install +- conan ```pip3 install --user conan``` +- CMake + +To generate the build files +```shell +cmake -S . -B build +``` +For build configuration use +```shell +ccmake build +``` +and to build the program +```shell +cmake --build build +``` diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 9e77618..3c3038a 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -4,8 +4,8 @@ find_package(Eigen3) # Generic test that uses conan libs add_library(integratoren SHARED LiveAgg.cpp Rod2d.cpp Simulation.cpp - Simulation.h Integratoren2d_forceless.cpp - Calculation.cpp Calculation.h Compute.cpp Compute.h) + Simulation.h Integratoren2d_forceless.cpp Integratoren2d_force.cpp + Calculation.cpp Compute.cpp) target_include_directories(integratoren PUBLIC .) target_link_libraries( integratoren diff --git a/src/Calculation.cpp b/src/Calculation.cpp index dc107e4..7689f6c 100644 --- a/src/Calculation.cpp +++ b/src/Calculation.cpp @@ -24,7 +24,7 @@ Calculation::Calculation(std::function t_integrator : sim(deltaT, seed), m_integrator( std::move(t_integrator)) { for (const auto &pair: t_computes) { - computes.emplace_back(rod, pair.first, pair.second); + computes.emplace_back(Compute(rod, pair.first, pair.second,sim)); } } diff --git a/src/Compute.cpp b/src/Compute.cpp index 28e6aab..85b609c 100644 --- a/src/Compute.cpp +++ b/src/Compute.cpp @@ -7,6 +7,7 @@ #include #include #include "Rod2d.hpp" +#include "Simulation.h" void Compute::evalMSD(const Rod2d &rod2D) { const auto &new_pos = rod2D.getPos(); @@ -32,7 +33,7 @@ void Compute::eval_empYY(const Rod2d & /*rod2D*/) { void Compute::eval(const Rod2d &rod2D) { time_step++; - if (time_step % every == 0){ + if (time_step % every == 0) { switch (type) { case Type::msd: evalMSD(rod2D); @@ -52,8 +53,23 @@ void Compute::eval(const Rod2d &rod2D) { } -Compute::Compute(Rod2d rod, Type t_type, size_t t_every) - : start_rod(std::move(rod)), agg({}), every(t_every), time_step(0), type(t_type) { +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); + switch (type) { + case Type::msd: + target = 4.0 * 0.5 * (rod.getDiff().trace()) * time; + break; + case Type::oaf: + target = std::exp(-rod.getDRot() * time); + break; + case Type::empxx: + target = 0; + break; + case Type::empyy: + target = 0; + break; + } } const LiveAgg &Compute::getAgg() const { @@ -67,3 +83,11 @@ Compute::Type Compute::getType() const { 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 diff --git a/src/Compute.h b/src/Compute.h index bdd24f2..ee1d1ac 100644 --- a/src/Compute.h +++ b/src/Compute.h @@ -8,6 +8,7 @@ #include "LiveAgg.hpp" #include "Rod2d.hpp" +#include "Simulation.h" class Compute { @@ -16,7 +17,7 @@ public: msd, oaf, empxx, empyy }; - explicit Compute(Rod2d rod, Type t_type, size_t t_every); + explicit Compute(Rod2d rod, Type t_type, size_t t_every, Simulation &sim); void eval(const Rod2d &rod2D); @@ -26,6 +27,8 @@ public: void eval_empXX(const Rod2d &rod2D); + double getTarget() const; + void eval_empYY(const Rod2d &rod2D); [[nodiscard]] const LiveAgg &getAgg() const; @@ -34,12 +37,17 @@ public: double getTime() const; + double getDifference() const; + + private: Rod2d start_rod; LiveAgg agg; size_t every; size_t time_step; Type type; + double target; + }; diff --git a/src/Integratoren2d_force.cpp b/src/Integratoren2d_force.cpp new file mode 100644 index 0000000..d92fd23 --- /dev/null +++ b/src/Integratoren2d_force.cpp @@ -0,0 +1,171 @@ +// +// Created by jholder on 21.10.21. +// + +#include "Integratoren2d_force.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; +} + +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::Set5_MBD(Rod2d &/*rod2D*/, Simulation &/*sim*/) { +} + +void Integratoren2d_force::Set1_Euler(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; + + //Force + Vector F_lab = force(rod2D.getPos(), rod2D.getE()); + Vector F_part = rod2D.getE_Base_matrix().inverse() * F_lab; + Vector F_trans = rod2D.getDiff_Sqrt() * F_part; + F_trans *= sim.getMDeltaT() / kBT; + + + //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(); + Vector new_e = e + (rot + M_rot) * ortho(e); + new_e.normalize(); + //apply + rod2D.setPos(rod2D.getPos() + trans_lab + F_trans); + rod2D.setE(new_e); +} + +Vector Integratoren2d_force::Heun_predictor_pos(const Rod2d &rod2D, Simulation &sim, + const std::function &force) { + auto standard_dev = sim.getSTD(); + Vector rand_pred = {sim.getNorm(standard_dev), sim.getNorm(standard_dev)}; //gaussian noise + 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 F_pred_lab = force(rod2D.getPos(), rod2D.getE()); + Vector F_pred_part = rod2D.getE_Base_matrix().inverse() * F_pred_lab; + Vector F_pred_trans = rod2D.getDiff_Sqrt() * F_pred_part; + F_pred_trans *= sim.getMDeltaT() / kBT; + + return F_pred_trans + trans_pred_lab; +} + +Vector Integratoren2d_force::Heun_predictor_rot(const Rod2d &rod2D, Simulation &sim, + const std::function &torque) { + auto std = sim.getSTD(); + double rot_predict = sim.getNorm(std) * rod2D.getDRot_Sqrt();// rotationsmatrix verwenden. + Vector e = rod2D.getE(); + double M_predict_rot = torque(rod2D.getPos(), rod2D.getE()) * rod2D.getDRot() / kBT * sim.getMDeltaT(); + Vector e_change_predict = (rot_predict + M_predict_rot) * ortho(e); + + return e_change_predict; +} + + +void Integratoren2d_force::Set2_Heun(Rod2d &rod2D, Simulation &sim, + const std::function &force, + const std::function &torque) { + + Vector delta_pos_predictor = Heun_predictor_pos(rod2D, sim, force); + 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(); + + + Rod2d pred_rod = rod2D; + pred_rod.setPos(pos_predictor); + pred_rod.setE(e_predict); + + 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.setE(e_integrated.normalized()); +} + +void Integratoren2d_force::Set3_Exact(Rod2d &rod2D, Simulation &sim, + std::function force, + 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 F_lab = force(rod2D.getPos(), rod2D.getE()); + Vector F_part = rod2D.getE_Base_matrix().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; + + //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(); + auto correction = -0.5 * pow(rod2D.getDRot(), 2) / pow(kBT, 2) * sim.getMDeltaT(); + Vector new_e = e + (rot + M_rot) * ortho(e) + correction * e; + new_e.normalize(); + //apply + rod2D.setPos(rod2D.getPos() + trans_lab + F_trans); + rod2D.setE(new_e); +} + +void Integratoren2d_force::Set4_BDAS(Rod2d &rod2D, Simulation &sim, + std::function force, + std::function /*torque*/) { + auto std = sim.getSTD(); + //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 F_lab = force(rod2D.getPos(), rod2D.getE()); + Vector F_part = rotMat.inverse() * F_lab; + + Vector F_trans = rod2D.getDiff_Sqrt() * F_part; + F_trans *= sim.getMDeltaT() / kBT; + + + auto rot = Eigen::Rotation2D(sim.getNorm(std) * rod2D.getDRot_Sqrt()); + Eigen::Rotation2Dd rotation(rot); + auto e_new = (rotation.toRotationMatrix() * rod2D.getE()).normalized(); + + // 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); +} + + diff --git a/src/Integratoren2d_force.h b/src/Integratoren2d_force.h new file mode 100644 index 0000000..db22fba --- /dev/null +++ b/src/Integratoren2d_force.h @@ -0,0 +1,48 @@ +// +// Created by jholder on 21.10.21. +// + +#pragma once + +#include "Rod2d.hpp" +#include "Simulation.h" + +class Integratoren2d_force { +public: + static void Set1_Euler(Rod2d &rod2D, Simulation &sim, + const std::function &force, + const std::function &torque); + + + static void Set2_Heun(Rod2d &rod2D, Simulation &sim, + const std::function &force, + const std::function &torque); + + static void Set3_Exact(Rod2d &rod2D, Simulation &sim, + std::function force, + std::function torque); + + static void Set4_BDAS(Rod2d &rod2D, Simulation &sim, + std::function force, + std::function torque); + + static const Eigen::Vector2d unitVec; + + static void Set5_MBD(Rod2d &rod2D, Simulation &sim); + + 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 &force); + + static Eigen::Vector2d + Heun_predictor_rot(const Rod2d &rod2D, Simulation &sim, + const std::function &torque); +}; + diff --git a/src/Integratoren2d_forceless.cpp b/src/Integratoren2d_forceless.cpp index 520785b..b77a215 100644 --- a/src/Integratoren2d_forceless.cpp +++ b/src/Integratoren2d_forceless.cpp @@ -4,17 +4,10 @@ #include "Integratoren2d_forceless.h" -Eigen::Vector2d ortho(Eigen::Vector2d e) { +Eigen::Vector2d Integratoren2d_forceless::ortho(Eigen::Vector2d e) { return Eigen::Vector2d{-e[1], e[0]}; } -Eigen::Matrix2d mat(Eigen::Vector2d e) { - Eigen::Matrix2d matrix; - matrix << e, ortho(e); - - return matrix; -} - void Integratoren2d_forceless::Set0_deterministic(Rod2d &rod2D, Simulation & /*sim*/) { auto trans_lab = Eigen::Vector2d({0, 0.01}); rod2D.setPos(rod2D.getPos() + trans_lab); @@ -57,6 +50,9 @@ void Integratoren2d_forceless::Set2_Heun(Rod2d &rod2D, Simulation &sim) { 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.setE(e_integrated.normalized()); diff --git a/src/Integratoren2d_forceless.h b/src/Integratoren2d_forceless.h index 5cf70b1..0033dd3 100644 --- a/src/Integratoren2d_forceless.h +++ b/src/Integratoren2d_forceless.h @@ -4,6 +4,7 @@ #ifndef MYPROJECT_INTEGRATOREN2D_FORCELESS_H #define MYPROJECT_INTEGRATOREN2D_FORCELESS_H + #include "Rod2d.hpp" #include "Simulation.h" @@ -19,11 +20,14 @@ public: static const Eigen::Vector2d unitVec; - static Eigen::Matrix2d e_2_matrix(Eigen::Vector2d m_e) ; + 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/main.cpp b/src/main.cpp index a070fda..8984f7d 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -5,36 +5,47 @@ #include #include "Calculation.h" #include "Integratoren2d_forceless.h" +#include "Integratoren2d_force.h" int main() { - constexpr int numStep = 1000000; - Calculation euler(Integratoren2d_forceless::Set4_BDAS, {{Compute::Type::msd, 1}, - {Compute::Type::msd, 2}, - {Compute::Type::msd, 4}, - {Compute::Type::msd, 8}, - {Compute::Type::msd, 16}, - {Compute::Type::oaf, 1}, - {Compute::Type::oaf, 2}, - {Compute::Type::oaf, 4}, - {Compute::Type::oaf, 8}, - {Compute::Type::oaf, 16}}, 0.01, 12345); + constexpr int numStep = 1000; + auto zero_Force = [](const Eigen::Vector2d & /*pos*/, + const Eigen::Vector2d & /*torque*/) -> Eigen::Vector2d { return {0.0, 0.0}; }; + 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::Set2_Heun(rod, sim, zero_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::oaf, 1}, + {Compute::Type::oaf, 2}, + {Compute::Type::oaf, 4}, + {Compute::Type::oaf, 8}, + {Compute::Type::oaf, 16}}, 0.01, 12345); euler.run(numStep); for (const auto &com: euler.getComputes()) { - if(com.getType()!=Compute::Type::msd) + if (com.getType() != Compute::Type::msd) continue; - auto time = euler.getSim().getMDeltaT() * com.getTime(); - auto targetMSD = 4.0 * 0.5 * (euler.getRod().getDiff().trace()) * time; - std::cout << "MSD: " << com.getAgg().getMean() - << " +- " - << com.getAgg().getSEM() - << " " << targetMSD << std::endl; + std::cout << "MSD: " << com.getDifference() << " " << com.getAgg().getMean() + << " <-> " << com.getTarget() << std::endl; } for (const auto &com: euler.getComputes()) { - if(com.getType()!=Compute::Type::oaf) + if (com.getType() != Compute::Type::oaf) continue; - auto time = euler.getSim().getMDeltaT() * com.getTime(); - auto targetOAF = std::exp(-euler.getRod().getDRot() * time); - std::cout <<"OAF: " << com.getAgg().getMean() - << " +- " << com.getAgg().getSEM() << " " << targetOAF << std::endl; + std::cout << "OAF: " << com.getDifference() << " " << com.getAgg().getMean() + << " <-> "<< com.getTarget() << std::endl; } } \ No newline at end of file diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 2d18e0c..4f528bb 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -11,7 +11,7 @@ target_link_libraries(catch_main PRIVATE project_options) set(TEST_FILES test_LiveAgg.cpp test_Rod2d.cpp - test_Calculation.cpp test_Compute.cpp test_Simulation.cpp) + test_Calculation.cpp test_Compute.cpp test_Simulation.cpp test_Integratoren.cpp) add_executable(tests ${TEST_FILES}) target_link_libraries(tests PRIVATE project_warnings project_options catch_main Eigen3::Eigen3 integratoren) diff --git a/test/test_Calculation.cpp b/test/test_Calculation.cpp index 8c24aac..5679651 100644 --- a/test/test_Calculation.cpp +++ b/test/test_Calculation.cpp @@ -9,13 +9,13 @@ TEST_CASE("Basic run of Calculation") { Calculation euler(Integratoren2d_forceless::Set1_Euler, {{Compute::Type::msd, 10}, - {Compute::Type::oaf, 10}}, 1, 1); + {Compute::Type::oaf, 10}}, 1, 1); Calculation heun(Integratoren2d_forceless::Set2_Heun, {{Compute::Type::msd, 10}, - {Compute::Type::oaf, 10}}, 1, 1); + {Compute::Type::oaf, 10}}, 1, 1); Calculation exact(Integratoren2d_forceless::Set3_Exact, {{Compute::Type::msd, 10}, - {Compute::Type::oaf, 10}}, 1, 1); + {Compute::Type::oaf, 10}}, 1, 1); Calculation bdas(Integratoren2d_forceless::Set4_BDAS, {{Compute::Type::msd, 10}, - {Compute::Type::oaf, 10}}, 1, 1); + {Compute::Type::oaf, 10}}, 1, 1); SECTION("Euler") { euler.run(100); CHECK(euler.getRod().getE().norm() == Catch::Detail::Approx(1.0)); diff --git a/test/test_Compute.cpp b/test/test_Compute.cpp index 12a951e..fbcceb3 100644 --- a/test/test_Compute.cpp +++ b/test/test_Compute.cpp @@ -7,8 +7,9 @@ #include "Integratoren2d_forceless.h" TEST_CASE("Compute") { - auto rod = Rod2d(1.0); - auto com = Compute(rod, Compute::Type::msd, 10); + Rod2d rod(1.0); + Simulation sim(0.1,1); + auto com = Compute(rod, Compute::Type::msd, 10, sim); SECTION("Mean of same values") { for (int i = 0; i < 100; ++i) { com.eval(rod); diff --git a/test/test_Integratoren.cpp b/test/test_Integratoren.cpp new file mode 100644 index 0000000..d71d141 --- /dev/null +++ b/test/test_Integratoren.cpp @@ -0,0 +1,76 @@ +// +// Created by jholder on 24.10.21. +// +#include +#include +#include +#include +#include + +TEST_CASE("Integratoren") { + auto zero_Force = [](const Eigen::Vector2d & /*pos*/, + const Eigen::Vector2d & /*torque*/) -> Eigen::Vector2d { return {0.0, 0.0}; }; + 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::Set2_Heun(rod, sim, zero_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); + }; + SECTION("Euler") { + Calculation euler(Integratoren2d_forceless::Set1_Euler, {{Compute::Type::msd, 1}, + {Compute::Type::oaf, 1}}, 0.01, 12345); + euler.run(10000); + for (auto &c: euler.getComputes()) { CHECK(c.getDifference() <= 0.005); } + Calculation euler_force(euler_Zero, {{Compute::Type::msd, 1}, + {Compute::Type::oaf, 1}}, 0.01, 12345); + euler_force.run(10000); + for (auto &c: euler_force.getComputes()) { CHECK(c.getDifference() <= 0.005); } + }SECTION("Heun") { + { + Calculation heun(Integratoren2d_forceless::Set2_Heun, {{Compute::Type::msd, 1}, + {Compute::Type::oaf, 1}}, 0.01, 12345); + heun.run(10000); + for (auto &c: heun.getComputes()) { CHECK(c.getDifference() <= 0.005); } + } + { + Calculation heun_force(heun_Zero, {{Compute::Type::msd, 1}, + {Compute::Type::oaf, 1}}, 0.01, 12345); + heun_force.run(10000); + for (auto &c: heun_force.getComputes()) { CHECK(c.getDifference() <= 0.005); } + } + }SECTION("Exact") { + { + Calculation exact(Integratoren2d_forceless::Set3_Exact, {{Compute::Type::msd, 1}, + {Compute::Type::oaf, 1}}, 0.01, 12345); + exact.run(10000); + for (auto &c: exact.getComputes()) { CHECK(c.getDifference() <= 0.005); } + } + { + Calculation exact_force(exact_Zero, {{Compute::Type::msd, 1}, + {Compute::Type::oaf, 1}}, 0.01, 12345); + exact_force.run(10000); + for (auto &c: exact_force.getComputes()) { CHECK(c.getDifference() <= 0.005); } + } + }SECTION("BDAS") { + { + Calculation bdas(Integratoren2d_forceless::Set4_BDAS, {{Compute::Type::msd, 1}, + {Compute::Type::oaf, 1}}, 0.01, 12345); + bdas.run(10000); + for (auto &c: bdas.getComputes()) { CHECK(c.getDifference() <= 0.005); } + } + { + Calculation bdas_zero(bdas_Zero, {{Compute::Type::msd, 1}, + {Compute::Type::oaf, 1}}, 0.01, 12345); + bdas_zero.run(10000); + for (auto &c: bdas_zero.getComputes()) { CHECK(c.getDifference() <= 0.005); } + } + } +}