jacob/BD_Integratoren
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Multiple cahnges :D

Browse Source
This commit is contained in:
Jacob Holder 2021-10-25 18:45:54 +02:00
parent 362ec38dd2
commit 8af7eee267
Signed by: jacob
GPG Key ID: 2194FC747048A7FD
22 changed files with 440 additions and 685 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
benchmark/bench.cpp
View File

@ -24,4 +24,7 @@ TEST_CASE("Euler - Baseline", "[benchmark]") {
BENCHMARK("BDAS without force") {
Integratoren2d_forceless::Set4_BDAS(rod, sim);
};
BENCHMARK("MBD without force") {
Integratoren2d_forceless::Set5_MBD(rod, sim);
};
}

8
src/Calculation.cpp
View File

@ -20,8 +20,8 @@ const Rod2d &Calculation::getRod() const { return rod; }
Calculation::Calculation(
std::function<void(Rod2d &, Simulation &)> t_integrator,
std::initializer_list<std::pair<Compute::Type, size_t>> t_computes,
double deltaT, size_t seed)
: sim(deltaT, seed), m_integrator(std::move(t_integrator)) {
double deltaT, size_t seed,double length)
: rod(length), sim(deltaT, seed),m_integrator(std::move(t_integrator)) {
for (const auto &pair : t_computes) {
computes.emplace_back(Compute(rod, pair.first, pair.second, sim));
}
@ -36,8 +36,8 @@ Calculation::Calculation(
std::initializer_list<std::pair<Compute::Type, size_t>> t_computes,
double deltaT, size_t seed,
std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)> force,
std::function<double(Eigen::Vector2d, Eigen::Vector2d)> torque)
: sim(deltaT, seed) {
std::function<double(Eigen::Vector2d, Eigen::Vector2d)> torque, double length)
:rod(length), sim(deltaT, seed) {
m_integrator = [&](Rod2d &t_rod, Simulation &t_sim) {
t_integrator(t_rod, t_sim, force, torque);
};

19
src/Calculation.h
View File

@ -11,29 +11,28 @@
class Calculation {
private:
Rod2d rod = Rod2d(1.0);
Rod2d rod;
Simulation sim;
std::function<void(Rod2d &, Simulation &)> m_integrator;
std::vector<Compute> computes;
public:
const Simulation &getSim() const;
[[nodiscard]] const Simulation &getSim() const;
using force_type = std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)>;
using torque_type = std::function<double(Eigen::Vector2d, Eigen::Vector2d)>;
using inte_force_type = std::function<void(Rod2d &, Simulation &, force_type, torque_type)>;
Calculation(
std::function<void(
Rod2d &, Simulation &,
std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)>,
std::function<double(Eigen::Vector2d, Eigen::Vector2d)>)>
inte_force_type
t_integrator,
std::initializer_list<std::pair<Compute::Type, size_t>> t_computes,
double deltaT, size_t seed,
std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)> force,
std::function<double(Eigen::Vector2d, Eigen::Vector2d)> torque);
force_type force,
torque_type torque, double length = 1.0);
explicit Calculation(
std::function<void(Rod2d &, Simulation &)> t_integrator,
std::initializer_list<std::pair<Compute::Type, size_t>> t_computes,
double deltaT, size_t seed);
double deltaT, size_t seed, double length = 1.0);
[[nodiscard]] const std::vector<Compute> &getComputes() const;

26
src/Compute.cpp
View File

@ -28,7 +28,7 @@ void Compute::eval_empXX(const Rod2d &rod2D) {
const auto &new_pos = rod2D.getPos();
auto old_pos = start_rod.getPos();
auto old_e = start_rod.getE();
double empxx = pow((new_pos - old_pos).dot(old_e), 2.0);
double empxx = pow((new_pos - old_pos).dot(old_e), 2.0)/2.0 / time;
agg.feed(empxx);
}
@ -37,7 +37,7 @@ void Compute::eval_empYY(const Rod2d &rod2D) {
auto old_pos = start_rod.getPos();
auto old_e = start_rod.getE();
Eigen::Vector2d old_e_ortho = {-old_e[1], old_e[0]};
double empxx = pow((new_pos - old_pos).dot(old_e_ortho), 2.0);
double empxx = pow((new_pos - old_pos).dot(old_e_ortho), 2.0)/2.0 / time;
agg.feed(empxx);
}
@ -64,15 +64,19 @@ 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) {
auto time = sim.getMDeltaT() * static_cast<double>(every);
time = sim.getMDeltaT() * static_cast<double>(every);
switch (type) {
case Type::msd:
case Type::msd:{
type_str = "msd";
target = 4.0 * 0.5 * (rod.getDiff().trace()) * time;
break;
case Type::oaf:
break;}
case Type::oaf: {
type_str = "oaf";
target = std::exp(-rod.getDRot() * time);
break;
}
case Type::empxx: {
type_str = "empxx";
const double Dmean = 0.5 * (rod.getDiff().trace());
const double u = 4.0;
const double deltaD = rod.getDiff()(1, 1) - rod.getDiff()(0, 0);
@ -81,12 +85,13 @@ Compute::Compute(Rod2d rod, Type t_type, size_t t_every, Simulation &sim)
rod.getDRot() / time;
} break;
case Type::empyy: {
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.0 *
(1 - exp(-u * rod.getDRot() * time)) / u /
rod.getDRot() / time;
rod.getDRot() / time ;
} break;
}
}
@ -99,4 +104,9 @@ double Compute::getTime() const { return static_cast<double>(every); }
double Compute::getTarget() const { return target; }
double Compute::getDifference() const { return abs(agg.getMean() - target); }
double Compute::getDifference() const { return abs(agg.getMean() - target); }
std::ostream &operator<<(std::ostream &os, const Compute &compute) {
os << "agg: " << compute.agg << " type: " << compute.type_str;
return os;
}

7
src/Compute.h
View File

@ -5,6 +5,7 @@
#ifndef MYPROJECT_COMPUTE_H
#define MYPROJECT_COMPUTE_H
#include <ostream>
#include "LiveAgg.hpp"
#include "Rod2d.hpp"
#include "Simulation.h"
@ -35,13 +36,17 @@ class Compute {
double getDifference() const;
private:
friend std::ostream &operator<<(std::ostream &os, const Compute &compute);
private:
Rod2d start_rod;
LiveAgg agg;
size_t every;
size_t time_step;
Type type;
double target;
double time;
std::string type_str;
};
#endif // MYPROJECT_COMPUTE_H

137
src/Integratoren2d_force.cpp
View File

