Smoke Simulation

홍정모님의 그래픽스 새싹코스 강의를 듣고 정리한 내용입니다.

Smoke Simulation

\[\frac{\partial u}{\Delta t} = - \left ( u \cdot \bigtriangledown \right )u - \frac{1}{\rho} p + \nu \bigtriangledown ^{2} u\]

void StableFluids::Update() {
    // ...

    Sourcing(context);
    Advection(context);
    Diffuse(context);
    Projection(context);
}

• Sourcing : External Forces 계산 (예제에서는 마우스의 속도, 그리드 셀의 컬러)
• Advection : Velocity Field 업데이트 (유체가 스스로를 이동시키며 내부 상태 변화)
• Diffuse : viscosity(점성) 효과를 적용해 값이 퍼짐
• Projection : pressure(압력) 계산 (유체가 비압축성 성질을 갖도록 압력 보정)

Stable Fluids예제를 간단히 복습하면 Navier-Stokes 방정식을 위 방식으로 풀었고, pixel shader에서는 단순히 해당 픽셀의 density값을 가져와 사용했었다. 여기서 조금만 더 응용하면 smoke simulation을 구현할 수 있다. 바로 Stable Fluid를 이용해 유체의 흐름을 계산하고, Cloud Rendering에서 사용한 lighting 계산을 이용해 color를 결정하는 방법이다. 단순 Density를 이용하는 것이 아닌 Beer’s law를 이용해 빛을 감쇠시키고, HenyeyGreenstein 함수를 적용해 color를 지정하면 된다.

void RealtimeSmoke::Update() {
    // ...

    DownSample();
    Sourcing();
    Projection();
    DiffUpSample();
    Advection();
}

Smoke Simulation은 위와 같은 단계로 수행된다. Stable Fluids예제와 달리 smoke simulation은 3차원에서 연산을 수행해야 하기 때문에 비용이 매우 비싸진다. 따라서 3D 텍스처를 DownSampling한 공간에서 연산을 수행하고 다시 UpSampling한 결과를 이용해 최종적으로 simulation을 구현한다.

No Object	Add Object	Boundary Conditions

유체 중간에 object를 두게되면 와류가 발생해 훨씬 자연스러운 연기 효과가 표현된다. 하지만 이를 구현하기 위해서는 추가적인 세팅(boundary condition)이 필요하다.

• Dirichlet Condition : 전체 경계 부분에서 사용
  • 속도 값을 상수로 고정 (압력을 0으로 세팅)
• Neumann Condition : 내부 오브젝트와의 경계에서 사용
  • 유체가 경계를 통과하지 않는 상황 (flux = 0)
  • 속도의 변화율을 상수(0)로 고정
  • 오브젝트 경계 외부와 내부의 압력 크기를 같게 세팅(부호만 반대로)

DownSampling

void DownSample()
{
  // Vorticity Confinement...

  // DownSampling...

  // Backup
}

• Vorticity Confinement : 저해상도 격자를 이용해 유체의 흐름을 계산하는 경우 국소적인 소용돌이 현상이 소멸되는 경향이 있다.
  • 이를 방지하기 위해 인위적으로 소용돌이 힘을 추가해주는 작업. (velocity field의 curl이용)
• DownSampling : DownSampling. 단순 interpolation.
• Backup : 뒤에 진행되는 Upsampling 단계에서 사용하기 위해 업데이트 전 density, velocity 데이터 저장
  • UpSampling을 진행할 때 저해상도 데이터의 변화율을 기반으로 interpolation하면 훨씬 부드러운 결과가 연출됨.

Sourcing

void main(uint3 dtID : SV_DispatchThreadID)
{
    uint width, height, depth;
    bc.GetDimensions(width, height, depth);
    
    bc[dtID] = 0;
    
    if (dtID.x == 0 || dtID.y == 0 || dtID.z == 0
        || dtID.x == width - 1 || dtID.y == height - 1 || dtID.z == depth - 1)
    {
        bc[dtID] = -1;
        density[dtID.xyz] = 0.0;
    }

    // Source
    float3 center = float3(0.02, 0.3, 0.5) / dxBase;
    int radius = 0.2 * height;

    float dist = length(float3(dtID.xyz) - center) / radius;
    
    if (dist < 1.0)
    {
        velocity[dtID.xyz] = float4(32 * sourceStrength, 0, 0, 0) / 64.0 * float(width); // scale up velocity
        density[dtID.xyz] = max(smootherstep(1.0 - dist), density[dtID.xyz]);
    }

    // Object
    center = float3(0.15, 0.3, 0.5) / dxBase;
    radius = 0.1 * height;
    
    dist = length(float3(dtID.xyz) - center) / radius;
    
    if (dist < 1.0)
    {
        velocity[dtID.xyz] = float4(0, 0, 0, 0) / 64.0 * width;
        bc[dtID.xyz] = -2; // Neumann
    }
}

• 가장자리 1픽셀은 Dirichlet Boundary Condition(-1) 지정
• Smoke가 시작되는 Source 위치, 초기 상태 지정
• 와류를 발생시키는 Object 세팅, 해당 부분 Neumann Boundary Condition(-2) 지정
  • 이번 예제에서 object는 움직이지 않기 때문에 velocity = 0.

Projection

void Projection()
{
  // Compute Divergence...

  // Jacobi Iteration...

  // Apply Pressure...
}

