Divergence incompresible flow

[Jorge Arturo Rodriguez Horcasitas]

Hi everyone,

First of all let me say thank you to the developers and everybody involved with mantaflow creation. It is awesome! 

I'm a computer science student and I'm trying to understand some fluid dynamics concepts, in particular the fact that in an incompresible flow the divergence is 0 at all points.

So I tested this out with the "plume_2d.py" program available in the "scenes" directory. I computed the divergence after applying the solvePressure function and this is what i got:

(plume_2d.py velocity divergence map for t=50)

I strongly believe that my divergence computation function is correct. Shouldn't I see in the image above a divergence close to 0 everywhere?

Thank you in advance. 

Jorge

[Nils Thuerey]

Hi Jorge, good to hear you find mantaflow useful! Regarding your question I was wondering, how to you compute the divergence? The pressure projection minimizes the divergence given as right hand side for the Poisson problem. That should work with the accuracy given as the "cgAccuracy" parameter of the solvePressure function. If you run the RHS computation right after the pressure solve (you can pass a "rhs" grid to view it), you should see the correspondingly reduced divergence values. 

However, if you compute the divergence differently, or at other points of the simulation loop, I can imagine that you'll see larger values, such as the ones in your picture...

Cheers, Nils

[Jorge Rodríguez]

Hi Nils,  thank you for your reply.

Currently I am computing the divergence right after calling the solvePressure function. I still have not found why I obtain these results. Now I think I'm computing the divergence

in a wrong way.

This is what I'm currently doing:

1. After calling the solvePressure function I copy the velocity field to a numpy array and then I obtain the corresponding X and Y components. (npVel is a numpy array previously defined as: "npVel = np.ones( [width, height, 3], dtype=np.float32 )" to contain the velocity field where width and heigh correspond to the grid resolution).

# Pass vel to npVel array

copyGridToArrayMAC( target=npVel, source=vel )

# Divergence Computation

npU = npVel[:, :, 0]

npV = npVel[:, :, 1]

2. I call my divergence function and print the colormap:

g = divergence([npU, npV])

# Graphics

plt.pcolormesh(g)

plt.colorbar()

plt.show()

The divegence function is pretty simple, it takes advantage of the np.gradient function in Numpy. Explicitly the funciton is:

def divergence(f):

    num_dims = len(f)

    return np.ufunc.reduce(np.add, [np.gradient(f[i], axis=i) for i in range(num_dims)])

The default delta x and delta y in the np.gradient computation is 1. It uses a first order one sided forwar/backward difference approximation at the boundaries and a second order accurate central difference in the rest of the points.

I will gladly explain in more details any of the previous points.

Finally, can you tell me what RHS stands for? I habe been searching in the web with no success.

Really grateful 

Jorge

[Nils Thuerey]

RHS stands for "right hand side" of the pressure system of equations. Sorry, it's a common abbreviation for iterative solvers, but not that clear in the physical context. The RHS in this case contains the negative divergence, so that's what you could look at. You can create a RealGrid, e.g. "divergence" and then call "solvePressure(... , retRhs=divergence)" in your python script. Then the grid should contain the divergence values the pressure field corrects.

Cheers, Nils

[Jorge Rodríguez]

Thanks Nils!

You are right , there must be something wrong with the way I was computing the divergence. Here is an image with both results:

On the right is the negative retRhs.