코딩 어드벤처: 연기 시뮬레이션

비디오: 코딩 어드벤처: 연기 시뮬레이션
채널: Sebastian Lague
길이: 41분 39초
출처: 자막 (수동, en-GB)

안녕하세요 여러분, 코딩 어드벤처(coding adventures)의 새로운 에피소드에 오신 것을 환영합니다! 오늘은 연기를 시뮬레이션(simulate)해보고자 합니다. 연기는 공기 중에 떠 있는 수많은 작은 입자(particles)라고 생각할 수 있으며, 공기는 유체(fluid)이기 때문에 유체 시뮬레이션(fluid simulation)으로 연기를 모델링하는 것이 타당합니다. 이것은 제가 얼마 전에 만든 입자 기반 유체 시뮬레이션(particle-based fluid simulation)인데, 오늘 우리의 목적에 맞게 조정할 수 있을 것이라 확신합니다. 예를 들어 다른 프로젝트에서는 오디오 신호를 포착하기 위해 공기 중의 조잡한 압력파(pressure waves)를 만드는 데 이것을 사용하기도 했습니다. 하지만 이번에는 — 새로운 기술을 배우기 위해 — 유체를 개별 입자가 아닌 그리드(grid) 상의 값으로 모델링하는 완전히 다른 접근 방식을 시도해보고 싶습니다. 이것이 어떻게 작동하는지 거의 전혀 모르지만, 우리를 안내해 줄 유용해 보이는 강의 노트(course notes)를 찾았습니다. 우선, 비압축성 유체(incompressible fluid)의 흐름을 설명하는 다소 단순화된 형태의 나비에-스토크스 방정식(navier-stokes equations)에 대해 잠시 생각해 봅시다. 실제로 모든 유체는 아주 조금이라도 압축될 수 있지만, 그렇지 않은 척함으로써 이런 방정식을 보지 않아도 된다면 저는 전적으로 찬성입니다. 여기서 우리가 가진 식은 유체 내부에서 시뮬레이션하는 모든 미세한 부피가 압력 구배(pressure gradient)를 따라 가속되어야 함을 말해줍니다. 즉, 균형을 맞추기 위해 압력이 높은 지역에서 낮은 지역으로 이동한다는 의미입니다. 하지만 관성(inertia) 때문에 유체의 밀도(density)가 높을수록 가속하기가 더 어려워지며, 이를 반영하기 위해 여기에 밀도로 나누는 과정이 포함되어 있습니다.

유체는 중력이나 바람, 혹은 그 외의 다른 것과 같이 존재할 수 있는 모든 외력 (external forces)에 의해서도 가속될 것입니다. 그리고 마지막으로 점성 (viscosity) 항이 있는데, 이는 일종의 내부 마찰력 (internal friction force) 때문에 인접한 영역의 속도들이 시간이 지남에 따라 서로 더 유사해지도록 만듭니다. 하지만 오늘은 이 항의 승수를 0이라고 가정하고, 계산을 조금 더 쉽게 만들기 위해 당분간은 외력이 없다고 가정해 봅시다. 좋습니다, 그러면 우리는 유체의 발산 (divergence)이 0이어야 함을 알려주는 아주 작은 두 번째 방정식도 가지고 있습니다. 이를 생각해 보기 위해, 우리가 만들 더 큰 그리드 (grid) 내의 단일 유체 셀 (fluid cell)을 상상해 봅시다. 예를 들어, 오른쪽 가장자리를 통해 유체가 흘러나오고 있다고 가정해 보겠습니다. 이 경우 유체는 셀로부터 발산 (diverging)하고 있는데, 방금 본 방정식에 따르면 이는 매우 불가능한(illegal) 일입니다. 따라서 상태를 일정하게 유지하려면, 흘러나가는 유체만큼 동일한 양의 유체가 흘러 들어와야 합니다. 이 유입 (inflow)이 완전히 왼쪽에서만 일어나는지, 아니면 다른 모든 가장자리를 통해 조금씩 흘러 들어오는지 여부는 중요하지 않습니다. 최종적으로 합계가 0이 되기만 한다면 모든 것이 순조롭게 진행될 것입니다. 우리가 볼 수 있듯이, 각 가장자리에서의 정확한 유체 속도(velocity)를 아는 것 — 또는 적어도 왼쪽과 오른쪽 가장자리에서의 수평 성분(horizontal component)과 위아래 가장자리에서의 수직 성분(vertical component)을 아는 것은 유체 흐름이 잘 균형을 이루고 있는지 확인하는 데 매우 편리합니다.

