Hybrid Simulation

Date: 2024.12.03 Updated: 2024.12.03

카테고리: Graphics

태그: DirectX

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

Hybrid Simulation

Particle Based	Grid Based

Particle Based Simulation
- 주변 입자들을 이용해 각 입자의 물리량 계산
- 장점) 유체의 경계, 복잡한 경계 구조를 잘 표현함
- 단점) 모든 particle을 돌며 주변 입자를 찾아야하기 때문에 연산이 매우 비쌈
- 비압축성(divergence=0) 성질을 주변 입자로만 계산하기 때문에 전체 시스템적으로 불안정함
Grid Based Simulation
- 고정된 격자마다 물리량을 계산
- 장점) 주변 격자의 위치를 찾을 필요가 없기 때문에 연산이 매우 효율적
- 단점) memory bandwidth 높음, 복잡한 경계 구조를 표현하기 어려움

앞서 Fluids Simulation을 구현하는 두 가지 방법에 대해 알아봤다. Hybrid Simulation은 입자, 격자 기반의 성질을 모두 이용한 방법이다. 시스템 전체의 비압축성을 적용하기 위해 Divergence를 포함한 물리량을 격자 기반으로 계산하고, 이를 통해 개별 입자의 속도를 최종적으로 업데이트 한다. 이때 Bitonic-Sort 알고리듬을 통해 개별 입자가 어느 그리드 셀에 포함되는지를 병렬적으로 계산한다. 따라서 주변 입자를 하나하나 찾을 필요 없이 이웃 그리드 셀에 포함된 입자들에 대해서만 연산을 진행하면 된다.

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

    Projection();
    ParticleStep();
}

Projection : 격자 기반 연산
- Boundary Condition 세팅
- 격자 기반 Divergence, Pressure, Velocity 연산
ParticleStep : 입자 기반 연산
- 입자가 속한 그리드 셀 기반으로 velocity, Position 업데이트
- 입자-그리드셀 매칭 정보 업데이트 & 정렬
- 입자기반 정보를 바탕으로 그리드셀 정보 업데이트

Projection

Compute Divergence

void main(uint3 dtID : SV_DispatchThreadID)
{
    if (dtID.x == 0 || dtID.y == 0 || dtID.z == 0
        || dtID.x == width - 1 || dtID.y == height - 1 || dtID.z == depth - 1
        || density[dtID] < 1e-6)
    {
        bc[dtID] = -1; // Dirichlet
    }
    else
    {
        bc[dtID] = 0;
        
        float div = 0.0;

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

            if (bc[dtID.xyz + offset[i]] == -2) // Neumann
                div += dot(2 * velocityTemp[dtID.xyz + offset[i]].xyz - velocityTemp[dtID.xyz].xyz, float3(offset[i]));

            else
                div += dot(velocityTemp[dtID.xyz + offset[i]].xyz, float3(offset[i]));
        }

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

Smoke Simulation에서 구현한 것과 마찬가지로, 3D 그리드 셀 기반으로 Divergence 계산
이웃 셀이 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}\]

Apply Pressure

void main(uint3 dtID : SV_DispatchThreadID)
{
    if (bc[dtID.xyz] >= 0)
    {
        float temp = 0.0;
        
        [unroll]
        for (int i = 0; i < 6; i++)
        {
            if (bc[dtID.xyz + offset[i]] == -1) // Dirichlet
            {
                temp += -pressure[dtID.xyz];
            }
            else if (bc[dtID.xyz + offset[i]] == -2) // Neumann
            {
                temp += pressure[dtID.xyz];
            }
            else
                temp += pressure[dtID.xyz + offset[i]];
        }

        pressure[dtID.xyz] = (-divergence[dtID.xyz] + temp) / 6.0;
    }
}

먼저 Jacobi Iteration 기법을 이용해 격자 기반 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);
    }
}

