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

test

Browse Source
This commit is contained in:
Jacob Holder 2021-12-22 13:45:26 +01:00
parent 82e81eb91e
commit b5e8ef350a
Signed by: jacob
GPG Key ID: 2194FC747048A7FD
24 changed files with 413 additions and 863 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

165
Auswertung/notebooks/kovergence.ipynb Normal file
View File

File diff suppressed because one or more lines are too long

3
Auswertung/requirements.txt
View File

@ -1,4 +1,5 @@
numpy~=1.19.5
seaborn~=0.11.2
matplotlib~=3.3.4
pandas~=1.1.5
pandas~=1.1.5
jupyter~=6.4.6

17
C++/src/CMakeLists.txt
View File

@ -3,11 +3,11 @@ find_package(spdlog)
find_package(Eigen3)
add_library(logging logging.cpp)
# Generic test that uses conan libs
add_library(integratoren SHARED LiveAgg.cpp Rod2d.cpp Simulation.cpp
Simulation.h Integratoren2d.cpp
Calculation.cpp Compute.cpp)
Compute.cpp konvergenz_runner.cpp konvergenz_runner.h)
target_include_directories(integratoren PUBLIC .)
target_link_libraries(
integratoren
@ -15,7 +15,17 @@ target_link_libraries(
project_warnings
fmt::fmt
spdlog::spdlog
Eigen3::Eigen3)
Eigen3::Eigen3
logging
)
target_include_directories(logging PUBLIC .)
target_link_libraries(
logging
PRIVATE project_options
project_warnings
fmt::fmt
spdlog::spdlog)
add_executable(main main.cpp)
find_package(OpenMP)
@ -28,4 +38,5 @@ target_link_libraries(
project_warnings
integratoren
Eigen3::Eigen3
fmt::fmt
)

94
C++/src/Calculation.cpp
View File

@ -1,94 +0,0 @@
//
// Created by jholder on 22.10.21.
//
#include "Calculation.h"
#include <fstream>
#include <utility>
void Calculation::run(size_t steps) {
std::ofstream file ("dump.xyz");
for (size_t step = 0; step < steps; ++step) {
m_integrator(rod, sim);
for (auto &comp : computes) {
comp.eval(rod);
}
file<<rod.getPos()[0]<<" "<<rod.getPos()[0]<<" 0\n";
}
file.close();
}
const Rod2d &Calculation::getRod() const { return rod; }
const std::vector<Compute> &Calculation::getComputes() const {
return computes;
}
const Simulation &Calculation::getSim() const { return sim; }
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),
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, Compute::Force::zero_F, pair.second, sim));
}
}
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, Compute::Force type_force,
torque_type torque, double length)
: 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, type_force, pair.second, 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),
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, Compute::Force::zero_F, pair.second, sim));
}
}
Calculation::Calculation(
Calculation::inte_force_type t_integrator,
const std::vector<std::pair<Compute::Type, size_t>>& t_computes, double deltaT,
size_t seed, Calculation::force_type force, Compute::Force type_force,
Calculation::torque_type torque, double length)
: 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, type_force, pair.second, sim));
}
}
void Calculation::reset() {
rod = start_rod;
for (auto &c : computes) {
c.reset(start_rod);
}
}

57
C++/src/Calculation.h
View File

@ -1,57 +0,0 @@
//
// Created by jholder on 22.10.21.
//
#pragma once
#include <functional>
#include "Compute.h"
#include "Rod2d.hpp"
#include "Simulation.h"
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)>;
using torque_type = std::function<double(Eigen::Vector2d, Eigen::Vector2d)>;
using inte_force_type =
std::function<void(Rod2d &, Simulation &, force_type, torque_type)>;
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, Compute::Force force_type,
torque_type torque, double length = 1.0);
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 = 1.0);
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 = 1.0);
Calculation(inte_force_type t_integrator,
const std::vector<std::pair<Compute::Type, size_t>>& t_computes,
double deltaT, size_t seed, force_type force,
Compute::Force force_type, torque_type torque,
double length = 1.0);
[[nodiscard]] const std::vector<Compute> &getComputes() const;
void run(size_t steps);
[[nodiscard]] const Rod2d &getRod() const;
};

45
C++/src/Compute.cpp
View File

