Ahhhhhh

Auswertung/harmonic.py

@@ -0,0 +1,40 @@
+import glob
+
+import matplotlib.pyplot as plt
+import pandas as pd
+import seaborn as sb
+import numpy as np
+
+
+def main():
+    time = np.arange(0.0, 4.0, 0.01)
+    xstart  = 1.0
+    ystart = 1.0
+    x_mean = xstart * np.exp(-time)
+    y_mean = ystart * np.exp(-time)
+    x_sq_mean = xstart*xstart * np.exp(-2*time) + (1 - np.exp(-2*time))
+    y_sq_mean = ystart*ystart * np.exp(-2 * time) + (1 - np.exp(-2 * time))
+
+    x = pd.read_csv("./data/out/harmonic_force_euler_L0_x.dat")
+    x_sq = pd.read_csv("./data/out/harmonic_force_euler_L0_x_squared.dat")
+    msd = pd.read_csv("./data/out/harmonic_force_euler_L0_msd.dat")
+
+    msd_mean =  (1-np.exp(-2*time))
+    plt.plot(time, x_mean, label="a <x>")
+    plt.plot(x["time"], x["val"], label="<x>")
+
+    plt.plot(time, x_sq_mean, label="a <x²>")
+    plt.plot(x_sq["time"], x_sq["val"], label="<x²>")
+
+    plt.plot(time, 2*msd_mean,label="a msd")
+
+    plt.plot(msd["time"], msd["val"], label="msd")
+
+    #plt.plot(x["time"], x_sq["val"] - np.power(x["val"], 2.0))
+    plt.plot(x["time"], 2.0*(x_sq["val"] - np.power(x["val"], 2.0)), label="<x²>-<x>²")
+    plt.legend(loc = 1)
+    plt.show()
+
+
+if __name__ == "__main__":
+    main()

Auswertung/main.py

@@ -2,30 +2,30 @@ import glob
 
 import matplotlib.pyplot as plt
 import pandas as pd
-import seaborn as sb
+import numpy as np
 
 
 def make_figure(length_index, length, forces, integrators, computes):
-    sb.set_theme()
     fig, axs = plt.subplots(nrows=len(forces), ncols=len(computes), figsize=(40, 24))
     fig.suptitle(length, fontsize=30)
     for index_compute, compute in enumerate(computes):
         for index_forces, force in enumerate(forces):
-
+            print("Reading file")
             legend = []
             plt.sca(axs[index_forces][index_compute])
             for index_integrator, integrator in enumerate(integrators):
                 for filename in glob.glob(f'data/out/*{force}*{integrator}*L{length_index}*_{compute}.dat'):
                     file = pd.read_csv(filename, delimiter=',')
-                    sb.lineplot(x="time", y="val", data=file)
+                    file.dropna("columns")
+                    ax = plt.plot(file["time"], file["val"])
+                    #ax.fill_between(file["time"], y1=file["val"] - file["std"], y2=file["val"] + file["std"], alpha=.5) 
                     plt.ylabel(compute)
                     legend.append(integrator)
-                    if index_integrator == 0:
+                    if index_integrator == 0 and not np.any(file["target"].isnull()):
-                        sb.lineplot(x="time", y="target", data=file)
+                        plt.plot(file["time"], file["target"])
                         legend.append("analytisch")
 
             plt.legend(title='iterators', labels=legend)
-            # plt.title("OAF without a force")
     pad = 5  # in points
 
     for ax, col in zip(axs[0], computes):
@@ -39,6 +39,7 @@ def make_figure(length_index, length, forces, integrators, computes):
                     size='large', ha='right', va='center', fontsize=25)
 
     fig.tight_layout()
+    fig.savefig(f"fig{length_index}.png")
 
 
 def main():
@@ -47,7 +48,6 @@ def main():
     computes = ["msd", "oaf", "empxx", "empyy", "x", "x_squared"]
     for l_i, l in enumerate(["sphere", "L = 1.5", "L = 2.0"]):
         make_figure(length=l, length_index=l_i, forces=forces, computes=computes, integrators=integrators)
-    plt.show()
 
 
 if __name__ == "__main__":

Auswertung/requirements.txt

@@ -1,4 +1,4 @@
-numpy~=1.21.3
+numpy~=1.19.5
 seaborn~=0.11.2
-matplotlib~=3.4.3
+matplotlib~=3.3.4
-pandas~=1.3.4
+pandas~=1.1.5

C++/src/Compute.cpp

