Barnes-Hut 시뮬레이션은 n체 운동 시뮬레이션을 $O(n \log(n))$ 의 시간 복잡도로 풀 수 있게 해주는 알고리즘임.
기존 n체 시뮬레이션은 한 입자와 다른 n개의 입자가 서로 끌어당기는 힘을 계산하므로 n개의 입자가 있다면 시간 복잡도는 $n^2$임. 태양계처럼 n이 9 정도라면 현대 컴퓨터로는 1초에 수만 번 계산을 할 수 있겠지만, 만약 n이 500 이상이라면 시간복잡도 $n^2$의 시간복잡도로는 수 초가 걸림. 이런 연산은 핵융합 시뮬레이션이나 분자동역학 시뮬레이션, 우주 시뮬레이션(은하나 토성의 고리 같은 거) 등에 사용되고 보통 슈퍼컴퓨터가 이런 연산을 수행함.
Barnes-Hut 시뮬레이션은 쿼드트리라는 자료구조를 사용해서 연산의 시간 복잡도를 줄이는데, 방법은 이럼:
일단 공간을 쿼드트리로 재귀적으로 사분할(3D일 경우 8분할)함.
그럼 위 사진처럼 입자들이 각각 박스에 들어감.
평범한 n체 시뮬레이션이라면 입자끼리 계산했겠지만, 여기서는 박스마다 CM(질량 중심)과 총 질량을 계산함.
이 알고리즘의 핵심은 바로 어느 정도의 박스까지 계산을 할 것인가? 임. 만약 모든 기본 단위 박스(입자가 1개짜리인 박스)에 대해 계산한다면 기존처럼 n^2의 시간 복잡도를 가질 것이고 일정 크기 이하의 박스를 무시한다면 1개의 입자에 대한 힘을 계산할 떄 트리 탐색 시간복잡도와 같은 $\log(n)$의 시간 복잡도를 갖고 n개의 입자를 계산해야 하므로 총 시간복잡도는 $O(n \log(n))$임.
계산할 박스의 크기는 다음과 같이 계산함.
위의 그림에서 입자와 박스 질량 중심 사이의 거리 r과 박스의 길이 d가 주어져 있으면 $\frac{d}{r}$ 이 일정 값(보통 0.5나 1) 이상이 되도록 하면 됨.
코드로 나타내면 다음과 같음:
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 - tmp->CenterOfMass;
double Dist = Dist_V.norm ();
if (tmp->width / Dist <= THETA || !tmp->HasLeaf)
if (!tmp->contains (leaf))
NetAcc += gravity (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;
}
코드는 트리를 BFS 탐색함. 트리 탐색 중에 THETA 이하라면(or 그 아래 노드가 비어 있다면) 해당 박스가 입자에 작용하는 가속도를 더해줌.
이대로 코드를 짜고 실행해 보면 엄청난 속도로 튕겨나가는 입자가 생기는데, 그건 입자들 사이의 간격이 지나치게 가까워져서 빠른 속도로 경계를 향해 날아가고 이렇게 날아간 입자들은 세그멘테이션 폴트를 일으키는데 방지하지 위해서는 경계에서 입자를 반사해주는 함수와 Force Softening이라는 게 필요함. 중력에 의한 힘 식 $F = G\frac{mM}{|r|^3} \vec{r}$ 에서 r^3이 0으로 수렴하고 힘은 무한대로 발산하게 되는데 $r^3$을 $(r^2 + \epsilon^2)^{1.5}$ 으로 두고 엡실론을 적당히 작은 수로 두면 힘이 무한대로 발산하는 걸 막을 수 있음.
class Gravitational
{
public:
vec2 operator()( const double RootMass, const vec2 &dist)
{
double norm = dist.norm ();
return -(GRAVITY_CONST * RootMass)
/ (std::pow ((norm * norm) + (EPSILON * EPSILON), 1.5)) * dist;
}
};
이전에 이중 진자 시뮬을 만들 때는 룽게-쿠타 방법으로 미분 방정식을 풀었는데 이런 시뮬레이션에 룽게-쿠타 방법을 사용하는 건 지나치게 많은 계산을 요구하는데다가, 에너지 보존 법칙이 성립하지 않기 때문에(부동소수점 오차 등으로 위상 공간 에서 그래프로 그려보면 나선을 그리면서 에너지가 손실되거나 증가함) 사용하면 안됨. 쓴다고 해도 4차가 아닌 훨씬 고차 방법이 필요함. 여기서는 가장 간단한 Symplectic Integrator인 Semi-Implicit Euler방법을 사용하였음.
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);
}
};
}
코드는 간단함. 단지 라그랑주 역학과 해밀토니안을 알아야 이해할 수 있을 뿐임.
실제 타원 은하의 항성 분포인 Plummer Model 등을 적용해야 하는데 이걸 어떻게 해야 할지 몰라서 거대 질량 천체 주변을 도는 작은 질량의 천체들로 시뮬레이션했음.
void DiskDistribution(
std::vector<Body> &bodies, int number, double rad_max, double rad_min,
double mass_min, double mass_max, double centermass)
{
std::random_device rand;
std::mt19937 gen (rand ());
std::uniform_real_distribution<> rand_angle (0, 2 * M_PI);
std::uniform_real_distribution<> rand_dist (rad_min, rad_max);
std::uniform_real_distribution<> rand_mass (mass_min, mass_max);
bodies.emplace_back(vec2(SimWidth / 2.0, SimHeight / 2.0), vec2(0.0, 0.0), centermass);
for (int i = 1; i < number; ++i)
{
double angle = rand_angle (gen);
double distance = rand_dist (gen);
double pos_x = std::cos (angle) * distance + (ViewWidth / 2.0);
double pos_y = std::sin (angle) * distance + (ViewHeight / 2.0);
double velocity = std::pow (centermass * GRAVITY_CONST / distance, 0.5);
double vel_x = std::sin (angle) * velocity;
double vel_y = -std::cos (angle) * velocity;
double mass = rand_mass (gen);
bodies.emplace_back (vec2 (pos_x, pos_y), vec2 (vel_x, vel_y), mass);
}
}How to randomly distribute N masses, such that they follow a Plummer density distribution in Python
I am working in Python. I have N stars, each of one solar mass M_0. I want to randomly distribute these stars in a volume of radius R, such that the density distribution follows a Plummer model, gi...
stackoverflow.com
여기에 나와있긴 함
위 설명은 대부분 이곳 과 이곳을 참고하였음.
코드의 성능은 Ryzen 7 6800HS CPU + 32GB 램 위에서 측정하였고, Visual Studio 2022에서 OpenMP 16스레드 멀티코어 프로세싱, 최대 최적화 옵션, SIMD AVX2 활성화해서 측정했고, 1000개의 파티클을 시뮬레이팅 하면 루프 당 15밀리세컨드가 필요하였음.
작업관리자 CPU 100% 찍어봄
GPGPU 버전 만들어봐야겠음. AWS에서 서버 빌려서 한 번 돌려보고싶음.
Plummer Model 만들어서 이쁘게 나오면 유튜브에 올릴것임.
코드 깃허브:
https://github.com/hydrogendeuteride/Barnes_Hut_Demo
GitHub - hydrogendeuteride/Barnes_Hut_Demo: 2D Barnes-Hut galaxy simulation
2D Barnes-Hut galaxy simulation. Contribute to hydrogendeuteride/Barnes_Hut_Demo development by creating an account on GitHub.
github.com
'Physics > 전산물리' 카테고리의 다른 글
| 1차원 슈뢰딩거 방정식 시뮬레이션 (0) | 2023.08.28 |
|---|---|
| 2차원 파동방정식 FDM (0) | 2023.07.01 |
| Barnes-Hut 시뮬레이션 진행상황 2 (0) | 2023.03.03 |
| 조금 더 깔끔한 이중진자 시뮬레이션 (1) | 2023.02.17 |
| Barnes-Hut 시뮬레이션 진행 상황 (0) | 2023.02.08 |
댓글