@ -2,26 +2,34 @@
// Created by jholder on 22.10.21.
//
//#define SPDLOG_ACTIVE_LEVEL SPDLOG_LEVEL_TRACE
#include "Compute.h"
#include <iostream>
#include <utility>
#include <spdlog/spdlog.h>
#include "Rod2d.hpp"
#include "Simulation.h"
void Compute::evalMSD(const Rod2d &rod2D) {
SPDLOG_TRACE("Evaluating MSD");
const auto &new_pos = rod2D.getPos();
auto old_pos = start_rod.getPos();
auto diff = new_pos - old_pos;
auto msd = diff.dot(diff);
SPDLOG_TRACE("Feeding {}", msd);
agg.feed(msd);
}
void Compute::evalX(const Rod2d &rod2D) {
SPDLOG_TRACE("Evaluating X");
const auto &new_pos = rod2D.getPos();
auto msd = (new_pos[0]);
agg.feed(msd);
SPDLOG_TRACE("Feeding {}", msd);
}
void Compute::evalX_squared(const Rod2d &rod2D) {
@ -56,7 +64,8 @@ void Compute::eval_empYY(const Rod2d &rod2D) {
void Compute::eval(const Rod2d &rod2D) {
time_step++;
if (time_step % every == 0) {
SPDLOG_TRACE("resetting: {} time_step {}", resetting, time_step);
if (time_step % every == 0){
if (resetting or time_step == every) {
switch (type) {
case Type::msd:
@ -78,10 +87,7 @@ void Compute::eval(const Rod2d &rod2D) {
evalX_squared(rod2D);
break;
}
if (resetting) {
start_rod = rod2D;
}
return;
start_rod = rod2D;
}
}
}
@ -94,7 +100,7 @@ Compute::Compute(Rod2d rod, Type t_type, Force force, size_t t_every,
time(sim.getMDeltaT() * static_cast<double>(every)) {
switch (type) {
case Type::msd: {
resetting = false;
resetting = true;
type_str = "msd";
switch (force) {
case Force::zero_F:
@ -102,22 +108,24 @@ Compute::Compute(Rod2d rod, Type t_type, Force force, size_t t_every,
break;
case Force::const_F:
target = 0.0;
break;
break; //TODO missing analytic
case Force::harmonic_F:
target = 2 * (1 - 2 * std::exp(- time) + std::exp(-2 * time)) + 2 * (1 - std::exp(-2 * time));
target =
2 * (1 - 2 * std::exp(-time) + std::exp(-2 * time)) +
2 * (1 - std::exp(-2 * time));
break;
}
break;
}
case Type::oaf: {
resetting = false;
resetting = true;
type_str = "oaf";
target = std::exp(-rod.getDRot() * time);
break;
}
case Type::empxx: {
resetting = false;
resetting = true;
type_str = "empxx";
const double Dmean = 0.5 * (rod.getDiff().trace());
const double u = 4.0;
@ -130,12 +138,12 @@ Compute::Compute(Rod2d rod, Type t_type, Force force, size_t t_every,
break;
case Force::const_F:
case Force::harmonic_F:
target = 0.0;
target = 0.0; //TODO missing analytic
break;
}
} break;
case Type::empyy: {
resetting = false;
resetting = true;
type_str = "empyy";
const double Dmean = 0.5 * (rod.getDiff().trace());
const double u = 4.0;
@ -148,7 +156,7 @@ Compute::Compute(Rod2d rod, Type t_type, Force force, size_t t_every,
break;
case Force::const_F:
case Force::harmonic_F:
target = 0.0;
target = 0.0; //TODO missing analytic
break;
}
} break;
@ -160,7 +168,7 @@ Compute::Compute(Rod2d rod, Type t_type, Force force, size_t t_every,
target = rod.getPos().sum();
break;
case Force::const_F:
target = 0.0;
target = 0.0; //TODO missing analytic
break;
case Force::harmonic_F:
target = rod.getPos()[0] * exp(-1 * time);
@ -174,10 +182,8 @@ Compute::Compute(Rod2d rod, Type t_type, Force force, size_t t_every,
switch (force) {
case Force::zero_F:
target = 0.0;
break;
case Force::const_F:
target = 0.0;
target = 0.0; //TODO missing analytic
break;
case Force::harmonic_F:
target = pow(rod.getPos()[0], 2) * exp(-2 * time) +
@ -185,10 +191,7 @@ Compute::Compute(Rod2d rod, Type t_type, Force force, size_t t_every,
break;
}
} break;
}
resetting = false;
}
const LiveAgg &Compute::getAgg() const { return agg; }
@ -207,6 +210,8 @@ std::ostream &operator<<(std::ostream &os, const Compute &compute) {
}
void Compute::reset(const Rod2d &rod) {
SPDLOG_TRACE("Resettign Compute for new run");
time_step = 0;
start_rod = rod;
}

47
C++/src/Integratoren2d.h
View File

@ -9,36 +9,49 @@
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);
public:
using integrator_t = std::function<void(
Rod2d &, Simulation &,
const std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)>
&,
const std::function<double(Eigen::Vector2d, Eigen::Vector2d)> &)>;
using force_t = const std::function<Vector(Vector, Vector)>;
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_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_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_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_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);
};

6
C++/src/LiveAgg.cpp
View File

@ -12,6 +12,9 @@ double LiveAgg::getSEM() const noexcept {
double LiveAgg::getMean() const noexcept { return mean; }
void LiveAgg::feed(double value) noexcept {
if (value != value) {
return;
}
num_of_data += 1;
auto delta = value - mean;
mean += delta / num_of_data;
@ -21,6 +24,7 @@ void LiveAgg::feed(double value) noexcept {
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;
os << "num_of_data: " << agg.num_of_data << " mean: " << agg.mean
<< " S: " << agg.S;
return os;
}

2
C++/src/force_lambdas.h
View File

@ -1,7 +1,7 @@
//
// Created by jholder on 27.10.21.
//
#include "Eigen/Dense"
#pragma once
[[maybe_unused]] const static auto harmonic_Force =
[](const Eigen::Vector2d &pos,

51
C++/src/konvergenz_runner.cpp Normal file
View File

@ -0,0 +1,51 @@
//
// Created by jholder on 16.12.21.
//
#include "konvergenz_runner.h"
#include <fstream>
#include <iostream>
#include "force_lambdas.h"
void Konvergenz::konvergenz(Compute::Force force,
const Integratoren2d::integrator_t& integrator,
size_t every, std::ofstream& file) {
Integratoren2d::force_t force_lambda = harmonic_Force;
switch (force) {
case Compute::Force::zero_F:
break;
case Compute::Force::const_F:
break;
case Compute::Force::harmonic_F:
break;
}
double delta_t = 1.0 / static_cast<double>(every);
Rod2d rod(1.0);
Simulation sim(delta_t, 1234);
Compute com(rod, Compute::Type::oaf, force, every, sim);
while (true) {
for (int i = 0; i < 100000; ++i) {
integrator(rod, sim, force_lambda, zero_Torque);
com.eval(rod);
}
std::cout << every << " " << com.getDifference() << " "
<< com.getAgg().getSEM() << " \r";
if (com.getDifference() >= com.getAgg().getSEM() * 10 and
com.getAgg().getSEM() < 0.01) {
std::cout << "\n";
file << every << " " << com.getDifference() << " "
<< com.getAgg().getSEM() << '\n';
file.flush();
break;
}
}
}
void Konvergenz::runner(Compute::Force force,
const Integratoren2d::integrator_t& integrator,
std::string filename) {
std::ofstream file;
file.open(filename);
for (size_t i = 1; i < 100; i *= 2) {
konvergenz(force, integrator, i, file);
}
}

18
C++/src/konvergenz_runner.h Normal file
View File

@ -0,0 +1,18 @@
//
// Created by jholder on 16.12.21.
//
#ifndef MYPROJECT_KONVERGENZ_RUNNER_H
#define MYPROJECT_KONVERGENZ_RUNNER_H
#include "Compute.h"
#include "Integratoren2d.h"
class Konvergenz {
public:
static void runner(Compute::Force force,
const Integratoren2d::integrator_t& integrator,
std::string filename);
static void konvergenz(Compute::Force force,
const Integratoren2d::integrator_t& integrator,
size_t every, std::ofstream& file);
};
#endif // MYPROJECT_KONVERGENZ_RUNNER_H

184
C++/src/legacy/Integratoren2d_force.cpp
View File

@ -1,184 +0,0 @@
//
// Created by jholder on 21.10.21.
//
#include "Integratoren2d_force.h"
#include "vec_trafos.h"
using Vector = Eigen::Vector2d;
using Matrix = Eigen::Matrix2d;
const Vector Integratoren2d_force::unitVec = {1, 0};
void Integratoren2d_force::Set1_Euler(
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)
auto e = rod2D.getE();
// translation
Vector rand = {sim.getNorm(std), sim.getNorm(std)}; // gaussian noise
Vector trans_part = rod2D.getDiff_Sqrt() * rand;
Vector trans_lab = rotation_Matrix(e) * trans_part;
// Force
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();
// rotation
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();
// apply
rod2D.setPos(rod2D.getPos() + trans_lab + F_trans);
rod2D.setE(new_e);
}
Vector Integratoren2d_force::Heun_predictor_pos(
const Rod2d &rod2D, Simulation &sim,
const std::function<Vector(Vector, Vector)> &force) {
auto standard_dev = sim.getSTD();
Vector rand_pred = {sim.getNorm(standard_dev),
sim.getNorm(standard_dev)}; // gaussian noise
Vector trans_pred_part =
rod2D.getDiff_Sqrt() *
rand_pred; // translation vector in Particle System
Vector trans_pred_lab = rotation_Matrix(rod2D.getE()) * trans_pred_part;
Vector F_pred_lab = force(rod2D.getPos(), rod2D.getE());
Vector F_pred_part = rotation_Matrix(rod2D.getE()).inverse() * F_pred_lab;
Vector F_pred_trans = rotation_Matrix(rod2D.getE())*rod2D.getDiff_Sqrt() * F_pred_part;
F_pred_trans *= sim.getMDeltaT() ;
return F_pred_trans + trans_pred_lab;
}
Vector Integratoren2d_force::Heun_predictor_rot(
const Rod2d &rod2D, Simulation &sim,
const std::function<double(Vector, Vector)> &torque) {
auto std = sim.getSTD();
double rot_predict =
sim.getNorm(std) * rod2D.getDRot_Sqrt(); // rotationsmatrix verwenden.
const Vector& e = rod2D.getE();
double M_predict_rot = torque(rod2D.getPos(), rod2D.getE()) *
rod2D.getDRot() * sim.getMDeltaT();
Vector e_change_predict = (rot_predict + M_predict_rot) * ortho(e);
return e_change_predict;
}
void Integratoren2d_force::Set2_Heun(
Rod2d &rod2D, Simulation &sim,
const std::function<Vector(Vector, Vector)> &force,
const std::function<double(Vector, Vector)> &torque) {
Vector delta_pos_predictor = Heun_predictor_pos(rod2D, sim, force);
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();
Rod2d pred_rod = rod2D;
pred_rod.setPos(pos_predictor);
pred_rod.setE(e_predictor);
Vector delta_e_future = Heun_predictor_rot(pred_rod, sim, torque);
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_force::Set3_Exact(
Rod2d &rod2D, Simulation &sim,
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
Vector trans_part =
rod2D.getDiff_Sqrt() * rand; // translation vector in Particle System
Vector F_lab = force(rod2D.getPos(), rod2D.getE());
Vector F_part = rotation_Matrix(rod2D.getE()).inverse() * F_lab;
Vector F_trans = rotation_Matrix(rod2D.getE())*rod2D.getDiff_Sqrt() * F_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() * sim.getMDeltaT();
auto correction =
-0.5 * pow(rod2D.getDRot(), 2) * sim.getMDeltaT();
Vector new_e = e + (rot + M_rot) * ortho(e) + correction * e;
new_e.normalize();
// apply
rod2D.setPos(rod2D.getPos() + trans_lab + F_trans);
rod2D.setE(new_e);
}
void Integratoren2d_force::Set4_BDAS(
Rod2d &rod2D, Simulation &sim,
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
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 = rotation_Matrix(rod2D.getE()).inverse() * F_lab;
Vector F_trans = rotation_Matrix(rod2D.getE())*rod2D.getDiff_Sqrt() * F_part;
F_trans *= sim.getMDeltaT();
// rotation
double rot =
sim.getNorm(std) * rod2D.getDRot_Sqrt(); // rotationsmatrix verwenden.
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());
}

52
C++/src/legacy/Integratoren2d_force.h
View File

@ -1,52 +0,0 @@
//
// Created by jholder on 21.10.21.
//
#pragma once
#include "Rod2d.hpp"
#include "Simulation.h"
class Integratoren2d_force {
public:
static void Set1_Euler(
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 Set2_Heun(
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,
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,
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,
const std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)>& force,
const std::function<double(Eigen::Vector2d, Eigen::Vector2d)>& torque);
private:
static Eigen::Vector2d Heun_predictor_pos(
const Rod2d &rod2D, Simulation &sim,
const std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)>
&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);
};

123
C++/src/legacy/Integratoren2d_forceless.cpp
View File

@ -1,123 +0,0 @@
//
// Created by jholder on 21.10.21.
//
#include "Integratoren2d_forceless.h"
#include "vec_trafos.h"
void Integratoren2d_forceless::Set0_deterministic(Rod2d &rod2D,
Simulation & /*sim*/) {
auto trans_lab = Eigen::Vector2d({0, 0.01});
rod2D.setPos(rod2D.getPos() + trans_lab);
Eigen::Rotation2D<double> rot(0.1);
rod2D.setE(rot * rod2D.getE());
}
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)};
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 new_e = e + rot * ortho(e);
// apply
rod2D.setPos(rod2D.getPos() + trans_lab);
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;
Eigen::Vector2d trans_lab = rotation_Matrix(e) * trans_part;
Eigen::Vector2d pos_predictor = rod2D.getPos() + trans_lab;
// rotation
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));
// apply
rod2D.setPos(pos_predictor); // pos with Euler
rod2D.setE(e_integrated.normalized());
}
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)};
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();
auto correction =
-0.5 * pow(rod2D.getDRot(), 2) * sim.getMDeltaT();
Eigen::Vector2d delta_e = correction * e + rot * ortho(e);
// apply
rod2D.setPos(rod2D.getPos() + trans_lab);
rod2D.setE((e + delta_e).normalized());
}
// 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::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() * e).normalized();
rod2D.setPos(rod2D.getPos() + trans_lab);
rod2D.setE(e_new);
}
void Integratoren2d_forceless::Set5_MBD(Rod2d & rod2D,
Simulation & sim) {
auto std = sim.getSTD();
Eigen::Vector2d e = rod2D.getE();
//rotation
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();
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_integrated);
rod2D.setE(e_integrated.normalized());
}