@ -3,17 +3,10 @@
//
#include "Integratoren2d_force.h"
#include "vec_trafos.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;
}
using Matrix = Eigen::Matrix2d;
void Integratoren2d_force::Set0_deterministic(Rod2d &rod2D,
Simulation & /*sim*/) {
@ -25,31 +18,27 @@ void Integratoren2d_force::Set0_deterministic(Rod2d &rod2D,
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) {
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;
Vector trans_part = rod2D.getDiff_Sqrt() * rand;
Vector trans_lab = rotation_Matrix(rod2D.getE()) * trans_part;
// Force
Vector F_lab = force(rod2D.getPos(), rod2D.getE());
Vector F_part = rod2D.getE_Base_matrix().inverse() * F_lab;
Vector F_part = rotation_Matrix(rod2D.getE()).inverse() * F_lab;
Vector F_trans = rod2D.getDiff_Sqrt() * F_part;
F_trans *= sim.getMDeltaT() / kBT;
F_trans *= sim.getMDeltaT();
// 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();
double M_rot = torque(rod2D.getPos(), rod2D.getE()) * rod2D.getDRot() * sim.getMDeltaT();
Vector new_e = e + (rot + M_rot) * ortho(e);
new_e.normalize();
// apply
@ -66,12 +55,12 @@ Vector Integratoren2d_force::Heun_predictor_pos(
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 trans_pred_lab = rotation_Matrix(rod2D.getE()) * 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_part = rotation_Matrix(rod2D.getE()).inverse() * F_pred_lab;
Vector F_pred_trans = rod2D.getDiff_Sqrt() * F_pred_part;
F_pred_trans *= sim.getMDeltaT() / kBT;
F_pred_trans *= sim.getMDeltaT() ;
return F_pred_trans + trans_pred_lab;
}
@ -82,9 +71,9 @@ Vector Integratoren2d_force::Heun_predictor_rot(
auto std = sim.getSTD();
double rot_predict =
sim.getNorm(std) * rod2D.getDRot_Sqrt(); // rotationsmatrix verwenden.
Vector e = rod2D.getE();
const Vector& e = rod2D.getE();
double M_predict_rot = torque(rod2D.getPos(), rod2D.getE()) *
rod2D.getDRot() / kBT * sim.getMDeltaT();
rod2D.getDRot() * sim.getMDeltaT();
Vector e_change_predict = (rot_predict + M_predict_rot) * ortho(e);
return e_change_predict;
@ -98,30 +87,27 @@ void Integratoren2d_force::Set2_Heun(
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();
Vector e_predictor = rod2D.getE() + delta_e_predictor;
e_predictor.normalize();
Rod2d pred_rod = rod2D;
pred_rod.setPos(pos_predictor);
pred_rod.setE(e_predict);
pred_rod.setE(e_predictor);
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.setPos(pos_predictor);
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) {
const std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)>& force,
const 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
@ -129,20 +115,19 @@ void Integratoren2d_force::Set3_Exact(
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_part = rotation_Matrix(rod2D.getE()).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;
F_trans *= sim.getMDeltaT();
Vector trans_lab = rotation_Matrix(rod2D.getE()) * 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();
double M_rot = torque(rod2D.getPos(), rod2D.getE()) * rod2D.getDRot() * sim.getMDeltaT();
auto correction =
-0.5 * pow(rod2D.getDRot(), 2) / pow(kBT, 2) * sim.getMDeltaT();
-0.5 * pow(rod2D.getDRot(), 2) * sim.getMDeltaT();
Vector new_e = e + (rot + M_rot) * ortho(e) + correction * e;
new_e.normalize();
// apply
@ -152,30 +137,58 @@ void Integratoren2d_force::Set3_Exact(
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();
const std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)>& force,
const std::function<double(Eigen::Vector2d, Eigen::Vector2d)>& torque) {
auto std = sim.getSTD(); // sqrt(2*delta_T)
Vector e = rod2D.getE();
// 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 rand = {sim.getNorm(std), sim.getNorm(std)}; // gaussian noise
Vector trans_part = rod2D.getDiff_Sqrt() * rand;
Vector trans_lab = rotation_Matrix(rod2D.getE()) * trans_part;
// Force
Vector F_lab = force(rod2D.getPos(), rod2D.getE());
Vector F_part = rotMat.inverse() * F_lab;
Vector F_part = rotation_Matrix(rod2D.getE()).inverse() * F_lab;
Vector F_trans = rod2D.getDiff_Sqrt() * F_part;
F_trans *= sim.getMDeltaT() / kBT;
F_trans *= sim.getMDeltaT();
auto rot =
Eigen::Rotation2D<double>(sim.getNorm(std) * rod2D.getDRot_Sqrt());
Eigen::Rotation2Dd rotation(rot);
auto e_new = (rotation.toRotationMatrix() * rod2D.getE()).normalized();
// rotation
double rot =
sim.getNorm(std) * rod2D.getDRot_Sqrt(); // rotationsmatrix verwenden.
// 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);
double M_rot = torque(rod2D.getPos(), rod2D.getE()) * rod2D.getDRot() * sim.getMDeltaT();
Vector new_e = Eigen::Rotation2D<double>(rot + M_rot) * e;
new_e.normalize();
// apply
rod2D.setPos(rod2D.getPos() + trans_lab + F_trans);
rod2D.setE(new_e);
}
void Integratoren2d_force::Set5_MBD(Rod2d &rod2D, Simulation &sim,
const std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)>& force,
const std::function<double(Eigen::Vector2d, Eigen::Vector2d)>& torque) {
Vector delta_pos_predictor = Heun_predictor_pos(rod2D, sim, force);
[[maybe_unused]] Vector pos_predictor = rod2D.getPos() + delta_pos_predictor;
Vector delta_e_predictor = Heun_predictor_rot(rod2D, sim, torque);
Vector e_predictor = rod2D.getE() + delta_e_predictor;
e_predictor.normalize();
auto std = sim.getSTD();
Eigen::Vector2d e = rod2D.getE();
Eigen::Matrix2d Diff_combined = 0.5*(rod2D.getDiff() + rotation_Matrix(e).inverse() * rotation_Matrix(e_predictor)*rod2D.getDiff());
Eigen::Matrix2d D_combined_sqrt = Diff_combined.llt().matrixL();
auto rand = Eigen::Vector2d(sim.getNorm(std), sim.getNorm(std));
Eigen::Vector2d pos_integrated =rod2D.getPos() + rotation_Matrix(e) * D_combined_sqrt * rand;
double rot = sim.getNorm(std) * rod2D.getDRot_Sqrt();
Eigen::Vector2d e_integrated = e + 0.5 * (rot * ortho(e) + rot * ortho(e_predictor));
//apply
rod2D.setPos(pos_integrated);
rod2D.setE(e_integrated.normalized());
}

16
src/Integratoren2d_force.h
View File

@ -23,25 +23,23 @@ class Integratoren2d_force {
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);
const std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)>& force,
const 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);
const std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)>& force,
const std::function<double(Eigen::Vector2d, Eigen::Vector2d)>& torque);
static const Eigen::Vector2d unitVec;
static void Set5_MBD(Rod2d &rod2D, Simulation &sim);
static void Set5_MBD(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 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)>

139
src/Integratoren2d_forceless.cpp
View File