@@ -13,19 +13,20 @@
 void Compute::evalMSD(const Rod2d &rod2D) {
     const auto &new_pos = rod2D.getPos();
     auto old_pos = start_rod.getPos();
-    auto msd = (new_pos - old_pos).squaredNorm();
+    auto diff = (new_pos - old_pos);
+    auto msd = diff.dot(diff);
     agg.feed(msd);
 }
 
 void Compute::evalX(const Rod2d &rod2D) {
     const auto &new_pos = rod2D.getPos();
-    auto msd = (new_pos).sum();
+    auto msd = (new_pos[0]);
     agg.feed(msd);
 }
 
 void Compute::evalX_squared(const Rod2d &rod2D) {
     const auto &new_pos = rod2D.getPos();
-    auto msd = (new_pos).squaredNorm();
+    auto msd = pow(new_pos[0],2.0);
     agg.feed(msd);
 }
 
@@ -103,7 +104,7 @@ Compute::Compute(Rod2d rod, Type t_type, Force force, size_t t_every,
                     target = 2 * (rod.getDiff().trace()) * time;
                     break;
                 case Force::const_F:
-                    target = static_cast<double>(NAN);
+                    target = 0.0;
                     break;
                 case Force::harmonic_F:
                     target = 2 * (1.0 - std::exp(-2 * time));
@@ -132,7 +133,7 @@ Compute::Compute(Rod2d rod, Type t_type, Force force, size_t t_every,
                     break;
                 case Force::const_F:
                 case Force::harmonic_F:
-                    target = static_cast<double>(NAN);
+                    target = 0.0;
                     break;
             }
         } break;
@@ -149,7 +150,7 @@ Compute::Compute(Rod2d rod, Type t_type, Force force, size_t t_every,
                     break;
                 case Force::const_F:
                 case Force::harmonic_F:
-                    target = static_cast<double>(NAN);
+                    target = 0.0;
                     break;
             }
         } break;
@@ -161,10 +162,10 @@ Compute::Compute(Rod2d rod, Type t_type, Force force, size_t t_every,
                     target = rod.getPos().sum();
                     break;
                 case Force::const_F:
-                    target = static_cast<double>(NAN);
+                    target = 0.0;
                     break;
                 case Force::harmonic_F:
-                    target = rod.getPos().sum() * exp(-1 * time);
+                    target = rod.getPos()[0] * exp(-1 * time);
                     break;
             }
 
@@ -177,14 +178,14 @@ Compute::Compute(Rod2d rod, Type t_type, Force force, size_t t_every,
 
             switch (force) {
                 case Force::zero_F:
-                    target = static_cast<double>(NAN);
+                    target = 0.0;
                     break;
                 case Force::const_F:
-                    target = static_cast<double>(NAN);
+                    target = 0.0;
                     break;
                 case Force::harmonic_F:
-                    target = rod.getPos().squaredNorm() * exp(-2 * time) +
+                    target = pow(rod.getPos()[0],2) * exp(-2 * time) +
-                             2 * (1 - exp(-2 * time));
+                             1 * (1 - exp(-2 * time));
                     break;
             }
         } break;
@@ -195,7 +196,7 @@ const LiveAgg &Compute::getAgg() const { return agg; }
 
 Compute::Type Compute::getType() const { return type; }
 
-double Compute::getTime() const { return static_cast<double>(every); }
+double Compute::getTime() const { return time; }
 
 double Compute::getTarget() const { return target; }
 

C++/src/Integratoren2d.cpp

@@ -8,8 +8,8 @@
 using Vector = Eigen::Vector2d;
 using Matrix = Eigen::Matrix2d;
 
