From 82e81eb91e8de5c37f51abc8defdee60ca2dcd81 Mon Sep 17 00:00:00 2001 From: Jacob Date: Thu, 16 Dec 2021 16:20:41 +0100 Subject: [PATCH] Konv added --- Auswertung/harmonic.py | 24 ++++---- Auswertung/main.py | 8 ++- C++/src/Calculation.cpp | 4 ++ C++/src/Compute.cpp | 73 ++++++++++++----------- C++/src/Rod2d.cpp | 5 +- C++/src/main.cpp | 126 +++++++++++++++++++++++++++------------- 6 files changed, 145 insertions(+), 95 deletions(-) diff --git a/Auswertung/harmonic.py b/Auswertung/harmonic.py index cc59825..fad46a7 100644 --- a/Auswertung/harmonic.py +++ b/Auswertung/harmonic.py @@ -1,38 +1,38 @@ -import glob - import matplotlib.pyplot as plt -import pandas as pd -import seaborn as sb import numpy as np +import pandas as pd def main(): time = np.arange(0.0, 4.0, 0.01) - xstart = 1.0 + 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_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)) + msd_mean = (1 - np.exp(-2 * time)) + msd_2_mean = 2 * (1 - 2 * np.exp(- time) + np.exp(-2 * time)) + 2 * (1 - np.exp(-2 * time)) + plt.plot(time, x_mean, label="a ") plt.plot(x["time"], x["val"], label="") plt.plot(time, x_sq_mean, label="a ") plt.plot(x_sq["time"], x_sq["val"], label="") - plt.plot(time, 2*msd_mean,label="a msd") + plt.plot(time, 2 * msd_mean, label="a msd") + plt.plot(time, msd_2_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="-²") - plt.legend(loc = 1) + # plt.plot(x["time"], x_sq["val"] - np.power(x["val"], 2.0)) + plt.plot(x["time"], 2 * (x_sq["val"] - np.power(x["val"], 2.0)), label="-²") + plt.legend(loc=1) plt.show() diff --git a/Auswertung/main.py b/Auswertung/main.py index acfe38f..348ebed 100644 --- a/Auswertung/main.py +++ b/Auswertung/main.py @@ -21,11 +21,15 @@ def make_figure(length_index, length, forces, integrators, computes): #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 and not np.any(file["target"].isnull()): plt.plot(file["time"], file["target"]) legend.append("analytisch") plt.legend(title='iterators', labels=legend) + for ax1 in axs: + for ax in ax1: + ax.axvline(1.0) pad = 5 # in points for ax, col in zip(axs[0], computes): @@ -40,13 +44,13 @@ def make_figure(length_index, length, forces, integrators, computes): fig.tight_layout() fig.savefig(f"fig{length_index}.png") - + plt.show() 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"]): + for l_i, l in enumerate(["sphere"]): #["sphere", "L = 1.5", "L = 2.0"] make_figure(length=l, length_index=l_i, forces=forces, computes=computes, integrators=integrators) diff --git a/C++/src/Calculation.cpp b/C++/src/Calculation.cpp index b858535..4f47d41 100644 --- a/C++/src/Calculation.cpp +++ b/C++/src/Calculation.cpp @@ -4,15 +4,19 @@ #include "Calculation.h" +#include #include void Calculation::run(size_t steps) { + std::ofstream file ("dump.xyz"); for (size_t step = 0; step < steps; ++step) { m_integrator(rod, sim); for (auto &comp : computes) { comp.eval(rod); } + file<(every)) { switch (type) { case Type::msd: { - resetting = true; + resetting = false; type_str = "msd"; switch (force) { case Force::zero_F: @@ -107,20 +104,20 @@ Compute::Compute(Rod2d rod, Type t_type, Force force, size_t t_every, target = 0.0; break; case Force::harmonic_F: - target = 2 * (1.0 - std::exp(-2 * time)); + target = 2 * (1 - 2 * std::exp(- time) + std::exp(-2 * time)) + 2 * (1 - std::exp(-2 * time)); break; } break; } case Type::oaf: { - resetting = true; + resetting = false; type_str = "oaf"; target = std::exp(-rod.getDRot() * time); break; } case Type::empxx: { - resetting = true; + resetting = false; type_str = "empxx"; const double Dmean = 0.5 * (rod.getDiff().trace()); const double u = 4.0; @@ -138,14 +135,15 @@ Compute::Compute(Rod2d rod, Type t_type, Force force, size_t t_every, } } break; case Type::empyy: { - resetting = true; + resetting = false; 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); switch (force) { case Force::zero_F: - target = Dmean + deltaD / 2 * (1 - exp(-u * rod.getDRot() * time)) / + target = Dmean + deltaD / 2 * + (1 - exp(-u * rod.getDRot() * time)) / u / rod.getDRot() / time; break; case Force::const_F: @@ -169,8 +167,6 @@ Compute::Compute(Rod2d rod, Type t_type, Force force, size_t t_every, break; } - - } break; case Type::x_squared: { resetting = false; @@ -184,12 +180,15 @@ Compute::Compute(Rod2d rod, Type t_type, Force force, size_t t_every, target = 0.0; break; case Force::harmonic_F: - target = pow(rod.getPos()[0],2) * exp(-2 * time) + + target = pow(rod.getPos()[0], 2) * exp(-2 * time) + 1 * (1 - exp(-2 * time)); break; } } break; + + } + resetting = false; } const LiveAgg &Compute::getAgg() const { return agg; } diff --git a/C++/src/Rod2d.cpp b/C++/src/Rod2d.cpp index f306612..92c2f32 100644 --- a/C++/src/Rod2d.cpp +++ b/C++/src/Rod2d.cpp @@ -20,8 +20,8 @@ Rod2d::Rod2d(double L) : m_pos({1, 1}), m_e({1, 0}) { const double p = L; auto D_para = D0 / (M_Pl * 2.0) * (log(p) - 0.1401 + 1.034 / p - 0.228 / (p * p)); - auto D_ortho = - D0 / (M_Pl * 4.0) * (log(p) + 0.8369 + 0.5551 / p - 0.06066 / (p * p)); + auto D_ortho = D0 / (M_Pl * 4.0) * + (log(p) + 0.8369 + 0.5551 / p - 0.06066 / (p * p)); m_D_rot = 3 * D0 / (M_Pl * L * L) * (log(p) - 0.3512 + 0.7804 / p - 0.09801 / (p * p)); m_Diff << D_para, 0, 0, D_ortho; @@ -46,4 +46,3 @@ double Rod2d::getDRot_Sqrt() const { return m_D_rot_sqrt; } const Eigen::Vector2d &Rod2d::getE() const { return m_e; } void Rod2d::setE(const Eigen::Vector2d &mE) { m_e = mE; } - diff --git a/C++/src/main.cpp b/C++/src/main.cpp index 953e95a..85f0d51 100644 --- a/C++/src/main.cpp +++ b/C++/src/main.cpp @@ -6,18 +6,21 @@ #include #include +#include #include "Calculation.h" #include "Integratoren2d.h" #include "force_lambdas.h" -constexpr size_t SEED = 1234; +#include +constexpr size_t SEED = 1364; constexpr double stepSize = 0.01; -constexpr size_t n_computes = 100; +constexpr size_t n_computes = 10; constexpr size_t num_integratoren = 5; constexpr size_t num_forces = 3; -constexpr size_t num_length = 3; -constexpr size_t delta_compute = 1; -constexpr size_t numStep = 10000; -inline const std::string header = {"time,val,target,std\n"}; +constexpr size_t num_length = 1; +constexpr size_t delta_compute = 10; +constexpr size_t numStep = 50000; +inline const std::string header = {"time,val,target,std,numPoints\n"}; + void run(size_t integrator_index, size_t force_index, size_t length_index) { using type = std::pair; std::vector vec({type(Integratoren2d::Set1_Euler_f, "euler"), @@ -41,7 +44,8 @@ void run(size_t integrator_index, size_t force_index, size_t length_index) { std::string folder = "out/" + force_name + name + "_L" + std::to_string(length_index) + "_"; fs::create_directories("out"); - double length = std::vector({1.0, 1.5, 2.0})[length_index]; + // double length = std::vector({1.0, 1.5, 2.0})[length_index]; + double length = 1.0; { std::vector> computes; for (size_t i = 1; i < n_computes; ++i) { @@ -49,6 +53,8 @@ void run(size_t integrator_index, size_t force_index, size_t length_index) { 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); + computes.emplace_back(Compute::Type::x, i * delta_compute); + computes.emplace_back(Compute::Type::x_squared, i * delta_compute); } Calculation euler(integrator, computes, stepSize, SEED, force, @@ -56,12 +62,18 @@ void run(size_t integrator_index, size_t force_index, size_t length_index) { for (size_t i = 0; i < numStep; ++i) { euler.run(n_computes * delta_compute); + euler.reset(); } + std::ofstream msdFile(folder + "msd.dat"); std::ofstream oafFile(folder + "oaf.dat"); std::ofstream empxxFile(folder + "empxx.dat"); std::ofstream empyyFile(folder + "empyy.dat"); + std::ofstream xFile(folder + "x.dat"); + std::ofstream xsqFile(folder + "x_squared.dat"); + xFile << header; + xsqFile << header; msdFile << header; oafFile << header; empxxFile << header; @@ -71,70 +83,102 @@ void run(size_t integrator_index, size_t force_index, size_t length_index) { if (com.getType() == Compute::Type::msd) { msdFile << com.getTime() << ", " << com.getAgg().getMean() << ", " << com.getTarget() << ", " - << com.getAgg().getSEM() << std::endl; + << com.getAgg().getSEM() << ", " + << com.getAgg().getNumPoints() << std::endl; } if (com.getType() == Compute::Type::oaf) { oafFile << com.getTime() << ", " << com.getAgg().getMean() << ", " << com.getTarget() << ", " - << com.getAgg().getSEM() << std::endl; + << com.getAgg().getSEM() << ", " + << com.getAgg().getNumPoints() << std::endl; } if (com.getType() == Compute::Type::empxx) { empxxFile << com.getTime() << ", " << com.getAgg().getMean() << ", " << com.getTarget() << ", " - << com.getAgg().getSEM() << std::endl; + << com.getAgg().getSEM() << ", " + << com.getAgg().getNumPoints() << std::endl; } if (com.getType() == Compute::Type::empyy) { empyyFile << com.getTime() << ", " << com.getAgg().getMean() << ", " << com.getTarget() << ", " - << com.getAgg().getSEM() << std::endl; + << com.getAgg().getSEM() << ", " + << com.getAgg().getNumPoints() << std::endl; } if (com.getType() == Compute::Type::msd) { msdFile << com.getTime() << ", " << com.getAgg().getMean() << ", " << com.getTarget() << ", " - << com.getAgg().getSEM() << std::endl; + << com.getAgg().getSEM() << ", " + << com.getAgg().getNumPoints() << std::endl; } - } - } - { - std::vector> repeating_computes; - for (size_t i = 1; i < n_computes; ++i) { - repeating_computes.emplace_back(Compute::Type::x, - i * delta_compute); - repeating_computes.emplace_back(Compute::Type::x_squared, - i * delta_compute); - } - Calculation calc_repeat(integrator, repeating_computes, stepSize, SEED, - force, force_type, zero_Torque, 1.0); - for (size_t i = 0; i < numStep; ++i) { - calc_repeat.run(n_computes * delta_compute); - calc_repeat.reset(); - } - - std::ofstream xFile(folder + "x.dat"); - std::ofstream xsqFile(folder + "x_squared.dat"); - - xFile << header; - xsqFile << header; - - for (const auto &com : calc_repeat.getComputes()) { if (com.getType() == Compute::Type::x) { xFile << com.getTime() << ", " << com.getAgg().getMean() << ", " << com.getTarget() << ", " << com.getAgg().getSEM() - << std::endl; + << ", " << com.getAgg().getNumPoints() << std::endl; } if (com.getType() == Compute::Type::x_squared) { xsqFile << com.getTime() << ", " << com.getAgg().getMean() << ", " << com.getTarget() << ", " - << com.getAgg().getSEM() << std::endl; + << com.getAgg().getSEM() << ", " + << com.getAgg().getNumPoints() << std::endl; } } } - std::cout << "Finished run " << folder << std::endl; } -int main() { - run(0, 2, 0); + +void konv(int delta, std::ofstream &file, + Calculation::inte_force_type integrator, size_t seed) { + std::vector> computes; + double step_Size = 1.0 / delta; + computes.emplace_back(Compute::Type::msd, delta); + Calculation euler(integrator, computes, step_Size, seed, harmonic_Force, + Compute::Force::harmonic_F, zero_Torque, 1.0); + + for (int j = 0; j < 1000000; ++j) { + for (int i = 0; i < 1000; ++i) { + euler.run(static_cast(delta)); + euler.reset(); + } + std::cout << euler.getComputes()[0].getAgg().getSEM() << " " + << euler.getComputes()[0].getDifference() << std::endl; + if (euler.getComputes()[0].getAgg().getSEM() * 10 < + euler.getComputes()[0].getDifference()) { + file << delta << " " << euler.getComputes()[0].getAgg().getMean() + << " " << euler.getComputes()[0].getAgg().getSEM() << " " + << euler.getComputes()[0].getTarget() << " " + << euler.getComputes()[0].getDifference() << std::endl; + return; + } + } + file << "timeout: " << delta << " " + << euler.getComputes()[0].getAgg().getMean() << " " + << euler.getComputes()[0].getAgg().getSEM() << " " + << euler.getComputes()[0].getTarget() << " " + << euler.getComputes()[0].getDifference() << std::endl; +} +Calculation::inte_force_type integrator(int index) { + if (index % 5 == 0) return Integratoren2d::Set1_Euler_f; + if (index % 5 == 1) return Integratoren2d::Set2_Heun_f; + if (index % 5 == 2) return Integratoren2d::Set3_Exact_f; + if (index % 5 == 3) return Integratoren2d::Set4_BDAS_f; + if (index % 5 == 4) return Integratoren2d::Set5_MBD_f; + return Integratoren2d::Set1_Euler_f; +} +void konv_runner(int index) { + std::ofstream file; + file.open(fmt::format("out/{}_konv.out", index)); + for (int i = 1; i < 256; i *= 2) { + konv(i, file, integrator(index), static_cast(index / 5)); + } + file.close(); +} + +int main(int argc, char *argv[]) { + int index = 0; + if (argc >1) index = std::stoi(argv[1]); + konv_runner(index); + return EXIT_SUCCESS; #pragma omp parallel for default(none) for (size_t j = 0; j < num_forces * num_length * num_integratoren; ++j) {