업데이트 된 Pressure를 이용해 격자 기반 Velocity 업데이트\
여기까지의 과정을 마무리하면 격자 기반 Divergence, Pressure, Velocity 연산 완료
이때 사용한 격자 Texture는 보통 DownSampling 되어 있으므로 업데이트 전 데이터도 저장해야 함
- Upsampling할 때 저해상도의 변화율을 적용해야 부드러운 결과 연출 가능하기 때문

ParticleStep

191421

Hybrid Simulation의 동작을 간략하게 이야기하면 위 그림과 같다. SortElemnts라는 구조체 배열에는 그리드 셀의 인덱스를 의미하는 key, 파티클 배열의 인덱스를 의미하는 value 값이 들어있다. 매 프레임 Key기준 Bitonic sort를 이용해 SortElements 배열을 정렬하면 특정 그리드 셀에 포함되는 파티클들을 빠르게 구할 수 있다. 그리드 기반으로 계산한 물리량을 해당 그리드에 포함된 파티클에 적용시키고, 주변 파티클들과 입자 기반 연산을 수행한다. 그 후 움직인 파티클들을 기준으로 다시 그리드 셀의 정보를 수정해주어야 한다.

Upsampling with Grid to Particle

RWStructuredBuffer<Particle> particles : register(u0);
RWStructuredBuffer<SortElement> sortElements : register(u1);

[numthreads(1024, 1, 1)]
void main(uint3 dtID : SV_DispatchThreadID)
{
    float coeff = 0.97;
    
    SortElement ref = sortElements[dtID.x];
    
    float3 uvw = particles[ref.value].pos.xyz * dxBase;
    float d = densityNew.SampleLevel(linearClampSampler, uvw, 0).r;
    float s = signedDistance.SampleLevel(linearClampSampler, uvw, 0).r;
    
    if (d > 1e-5 && s < -0.1 && ref.key != INACTIVE) // Active particle
    {
        float3 velOld = velocityOld.SampleLevel(linearClampSampler, uvw, 0).xyz;
        float3 velNew = velocityNew.SampleLevel(linearClampSampler, uvw, 0).xyz;
        
        particles[ref.value].vel = lerp(velNew, particles[ref.value].vel + velNew - velOld, coeff);
    }
}

DownSampling된 공간에서 계산한 Grid 기반 정보를 Upsampling
이때 저해상도 데이터의 변화율을 기반으로 업데이트하면 부드러운 결과 연출
0.97 Coefficient를 곱해 속도가 점점 줄어들도록 구현

Update Particles

void main(uint3 gId : SV_GroupID, uint3 dtID : SV_DispatchThreadID)
{
    SortElement ref = sortElements[dtID.x];
    
    if (ref.key != INACTIVE)
    {
        Particle p = particles[ref.value];
   
        float3 uvw = p.pos.xyz * dxBase;

        if (IsInside(uvw))
        {
            p.vel += float3(0, -0.2, 0) * dt * float3(width, height, depth); // 중력
            p.pos += p.vel * dt;
            sortElements[dtID.x].key = MinCornerCell(particles[ref.value].pos);
            particles[ref.value] = p;
        }

        else
        {
            sortElements[dtID.x].key = INACTIVE;
        }
    }

    else if (dtID.x < numActiveParticles[0] + numNewParticles)
    {
        // Add a new particle
    }
}

경계를 벗어난 파티클에 대응되는 SortElements의 key값은 INACTIVE(int 자료형의 최대값)으로 설정
- 정렬시 알아서 맨 뒤로 밀려나도록
파티클별로 External Forces(중력, …) 추가
- 이때 파티클 좌표계에 따라 UpSampling 기준 Scale을 곱해줘야 함
경계를 벗어나는 입자들은 INACTIVE로 꺼주기
필요한 경우 조건에 따라 새로운 파티클 생성

Find First Index

210923

uint3 Index3(in uint i1)
{
    uint width, height, depth;
    firstIndex.GetDimensions(width, height, depth);
    
    uint3 i3;
    i3.z = i1 / (width * height);
    i1 -= i3.z * (width * height);
    i3.y = i1 / width;
    i1 -= i3.y * width;
    i3.x = i1;
    
    return i3;
}

