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

OpenGL OpenGL: Advanced Coding
1

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.

2

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

3

its the gradient of the divergence of U.
i know, what you wrote is what used to start with.

4

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?

5

He wants the DIVERGENCE OF the gradient of U. This is a vector.

6

bakery is right.

7

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

8

ok, so why in the paper it is 1/2?
and what should y/z look like?

9

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.

10

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)]/.

?

11

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.

12

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).]

13

by saying stronger incompressibility you mean doing viscosity and pressure seperatly, each with a poission iterative solver?

14

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?