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

A lot of refactoring

Browse Source
This commit is contained in:
Jacob Holder 2021-10-27 16:11:55 +02:00
parent 842957bbba
commit 1997edeae7
Signed by: jacob
GPG Key ID: 2194FC747048A7FD
25 changed files with 510 additions and 486 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
.clang-tidy
View File

@ -1,5 +1,5 @@
---
Checks: '*,-fuchsia-*,-google-*,-zircon-*,-abseil-*,-modernize-use-trailing-return-type,-llvm*,-readability-magic-numbers,-cppcoreguidelines-avoid-magic-numbers'
Checks: '*,-fuchsia-*,-google-*,-zircon-*,-abseil-*,-modernize-use-trailing-return-type,-llvm*'
WarningsAsErrors: ''
HeaderFilterRegex: ''
FormatStyle: none

2
benchmark/bench.cpp
View File

@ -6,7 +6,7 @@
#define CATCH_CONFIG_ENABLE_BENCHMARKING
#include <catch2/catch.hpp>
#include "Integratoren2d_forceless.h"
#include "legacy/Integratoren2d_forceless.h"
#include "Rod2d.hpp"
TEST_CASE("Euler - Baseline", "[benchmark]") {

2
benchmark/bench_euler.cpp
View File

@ -3,7 +3,7 @@
//
#include "Rod2d.hpp"
#include "Simulation.h"
#include "Integratoren2d_forceless.h"
#include "legacy/Integratoren2d_forceless.h"
int main(){
Rod2d rod(1.0);
Simulation sim(0.01, 1);

2
benchmark/bench_mbd.cpp
View File

@ -4,7 +4,7 @@
#include "Rod2d.hpp"
#include "Simulation.h"
#include "Integratoren2d_forceless.h"
#include "legacy/Integratoren2d_forceless.h"
int main(){
Rod2d rod(1.0);
Simulation sim(0.01, 1);

2
src/CMakeLists.txt
View File

@ -4,7 +4,7 @@ find_package(Eigen3)
# Generic test that uses conan libs
add_library(integratoren SHARED LiveAgg.cpp Rod2d.cpp Simulation.cpp
Simulation.h Integratoren2d_forceless.cpp Integratoren2d_force.cpp
Simulation.h Integratoren2d.cpp
Calculation.cpp Compute.cpp)
target_include_directories(integratoren PUBLIC .)
target_link_libraries(

17
src/Calculation.cpp
View File

@ -21,7 +21,7 @@ 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, double length)
: rod(length), sim(deltaT, seed), m_integrator(std::move(t_integrator)) {
: rod(length), start_rod(rod), 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));
}
@ -30,7 +30,7 @@ Calculation::Calculation(
Calculation::Calculation(inte_force_type t_integrator,
std::initializer_list<std::pair<Compute::Type, size_t>> t_computes,
double deltaT, size_t seed, force_type force, torque_type torque, double length)
: rod(length), sim(deltaT, seed) {
: rod(length), start_rod(rod), sim(deltaT, seed) {
m_integrator = [&](Rod2d &t_rod, Simulation &t_sim) {
t_integrator(t_rod, t_sim, force, torque);
};
@ -46,7 +46,7 @@ const std::vector<Compute> &Calculation::getComputes() const {
const Simulation &Calculation::getSim() const { return sim; }
Calculation::Calculation(std::function<void(Rod2d &, Simulation &)> t_integrator, const std::vector<std::pair<Compute::Type,size_t>>& t_computes,
double deltaT, size_t seed, double length) : rod(length), sim(deltaT, seed),
double deltaT, size_t seed, double length) : rod(length), start_rod(rod), 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));
@ -55,13 +55,18 @@ Calculation::Calculation(std::function<void(Rod2d &, Simulation &)> t_integrator
Calculation::Calculation(Calculation::inte_force_type t_integrator, std::vector<std::pair<Compute::Type,size_t>> t_computes, double deltaT,
size_t seed, Calculation::force_type force, Calculation::torque_type torque, double length)
: rod(length), sim(deltaT, seed) {
: rod(length), start_rod(rod), sim(deltaT, seed) {
m_integrator = [&](Rod2d &t_rod, Simulation &t_sim) {
t_integrator(t_rod, t_sim, force, torque);
};
for (const auto &pair: t_computes) {
computes.emplace_back(Compute(rod, pair.first, pair.second, sim));
}
rod.setPos({1000,1000});
}
void Calculation::reset() {
rod =start_rod;
for(auto & c:computes){
c.reset(start_rod);
}
}

4
src/Calculation.h
View File

@ -12,11 +12,15 @@
class Calculation {
private:
Rod2d rod;
const Rod2d start_rod;
Simulation sim;
std::function<void(Rod2d &, Simulation &)> m_integrator;
std::vector<Compute> computes;
public:
void reset();
[[nodiscard]] const Simulation &getSim() const;
using force_type = std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)>;

65
src/Compute.cpp
View File

@ -17,6 +17,13 @@ void Compute::evalMSD(const Rod2d &rod2D) {
agg.feed(msd);
}
void Compute::evalX(const Rod2d &rod2D) {
const auto &new_pos = rod2D.getPos();
auto old_pos = start_rod.getPos();
auto msd = (new_pos - old_pos).norm();
agg.feed(msd);
}
void Compute::evalOAF(const Rod2d &rod2D) {
const auto &new_e = rod2D.getE();
auto old_e = start_rod.getE();
@ -57,51 +64,86 @@ void Compute::eval(const Rod2d &rod2D) {
case Type::empyy:
eval_empYY(rod2D);
break;
case Type::msd_harmonic_k_10:
case Type::msd_harmonic_k_1:
evalMSD(rod2D);
break;
case Type::msd_const:
break;
case Type::x:
evalX(rod2D);
break;
case Type::x_squared:
evalMSD(rod2D);
break;
}
start_rod = rod2D;
if (resetting) {
start_rod = rod2D;
}
return;
}
}
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;
type_str = "msd";
target = 4.0 * 0.5 * (rod.getDiff().trace()) * time;
target = 2 * (rod.getDiff().trace()) * time;
break;
}
case Type::oaf: {
resetting = true;
type_str = "oaf";
target = std::exp(-rod.getDRot() * time);
break;
}
case Type::empxx: {
resetting = true;
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);
target = Dmean - deltaD / 2.0 *
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.0 *
target = Dmean + deltaD / 2 *
(1 - exp(-u * rod.getDRot() * time)) / u /
rod.getDRot() / time;
}
break;
case Type::msd_harmonic_k_10: {
case Type::msd_harmonic_k_1: {
resetting = true;
type_str = "msd_harmonic";
target = 2.0/10.0*(1.0 - std::exp(-20 * time));
target = 2 * (1.0 - std::exp(-2 * time));
}
break;
case Type::msd_const: {
resetting = true;
type_str = "msd_const";
target = -1.0;
}
break;
case Type::x: {
resetting = false;
type_str = "<x>";
target = -1.0;
}
break;
case Type::x_squared: {
resetting = false;
type_str = "<x_squared>";
target = -1.0;
}
break;
}
@ -121,3 +163,8 @@ std::ostream &operator<<(std::ostream &os, const Compute &compute) {
os << "agg: " << compute.agg << " type: " << compute.type_str;
return os;
}
void Compute::reset(const Rod2d &rod) {
time_step = 0;
start_rod = rod;
}

15
src/Compute.h
View File

@ -12,7 +12,7 @@
class Compute {
public:
enum class Type { msd, oaf, empxx, empyy, msd_harmonic_k_10 };
enum class Type { msd, oaf, empxx, empyy, msd_harmonic_k_1,msd_const, x,x_squared };
explicit Compute(Rod2d rod, Type t_type, size_t t_every, Simulation &sim);
@ -24,7 +24,7 @@ class Compute {
void eval_empXX(const Rod2d &rod2D);
double getTarget() const;
[[nodiscard]] double getTarget() const;
void eval_empYY(const Rod2d &rod2D);
@ -32,12 +32,14 @@ class Compute {
[[nodiscard]] Type getType() const;
double getTime() const;
[[nodiscard]] double getTime() const;
double getDifference() const;
[[nodiscard]] double getDifference() const;
friend std::ostream &operator<<(std::ostream &os, const Compute &compute);
void reset(const Rod2d &rod);
private:
Rod2d start_rod;
LiveAgg agg;
@ -47,6 +49,11 @@ private:
double target;
double time;
std::string type_str;
void evalX(const Rod2d &rod2D);
bool resetting;
};
#endif // MYPROJECT_COMPUTE_H

205
src/Integratoren2d.cpp Normal file
View File

@ -0,0 +1,205 @@
//
// Created by jholder on 21.10.21.
//
#include "Integratoren2d.h"
#include "vec_trafos.h"
using Vector = Eigen::Vector2d;
using Matrix = Eigen::Matrix2d;
Vector
force_euler(const Rod2d &rod2D, const Simulation &sim, const std::function<Vector(Vector, Vector)> &force) {
if (not force) {
return {0.0, 0.0};
}
const auto &e = rod2D.getE();
Vector F_lab = force(rod2D.getPos(), e);
Vector F_part = rotation_Matrix(e).inverse() * F_lab;
Vector F_trans = rotation_Matrix(e) * rod2D.getDiff_Sqrt() * F_part;
F_trans *= sim.getMDeltaT();
return F_trans;
}
Vector torque_euler(const Rod2d &rod2D, const Simulation &sim, const std::function<double(Vector, Vector)> &torque) {
if (not torque) {
return {0.0, 0.0};
}
const auto &e = rod2D.getE();
return torque(rod2D.getPos(), e) * rod2D.getDRot() * sim.getMDeltaT() * ortho(rod2D.getE());
}
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
Vector trans_part = rod2D.getDiff_Sqrt() * rand;
return rotation_Matrix(e) * trans_part;
}
Vector rand_rot_euler(const Rod2d &rod2D, Simulation &sim) {
const auto std = sim.getSTD();
return sim.getNorm(std) * rod2D.getDRot_Sqrt() * ortho(rod2D.getE());
}
void Integratoren2d::Set1_Euler_f(
Rod2d &rod2D, Simulation &sim,
const std::function<Vector(Vector, Vector)> &force,
const std::function<double(Vector, Vector)> &torque) {
// Gaussian - translation
auto trans = rand_trans_euler(rod2D, sim);
// Force
auto F_trans = force_euler(rod2D, sim, force);
// Gaussian rotation
auto rot = rand_rot_euler(rod2D, sim);
// torque
auto M_rot = torque_euler(rod2D, sim, torque);
Vector integrated_e = rod2D.getE() + rot + M_rot;
integrated_e.normalize();
Vector integrated_pos = rod2D.getPos() + trans + F_trans;
// apply
rod2D.setPos(integrated_pos);
rod2D.setE(integrated_e);
}
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<Vector(Vector, Vector)> &force,
const std::function<double(Vector, Vector)> &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
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(rod2D, sim) + torque_euler(rod2D, sim, torque);
//Heun integration
Vector e_integrated =
rod2D.getE() + 0.5 * (delta_e_predictor + delta_e_future);
// apply
rod2D.setPos(pos_predictor);
rod2D.setE(e_integrated.normalized());
}
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<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)> &force,
const std::function<double(Eigen::Vector2d, Eigen::Vector2d)> &torque) {
// Gaussian - translation
auto trans = rand_trans_euler(rod2D, sim);
// Force
auto F_trans = force_euler(rod2D, sim, force);
Vector integrated_pos = rod2D.getPos() + trans + F_trans;
// Gaussian rotation
auto rot = rand_rot_euler(rod2D, sim);
// torque
auto M_rot = torque_euler(rod2D, sim, torque);
auto correction =
-0.5 * pow(rod2D.getDRot(), 2) * sim.getMDeltaT() * rod2D.getE();
Vector integrated_e = rod2D.getE() + rot + M_rot + correction;
integrated_e.normalize();
// apply
rod2D.setPos(integrated_pos);
rod2D.setE(integrated_e);
}
void Integratoren2d::Set3_Exact(Rod2d &rod2D, Simulation &sim) {
Set3_Exact_f(rod2D, sim, {}, {});
}
void Integratoren2d::Set4_BDAS_f(
Rod2d &rod2D, Simulation &sim,
const std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)> &force,
const std::function<double(Eigen::Vector2d, Eigen::Vector2d)> &torque) {
// Gaussian - translation
auto trans = rand_trans_euler(rod2D, sim);
// Force
auto F_trans = force_euler(rod2D, sim, force);
Vector integrated_pos = rod2D.getPos() + trans + F_trans;
// 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;
Vector integrated_e = Eigen::Rotation2D<double>(rot + M_rot) * rod2D.getE();
integrated_e.normalize();
// apply
rod2D.setPos(integrated_pos);
rod2D.setE(integrated_e);
}
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<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)> &force,
const std::function<double(Eigen::Vector2d, Eigen::Vector2d)> &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
Rod2d pred_rod = rod2D;
pred_rod.setPos(pos_predictor);
pred_rod.setE(e_predictor);
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)));
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 pos_integrated = rod2D.getPos() + D_combined_sqrt * rand + 0.5*(force_euler(rod2D,sim,force)+force_euler(pred_rod,sim,force));
//apply
rod2D.setPos(pos_integrated);
rod2D.setE(e_integrated.normalized());
}
void Integratoren2d::Set5_MBD(Rod2d &rod2D, Simulation &sim) {
Set5_MBD_f(rod2D, sim, {}, {});
}

