이승민의 개발 블로그

SPH 탄생 배경

SPH의 탄생 배경은 SPH(Smoothed Particle Hydrodynamics)의 시초로 꼽히는 논문( Gingold & Monaghan 1977)에 잘 나와있습니다.

SPH는 1977년경 천문학 및 우주론 분야에서 처음 제안되었습니다.

당시 사용되던 오일러 방식(Eulerian)의 격자 기반 수치 해석 기법은, 고정된 격자(공간의 한 지점)를 기준으로 유체의 물리량 변화를 관찰하는 방식인데, 3차원 문제의 경우 격자 밀도가 급격히 증가하여 계산량이 기하급수적으로 늘어나 시뮬레이션하기 어려웠습니다.

이에 반해, 라그랑지안 방식(Lagrangian)은 유체를 구성하는 개별 입자 혹은 물질 덩어리를 따라가며 그 변화를 추적하는 방식으로, 자유 표면이나 큰 변형을 자연스럽게 처리할 수 있으면서 계산량을 줄일 수 있다는 장점이 있었습니다.

따라서, SPH는 계산 효율성을 크게 개선하고 복잡한 천체 물리 문제를 다루기 위해 이 라그랑지안 방식을 도입한 것입니다

그러면 Smoothing 함수가 무엇인가?

여기서의 스무딩 함수는 SPH에서는 각 입자가 고유의 질량, 밀도, 속도 등의 값을 갖고 있지만, 해당 값들이 하나의 점에 국한되지 않고, "주변 영역에 걸쳐 영향을 미치도록" 하기위해서 Smoothing 함수 정의하여 사용한다.

이 Smoothing 함수를 Kernel(커널) 함수라고 부르고, 이 함수는 보통 입자 사이의 거리와 Smoothing(평활화) 반경(h)에 따라서 가중치를 부여하고, 하나의 입자에 의한 물리량을 주변으로 "환산"시켜 마치 블러 효과를 주듯이, 마치 하나의 점의 영향력을 표현하자면, 하나의 점이 주변으로 스무딩되듯 부드럽게 연결한다.

위 사진은 커널 함수의 원리에 대한 사진이다. 한 입자에서 이웃 입자 사이의 거리와 평활화 반경(h) 매개변수로 kernel 입자에 전달해주게 되면, 얼마만큼의 가중치가 발생하는지를 알수있다.

여러 SPH Kernel Function 소개 : https://interactivecomputergraphics.github.io/physics-simulation/examples/sph_kernel.html

그렇다면 이 커널 함수의 결과를 어디에 사용하는가?

직전에 말했듯, 커널 함수를 통해서 한 입자와 다른 입자 사이의 거리와 평활화 반경을 통해서 얼마만큼의 영향력이 발생하는지를 알수있다고 했다.

그렇다면 평활화 반경을 통해 한 입자에서 평활화 반경 내에 있는 모든 입자들까지의 영향력을 알아낼수있다.

즉, 한 입자에서 평활화 반경 내에 존재하는 입자들과의 영향력의 총합을 통해, 해당 입자의 물리적 특성을 근사할수있다.

그렇게 해서 근사해낸 영향력으로 한 입자의 밀도를 구하고, 밀도로 압력을 구하고, 압력으로 힘을 구해서, 해당 입자가 다음 타임 스텝에 어느 방향으로 움직여야하는지를 근사해, 다음 위치를 계산하는 것이다.

그럼 먼저 사용할 커널 함수 방정식을 세운다.

평활화 반경(h)와 입자 사이의 거리(x)의 차에 세제곱(음수가 될수있으니 미리 0이상으로 설정)을 하면, 왼쪽 사진의 커널 그래프 모습과 비슷하게 오른쪽과 같이 나온다.

그런데 이 식을 그대로 쓰면안되고, 위 식을 통해서 나오는 값은 영향력 분포 형태를 정의하는것이다.

이걸 단위부피당 값으로 표현해서, 전체 부피를 1로 만들어야하기때문에 해당 위치에서의 부피로 나눠줘야한다.

그래서 위 적분식의 결과를 통해서 오른쪽의 식으로 부피를 구한다음 이걸로 나눠준 결과를 최종 영향력으로 사용하면된다.

이렇게 커널 함수를 설정하고, 이 결과 값을 가중치로 사용한다.

위는 입자를 64 x 64로 배치한 모습이다.

그렇다면 일반적으로 이 상황에서 "한 입자"에서 "평활화 반경 내의 밀도"는 어떻게 측정되어야할까?

모든 입자들이 동일한 간격으로 배치되어있으니, 반경 내에 입자들의 간격이 동일하다면, 입자 배치의 중심에서 평활화 반경을 증가시키거나 줄였을때, 측정되는 밀도는 항상 동일할것이다.

왼쪽은 평활화 반경이 1일때이고, 오른쪽은 3일때 측정한것인데, 수치가 비슷하지만 오차가 발생하는 이유는 부동소수점 연산때문에 그렇다.

위처럼 특정 평활화 반경에서 그 주변 입자의 몇몇 입자들과 밀도를 비교해서 거의 같으면된다

밀도 구하기

어떤 입자의 물리량을 다른 입자들과의 연관성을 통해서 알고싶을때 위 보간식을 사용해서 계산할수있다.

위 식은 위치 x에서의 물리량 A(x) 는 주변 입자 i 들의 물리량 A_i에 각 입자의 크기(부피)와 위치 x에 대한 영향력(커널 함수)를 곱한 값들을 모두 더해서 얻어진다. 라는 의미이다.

