Wokring Without Heunwith force and MBD

2021-10-24 17:41:05 +02:00 · 2021-10-24 17:41:05 +02:00 · 29779754e9
commit 29779754e9
parent e63239ed47
14 changed files with 417 additions and 46 deletions
--- a/README.md
+++ 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
+```
--- 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
--- a/src/Calculation.cpp
+++ b/src/Calculation.cpp
@ -24,7 +24,7 @@ Calculation::Calculation(std::function<void(Rod2d &, Simulation &)> 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));
    }
 }

--- a/src/Compute.cpp
+++ b/src/Compute.cpp
@ -7,6 +7,7 @@
 #include <utility>
 #include <iostream>
 #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<double>(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<double>(every);
 }
+
+double Compute::getTarget() const {
+    return target;
+}
+
+double Compute::getDifference() const {
+    return abs(agg.getMean() -target);
+}
--- 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;
+
 };


--- a/src/Integratoren2d_force.cpp
+++ 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<double> 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<Vector(Vector, Vector)> &force,
+                                      const std::function<double(Vector, Vector)> &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<Vector(Vector, Vector)> &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<double(Vector, Vector)> &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<Vector(Vector, Vector)> &force,
+                                     const std::function<double(Vector, Vector)> &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<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)> force,
+                                      std::function<double(Eigen::Vector2d, Eigen::Vector2d)> 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<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)> force,
+                                     std::function<double(Eigen::Vector2d, Eigen::Vector2d)> /*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<double>(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);
+}
+
+
--- a/src/Integratoren2d_force.h
+++ 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<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)> &force,
+                           const std::function<double(Eigen::Vector2d, Eigen::Vector2d)> &torque);
+
+
+    static void Set2_Heun(Rod2d &rod2D, Simulation &sim,
+                          const std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)> &force,
+                          const std::function<double(Eigen::Vector2d, Eigen::Vector2d)> &torque);
+
+    static void Set3_Exact(Rod2d &rod2D, Simulation &sim,
+                           std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)> force,
+                           std::function<double(Eigen::Vector2d, Eigen::Vector2d)> torque);
+
+    static void Set4_BDAS(Rod2d &rod2D, Simulation &sim,
+                          std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)> force,
+                          std::function<double(Eigen::Vector2d, Eigen::Vector2d)> 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<double, 2, 1> e);
+
+    static Eigen::Vector2d Heun_predictor_pos(const Rod2d &rod2D, Simulation &sim,
+                                              const std::function<Eigen::Vector2d(Eigen::Vector2d,
+                                                                                  Eigen::Vector2d)> &force);
+
+    static Eigen::Vector2d
+    Heun_predictor_rot(const Rod2d &rod2D, Simulation &sim,
+                       const std::function<double(Eigen::Vector2d, Eigen::Vector2d)> &torque);
+};
+
--- 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());
--- 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);
 };


--- a/src/main.cpp
+++ b/src/main.cpp
@ -5,36 +5,47 @@
 #include <iostream>
 #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;
    }
 }
--- 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)
--- 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));
--- 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);
--- a/test/test_Integratoren.cpp
+++ b/test/test_Integratoren.cpp
@ -0,0 +1,76 @@
+//
+// Created by jholder on 24.10.21.
+//
+#include <catch2/catch.hpp>
+#include <Integratoren2d_forceless.h>
+#include <Calculation.h>
+#include <Compute.h>
+#include <Integratoren2d_force.h>
+
+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); }
+        }
+    }
+}