three.js GPGPU SPH - Fluid Simulation 구현

 

이번에는 integration에서 넘어가서 density & pressure에 대한 연산을 다루어보자.

밀도와 압력은 particle 시뮬레이션에서 유체의 느낌을 주기 위해서 필요하다. 

 

Density Pressure Function

밀도 압력 함수를 다룰 때,

유체 시뮬레이션에서는 다양한 커널이 존재한다.

이건 논문을 뒤져봐야하는 어려운 내용이므로 개념을 다루지는 않고, 어떻게 코드로 짤 것인가에 집중한다.

 

커널의 종류도 많아서 이걸 다 glsl에다 구현해야할 것 같다. 

위 세 함수는 결과적으로 compute force에서 사용하게 될 것이다.

 

Compute Force Function

또한 입자간 서로 밀어내는 힘에 대한 함수도 shader에서 필요하다. 

 

이렇게 위 크게 두 부류의 density & pressure, compute force shader를 만들고, 

integration과 같이 매 프레임마다 compute 해주는 코드를 구현하도록 하자.

 

Density Pressure Shader 구현하기 

일단 위의 함수에서 density pressure 두 개의 변수가 결과로 저장된다.

그러면 1 버전의 glsl에서는 둘 다 꺼낼 수 있도록 하나의 vector에 저장해서 출력해야한다. 

 

precision mediump float;

uniform int particleLength
uniform float particleMass
uniform float viscosity
uniform float gasConstant
uniform float restDensity
uniform float boundDamping
uniform float pi
uniform vec3 boxSize
uniform float radius
uniform float radius2
uniform float radius3
uniform float radius4
uniform float radius5
uniform float timestep

float StdKernel(float distanceSquared) {
  float x = 1.0f - distanceSquared / radius2;
  return 315.0 / (64.0 * pi * radius3) * pow(x, 3.0);
}

void main() {
    int id = int(gl_FragCoord.x);

    // 파티클 위치 읽기
    vec3 origin = texture2D(positionTexture, vec2(float(id) / float(particleLength), 0.0)).xyz;
    float sum = 0.0;

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

        vec3 otherPos = texture2D(positionTexture, vec2(float(i) / float(particleLength), 0.0)).xyz;
        vec3 diff = origin - otherPos;
        float distanceSquared = dot(diff, diff);
        
        if (radius2 * 0.004 > distanceSquared * 0.004) {
            sum += StdKernel(distanceSquared * 0.004);
        }
    }

    // 파티클 밀도 및 압력 계산
    float density = sum * particleMass + 0.000001;
    float pressure = gasConstant * (density - restDensity);

    // 결과 출력 (압력, 밀도를 텍스처에 쓰는 방법은 WebGL의 FBO를 통해 처리해야 함)
    // gl_FragColor로 간단히 출력
    gl_FragColor = vec4(density, pressure, 0.0, 1.0); 
}

이렇게 glsl shader를 작성하고 vector에 density와 pressure를 둘 다 담아서 FragColor로 return한다. 

 

particle texture의 uv 좌표

보면서 헷갈렸던 점이 있다.

이전에 integrate function을 작성할 때는

 

texture2D(positionTexture, vec2(gl_FragCoord.xy / resolution.xy)).xyz;

이런식으로 texture의 texel값으로 파티클 postion을 가져왔다.

 

vec3 otherPosition = texture2D(positionTexture, vec2(float(i) / float(particleLength), 0.0)).xyz;

다른 particle position에 접근하기 위해서 현재 이런식으로 접근하고 있는데,

이렇게 하는게 맞는 건가?

 

this.gpuCompute = new GPUComputationRenderer(this.totalParticles, 1, this.renderer);

texture 에 uv로 접근한다는 점에서 현재 해상도를 y는 1로 넣고, x로 크기를 설정하고 있다.

 

resolution의 x좌표는 particle의 총 개수로 들어가 있는 상태이고,

y좌표는 1로 고정해놓기 때문에 결국 particleLength와 resolution.x는 같다.

 

gl_FragCoord.x의 경우도 마찬가지로 x좌표 중 몇 번째인지 index를 담고 있어서