그렇다면 코드에서도 그런 방식으로 저장해 봅시다. 저는 원하는 셀(cell)의 개수와 월드에서의 물리적 크기를 입력받는 작은 유체 그리드(fluid grid) 클래스를 만들었습니다. 이 클래스는 수평 및 수직 속도(velocity) 배열을 초기화하는데, 이 배열들의 크기는 약간 다르게 설정되어 조금 어색할 수 있습니다. 예를 들어 3x2 그리드를 상상해 본다면, 이는 일종의 4x2 좌우 가장자리(edge) 그리드와 3x3 상하 가장자리 그리드를 포함하기 때문입니다. 어쨌든, 우리는 셀을 가로 방향과 세로 방향으로 속도가 얼마나 빠르게 변하는지를 살펴보고 이를 모두 더함으로써 어떤 셀에 대해서든 발산(divergence) 값을 계산할 수 있습니다. 이를 통해 유체가 주로 흘러나가면 음수(-), 주로 흘러 들어오면 양수(+), 그리고 완벽하게 균형을 이루면 이상적으로 0인 전체적인 측정값을 얻을 수 있습니다. 또한 저는 이를 위한 작은 시각화(visualization) 도구를 만들기 시작하고 싶습니다. 그렇지 않으면 코드가 무엇을 하고 있는지에 대한 저의 정신적 모델(mental model)이 실제와 매우 빠르게 어긋나는 것을 발견하곤 하기 때문입니다. 그래서 그리드를 입력받아 전체 크기 등을 파악하는 작은 스크립트를 만들었습니다. 이를 통해 어떤 셀이나 가장자리의 위치도 쉽게 찾아볼 수 있습니다. 그런 다음 단순히 루프를 돌며 각 셀을 나타내는 작은 상자와 각 속도 성분(velocity component)을 나타내는 작은 화살표를 그릴 수 있습니다. 자, 그럼 아마도 아주 작은 5x3 그리드로 이것을 테스트해 봅시다. 그리고 여기서 시각화 파라미터(parameter)를 약간 조정하겠습니다. 현재는 속도가 모두 0이라서 화살표가 나타나지 않지만, 어떤 모습인지 볼 수 있도록 속도를 무작위로 설정하는 키보드 단축키를 만들어 두었습니다.

이제 이 모든 셀(cell)들의 발산(divergence)을 확인해 보고 싶으므로, 다시 해당 코드로 돌아가서 각 셀에 대해 우리의 작은 발산 함수를 실행하고, 값이 음수인지 양수인지에 따라 색상을 빨간색 또는 파란색으로 혼합해 보겠습니다. 좋습니다, 한번 살펴봅시다. 다시 속도(velocities)를 무작위로 설정하겠습니다... 이제 예를 들어, 많은 유체가 여기 있는 이 셀로 수렴(converging)하고 압축(compressing)되고 있는 것을 볼 수 있으며, 그래서 밝은 빨간색으로 표시됩니다. 그리고 바로 옆 칸에는 유체가 발산(diverging)하며 멀어지는 반대의 상황이 밝은 파란색으로 나타납니다. 이 속도 성분(velocity components)들을 이리저리 드래그하여 수동으로 이를 수정해 볼 수도 있습니다. 비록 각 조정이 우리가 수정하려는 셀 하나만이 아니라 두 개의 셀에 영향을 미친다는 사실 때문에 조금 까다롭기는 하지만 말입니다. 모든 속도를 0으로 설정하는 것이 발산을 0으로 보장하는 쉽고 좋은 방법이 될 수도 있겠지만, 당연히 우리가 이 화살표들을 마음대로 바꿀 수는 없습니다. 유체가 가속되는 방식에 대한 첫 번째 방정식을 따라야 합니다. 그러니 하나의 셀에 집중해 봅시다. 이를 중앙 셀 $C$라고 부르고, 왼쪽, 오른쪽, 위, 아래에 이웃(neighbours)이 있다고 가정하겠습니다. 그리고 이들 각각은 자신만의 압력(pressure) 값을 가질 것입니다. 또한 중앙 셀의 가장자리 속도(edge velocities)를 $u_l, u_r, u_t, u_b$라고 표시해 보겠습니다. 참고로 우리가 속도에 $u$를 사용하는 이유는, 이제 아주 작은 타임스텝(timestep) 후 미래의 유체 새로운 속도에 대한 예측값으로 $v$를 사용할 것이기 때문입니다. 예를 들어, 오른쪽 가장자리의 새로운 속도는 현재 값에 그곳에서 유체가 겪는 가속도를 더하고, 여기에 아주 작은 시간 단계를 곱한 값과 같다고 말할 수 있을 것입니다.