[numthreads(1024, 1, 1)]
void main(uint3 dtID : SV_DispatchThreadID)
{
    uint numStructs, stride;
    sortElements.GetDimensions(numStructs, stride);
    
    if (dtID.x < numStructs && sortElements[dtID.x].key != INACTIVE)
    {
        if (dtID.x == numStructs - 1)
        {
            numActiveParticles[0] = numStructs;
        }

        else if (sortElements[dtID.x + 1].key == INACTIVE)
        {
            numActiveParticles[0] = dtID.x + 1;
        }
        
        if (dtID.x == 0 || sortElements[dtID.x - 1].key != sortElements[dtID.x].key)
        {
            uint3 i3 = Index3(sortElements[dtID.x].key);
            firstIndex[i3] = dtID.x;
        }
    }
}

특정 그리드 셀에 포함된 파티클들을 빠르게 접근하기 위해 firstIndex Buffer 세팅
- first index부터 key값이 바뀌는 구간까지가 해당 그리드 셀에 포함된 파티클들
sortElements 한 바퀴 도는 김에 활성화된 파티클 갯수도 구함
- dtID.x == numStructs - 1 : 모든 파티클 ACTIVE 상태
- sortElements[dtID.x + 1].key == INACTIVE : INACTIVE 상태의 파티클 존재
인덱스를 다시 3D 좌표로 변환해주는 Index3 함수 이용

Particle to Grid Cell

void main(uint3 dtID : SV_DispatchThreadID)
{
    uint width, height, depth;
    firstIndex.GetDimensions(width, height, depth);
    
    uint numStructs, stride;
    sortElements.GetDimensions(numStructs, stride);
    
    density[dtID] = 0.0;
    velocity[dtID] = float4(0, 0, 0, 0);
    
    float3 posTemp = float3(0, 0, 0);
    
    if (dtID.x > 0 && dtID.y > 0 && dtID.z > 0
        && dtID.x < width - 1 && dtID.y < height - 1 && dtID.z < depth - 1)
    {
        float3 cellCenter = dtID + 0.5;
        
        for (int k = -1; k < 1; k++)
            for (int j = -1; j < 1; j++)
                for (int i = -1; i < 1; i++)
                {
                    uint firstId = firstIndex[dtID + int3(i, j, k)];
                    
                    if (firstId != INACTIVE)
                    {
                        int p = 0;

                        [loop]
                        while (true)
                        {
                            SortElement ref = sortElements[firstId + p];
                        
                            if (ref.key == INACTIVE)
                                break;

                            if (ref.key != sortElements[firstId].key)
                                break;
                        
                            float dist = length(particles[ref.value].pos - cellCenter);
                            float weight = CubicSpline(dist * 2.0);
                            density[dtID] += weight;
                            velocity[dtID] += weight * float4(particles[ref.value].vel, 0.0);
                            
                            p++;
                        
                            if (firstId + p == numStructs)
                                break;
                        }
                    }
                }

        if (density[dtID] > 1e-6)
        {
            velocity[dtID] /= density[dtID];
        }
    }
}

업데이트된 파티클들의 정보를 다시 그리드셀로 적용
이웃 그리드 셀에 포함된 파티클들을 순회하며 밀도, 속도 wieghted sum
- 각 이웃 그리드 셀의 인덱스는 정육면체 형태의 꼭지점 8개에 저장
weighted sum 결과를 바탕으로 그리드 셀 정보 업데이트

Results

HongLabGraphicsExample2024-12-0422-01-11-ezgif com-crop

맨 위로 이동하기

등딱지빌런

Hybrid Simulation

Hybrid Simulation

Projection

Compute Divergence

Apply Pressure

ParticleStep

Upsampling with Grid to Particle

Update Particles

Find First Index

Particle to Grid Cell

Results

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

댓글 남기기

최근 글 10 개 :)