Stable Fluids

Date: 2024.10.05 Updated: 2024.10.05

카테고리: Graphics

태그: DirectX

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

Other References:
Stable Fluids Paper
Ten Minute Physics
Projection method (fluid dynamics)
Stable Fluids with three.js

Chorin’s Projection Method

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

Stable Fluid 역시 SPH와 동일하게 Incompressible한 Navier-Stokes 방정식을 푸는 컨셉이다. 차이점이라 한다면 SPH는 Particle의 속도와 압력을 직접 계산해 수정하는 방식이었지만, Stable Fluid에서는 공간의 그리드 셀 마다 속도와 압력을 가지고 있으며 해당 공간의 Vector Field를 업데이트하는 방식이다. 또한 방정식을 푸는 방법도 다른데 Stable Fluid에서는 Chorin이라는 수학자가 제시한 Projection Method라는 방식을 이용한다. (수학과 만세..)

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

n번째 스텝에서 바로 n+1번째 스텝으로 가는 것이 아니라 중간 스텝(u*)을 먼저 구함
이때 u*를 구하는 과정은 Navier-Stokes식에서 압력 p 텀을 무시한 형태
(1)에서 구한 u*를 이용해 n+1번째 스텝 계산
여기서 문제는 압력 p텀을 알아야 하는데, Incompressible 조건을 이용해 아래와 같이 구할 수 있음

\[u^{n+1} - u^{*} = - \frac{\Delta t}{\rho} \bigtriangledown p^{n+1}\] \[\bigtriangledown \left (u^{n+1} - u^{*} \right ) = - \frac{\Delta t}{\rho} \bigtriangledown^{2} p^{n+1}\] \[- \bigtriangledown u^{*} = - \frac{\Delta t}{\rho} \bigtriangledown^{2} p^{n+1} \quad \left ( \because \bigtriangledown u^{n+1} = 0 \right )\] \[\bigtriangledown^{2} p^{n+1} = \frac{\rho}{\Delta t} \bigtriangledown \cdot u^{*}\]

Incompressible 하다는 것은 Divergence가 0이라는 의미!
- (중간 단계 u*에서는 0이 아니고, 최종 단계에서 0)
이렇게 구한 p를 (2)식에 대입해 u의 n+1 스텝을 계산

Grid Texture	FDM

참고로 미분 관련 연산은 모두 FDM 방식을 이용해 계산
이때 각 그리드 셀의 속도, 압력 벡터들을 텍스처 형태로 들고다니면 매우매우 편리해짐
- 속도의 경우 그리드 셀 가운데 하나만 정의해도 괜찮음
- 반면 압력의 경우 그리드 셀의 모서리마다 벡터가 정의되어야 내부의 변화량을 측정하기 편함

Stable Fluid Algorithm

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

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

Sourcing : External Forces 계산 (예제에서는 마우스의 속도, 그리드 셀의 컬러)
Advection : Velocity Field 업데이트
Diffuse : viscosity(점성) 계산 후 속도의 확산
- 여기까지 계산한 결과가 u*
Projection : pressure(압력) 계산 & 최종 속도 업데이트

Externel Forces

[numthreads(32, 32, 1)]
void main(int3 gID : SV_GroupID, int3 gtID : SV_GroupThreadID,
          uint3 dtID : SV_DispatchThreadID)
{
    int radius = 50;
    float dist = length(float2(dtID.xy) - float2(i, j)) / radius;
    float scale = smootherstep(1.0 - dist);

    velocity[dtID.xy] += sourcingVelocity * scale;
    density[dtID.xy] += sourcingDensity * scale;
    
}

Sourcing 단계에서는 간단하게 sourcingVelocity와 sourcingDensity를 넘겨주면 됨
sourcingVelocity의 경우 마우스의 속도
sourcingDensity는 이것저것 바꿔가며 멋진 걸로..

Advection

fluid-images11-en