25
C++/src/legacy/Integratoren2d_forceless.h
View File

@ -1,25 +0,0 @@
//
// Created by jholder on 21.10.21.
//
#pragma once
#include "Rod2d.hpp"
#include "Simulation.h"
class Integratoren2d_forceless {
public:
static void Set1_Euler(Rod2d &rod2D, Simulation &sim);
static void Set2_Heun(Rod2d &rod2D, Simulation &sim);
static void Set3_Exact(Rod2d &rod2D, Simulation &sim);
static void Set4_BDAS(Rod2d &rod2D, Simulation &sim);
static const Eigen::Vector2d unitVec;
static void Set5_MBD(Rod2d &rod2D, Simulation &sim);
static void Set0_deterministic(Rod2d &rod2D, Simulation &sim);
};

20
C++/src/logging.cpp Normal file
View File

@ -0,0 +1,20 @@
//
// Created by jholder on 16.12.21.
//
#include "logging.hpp"
#include <spdlog/spdlog.h>
void LOGGER::setLogging(LOG_LEVEL level) {
switch (level) {
case LOG_LEVEL::WARN:
spdlog::set_level(spdlog::level::warn);
break;
case LOG_LEVEL::DEBUG:
spdlog::set_level(spdlog::level::debug);
break;
case LOG_LEVEL::TRACE:
spdlog::set_level(spdlog::level::trace);
break;
}
}

