Big restructure

.gitignore

@@ -2,6 +2,7 @@
 build/
 out/
 cmake-build-*/
+data/
 
 # IDE files
 .vs/
@@ -25,4 +26,4 @@ cmake-build-*/
 $RECYCLE.BIN/
 .TemporaryItems
 ehthumbs.db
-Thumbs.db
+Thumbs.db

Auswertung/main.py

@@ -0,0 +1,54 @@
+import glob
+
+import matplotlib.pyplot as plt
+import pandas as pd
+import seaborn as sb
+
+
+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):
+
+            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)
+                    plt.ylabel(compute)
+                    legend.append(integrator)
+                    if index_integrator == 0:
+                        sb.lineplot(x="time", y="target", data=file)
+                        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):
+        ax.annotate(col, xy=(0.5, 1), xytext=(0, pad),
+                    xycoords='axes fraction', textcoords='offset points',
+                    size='large', ha='center', va='baseline', fontsize=25)
+
+    for ax, row in zip(axs[:, 0], forces):
+        ax.annotate(row, xy=(0, 0.5), xytext=(-ax.yaxis.labelpad - pad, 0),
+                    xycoords=ax.yaxis.label, textcoords='offset points',
+                    size='large', ha='right', va='center', fontsize=25)
+
+    fig.tight_layout()
+
+
+def main():
+    forces = ["zero", "const", "harmonic"]
+    integrators = ["euler", "heun", "exact", "bdas", "mbd"]
+    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__":
+    main()

Auswertung/requirements.txt

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

C++/src/Compute.cpp

@@ -70,10 +70,14 @@ void Compute::eval(const Rod2d &rod2D) {
             case Type::msd_const:
                 break;
             case Type::x:
-                evalX(rod2D);
+                if (time_step == every) {
+                    evalX(rod2D);
+                }
                 break;
             case Type::x_squared:
-                evalMSD(rod2D);
+                if (time_step == every) {
+                    evalMSD(rod2D);
+                }
                 break;
         }
         if (resetting) {
@@ -84,9 +88,11 @@ void Compute::eval(const Rod2d &rod2D) {
 }
 
 Compute::Compute(Rod2d rod, Type t_type, size_t t_every, Simulation &sim)
-        : start_rod(std::move(rod)), every(t_every), time_step(0), type(t_type),
-          time(sim.getMDeltaT() * static_cast<double>(every)) {
-
+    : start_rod(std::move(rod)),
+      every(t_every),
+      time_step(0),
+      type(t_type),
+      time(sim.getMDeltaT() * static_cast<double>(every)) {
     switch (type) {
         case Type::msd: {
             resetting = true;
@@ -106,46 +112,39 @@ Compute::Compute(Rod2d rod, Type t_type, size_t t_every, Simulation &sim)
             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 *
-                             (1 - exp(-u * rod.getDRot() * time)) / u /
-                             rod.getDRot() / time;
-        }
-            break;
+            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 *
-                             (1 - exp(-u * rod.getDRot() * time)) / u /
-                             rod.getDRot() / time;
-        }
-            break;
+            target = Dmean + deltaD / 2 * (1 - exp(-u * rod.getDRot() * time)) /
+                                 u / rod.getDRot() / time;
+        } break;
         case Type::msd_harmonic_k_1: {
             resetting = true;
             type_str = "msd_harmonic";
             target = 2 * (1.0 - std::exp(-2 * time));
-        }
-            break;
+        } break;
         case Type::msd_const: {
             resetting = true;
             type_str = "msd_const";
             target = -1.0;
-        }
-            break;
+        } break;
         case Type::x: {
             resetting = false;
             type_str = "<x>";
-            target = -1.0;
-        }
-            break;
+            target = rod.getPos().norm()*exp(-1.0*time);
+        } break;
         case Type::x_squared: {
             resetting = false;
             type_str = "<x_squared>";
-            target = -1.0;
-        }
-            break;
+            target = rod.getPos().squaredNorm()*exp(-2.0*time)+
+                2.0*(1-exp(-2.0*time));
+        } break;
     }
 }
 

C++/src/Integratoren2d.cpp

