Lattice Boltzmann WebGL 翻译结果到处都是 NaN

Question

我遇到了 this implementation of Lattice Boltzmann fluid solver from this blog，其中介绍了它的实现。我决定用 webgl 将其转换为 ShaderToy。着色器玩具的局限性迫使我使用多个纹理帧，而不是计算平衡所需的每个步骤。

我的算法基本上是这样的：

从纹理 B、C 和 D 计算纹理 A 中每个单元格的体积速度和密度，其中纹理 B、C 和 D 包含 9 个方向分量（B、E 中每个 vec4 中的 NW、N、NE、W ,SW,S,SE 在 C 中的每个 vec4 中，在 D 中的每个 vec4 中居中）。
计算完这些分量后，为每个纹理帧B、C、D的每个必要单元格重新计算stream/advection，并从A中取出密度和速度分量来计算平衡.
将每个方向的最终值设置为 new_direction - (new_direction - new_direction_equilibrium), ie(north_west - (north_west - north_west_eq))`。这与参考代码没有任何不同。

另外我把坐标环绕起来，所以没有边界条件，以避免处理边界逻辑，用户点击会引起格子的扰动，从而将"not moving"方向值设置为数字.

然而，在我的代码中，我最终得到了很多 NaN（此处为白色，红色代表密度），我不确定为什么。我在我的代码中的位置放置了安全防护以避免零密度导致问题，但这似乎没有做任何事情。

可以test this out on shadertoy，但是我的代码如下：

//COMMON functions
const int DIRECTION_COUNT = 9;
const int DIMENSION_COUNT = 2;
const float LATTICE_SPEED = 0.1;
const float TAU = 0.9;


const vec2 north_offset = vec2(0.0,1.0);
const vec2 north_west_offset = vec2(-1.0,1.0);
const vec2 north_east_offset = vec2(1.0,1.0);
const vec2 west_offset = vec2(-1.0,0.0);
const vec2 east_offset = vec2(1.0,0.0);
const vec2 south_offset = vec2(0.0,-1.0);
const vec2 south_west_offset = vec2(-1.0,-1.0);
const vec2 south_east_offset = vec2(1.0,-1.0);
const vec2 center_offset = vec2(0.0,0.0);


const vec2 offsets[DIRECTION_COUNT] = vec2[DIRECTION_COUNT](
    north_west_offset, 
    north_offset, 
    north_east_offset,
    west_offset, 
    center_offset,
    east_offset,
    south_west_offset,
    south_offset,
    south_east_offset);

const int north_west_tex_idx = 0;
const int north_tex_idx = 1;
const int north_east_tex_idx = 2;
const int west_tex_idx = 3;

const int east_tex_idx = 0;
const int south_west_tex_idx = 1;
const int south_tex_idx = 2;
const int south_east_tex_idx = 3;

const int center_tex_idx = 0;

float textureN(sampler2D NW_N_NE_W_channel, vec2 coord, vec2 resolution){
    vec2 offset_coord = coord + south_offset;
    return texture(NW_N_NE_W_channel, offset_coord/resolution)[north_tex_idx];
}

float textureNW(sampler2D NW_N_NE_W_channel, vec2 coord, vec2 resolution){
    vec2 offset_coord = coord + south_east_offset;
    return texture(NW_N_NE_W_channel, offset_coord/resolution)[north_west_tex_idx];
}

float textureNE(sampler2D NW_N_NE_W_channel, vec2 coord, vec2 resolution){
    vec2 offset_coord = coord + south_west_offset;
    return texture(NW_N_NE_W_channel, offset_coord/resolution)[north_east_tex_idx];
}

float textureW(sampler2D NW_N_NE_W_channel, vec2 coord, vec2 resolution){
    vec2 offset_coord = coord + east_offset;
    return texture(NW_N_NE_W_channel, offset_coord/resolution)[west_tex_idx];
}

float textureS(sampler2D E_SW_S_SE_channel, vec2 coord, vec2 resolution){
    vec2 offset_coord = coord + north_offset;
    return texture(E_SW_S_SE_channel, offset_coord/resolution)[south_tex_idx];
}

float textureSW(sampler2D E_SW_S_SE_channel, vec2 coord, vec2 resolution){
    vec2 offset_coord = coord + north_east_offset;
    return texture(E_SW_S_SE_channel, offset_coord/resolution)[south_west_tex_idx];
}

float textureSE(sampler2D E_SW_S_SE_channel, vec2 coord, vec2 resolution){
    vec2 offset_coord = coord + north_west_offset;
    return texture(E_SW_S_SE_channel, offset_coord/resolution)[south_east_tex_idx];
}

float textureE(sampler2D E_SW_S_SE_channel, vec2 coord, vec2 resolution){
    vec2 offset_coord = coord + west_offset;
    return texture(E_SW_S_SE_channel, offset_coord/resolution)[east_tex_idx];
}

float textureC(sampler2D C_channel, vec2 coord, vec2 resolution){
    vec2 offset_coord = coord + center_offset;
    return texture(C_channel, offset_coord/resolution)[center_tex_idx];
}


float calc_equilibrium(const in float density, 
                 const in vec2 velocity, 
                 const in ivec2 ij) {

    int i = ij.x;
    int j = ij.y;
    // u . u
    float velmag = dot(velocity, velocity);
    // Compute the weight.
    float weight;
    if(i == 0 && j == 0) {
        weight = 4.0 / 9.0;
    } else if(i == 0 || j == 0) {
        weight = 1.0 / 9.0;
    } else {
        weight = 1.0 / 36.0;
    }

    // e_i . u
    float dotprod = float(i) * velocity.x + float(j) * velocity.y;

    float sum = 1.0;
    sum += (3.0 / LATTICE_SPEED) * dotprod;
    sum += (4.5 / (LATTICE_SPEED * LATTICE_SPEED)) * dotprod * dotprod;
    sum -= (1.5 / (LATTICE_SPEED * LATTICE_SPEED)) * velmag;
    if(density == 0.0){
        return 0.0;
    }
    return  weight * density * sum;
}

|

//Buffer A, takes in B, C, and D as in put in that order
float[DIRECTION_COUNT] stream_all(
    sampler2D NW_N_NE_W_channel, 
    sampler2D E_SW_S_SE_channel,
    sampler2D C_channel, 
    in vec2 ifragCoord){

    float north_west = textureNW(NW_N_NE_W_channel, ifragCoord, iResolution.xy);
    float north = textureN(NW_N_NE_W_channel, ifragCoord, iResolution.xy);
    float north_east = textureNE(NW_N_NE_W_channel, ifragCoord, iResolution.xy);
    float west = textureW(NW_N_NE_W_channel, ifragCoord, iResolution.xy);

    float east = textureE(E_SW_S_SE_channel, ifragCoord, iResolution.xy);
    float south_west = textureSW(E_SW_S_SE_channel, ifragCoord, iResolution.xy);
    float south = textureS(E_SW_S_SE_channel, ifragCoord, iResolution.xy);
    float south_east = textureSE(E_SW_S_SE_channel, ifragCoord, iResolution.xy);

    float center = textureC(C_channel, ifragCoord, iResolution.xy);
    return float[DIRECTION_COUNT](
        north_west, north, north_east, west, center, east, south_west, south, south_east
    );

}


float calc_density(const in float new_directions[DIRECTION_COUNT]) {
    float density; 
    for(int i = 0; i < DIRECTION_COUNT; ++i){
        density += new_directions[i];
    }
    return density;
}

vec2 calc_velocity(const in float new_directions[DIRECTION_COUNT], const in float density) {

    if(density == 0.0){
        return vec2(0.0);
    }
    if(isinf(density)){
        return vec2(0.0);
    }
    // Compute target indices.
    vec2 velocity = vec2(0.0);
    for(int idx = 0; idx < DIRECTION_COUNT; ++idx){
        vec2 ij = offsets[idx];
        float i = ij.x;
        float j = ij.y;
        velocity.x += new_directions[idx] * (i);
        velocity.y += new_directions[idx] * (j);
    }

    return velocity * (LATTICE_SPEED/density);
}



void mainImage( out vec4 fragColor, in vec2 fragCoord )
{

    ivec2 ifragCoord = ivec2(fragCoord);
    float new_directions[DIRECTION_COUNT] = stream_all(iChannel0, iChannel1, iChannel2, fragCoord);
    float density = calc_density(new_directions);
    vec2 velocity = calc_velocity(new_directions, density);
    fragColor = vec4(density,velocity.x,velocity.y,0.0);
    float center = textureC(iChannel2, fragCoord, iResolution.xy);
    float debug = center;
    if(isnan(density)){
        debug = 1.0;
        fragColor.w = debug;
    }

    //fragColor = vec4(1.0);
}

|

//Buffer B, takes in B, and A in that order
void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    if(iFrame < 10){
        fragColor = vec4(0.0);
        return;
    }
    ivec2 ifragCoord = ivec2(fragCoord - 0.5);
    float north_west = textureNW(iChannel0, fragCoord, iResolution.xy);
    float north = textureN(iChannel0, fragCoord, iResolution.xy);
    float north_east = textureNE(iChannel0, fragCoord, iResolution.xy);
    float west = textureW(iChannel0, fragCoord, iResolution.xy);

    vec4 density_velocity = texelFetch(iChannel1, ifragCoord, 0);
    float density = density_velocity.x;
    vec2 velocity = density_velocity.yz;

    float north_west_eq = calc_equilibrium(density, velocity, ivec2(north_west_offset));
    float north_eq = calc_equilibrium(density, velocity, ivec2(north_offset));
    float north_east_eq = calc_equilibrium(density, velocity, ivec2(north_east_offset));
    float west_eq = calc_equilibrium(density, velocity, ivec2(west_offset));


    fragColor = vec4((north_west - (north_west - north_west_eq) / TAU),
                     (north - (north - north_eq) / TAU),
                     (north_east - (north_east - north_east_eq) / TAU),
                     (west - (west - west_eq) / TAU));
}

|

//Buffer C, takes in C and A in that order. 
void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    if(iFrame < 10){
        fragColor = vec4(0.0);
        return;
    }
    if(iFrame < 30 && fragCoord.y < -1.0){
        fragColor = vec4(10.0, 0.0,10.0,0.0);
        return;
    }
    ivec2 ifragCoord = ivec2(fragCoord - 0.5);
    float east = textureE(iChannel0, fragCoord, iResolution.xy);
    float south_west = textureSW(iChannel0, fragCoord, iResolution.xy);
    float south = textureS(iChannel0, fragCoord, iResolution.xy);
    float south_east = textureSE(iChannel0, fragCoord, iResolution.xy);

    vec4 density_velocity = texelFetch(iChannel1, ifragCoord, 0);
    float density = density_velocity.x;
    vec2 velocity = density_velocity.yz;

    float east_eq = calc_equilibrium(density, velocity, ivec2(east_offset));
    float south_west_eq = calc_equilibrium(density, velocity, ivec2(south_west_offset));
    float south_eq = calc_equilibrium(density, velocity, ivec2(south_offset));
    float south_east_eq = calc_equilibrium(density, velocity, ivec2(south_east_offset));


    fragColor = vec4((east - (east - east_eq) / TAU),
                     (south_west - (south_west - south_west_eq) / TAU),
                     (south - (south - south_eq) / TAU),
                     (south_east - (south_east - south_east_eq) / TAU));

}

|

//Buffer D takes in D and A in that order
void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    if(iFrame < 10){
        fragColor = vec4(1, 0.0,0.0,0.0);
        return;
    }
    ivec2 ifragCoord = ivec2(fragCoord - 0.5);
    float center = textureC(iChannel0, fragCoord, iResolution.xy);


    vec4 density_velocity = texelFetch(iChannel1, ifragCoord, 0);
    float density = density_velocity.x;
    vec2 velocity = density_velocity.yz;

    float center_eq = calc_equilibrium(density, velocity, ivec2(center_offset));



    fragColor = vec4((center - (center - center_eq) / TAU),
                     0.0,
                     0.0,
                     0.0);

    vec2 mouse = vec2(iMouse.zw);
    if(mouse.x > 0.0 && mouse.y > 0.0){
        vec2 current_mouse = vec2(iMouse.xy);
        if(distance(fragCoord, current_mouse) < 3.0){
            fragColor.r = vec4(10.0).r;
        }
    }
}

|

//main image output, only takes in A as an iChannel
void mainImage( out vec4 fragColor, in vec2 fragCoord )
{
    ivec2 ifragCoord = ivec2(fragCoord-0.5);
    vec4 density_velocity = texelFetch(iChannel0, ifragCoord, 0);
    float density = density_velocity.r;
    vec2 velocity = density_velocity.gb;
    float vel_length = length(velocity);
    velocity = normalize(velocity);

    //Output to screen
    //fragColor = vec4(abs(velocity),density/100.0,vel_length/100.0);
    //fragColor = vec4(abs(velocity),0.0,1.0);
    fragColor = vec4(density/10.0,0.0,0.0,1.0);
    //
    if(density_velocity.w == 1.0){
        fragColor = vec4(1.0);
    }
}

我做错了什么导致所有这些 Nans？有办法阻止他们吗？

Answer 1

从 calc_equilibrium 钳制 return 值应该可以避免白色 NaN 溢出。

return clamp(weight * density * sum, -1000.0, 1000.0);

防止 red/black 噪音泛滥似乎并不那么简单。对于按住鼠标按钮时发生的每一帧，都会向系统添加大量能量，并且在某些时候它必然会沸腾。

Lattice Boltzmann WebGL 翻译结果到处都是 NaN

Lattice Boltzman WebGL translation results in NaNs everywhere

glsl

webgl

scientific-computing

fluid-dynamics