I am trying to implement a 3D fluid simulator based on SDK’s fluidGL, but i don’t know how to translate the projection part to 3D, here is the code
// This method performs velocity diffusion and forces mass conservation
// in the frequency domain. The inputs 'vx' and 'vy' are complex-valued
// arrays holding the Fourier coefficients of the velocity field in
// X and Y. Diffusion in this space takes a simple form described as:
// v(k,t) = v(k,t) / (1 + visc * dt * k^2), where visc is the viscosity,
// and k is the wavenumber. The projection step forces the Fourier
// velocity vectors to be orthogonal to the vectors for each
// wavenumber: v(k,t) = v(k,t) - ((k dot v(k,t) * k) / k^2.
__global__ void
diffuseProject_k(cData *vx, cData *vy, int dx, int dy, float dt,
float visc, int lb) {
int gtidx = blockIdx.x * blockDim.x + threadIdx.x;
int gtidy = blockIdx.y * (lb * blockDim.y) + threadIdx.y * lb;
int p;
cData xterm, yterm;
// gtidx is the domain location in x for this thread
if (gtidx < dx) {
for (p = 0; p < lb; p++) {
// fi is the domain location in y for this thread
int fi = gtidy + p;
if (fi < dy) {
int fj = fi * dx + gtidx;
xterm = vx[fj];
yterm = vy[fj];
// Compute the index of the wavenumber based on the
// data order produced by a standard NN FFT.
int iix = gtidx;
int iiy = (fi>dy/2)?(fi-(dy)):fi;
// Velocity diffusion
float kk = (float)(iix * iix + iiy * iiy); // k^2
float diff = 1.f / (1.f + visc * dt * kk);
xterm.x *= diff; xterm.y *= diff;
yterm.x *= diff; yterm.y *= diff;
// Velocity projection
if (kk > 0.f) {
float rkk = 1.f / kk;
// Real portion of velocity projection
float rkp = (iix * xterm.x + iiy * yterm.x);
// Imaginary portion of velocity projection
float ikp = (iix * xterm.y + iiy * yterm.y);
xterm.x -= rkk * rkp * iix;
xterm.y -= rkk * ikp * iix;
yterm.x -= rkk * rkp * iiy;
yterm.y -= rkk * ikp * iiy;
}
vx[fj] = xterm;
vy[fj] = yterm;
}
}
}
}
Anyone give me some ideas?
I don’t know the wavenumber in Z. And i think i need to forces the Fourier velocity vectors to be orthogonal to the plane for each wavenumber, but i have no idea how to start.