gl_FragCoord.x / resolution.x가 0~1 사이의 uv 좌표가 되는 것이다.

 

만약 다른 particle의 위치를 원한다면 이 index를 다른 것으로 바꾸면 된다.

따라서 i / particleLength의 방법이 충분히 가능하다는 이론이다.

 

로직상으로는 그렇게 이해가 되는데 실제로 동작해서 문제가 되지 않는지 확인해보자. 

 

ComputeForce

precision highp float;

uniform int particleLength
uniform float particleMass
uniform float viscosity
uniform float gasConstant
uniform float restDensity
uniform float boundDamping
uniform float pi
uniform vec3 boxSize
uniform float radius
uniform float radius2
uniform float radius3
uniform float radius4
uniform float radius5
uniform float timestep

// uniform sampler2D forceTexture;     // 힘을 입력으로 받음

// 첫 번째 도함수
float SpikyKernelFirstDerivative(float distance) {
    float x = 1.0 - distance / radius;
    return -45.0 / (pi * radius4) * x * x;
}

// 두 번째 도함수
float SpikyKernelSecondDerivative(float distance) {
    float x = 1.0 - distance / radius;
    return 90.0 / (pi * radius5) * x;
}

// 그라디언트 계산 함수
vec3 SpikyKernelGradient(float distance, vec3 direction) {
    return SpikyKernelFirstDerivative(distance) * direction;
}

void main() {
    int id = int(gl_FragCoord.x);  // 각 픽셀은 파티클 ID와 매핑

    if (id >= particleLength) {
        gl_FragColor = vec4(0.0);  // 범위 밖일 경우 아무런 작업도 하지 않음
        return;
    }

    vec3 origin = texture2D(positionTexture, vec2(float(id) / float(particleLength), 0.0)).xyz;
    float density = texture2D(densityPressureTexture, vec2(float(id) / float(particleLength), 0.0)).x;
    float pressure = texture2D(densityPressureTexture, vec2(float(id) / float(particleLength), 0.0)).y
    float density2 = density * density;
    float mass2 = particleMass * particleMass;
    vec3 resultPressure = vec3(0.0);
    vec3 visc = vec3(0.0);

    for (int i = 0; i < particleLength; i++) {
        vec3 otherPosition = texture2D(positionTexture, vec2(float(i) / float(particleLength), 0.0)).xyz;
        
        // 자기 자신과 비교하지 않기
        if (origin == otherPosition) continue;

        float dist = distance(origin, otherPosition);
        if (dist < radius * 2.0) {
            // 압력 계산
            vec3 pressureGradientDirection = normalize(origin - otherPosition);

            vec3 pressureContribution = mass2 * SpikyKernelGradient(dist, pressureGradientDirection);
            float otherDensity = texture2D(densityPressureTexture, vec2(float(i) / float(particleLength), 0.0)).x;
            float otherPressure = texture2D(densityPressureTexture, vec2(float(i) / float(particleLength), 0.0)).y;
            pressureContribution *= (pressure / density2) + otherPressure / (otherDensity * otherDensity);

            // 점성 계산
            vec3 velocity = texture2D(velocityTexture, vec2(float(id) / float(particleLength), 0.0)).xyz;
            vec3 otherVelocity = texture2D(velocityTexture, vec2(float(i) / float(particleLength), 0.0)).xyz;
            vec3 viscosityContribution = viscosity * mass2 * (otherVelocity - velocity) / otherDensity;
            viscosityContribution *= SpikyKernelSecondDerivative(dist);
            
            resultPressure += pressureContribution;
            visc += viscosityContribution;
        }
    }

    // 중력 힘 적용
    vec3 currentForce = vec3(0, -9.81 * particleMass, 0) - resultPressure + visc;
    gl_FragColor = vec4(currentForce, 1.0);
}

일단 두 셰이더를 glsl로 포워딩하였고, 

