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

Minor fixes should run

Browse Source
This commit is contained in:
Jacob Holder 2021-10-28 13:43:59 +02:00
parent 7f947ff1d1
commit 2c16c36fdc
Signed by: jacob
GPG Key ID: 2194FC747048A7FD
7 changed files with 228 additions and 134 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -2,6 +2,8 @@ find_package(fmt)
find_package(spdlog)
find_package(Eigen3)
# Generic test that uses conan libs
add_library(integratoren SHARED LiveAgg.cpp Rod2d.cpp Simulation.cpp
Simulation.h Integratoren2d.cpp
@ -16,6 +18,10 @@ target_link_libraries(
Eigen3::Eigen3)
add_executable(main main.cpp)
find_package(OpenMP)
if(OpenMP_CXX_FOUND)
target_link_libraries(main PUBLIC OpenMP::OpenMP_CXX)
endif()
target_link_libraries(
main
PRIVATE project_options

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

@ -9,7 +9,7 @@
void Calculation::run(size_t steps) {
for (size_t step = 0; step < steps; ++step) {
m_integrator(rod, sim);
for (auto &comp: computes) {
for (auto &comp : computes) {
comp.eval(rod);
}
}
@ -17,56 +17,74 @@ void Calculation::run(size_t steps) {
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, 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));
}
}
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), 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));
}
}
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, 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, pair.second, 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(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), start_rod(rod), sim(deltaT, seed) {
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, pair.second, sim));
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){
rod = start_rod;
for (auto &c : computes) {
c.reset(start_rod);
}
}

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

@ -10,35 +10,44 @@
#include "Simulation.h"
class Calculation {
private:
private:
Rod2d rod;
const Rod2d start_rod;
Simulation sim;
std::function<void(Rod2d &, Simulation &)> m_integrator;
std::vector<Compute> computes;
public:
public:
void reset();
[[nodiscard]] const Simulation &getSim() const;
using force_type = std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)>;
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)>;
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, torque_type torque, double length = 1.0);
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,
std::initializer_list<std::pair<Compute::Type, size_t>> t_computes,
const std::vector<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, std::vector<std::pair<Compute::Type,size_t>> t_computes,
double deltaT, size_t seed, force_type force, torque_type torque, 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;

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

@ -19,8 +19,13 @@ void Compute::evalMSD(const Rod2d &rod2D) {
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();
auto msd = (new_pos).sum();
agg.feed(msd);
}
void Compute::evalX_squared(const Rod2d &rod2D) {
const auto &new_pos = rod2D.getPos();
auto msd = (new_pos).squaredNorm();
agg.feed(msd);
}
@ -35,7 +40,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) / 2.0 / time;
double empxx = pow((new_pos - old_pos).dot(old_e), 2) / 2 / time;
agg.feed(empxx);
}
@ -44,7 +49,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) / 2.0 / time;
double empxx = pow((new_pos - old_pos).dot(old_e_ortho), 2) / 2 / time;
agg.feed(empxx);
}
@ -64,11 +69,6 @@ void Compute::eval(const Rod2d &rod2D) {
case Type::empyy:
eval_empYY(rod2D);
break;
case Type::msd_harmonic_k_1:
evalMSD(rod2D);
break;
case Type::msd_const:
break;
case Type::x:
if (time_step == every) {
evalX(rod2D);
@ -76,7 +76,7 @@ void Compute::eval(const Rod2d &rod2D) {
break;
case Type::x_squared:
if (time_step == every) {
evalMSD(rod2D);
evalX_squared(rod2D);
}
break;
}
@ -87,7 +87,8 @@ void Compute::eval(const Rod2d &rod2D) {
}
}
Compute::Compute(Rod2d rod, Type t_type, size_t t_every, Simulation &sim)
Compute::Compute(Rod2d rod, Type t_type, Force force, size_t t_every,
Simulation &sim)
: start_rod(std::move(rod)),
every(t_every),
time_step(0),
@ -97,9 +98,20 @@ Compute::Compute(Rod2d rod, Type t_type, size_t t_every, Simulation &sim)
case Type::msd: {
resetting = true;
type_str = "msd";
target = 2 * (rod.getDiff().trace()) * time;
switch (force) {
case Force::zero_F:
target = 2 * (rod.getDiff().trace()) * time;
break;
case Force::const_F:
target = static_cast<double>(NAN);
break;
case Force::harmonic_F:
target = 2 * (1.0 - std::exp(-2 * time));
break;
}
break;
}
case Type::oaf: {
resetting = true;
type_str = "oaf";
@ -112,8 +124,17 @@ Compute::Compute(Rod2d rod, Type t_type, size_t t_every, Simulation &sim)
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 * (1 - exp(-u * rod.getDRot() * time)) /
u / rod.getDRot() / time;
switch (force) {
case Force::zero_F:
target = Dmean - deltaD / 2 *
(1 - exp(-u * rod.getDRot() * time)) /
u / rod.getDRot() / time;
break;
case Force::const_F:
case Force::harmonic_F:
target = static_cast<double>(NAN);
break;
}
} break;
case Type::empyy: {
resetting = true;
@ -121,29 +142,51 @@ Compute::Compute(Rod2d rod, Type t_type, size_t t_every, Simulation &sim)
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 * (1 - exp(-u * rod.getDRot() * time)) /
u / rod.getDRot() / time;
} break;
case Type::msd_harmonic_k_1: {
resetting = true;
type_str = "msd_harmonic";
target = 2 * (1.0 - std::exp(-2 * time));
} break;
case Type::msd_const: {
resetting = true;
type_str = "msd_const";
target = -1.0;
switch (force) {
case Force::zero_F:
target = Dmean + deltaD / 2 * (1 - exp(-u * rod.getDRot() * time)) /
u / rod.getDRot() / time;
break;
case Force::const_F:
case Force::harmonic_F:
target = static_cast<double>(NAN);
break;
}
} break;
case Type::x: {
resetting = false;
type_str = "<x>";
target = rod.getPos().norm()*exp(-1.0*time);
switch (force) {
case Force::zero_F:
target = rod.getPos().sum();
break;
case Force::const_F:
target = static_cast<double>(NAN);
break;
case Force::harmonic_F:
target = rod.getPos().sum() * exp(-1 * time);
break;
}
} break;
case Type::x_squared: {
resetting = false;
type_str = "<x_squared>";
target = rod.getPos().squaredNorm()*exp(-2.0*time)+
2.0*(1-exp(-2.0*time));
switch (force) {
case Force::zero_F:
target = static_cast<double>(NAN);
break;
case Force::const_F:
target = static_cast<double>(NAN);
break;
case Force::harmonic_F:
target = rod.getPos().squaredNorm() * exp(-2 * time) +
2 * (1 - exp(-2 * time));
break;
}
} break;
}
}