15
C++/src/logging.hpp Normal file
View File

@ -0,0 +1,15 @@
//
// Created by jholder on 16.12.21.
//
#ifndef MYPROJECT_LOGGING_H
#define MYPROJECT_LOGGING_H
enum LOG_LEVEL{WARN,DEBUG, TRACE};
class LOGGER{
public:
static void setLogging(LOG_LEVEL level);
};
#endif // MYPROJECT_LOGGING_H

210
C++/src/main.cpp
View File

@ -1,189 +1,29 @@
//
// Created by jholder on 22.10.21.
//
#include <fstream>
#include <iostream>
#include <filesystem>
#include <utility>
#include "Calculation.h"
#include "Integratoren2d.h"
#include "force_lambdas.h"
#include <fmt/core.h>
constexpr size_t SEED = 1364;
constexpr double stepSize = 0.01;
constexpr size_t n_computes = 10;
constexpr size_t num_integratoren = 5;
constexpr size_t num_forces = 3;
constexpr size_t num_length = 1;
constexpr size_t delta_compute = 10;
constexpr size_t numStep = 50000;
inline const std::string header = {"time,val,target,std,numPoints\n"};
void run(size_t integrator_index, size_t force_index, size_t length_index) {
using type = std::pair<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];
auto force_type = std::vector<Compute::Force>(
{Compute::Force::zero_F, Compute::Force::const_F,
Compute::Force::harmonic_F})[force_index];
namespace fs = std::filesystem;
std::string folder =
"out/" + force_name + name + "_L" + std::to_string(length_index) + "_";
fs::create_directories("out");
// double length = std::vector<double>({1.0, 1.5, 2.0})[length_index];
double length = 1.0;
{
std::vector<std::pair<Compute::Type, size_t>> computes;
for (size_t i = 1; i < n_computes; ++i) {
computes.emplace_back(Compute::Type::msd, i * delta_compute);
computes.emplace_back(Compute::Type::oaf, i * delta_compute);
computes.emplace_back(Compute::Type::empxx, i * delta_compute);
computes.emplace_back(Compute::Type::empyy, i * delta_compute);
computes.emplace_back(Compute::Type::x, i * delta_compute);
computes.emplace_back(Compute::Type::x_squared, i * delta_compute);
}
Calculation euler(integrator, computes, stepSize, SEED, force,
force_type, zero_Torque, length);
for (size_t i = 0; i < numStep; ++i) {
euler.run(n_computes * delta_compute);
euler.reset();
}
std::ofstream msdFile(folder + "msd.dat");
std::ofstream oafFile(folder + "oaf.dat");
std::ofstream empxxFile(folder + "empxx.dat");
std::ofstream empyyFile(folder + "empyy.dat");
std::ofstream xFile(folder + "x.dat");
std::ofstream xsqFile(folder + "x_squared.dat");
xFile << header;
xsqFile << header;
msdFile << header;
oafFile << header;
empxxFile << header;
empyyFile << header;
for (const auto &com : euler.getComputes()) {
if (com.getType() == Compute::Type::msd) {
msdFile << com.getTime() << ", " << com.getAgg().getMean()
<< ", " << com.getTarget() << ", "
<< com.getAgg().getSEM() << ", "
<< com.getAgg().getNumPoints() << std::endl;
}
if (com.getType() == Compute::Type::oaf) {
oafFile << com.getTime() << ", " << com.getAgg().getMean()
<< ", " << com.getTarget() << ", "
<< com.getAgg().getSEM() << ", "
<< com.getAgg().getNumPoints() << std::endl;
}
if (com.getType() == Compute::Type::empxx) {
empxxFile << com.getTime() << ", " << com.getAgg().getMean()
<< ", " << com.getTarget() << ", "
<< com.getAgg().getSEM() << ", "
<< com.getAgg().getNumPoints() << std::endl;
}
if (com.getType() == Compute::Type::empyy) {
empyyFile << com.getTime() << ", " << com.getAgg().getMean()
<< ", " << com.getTarget() << ", "
<< com.getAgg().getSEM() << ", "
<< com.getAgg().getNumPoints() << std::endl;
}
if (com.getType() == Compute::Type::msd) {
msdFile << com.getTime() << ", " << com.getAgg().getMean()
<< ", " << com.getTarget() << ", "
<< com.getAgg().getSEM() << ", "
<< com.getAgg().getNumPoints() << std::endl;
}
if (com.getType() == Compute::Type::x) {
xFile << com.getTime() << ", " << com.getAgg().getMean() << ", "
<< com.getTarget() << ", " << com.getAgg().getSEM()
<< ", " << com.getAgg().getNumPoints() << std::endl;
}
if (com.getType() == Compute::Type::x_squared) {
xsqFile << com.getTime() << ", " << com.getAgg().getMean()
<< ", " << com.getTarget() << ", "
<< com.getAgg().getSEM() << ", "
<< com.getAgg().getNumPoints() << std::endl;
}
}
}
std::cout << "Finished run " << folder << std::endl;
#include "konvergenz_runner.h"
void msd(){
Konvergenz::runner(Compute::Force::harmonic_F, Integratoren2d::Set1_Euler_f,
"1Euler_msd_harmonic.dat");
Konvergenz::runner(Compute::Force::harmonic_F, Integratoren2d::Set2_Heun_f,
"2Heun_msd_harmonic.dat");
Konvergenz::runner(Compute::Force::harmonic_F, Integratoren2d::Set3_Exact_f,
"3Exact_msd_harmonic.dat");
Konvergenz::runner(Compute::Force::harmonic_F, Integratoren2d::Set4_BDAS_f,
"4BDAS_msd_harmonic.dat");
Konvergenz::runner(Compute::Force::harmonic_F, Integratoren2d::Set5_MBD_f,
"5MBD_msd_harmonic.dat");
}
void konv(int delta, std::ofstream &file,
Calculation::inte_force_type integrator, size_t seed) {
std::vector<std::pair<Compute::Type, size_t>> computes;
double step_Size = 1.0 / delta;
computes.emplace_back(Compute::Type::msd, delta);
Calculation euler(integrator, computes, step_Size, seed, harmonic_Force,
Compute::Force::harmonic_F, zero_Torque, 1.0);
for (int j = 0; j < 1000000; ++j) {
for (int i = 0; i < 1000; ++i) {
euler.run(static_cast<size_t>(delta));
euler.reset();
}
std::cout << euler.getComputes()[0].getAgg().getSEM() << " "
<< euler.getComputes()[0].getDifference() << std::endl;
if (euler.getComputes()[0].getAgg().getSEM() * 10 <
euler.getComputes()[0].getDifference()) {
file << delta << " " << euler.getComputes()[0].getAgg().getMean()
<< " " << euler.getComputes()[0].getAgg().getSEM() << " "
<< euler.getComputes()[0].getTarget() << " "
<< euler.getComputes()[0].getDifference() << std::endl;
return;
}
}
file << "timeout: " << delta << " "
<< euler.getComputes()[0].getAgg().getMean() << " "
<< euler.getComputes()[0].getAgg().getSEM() << " "
<< euler.getComputes()[0].getTarget() << " "
<< euler.getComputes()[0].getDifference() << std::endl;
void oaf(){
Konvergenz::runner(Compute::Force::harmonic_F, Integratoren2d::Set1_Euler_f,
"1Euler_oaf_harmonic.dat");
Konvergenz::runner(Compute::Force::harmonic_F, Integratoren2d::Set2_Heun_f,
"2Heun_oaf_harmonic.dat");
// Konvergenz::runner(Compute::Force::harmonic_F, Integratoren2d::Set3_Exact_f,
// "3Exact_oaf_harmonic.dat");
//Konvergenz::runner(Compute::Force::harmonic_F, Integratoren2d::Set4_BDAS_f,
// "4BDAS_oaf_harmonic.dat");
Konvergenz::runner(Compute::Force::harmonic_F, Integratoren2d::Set5_MBD_f,
"5MBD_oaf_harmonic.dat");
}
Calculation::inte_force_type integrator(int index) {
if (index % 5 == 0) return Integratoren2d::Set1_Euler_f;
if (index % 5 == 1) return Integratoren2d::Set2_Heun_f;
if (index % 5 == 2) return Integratoren2d::Set3_Exact_f;
if (index % 5 == 3) return Integratoren2d::Set4_BDAS_f;
if (index % 5 == 4) return Integratoren2d::Set5_MBD_f;
return Integratoren2d::Set1_Euler_f;
}
void konv_runner(int index) {
std::ofstream file;
file.open(fmt::format("out/{}_konv.out", index));
for (int i = 1; i < 256; i *= 2) {
konv(i, file, integrator(index), static_cast<size_t>(index / 5));
}
file.close();
}
int main(int argc, char *argv[]) {
int index = 0;
if (argc >1) index = std::stoi(argv[1]);
konv_runner(index);
int main(){
oaf();
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),
(j % (num_integratoren * num_forces)) / num_integratoren,
j / (num_integratoren * num_forces));
}
}

