struct ParticleHash
{
uint particleID; // 원래 파티클 인덱스
uint hashValue; // 계산된 해시 값
uint flag; // 그룹 Flag
};
1. 위 구조체를 담는 배열을 생성하고, 아래와 같이 해시를 간단히 설정 -> 이후에 해시 생성 방식 변경
// 파티클 위치 p.position 을 그리드 최소 바운더리 위치를 기준으로 상대 위치로 변환
float3 relativePos = p.position - minBounds;
// 상대 위치를 셀 크기로 나누어, 각 축별로 몇 번째 셀에 해당하는지 실수 값으로 계산
float3 normalizedPos = relativePos / cellSize;
// floor 연산으로 소수점을 버려 정수 인덱스를 얻음 (이때 음수 인덱스 가능성 있음)
int3 cellID = int3(floor(normalizedPos));
uint hashValue = uint(cellID.x) +
uint(cellID.y) * gridDimX +
uint(cellID.z) * gridDimX * gridDimY;
ParticleHashesOutput[index].particleID = index;
ParticleHashesOutput[index].hashValue = hashValue;
ParticleHashesOutput[index].flag = 0;
2. 생성된 해시를 BitonicSort로 정렬하고(이후에 입자 개수가 많아질경우 Radix Sort 적용), 정렬된 해시 배열을 순회하면서, 이전 위치의 해시값과 지금 위치의 해시값이 다르면, 새로운 그룹으로 설정(flag = 1)
3. flag가 설정된 정렬된 해시 배열을 병렬 스캔 진행해서 전역 Prefix Sum 생성
효율적인 Belloch 병렬 스캔 알고리즘 적용(CUDA 예제 참고)
1. Pass 1 결과로 LocalScanSum이랑 PartialSum을 얻습니다.
2. Pass2에서는 PartialSum을 인풋으로, 결과를 blockOffsets를 얻습니다.
3. LocalScanSum이랑 BlockOffsets를 이용해서, 최종 PrefixSum을 얻습니다.
해당 논문은 알고리즘이나 수학, 물리에 대한 깊은 지식을 기술한것이 아니라 이 당시(2014)에 좋은 성능을 끌어내던 기술들에 대한 내용을 기술하고있다.
지배 방정식을 가장 먼저 알아야한다. 입자가 존재하고, 이 입자 i는 고유한 질량(m_i)을 가지며, 밀도(ρ_i, rho), 압력(p_i), 부피(V_i), 위치(x_i), 속도(v_i) 등의 특성을 가진다. 결국 입자 시뮬레이션에서는 이 입자를 움직여야한다. 입자가 움직임, 즉, 위치 변화율(dx_i / dt = v_i)을 구해야하고, 위치 변화율은 그 입자의 속도(v_i)와 같다. 그러므로, 새로운 위치 = 현재 위치 + 속도 * dt와 같이 게산된다. 이것을 이류(Advection)이라고한다. 이류는 움직이는 유체와 함께 물질 또는 보존량이 이송되는 것을 의미한다.
다음으로 위치를 알기위한, 속도를 얻어야하는데, 이것을 운동량 방정식(Momentum Equation)인 Navier-Stokes를 이용한다. 이것의 의미는 입자 i의 속도 변화율(dv_i / dt), 즉 가속도(acceleration)는 여러 힘들의 합에 의해 결정된다. 이 방정식은 유체 동역학의 핵심인 나비에르-스토크스 방정식을 라그랑주 관점(Lagrangian perspective)에서 본것이다.(즉, 개별 입자를 따라가면서 보는 관점) 이 가속도를 계산해서(새로운 속도 = 현재 속도 + 가속도 *dt), 다음 시간 스텝의 속도를 얻는다.
라그랑주 관점에서의 나비에 스토크스 방정식의 각 항을 분석한다.
각 항은 가속도에 해당한다. 속도를 구하기 위해서 가속도를 알아야하기때문에 각 항은 가속도의 값을 가진다는 것을 알아야한다.
압력 항 유체 내의 압력 차이로 인해 발생하는 가속도이다. ∇p_i는 압력의 gradient로, 압력이 가장 급격하게 증가하는 방향과 크기를 나타낸다. 사실 압력 값 자체는 스칼라 값이지만, 스칼라의 gradient는 벡터(방향과 크기)이다. 유체는 압력이 높은 곳에서 낮은 곳으로 밀려나려는 경향이 있고, ∇p_i는 압력이 증가하는 방향을 가리키므로, 앞에 붙은 음수 부호(-)는 힘이 압력이 감소하는 방향, 즉 고압에서 저압방향으로 작용하게만든다. 1/p_i는 밀도가 높을수록 같은 압력 차이(gradient)에 대해 덜 민감하게 반응(덜 가속)함을 의미한다.(a = F/m)
점성 항 서로 다른 속도로 움직이는 입자들 사이의 마찰력(friction), 즉 점성으로 인한 가속도를 의미한다. ν(뉴, nu)는 동점성 계수(kinematic viscosity)로, 유체의 끈적이는 정도를 나타낸다.(물이 꿀보다 ν가 훨씬 작음) ∇²v_i는 속도의 라플라시안(Laplacian)으로, 주변 입자들과의 속도 차이를 측정하는 역할을 한다. 어떤 입자의 속도가 주변보다 훨씬 빠르거나 느리면 점성항이 커져서, 속도 차이를 줄이는 방향(즉, 주변 속도와 비슷해지도록)으로 힘이 작용한다. 즉, 속도를 부드럽게 퍼뜨리는 역할을 한다. 유체의 자연스러운 감쇠 효과(흐름이 점점 느려지는 것)을 만들고, 시뮬레이션의 수치적 안정성(stability)을 향상시키는 데 도움을 준다. 텍스트에서는 실제 물의 점성 계수( ν ≈ 10⁻⁶ m²/s)보다 더 큰 값을 SPH 안정성을 위해 사용하거나, 인공 점성(Artificial Viscosity, Mon92) 또는 XSPH(Mon89) 같은 대체 기법을 사용하기도 한다. (점성 “계수”임 즉, 스칼라값)
기타 힘 항 압력과 점성 외에 입자에 작용하는 다른 모든 외부 힘 또는 체적력(Body Force)로 인한 가속도이다. F_i ^ other 는 임자 i에 작용하는 외부 힘의 합이고, 이를 질량 m_i로 나눈 것이 가속도이다. ( 뉴턴의 제 2법칙, a = F/m) 가장 대표적인 것이 중력이다. 그 외에 사용자가 인위적으로 가하는 힘(예: 파도 생성), 표면 장력(surface tension, 때로는 별도 항으로 처리), 전자기력 등이 포함될 수 있다.
논문에서의 1.2 SPH절에서는 이전 지배 방정식을 실제로 입자 시스템에서 어떻게 계산하는지 알아본다. 수학적 도구와 근사 방법을 설명해준다. 핵심은 “이산적인 입자 정보로부터 연속적인 유체 정보(값과 그 미분)를 어떻게 보간할 것인가?” 이다. 논문에서는 SPH는 임의의 위치에 유체량을 보간하고, 유한한 수의 샘플 위치(인접 입자)로의 공간 미분을 근사하는데 사용된다고한다. 임의의 위치에 유체 물리량(밀도, 압력 등)과 공간적 미분(gradient, divergence, Laplacian)을 주변에 있는 유한한 개수의 인접 입자들의 정보를 이용하여 보간(interpolate)하고, 근사(approximate)하는 방법이다.
먼저 보간의 공식은 위와 같다. 임의의 입자 i에서의 물리량 A의 값(A_i)은 그 주변에 있는 모든 인접 입자 j 들의 물리량 값(A_j)에 가중치를 곱하여 합산함으로써 근사적으로 계산할수있다는 것이고, 가중치(W_ij)는 커널 함수라고하며, 이것이 SPH의 핵심이라고할수있다. 입자 i와 j 사이의 거리에 따라 가중치를 부여하는 함수이다. 가까울수록 큰 가중치, 멀수록 작은 가중치를 부여한다.
커널함수는 위와 같다. 이름처럼 함수이므로, 괄호안에 들어가는 q에 대해서 알면된다. ||x_i - x_j||는 입자 i와 j 사이의 거리이고, h는 스무딩 길이로 커널함수의 영향력을 미칠 반경을 결정하는 값이다. 즉, 임의의 입자에서 h 값 범위 내에 있는 입자들을 “인접 입자”로 간주할것이다. 그러기에 아주 중요하다. 그래서 입자 사이의 거리(r) / 스무딩 길이(h) 이므로, 정규화된 거리이다. 보통 f(q)를 r / h로 표현하고, 1 / h^d를 곱해주는 이유는 모든 방향에 대한 커널 함수의 적분 값이 1이 되도록 보장하는 역할을 해준다. 그래서 d의 값은 공간의 차원에 해당한다(2차원이면 2, 3차원이면 3).
그리고 커널함수의 종류는 다양한데 여기서는 3D SPH에서 흔히 사용되는 큐빅 스플라인 커널을 사용한것이다. 위 공식 처럼 q값의 범위에 따라 다르게 정의되는 구간별 함수이다. q가 2이상이면 값이 0이 된다.
그다음으로 중요한 위의 공간미분(Spatial Derivatives)이 있는데, 순서대로 gradient, divergence, laplacian 이다.
이전에 소개하지않는 것은 x_ij인데 이것은 입자 i와 j 간의 벡터 차이다.
논문의 1.3. Concept of an SPH-based Fluid Solver절에서는 이제 근사 방법에 대해서 알았으니, 실제 구현하기 위한 구조와 절차를 설명한다. SPH 기반 유체 솔버의 기본 구성 요소는 다음과 같다.
이웃 탐색(Neighborhood search)는 각 입자 i에 대해서 주변 이웃 입자들을 효율적으로 찾는 과정이다. 한 입자에 대해서 다른 입자 모두를 비교하게되면, O(N^2)이 되어 너무 느리기 때문에 보통은 공간 분할 자료구조(Spatial access structure), 특히 균일 격자(Uniform Grid)를 사용해서 가속한다. 격자의 셀 크기는 보통 커널 함수의 영향 반경(ex: 2h)와 비슷하게 설정하여, 해당 셀과 인접 셀들만 탐색하면 되도록한다.
압력 계산(Pressure computation)은 지배 방정식의 압력항에서 압력 gradient(∇p_i)를 구하기위해서, 먼저 각 입자의 압력 p_i를 알아내야한다. SPH에서는 보통 밀도 ρ_i로부터 압력 p_i를 계산한다.
시간 적분(Time integration)은 계산된 힘(가속도)을 이용해 시간에 따라 입자의 속도와 위치를 업데이트한다.
압력 계산하는 방법으로 위의 상태 방정식을 이용한다(Equation of state, EOS) 밀도 ρ_i를 압력 p_i로 변환하는 함수로, SPH는 약간 압축 가능한 유체를 시뮬레이션하는 방식이기에 이 EOS는 밀도가 기준 밀도(ρ_0)에서 벗어날 때 얼마나 강하게 저항(압력 발생)하는지를 모델링한다. ρ_0는 유체의 기준 밀도이고, 유체가 안정 상태일 때 가지려는 목표 밀도이다. k는 강성 계수(Stiffness constant)로, 이 값이 클수록 밀도 변화에 대한 압력 반응이 커져서 유체가 더 비압축성처럼 행동하게된다. 하지만 k가 너무 크면 수치적 불안정성을 유발하므로, 매우 작은 시간 간격을 사용해야한다. 7은 상태 방정식의 지수인데 보통 7을 사용한다고함.
이후부터는 SPH의 여러 알고리즘에 대한 수도코드, Neighborhood 찾는 방법, EOS, Boundary Handling, Adaptive Time/Space Discretization, Multiphase Fluids, 입자 기반 인터페이스 표현 및 렌더링, 디테일 추가(거품, 스프레이, 공기방울 등), 향후 연구 방향에 대해 소개한다.
가장 먼저 Algorithm 1에 해당하는 SPH with state equation(SESPH)을 구현한다. 구현하기에 앞서 수도 코드에 대해 설명한다.
먼저 어떤 입자에 대해서 이웃 입자들을 찾는다. 여기서는 앞서 언급한 Uniform Grid를 사용한다.
어떤 입자의 밀도를 계산하는데, 이웃 입자들의 질량과 커널 가중치의 기여도의 합산으로 계산하고, 그 밀도를 이용해서 아래식으로 압력을 구한다.
힘을 계산한다. 압력 gradient를 구해서 압력힘을 구하고, 속도 laplacian으로 점성힘을 구함. 이후에 파도 외부힘 구현할때, 이 항을 더 사용한다. 기본적으로 중력(g)를 설정한다. 이제 구한 힘을 다 더해 입자 i의 총 힘을 계산한다.
마지막으로 시간 적분단계를 진행한다. 맨처음에 설명했듯이, 새로운 위치 = 현재 위치 + 새로운 속도 * dt이고, 새로운 속도 = 현재 속도 + 가속도 * dt이다. 최종적으로 구하고싶은 새로운 위치를 구해야하고, 현재는 a = F / m로 가속도를 구할수있다. 구한 가속도로 새로운 속도를 구하고, 새로운 속도로 새로운 위치를 구하면된다.
당장의 구현 방향은 다음과 같다. 앞으로 진행될 과정에 대해서는 최대한 작은 단위로 나눠서 세분화하도록한다.
2D 파티클 시스템 구현 및 생명 주기를 가진 입자를 주기적으로 생성하고, 간단한 경계 충돌 테스트
Compute Shader에서 구현
2D 상에서 SESPH를 구현 및 입자 생성, 소멸, 간단한 경계 충돌 테스트
입자를 특정 위치에서 특정 속도로 생성. 인접 인자 확인은 Uniform Grid로 구현. 여전히 간단한 경계 충돌
