PDA

View Full Version : A bit OT: a discrete form of grad(divU).



okapota
04-27-2003, 10:33 AM
hello. i just read a paper about modelling clouds, and im implementing it's algorithm. in order to do that, i need the discrete form of W = grad(divU), where U is the velocity vector. so, in the paper, the x component of W is presented, and it is noted that the y and z are computed in a similiar way, so i sat and wrote the equations, and ended up with a slightly different solution, so maybe im wrong, or maybe there is a mistake in the paper. so if you can take a look at that-

my solution -
Wx = 1/4[Ux(x+2,y,z)+Ux(x-2,y,z)-2Ux(x,y,z)]
+ 1/4[Uy(x+1,y+1,z)+Uy(x-1,y-1,z)-Uy(x-1,y+1,z)-Uy(x+1,y-1,z)+Uz(x+1,y,z+1)+Uz(x-1,y,z-1)-Uz(x-1,y,z+1)-Uz(x+1,y,z-1)].

the papers solution, the second term is identical-

Wx = 1/2[Ux(x+1,y,z)+Ux(x-1,y,z)-2Ux(x,y,z)]
+ 1/4[Uy(x+1,y+1,z)+Uy(x-1,y-1,z)-Uy(x-1,y+1,z)-Uy(x+1,y-1,z)+Uz(x+1,y,z+1)+Uz(x-1,y,z-1)-Uz(x-1,y,z+1)-Uz(x+1,y,z-1)].

my question is who is right and what should be Wy,Wz?
i will appriciate any answer, thnx.

graphicsMan
04-27-2003, 01:15 PM
Is grad(divU) the divergence of U?

If so, then the grad dot U =

d(Ux)/dx + d(Uy)/dy + d(Uz)/dz

which ends up being estimated by

(1/2)*(Ux(x+1,y,z)-Ux(x-1,y,z)) etc... (by central differencing)

Brian

okapota
04-27-2003, 11:54 PM
its the gradient of the divergence of U.
i know, what you wrote is what used to start with.

graphicsMan
04-28-2003, 07:47 AM
Divergence is a scalar quantity... so the gradient should be a single-valued function. I don't get how you are getting Wx,Wy, and Wz...

Maybe I'm missing something?

bakery2k
04-28-2003, 08:55 AM
He wants the DIVERGENCE OF the gradient of U. This is a vector.

okapota
04-28-2003, 09:12 AM
bakery is right.

graphicsMan
04-28-2003, 07:23 PM
Okay, the answer I get goes:

1/4[Ux(x+1,y,z)+Ux(x-1,y,z)-2Ux(x,y,z)]...

This is because the 1/4 comes from 1/(deltaX^2) where your deltaX is 2. With your formula, you should get 1/16...

Brian

okapota
04-28-2003, 11:20 PM
ok, so why in the paper it is 1/2?
and what should y/z look like?

graphicsMan
05-01-2003, 08:51 PM
Originally posted by okapota:
ok, so why in the paper it is 1/2?
and what should y/z look like?

Good question! I don't know why the paper has 1/2. I believe that y/z should be analogous.

okapota
05-01-2003, 11:05 PM
by analogous you mean like that -

Wy = 1/2[Uy(x,y+1,z)+Uy(x,y-1,z)-2Uy(x,y,z)]
+ 1/4[Ux(x+1,y+1,z)+Ux(x-1,y-1,z)-Ux(x-1,y+1,z)-Ux(x+1,y-1,z)+Uz(x,y+1,z+1)+Uz(x,y-1,z-1)-Uz(x,y-1,z+1)-Uz(x,y+1,z-1)]/.

?

graphicsMan
05-02-2003, 09:59 AM
Yes... except that I think the 1/2 should be a 1/4?

Basically, you're just taking a different partial derivative, so it should be the same computation just with respect to a different variable.

tie
05-02-2003, 02:39 PM
The paper is incorrect. They copied their formula off Yanagita & Kaneko's paper, which had the incorrect formula.

PS You might find vorticity confinement and stronger incompressibility helpful for implementing clouds.

[This message has been edited by tie (edited 05-02-2003).]

okapota
05-02-2003, 11:09 PM
by saying stronger incompressibility you mean doing viscosity and pressure seperatly, each with a poission iterative solver?

okapota
05-04-2003, 09:48 AM
ok, nevermind. what about this, in a slide from gdc 2003 it is said that the decrete version of a poission equation in 2D is

P'i,j = 0.25*(Pi+1,j + Pi-1,j + Pi,j+1 + Pi,j-1 - Delta^2(DivU).

i understand this, but i wont to implement such a solver for a 3d grid, and i dont know how to derive the formula, anybody knows?