35
C++/src/Compute.h
View File

@ -12,22 +12,26 @@
class Compute {
public:
enum class Type { msd, oaf, empxx, empyy, msd_harmonic_k_1,msd_const, x,x_squared };
enum class Type {
msd,
oaf,
empxx,
empyy,
x,
x_squared
};
enum class Force {
zero_F,
const_F,
harmonic_F
};
explicit Compute(Rod2d rod, Type t_type, size_t t_every, Simulation &sim);
explicit Compute(Rod2d rod, Type t_type, Force force,size_t t_every, Simulation &sim);
void eval(const Rod2d &rod2D);
void evalMSD(const Rod2d &rod2D);
void evalOAF(const Rod2d &rod2D);
void eval_empXX(const Rod2d &rod2D);
[[nodiscard]] double getTarget() const;
void eval_empYY(const Rod2d &rod2D);
[[nodiscard]] const LiveAgg &getAgg() const;
[[nodiscard]] Type getType() const;
@ -40,7 +44,7 @@ class Compute {
void reset(const Rod2d &rod);
private:
private:
Rod2d start_rod;
LiveAgg agg;
size_t every;
@ -49,11 +53,14 @@ private:
double target;
double time;
std::string type_str;
bool resetting;
void evalX(const Rod2d &rod2D);
bool resetting;
void evalMSD(const Rod2d &rod2D);
void evalOAF(const Rod2d &rod2D);
void eval_empXX(const Rod2d &rod2D);
void eval_empYY(const Rod2d &rod2D);
void evalX_squared(const Rod2d &rod2D);
};
#endif // MYPROJECT_COMPUTE_H

2
C++/src/Rod2d.cpp
View File

@ -10,7 +10,7 @@ void sqrt(Eigen::Matrix2d &mat) {
mat.data()[i] = sqrt(mat.data()[i]);
}
}
Rod2d::Rod2d(double L) : m_pos({0, 0}), m_e({1, 0}) {
Rod2d::Rod2d(double L) : m_pos({1, 1}), m_e({1, 0}) {
assert(L >= 1.0);
if (L == 1.0) {
m_D_rot = 3.0;

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

@ -11,12 +11,14 @@
#include "force_lambdas.h"
constexpr size_t SEED = 1234;
constexpr double stepSize = 0.001;
constexpr int n_computes = 100;
constexpr size_t num_runs = 5;
constexpr int numStep = 1000000;
constexpr int delta_compute = 4;
void run_no_force(size_t integrator_index, size_t force_index,
size_t length_index) {
constexpr size_t n_computes = 100;
constexpr size_t num_integratoren = 5;
constexpr size_t num_forces = 3;
constexpr size_t num_length = 3;
constexpr size_t delta_compute = 4;
constexpr size_t numStep = 10000;
inline const std::string header = {"time,val,target,std\n"};
void run(size_t integrator_index, size_t force_index, size_t length_index) {
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"),
@ -30,90 +32,100 @@ void run_no_force(size_t integrator_index, size_t force_index,
{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)+"_";
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];
{
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 * delta_compute);
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);
}
Calculation euler(integrator, computes, stepSize, SEED, force,
zero_Torque, length);
force_type, zero_Torque, length);
euler.run(numStep);
for (size_t i = 0; i < numStep; ++i) {
euler.run(n_computes * delta_compute);
}
std::ofstream msdFile(folder + "msd.dat");
std::ofstream oafFile(folder + "oaf.dat");
std::ofstream empxxFile(folder + "empxx.dat");
std::ofstream empyyFile(folder + "empyy.dat");
msdFile << "time,val,target\n";
oafFile << "time,val,target\n";
empxxFile << "time,val,target\n";
empyyFile << "time,val,target\n";
msdFile << header;
oafFile << header;
empxxFile << header;
empyyFile << header;
for (const auto &com : euler.getComputes()) {
if (com.getType() == MSD_VARIANT) {
if (com.getType() == Compute::Type::msd) {
msdFile << com.getTime() << ", " << com.getAgg().getMean()
<< ", " << com.getTarget() << std::endl;
<< ", " << com.getTarget() << ", "
<< com.getAgg().getSEM() << std::endl;
}
if (com.getType() == Compute::Type::oaf) {
oafFile << com.getTime() << ", " << com.getAgg().getMean()
<< ", " << com.getTarget() << std::endl;
<< ", " << com.getTarget() << ", "
<< com.getAgg().getSEM() << std::endl;
}
if (com.getType() == Compute::Type::empxx) {
empxxFile << com.getTime() << ", " << com.getAgg().getMean()
<< ", " << com.getTarget() << std::endl;
<< ", " << com.getTarget() << ", "
<< com.getAgg().getSEM() << std::endl;
}
if (com.getType() == Compute::Type::empyy) {
empyyFile << com.getTime() << ", " << com.getAgg().getMean()
<< ", " << com.getTarget() << std::endl;
<< ", " << com.getTarget() << ", "
<< com.getAgg().getSEM() << std::endl;
}
if (com.getType() == Compute::Type::msd) {
msdFile << com.getTime() << ", " << com.getAgg().getMean()
<< ", " << com.getTarget() << std::endl;
<< ", " << com.getTarget() << ", "
<< com.getAgg().getSEM() << 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 * delta_compute);
repeating_computes.emplace_back(Compute::Type::x_squared, i * delta_compute);
for (size_t i = 1; i < n_computes; ++i) {
repeating_computes.emplace_back(Compute::Type::x,
i * delta_compute);
repeating_computes.emplace_back(Compute::Type::x_squared,
i * delta_compute);
}
Calculation calc_repeat(integrator, repeating_computes, stepSize, SEED,
force, zero_Torque, 1.0);
for (int i = 0; i < numStep / n_computes; ++i) {
calc_repeat.run(n_computes);
force, force_type, zero_Torque, 1.0);
for (size_t i = 0; i < numStep; ++i) {
calc_repeat.run(n_computes * delta_compute);
calc_repeat.reset();
}
std::ofstream xFile(folder + "x.dat");
std::ofstream xsqFile(folder + "x_squared.dat");
xFile << "time,val,target\n";
xsqFile << "time,val,target\n";
xFile << header;
xsqFile << header;
for (const auto &com : calc_repeat.getComputes()) {
if (com.getType() == Compute::Type::x) {
xFile << com.getTime() << ", " << com.getAgg().getMean() << ", "
<< com.getTarget() << std::endl;
<< com.getTarget() << ", " << com.getAgg().getSEM()
<< std::endl;
}
if (com.getType() == Compute::Type::x_squared) {
xsqFile << com.getTime() << ", " << com.getAgg().getMean()
<< ", " << com.getTarget() << std::endl;
<< ", " << com.getTarget() << ", "
<< com.getAgg().getSEM() << std::endl;
}
}
}
@ -122,11 +134,10 @@ void run_no_force(size_t integrator_index, size_t force_index,
}
int main() {
for (size_t j = 0; j < 3; ++j) {
for (size_t i = 0; i < num_runs; ++i) {
for (size_t k = 0; k < 1; ++k) {
run_no_force(i, j, k);
}
}
#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));
}
}