4
C++/test/CMakeLists.txt
View File

@ -5,13 +5,13 @@ include(CTest)
include(Catch)
add_library(catch_main STATIC catch_main.cpp)
target_link_libraries(catch_main PUBLIC Catch2::Catch2)
target_link_libraries(catch_main PUBLIC Catch2::Catch2 logging)
target_link_libraries(catch_main PRIVATE project_options)
set(TEST_FILES
test_LiveAgg.cpp
test_Rod2d.cpp
test_Compute.cpp test_Simulation.cpp test_Integratoren.cpp)
test_Compute.cpp)
add_executable(tests ${TEST_FILES})
target_link_libraries(tests PRIVATE project_warnings project_options catch_main Eigen3::Eigen3 integratoren)

10
C++/test/catch_main.cpp
View File

@ -1,3 +1,11 @@
#define CATCH_CONFIG_MAIN // This tells the catch header to generate a main
#define CATCH_CONFIG_RUNNER
#include <catch2/catch.hpp>
#include "logging.hpp"
int main( int argc, char* argv[] ) {
LOGGER::setLogging(LOG_LEVEL::TRACE);
auto result = Catch::Session().run( argc, argv );
return result;
}

36
C++/test/test_Compute.cpp
View File

@ -2,19 +2,41 @@
// Created by jholder on 22.10.21.
//
#include "Compute.h"
#include <catch2/catch.hpp>
#include "Compute.h"
TEST_CASE("Compute") {
Rod2d rod(1.0);
Rod2d rod_old(1.0);
Simulation sim(0.1, 1);
auto com = Compute(rod, Compute::Type::msd, 10, sim);
SECTION("Mean of same values") {
for (int i = 0; i < 100; ++i) {
SECTION("MSD") {
Compute com(rod, Compute::Type::msd, Compute::Force::zero_F, 2, sim);
for (int i = 0; i < 11; ++i) {
rod.setPos(rod.getPos() + Eigen::Vector2d{1, 1});
com.eval(rod);
}
CHECK(com.getAgg().getNumPoints() == 10);
CHECK(com.getAgg().getMean() == 0.0);
CHECK(com.getAgg().getNumPoints() == 5);
CHECK(com.getAgg().getMean() == 8);
}
SECTION("X") {
Compute com(rod, Compute::Type::x, Compute::Force::zero_F, 2, sim);
for (int i = 0; i < 11; ++i) {
rod.setPos(rod.getPos() + Eigen::Vector2d{1, 1});
com.eval(rod);
}
CHECK(com.getAgg().getNumPoints() == 1);
CHECK(com.getAgg().getMean() == 3);
for (int i = 0; i < 11; ++i) {
rod.setPos(rod.getPos() + Eigen::Vector2d{1, 1});
com.eval(rod);
if ((i + 1) % 2 == 0) {
rod = rod_old;
com.reset(rod_old);
}
}
CHECK(com.getAgg().getNumPoints() == 5);
CHECK(com.getAgg().getMean() == 3);
}
}

62
C++/test/test_Integratoren.cpp
View File

@ -1,62 +0,0 @@
//
// Created by jholder on 24.10.21.
//
#include "Calculation.h"
#include "Compute.h"
#include "Integratoren2d.h"
#include "force_lambdas.h"
#include <catch2/catch.hpp>
#include <utility>
TEST_CASE("Integrator Check") {
const size_t SEED = Catch::rngSeed();
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 calc(integrator,
{{Compute::Type::msd, 1},
{Compute::Type::oaf, 1},
{Compute::Type::empxx, 1},
{Compute::Type::empyy, 1}},
0.01, SEED, length);
calc.run(steps);
SECTION("ForceLess") {
CAPTURE(name);
CAPTURE(length);
size_t i = 0;
for (const auto &c: calc.getComputes()) {
CAPTURE(c);
CHECK(c.getDifference() <= delta * c.getAgg().getMean());
++i;
}
}
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);
CAPTURE(name);
for (const auto &c: calc_f.getComputes()) {
CAPTURE(c);
CHECK(c.getDifference() <= delta * c.getAgg().getMean());
}
}
}

