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

Added EMPXX and EMPYY

Browse Source
This commit is contained in:
Jacob Holder 2021-10-25 15:34:34 +02:00
parent 81f4796eeb
commit 7eac132a5f
Signed by: jacob
GPG Key ID: 2194FC747048A7FD
4 changed files with 143 additions and 91 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

142
src/Compute.cpp
View File

@ -9,85 +9,109 @@
#include "Rod2d.hpp" #include "Rod2d.hpp"
#include "Simulation.h" #include "Simulation.h"
void Compute::evalMSD(const Rod2d &rod2D) { void Compute::evalMSD(const Rod2d &rod2D)
const auto &new_pos = rod2D.getPos(); {
auto old_pos = start_rod.getPos(); const auto &new_pos = rod2D.getPos();
auto msd = (new_pos - old_pos).squaredNorm(); auto old_pos = start_rod.getPos();
agg.feed(msd); auto msd = (new_pos - old_pos).squaredNorm();
agg.feed(msd);
} }
void Compute::evalOAF(const Rod2d &rod2D) { void Compute::evalOAF(const Rod2d &rod2D)
const auto &new_e = rod2D.getE(); {
auto old_e = start_rod.getE(); const auto &new_e = rod2D.getE();
auto oaf = old_e.dot(new_e); auto old_e = start_rod.getE();
agg.feed(oaf); auto oaf = old_e.dot(new_e);
agg.feed(oaf);
} }
void Compute::eval_empXX(const Rod2d & /*rod2D*/) { void Compute::eval_empXX(const Rod2d &rod2D)
//TODO {
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);
agg.feed(empxx);
} }
void Compute::eval_empYY(const Rod2d & /*rod2D*/) { void Compute::eval_empYY(const Rod2d &rod2D)
//TODO {
const auto &new_pos = rod2D.getPos();
auto old_pos = start_rod.getPos();
auto old_e_ortho = start_rod.getE_ortho();
double empxx = pow((new_pos - old_pos).dot(old_e_ortho), 2.0);
agg.feed(empxx);
} }
void Compute::eval(const Rod2d &rod2D) { void Compute::eval(const Rod2d &rod2D)
time_step++; {
if (time_step % every == 0) { time_step++;
switch (type) { if (time_step % every == 0) {
case Type::msd: switch (type) {
evalMSD(rod2D); case Type::msd:
break; evalMSD(rod2D);
case Type::oaf: break;
evalOAF(rod2D); case Type::oaf:
break; evalOAF(rod2D);
case Type::empxx: break;
eval_empXX(rod2D); case Type::empxx:
break; eval_empXX(rod2D);
case Type::empyy: break;
eval_empYY(rod2D); case Type::empyy:
break; eval_empYY(rod2D);
} break;
start_rod = rod2D;
} }
start_rod = rod2D;
}
} }
Compute::Compute(Rod2d rod, Type t_type, size_t t_every, Simulation &sim) 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) { : start_rod(std::move(rod)), every(t_every), time_step(0), type(t_type)
auto time = sim.getMDeltaT() * static_cast<double>(every); {
switch (type) { auto time = sim.getMDeltaT() * static_cast<double>(every);
case Type::msd: switch (type) {
target = 4.0 * 0.5 * (rod.getDiff().trace()) * time; case Type::msd:
break; target = 4.0 * 0.5 * (rod.getDiff().trace()) * time;
case Type::oaf: break;
target = std::exp(-rod.getDRot() * time); case Type::oaf:
break; target = std::exp(-rod.getDRot() * time);
case Type::empxx: break;
target = 0; case Type::empxx: {
break; const double Dmean = 0.5 * (rod.getDiff().trace());
case Type::empyy: const double u = 4.0;
target = 0; const double deltaD = rod.getDiff()(1, 1) - rod.getDiff()(0, 0);
break; target = Dmean - deltaD / 2.0 * (1 - exp(-u * rod.getDRot() * time)) / u / rod.getDRot() / time;
} } break;
case Type::empyy: {
const double Dmean = 0.5 * (rod.getDiff().trace());
const double u = 4.0;
const double deltaD = rod.getDiff()(1, 1) - rod.getDiff()(0, 0);
target = Dmean + deltaD / 2.0 * (1 - exp(-u * rod.getDRot() * time)) / u / rod.getDRot() / time;
} break;
}
} }
const LiveAgg &Compute::getAgg() const { const LiveAgg &Compute::getAgg() const
return agg; {
return agg;
} }
Compute::Type Compute::getType() const { Compute::Type Compute::getType() const
return type; {
return type;
} }
double Compute::getTime() const { double Compute::getTime() const
return static_cast<double>(every); {
return static_cast<double>(every);
} }
double Compute::getTarget() const { double Compute::getTarget() const
return target; {
return target;
} }
double Compute::getDifference() const { double Compute::getDifference() const
return abs(agg.getMean() -target); {
return abs(agg.getMean() - target);
} }