이제 shader를 스크립트에 포함시키는 코드를 만든다. 

 

    this.initialPosition = this.gpuCompute.createTexture();
    this.initialVelocity = this.gpuCompute.createTexture();
    this.initialForce = this.gpuCompute.createTexture();
    this.initParticlePosition();
    this.initialDenstPress = this.gpuCompute.createTexture();

    this.positionVariable = this.gpuCompute.addVariable('positionTexture', computeIntegratePosition, this.initialPosition);
    this.velocityVariable = this.gpuCompute.addVariable('velocityTexture', computeIntegrateVelocity, this.initialVelocity);
    this.forceVariable = this.gpuCompute.addVariable('forceTexture', computeForce, this.initialForce);
    this.denstPressVariable = this.gpuCompute.addVariable('densityPressureTexture', computeDensityPressure, this.initialDenstPress);

    this.gpuCompute.setVariableDependencies(this.positionVariable, [this.positionVariable, this.velocityVariable]);
    this.gpuCompute.setVariableDependencies(this.velocityVariable, [this.positionVariable, this.velocityVariable, this.forceVariable]);
    this.gpuCompute.setVariableDependencies(this.forceVariable, [this.positionVariable, this.velocityVariable, this.forceVariable, this.denstPressVariable]);
    this.gpuCompute.setVariableDependencies(this.denstPressVariable, [this.positionVariable]);

    //#region position uniform
    this.positionVariable.material.uniforms.particleLength = { value: this.totalParticles };
    this.positionVariable.material.uniforms.particleMass = { value: this.particleMass };
    this.positionVariable.material.uniforms.viscosity = { value: this.viscosity };
    this.positionVariable.material.uniforms.gasConstant = { value: this.gasConstant };
    this.positionVariable.material.uniforms.restDensity = { value: this.restingDensity };
    this.positionVariable.material.uniforms.boundDamping = { value: this.boundDamping };
    this.positionVariable.material.uniforms.pi = { value: Math.PI };
    this.positionVariable.material.uniforms.boxSize = { value: this.boxSize };

    this.positionVariable.material.uniforms.radius = { value: this.particleRadius };
    this.positionVariable.material.uniforms.radius2 = { value: Math.pow(this.particleRadius, 2) };
    this.positionVariable.material.uniforms.radius3 = { value: Math.pow(this.particleRadius, 3) };
    this.positionVariable.material.uniforms.radius4 = { value: Math.pow(this.particleRadius, 4) };
    this.positionVariable.material.uniforms.radius5 = { value: Math.pow(this.particleRadius, 5) };
    
    this.positionVariable.material.uniforms.timestep = { value: this.timestep };
    //#endregion

    //#region velocity uniform
    this.velocityVariable.material.uniforms.particleLength = { value: this.totalParticles };
    this.velocityVariable.material.uniforms.particleMass = { value: this.particleMass };
    this.velocityVariable.material.uniforms.viscosity = { value: this.viscosity };
    this.velocityVariable.material.uniforms.gasConstant = { value: this.gasConstant };
    this.velocityVariable.material.uniforms.restDensity = { value: this.restingDensity };
    this.velocityVariable.material.uniforms.boundDamping = { value: this.boundDamping };
    this.velocityVariable.material.uniforms.pi = { value: Math.PI };
    this.velocityVariable.material.uniforms.boxSize = { value: this.boxSize };
    
    this.velocityVariable.material.uniforms.radius = { value: this.particleRadius };
    this.velocityVariable.material.uniforms.radius2 = { value: Math.pow(this.particleRadius, 2) };
    this.velocityVariable.material.uniforms.radius3 = { value: Math.pow(this.particleRadius, 3) };
    this.velocityVariable.material.uniforms.radius4 = { value: Math.pow(this.particleRadius, 4) };
    this.velocityVariable.material.uniforms.radius5 = { value: Math.pow(this.particleRadius, 5) };
    
    this.velocityVariable.material.uniforms.timestep = { value: this.timestep };
    //#endregion

    //#region force uniform
    this.forceVariable.material.uniforms.particleLength = { value: this.totalParticles };
    this.forceVariable.material.uniforms.particleMass = { value: this.particleMass };
    this.forceVariable.material.uniforms.viscosity = { value: this.viscosity };
    this.forceVariable.material.uniforms.gasConstant = { value: this.gasConstant };
    this.forceVariable.material.uniforms.restDensity = { value: this.restingDensity };
    this.forceVariable.material.uniforms.boundDamping = { value: this.boundDamping };
    this.forceVariable.material.uniforms.pi = { value: Math.PI };
    this.forceVariable.material.uniforms.boxSize = { value: this.boxSize };
    
    this.forceVariable.material.uniforms.radius = { value: this.particleRadius };
    this.forceVariable.material.uniforms.radius2 = { value: Math.pow(this.particleRadius, 2) };
    this.forceVariable.material.uniforms.radius3 = { value: Math.pow(this.particleRadius, 3) };
    this.forceVariable.material.uniforms.radius4 = { value: Math.pow(this.particleRadius, 4) };
    this.forceVariable.material.uniforms.radius5 = { value: Math.pow(this.particleRadius, 5) };
    
    this.forceVariable.material.uniforms.timestep = { value: this.timestep };
    //#endregion

    //#region density & pressure uniform
    this.denstPressVariable.material.uniforms.particleLength = { value: this.totalParticles };
    this.denstPressVariable.material.uniforms.particleMass = { value: this.particleMass };
    this.denstPressVariable.material.uniforms.viscosity = { value: this.viscosity };
    this.denstPressVariable.material.uniforms.gasConstant = { value: this.gasConstant };
    this.denstPressVariable.material.uniforms.restDensity = { value: this.restingDensity };
    this.denstPressVariable.material.uniforms.boundDamping = { value: this.boundDamping };
    this.denstPressVariable.material.uniforms.pi = { value: Math.PI };
    this.denstPressVariable.material.uniforms.boxSize = { value: this.boxSize };
    
    this.denstPressVariable.material.uniforms.radius = { value: this.particleRadius };
    this.denstPressVariable.material.uniforms.radius2 = { value: Math.pow(this.particleRadius, 2) };
    this.denstPressVariable.material.uniforms.radius3 = { value: Math.pow(this.particleRadius, 3) };
    this.denstPressVariable.material.uniforms.radius4 = { value: Math.pow(this.particleRadius, 4) };
    this.denstPressVariable.material.uniforms.radius5 = { value: Math.pow(this.particleRadius, 5) };
    
    this.denstPressVariable.material.uniforms.timestep = { value: this.timestep };
    //#endregion