44
src/Integratoren2d.h Normal file
View File

@ -0,0 +1,44 @@
//
// Created by jholder on 21.10.21.
//
#pragma once
#include "Rod2d.hpp"
#include "Simulation.h"
class Integratoren2d {
using Vector = Eigen::Vector2d;
public:
static void Set1_Euler_f(Rod2d &rod2D, Simulation &sim, const std::function<Vector(Vector, Vector)> &force,
const std::function<double(Vector, Vector)> &torque);
static void Set1_Euler(Rod2d &rod2D, Simulation &sim);
static void Set2_Heun_f(Rod2d &rod2D, Simulation &sim, const std::function<Vector(Vector, Vector)> &force,
const std::function<double(Vector, Vector)> &torque);
static void Set2_Heun(Rod2d &rod2D, Simulation &sim);
static void
Set3_Exact_f(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 Set3_Exact(Rod2d &rod2D, Simulation &sim);
static void
Set4_BDAS_f(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 Set4_BDAS(Rod2d &rod2D, Simulation &sim);
static void
Set5_MBD_f(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 Set5_MBD(Rod2d &rod2D, Simulation &sim);
};

5
src/Rod2d.cpp
View File

@ -31,11 +31,6 @@ Rod2d::Rod2d(double L) : m_pos({0, 0}), m_e({1, 0}) {
sqrt(m_Diff_sqrt);
}
void Rod2d::reset() {
m_pos = {0, 0};
m_e = {1, 0};
}
void Rod2d::setPos(const Eigen::Vector2d &Pos) { m_pos = Pos; }
double Rod2d::getDRot() const { return m_D_rot; }

1
src/Rod2d.hpp
View File

@ -4,7 +4,6 @@
class Rod2d {
public:
explicit Rod2d(double L);
void reset();
[[nodiscard]] const Eigen::Matrix2d &getDiff() const;

25
src/force_lambdas.h Normal file
View File

@ -0,0 +1,25 @@
//
// Created by jholder on 27.10.21.
//
#pragma once
[[maybe_unused]] const static auto harmonic_Force =
[](const Eigen::Vector2d &pos,
const Eigen::Vector2d & /*torque*/) -> Eigen::Vector2d {
return -1.0 * pos;
};
[[maybe_unused]] const static auto const_Force =
[](const Eigen::Vector2d & /*pos*/,
const Eigen::Vector2d & /*torque*/) -> Eigen::Vector2d {
return {1,1};
};
[[maybe_unused]] const static auto zero_Force =
[](const Eigen::Vector2d & /*pos*/,
const Eigen::Vector2d & /*torque*/) -> Eigen::Vector2d {
return {0.0, 0.0};
};
[[maybe_unused]] const static auto zero_Torque = [](const Eigen::Vector2d & /*pos*/,
const Eigen::Vector2d & /*torque*/) -> double {
return 0.0;
};

11
src/Integratoren2d_force.cpp → src/legacy/Integratoren2d_force.cpp
View File

@ -8,14 +8,6 @@
using Vector = Eigen::Vector2d;
using Matrix = Eigen::Matrix2d;
void Integratoren2d_force::Set0_deterministic(Rod2d &rod2D,
Simulation & /*sim*/) {
auto trans_lab = Vector({0, 0.01});
rod2D.setPos(rod2D.getPos() + trans_lab);
Eigen::Rotation2D<double> rot(0.1);
rod2D.setE(rot * rod2D.getE());
}
const Vector Integratoren2d_force::unitVec = {1, 0};
void Integratoren2d_force::Set1_Euler(
@ -36,8 +28,7 @@ void Integratoren2d_force::Set1_Euler(
F_trans *= sim.getMDeltaT();
// rotation
double rot =
sim.getNorm(std) * rod2D.getDRot_Sqrt(); // rotationsmatrix verwenden.
double rot = sim.getNorm(std) * rod2D.getDRot_Sqrt();
double M_rot = torque(rod2D.getPos(), e) * rod2D.getDRot() * sim.getMDeltaT();
Vector new_e = e + (rot + M_rot) * ortho(e);
new_e.normalize();

7
src/Integratoren2d_force.h → src/legacy/Integratoren2d_force.h
View File

@ -37,9 +37,7 @@ class Integratoren2d_force {
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:
private:
static Eigen::Vector2d Heun_predictor_pos(
const Rod2d &rod2D, Simulation &sim,
const std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)>
@ -48,4 +46,7 @@ class Integratoren2d_force {
static Eigen::Vector2d Heun_predictor_rot(
const Rod2d &rod2D, Simulation &sim,
const std::function<double(Eigen::Vector2d, Eigen::Vector2d)> &torque);
void Set1_Euler(Rod2d &rod2D, Simulation &sim, const std::function<Vector(Vector, Vector)> &force,
const std::function<double(Vector, Vector)> &torque);
};

0
src/Integratoren2d_forceless.cpp → src/legacy/Integratoren2d_forceless.cpp
View File

0
src/Integratoren2d_forceless.h → src/legacy/Integratoren2d_forceless.h
View File

152
src/main.cpp
View File

@ -6,60 +6,114 @@
#include <fstream>
#include "Calculation.h"
#include "Integratoren2d_force.h"
#include "Integratoren2d_forceless.h"
#include "Integratoren2d.h"
#include "force_lambdas.h"
int main() {
constexpr int numStep = 10000000;
constexpr double k = 10;
[[maybe_unused]] auto harmonic_Force =
[](const Eigen::Vector2d &pos,
const Eigen::Vector2d & /*torque*/) -> Eigen::Vector2d {
return -k * pos;
};
[[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;
};
std::vector<std::pair<Compute::Type, size_t>> computes;
for (int i = 1; i < 100; ++i) {
computes.emplace_back(Compute::Type::msd_harmonic_k_10, 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);
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;
}
}
}
Calculation euler(Integratoren2d_force::Set4_BDAS,
computes, 0.01, 12345, harmonic_Force, zero_Torque, 1.0);
std::cout << "Finished run " << folder << std::endl;
}
euler.run(numStep);
std::ofstream msdFile("msd.dat");
std::ofstream msd_harmFile("msd_harm.dat");
std::ofstream oafFile("oaf.dat");
std::ofstream empxxFile("empxx.dat");
std::ofstream empyyFile("empyy.dat");
for (const auto &com: euler.getComputes()) {
if (com.getType() == Compute::Type::msd) {
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_harmonic_k_10) {
msd_harmFile << com.getTime() << ", " << com.getAgg().getMean() << ", " << com.getTarget() << 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);
}
}
}

2
src/vec_trafos.h
View File

@ -10,4 +10,4 @@ inline static Eigen::Matrix2d rotation_Matrix(Eigen::Vector2d e) {
Eigen::Matrix2d mat;
mat << e[0], -e[1], e[1], e[0];
return mat;
}
}

7
test/.clang-tidy Normal file
View File

@ -0,0 +1,7 @@
---
Checks: '*,-fuchsia-*,-google-*,-zircon-*,-abseil-*,-modernize-use-trailing-return-type,-llvm*,-readability-magic-numbers,-cppcoreguidelines-avoid-magic-numbers'
WarningsAsErrors: ''
HeaderFilterRegex: ''
FormatStyle: none

2
test/CMakeLists.txt
View File

@ -11,7 +11,7 @@ target_link_libraries(catch_main PRIVATE project_options)
set(TEST_FILES
test_LiveAgg.cpp
test_Rod2d.cpp
test_Calculation.cpp test_Compute.cpp test_Simulation.cpp test_Integratoren.cpp)
test_Compute.cpp test_Simulation.cpp test_Integratoren.cpp)
add_executable(tests ${TEST_FILES})
target_link_libraries(tests PRIVATE project_warnings project_options catch_main Eigen3::Eigen3 integratoren)

64
test/test_Calculation.cpp
View File

@ -1,64 +0,0 @@
//
// Created by jholder on 22.10.21.
//
#include "Calculation.h"
#include <catch2/catch.hpp>
#include "Integratoren2d_forceless.h"
#include "LiveAgg.hpp"
TEST_CASE("Basic run of Calculation") {
Calculation euler(Integratoren2d_forceless::Set1_Euler,
{{Compute::Type::msd, 10}, {Compute::Type::oaf, 10}}, 1,
1);
Calculation heun(Integratoren2d_forceless::Set2_Heun,
{{Compute::Type::msd, 10}, {Compute::Type::oaf, 10}}, 1,
1);
Calculation exact(Integratoren2d_forceless::Set3_Exact,
{{Compute::Type::msd, 10}, {Compute::Type::oaf, 10}}, 1,
1);
Calculation bdas(Integratoren2d_forceless::Set4_BDAS,
{{Compute::Type::msd, 10}, {Compute::Type::oaf, 10}}, 1,
1);
SECTION("Euler") {
euler.run(100);
CHECK(euler.getRod().getE().norm() == Catch::Detail::Approx(1.0));
CHECK(euler.getComputes()[0].getType() == Compute::Type::msd);
CHECK(euler.getComputes()[1].getType() == Compute::Type::oaf);
}
SECTION("Heun") {
heun.run(100);
CHECK(heun.getRod().getE().norm() == Catch::Detail::Approx(1.0));
}
SECTION("Exact") {
exact.run(100);
CHECK(exact.getRod().getE().norm() == Catch::Detail::Approx(1.0));
}
SECTION("Euler") {
bdas.run(100);
CHECK(bdas.getRod().getE().norm() == Catch::Detail::Approx(1.0));
}
SECTION("MSD") {
euler.run(100);
CHECK(euler.getComputes()[0].getAgg().getNumPoints() == 10);
CHECK(euler.getComputes()[1].getAgg().getNumPoints() == 10);
}
SECTION("Deterministic") {
Calculation determ(Integratoren2d_forceless::Set0_deterministic,
{{Compute::Type::msd, 1}, {Compute::Type::oaf, 1}},
1, 1);
determ.run(10);
auto time = 1;
auto targetMSD = pow(0.01, 2) * time;
auto targetOAF = cos(0.1);
CHECK(determ.getComputes()[0].getAgg().getMean() ==
Catch::Detail::Approx(targetMSD));
CHECK(determ.getComputes()[1].getAgg().getMean() ==
Catch::Detail::Approx(targetOAF));
}
}

1
test/test_Compute.cpp
View File

@ -5,7 +5,6 @@
#include "Compute.h"
#include <catch2/catch.hpp>
#include "Integratoren2d_forceless.h"
TEST_CASE("Compute") {
Rod2d rod(1.0);

359
test/test_Integratoren.cpp
View File

@ -1,357 +1,62 @@
//
// Created by jholder on 24.10.21.
//
#include "Calculation.h"
#include "Compute.h"
#include "Integratoren2d_force.h"
#include "Integratoren2d_forceless.h"
#include "Integratoren2d.h"
#include "force_lambdas.h"
#include <catch2/catch.hpp>
#include <utility>
TEST_CASE("Forceless Integratoren") {
TEST_CASE("Integrator Check") {
const size_t SEED = Catch::rngSeed();
const double length = GENERATE(1.0,1.4,2.0);
const double length = GENERATE(1.0, 1.4, 2.0);
constexpr int steps = 10000;
constexpr double delta = 0.1;
using type = std::tuple<std::function<void(Rod2d &, Simulation &)>, Calculation::inte_force_type, std::string>;
auto[integrator, integrator_force, name] = GENERATE(
type(Integratoren2d::Set1_Euler, Integratoren2d::Set1_Euler_f, "Euler"),
type(Integratoren2d::Set2_Heun, Integratoren2d::Set2_Heun_f, "Heun"),
type(Integratoren2d::Set3_Exact, Integratoren2d::Set3_Exact_f, "Exact"),
type(Integratoren2d::Set4_BDAS, Integratoren2d::Set4_BDAS_f, "BDAS"),
type(Integratoren2d::Set5_MBD, Integratoren2d::Set5_MBD_f, "MBD"));
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,
Calculation calc(integrator,
{{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") {
0.01, SEED, length);
calc.run(steps);
SECTION("ForceLess") {
CAPTURE(name);
CAPTURE(length);
size_t i = 0;
for (auto &c: euler_zero.getComputes()) {
for (const auto &c: calc.getComputes()) {
CAPTURE(c);
CHECK(c.getDifference() <= delta*c.getAgg().getMean());
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") {
Calculation calc_f(integrator_force,
{{Compute::Type::msd, 1},
{Compute::Type::oaf, 1},
{Compute::Type::empxx, 1},
{Compute::Type::empyy, 1}},
0.01, SEED, zero_Force, zero_Torque, length);
calc_f.run(steps);
SECTION("ForceLess") {
CAPTURE(length);
size_t i = 0;
for (auto &c: exact_zero.getComputes()) {
CAPTURE(name);
for (const auto &c: calc_f.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;
CHECK(c.getDifference() <= delta * c.getAgg().getMean());
}
}
}
/*
SECTION("Euler") {
{
}
Calculation euler_force(
Integratoren2d_force::Set1_Euler,
{{Compute::Type::msd, 1}, {Compute::Type::oaf, 1}}, 0.01, SEED,
force, torque);
{
euler_force.run(10000);
for (auto &c : euler_force.getComputes()) {
CHECK(c.getDifference() <= 0.005);
}
}
CHECK(euler.getComputes()[0].getAgg().getMean() ==
euler_force.getComputes()[0].getAgg().getMean());
CHECK(euler.getComputes()[1].getAgg().getMean() ==
euler_force.getComputes()[1].getAgg().getMean());
}
SECTION("Heun") {
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);
{
heun.run(10000);
for (auto &c : heun.getComputes()) {
CHECK(c.getDifference() <= 0.005);
}
}
Calculation heun_force(Integratoren2d_force::Set2_Heun,
{{Compute::Type::msd, 1},
{Compute::Type::oaf, 1},
{Compute::Type::empxx, 1},
{Compute::Type::empyy, 1}},
0.01, SEED, force, torque);
{
heun_force.run(10000);
for (auto &c : heun_force.getComputes()) {
CHECK(c.getDifference() <= 0.005);
}
}
CHECK(heun.getComputes()[0].getAgg().getMean() ==
heun_force.getComputes()[0].getAgg().getMean());
CHECK(heun.getComputes()[1].getAgg().getMean() ==
heun_force.getComputes()[1].getAgg().getMean());
}
SECTION("Exact") {
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);
{
exact.run(10000);
for (auto &c : exact.getComputes()) {
CHECK(c.getDifference() <= 0.005);
}
}
Calculation exact_force(Integratoren2d_force::Set3_Exact,
{{Compute::Type::msd, 1},
{Compute::Type::oaf, 1},
{Compute::Type::empxx, 1},
{Compute::Type::empyy, 1}},
0.01, SEED, force, torque);
{
exact_force.run(10000);
for (auto &c : exact_force.getComputes()) {
CHECK(c.getDifference() <= 0.005);
}
}
CHECK(exact.getComputes()[0].getAgg().getMean() ==
Approx(exact_force.getComputes()[0].getAgg().getMean())
.epsilon(1e-10));
CHECK(exact.getComputes()[1].getAgg().getMean() ==
exact_force.getComputes()[1].getAgg().getMean());
}
SECTION("BDAS") {
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);
{
bdas.run(10000);
for (auto &c : bdas.getComputes()) {
CHECK(c.getDifference() <= 0.005);
}
}
Calculation bdas_force(Integratoren2d_force::Set4_BDAS,
{{Compute::Type::msd, 1},
{Compute::Type::oaf, 1},
{Compute::Type::empxx, 1},
{Compute::Type::empyy, 1}},
0.01, SEED, force, torque);
{
bdas_force.run(10000);
for (auto &c : bdas_force.getComputes()) {
CHECK(c.getDifference() <= 0.005);
}
}
CHECK(bdas.getComputes()[0].getAgg().getMean() ==
Approx(bdas_force.getComputes()[0].getAgg().getMean())
.epsilon(10e-10));
CHECK(bdas.getComputes()[1].getAgg().getMean() ==
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));
}
}
*/