Add new method IVOCK about fluid-engine-dev

[zhxx1987]

@igbondar @doyubkim
here are some helpful instructions for ivock implementation:

1. the monopole tree code is not necessary, it is just used to compute the boundary values for the stream function poisson problem so that a Multigrid solver then performs more accurate asif the in an open-space. you can just use 0 values instead. while a slightly better solution could be to simply extend the solver domain a little bit when solving the vorticity stream function, with 0 boundary values.
2. the vorticity stretching is important, and we approximated it from a geometric point of view:
   which is, when you look at a conceptually "vortex tube", the velocity at its two ends are usually different, as a result, the vortex tube becomes longer or shorter by the velocity advection,
   now:
   a). circulation preservation states that the integral of vorticity flux across the cut plane of this infinitesimal tube stay unchanged,
   b) incompressibility states the volume of the tube unchanged,
   so, the cut area of the tube became smaller when it is stretched by the velocity field, which requires larger voriticity flux across the cut plan to maintain circulation.
   hence, the vorticity(circulation per unit area of this vortex tube), becomes larger when the tube is stretched and smaller when the tube shrinks, when we approximate the vortex tube by vortex segments, this means the vorticity represented by a vortex segment enlarges at the same rate the segment stretches, which is the new_length/old_length. and vice versa.
   ||w||2 is the 2-norm, or say, length(w) of the vorticity, and you can verify the fact that we were in fact computing the length ratio before and after vorticity advection. in a backward way.

[doyubkim]

Thanks for the suggestion! Yes, this paper shows very nice results and would be nice to have for Jet codebase as well. I will add this to my feature wishlist.

[zhxx1987]

feel free to let me know regarding algorithm details and implementation

[doyubkim]

Thanks, @zhxx1987 ! I haven't started developing IVOCK, yet, but will definitely ask you if I face any problems during the implementation.

[igbondar]

Hi Doyub! What is your plan on implementing this?
The part with getting velocity from vorticity bit hard to transfer to code(.

[doyubkim]

Hi @igbondar!

At the moment, my general plan for Jet Framework is to

1. Add/fix minor improvements/issues
2. Consolidate existing branches for v1.4 and
3. Refactor the API and some of the underlying code structure for v2.0

This means I won't be working on adding "features" like new solvers in the time frame above. If anybody can take a ticket like this one, I am more than happy to help.

Regarding your issue:

The part with getting velocity from vorticity bit hard to transfer to code(.

Feel free to ask me about the code. We can figure out what are needed to support your work.

[igbondar]

Hi Doyub! Thanks for answer.

For now my main interest to test IVOCK in Houdini for real-world scenarios. Especially since it can preserve shape stability overtime, compared to existing vortex confinement.

As I understand to obtain u from v I need to almost re-create MultiGrid(same as in source code to IVOCK pdf) with monopole tree as described. I will try to readopt the code, that they have shared but for sure it is much better to understand all theory they mentioned in explanation.

On thing I don't get in ivock pdf - after equation 11 in text E=h/||w||2 - what 2 meaning in this case?

And thanks for your book and source code really helpful!

[doyubkim]

For this:

On thing I don't get in ivock pdf - after equation 11 in text E=h/||w||2 - what 2 meaning in this case?

I think @zhxx1987 can chime in to share some insights.

[igbondar]

Hi @zhxx1987!
Thanks for help. I tried implement 2d case of IVOCK since it simpler and doesn't require vortex stretching.
The problem it is quickly become unstable even with a small increasing of correction values.
In picture - from left to right - stable fluids, ivock with 0.015 multiplier, ivock with 0.018 multiplier.
I think I did somewhere logic mistake or didn't use some clamping. If you have some time to answer I have couple questions.
I follow your pdf guide for 2d case w=(0,0,value):

1. Get vorticity from vel
2. advect it (BFECC - advection method Trace)
3. advect velocity (BFECC - advection method Trace)
4. get new vorticity from adv vel
5. get vort difference
6. For solving poisson equation I'm using this discretized formula
   

Xnn - sampled values of vorticity correction near examined cell - for 2d case (x-1,y), (x+1,y), (x,y-1),(x,y+1)
denominator b =4, a = -pow(2,grid_cell_size), numerator b= divergence of vorticity correction.
7. get curl from stream function and add it to velocity

1. Can you suggest where could be possible error and how to control simulation so it isn't explode over time and be more predictable?
2. In your C++ code you also using decay_vortices - what purpose of it? Or it is just for pyro model - I see only burn sampling.

[zhxx1987]

Hi @zhxx1987!
Thanks for help. I tried implement 2d case of IVOCK since it simpler and doesn't require vortex stretching.
The problem it is quickly become unstable even with a small increasing of correction values.
In picture - from left to right - stable fluids, ivock with 0.015 multiplier, ivock with 0.018 multiplier.
I think I did somewhere logic mistake or didn't use some clamping. If you have some time to answer I have couple questions.
I follow your pdf guide for 2d case w=(0,0,value):

1. Get vorticity from vel
2. advect it (BFECC - advection method Trace)
3. advect velocity (BFECC - advection method Trace)
4. get new vorticity from adv vel
5. get vort difference
6. For solving poisson equation I'm using this discretized formula
   

Xnn - sampled values of vorticity correction near examined cell - for 2d case (x-1,y), (x+1,y), (x,y-1),(x,y+1)
denominator b =4, a = -pow(2,grid_cell_size), numerator b= divergence of vorticity correction.
7. get curl from stream function and add it to velocity

1. Can you suggest where could be possible error and how to control simulation so it isn't explode over time and be more predictable?
2. In your C++ code you also using decay_vortices - what purpose of it? Or it is just for pyro model - I see only burn sampling.

Hi igbondar,
here are couple things to check:

1. BFECC itself is not an unconditionally stable advection scheme, it will introduce new extrema after advection. so, make sure you do the extrema clamping as in many real fluid solvers:
   let's say the grid value at the grid point x of any field is phi_new after bfecc advection.
   do a back trace of x using the divergence-free velocity field again and sample the 4 corners of the un-advected field.
   clamp phi_new by the max and min of the 4 corner values.
2. there would be un-physical vorticity near the domain boundary and object boundary, make-sure to mask them out. after you computed delta_omega, you can zerofy those values close to domain or object boundaries.i mask out 10 cells away from the boundaries.
3. the way you are solving the stream function equation looks like Jacobi iteration to me. which doesn't give best visual quality at all, for the vorticity stream function. you might be able to find 2D multigrid solvers from one of my repo and please use it instead.
4. I believe the figures in my paper for 2D simulations was obtained with Semi-Lagrangian advectors for both velocity and vorticity fields. So it could be much less troublesome if you use this most robust advector instead.
5. numerator b= divergence of vorticity correction. this is also incorrect, the 2D stream function for vorticity is simply
   Laplace p = - omega
   with omega = dv/dx - du/dy
   and you shall recover (u,v) = (dp/dy, -dp/dx)

[zhxx1987]

@igbondar after you done my check list, you shall be able to add back 1.0 of the missing vortices with no problem in 2D.
although I used 0.95 in my 3D cases due to some extra errors introduced by the vortex stretching.

[zhxx1987]

for using multigrid solvers in 2D, please check the code at https://github.com/ziyinq/Bimocq/blob/master/src/bimocq2D/BimocqSolver2D.cpp
"buildMultiGrid" and "projection". for vorticity stream functions, PURE_NEUMANN = false.