   CFD - Vorticity Confinement
  I just tried to add back any curls, that are lost during my numerical fluid dynamics simulation. I use a 3D grid and the stable fluids approach from Stam to solve the incompressible Navier-Stokes equations. I assume the stable fluid convection method damps out those interesting curls.So I tried vorticity confinement, but it didn''t give the results I expected. Curls behave rather chaotic.I computed for each cell: W = Nabla cross U   // U = velocity vectorN = Nabla |W| / |Nabla |W||   // W = Vorticity, F = N cross W // N = Vector pointing to vortex-source// F = vorticity confinement forceF is then Euler-integrated and added as source term to the velocity.Anyone tried this? Did I miss anything?Cheers,   Oliver 
  I see nothing wrong in your formulas. Perhaps one source of chaotism might be errors in the cross product formulas? I had some in the first place.Another (obvious) cause of chaotism is that you''re not (at least in your post) scaling your confinement force. It should always be scaled with some confinement scalar - otherwise the velocity field diverges (Stable Fluids won''t help here as more and more mad turbulence is added every frame). If it diverges, use a smaller scalar.Another idea could be to negate the confinement force. It might have the wrong sign. (If this is the case, it should be seen as increased dissipation anyway.)So all in all, I bet on scaling the confinement force and possibly negating it  - Mikko Kauppila ;  ..Yes, the cross product formulas might be faulty.Here is how I implemented the vorticity confinement.vec3 getVortex(int i, int j, int k) {  float DwDy = (W[i,j+1,k]-W[i,j-1,k])*0.5f;  // central difference, with cell length of 1  float DvDz = (V[i,j,k+1]-V[i,j,k-1])*0.5f;  float DuDz = (U[i,j,k+1]-U[i,j,k-1])*0.5f;  float DwDx = (W[i+1,j,k]-W[i-1,j,k])*0.5f;  float DvDx = (V[i+1,j,k]-V[i-1,j,k])*0.5f;  float DuDy = (U[i,j+1,k]-U[i,j-1,k])*0.5f;  return vec3(DwDy-DvDz,DuDz-DwDx,DvDx-DuDy);}// compute vorticity confinement force (for every cell)for every cell (i,j,k) {  // float nrm2(const vec3& v)  computes the vectorlength as : return sqrt(v.x*v.x + v.y*v.y + v.z*v.z);  float DoDx = (nrm2(getVortex(i+1,j,k)) - nrm2(getVortex(i-1,j,k))) * 0.5f;    float DoDy = (nrm2(getVortex(i,j+1,k)) - nrm2(getVortex(i,j-1,k))) * 0.5f;  float DoDz = (nrm2(getVortex(i,j,k+1)) - nrm2(getVortex(i,j,k-1))) * 0.5f;  float len = nrm2(vec3(DoDx,DoDy,DoDz));  DoDx/=len;  DoDy/=len;  DoDz/=len;  vec3 v=getVortex(i,j,k);  curl_x[i,j,k] = (DoDy*v.z - v.y*DoDz) * 60.0f; // curl_x is the force field''s x-component  curl_y[i,j,k] = (DoDz*v.x - v.z*DoDx) * 60.0f;  curl_z[i,j,k] = (DoDx*v.y - v.x*DoDy) * 60.0f; }Any errors?...Negating the function did not change the chaotic behavior.Cheers,  Oliver ;  Mikko, Please tell me, if you implemented that, and had successful (or non-successful) results.I am hanging around with this for a while, and I need these curls.Cheers, Oliver ;  I have implemented vorticity confinement successfully in both 2d and 3d, with cell-centered and staggered grids.My idea number one is to not use the constant 60.0f to scale the confinement force, but try a smaller number. If it still diverges, try a much smaller number. Simple.I checked your code, and I see nothing wrong in there. Just change the value 60.0f to a lower one.If you get artifical dissipation instead of artificial turbulence (AFTER lowering the constant 60.0f), try negating the confinement force. - Mikko Kauppila   [edited by - uutee on August 2, 2003 6:20:12 AM];  hi Mikko.tell me, how large is your 3D grid? did vorticity confinement produced good results? i mean, im sure its doing its job, but can help with creating curls in small grids, like 32*32 or 64*64? ;  pickups: vorticity confinement definitely helps in (practical sized) 32x32/64x64 grids. On extremely coarse grids (8x8 for example) it doesn''t help, but so does nothing else really  My 3D grids are at the moment 20x20x20 and 20x24x20.About your projection step problem (Ap=b). p (the pressure) and b (gradient of velocity) are the "stretched out" (as you called it) vectors in which each node corresponds to one grid point.In the matrix A, each row encaptures one equation of the formP(x+1,y,z) + P(x,y+1,z) + P(x,y,z+1) - 6*P(x,y,z) + P(x-1,y,z) + P(x,y-1,z) + P(x,y,z-1) = (h/2)*[ u(x+1,y,z) - u(x-1,y,z) + v(x,y+1,z) - v(x,y-1,z) + w(x,y,z+1) - w(x,y,z-1)]to a single grid point (x,y,z). On that row in the matrix you store the coefficients of pressures (as you can see, those are only 1 on off-diagonal entries and -6 diagonal column entries).I don''t know how it can be explained with more intuition. For me, the talks of Jos Stam http://www.dgp.toronto.edu/people/stam/reality/Talks/index.html on Stable Fluids and Visual Simulation of Smoke gave good insights into forming the linear system.You don''t have to understand the meaning of a band of a matrix for implementing the method. If you write some practice matrices on a paper (helped for me), you should see that the matrix has clear diagonally oriented lines, which are called bands. As a little tip for the future I can reveal that there exists a great spot for a speedup in exploiting the banded structure... - Mikko Kauppila ;  Ok. thnks. i didnt know the links for the talks, ill go over it. and just try to understand. thanks. ;  Allright! i think i got it. each row in the matrix is the equation for this cell and we have N*N*N rows. in ech column there is 0 or 1 depends on whether this cell is adjacent to the center grid. and the diaognal is always -6. right?so this matrix can be precalculated, which is a good thing. ;  Hi Mikka,I guess that was it. I now took a vorticity confinement factor of 15, which gives decent results. :-)Which factor do you take?I just have an asymmetrie left, where the confinement force tends to strengthen curls on one side of my computing grid.The odd thing is, when I only negate the y component of my confinement force, this asymmetrie get''s more or less away.(I have an inflow-BC at the bottom)But I have no reason to do so... I mean this should not be.Anyway, thx, that you opened my eyes on this.Cheers,  Oliver 