49
src/Compute.h
View File

@ -10,45 +10,48 @@
#include "Rod2d.hpp" #include "Rod2d.hpp"
#include "Simulation.h" #include "Simulation.h"
class Compute { class Compute
{
public: public:
enum class Type { enum class Type {
msd, oaf, empxx, empyy msd,
}; oaf,
empxx,
empyy
};
explicit Compute(Rod2d rod, Type t_type, size_t t_every, Simulation &sim); explicit Compute(Rod2d rod, Type t_type, size_t t_every, Simulation &sim);
void eval(const Rod2d &rod2D); void eval(const Rod2d &rod2D);
void evalMSD(const Rod2d &rod2D); void evalMSD(const Rod2d &rod2D);
void evalOAF(const Rod2d &rod2D); void evalOAF(const Rod2d &rod2D);
void eval_empXX(const Rod2d &rod2D); void eval_empXX(const Rod2d &rod2D);
double getTarget() const; double getTarget() const;
void eval_empYY(const Rod2d &rod2D); void eval_empYY(const Rod2d &rod2D);
[[nodiscard]] const LiveAgg &getAgg() const; [[nodiscard]] const LiveAgg &getAgg() const;
[[nodiscard]] Type getType() const; [[nodiscard]] Type getType() const;
double getTime() const; double getTime() const;
double getDifference() const; double getDifference() const;
private: private:
Rod2d start_rod; Rod2d start_rod;
LiveAgg agg; LiveAgg agg;
size_t every; size_t every;
size_t time_step; size_t time_step;
Type type; Type type;
double target; double target;
}; };
#endif //MYPROJECT_COMPUTE_H #endif// MYPROJECT_COMPUTE_H

32
test/test_Integratoren.cpp
View File

