我从这个博客中遇到了格子玻尔兹曼流体求解器的实现，该博客介绍了它的实现。我决定用 webgl 把它翻译成 ShaderToy。着色器玩具的局限性迫使我使用多个纹理帧，而不是计算平衡所需的每个步骤。
我的算法基本上是这样的：
从纹理 B、C 和 D 计算纹理 A 中每个单元的体积速度和密度，这些纹理 B、C 和 D 包含 9 个方向分量（B 中的每个 vec4 中的 NW、N、NE、W，每个 vec4 中的 E、SW、S、SE C，在 D 中的每个 vec4 中居中）。
计算完这些分量后，重新计算每个纹理帧 B、C、D 的每个必要单元的流/平流，并从 A 中获取密度和速度分量来计算平衡。
将每个方向的最终值设置为
new_direction - (new_direction - new_direction_equilibrium), ie(north_west - (north_west - north_west_eq))`。这与参考代码没有任何不同。
另外我环绕坐标，所以没有边界条件以避免处理边界逻辑，并且用户单击以在晶格中引起干扰，这将“不移动”方向值设置为一个数字。
然而，在我的代码中，我最终得到了很多 NaN（这里用白色着色，红色代表密度），我不知道为什么。我在代码中的位置放置了安全卫士，以避免零密度导致问题，但这似乎没有任何作用。
您可以在 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？有没有办法阻止他们？