하지만 우리는 가속도 (acceleration)가 밀도 (density)에 대한 음의 압력 구배 (negative pressure gradient)와 같아야 한다는 것을 알고 있으므로, 이를 대입할 수 있습니다. 그런데 압력 구배 (pressure gradient)란 공간에 따라 압력이 얼마나 빠르게 증가하는지를 측정하는 값일 뿐입니다. 따라서 현재 우리가 관심을 두고 있는 격자 (grid) 상의 수평 성분 (horizontal component)은 단순히 오른쪽 셀과 중앙 셀 사이의 압력 차이를 그 사이의 거리로 나눈 값이 될 것입니다. 참고로, 속도 (velocities)를 가장자리를 따라 이 약간은 생소한 방식으로 저장하는 것이 매우 편리한 이유가 바로 이것입니다. 압력 구배 (pressure gradient)를 계산할 때 양옆의 압력 값과 아주 깔끔하게 맞아떨어지기 때문입니다. 어쨌든, 이제 우리는 이 방정식을 얻게 됩니다. 저는 여기서 격자 전체에서 일정한 값들을 초록색으로 표시했습니다. 물론 시간 단계 (time step)는 어디서나 동일하며, 셀의 크기 (size of the cell)와 중요한 요소인 밀도 (density) 또한 마찬가지입니다. 우리는 발산 제약 조건 (divergence constraint) 때문에 이것이 사실임을 알고 있습니다. 왜냐하면 공간의 모든 영역에서 유입되는 양만큼의 유체가 유출되어야 하므로, 각 위치의 질량 (mass)은 변할 수 없기 때문입니다. 좋습니다, 이 작은 방정식을 아주 조금 정리하기 위해 상수 (constants)들을 문자 k로 묶어서 다음과 같이 나타내겠습니다. 그런 다음 왼쪽, 위, 아래에 대해서도 정확히 동일한 과정을 반복하여, 우리 중앙 셀의 모든 가장자리에 대한 새로운 속도 (velocities) 방정식을 도출할 수 있습니다. 불행히도, 압력 값 (pressure values)이 무엇인지 모른다면 이것은 다소 쓸모가 없습니다. 따라서 우리의 목표는 실제로 중앙 셀의 압력을 구하는 것이며, 그렇게 되면 모든 곳의 압력을 계산할 수 있게 될 것입니다.

이 4개의 작은 방정식들을 하나의 큰 방정식으로 결합하는 것부터 시작해 봅시다. 그러면 실제로 4개의 가장자리 속도 (edge velocities)가 포함된 편리한 식을 얻게 되는데, 이것은 다시 한번 발산 제약 조건 (divergence constraint)입니다. 왜냐하면 이 계산 방식은 각 축을 따라 셀 전체에서 속도가 어떻게 변하는지를 살펴보고, 그것이 모두 합쳐졌을 때 어떤 값이 되는지를 확인하는 것이기 때문입니다. 따라서 우리는 새로운 속도들의 발산 (divergence)이 0과 같아야 한다고 쓸 수 있으며, 그다음 각 속도에 대한 실제 계산값들을 대입하면 됩니다. 이 과정에서 식이 꽤 복잡해지긴 하지만, 이제 중앙 압력 (central pressure)을 구하는 것은 단순한 대수학 (algebra)의 문제입니다. 자, 결국 우리가 발견한 것은 중앙 셀의 압력이 이웃 셀들의 평균 압력에서, 셀의 크기와 밀도를 곱한 값을 빼고, 이를 속도 변화량 (velocity deltas)의 합을 4배의 타임스텝 (timestep)으로 나눈 값으로 스케일링한 것과 같다는 사실입니다. 이제 우리는 이 방정식을 실제로 풀기 위한 코드를 작성해야 합니다. 각 셀의 압력이 주변 셀들의 미지수 압력을 포함하고 있기 때문에 이 과정은 약간 까다롭습니다. 노트에서는 "수정된 불완전 촐레스키 공액 기울기법 (modified incomplete Cholesky conjugate gradient)"이라는 알고리즘을 추천하고 있는데, 이름만큼 무서운 것은 아니라고 합니다. 하지만 저는 오늘 일단 익숙한 방식인 Gauss-Seidel 반복법 (Gauss-Seidel iteration)으로 시작하고 싶습니다. 저는 빨리 유체가 흐르는 것을 보고 싶어 조바심이 나거든요. 더 효율적인 함수들은 나중에 언제든지 만져볼 수 있으니까요. 그래서 코드에서는 압력 값들의 그리드 (grid)를 정의하고, 타임스텝과 밀도에 대한 상수들을 설정한 다음, 우리가 찾아낸 방정식을 사용하여 특정 셀의 압력을 업데이트하는 작은 함수를 만들었습니다.

이웃한 압력 (neighbouring pressures) 값들에 대해서는, 현재 그리드 (grid)에 저장되어 있는 값을 그대로 사용합니다. 따라서 시작 단계에서는 모두 0일 것입니다. 예를 들어, 여기 있는 이 셀의 압력을 계산한다고 가정해 봅시다. 현재 이 셀은 12.5라는 값을 가지고 있지만, 만약 우리가 바로 옆에 있는 셀의 압력 값을 계산한다면, 그 값은 이제 구식이 되어버릴 것입니다. 이웃 셀에 대한 최신 정보를 사용하지 않았기 때문입니다. 이를 해결하기 위해, 우리는 아주 놀라운 트릭을 사용할 수 있습니다... 바로 다시 계산하는 것이죠. 하지만 물론 이제는 이웃 셀의 값이 구식이 되었으므로 그것 역시 수정해야 하며, 이런 식으로 계속 반복됩니다. 하지만 꽤 빠르게 값들이 안정적인 상태 (stable state)로 수렴하는 것을 볼 수 있습니다. 그러니 단순히 전체 그리드를 훑으며 진행하면서 각 압력 값을 업데이트하는 작은 함수를 만들어 봅시다. 그런 다음, 우리가 희망하는 결과에 도달할 때까지 그 함수를 여러 번 실행하면 됩니다...