@ -16,7 +16,9 @@ TEST_CASE("Forceless Integratoren") {
SECTION("Euler") { SECTION("Euler") {
Calculation euler(Integratoren2d_forceless::Set1_Euler, {{Compute::Type::msd, 1}, Calculation euler(Integratoren2d_forceless::Set1_Euler, {{Compute::Type::msd, 1},
{Compute::Type::oaf, 1}}, 0.01, SEED); {Compute::Type::oaf, 1},
{Compute::Type::empxx, 1},
{Compute::Type::empyy, 1}}, 0.01, SEED);
{ {
euler.run(10000); euler.run(10000);
for (auto &c: euler.getComputes()) { CHECK(c.getDifference() <= 0.005); } for (auto &c: euler.getComputes()) { CHECK(c.getDifference() <= 0.005); }
@ -34,14 +36,18 @@ TEST_CASE("Forceless Integratoren") {
}SECTION("Heun") { }SECTION("Heun") {
Calculation heun(Integratoren2d_forceless::Set2_Heun, {{Compute::Type::msd, 1}, Calculation heun(Integratoren2d_forceless::Set2_Heun, {{Compute::Type::msd, 1},
{Compute::Type::oaf, 1}}, 0.01, SEED); {Compute::Type::oaf, 1},
{Compute::Type::empxx, 1},
{Compute::Type::empyy, 1}}, 0.01, SEED);
{ {
heun.run(10000); heun.run(10000);
for (auto &c: heun.getComputes()) { CHECK(c.getDifference() <= 0.005); } for (auto &c: heun.getComputes()) { CHECK(c.getDifference() <= 0.005); }
} }
Calculation heun_force(Integratoren2d_force::Set2_Heun, {{Compute::Type::msd, 1}, Calculation heun_force(Integratoren2d_force::Set2_Heun, {{Compute::Type::msd, 1},
{Compute::Type::oaf, 1}}, 0.01, SEED, force, {Compute::Type::oaf, 1},
{Compute::Type::empxx, 1},
{Compute::Type::empyy, 1}}, 0.01, SEED, force,
torque); torque);
{ {
heun_force.run(10000); heun_force.run(10000);
@ -52,14 +58,18 @@ TEST_CASE("Forceless Integratoren") {
}SECTION("Exact") { }SECTION("Exact") {
Calculation exact(Integratoren2d_forceless::Set3_Exact, {{Compute::Type::msd, 1}, Calculation exact(Integratoren2d_forceless::Set3_Exact, {{Compute::Type::msd, 1},
{Compute::Type::oaf, 1}}, 0.01, SEED); {Compute::Type::oaf, 1},
{Compute::Type::empxx, 1},
{Compute::Type::empyy, 1}}, 0.01, SEED);
{ {
exact.run(10000); exact.run(10000);
for (auto &c: exact.getComputes()) { CHECK(c.getDifference() <= 0.005); } for (auto &c: exact.getComputes()) { CHECK(c.getDifference() <= 0.005); }
} }
Calculation exact_force(Integratoren2d_force::Set3_Exact, {{Compute::Type::msd, 1}, Calculation exact_force(Integratoren2d_force::Set3_Exact, {{Compute::Type::msd, 1},
{Compute::Type::oaf, 1}}, 0.01, SEED, force, {Compute::Type::oaf, 1},
{Compute::Type::empxx, 1},
{Compute::Type::empyy, 1}}, 0.01, SEED, force,
torque); torque);
{ {
exact_force.run(10000); exact_force.run(10000);
@ -70,20 +80,24 @@ TEST_CASE("Forceless Integratoren") {
}SECTION("BDAS") { }SECTION("BDAS") {
Calculation bdas(Integratoren2d_forceless::Set4_BDAS, {{Compute::Type::msd, 1}, Calculation bdas(Integratoren2d_forceless::Set4_BDAS, {{Compute::Type::msd, 1},
{Compute::Type::oaf, 1}}, 0.01, SEED); {Compute::Type::oaf, 1},
{Compute::Type::empxx, 1},
{Compute::Type::empyy, 1}}, 0.01, SEED);
{ {
bdas.run(10000); bdas.run(10000);
for (auto &c: bdas.getComputes()) { CHECK(c.getDifference() <= 0.005); } for (auto &c: bdas.getComputes()) { CHECK(c.getDifference() <= 0.005); }
} }
Calculation bdas_force(Integratoren2d_force::Set4_BDAS, {{Compute::Type::msd, 1}, Calculation bdas_force(Integratoren2d_force::Set4_BDAS, {{Compute::Type::msd, 1},
{Compute::Type::oaf, 1}}, 0.01, SEED, force, {Compute::Type::oaf, 1},
{Compute::Type::empxx, 1},
{Compute::Type::empyy, 1}}, 0.01, SEED, force,
torque); torque);
{ {
bdas_force.run(10000); bdas_force.run(10000);
for (auto &c: bdas_force.getComputes()) { CHECK(c.getDifference() <= 0.005); } for (auto &c: bdas_force.getComputes()) { CHECK(c.getDifference() <= 0.005); }
} }
CHECK(bdas.getComputes()[0].getAgg().getMean() == bdas_force.getComputes()[0].getAgg().getMean()); CHECK(bdas.getComputes()[0].getAgg().getMean() == Approx(bdas_force.getComputes()[0].getAgg().getMean()).epsilon(10e-10));
CHECK(bdas.getComputes()[1].getAgg().getMean() == bdas_force.getComputes()[1].getAgg().getMean()); CHECK(bdas.getComputes()[1].getAgg().getMean() == Approx(bdas_force.getComputes()[1].getAgg().getMean()).epsilon(10e-10));
} }
} }

11
test/test_Rod2d.cpp
View File

@ -43,5 +43,16 @@ TEST_CASE("Rods") {
} }
SECTION("Angle Stuff"){
const Eigen::Vector2d unitVec = {1, 0};
double f1 = GENERATE(take(10,random(-100,100)));
double f2 = GENERATE(take(10,random(0,100)));
sphere.setE(Eigen::Vector2d({f1,f2}).normalized());
Eigen::Rotation2D rotMat(acos(unitVec.dot(sphere.getE())));
CHECK(rotMat.toRotationMatrix()(0,0) == Approx(sphere.getE_Base_matrix()(0,0)).epsilon(10e-10));
CHECK(rotMat.toRotationMatrix()(1,0) == Approx(sphere.getE_Base_matrix()(1,0)).epsilon(10e-10));
CHECK(rotMat.toRotationMatrix()(0,1) == Approx(sphere.getE_Base_matrix()(0,1)).epsilon(10e-10));
CHECK(rotMat.toRotationMatrix()(1,1) == Approx(sphere.getE_Base_matrix()(1,1)).epsilon(10e-10));
}
} }