Barnes Hut 시뮬레이션은 은하나 우주 거대 구조, 분자동역학처럼 매우 많은(대충 1000이상) 입자의 n체 시뮬레이션을 하려고 만든 알고리즘임.
이제 GUI, 멀티코어, 올바른 상수값 입력 정도만 완성하면 되는데 그거 만들기 전에 기술실증용으로 2중 진자 시뮬레이션 먼저 만들거임.
저 세개중에서는 멀티코어 프로세싱이 경험도 없고 제일 어려운건데 왜 이 시뮬레이션에 멀티코어가 필수인지는 아래 설명하겠음.
일단 코드(Eigen 라이브러리 필요, C++17이상)
#ifndef BARNES_HUT_DEMO_QUADTREE_H
#define BARNES_HUT_DEMO_QUADTREE_H
#include <memory>
#include <eigen3/Eigen/Dense>
#include <utility>
#include <vector>
#include <queue>
#include <algorithm>
#include <stack>
#include <numeric>
typedef Eigen::Vector2d vec2;
constexpr double PRECISION = 1e-15;
struct Body
{
vec2 pos;
vec2 vel;
double mass;
Body(vec2 p, vec2 v, double m) : pos (std::move(p)), vel (std::move(v)), mass (m)
{};
bool operator==(const Body& x)const
{
return (pos.isApprox (x.pos, PRECISION)) &&
(vel.isApprox (x.vel, PRECISION)) &&
(abs((mass - x.mass)) <= PRECISION);
}
};
class Node
{
public:
std::vector<Body> data;
std::shared_ptr<Node> q1;
std::shared_ptr<Node> q2;
std::shared_ptr<Node> q3;
std::shared_ptr<Node> q4;
vec2 CenterOfMass;
double TotalMass = 0;
double width = 0, height = 0;
double BoxPosX = 0, BoxPosY = 0;
bool HasLeaf = false;
Node(std::vector<Body> data, double w, double h);
Node(std::vector<Body> data, double w, double h, double BoxX, double BoxY);
bool contains(const Body& body);
};
class QuadTree
{
public:
QuadTree() : root (nullptr)
{};
void addNodeIterative(std::vector<Body> &data, double w, double h);
void reset();
//private:
std::shared_ptr<Node> root;
};
#endif //BARNES_HUT_DEMO_QUADTREE_Hquadtree.h
#include "quadtree.h"
Node::Node(std::vector<Body> data, double w, double h)
: data (std::move (data)), width (w), height (h)
{
TotalMass = std::accumulate (this->data.begin (), this->data.end (), 0,
[](double sum, const auto &x) { return sum + x.mass; });
CenterOfMass = std::accumulate (this->data.begin (), this->data.end (), vec2 (0.0, 0.0),
[](const vec2 &sum, const auto &x) { return sum + (x.mass * x.pos); }) / TotalMass;
}
Node::Node(std::vector<Body> data, double w, double h, double BoxX, double BoxY)
: data (std::move (data)), width (w), height (h), BoxPosX (BoxX), BoxPosY (BoxY)
{
TotalMass = std::accumulate (this->data.begin (), this->data.end (), 0,
[](double sum, const auto &x) { return sum + x.mass; });
CenterOfMass = std::accumulate (this->data.begin (), this->data.end (), vec2 (0.0, 0.0),
[](const vec2 &sum, const auto &x) { return sum + (x.mass * x.pos); }) / TotalMass;
}
bool Node::contains(const Body &body)
{
if (std::find (this->data.begin (), this->data.end (), body) == this->data.end ())
return false;
else
return true;
}
void QuadTree::reset()
{
std::stack<std::shared_ptr<Node>> stack1, stack2;
std::shared_ptr<Node> tmp;
stack1.push (root);
/*
* leaf4 -> pop
* leaf3 (leaf3)
* leaf2 (leaf4)
* leaf1 -(leaf4)-> root
* ______ ______
* stack1 stack2
*/
while (!stack1.empty ())
{
tmp = stack1.top ();
stack1.pop ();
stack2.push (tmp);
if (tmp->q1 != nullptr)
stack1.push (tmp->q1);
if (tmp->q2 != nullptr)
stack1.push (tmp->q2);
if (tmp->q3 != nullptr)
stack1.push (tmp->q3);
if (tmp->q4 != nullptr)
stack1.push (tmp->q4);
}
while (!stack2.empty())
{
tmp = stack2.top();
stack2.pop();
tmp->HasLeaf = false;
tmp.reset();
}
root.reset();
}
void QuadTree::addNodeIterative(std::vector<Body> &data, double w, double h)
{
std::queue<std::shared_ptr<Node>> queue;
if (root != nullptr)
{
root->HasLeaf = true;
queue.push (root);
}
else
{
root = std::make_shared<Node> (data, w, h, 0, 0);
root->HasLeaf = true;
queue.push (root);
}
while (!queue.empty ())
{
std::shared_ptr<Node> tmp = queue.front ();
queue.pop ();
std::vector<Body> q1, q2, q3, q4;
for (auto &body1: tmp->data)
{
if (body1.pos (0, 0) < (tmp->BoxPosX + (tmp->width / 2.0)))
{
if (body1.pos (1, 0) < (tmp->BoxPosY) + (tmp->height / 2.0))
q1.push_back (body1);
else
q3.push_back (body1);
}
else
{
if (body1.pos (1, 0) < (tmp->BoxPosY) + (tmp->height / 2.0))
q2.push_back (body1);
else
q4.push_back (body1);
}
}
/*
*0 1 2 3 4 5 6 7
*1 |
*2 q1 | q2
*3 |
*4----------|--------------
*5 |
*6 q3 | q4
*7 |
*/
if (!q1.empty ())
{
tmp->q1 = std::make_shared<Node> (q1, tmp->width / 2.0, tmp->height / 2.0,
tmp->BoxPosX, tmp->BoxPosY);
if (q1.size () > 1)
{
queue.push (tmp->q1);
tmp->HasLeaf = true;
}
}
if (!q2.empty ())
{
tmp->q2 = std::make_shared<Node> (q2, tmp->width / 2.0, tmp->height / 2.0,
tmp->BoxPosX + tmp->width / 2.0, tmp->BoxPosY);
if (q2.size () > 1)
{
queue.push (tmp->q2);
tmp->HasLeaf = true;
}
}
if (!q3.empty ())
{
tmp->q3 = std::make_shared<Node> (q3, tmp->width / 2.0, tmp->height / 2.0,
tmp->BoxPosX, tmp->BoxPosY + tmp->height / 2.0);
if (q3.size () > 1)
{
queue.push (tmp->q3);
tmp->HasLeaf = true;
}
}
if (!q4.empty ())
{
tmp->q4 = std::make_shared<Node> (q4, tmp->width / 2.0, tmp->height / 2.0,
tmp->BoxPosX + tmp->width / 2.0, tmp->BoxPosY + tmp->height / 2.0);
if (q4.size () > 1)
{
queue.push (tmp->q4);
tmp->HasLeaf = true;
}
}
}
}quadtree.cc
#include <iostream>
#include <random>
#include <tuple>
#include <chrono>
#include "quadtree.h"
void postorder(const std::shared_ptr<Node> &Node)
{
if (Node == nullptr)
return;
postorder (Node->q1);
postorder (Node->q2);
postorder (Node->q3);
postorder (Node->q4);
std::cout << Node->width << " " << Node->height << " " << Node->data.size () << std::endl;
}
constexpr double GRAVITY_CONST = 0.1;
constexpr double SOFTENER = 5.0;
constexpr double THETA = 0.5;
namespace Acceleration
{
class Gravitational
{
public:
vec2 operator()(const double LeafMass, const double RootMass, const vec2 &dist)
{
return -(GRAVITY_CONST * LeafMass * RootMass) / (std::pow (dist.norm (), 3)) * dist;
}
};
}
namespace Integrator
{
class Semi_Implicit_Euler
{
public:
std::tuple<vec2, vec2> operator()(
const std::tuple<vec2, vec2> &Pos_Vel,
const vec2 &accel, const double timestep)
{
auto [x, v] = Pos_Vel;
vec2 x_1 (0.0, 0.0), v_1 (0.0, 0.0);
v_1 = v + accel * timestep;
x_1 = x + v_1 * timestep;
return std::make_tuple (x_1, v_1);
}
};
}
class BarnesHutTree : public QuadTree
{
private:
vec2 NetAcceleration(Body &leaf);
Acceleration::Gravitational gravity;
public:
void CalcMovement(Body &body, double timestep);
};
vec2 BarnesHutTree::NetAcceleration(Body &leaf)
{
vec2 NetAcc = vec2 (0.0, 0.0);
std::queue<std::shared_ptr<Node>> queue;
queue.push (root);
while (!queue.empty ())
{
std::shared_ptr<Node> tmp = queue.front ();
queue.pop ();
vec2 Dist_V = leaf.pos - root->CenterOfMass;
double Dist = Dist_V.norm ();
if (tmp->width / Dist <= THETA || !tmp->HasLeaf)
if (!tmp->contains (leaf))
NetAcc += gravity (leaf.mass, tmp->TotalMass, Dist_V);
if (tmp->q1 != nullptr)
queue.push (tmp->q1);
if (tmp->q2 != nullptr)
queue.push (tmp->q2);
if (tmp->q3 != nullptr)
queue.push (tmp->q3);
if (tmp->q4 != nullptr)
queue.push (tmp->q4);
}
return NetAcc;
}
void BarnesHutTree::CalcMovement(Body &body, double timestep)
{
Integrator::Semi_Implicit_Euler euler;
auto [x, v] =
euler (std::make_tuple (body.pos, body.vel),
NetAcceleration (body), timestep);
body.pos = x, body.vel = v;
}
int main()
{
std::random_device rand;
std::mt19937 gen (rand ());
std::uniform_real_distribution<> dist (0.0, 10.0);
std::uniform_real_distribution<> mass (1.0, 2.0);
std::vector<Body> bodies;
for (int i = 0; i < 1000; ++i)
bodies.emplace_back (vec2 (dist (gen), dist (gen)),
vec2 (dist (gen), dist (gen)),
mass (gen));
BarnesHutTree tree = BarnesHutTree ();
for (int i = 0; i < 20; ++i)
{
auto start = std::chrono::high_resolution_clock::now();
tree.addNodeIterative (bodies, 50, 50);
//postorder (tree.root);
for (auto& x:bodies)
{
tree.CalcMovement (x, 1.0);
}
auto stop = std::chrono::high_resolution_clock::now();
auto duration = std::chrono::duration_cast<std::chrono::milliseconds>(stop - start);
std::cout<<i<<"\t"<<duration.count()<<"\n";
//postorder (tree.root);
//tree.reset();
}
return 0;
}main.cc
여기서 조금 더 손봐야 될 요소들이 있긴 한데 일단은 코드는 이렇고 아직 깃허브에는 안올림.
코드를 실행시켜보면
0 3384
1 4469
2 4883
3 5028
4 5125
5 5170
6 5357
7 5343
8 5299
9 5576
10 5534
11 5360
12 5605
13 5524
14 5544
15 5409
16 5453
17 5404
18 5501
19 5485
이런 식으로 출력이 되는데 저기 찍힌 숫자가 밀리세컨드임.
컴퓨터 환경은 R7 6800HS, 1코어, 32G DDR5, Windows 11 + WSL2임.
1000개의 물체를 계산하는데 5.5초가 걸리므로 멀티코어 프로세싱이 잘 되서 16으로 나눈다 쳐도 한 프레임 당 0.3초 정도 소모됨.
프레임을 미리 계산해놓고 그 데이터로 GUI를 굴릴 예정인데 1초당 60프레임이라고 가정하면 1분에 3600장이 필요하고 계산하는 데 1080초(18분) 이상이 걸림. 데이터 크기는 48바이트 * 1000개 * 3600장 = 172800000바이트 = 대충 180메가바이트임.
아마 컴파일도 릴리즈로 하고 O2나 O3 같은 옵션 주면 10분 내로 할 수도 있고 아닐수도 있고.
코드는 이미 최적화할 수 있는만큼 해놔서 건드려봐도 속도 많이 안올라갈거임. std::move 쓸 수 있는 만큼 쓰고 재귀도 테스트 코드 빼면 하나도 안씀.
최적화 해볼만한 건 shared_ptr을 그냥 포인터로 쓰고 vector를 list로 바꾸는 건데 일단 완성하고 해볼거임. 첫번재 건 CUDA 버전 만들 때 할 건데 두번째 건 멀티코어 프로그래밍이 불가능할 것 갈음.
결과 저장할때는 Memory Mapped File 같은 거 써서 기록해야될듯함.
'Physics > 전산물리' 카테고리의 다른 글
| Barnes-Hut 시뮬레이션 진행상황 2 (0) | 2023.03.03 |
|---|---|
| 조금 더 깔끔한 이중진자 시뮬레이션 (1) | 2023.02.17 |
| 파이썬 궤도역학 연재 예정 (0) | 2021.03.22 |
| 파이썬 궤도 시뮬레이션 -1 원궤도 (1) | 2021.03.15 |
| N체 시뮬레이션 (0) | 2021.03.04 |
댓글