이렇게 density pressure shader에 대해서도 compute 과정을 완료하였다.

 

그런데 뭔가 이상해졌다.. 

 

NaN?

NaN not a number.

position이 NaN이 되어버린 particle들은 사라진 것 같다. 

왜 연산 과정에서 NaN이 나오지? 

 

역시 pressure와 viscosity를 계산하는 과정에서 문제가 생겼을 것이다. 

 

    // 중력 힘 적용
    vec3 currentForce = vec3(0, -9.81 * particleMass, 0); - resultPressure + visc;
    gl_FragColor = vec4(currentForce, 1.0);

아 혹시 이 부분인가.

compute shader에서 이전 force의 결과를 토대로 update해야하는데 그러지 않았다. 

 

흠... 근데 그건 아니다.

현재 force는 이전 force의 영향을 받지 않는다. 

 

디버깅

보니 dependency 사이에 순서가 있긴 하다. 

 

디버깅으로 첫 번째 compute의 결과 force texture를 뽑아보았더니. 

이미 안에 NaN들이 들어있다. 

 

원인파악을 좀 해보자.

 

    gl_FragColor = vec4(gl_FragCoord.x, 0.0, 0.0, 1.0);

우선 position shader에서 gl_FragCoord가 float이어야만 하는 이유를 찾을려고 

그대로 출력해보았다.

 

보면 차이는 1씩 나지만, .5의 정수로 이루어져 있다.

범위는 0.5, 1.5 ... 999.5 까지

 

이걸 int로 하면 안되나? 

int로 변경하면 0, 1, ... 999 까지가 된다.

솔직히 뭐가 별반 다른가 싶기도 하다...

 

gl_FragCoord에 이렇게 0.5 부터 1단위로 들어간 이유는

fragment 중심좌표를 사용하기 때문이라고 한다. 

이것 때문일 가능성이 좀 있어 보인다...

 