-Vector
+Vector force_euler(const Rod2d &rod2D, const Simulation &sim,
-force_euler(const Rod2d &rod2D, const Simulation &sim, const std::function<Vector(Vector, Vector)> &force) {
+                   const std::function<Vector(Vector, Vector)> &force) {
     if (not force) {
         return {0.0, 0.0};
     }
@@ -21,16 +21,17 @@ force_euler(const Rod2d &rod2D, const Simulation &sim, const std::function<Vecto
     return F_trans;
 }
 
-Vector torque_euler(const Rod2d &rod2D, const Simulation &sim, const std::function<double(Vector, Vector)> &torque) {
+Vector torque_euler(const Rod2d &rod2D, const Simulation &sim,
+                    const std::function<double(Vector, Vector)> &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());
+    return torque(rod2D.getPos(), e) * rod2D.getDRot() * sim.getMDeltaT() *
+           ortho(rod2D.getE());
 }
 
-Vector
+Vector rand_trans_euler(const Rod2d &rod2D, Simulation &sim) {
-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
@@ -44,10 +45,9 @@ Vector rand_rot_euler(const Rod2d &rod2D, Simulation &sim) {
 }
 
 void Integratoren2d::Set1_Euler_f(
-        Rod2d &rod2D, Simulation &sim,
+    Rod2d &rod2D, Simulation &sim,
-        const std::function<Vector(Vector, Vector)> &force,
+    const std::function<Vector(Vector, Vector)> &force,
-        const std::function<double(Vector, Vector)> &torque) {
+    const std::function<double(Vector, Vector)> &torque) {
-
     // Gaussian - translation
     auto trans = rand_trans_euler(rod2D, sim);
 
@@ -69,48 +69,49 @@ void Integratoren2d::Set1_Euler_f(
     rod2D.setE(integrated_e);
 }
 
-void Integratoren2d::Set1_Euler(
+void Integratoren2d::Set1_Euler(Rod2d &rod2D, Simulation &sim) {
-        Rod2d &rod2D, Simulation &sim) {
     Set1_Euler_f(rod2D, sim, {}, {});
 }
 
-
 void Integratoren2d::Set2_Heun_f(
-        Rod2d &rod2D, Simulation &sim,
+    Rod2d &rod2D, Simulation &sim,
-        const std::function<Vector(Vector, Vector)> &force,
+    const std::function<Vector(Vector, Vector)> &force,
-        const std::function<double(Vector, Vector)> &torque) {
+    const std::function<double(Vector, Vector)> &torque) {
-    //position
+    // position
-    Vector pos_predictor = rod2D.getPos() + rand_trans_euler(rod2D, sim) + force_euler(rod2D, sim, force);
+    Vector pos_predictor = rod2D.getPos() + rand_trans_euler(rod2D, sim) +
-    //rotation
+                           force_euler(rod2D, sim, force);
-    Vector delta_e_predictor = rand_rot_euler(rod2D, sim) + torque_euler(rod2D, sim, torque);
+    // 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
+    // make predicted Particle
     Rod2d pred_rod = rod2D;
     pred_rod.setPos(pos_predictor);
     pred_rod.setE(e_predictor);
 
-    //integrate euler for future
+    // Heun integration
-    Vector delta_e_future = rand_rot_euler(pred_rod, sim) + torque_euler(pred_rod, sim, torque);
+    Vector e_integrated = rod2D.getE() +
-
+                          0.5 * (ortho(rod2D.getE()) + ortho(e_predictor)) *
-    //Heun integration
+                              sim.getNorm(sim.getSTD()) * rod2D.getDRot_Sqrt() +
-    Vector e_integrated =
+                          0.5 * (torque_euler(rod2D, sim, torque) +
-            rod2D.getE() + 0.5 * (delta_e_predictor + delta_e_future);
+                                 torque_euler(pred_rod, sim, torque));
+    ;
     // apply
     rod2D.setPos(pos_predictor);
     rod2D.setE(e_integrated.normalized());
 }
 
-void Integratoren2d::Set2_Heun(
+void Integratoren2d::Set2_Heun(Rod2d &rod2D, Simulation &sim) {
-        Rod2d &rod2D, Simulation &sim) {
     Set2_Heun_f(rod2D, sim, {}, {});
 }
 
 void Integratoren2d::Set3_Exact_f(
-        Rod2d &rod2D, Simulation &sim,
+    Rod2d &rod2D, Simulation &sim,
-        const std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)> &force,
+    const std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)>
-        const std::function<double(Eigen::Vector2d, Eigen::Vector2d)> &torque) {
+        &force,
+    const std::function<double(Eigen::Vector2d, Eigen::Vector2d)> &torque) {
     // Gaussian - translation
     auto trans = rand_trans_euler(rod2D, sim);
 
@@ -125,8 +126,7 @@ void Integratoren2d::Set3_Exact_f(
     // torque
     auto M_rot = torque_euler(rod2D, sim, torque);
 
-    auto correction =
+    auto correction = rod2D.getDRot() * sim.getMDeltaT() * rod2D.getE();
-            -0.5 * pow(rod2D.getDRot(), 2) * sim.getMDeltaT() * rod2D.getE();
 
     Vector integrated_e = rod2D.getE() + rot + M_rot + correction;
     integrated_e.normalize();
@@ -140,9 +140,10 @@ void Integratoren2d::Set3_Exact(Rod2d &rod2D, Simulation &sim) {
 }
 
 void Integratoren2d::Set4_BDAS_f(
-        Rod2d &rod2D, Simulation &sim,
+    Rod2d &rod2D, Simulation &sim,
-        const std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)> &force,
+    const std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)>
-        const std::function<double(Eigen::Vector2d, Eigen::Vector2d)> &torque) {
+        &force,
+    const std::function<double(Eigen::Vector2d, Eigen::Vector2d)> &torque) {
     // Gaussian - translation
     auto trans = rand_trans_euler(rod2D, sim);
     // Force
@@ -152,7 +153,9 @@ void Integratoren2d::Set4_BDAS_f(
     // 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;
+    double M_rot = torque ? torque(rod2D.getPos(), rod2D.getE()) *
+                                rod2D.getDRot() * sim.getMDeltaT()
+                          : 0.0;
 
     Vector integrated_e = Eigen::Rotation2D<double>(rot + M_rot) * rod2D.getE();
     integrated_e.normalize();
@@ -165,19 +168,21 @@ void Integratoren2d::Set4_BDAS(Rod2d &rod2D, Simulation &sim) {
     Set4_BDAS_f(rod2D, sim, {}, {});
 }
 
-
+void Integratoren2d::Set5_MBD_f(
-void Integratoren2d::Set5_MBD_f(Rod2d &rod2D, Simulation &sim,
+    Rod2d &rod2D, Simulation &sim,
-                                const std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)> &force,
+    const std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)>
-                                const std::function<double(Eigen::Vector2d, Eigen::Vector2d)> &torque) {
+        &force,
-
+    const std::function<double(Eigen::Vector2d, Eigen::Vector2d)> &torque) {
-    //position
+    // position
-    Vector pos_predictor = rod2D.getPos() + rand_trans_euler(rod2D, sim) + force_euler(rod2D, sim, force);
+    Vector pos_predictor = rod2D.getPos() + rand_trans_euler(rod2D, sim) +
-    //rotation
+                           force_euler(rod2D, sim, force);
-    Vector delta_e_predictor = rand_rot_euler(rod2D, sim) + torque_euler(rod2D, sim, torque);
+    // 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
+    // make predicted Particle
     Rod2d pred_rod = rod2D;
     pred_rod.setPos(pos_predictor);
     pred_rod.setE(e_predictor);
@@ -185,17 +190,23 @@ void Integratoren2d::Set5_MBD_f(Rod2d &rod2D, Simulation &sim,
     auto std = sim.getSTD();
     Vector e = rod2D.getE();
 
-    auto rot_Mat = [](const Matrix& mat,const Matrix& rotMat) { return rotMat * mat * rotMat.transpose(); };
+    auto rot_Mat = [](const Matrix &mat, const Matrix &rotMat) {
-    Eigen::Matrix2d Diff_combined = 0.5 * (rot_Mat(rod2D.getDiff(), rotation_Matrix(e)) +
+        return rotMat * mat * rotMat.transpose();
-                                           rot_Mat(rod2D.getDiff(), rotation_Matrix(e_predictor)));
+    };
+    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 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));
+    Eigen::Vector2d pos_integrated = rod2D.getPos() + D_combined_sqrt * rand +
+                                     0.5 * (force_euler(rod2D, sim, force) +
+                                            force_euler(pred_rod, sim, force));
 
-    //apply
+    // apply
     rod2D.setPos(pos_integrated);
     rod2D.setE(e_integrated.normalized());
 }

C++/src/main.cpp

@@ -10,12 +10,12 @@
 #include "Integratoren2d.h"
 #include "force_lambdas.h"
 constexpr size_t SEED = 1234;
-constexpr double stepSize = 0.001;
+constexpr double stepSize = 0.01;
 constexpr size_t n_computes = 100;
 constexpr size_t num_integratoren = 5;
 constexpr size_t num_forces = 3;
 constexpr size_t num_length = 3;
-constexpr size_t delta_compute = 4;
+constexpr size_t delta_compute = 1;
 constexpr size_t numStep = 10000;
 inline const std::string header = {"time,val,target,std\n"};
 void run(size_t integrator_index, size_t force_index, size_t length_index) {
@@ -134,6 +134,8 @@ void run(size_t integrator_index, size_t force_index, size_t length_index) {
 }
 
 int main() {
+    run(0, 2, 0);
+    return EXIT_SUCCESS;
 #pragma omp parallel for default(none)
     for (size_t j = 0; j < num_forces * num_length * num_integratoren; ++j) {
         run(j % (num_integratoren),