2
C++/test/test_Rod2d.cpp
View File

@ -7,10 +7,12 @@
#include "Rod2d.hpp"
#include "vec_trafos.h"
TEST_CASE("Rods") {
Rod2d sphere(1);
Rod2d rod(2);
SECTION("Checking translation") {
rod.setPos({0,0});
rod.setPos(rod.getPos() + Eigen::Vector2d(1, -1));
auto newPos = rod.getPos();
REQUIRE(newPos[0] == 1);

28
C++/test/test_Simulation.cpp
View File

@ -1,28 +0,0 @@
//
// Created by jholder on 22.10.21.
//
#include <catch2/catch.hpp>
#include "Simulation.h"
TEST_CASE("Simulation") {
auto sim = Simulation(1, 0);
REQUIRE(sim.getSTD() == sqrt(1 * 2.0));
SECTION("Symmetrie") {
constexpr int num = 10000000;
double sum = 0.0;
for (int i = 0; i < num; ++i) {
sum += sim.getNorm(1);
}
CHECK(sum / num == Catch::Detail::Approx(0.0).margin(0.001));
}
SECTION("STD") {
constexpr int num = 10000000;
double sum = 0.0;
for (int i = 0; i < num; ++i) {
sum += pow(sim.getNorm(1), 2);
}
CHECK(sum / num == Catch::Detail::Approx(1.0).margin(0.001));
}
}