int로 지정된 id를 싹 다 float형태로 바꿔보자. 

뭔가 서로 튕겨내는 힘이 생기긴 했는데,

갑자기 3x3 형태로 바뀌었고 무슨 트램펄린처럼 튀어버린다.

 

이번에는 그래도 density pressure 값에서 NaN이 나오지 않았는데 

 

position을 뽑아봤더니 2~3 프레임 정도에서 갑자기 

position이 이상하게 변했다. 

 

위 그림에서 본 것 처럼

0~99 / 400~499 / 800~899 번째의 particle 9개씩만 position이 정상으로 처리되고 

나머지는 전부 NaN이다. 

 

왜 이러지...

 

    vec3 currentForce = vec3(0, -9.81 * particleMass, 0); - resultPressure + visc;
    gl_FragColor = vec4(visc, 1.0);

compute force shader에서 pressure와 visc를 뽑아본 결과.

첫 계산부터 NaN이 나오는 것을 확인..

 

shader 순서 문제라면

만약 density pressure보다 force계산이 앞서 계산되어서 생기는 문제라면, 

정말 이 프로젝트를 접어야 될 수도 있는 큰 문제이다.

 

        shader.Dispatch(computeForceKernel, totalParticles / 256, 1, 1);
        shader.Dispatch(densityPressureKernel, totalParticles / 256, 1, 1);
        shader.Dispatch(integrateKernel, totalParticles / 256, 1, 1);

타당한 이유일 수 있는게, unity SPH 프로젝트에서도

density pressure kernel보다 compute force가 맞서 연산되면 오류가 발생한다.

 

그럼 셰이더를 나누지 않고 합치면 안될려나?

아.. 보니까 다른 particle의 density와 pressure를 사용하기 때문에 어렵다. 

 

Fluid Simulation 구현

아 꼬여서 이것저것 하고 있는데... 어라?

오 나름.. 파티클 유체 효과를 내도록 바뀌었다. 

이게 어떻게 된 거냐면,

 

precision highp float;

uniform int particleLength;
uniform float particleMass;
uniform float viscosity;
uniform float gasConstant;
uniform float restDensity;
uniform float boundDamping;
uniform float pi;
uniform vec3 boxSize;
uniform float radius;
uniform float radius2;
uniform float radius3;
uniform float radius4;
uniform float radius5;
uniform float timestep;

// uniform sampler2D positionTexture;
// uniform sampler2D velocityTexture;
// uniform sampler2D forceTexture;
// uniform sampler2D densityPressureTexture;

// 첫 번째 도함수
float SpikyKernelFirstDerivative(float distance) {
    float x = 1.0 - distance / radius;
    return -45.0 / (pi * radius4) * x * x;
}

// 두 번째 도함수
float SpikyKernelSecondDerivative(float distance) {
    float x = 1.0 - distance / radius;
    return 90.0 / (pi * radius5) * x;
}

// 그라디언트 계산 함수
vec3 SpikyKernelGradient(float distance, vec3 direction) {
    return SpikyKernelFirstDerivative(distance) * direction;
}

float StdKernel(float distanceSquared) {
  float x = 1.0f - distanceSquared / radius2;
  return 315.0 / (64.0 * pi * radius3) * x * x * x;
}

vec2 getDensityPressure(float index){
    // 파티클 위치 읽기
    vec3 origin = texture2D(positionTexture, vec2(index / float(particleLength))).xyz;
    float sum = 0.0;

    for (float i = 0.5; int(i) < particleLength; i++) {

        vec3 otherPos = texture2D(positionTexture, vec2(i / float(particleLength))).xyz;
        vec3 diff = origin - otherPos;
        float distanceSquared = dot(diff, diff);
        
        if (radius2 * 0.004 > distanceSquared * 0.004) {
            sum += StdKernel(distanceSquared * 0.004);
        }
    }

    // 파티클 밀도 및 압력 계산
    float density = sum * particleMass + 0.000001;
    float pressure = gasConstant * (density - restDensity);
    vec2 densityPressure = vec2(density, pressure);

    return densityPressure;
}