[numthreads(32, 32, 1)]
void main(int3 gID : SV_GroupID, int3 gtID : SV_GroupThreadID,
          uint3 dtID : SV_DispatchThreadID)
{
    uint width, height;
    velocity.GetDimensions(width, height);
    float2 dx = float2(1.0 / width, 1.0 / height);
    float2 pos = (dtID.xy + 0.5) * dx; // point a
    
    float2 vel = velocityTemp.SampleLevel(pointWrapSS, pos, 0).xy; // velocity at point a
    
    float2 posBack = pos - vel * dt;

    velocity[dtID.xy] = velocityTemp.SampleLevel(linearWrapSS, posBack, 0);
    density[dtID.xy] = densityTemp.SampleLevel(linearWrapSS, posBack, 0);
}

Advection 단계에서는 가장 기본적인 작업인 Velocity Field를 업데이트 해야한다. 이때 비교적 간단하고 안정적이라고 알려진 Semi-Lagrangian Method를 이용할 수 있다. 컨셉을 한 마디로 정리하면, 현재 위치에서의 속도 벡터를 역으로 추적한 곳의 속도 벡터를 가져오는 것이다. 위 그림으로 보면 b의 속도 벡터를 a위치로 가져오는 것이다. 다만 이 방법은 advection이 선형이라는 가정이 필요하다. 따라서 역추적한 곳의 속도 벡터가 선형이 아니라면 오차가 발생하게 되는데 이를 줄이기 위해 BFECC(Back and Forth Error Compensation and Correction) 이라는 방법이 제안되었다.

BFECC Algorithm	Result

[numthreads(32, 32, 1)]
void main(int3 gID : SV_GroupID, int3 gtID : SV_GroupThreadID,
          uint3 dtID : SV_DispatchThreadID)
{
    uint width, height;
    velocity.GetDimensions(width, height);
    float2 dx = float2(1.0 / width, 1.0 / height);
    float2 pos = (dtID.xy + 0.5) * dx; 
    
    float2 vel = velocityTemp.SampleLevel(pointWrapSS, pos, 0).xy;
    
    // back trace
    float2 posBack1 = pos - vel * dt;
    vel = velocityTemp.SampleLevel(pointWrapSS, posBack1, 0).xy;
    
    // forward trace
    float2 posNew1 = posBack1 + vel * dt;
    float2 error = posNew1 - pos;
    
    float2 posNew2 = pos - error / 2.0;
    vel = velocityTemp.SampleLevel(pointWrapSS, posNew2, 0).xy;
    
    // back trace 2
    float2 posBack2 = posNew2 - vel * dt;
    
    velocity[dtID.xy] = velocityTemp.SampleLevel(pointWrapSS, posBack2, 0);
    density[dtID.xy] = densityTemp.SampleLevel(linearWrapSS, posBack2, 0);
}

Diffuse

\[\frac{\partial u}{\partial t} = \nu \bigtriangledown ^{2} u\] \[\frac{u_{i,j}^{*} - u_{i,j}^{n}}{\Delta t} = \nu \bigtriangledown ^{2} u_{i,j}^{n}\] \[Note \; that \; FDM \; of \; Laplacian \quad \bigtriangledown ^{2} u_{i,j}^{n} = \frac{u_{i+1,j}^{n} + u_{i-1,j}^{n} + u_{i,j+1}^{n} + u_{i,j-1}^{n} - 4 u_{i,j}^{n}}{\Delta x^{2}}\] \[\therefore \quad u_{i,j}^{*} = u_{i,j}^{n} + \nu \Delta t \left ( \frac{u_{i+1,j}^{n} + u_{i-1,j}^{n} + u_{i,j+1}^{n} + u_{i,j-1}^{n} - 4 u_{i,j}^{n}}{\Delta x^{2}} \right )\]