그랬을때, A_i가 밀도라고 하면, 이미 밀도로 나눠주기때문에 상쇄되어 사라지게된다.

그럼 최종적으로 밀도를 구하는 식은 위와 같다.

출처

https://www.researchgate.net/figure/The-principle-of-the-SPH-kernel-function-modified-from-90_fig6_334481027

https://www.youtube.com/watch?v=rSKMYc1CQHE&t=287s&ab_channel=SebastianLague

https://interactivecomputergraphics.github.io/physics-simulation/examples/sph_kernel.html

https://matthias-research.github.io/pages/publications/sca03.pdf

동작 과정

1. 파티클 해시 계산

파티클의 위치값과 셀크기로 해당 파티클의 해시 값 계산

 Particle p = Particles[index];

    // 파티클 위치 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;

입력 : 파티클 버퍼

struct Particle {
		XMFLOAT3 position = XMFLOAT3(0.0f, 0.0f, 0.0f);
		float radius = 0.0f;
		XMFLOAT3 velocity = XMFLOAT3(0.0f, 0.0f, 0.0f);
		float life = -1.0f;
		XMFLOAT3 color = XMFLOAT3(1.0f, 1.0f, 1.0f);
		float density = 0.0f;
		XMFLOAT3 force = XMFLOAT3(0.0f, 0.0f, 0.0f);
		float pressure = 0.0f;
	};

출력 : 파티클 해시 버퍼

struct ParticleHash
	{
		UINT particleID = 0; // 원래 파티클 인덱스
		UINT hashValue = 0;  // 계산된 해시 값
		UINT flag = 0; // 그룹 시작 플래그
	};

2. 해시 값 기준 정렬

BitonicSort로 해시 값을 기준으로 오름차순으로 정렬 - Bitonic Sort 사용.

입력: ParticleHash

출력 : 정렬된 ParticleHash

3. 플래그 생성

정렬된 ParticleHash에서 i와 i - 1에서의 hashValue가 다르면 i 위치의 flag를 1로 설정

입력 : 정렬된 ParticleHash

출력 : 플래그 설정된 정렬된 ParticleHash

4. 플래그 스캔으로 그룹 생성

4-1. 지역스캔

스레드 그룹들끼리 공유 메모리를 이용해서 지역적으로 자기 위치 이전까지 몇개의 그룹이 존재하는지 측정

그룹 사이즈 - 1 위치의 값은 그 위치 이전에 지역적으로 존재하는 최대 그룹 수

그룹 사이즈 - 1 위치"까지"의 최대 그룹수는 그룹 사이즈 - 1 위치의 flag의 값이 1이면, 이전 최대 그룹 수에 1을 더해주면됨

4-2. 지역 스캔 결과로 나온 BlockSum들의 지역 스캔

4-1의 결과로 나온 그룹당 최대 그룹 수들의 모음을 가지고 4-1의 동작을 한번더 실행하면, 4-1의 각 그룹에 해당하는 그룹 Offset 값들( BlockSum )이 나오게 되고, 4-1의 지역 스캔결과와 4-2의 BlockSum을 연산해주면, 전체 정렬된 ParticleHash 의 각 객체에 해당하는 그룹들이 설정됨

입력 : 정렬된 ParticleHash

출력 : 전체 Scan

5. CompactCell 및 CellMap 설정

struct CompactCell
	{
		UINT hashValue = 0;
		UINT startIndex = 0;
		UINT endIndex = 0;
	};

4의 결과를 통해서, 각 그룹이 정렬된 ParticleHash에서 어느 범위에 해당하는지를 설정

스캔 결과의 현재 index 위치의 값이 그룹 값이 된다

5-1. 현재 index의 정렬된 ParticleHash의 flag가 1이라면, 현재 그룹의 startIndex에 해당함

5-2. 현재 그룹이 최대 그룹 개수 - 1이면, ParticleHash의 최대개수- 1이 endIndex가 되고, 그 외에는 CompactCell에서 현재 그룹 + 1인 다음 그룹의 startIndex가 endIndex가 된다.

5-3. 이후 uint로 이루어진 CellMap을 설정할건데, 현재 그룹의 해시값 위치의 CellMap값이 현재 그룹 값으로 설정된다.

이전에 CellMap의 각 값은 -1로 초기화된다.

6. 파티클 시스템

https://minsdevblog.tistory.com/28

참조

https://cg.informatik.uni-freiburg.de/publications/2014_EG_SPH_STAR.pdf

https://github.com/steampower33/LSMEngine

입자의 반지름(r) : 1.0

커널 크기(h) : 반지름 * 4

질량(m) : 0.5

압력 계수(pressureCoeff) : 5.0

기준 밀도(density0) : 1.0

점성 계수(viscosityCoeff) : 0.5

으로 설정해서 시뮬레이션을 진행했습니다.

위 과정을 진행하여 밀도, 압력을 구하고, 힘을 구해 속도와 위치를 업데이트합니다.

위 방정식으로 이전에 구한 밀도를 통해 압력을 구합니다.

그리고 위 방정식들을 통해서 각 입자의 힘을 구합니다.

커널 함수의 경우 위 식을 사용하는데, 현재 h^d를 f(q)에 나눠주게 되면 시뮬레이션이 불안정해져서 그 원인을 찾고있습니다.

입자개수 8192개를 렌더한 결과입니다.

참조

https://cg.informatik.uni-freiburg.de/publications/2014_EG_SPH_STAR.pdf