void main() {
    float id = float(gl_FragCoord.x);  // 각 픽셀은 파티클 ID와 매핑

    if (id >= float(particleLength)) {
        gl_FragColor = vec4(0.0);  // 범위 밖일 경우 아무런 작업도 하지 않음
        return;
    }

    // 파티클 위치 읽기
    vec3 origin = texture2D(positionTexture, vec2(id / float(particleLength))).xyz;
    float sum = 0.0;

    for (float i = 0.5; int(i) < particleLength; i++) {

        vec3 otherPos = texture2D(positionTexture, vec2(i / float(particleLength))).xyz;
        vec3 diff = origin - otherPos;
        float distanceSquared = dot(diff, diff);
        
        if (radius2 * 0.004 > distanceSquared * 0.004) {
            sum += StdKernel(distanceSquared * 0.004);
        }
    }

    // 파티클 밀도 및 압력 계산
    float density = sum * particleMass + 0.000001;
    float pressure = gasConstant * (density - restDensity);

    // ---------------

    // float density = texture2D(densityPressureTexture, vec2(id / float(particleLength))).x;
    // float pressure = texture2D(densityPressureTexture, vec2(id / float(particleLength))).y;
    float density2 = density * density;
    float mass2 = particleMass * particleMass;
    vec3 resultPressure = vec3(0.0);
    vec3 visc = vec3(0.0);

    // if (density != density || pressure != pressure){
    //     gl_FragColor = vec4(0.0);
    //     return;
    // }

    for (float i = 0.5; int(i) < particleLength; i++) {
        vec3 otherPosition = texture2D(positionTexture, vec2(i / float(particleLength))).xyz;
        
        // 자기 자신과 비교하지 않기
        if (id == i) continue;

        float dist = distance(otherPosition, origin);
        if (dist < radius * 2.0) {
            // 압력 계산
            vec3 pressureGradientDirection = normalize(origin - otherPosition);

            vec3 pressureContribution = mass2 * SpikyKernelGradient(dist, pressureGradientDirection);
            // float otherDensity = texture2D(densityPressureTexture, vec2(i / float(particleLength))).x;
            // float otherPressure = texture2D(densityPressureTexture, vec2(i / float(particleLength))).y;
            float otherDensity = density;
            float otherPressure = pressure;
            pressureContribution *= (pressure / density2) + otherPressure / (otherDensity * otherDensity);

            // 점성 계산
            vec3 velocity = texture2D(velocityTexture, vec2(id / float(particleLength))).xyz;
            vec3 otherVelocity = texture2D(velocityTexture, vec2(i / float(particleLength))).xyz;
            vec3 viscosityContribution = viscosity * mass2 * (otherVelocity - velocity) / otherDensity;
            viscosityContribution *= SpikyKernelSecondDerivative(dist);
            
            resultPressure += pressureContribution;
            visc += viscosityContribution;
        }
    }

    // 중력 힘 적용
    vec3 currentForce = vec3(0, -9.81 * particleMass, 0) - resultPressure + visc;
    gl_FragColor = vec4(currentForce, 1.0);
}

이렇게 결국은 합쳐 버렸다.

 

            vec3 pressureGradientDirection = normalize(origin - otherPosition);

            vec3 pressureContribution = mass2 * SpikyKernelGradient(dist, pressureGradientDirection);
            // float otherDensity = texture2D(densityPressureTexture, vec2(i / float(particleLength))).x;
            // float otherPressure = texture2D(densityPressureTexture, vec2(i / float(particleLength))).y;
            float otherDensity = density;
            float otherPressure = pressure;
            pressureContribution *= (pressure / density2) + otherPressure / (otherDensity * otherDensity);

그 와중에 다른 particle의 density와 pressure가 필요한 부분은

대체법으로 현재 파티클의 density와 pressure로 넣었다.

 

즉 오류는 해결했지만 완벽한 로직은 아닌 셈이다. 

 

(부록) WebGL 과부하 

render-setting.ts:5 THREE.WebGLRenderer: A WebGL context could not be created. Reason: Web page caused context loss and was blocked

그렇다고 너무 shader에서 과하게 연산을 하게 되면 이렇게 잠시동안 webGL이 과부하가 걸리면서 렌더링이 안된다.