이때 delta x = 1이고, 다음 스텝의 기울기를 사용하는 Implicit Integration 방법을 사용할 것이기 때문에 Laplacian의 마지막 4u^n는 중간 스텝(u)으로 변경한다. (실험적으로 더 안정적이어서..?)

\[u_{i,j}^{*} = u_{i,j}^{n} + \nu \Delta t \left ( u_{i+1,j}^{n} + u_{i-1,j}^{n} + u_{i,j+1}^{n} + u_{i,j-1}^{n} - 4 u_{i,j}^{*} \right )\] \[u_{i,j}^{*} \left ( 1 - 4 \nu \Delta t \right ) = u_{i,j}^{n} + \nu \Delta t \left ( u_{i+1,j}^{n} + u_{i-1,j}^{n} + u_{i,j+1}^{n} + u_{i,j-1}^{n} \right )\] \[u_{i,j}^{*} = \frac{ u_{i,j}^{n} + \nu \Delta t \left ( u_{i+1,j}^{n} + u_{i-1,j}^{n} + u_{i,j+1}^{n} + u_{i,j-1}^{n} \right ) } {1 - 4 \nu \Delta t}\]

[numthreads(32, 32, 1)]
void main(int3 gID : SV_GroupID, int3 gtID : SV_GroupThreadID,
          uint3 dtID : SV_DispatchThreadID)
{
    uint width, height;
    velocity.GetDimensions(width, height);
    
    uint2 left = uint2(dtID.x == 0 ? width - 1 : dtID.x - 1, dtID.y);
    uint2 right = uint2(dtID.x == width - 1 ? 0 : dtID.x + 1, dtID.y);
    uint2 up = uint2(dtID.x, dtID.y == height - 1 ? 0 : dtID.y + 1);
    uint2 down = uint2(dtID.x, dtID.y == 0 ? height - 1 : dtID.y - 1);
    
    // Explicit integration
    velocity[dtID.xy] = (velocityTemp[dtID.xy] + viscosity * dt * (velocityTemp[left] + 
                          velocityTemp[right] + velocityTemp[up] + velocityTemp[down])) / (1.0 + 4 * viscosity * dt);
}

Viscosity를 구한 후 Diffuse하는 과정에서 주의해야 할 점이 있다. 단 한번의 연산 만으로는 정확한 Diffuse가 표현되지 않는다는 것이다. 따라서 정확하고 안정적인 결과을 얻기 위해서는 Jacobi Iteration 방법을 사용한 반복적인 계산을 진행해야 한다. 복잡한 내용은 아니고 UAV 텍스처 두개를 준비한 후 홀짝 번갈아가며 Diffuse Compute Shader를 실행하면 된다.

// ...
context->CSSetShader(m_diffuseCS.Get(), 0, 0);

for (int i = 0; i < 10; i++) {

    if (i % 2 == 0) {
        context->CSSetShaderResources(0, 2, evenSRVs);
        context->CSSetUnorderedAccessViews(0, 2, evenUAVs, NULL);
    } else {
        context->CSSetShaderResources(0, 2, oddSRVs);
        context->CSSetUnorderedAccessViews(0, 2, oddUAVs, NULL);
    }
      
    context->Dispatch(UINT(ceil(m_width / 32.0f)),
                      UINT(ceil(m_height / 32.0f)), 1);
      
    ComputeShaderBarrier(context);
}

Projection

