I am making a 3D fluid simulation in JS/WebGL based on Jos Stam's papers. I believe everything is working except for either the boundary conditions, advection, or something above my expertise.
My simulation diffuses, responds to changes in velocity, but continues beyond its container after reaching the boundary.
Transparent parts are empty cells and container boundary in those pictures.
I am concerned there is something fundamentally wrong with what I am doing. I want to see fluid spawn, move, and interact with its container. Maybe the nature of the simulation or my mesh is flawed. Please let me know if I can provide any more details.
Below are my boundary, diffusion, and advection routines. My mesh is rectangular with dimensions N + extra * N * N instead of (N + 1) ^3 like in Stam's papers. Thus, the boundaries are at:
0, j, k // horizontal parallel to page
N + extra - 1, j, k
i, 0, k // vertical parallel to page
i, N - 1, k
i, j, 0 // perpendicular to page
i, j, N - 1
boundary(value, flag) {
for (let i = 0; i < N; i++) {
for (let j = 0; j < N; j++) {
// yz parallel horizontal plane
if (flag == 1) {
value[idx(0, i, j)] = -value[idx(1, i, j)];
value[idx(N + extra - 1, i, j)] = -value[idx(N + extra - 2, i, j)];
} else {
value[idx(0, i, j)] = value[idx(1, i, j)];
value[idx(N + extra - 1, i, j)] = value[idx(N + extra - 2, i, j)];
}
}
}
for (let i = 0; i < N + extra; i++) {
for (let j = 0; j < N; j++) {
// xz parallel vertical plane
if (flag == 2) {
value[idx(i, 0, j)] = -value[idx(i, 1, j)];
value[idx(i, N - 1, j)] = -value[idx(i, N - 2, j)];
} else {
value[idx(i, 0, j)] = value[idx(i, 1, j)];
value[idx(i, N - 1, j)] = value[idx(i, N - 2, j)];
}
// xy perpendicular plane
if (flag == 3) {
value[idx(i, j, 0)] = -value[idx(i, j, 1)];
value[idx(i, j, N - 1)] = -value[idx(i, j, N - 2)];
} else {
value[idx(i, j, 0)] = value[idx(i, j, 1)];
value[idx(i, j, N - 1)] = value[idx(i, j, N - 2)];
}
}
}
value[idx(0, 0, 0)] = .33 * (value[idx(1, 0, 0)] + value[idx(0, 1, 0)] + value[idx(0, 0, 1)]);
value[idx(0, N - 1, 0)] = .33 * (value[idx(1, N - 1, 0)] + value[idx(0, N - 2, 0)] + value[idx(0, N - 1, 1)]);
value[idx(N + extra - 1, 0, 0)] = .33 * (value[idx(N + extra - 2, 0, 0)] + value[idx(N + extra - 1, 1, 0)] + value[idx(N + extra - 1, 0, 1)]);
value[idx(N + extra - 1, N - 1, 0)] = .33 * (value[idx(N + extra - 2, N - 1, 0)] + value[idx(N + extra - 1, N - 2, 0)] + value[idx(N + extra - 1, N - 1, 1)]);
value[idx(0, 0, N - 1)] = .33 * (value[idx(1, 0, N - 1)] + value[idx(0, 1, N - 1)] + value[idx(0, 0, N - 2)]);
value[idx(0, N - 1, N - 1)] = .33 * (value[idx(1, N - 1, N - 1)] + value[idx(0, N - 2, N - 1)] + value[idx(0, N - 1, N - 2)]);
value[idx(N + extra - 1, 0, N - 1)] = .33 * (value[idx(N + extra - 2, 0, N - 1)] + value[idx(N + extra - 1, 1, N - 1)] + value[idx(N + extra - 1, 0, N - 2)]);
value[idx(N + extra - 1, N - 1, N - 1)] = .33 * (value[idx(N + extra - 2, N - 1, N - 1)] + value[idx(N + extra - 1, N - 2, N - 1)] + value[idx(N + extra - 1, N - 1, N - 2)]);
return value;
}
linear_solve(flag, value, value_old, rate, c) {
for (let k = 0; k < 10; k++) {
for (let x = 1; x < N + extra - 1; x++) {
for (let y = 1; y < N - 1; y++) {
for (let z = 1; z < N - 1; z++) {
let sum = value[idx(x - 1, y, z)] + value[idx(x + 1, y, z)] +
value[idx(x, y - 1, z)] + value[idx(x, y + 1, z)] +
value[idx(x, y, z - 1)] + value[idx(x, y, z + 1)];
value[idx(x, y, z)] = (value_old[idx(x, y, z)] + rate * sum) / c;
}
}
}
}
value = this.boundary(value, flag);
return value;
}
diffuse(flag, value, value_old, viscosity, dt) {
let rate = dt * viscosity * (N + extra) * N**2;
let c = 1 + 6 * rate;
value = this.linear_solve(flag, value, value_old, rate, c)
return value;
}
advect(flag, value, value_old, dt) {
let dtOX = dt * (N + extra);
let dtO = dt * N;
for (let i = 1; i < N + extra - 1; i++) {
for (let j = 1; j < N - 1; j++) {
for (let k = 1; k < N - 1; k++) {
let x = i - dtOX * this.u[idx(i, j, k)];
let y = j - dtO * this.v[idx(i, j, k)];
let z = k - dtO * this.w[idx(i, j, k)];
if (x < 0.5) {
x = 0.5;
}
if (y < 0.5) {
y = 0.5;
}
if (z < 0.5) {
z = 0.5;
}
if (x > N + extra - 1.5) {
x = N + extra - 1.5;
}
if (y > N - 1.5) {
y = N - 1.5;
}
if (z > N - 1.5) {
z = N - 1.5;
}
let i0 = Math.floor(x);
let j0 = Math.floor(y);
let k0 = Math.floor(z);
let i1 = i0 + 1;
let j1 = j0 + 1;
let k1 = k0 + 1;
let s1 = x - i0;
let t1 = y - j0;
let u1 = z - k0;
let s0 = 1 - s1;
let t0 = 1 - t1;
let u0 = 1 - u1;
value[idx(i, j, k)] =
s0 * (
t0 * u0 * value_old[idx(i0, j0, k0)] +
t1 * u0 * value_old[idx(i0, j1, k0)] +
t0 * u1 * value_old[idx(i0, j0, k1)] +
t1 * u1 * value_old[idx(i0, j1, k1)]
) +
s1 * (
t0 * u0 * value_old[idx(i1, j0, k0)] +
t1 * u0 * value_old[idx(i1, j1, k0)] +
t0 * u1 * value_old[idx(i1, j0, k1)] +
t1 * u1 * value_old[idx(i1, j1, k1)]
);
}
}
}
value = this.boundary(value, flag);
return value;
}