@@ -92,7 +92,7 @@ void Integratoren2d::Set2_Heun_f(
     pred_rod.setE(e_predictor);
 
     //integrate euler for future
-    Vector delta_e_future = rand_rot_euler(rod2D, sim) + torque_euler(rod2D, sim, torque);
+    Vector delta_e_future = rand_rot_euler(pred_rod, sim) + torque_euler(pred_rod, sim, torque);
 
     //Heun integration
     Vector e_integrated =

C++/src/main.cpp

@@ -0,0 +1,132 @@
+//
+// Created by jholder on 22.10.21.
+//
+
+#include <fstream>
+#include <iostream>
+
+#include <filesystem>
+#include "Calculation.h"
+#include "Integratoren2d.h"
+#include "force_lambdas.h"
+constexpr size_t SEED = 1234;
+constexpr double stepSize = 0.001;
+constexpr int n_computes = 100;
+constexpr size_t num_runs = 5;
+constexpr int numStep = 1000000;
+constexpr int delta_compute = 4;
+void run_no_force(size_t integrator_index, size_t force_index,
+                  size_t length_index) {
+    using type = std::pair<Calculation::inte_force_type, std::string>;
+    std::vector<type> 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<
+        std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)>>(
+        {zero_Force, const_Force, harmonic_Force})[force_index];
+    auto force_name = std::vector<std::string>(
+        {"zero_force_", "const_force_", "harmonic_force_"})[force_index];
+
+    namespace fs = std::filesystem;
+
+    std::string folder = "out/" + force_name + name + "_L"+std::to_string(length_index)+"_";
+    fs::create_directories("out");
+    double length = std::vector<double>({1.0, 1.5, 2.0})[length_index];
+    {
+        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<std::pair<Compute::Type, size_t>> computes;
+        for (int i = 1; i < n_computes; ++i) {
+            computes.emplace_back(MSD_VARIANT, i * delta_compute);
+            computes.emplace_back(Compute::Type::oaf, i * delta_compute);
+            computes.emplace_back(Compute::Type::empxx, i * delta_compute);
+            computes.emplace_back(Compute::Type::empyy, i * delta_compute);
+        }
+
+        Calculation euler(integrator, computes, stepSize, SEED, force,
+                          zero_Torque, length);
+
+        euler.run(numStep);
+        std::ofstream msdFile(folder + "msd.dat");
+        std::ofstream oafFile(folder + "oaf.dat");
+        std::ofstream empxxFile(folder + "empxx.dat");
+        std::ofstream empyyFile(folder + "empyy.dat");
+
+        msdFile << "time,val,target\n";
+        oafFile << "time,val,target\n";
+        empxxFile << "time,val,target\n";
+        empyyFile << "time,val,target\n";
+
+        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<std::pair<Compute::Type, size_t>> 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 < numStep / n_computes; ++i) {
+            calc_repeat.run(n_computes);
+            calc_repeat.reset();
+        }
+
+        std::ofstream xFile(folder + "x.dat");
+        std::ofstream xsqFile(folder + "x_squared.dat");
+
+        xFile << "time,val,target\n";
+        xsqFile << "time,val,target\n";
+
+        for (const auto &com : calc_repeat.getComputes()) {
+            if (com.getType() == Compute::Type::x) {
+                xFile << com.getTime() << ", " << com.getAgg().getMean() << ", "
+                      << com.getTarget() << std::endl;
+            }
+            if (com.getType() == Compute::Type::x_squared) {
+                xsqFile << com.getTime() << ", " << com.getAgg().getMean()
+                        << ", " << com.getTarget() << std::endl;
+            }
+        }
+    }
+
+    std::cout << "Finished run " << folder << std::endl;
+}
+
+int main() {
+    for (size_t j = 0; j < 3; ++j) {
+        for (size_t i = 0; i < num_runs; ++i) {
+            for (size_t k = 0; k < 1; ++k) {
+                run_no_force(i, j, k);
+            }
+        }
+    }
+}

src/main.cpp

@@ -1,119 +0,0 @@
-//
-// Created by jholder on 22.10.21.
-//
-
-#include <iostream>
-#include <fstream>
-
-#include "Calculation.h"
-#include "Integratoren2d.h"
-#include "force_lambdas.h"
-
-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<Calculation::inte_force_type, std::string>;
-    std::vector<type> 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<std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)>>({zero_Force, const_Force,
-                                                                                                harmonic_Force})[force_index];
-    auto force_name = std::vector<std::string>({"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<std::pair<Compute::Type, size_t>> 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<std::pair<Compute::Type, size_t>> 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 <numStep/n_computes; ++i) {
-            calc_repeat.run(n_computes);
-            calc_repeat.reset();
-        }
-
-
-        std::ofstream xFile(folder + "x.dat");
-        std::ofstream xsqFile(folder + "x_squared.dat");
-
-        for (const auto &com: calc_repeat.getComputes()) {
-            if (com.getType() == Compute::Type::x) {
-                xFile << com.getTime() << ", " << com.getAgg().getMean() << ", " << com.getTarget() << std::endl;
-            }
-            if (com.getType() == Compute::Type::x_squared) {
-                xsqFile << com.getTime() << ", " << com.getAgg().getMean() << ", " << com.getTarget() << std::endl;
-            }
-        }
-    }
-
-    std::cout << "Finished run " << folder << std::endl;
-}
-
-
-int main() {
-    for (size_t j = 0; j < 3; ++j) {
-        for (size_t i = 0; i < num_runs; ++i) {
-            run_no_force(i, j);
-        }
-    }
-}