void StableFluids::Projection(ComPtr<ID3D11DeviceContext> &context) {
    // Compute divergence
    context->CSSetShaderResources(0, 1, m_velocity.GetAddressOfSRV());

    ID3D11UnorderedAccessView *uavs[3] = {
        m_divergence.GetUAV(), m_pressure.GetUAV(), m_pressureTemp.GetUAV()};

    context->CSSetUnorderedAccessViews(0, 3, uavs, NULL);
    context->CSSetShader(m_divergenceCS.Get(), 0, 0);
    context->Dispatch(UINT(ceil(m_width / 32.0f)), UINT(ceil(m_height / 32.0f)),
                      1);
    ComputeShaderBarrier(context);

    // Jacobi iteration
    context->CSSetShader(m_jacobiCS.Get(), 0, 0);

    for (int i = 0; i < 100; i++) {
        if (i % 2 == 0) {
            context->CSSetShaderResources(0, 1, m_pressure.GetAddressOfSRV());
            context->CSSetUnorderedAccessViews(
                0, 1, m_pressureTemp.GetAddressOfUAV(), NULL);
        } else {
            context->CSSetShaderResources(0, 1,
                                          m_pressureTemp.GetAddressOfSRV());
            context->CSSetUnorderedAccessViews(
                0, 1, m_pressure.GetAddressOfUAV(), NULL);
        }
        context->CSSetShaderResources(1, 1, m_divergence.GetAddressOfSRV());
        context->Dispatch(UINT(ceil(m_width / 32.0f)),
                          UINT(ceil(m_height / 32.0f)), 1);
        ComputeShaderBarrier(context);
    }

    // Apply pressure
    context->CSSetShaderResources(0, 1, m_pressure.GetAddressOfSRV());
    context->CSSetUnorderedAccessViews(0, 1, m_velocity.GetAddressOfUAV(),
                                       NULL);
    context->CSSetShader(m_applyPressureCS.Get(), 0, 0);
    context->Dispatch(UINT(ceil(m_width / 32.0f)), UINT(ceil(m_height / 32.0f)),
                      1);
    ComputeShaderBarrier(context);
}

u*의 Divergence 먼저 계산
Jacobi Iteration을 이용해 n+1 스텝에서의 p 계산
- 압력은 점성에비해 중요도가 매우 크기때문에 반복 횟수도 훨씬 많이 지정
계산한 p를 기반으로 최종 속도 업데이트

\[\bigtriangledown^{2} p^{n+1} = \frac{\rho}{\Delta t} \bigtriangledown \cdot u^{*}\] \[\frac{p_{i+1,j}^{n} + p_{i-1,j}^{n} + p_{i,j+1}^{n} + p_{i,j-1}^{n} - 4 p_{i,j}^{n+1}}{\Delta t} = \frac{\rho}{\Delta t} \bigtriangledown \cdot u_{i,j}^{*}\] \[p_{i,j}^{n+1} = \frac{\rho}{4} \left ( - \bigtriangledown \cdot u^{*} + p_{i+1,j}^{n} + p_{i-1,j}^{n} + p_{i,j+1}^{n} + p_{i,j-1}^{n} \right )\]

[numthreads(32, 32, 1)]
void main(int3 gID : SV_GroupID, int3 gtID : SV_GroupThreadID,
          uint3 dtID : SV_DispatchThreadID)
{
    // Dirichlet boundary condition
    if (dtID.x == 0 && dtID.y == 0)
    {
        pressureOut[dtID.xy] = 0.0;
        return;
    }
    
    uint width, height;
    pressureOut.GetDimensions(width, height);
    float2 dx = float2(1.0 / width, 1.0 / height);
        
    uint2 left = uint2(dtID.x == 0 ? width - 1 : dtID.x - 1, dtID.y);
    uint2 right = uint2(dtID.x == width - 1 ? 0 : dtID.x + 1, dtID.y);
    uint2 up = uint2(dtID.x, dtID.y == height - 1 ? 0 : dtID.y + 1);
    uint2 down = uint2(dtID.x, dtID.y == 0 ? height - 1 : dtID.y - 1);
    
    pressureOut[dtID.xy] = (-divergence[dtID.xy] + pressureTemp[left] + pressureTemp[right] +
                             pressureTemp[up] + pressureTemp[down]) * 0.25;
}

Result

HongLabGraphicsExample2024-10-0517-58-35-ezgif com-crop

HongLabGraphicsExample2024-10-0522-01-57-ezgif com-crop

맨 위로 이동하기

등딱지빌런