@ -3,10 +3,7 @@
//
#include "Integratoren2d_forceless.h"
Eigen::Vector2d Integratoren2d_forceless::ortho(Eigen::Vector2d e) {
return Eigen::Vector2d{-e[1], e[0]};
}
#include "vec_trafos.h"
void Integratoren2d_forceless::Set0_deterministic(Rod2d &rod2D,
Simulation & /*sim*/) {
@ -18,165 +15,109 @@ void Integratoren2d_forceless::Set0_deterministic(Rod2d &rod2D,
void Integratoren2d_forceless::Set1_Euler(Rod2d &rod2D, Simulation &sim) {
auto std = sim.getSTD(); // sqrt(2*delta_T)
Eigen::Vector2d e = rod2D.getE();
// translation
Eigen::Vector2d rand = {sim.getNorm(std),
sim.getNorm(std)}; // gaussian noise
Eigen::Vector2d trans_part =
rod2D.getDiff_Sqrt() * rand; // translation vector in Particle System
Eigen::Vector2d trans_lab = rod2D.getE_Base_matrix() * trans_part;
/* TODO
Statt Zeile 22 und 23 kannst du trans_lab = paralleler_matrix_eintrag * e
plus senkrechter_matrix_eintrag * ortho(e);
*/
Eigen::Vector2d rand = {sim.getNorm(std), sim.getNorm(std)};
Eigen::Vector2d trans_part = rod2D.getDiff_Sqrt() * rand;
Eigen::Vector2d trans_lab = rotation_Matrix(e) * trans_part;
// rotation
double rot =
sim.getNorm(std) * rod2D.getDRot_Sqrt(); // rotationsmatrix verwenden.
Eigen::Vector2d e = rod2D.getE();
double rot = sim.getNorm(std) * rod2D.getDRot_Sqrt();
Eigen::Vector2d new_e = e + rot * ortho(e);
new_e.normalize();
// apply
rod2D.setPos(rod2D.getPos() + trans_lab);
rod2D.setE(new_e);
rod2D.setE(new_e.normalized());
}
void Integratoren2d_forceless::Set2_Heun(Rod2d &rod2D, Simulation &sim) {
// Predict with Euler
auto std = sim.getSTD();
Eigen::Vector2d e = rod2D.getE();
Eigen::Vector2d rand = Eigen::Vector2d(sim.getNorm(std), sim.getNorm(std));
Eigen::Vector2d trans_part =
rod2D.getDiff_Sqrt() * rand; // translation vector in Particle System
Eigen::Vector2d trans_lab = rod2D.getE_Base_matrix() * trans_part;
rod2D.getDiff_Sqrt() * rand;
Eigen::Vector2d trans_lab = rotation_Matrix(e) * trans_part;
Eigen::Vector2d pos_predictor = rod2D.getPos() + trans_lab;
/*
TODO
Siehe Set_1: Hier auch über e und orth(e)
*/
// rotation
double rot_predict =
sim.getNorm(std) * rod2D.getDRot_Sqrt(); // rotationsmatrix verwenden.
Eigen::Vector2d e = rod2D.getE();
double rot_predict = sim.getNorm(std) * rod2D.getDRot_Sqrt();
Eigen::Vector2d delta_e_predict = rot_predict * ortho(e);
Eigen::Vector2d e_predict = (e + delta_e_predict).normalized();
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.setPos(pos_predictor); // pos with Euler
rod2D.setE(e_integrated.normalized());
}
constexpr double kBT = 1.0;
void Integratoren2d_forceless::Set3_Exact(Rod2d &rod2D, Simulation &sim) {
auto std = sim.getSTD();
Eigen::Vector2d e = rod2D.getE();
// translation
Eigen::Vector2d rand = {sim.getNorm(std),
sim.getNorm(std)}; // gaussian noise
Eigen::Vector2d trans_part =
rod2D.getDiff_Sqrt() * rand; // translation vector in Particle System
Eigen::Vector2d trans_lab = rod2D.getE_Base_matrix() * trans_part;
/*
TODO: Siehe Set1
*/
Eigen::Vector2d rand = {sim.getNorm(std), sim.getNorm(std)};
Eigen::Vector2d trans_part = rod2D.getDiff_Sqrt() * rand;
Eigen::Vector2d trans_lab = rotation_Matrix(e) * trans_part;
// rotation
double rot = sim.getNorm(std) * rod2D.getDRot_Sqrt();
Eigen::Vector2d e = rod2D.getE();
auto correction =
-0.5 * pow(rod2D.getDRot(), 2) / pow(kBT, 2) * sim.getMDeltaT();
-0.5 * pow(rod2D.getDRot(), 2) * sim.getMDeltaT();
Eigen::Vector2d delta_e = correction * e + rot * ortho(e);
/*
TODO Warum hier kBT definiert?
*/
// apply
rod2D.setPos(rod2D.getPos() + trans_lab);
rod2D.setE((e + delta_e).normalized());
}
const Eigen::Vector2d Integratoren2d_forceless::unitVec = {1, 0};
// this is a slow implementation to keep the same data structure for performance
// analysis this must be altered
// Normalisation should not be necessary if a proper angular representation
// is used. But with vector e it is done in case of numerical errors
void Integratoren2d_forceless::Set4_BDAS(Rod2d &rod2D, Simulation &sim) {
auto std = sim.getSTD();
Eigen::Vector2d e = rod2D.getE();
// translation
auto rand = Eigen::Vector2d(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;
/*
TODO: Warum hier nicht die Matrix von oben?
*/
Eigen::Vector2d trans_lab = rotation_Matrix(e) * trans_part;
auto rot =
Eigen::Rotation2D<double>(sim.getNorm(std) * rod2D.getDRot_Sqrt());
Eigen::Rotation2Dd rotation(rot);
auto e_new = (rotation.toRotationMatrix() * rod2D.getE()).normalized();
/*
TODO Warum Normierung? Ist der Fehler so gros? Wie siehts mit einer
Winkelvariable aus, dann musst du den Winkel nicht jedesmal berechnen?
*/
// Normalisation should not be necessary if a proper angular representation
// is used. But with vector e it is done in case of numerical errors
auto e_new = (rotation.toRotationMatrix() * e).normalized();
rod2D.setPos(rod2D.getPos() + trans_lab);
rod2D.setE(e_new);
}
Eigen::Matrix2d Integratoren2d_forceless::e_2_matrix(Eigen::Vector2d m_e) {
Eigen::Matrix2d mat;
mat << m_e[0], -m_e[1], m_e[1], m_e[0];
return mat;
}
void Integratoren2d_forceless::Set5_MBD(Rod2d & rod2D,
Simulation & sim) {
void Integratoren2d_forceless::Set5_MBD(Rod2d & /*rod2D*/,
Simulation & /*sim*/) {
/*
auto std = sim.getSTD();
auto rand = Eigen::Vector2d(sim.getNorm(std), sim.getNorm(std));
auto trans_part = rod2D.getDiff_Sqrt() * rand; //translation vector in
Particle System auto trans_lab = rod2D.getE_Base_matrix() * trans_part; auto
pos_predictor = rod2D.getPos() + trans_lab;
Eigen::Vector2d e = rod2D.getE();
//rotation
auto rot_predict = Eigen::Rotation2D<double>(sim.getNorm(std) *
rod2D.getDRot_Sqrt());// rotationsmatrix verwenden. auto e = rod2D.getE();
auto delta_e_predict = rot_predict * ortho(e); auto e_predict = (e +
delta_e_predict).normalized();
double rot_predict = sim.getNorm(std) *
rod2D.getDRot_Sqrt();
auto delta_e_predict = rot_predict * ortho(e);
auto e_predict = (e + delta_e_predict).normalized();
auto std_combined = rod2D.getDiff() + rod2D.getE_Base_matrix() *
rod2D.getDiff(); //old + new diff
auto rot = Eigen::Rotation2D<double>(sim.getNorm(std) * rod2D.getDRot_Sqrt());
auto e_integrated = e + 0.5 * (rot * ortho(e) + rot * ortho(e_predict));
Eigen::Matrix2d Diff_combined = 0.5*(rod2D.getDiff() + rotation_Matrix(e).inverse() * rotation_Matrix(e_predict)*rod2D.getDiff());
Eigen::Matrix2d D_combined_sqrt = Diff_combined.llt().matrixL();
auto rand = Eigen::Vector2d(sim.getNorm(std), sim.getNorm(std));
Eigen::Vector2d pos_integrated =rod2D.getPos() + rotation_Matrix(e) * D_combined_sqrt * rand;
double rot = sim.getNorm(std) * rod2D.getDRot_Sqrt();
Eigen::Vector2d e_integrated = e + 0.5 * (rot * ortho(e) + rot * ortho(e_predict));
//apply
rod2D.setPos(pos_predictor);
rod2D.setPos(pos_integrated);
rod2D.setE(e_integrated.normalized());
*/
}

5
src/Integratoren2d_forceless.h
View File

@ -19,12 +19,7 @@ class Integratoren2d_forceless {
static const Eigen::Vector2d unitVec;
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);
};

5
src/LiveAgg.cpp
View File

@ -19,3 +19,8 @@ void LiveAgg::feed(double value) noexcept {
S += delta * delta2;
}
int LiveAgg::getNumPoints() const noexcept { return num_of_data; }
std::ostream &operator<<(std::ostream &os, const LiveAgg &agg) {
os << "num_of_data: " << agg.num_of_data << " mean: " << agg.mean << " S: " << agg.S;
return os;
}

7
src/LiveAgg.hpp
View File

@ -1,4 +1,7 @@
#pragma once
#include <ostream>
class LiveAgg {
public:
LiveAgg();
@ -8,7 +11,9 @@ class LiveAgg {
[[nodiscard]] int getNumPoints() const noexcept;
void feed(double value) noexcept;
private:
friend std::ostream &operator<<(std::ostream &os, const LiveAgg &agg);
private:
int num_of_data;
double mean;
double S;

5
src/Rod2d.cpp
View File

@ -52,8 +52,3 @@ const Eigen::Vector2d &Rod2d::getE() const { return m_e; }
void Rod2d::setE(const Eigen::Vector2d &mE) { m_e = mE; }
Eigen::Matrix2d Rod2d::getE_Base_matrix() const {
Eigen::Matrix2d mat;
mat << m_e[0], -m_e[1], m_e[1], m_e[0];
return mat;
}

2
src/Rod2d.hpp
View File

@ -22,8 +22,6 @@ class Rod2d {
[[nodiscard]] double getDRot_Sqrt() const;
Eigen::Matrix2d getE_Base_matrix() const;
private:
Eigen::Matrix2d m_Diff;
Eigen::Matrix2d m_Diff_sqrt;

130
src/legacy/Integratoren.cpp
View File

@ -1,130 +0,0 @@
#include "Integratoren.hpp"
constexpr double kBT = 1.0;
void Integratoren::integrateITOEuler_1(Rod2d &rod)
{
//DLOG_SCOPE_F(1,"Integrateing Rod2d:");
double movePara = std::sqrt(2.0 * rod.Dpara * rod.deltaT) * rod.dis(rod.generator);
double moveOrtho = std::sqrt(2.0 * rod.Dortho * rod.deltaT) * rod.dis(rod.generator);
double tmpx1 = rod.x1 + movePara * rod.e1 - moveOrtho * rod.e2;
double tmpx2 = rod.x2 + movePara * rod.e2 + moveOrtho * rod.e1;
double rodOrtho = std::sqrt(2.0 * rod.Drot * rod.deltaT) * rod.dis(rod.generator);
double tmpe1 = rod.e1 - rodOrtho * rod.e2;
double tmpe2 = rod.e2 + rodOrtho * rod.e1;
rod.x1 = tmpx1;
rod.x2 = tmpx2;
rod.e1 = tmpe1;
rod.e2 = tmpe2;
rod.normalize();
}
void Integratoren::integrateITOHeun_2(Rod2d &rod)
{
//DLOG_SCOPE_F(1,"Integrateing Rod2d:");
double movePara = std::sqrt(2.0 * rod.Dpara * rod.deltaT) * rod.dis(rod.generator);
double moveOrtho = std::sqrt(2.0 * rod.Dortho * rod.deltaT) * rod.dis(rod.generator);
double tmpx1 = rod.x1 + movePara * rod.e1 - moveOrtho * rod.e2;
double tmpx2 = rod.x2 + movePara * rod.e2 + moveOrtho * rod.e1;
double rodOrtho = std::sqrt(2.0 * rod.Drot * rod.deltaT) * rod.dis(rod.generator);
double e1 = rod.e1 - rodOrtho * rod.e1;
double e2 = rod.e2 + rodOrtho * rod.e2;
double norm = std::sqrt(e1 * e1 + e2 * e2);
e1 /= norm;
e2 /= norm;
double rodOrtho1 = std::sqrt(2.0 * rod.Drot * rod.deltaT);
double rodOrthoRand = rod.dis(rod.generator);
double tmpe1 = rod.e1 - 0.5 * (rodOrtho1 * rod.e2 + rodOrtho1 * e2) * rodOrthoRand;//TODO Check if 0.5 faktor
double tmpe2 = rod.e2 + 0.5 * (rodOrtho1 * rod.e1 + rodOrtho1 * e1) * rodOrthoRand;
rod.x1 = tmpx1;
rod.x2 = tmpx2;
rod.e1 = tmpe1;
rod.e2 = tmpe2;
rod.normalize();
}
void Integratoren::integrateExact_3(Rod2d &rod)
{
//DLOG_SCOPE_F(1,"Integrateing Rod2d:");
double movePara = std::sqrt(2.0 * rod.Dpara * rod.deltaT) * rod.dis(rod.generator);
double moveOrtho = std::sqrt(2.0 * rod.Dortho * rod.deltaT) * rod.dis(rod.generator);
double tmpx1 = rod.x1 + movePara * rod.e1 - moveOrtho * rod.e2;
double tmpx2 = rod.x2 + movePara * rod.e2 + moveOrtho * rod.e1;
double rodOrtho = std::sqrt(2.0 * rod.Drot * rod.deltaT) * rod.dis(rod.generator);
double korr = -0.5 * rod.Drot * rod.Drot / kBT / kBT * rod.deltaT;
double tmpe1 = rod.e1 - rodOrtho * rod.e2 + korr * rod.e1;
double tmpe2 = rod.e2 + rodOrtho * rod.e1 + korr * rod.e2;
rod.x1 = tmpx1;
rod.x2 = tmpx2;
rod.e1 = tmpe1;
rod.e2 = tmpe2;
rod.normalize();
}
void Integratoren::integrateBDAS_4(Rod2d &rod)
{
double deltaRTpara = std::sqrt(2.0 * rod.Dpara * rod.deltaT) * rod.dis(rod.generator);
double deltaRTortho = std::sqrt(2.0 * rod.Dortho * rod.deltaT) * rod.dis(rod.generator);
double tmpx1 = rod.x1 + deltaRTpara * rod.e1 - deltaRTortho * rod.e2;
double tmpx2 = rod.x2 + deltaRTpara * rod.e2 + deltaRTortho * rod.e1;
double deltaPhi = std::sqrt(2.0 * rod.Drot * rod.deltaT) * rod.dis(rod.generator);
double tmpe1 = rod.e1 * std::cos(deltaPhi) - rod.e2 * std::sin(deltaPhi);
double tmpe2 = rod.e1 * std::sin(deltaPhi) + rod.e2 * std::cos(deltaPhi);
rod.x1 = tmpx1;
rod.x2 = tmpx2;
rod.e1 = tmpe1;
rod.e2 = tmpe2;
rod.normalize();
}
void Integratoren::integrateMBD_5(Rod2d &rod)
{
/*
double movePara = std::sqrt(2.0*rod.Dpara*rod.deltaT)*rod.dis(rod.generator); // Taken from euler
double moveOrtho = std::sqrt(2.0*rod.Dortho*rod.deltaT)*rod.dis(rod.generator); // taken from euler
double x1 = rod.x1 + movePara*rod.e1- moveOrtho*rod.e2;
double x2 = rod.x2 + movePara*rod.e2+ moveOrtho*rod.e1;
*/
// unused without force
double rodOrtho = std::sqrt(2.0 * rod.Drot * rod.deltaT) * rod.dis(rod.generator);
double e1 = rod.e1 - rodOrtho * rod.e2;
double e2 = rod.e2 + rodOrtho * rod.e1;
double norm = std::sqrt(e1 * e1 + e2 * e2);
e1 /= norm;
e2 /= norm;
double movePara1 = std::sqrt(1.0 * rod.Dpara * rod.deltaT) * rod.dis(rod.generator);// If difussion would be time dependend, average of the predicted D and now D
double moveOrtho1 = std::sqrt(1.0 * rod.Dortho * rod.deltaT) * rod.dis(rod.generator);
double movePara2 = std::sqrt(1.0 * rod.Dpara * rod.deltaT) * rod.dis(rod.generator);// If difussion would be time dependend, average of the predicted D and now D
double moveOrtho2 = std::sqrt(1.0 * rod.Dortho * rod.deltaT) * rod.dis(rod.generator);
double tmpx1 = rod.x1 + (movePara1 * rod.e1 - moveOrtho1 * rod.e2) + (movePara2 * e1 - moveOrtho2 * e2);
double tmpx2 = rod.x2 + (movePara1 * rod.e2 + moveOrtho1 * rod.e1) + (movePara2 * e2 + moveOrtho2 * e1);
double rodOrtho1 = std::sqrt(2.0 * rod.Drot * rod.deltaT);
double rodOrtho2 = rod.dis(rod.generator);
double tmpe1 = rod.e1 - 0.5 * (rodOrtho1 * rod.e2 + rodOrtho1 * e2) * rodOrtho2;
double tmpe2 = rod.e2 + 0.5 * (rodOrtho1 * rod.e1 + rodOrtho1 * e1) * rodOrtho2;
rod.x1 = tmpx1;
rod.x2 = tmpx2;
rod.e1 = tmpe1;
rod.e2 = tmpe2;
rod.normalize();
}

19
src/legacy/Integratoren.hpp
View File

@ -1,19 +0,0 @@
#ifndef INTEGRATOREN_HPP
#define INTEGRATOREN_HPP
#include "Rod2d.hpp"
class Integratoren
{
public:
static void integrateITOEuler_1(Rod2d &rod);
static void integrateITOHeun_2(Rod2d &rod);
static void integrateExact_3(Rod2d &rod);
static void integrateBDAS_4(Rod2d &rod);
static void integrateMBD_5(Rod2d &rod);
};
#endif// INTEGRATOREN_HPP

6
src/legacy/Simulation2d.cpp
View File

@ -1,6 +0,0 @@
#include "Simulation2d.hpp"
#include <vector>
#include<iostream>
#include "LiveAgg.hpp"
template <void (*T)(Rod2d &)>
std::vector<double> simulate(double deltaT, int seed, int meanSteps);

157
src/legacy/Simulation2d.hpp
View File

@ -1,157 +0,0 @@
#ifndef SIMULATION2D_HPP
#define SIMULATION2D_HPP
#include "Rod2d.hpp"
#include "LiveAgg.hpp"
#include <fstream>
#include <iostream>
/* This is used to calculate the convergence order of the particles
*/
inline void printProgress(double progress)
{
const int barWidth = 70;
std::cout << "[";
const int pos = static_cast<int>(barWidth * progress);
for (int i = 0; i < barWidth; ++i) {
if (i < pos) {
std::cout << "=";
} else if (i == pos) {
std::cout << ">";
} else {
std::cout << " ";
}
}
std::cout << "] " << int(progress * 100.0) << " % \r";
std::cout.flush();
}
template<void (*T)(Rod2d &)>
void simulate(double L, int seed, const std::string &filename)
{
// we now simulate the msd after 10 steps
Rod2d rod(L, 1e-5, seed);
const int sampleWidth = 1000;
std::vector<LiveAgg> msd(sampleWidth, LiveAgg());
std::vector<LiveAgg> oaf(sampleWidth, LiveAgg());
std::vector<LiveAgg> empxx(sampleWidth, LiveAgg());
std::vector<LiveAgg> empyy(sampleWidth, LiveAgg());
std::vector<double> xPos(sampleWidth, 0.0);
std::vector<double> yPos(sampleWidth, 0.0);
std::vector<double> e1(sampleWidth, 0.0);
std::vector<double> e2(sampleWidth, 0.0);
std::vector<double> targetMSD(sampleWidth, 0.0);
std::vector<double> targetOAF(sampleWidth, 0.0);
std::vector<double> targetEMPXX(sampleWidth, 0.0);
std::vector<double> targetEMPYY(sampleWidth, 0.0);
std::vector<size_t> samplePoints(sampleWidth, 1);
size_t temp = 1;
const double Dmean = 0.5 * (rod.Dpara + rod.Dortho);
const double deltaD = rod.Dortho - rod.Dpara;
for (size_t i = 0; i < samplePoints.size(); ++i) {
samplePoints[i] = temp;
temp += 100;
const double time = static_cast<double>(temp) * rod.deltaT;
targetMSD[i] = 4.0 * 0.5 * (rod.Dpara + rod.Dortho) * time;
targetOAF[i] = std::exp(-rod.Drot * time);
const double u = 4.0;
targetEMPXX[i] = Dmean - deltaD / 2.0 * (1 - exp(-u * rod.Drot * time)) / u / rod.Drot / time;
targetEMPYY[i] = Dmean + deltaD / 2.0 * (1 - exp(-u * rod.Drot * time)) / u / rod.Drot / time;
}
unsigned long long int iteration = 0;
//const unsigned long long int maxItr = ((long int)10000)*1000000;
const unsigned long long int maxItr = 10000000;
const unsigned long long int updateInter = 10000;
while (iteration < maxItr) {
T(rod);
iteration++;
for (size_t i = 0; i < samplePoints.size(); ++i) {
const size_t sample = samplePoints[i];
if (iteration % sample == 0) {
msd[i].feed((rod.x1 - xPos[i]) * (rod.x1 - xPos[i]) + (rod.x2 - yPos[i]) * (rod.x2 - yPos[i]));
oaf[i].feed(rod.e1 * e1[i] + rod.e2 * e2[i]);
empxx[i].feed(std::pow((rod.x1 - xPos[i]) * e1[i] + (rod.x2 - yPos[i]) * e2[i], 2));
empyy[i].feed(std::pow((rod.x1 - xPos[i]) * e2[i] - (rod.x2 - yPos[i]) * e1[i], 2));
xPos[i] = rod.x1;
yPos[i] = rod.x2;
e1[i] = rod.e1;
e2[i] = rod.e2;
}
}
if (iteration % updateInter == 0) {
printProgress(static_cast<double>(iteration)/ maxItr);
}
}
std::cout << "\n";
std::ofstream fileMSD(filename + "_" + std::to_string(L) + "_MSD.out");
std::ofstream fileOAF(filename + "_" + std::to_string(L) + "_OAF.out");
std::ofstream fileXX(filename + "_" + std::to_string(L) + "_XX.out");
std::ofstream fileYY(filename + "_" + std::to_string(L) + "_YY.out");
for (size_t i = 0; i < samplePoints.size(); ++i) {
fileMSD << static_cast<double>(samplePoints[i]) * rod.deltaT << "\t" << targetMSD[i] << "\t" << msd[i].getMean() << "\t" << msd[i].getSEM() << "\t" << msd[i].getSD() << "\t" << msd[i].getNumPoints() << "\n";
fileOAF << static_cast<double>(samplePoints[i]) * rod.deltaT << "\t" << targetOAF[i] << "\t" << oaf[i].getMean() << "\t" << oaf[i].getSEM() << "\t" << oaf[i].getSD() << "\t" << oaf[i].getNumPoints() << "\n";
fileXX << static_cast<double>(samplePoints[i]) * rod.deltaT << "\t" << targetEMPXX[i] << "\t" << empxx[i].getMean() << "\t" << empxx[i].getSEM() << "\t" << empxx[i].getSD() << "\t" << empxx[i].getNumPoints() << "\n";
fileYY << static_cast<double>(samplePoints[i]) * rod.deltaT << "\t" << targetEMPYY[i] << "\t" << empyy[i].getMean() << "\t" << empyy[i].getSEM() << "\t" << empyy[i].getSD() << "\t" << empyy[i].getNumPoints() << "\n";
}
fileMSD.close();
fileOAF.close();
fileXX.close();
fileYY.close();
}
template<void (*T)(Rod2d &)>
void konvergenz(double L, int seed, double deltaT, const std::string &filename)
{
// we now simulate the msd after 10 steps
Rod2d rod(L, deltaT, seed);
const double time = 2;
const int steps = static_cast<int>(time / deltaT);
const int minSteps = 1000;
//const double targetMSD = 4.0*0.5*(rod.Dpara+rod.Dortho)*time;
const double targetOAF = std::exp(-rod.Drot * time);
/*const double u = 4.0;
const double Dmean = 0.5*(rod.Dpara+rod.Dortho);
const double deltaD = rod.Dortho - rod.Dpara;
const double targetEMPXX = Dmean-deltaD/2.0*(1 - exp(-u*rod.Drot*time))/u/rod.Drot / time;
const double targetEMPYY = Dmean+deltaD/2.0*(1 - exp(-u*rod.Drot*time))/u/rod.Drot / time;*/
//LiveAgg msd;
LiveAgg oaf;
//LiveAgg empxx;
//LiveAgg empyy;
int meanCount = 0;
while (meanCount < minSteps or std::abs(targetOAF - oaf.getMean()) < oaf.getSEM() * 10) {
//reset the position
rod.reset();
for (int i = 0; i < steps; ++i) {
T(rod);
}
//msd.feed(rod.x1*rod.x1+rod.x2*rod.x2);
oaf.feed(rod.e1);
//.feed(rod.x1*rod.x1/2/time);
//empyy.feed(rod.x2*rod.x2/2/time);
if (meanCount % 100000 == 0) {
//std::cout<<msd.getSEM()<<" "<<std::abs(targetMSD-msd.getMean());
std::cout << "\t" << oaf.getSEM() << " " << std::abs(targetOAF - oaf.getMean()) << std::endl;
//std::cout<<"\t"<<empxx.getSEM()<<" "<<std::abs(targetEMPXX-empxx.getMean());
//std::cout<<"\t"<<empyy.getSEM()<<" "<<std::abs(targetEMPYY-empyy.getMean())<<std::endl;;
}
if (meanCount % 1000000 == 0) {
std::ofstream fileOAF(filename + "_length=" + std::to_string(L) + "_dt=" + std::to_string(rod.deltaT) + "_OAF.out");
fileOAF << rod.deltaT << "\t" << std::abs(targetOAF - oaf.getMean()) << "\t" << oaf.getSEM() << "\n";
fileOAF.close();
}
meanCount++;
}
std::ofstream fileOAF(filename + "_length=" + std::to_string(L) + "_dt=" + std::to_string(rod.deltaT) + "_OAF.out");
fileOAF << rod.deltaT << "\t" << std::abs(targetOAF - oaf.getMean()) << "\t" << oaf.getSEM() << "\n";
fileOAF.close();
}
#endif// SIMULATION2D_HPP

72
src/legacy/main.cpp
View File

@ -1,72 +0,0 @@
#include <random>
#include <span>
#include <string>
#include "legacy/Integratoren.hpp"
#include "legacy/Simulation2d.hpp"
const int seed = 12345;
[[maybe_unused]] void overTime(int sim)
{
const std::vector<double> Length = { 1.0, 1.3974, 20.0 };
const int numCases = 5;
int i = 0;
for (auto L : Length) {
if (sim == numCases * i) {
simulate<Integratoren::integrateITOEuler_1>(L, seed, "1");
}
if (sim == numCases * i + 1) {
simulate<Integratoren::integrateITOHeun_2>(L, seed, "2");
}
if (sim == numCases * i + 2) {
simulate<Integratoren::integrateExact_3>(L, seed, "3");
}
if (sim == numCases * i + 3) {
simulate<Integratoren::integrateBDAS_4>(L, seed, "4");
}
if (sim == numCases * i + 4) {
simulate<Integratoren::integrateMBD_5>(L, seed, "5");
}
i++;
}
}
[[maybe_unused]]void konvergenz(int sim)
{
const std::vector<double> Length = { 1.0, 1.3974, 20.0 };
const std::vector<double> DELTAT = { 2.0, 1.0, 0.5, 0.25, 0.2, 0.1, 0.05, 0.025, 0.02, 0.01, 0.005, 0.0025, 0.200, 0.001, 0.0005, 0.00025 };
int i = 0;
for (auto dt : DELTAT) {
for (auto L : Length) {
if (sim == i++) {
konvergenz<Integratoren::integrateITOEuler_1>(L, seed, dt, "1");
}
if (sim == i++) {
konvergenz<Integratoren::integrateITOHeun_2>(L, seed, dt, "2");
}
if (sim == i++) {
konvergenz<Integratoren::integrateExact_3>(L, seed, dt, "3");
}
if (sim == i++) {
konvergenz<Integratoren::integrateBDAS_4>(L, seed, dt, "4");
}
if (sim == i++) {
konvergenz<Integratoren::integrateMBD_5>(L, seed, dt, "5");
}
}
}
}
int main(int argc, char *argv[])
{
int sim;
auto args = std::span(argv, size_t(argc));
if (argc > 1) {
sim = std::stoi(args[1]) - 1;
} else {
sim = 0;
}
overTime(sim);
//konvergenz(sim);
return 0;
}

66
src/main.cpp
View File

@ -9,7 +9,7 @@
#include "Integratoren2d_forceless.h"
int main() {
constexpr int numStep = 1000000;
constexpr int numStep = 10000000;
constexpr double k = 0.01;
[[maybe_unused]] auto harmonic_Force =
[](const Eigen::Vector2d &pos,
@ -21,55 +21,31 @@ int main() {
const Eigen::Vector2d & /*torque*/) -> Eigen::Vector2d {
return {0.0, 0.0};
};
auto zero_Torque = [](const Eigen::Vector2d & /*pos*/,
[[maybe_unused]] 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::Set1_Euler(rod, sim, harmonic_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::msd, 32},
{Compute::Type::msd, 64},
{Compute::Type::msd, 128},
{Compute::Type::msd, 256},
{Compute::Type::msd, 512},
{Compute::Type::msd, 1024},
{Compute::Type::oaf, 1},
{Compute::Type::oaf, 2},
{Compute::Type::oaf, 4},
{Compute::Type::oaf, 8},
{Compute::Type::oaf, 16}},
0.01, 12345);
constexpr Compute::Type compute = Compute::Type::oaf;
Calculation euler(Integratoren2d_force::Set2_Heun,
{{compute, 1},
{compute, 2},
{compute, 4},
{compute, 8},
{compute, 16},
{compute, 32},
{compute, 64},
{compute, 128},
{compute, 256},
{compute, 512},
{compute, 1024}},
0.01, 12345, zero_Force,zero_Torque,1.0);
euler.run(numStep);
for (const auto &com : euler.getComputes()) {
if (com.getType() != Compute::Type::msd) continue;
std::cout << "MSD: " << com.getDifference() << " "
<< com.getAgg().getMean() << " <-> " << com.getTarget()
<< std::endl;
}
for (const auto &com : euler.getComputes()) {
if (com.getType() != Compute::Type::oaf) continue;
std::cout << "OAF: " << com.getDifference() << " "
<< com.getAgg().getMean() << " <-> " << com.getTarget()
if (com.getType() != compute) { continue;
}
std::cout
<< com.getAgg().getMean() << " , " << com.getTarget()
<< std::endl;
}
}

13
src/vec_trafos.h Normal file
View File

@ -0,0 +1,13 @@
//
// Created by jholder on 10/25/21.
//
#pragma once
#include <Eigen/Dense>
inline static Eigen::Vector2d ortho(Eigen::Vector2d e) {
return Eigen::Vector2d{-e[1], e[0]};
}
inline static Eigen::Matrix2d rotation_Matrix(Eigen::Vector2d e) {
Eigen::Matrix2d mat;
mat << e[0], -e[1], e[1], e[0];
return mat;
}

243
test/test_Integratoren.cpp
View File

@ -10,26 +10,202 @@
TEST_CASE("Forceless Integratoren") {
const size_t SEED = Catch::rngSeed();
auto force = [](const Eigen::Vector2d & /*pos*/,
const Eigen::Vector2d & /*torque*/) -> Eigen::Vector2d {
return {0.0, 0.0};
};
auto torque = [](const Eigen::Vector2d & /*pos*/,
const Eigen::Vector2d & /*torque*/) -> double {
const double length = GENERATE(1.0,1.4,2.0);
constexpr int steps = 10000;
constexpr double delta = 0.1;
Calculation euler(Integratoren2d_forceless::Set1_Euler,
{{Compute::Type::msd, 1},
{Compute::Type::oaf, 1},
{Compute::Type::empxx, 1},
{Compute::Type::empyy, 1}},
0.01, SEED,length);
euler.run(steps);
Calculation heun(Integratoren2d_forceless::Set2_Heun,
{{Compute::Type::msd, 1},
{Compute::Type::oaf, 1},
{Compute::Type::empxx, 1},
{Compute::Type::empyy, 1}},
0.01, SEED,length);
heun.run(steps);
Calculation exact(Integratoren2d_forceless::Set3_Exact,
{{Compute::Type::msd, 1},
{Compute::Type::oaf, 1},
{Compute::Type::empxx, 1},
{Compute::Type::empyy, 1}},
0.01, SEED,length);
exact.run(steps);
Calculation bdas(Integratoren2d_forceless::Set4_BDAS,
{{Compute::Type::msd, 1},
{Compute::Type::oaf, 1},
{Compute::Type::empxx, 1},
{Compute::Type::empyy, 1}},
0.01, SEED,length);
bdas.run(steps);
Calculation mbd(Integratoren2d_forceless::Set5_MBD,
{{Compute::Type::msd, 1},
{Compute::Type::oaf, 1},
{Compute::Type::empxx, 1},
{Compute::Type::empyy, 1}},
0.01, SEED,length);
mbd.run(steps);
SECTION("ForceLess Euler") {
CAPTURE(length);
size_t i = 0;
for (auto &c: euler.getComputes()) {
CAPTURE(c);
CHECK(c.getDifference() <= delta*c.getAgg().getMean());
++i;
}
}
SECTION("ForceLess Heun") {
CAPTURE(length);
size_t i = 0;
for (auto &c: heun.getComputes()) {
CAPTURE(c);
CHECK(c.getDifference() <= delta*c.getAgg().getMean());
++i;
}
}
SECTION("ForceLess Exact") {
CAPTURE(length);
size_t i = 0;
for (auto &c: exact.getComputes()) {
CAPTURE(c);
CHECK(c.getDifference() <= delta*c.getAgg().getMean());
++i;
}
}
SECTION("ForceLess BDAS") {
CAPTURE(length);
size_t i = 0;
for (auto &c: bdas.getComputes()) {
CAPTURE(c);
CHECK(c.getDifference() <= delta*c.getAgg().getMean());
++i;
}
}
SECTION("ForceLess MBD") {
CAPTURE(length);
size_t i = 0;
for (auto &c: mbd.getComputes()) {
CAPTURE(c);
CHECK(c.getDifference() <= delta*c.getAgg().getMean());
++i;
}
}
[[maybe_unused]] auto zero_Force =
[](const Eigen::Vector2d & /*pos*/,
const Eigen::Vector2d & /*torque*/) -> Eigen::Vector2d {
return {0.0, 0.0};
};
[[maybe_unused]] auto zero_Torque = [](const Eigen::Vector2d & /*pos*/,
const Eigen::Vector2d & /*torque*/) -> double {
return 0.0;
};
Calculation euler_zero(Integratoren2d_force::Set1_Euler,
{{Compute::Type::msd, 1},
{Compute::Type::oaf, 1},
{Compute::Type::empxx, 1},
{Compute::Type::empyy, 1}},
0.01, SEED,zero_Force,zero_Torque,length);
euler_zero.run(steps);
Calculation heun_zero(Integratoren2d_force::Set2_Heun,
{{Compute::Type::msd, 1},
{Compute::Type::oaf, 1},
{Compute::Type::empxx, 1},
{Compute::Type::empyy, 1}},
0.01, SEED,zero_Force,zero_Torque,length);
heun_zero.run(steps);
Calculation exact_zero(Integratoren2d_force::Set3_Exact,
{{Compute::Type::msd, 1},
{Compute::Type::oaf, 1},
{Compute::Type::empxx, 1},
{Compute::Type::empyy, 1}},
0.01, SEED,zero_Force,zero_Torque,length);
exact_zero.run(steps);
Calculation bdas_zero(Integratoren2d_force::Set4_BDAS,
{{Compute::Type::msd, 1},
{Compute::Type::oaf, 1},
{Compute::Type::empxx, 1},
{Compute::Type::empyy, 1}},
0.01, SEED,zero_Force,zero_Torque,length);
bdas_zero.run(steps);
Calculation mbd_zero(Integratoren2d_force::Set5_MBD,
{{Compute::Type::msd, 1},
{Compute::Type::oaf, 1},
{Compute::Type::empxx, 1},
{Compute::Type::empyy, 1}},
0.01, SEED,zero_Force,zero_Torque,length);
mbd_zero.run(steps);
SECTION("Force Euler") {
CAPTURE(length);
size_t i = 0;
for (auto &c: euler_zero.getComputes()) {
CAPTURE(c);
CHECK(c.getDifference() <= delta*c.getAgg().getMean());
++i;
}
}
SECTION("Force Heun") {
CAPTURE(length);
size_t i = 0;
for (auto &c: heun_zero.getComputes()) {
CAPTURE(c);
CHECK(c.getDifference() <= delta*c.getAgg().getMean());
++i;
}
}
SECTION("Force Exact") {
CAPTURE(length);
size_t i = 0;
for (auto &c: exact_zero.getComputes()) {
CAPTURE(c);
CHECK(c.getDifference() <= delta*c.getAgg().getMean());
++i;
}
}
SECTION("Force BDAS") {
CAPTURE(length);
size_t i = 0;
for (auto &c: bdas_zero.getComputes()) {
CAPTURE(c);
CHECK(c.getDifference() <= delta*c.getAgg().getMean());
++i;
}
}
SECTION("Force MBD") {
CAPTURE(length);
size_t i = 0;
for (auto &c: mbd_zero.getComputes()) {
CAPTURE(c);
CHECK(c.getDifference() <= delta*c.getAgg().getMean());
++i;
}
}
}
/*
SECTION("Euler") {
Calculation euler(Integratoren2d_forceless::Set1_Euler,
{{Compute::Type::msd, 1},
{Compute::Type::oaf, 1},
{Compute::Type::empxx, 1},
{Compute::Type::empyy, 1}},
0.01, SEED);
{
euler.run(10000);
for (auto &c : euler.getComputes()) {
CHECK(c.getDifference() <= 0.005);
}
}
Calculation euler_force(
@ -143,4 +319,39 @@ TEST_CASE("Forceless Integratoren") {
Approx(bdas_force.getComputes()[1].getAgg().getMean())
.epsilon(10e-10));
}
SECTION("MBD") {
Calculation mbd(Integratoren2d_forceless::Set5_MBD,
{{Compute::Type::msd, 1},
{Compute::Type::oaf, 1},
{Compute::Type::empxx, 1},
{Compute::Type::empyy, 1}},
0.01, SEED);
{
mbd.run(10000);
for (auto &c : mbd.getComputes()) {
CHECK(c.getDifference() <= 0.005);
}
}
Calculation mbd_force(Integratoren2d_force::Set5_MBD,
{{Compute::Type::msd, 1},
{Compute::Type::oaf, 1},
{Compute::Type::empxx, 1},
{Compute::Type::empyy, 1}},
0.01, SEED, force, torque);
{
mbd_force.run(10000);
for (auto &c : mbd_force.getComputes()) {
CHECK(c.getDifference() <= 0.005);
}
}
CHECK(mbd.getComputes()[0].getAgg().getMean() ==
Approx(mbd_force.getComputes()[0].getAgg().getMean())
.epsilon(10e-10));
CHECK(mbd.getComputes()[1].getAgg().getMean() ==
Approx(mbd_force.getComputes()[1].getAgg().getMean())
.epsilon(10e-10));
}
}
*/

40
test/test_Rod2d.cpp
View File

@ -6,7 +6,7 @@
#include <catch2/catch.hpp>
#include "Rod2d.hpp"
#include "vec_trafos.h"
TEST_CASE("Rods") {
Rod2d sphere(1);
Rod2d rod(2);
@ -21,41 +21,13 @@ TEST_CASE("Rods") {
REQUIRE(sphere.getDiff_Sqrt()(1, 1) == 1.0);
}
SECTION("Checking rotation") {
SECTION("Forwards == e 90") {
//
// =====
//
rod.setE({0, 1}); // Lab system is 90 degree rotated
double f1 = GENERATE(take(10, random(-100, 100)));
double f2 = GENERATE(take(10, random(-100, 100)));
rod.setE(Eigen::Vector2d({f1, f2}).normalized());
auto test = Eigen::Vector2d({1, 0});
auto erg = (rod.getE_Base_matrix() * test).normalized();
auto e = rod.getE();
auto erg = (rotation_Matrix(e) * test).normalized();
REQUIRE(erg.isApprox(e));
}
SECTION("Forwards == e 45") {
// o
// o
// o
rod.setE(Eigen::Vector2d({1, 1})
.normalized()); // Lab system is 90 degree rotated
auto test = Eigen::Vector2d({1, 0});
auto erg = (rod.getE_Base_matrix() * test).normalized();
auto e = rod.getE();
REQUIRE(erg.isApprox(e));
}
}
SECTION("Angle Stuff") {
const Eigen::Vector2d unitVec = {1, 0};
double f1 = GENERATE(take(10, random(-100, 100)));
double f2 = GENERATE(take(10, random(0, 100)));
sphere.setE(Eigen::Vector2d({f1, f2}).normalized());
Eigen::Rotation2D rotMat(acos(unitVec.dot(sphere.getE())));
CHECK(rotMat.toRotationMatrix()(0, 0) ==
Approx(sphere.getE_Base_matrix()(0, 0)).epsilon(10e-10));
CHECK(rotMat.toRotationMatrix()(1, 0) ==
Approx(sphere.getE_Base_matrix()(1, 0)).epsilon(10e-10));
CHECK(rotMat.toRotationMatrix()(0, 1) ==
Approx(sphere.getE_Base_matrix()(0, 1)).epsilon(10e-10));
CHECK(rotMat.toRotationMatrix()(1, 1) ==
Approx(sphere.getE_Base_matrix()(1, 1)).epsilon(10e-10));
}
}