• Compute Divergence : Boundary Condition을 고려해 각 셀의 Divergence 계산
• Jacobi Iteration : Divergence를 이용해 Pressure 업데이트 (안정성을 위해 jacobi iteration)
• Apply Pressure : 업데이트 된 Pressure로 Velocity 업데이트

Divergence

static int3 offset[6] =
{
    int3(1, 0, 0), // right
    int3(-1, 0, 0), // left
    int3(0, 1, 0), // up
    int3(0, -1, 0), // down
    int3(0, 0, 1), // back
    int3(0, 0, -1) // front
};

void main(uint3 dtID : SV_DispatchThreadID)
{
    if (bc[dtID.xyz] >= 0)
    {
        float div = 0.0;

        [unroll]
        for (int i = 0; i < 6; i++)
        {
            if (bc[dtID.xyz + offset[i]] == -1) // Dirichlet
            {
                div += dot(velocity[dtID.xyz].xyz, float3(offset[i]));
            }
            if (bc[dtID.xyz + offset[i]] == -2) // Neumann
            {
                div += dot(2 * velocity[dtID.xyz + offset[i]].xyz - velocity[dtID.xyz].xyz, float3(offset[i]));
            }
            else
            {
                div += dot(velocity[dtID.xyz + offset[i]].xyz, float3(offset[i]));
            }
        }

        divergence[dtID.xyz] = 0.5 * div;
    }
}

• 이웃 셀이 Dirichlet Condition인 경우
  • 이웃 셀의 압력이 0이기 때문에 현재 셀의 성분이 그대로 발산
• 이웃 셀이 Neumann Condition인 경우
  • 경계 부분에서 속도의 변화율을 0으로 고정
  • \[\frac{v_{bdr} - v_{nxt}}{\Delta x} = \frac{v_{nxt} - v_{cur}}{\Delta x}\]
  • \[v_{bdr} = 2v_{nxt} - v_{cur}\]
  • v_bdr : 경계 부분의 속도, v_cur : 해당 셀의 속도, v_nxt : 이웃 셀의 속도

Apply Pressure

void main(uint3 dtID : SV_DispatchThreadID)
{
    if (bc[dtID.xyz] >= 0)
    {
        float p[6];

        [unroll]
        for (int i = 0; i < 6; i++)
        {
            if (bc[dtID.xyz + offset[i]] == -1) // Dirichlet
            {
                p[i] = -pressure[dtID.xyz];
            }
            else if (bc[dtID.xyz + offset[i]] == -2) // Neumann
            {
                p[i] = pressure[dtID.xyz];
            }
            else
                p[i] = pressure[dtID.xyz + offset[i]];
        }

        velocity[dtID.xyz] -= 0.5 * float4(p[0] - p[1], p[2] - p[3], p[4] - p[5], 0);
    }
}

• 이웃 셀이 Dirichlet Condition인 경우
  • 유체의 흐름이 유지되려면 압력의 평균이 0이어야 함
  • P[i] = -pressure[cur]
• 이웃 셀이 Neumann Condition인 경우
  • 유체가 경계를 통과하지 않으려면 압력이 동일해야 함
  • P[i] = pressure[cur]

Diffuse & UpSampling

void main(uint3 dtID : SV_DispatchThreadID)
{
    float3 uvw = (dtID + 0.5) * dxUp;
    
    float coeff = 0.99;
    
    float4 velOld = velocityOld.SampleLevel(linearClampSS, uvw, 0);
    float4 velNew = velocityNew.SampleLevel(linearClampSS, uvw, 0);
    velocityUp[dtID] = lerp(velNew, velocityUp[dtID] + velNew - velOld, coeff);
    
    float denOld = densityOld.SampleLevel(linearClampSS, uvw, 0);
    float denNew = densityNew.SampleLevel(linearClampSS, uvw, 0);
    densityUp[dtID] = lerp(denNew, densityUp[dtID] + denNew - denOld, coeff);
}

• velocityUp[dtID] = lerp(velNew, velocityUp[dtID] + velNew - velOld, coeff)
  • velNew : 저해상도에서 업데이트 된 속도
  • velOld : 저해상도에서 업데이트 되기 전 속도
  • velocityUp : 고해상도에서 업데이트 되기 전 속도
• UpSampling을 진행할 때 저해상도 데이터의 변화율을 기반으로 interpolation하면 훨씬 부드러운 결과 연출
• 유체 시뮬레이션 분야에서 Upsampling할 때 많이 사용하는 기법이라고 함

Advection

void main(uint3 dtID : SV_DispatchThreadID) 
{
    uint width, height, depth;
    velocity.GetDimensions(width, height, depth);
    float3 dx = float3(1.0 / width, 1.0 / height, 1.0 / depth);
    float3 uvw = (dtID.xyz + 0.5) * dx;
    
    float3 vel = velocityTemp[dtID.xyz].xyz * dx;
    float3 uvwBack = uvw - vel * dt * upScale;
    
    density[dtID.xyz] = densityTemp.SampleLevel(linearClampSS, uvwBack, 0) * 0.999;
    velocity[dtID.xyz] = velocityTemp.SampleLevel(linearClampSS, uvwBack, 0) * 0.999;

}

• Advection의 경우 Stable Fluids 예제와 동일
• density, velocity를 업데이트 할 때 0.999를 곱해주어 점점 약해지도록 구현
  • 수치적으로 안정화 시키는 역할도 있다고 함

Results

맨 위로 이동하기

Graphics 카테고리 내 다른 글 보러가기