diff --git a/Auswertung/harmonic.py b/Auswertung/harmonic.py new file mode 100644 index 0000000..cc59825 --- /dev/null +++ b/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 ") + 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(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.show() + + +if __name__ == "__main__": + main() diff --git a/Auswertung/main.py b/Auswertung/main.py index 4a1d53e..acfe38f 100644 --- a/Auswertung/main.py +++ b/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: - sb.lineplot(x="time", y="target", data=file) + 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) - # 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__": diff --git a/Auswertung/requirements.txt b/Auswertung/requirements.txt index c8628cc..9ee03f4 100644 --- a/Auswertung/requirements.txt +++ b/Auswertung/requirements.txt @@ -1,4 +1,4 @@ -numpy~=1.21.3 +numpy~=1.19.5 seaborn~=0.11.2 -matplotlib~=3.4.3 -pandas~=1.3.4 \ No newline at end of file +matplotlib~=3.3.4 +pandas~=1.1.5 \ No newline at end of file diff --git a/C++/src/Compute.cpp b/C++/src/Compute.cpp index 68eab69..c138809 100644 --- a/C++/src/Compute.cpp +++ b/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(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(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(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(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(NAN); + target = 0.0; break; case Force::const_F: - target = static_cast(NAN); + target = 0.0; break; case Force::harmonic_F: - target = rod.getPos().squaredNorm() * exp(-2 * time) + - 2 * (1 - exp(-2 * time)); + target = pow(rod.getPos()[0],2) * 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(every); } +double Compute::getTime() const { return time; } double Compute::getTarget() const { return target; } diff --git a/C++/src/Integratoren2d.cpp b/C++/src/Integratoren2d.cpp index c664a34..8714001 100644 --- a/C++/src/Integratoren2d.cpp +++ b/C++/src/Integratoren2d.cpp @@ -8,8 +8,8 @@ using Vector = Eigen::Vector2d; using Matrix = Eigen::Matrix2d; -Vector -force_euler(const Rod2d &rod2D, const Simulation &sim, const std::function &force) { +Vector force_euler(const Rod2d &rod2D, const Simulation &sim, + const std::function &force) { if (not force) { return {0.0, 0.0}; } @@ -21,16 +21,17 @@ force_euler(const Rod2d &rod2D, const Simulation &sim, const std::function &torque) { +Vector torque_euler(const Rod2d &rod2D, const Simulation &sim, + const std::function &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 -rand_trans_euler(const Rod2d &rod2D, Simulation &sim) { +Vector 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, - const std::function &force, - const std::function &torque) { - + Rod2d &rod2D, Simulation &sim, + const std::function &force, + const std::function &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( - Rod2d &rod2D, Simulation &sim) { +void Integratoren2d::Set1_Euler(Rod2d &rod2D, Simulation &sim) { Set1_Euler_f(rod2D, sim, {}, {}); } - void Integratoren2d::Set2_Heun_f( - Rod2d &rod2D, Simulation &sim, - const std::function &force, - const std::function &torque) { - //position - Vector pos_predictor = rod2D.getPos() + rand_trans_euler(rod2D, sim) + force_euler(rod2D, sim, force); - //rotation - Vector delta_e_predictor = rand_rot_euler(rod2D, sim) + torque_euler(rod2D, sim, torque); + Rod2d &rod2D, Simulation &sim, + const std::function &force, + const std::function &torque) { + // position + Vector pos_predictor = rod2D.getPos() + rand_trans_euler(rod2D, sim) + + force_euler(rod2D, sim, force); + // 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 - Vector delta_e_future = rand_rot_euler(pred_rod, sim) + torque_euler(pred_rod, sim, torque); - - //Heun integration - Vector e_integrated = - rod2D.getE() + 0.5 * (delta_e_predictor + delta_e_future); + // Heun integration + Vector e_integrated = rod2D.getE() + + 0.5 * (ortho(rod2D.getE()) + ortho(e_predictor)) * + sim.getNorm(sim.getSTD()) * rod2D.getDRot_Sqrt() + + 0.5 * (torque_euler(rod2D, sim, torque) + + torque_euler(pred_rod, sim, torque)); + ; // apply rod2D.setPos(pos_predictor); rod2D.setE(e_integrated.normalized()); } -void Integratoren2d::Set2_Heun( - Rod2d &rod2D, Simulation &sim) { +void Integratoren2d::Set2_Heun(Rod2d &rod2D, Simulation &sim) { Set2_Heun_f(rod2D, sim, {}, {}); } void Integratoren2d::Set3_Exact_f( - Rod2d &rod2D, Simulation &sim, - const std::function &force, - const std::function &torque) { + Rod2d &rod2D, Simulation &sim, + const std::function + &force, + const std::function &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 = - -0.5 * pow(rod2D.getDRot(), 2) * sim.getMDeltaT() * rod2D.getE(); + auto correction = rod2D.getDRot() * 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, - const std::function &force, - const std::function &torque) { + Rod2d &rod2D, Simulation &sim, + const std::function + &force, + const std::function &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(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(Rod2d &rod2D, Simulation &sim, - const std::function &force, - const std::function &torque) { - - //position - Vector pos_predictor = rod2D.getPos() + rand_trans_euler(rod2D, sim) + force_euler(rod2D, sim, force); - //rotation - Vector delta_e_predictor = rand_rot_euler(rod2D, sim) + torque_euler(rod2D, sim, torque); +void Integratoren2d::Set5_MBD_f( + Rod2d &rod2D, Simulation &sim, + const std::function + &force, + const std::function &torque) { + // position + Vector pos_predictor = rod2D.getPos() + rand_trans_euler(rod2D, sim) + + force_euler(rod2D, sim, force); + // 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(); }; - Eigen::Matrix2d Diff_combined = 0.5 * (rot_Mat(rod2D.getDiff(), rotation_Matrix(e)) + - rot_Mat(rod2D.getDiff(), rotation_Matrix(e_predictor))); + auto rot_Mat = [](const Matrix &mat, const Matrix &rotMat) { + return rotMat * mat * rotMat.transpose(); + }; + 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()); } diff --git a/C++/src/main.cpp b/C++/src/main.cpp index f4c7032..953e95a 100644 --- a/C++/src/main.cpp +++ b/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),