Eulerian fluid simulations simulate the flow of fluids by tracking fluid velocity and density over a set of individual (discreet) evenly spaced grid locations. One downside to this approach is that the finer details in the fluid can be smoothed out, so you lose those little swirls and vortices.

A simple fix for this is to add Vorticity Confinement. If you read the Wikipedia page on Vorticity Confinement you may be no wiser on what it is or how to add it into your fluid simulations.

My explanation of vorticity confinement is that it looks for curls (vortices) in the fluid and adds in velocity to help boost the swirling motion of the fluid. Adding vorticity confinement can also give more turbulent looking fluid simulations which tend to be more aesthetically pleasing in simulations (unless you are a member of team laminar flow).

The code for implementing vorticity confinement is relatively simple. For 2D I used the snippet provided by Iam0x539 in this video.

function Curl(x,y:integer):double;
begin
Curl:=xvelocity[x,y+1]-xvelocity[x,y-1] + yvelocity[x-1,y]-yvelocity[x+1, y];
end;
procedure VorticityConfinement(vorticity:double);
var dx,dy,len:double;
x,y:integer;
begin
for y:=2 to _h-3 do
begin
for x:=2 to _w-3 do
begin
dx:=abs(curl(x + 0, y - 1)) - abs(curl(x + 0, y + 1));
dy:=abs(curl(x + 1, y + 0)) - abs(curl(x - 1, y + 0));
len:=sqrt(sqr(dx)+sqr(dy))+1e-5;
dx:=vorticity/len*dx;
dy:=vorticity/len*dy;
xvelocity[x,y]:=xvelocity[x,y]+timestep*curl(x,y)*dx);
yvelocity[x,y]:=yvelocity[x,y]+timestep*curl(x,y)*dy);
end;
end;
end;

The VorticityConfinement procedure is called once per simulation step. It looks for local curl at each fluid grid point and then increases the local x and y velocities using the curl. This is what helps preserve the little vortices and helps reduce the smoothing out of the fluid.

To demonstrate how vorticity confinement changes a fluid simulation, the images within this post and the following movie add vorticity confinement to my previous Eulerian MAC Fluid Simulations code.

Eulerian MAC Fluid Simulations with Vorticity Confinement is now included in the latest version of Visions of Chaos.

He generously shares the source code to a series of programs on his Incremental Fluids GitHub that cover implementing a 2D fluid simulation. His code is based on Robert Bridson’s book, “Fluid Simulation for Computer Graphics”. I have seen that book mentioned all over the place and almost bought a copy, but reviews say it is focused more on the math (not so helpful to me) and not on the code (which I can follow much easier than math formulas).

So far, I have converted Benedikt’s first and second programs for inclusion in Visions of Chaos. Calculations at 4K resolution were originally taking up to 10 minutes per frame, but with some multi-threading and code optimizations I got it down to around 10 seconds per 4K resolution frame on a relatively modern i7 CPU.

The results so far are really nice. The resulting flows show very high detailed vortices and fluid behavior.

Here is a sample 4K resolution movie showing these fluids in motion.

Eulerian Marker-and-Cell Fluid Simulations are now available in the latest version of Visions of Chaos.

I am a huge fan of anything related to fluid simulations. This interest was once again sparked when Daniel Shiffman covered fluid simulation in his latest Coding Train video.

As always, I recommend Coding Train to any developer. Dan has a unique way of making his programming topics both interesting and entertaining.

The code converted during the Coding Train video (Mike Ash’s Fluid Simulation For Dummies) is based on Jos Stam’s Stable Fluids code.

Jos Stam

Jos Stam wrote his seminal paper Real-Time Fluid Dynamics for Games back in March 2003. I have lost count of the times I seen the paper cited or linked to. It has had a huge influence on fluid simulation.

The original source code from the paper is provided here and here.

Going back 11 years, one of my first YouTube videos was this super low resolution example of 2D fluid simulation using Stam’s methods.

3D Fluid

My main objective when revisiting Stam’s stable fluids was to get a 3D version going.

A quick Google search led me to Blain Maguire’s implementation which I was able to translate into a working 3D fluid simulation.

By default Stam’s fluid method generates a fairly smoothed fluid/smoke flow. To make things more interesting you want to add a bit of turbulence. The key term here is “vorticity confinement”. Vorticity confinemement was first described in the paper Visual Simulation of Smoke by Fedkiw et al.

I found some 3D source code for adding vorticity confinement to Stam’s stable fluids here. See the fluid.cpp file inside the fire32.tar.gz archive.

At this stage I had 3D fluid working. Now comes the fun part, how to display the fluid.

Displaying 3D Fluid

Once you get the code working, the main issue becomes how to display the fluid.

Stam’s fluid uses what is known as an Eulerian approach to fluid simulation. Rather than track individual fluid particles (as in SPH simulations) the fluid is simulated by tracking velocity and density at fixed grid based locations in 3D space. This means that as the simulation runs you have a 3D grid of fluid properties that need to be displayed.

The method I use is to hide the cells that contain a fluid velocity lower than a threshold value. 0.01 seems to work OK for my tests. Then render the cells as spheres or cubes. Color shaded based on fluid velocities.

For a fake volume rendering like approach you can render the fluid using translucent billboard quads. If you run a pre-render pass over the array and strip any cell that is surrounded by other cells it results in only having the “shell” of the fluids rendering. Rendering these shells with the billboard particles gives a nice result.

Back in the 1980’s Fred Mitchell designed a method to visualize Newtonian gravity he called the Mitchell Gravity Set or MGS. See here for his page explaining MGS.

This process of plotting dynamic systems is usually referred to as basins of attraction. See this paper and this paper for other versions of plotting gravitational basins of attraction. Similar methods are used when plotting Root-Finding Fractals and the basins of a Magnetic Pendulum.

The 2D images in this post are from when I first experimented with MGS back in 2005.

The basics of MGS is that you have a number of gravity objects (stars or gravity wells) in 2D space. Each pixel of a 2D image becomes the starting coordinate for a particle in that 2D space. The particle is attracted to the stars through Newtonian gravity and moves. You follow the path the particle takes for a number of steps and color the pixel depending on how long it took to be flung out of the simulation area.

I ended up rewriting/translating my old code into GLSL for speed. Click here to see the shader source code that generates these 2D Gravity Set images.

By playing with the various star positions and masses you can get a variety of different images.

The following image shows the current Gravity Set Options dialog in Visions of Chaos showing the various settings that can be tweaked.

One thing I have not tried yet that is mentioned in Fred’s YouTube video is that you can use negative gravity constants and negative masses.

Revisiting and Extending Into The 3rd Dimension

I was recently contacted by Fred again and he mentioned he was working on getting a 3D version of MGS going. Getting the 2D code going in 3D is simple enough. Just a few extra lines for adding the Z dimension, so I was inspired to have a go at seeing what the 3D versions could look like. The settings dialog is extended as follows.

Similar to the other 3D grid based displays like 3D Cellular Automata I have worked with, displaying a 3D grid of values can be tricky.

The following examples use a 500x500x500 resolution space.

The most obvious is a 3D Grid of little cubes. Each smaller cube is a starting point for the gravity calculations. Color the cubes depending on which gravity well they come closest to during their journey.

That gives you a general idea of what the shapes of the 3D Gravity Set will be, but the interior is completely hidden.

You can slice an eighth of the cubes away….

or even slice a half…

Another option is to only display the cubes for only half of the gravity wells….

or even for only a single well.

Then came the idea for “volumetric rendering”. I already had code for rendering transparent billboard quads so I used that. I am only rendering the edges of each gravity well area. If you render all points the image is just a white blob mess.

Another method. Export the visible cubes as a OBJ point cloud. Import the OBJ into MeshLab. Use MeshLab to generate normals for the points and then marching cubes to generate a mesh from the point cloud. Export the mesh as a PLY format file. Import the PLY into Blender and render it. These are rough results. I am sure someone who knows more about MeshLab and Blender could clean up the mesh and make a much better looking render.

Here is a sample movie of rotating gravity sets.

Magnetic Pendulum Alternative to MGS

I also tried using Magnetic Pendulum formulas in 3D to do similar plots. The formula is a 3D extension of the 2D code I used here. Otherwise the rest of the code is the same as for the 3D MGS above.

Availability

Both 2D and 3D Gravity Set Fractals are now included in the latest version of Visions of Chaos.

Another year gone. Over the past year Visions of Chaos has had many new features added, bugs fixed and loads of smaller enhancements. In this blog I managed to cover and experiment with a wide variety of topics over the last 12 months.

Post Summary

Here is a list of posts I added to this blog in 2018;

I have tried to include more useful information as this blog has gone on over the years to hopefully help others, but mainly blogging is most helpful to me. There is an old adage saying that if you cannot explain something clearly then you do not understand it well enough. Writing blog posts has helped me clarify many topics in my mind.

Another big advantage of blogging is that you get contacted by like minded people. Many of the interests I have blogged about have been greatly enhanced by someone emailing me a new idea or fix for an issue I am having.

I really do recommend that everyone blog. Everyone has a “something” they are good at or knowledgeable about. Share that knowledge. You may be surprised how much you benefit from the process.

Looking Forward

I have no intentions of stopping development on Visions of Chaos any time soon. I still have loads of ideas for new features. That also means this blog will continue to grow for the foreseeable future. Onward to 2019 and beyond.

3D Gravity simulations are something I have been interested in for many years now. Some worked, some didn’t, some were more realistic than others.

The first 3D Gravity simulation movie I still have on my YouTube channel is this one from way back in May of 2007. Low res with only a bunch of blurry objects.

Since then I have increased the details and object counts. I also started experimenting with OpenCL for big speedups that allowed many more objects to be simulated in a reasonable time frame.

Moving forward to now

For this latest post I went back and rewrote my code and the OpenCL kernel code to correctly compare every object to all other objects in the gravity calculations. The simulation is using Newton’s law of universal gravitation.

Every point mass attracts every single other point mass by a force acting along the line intersecting both points. The force is proportional to the product of the two masses and inversely proportional to the square of the distance between them.

How this simulation works

I am using software OpenGL on the CPU for all the rendering of the visuals and CPU and OpenCL for the gravity calculations. OpenCL code runs on your graphics card GPU and GPUs are great at running lots of small bits of code fast at the same time. The gravity formula maths is perfect for multi threading. Every objects velocity and acceleration can be calculated at the same time as the other objects.

The basics of using OpenCL is you fill arrays with the information you want the OpenCL code to use (for 3D gravity I am passing position, velocity, acceleration and mass of the objects), pass it to OpenCL, run the code on the GPU, and then read back the results from the GPU when it is done.

This is the current OpenCL code I am using for these latest simulations. Each of the arrays passed (posx, posy, etc) contain all the current objects, ie for a 1 million object simulation the posx array has 1,000,000 floating point values to cover every object’s X position in 3D space.

__kernel void Gravity3DKernel( __global float * posx,
__global float * posy
__global float * posz,
__global float * velx,
__global float * vely,
__global float * velz,
__global float * accx,
__global float * accy,
__global float * accz,
__global float * mass,
__global float * mingravdist)
{
int index=get_global_id(0);
float dx,dy,dz,distance,force;
float positionx=posx[index];
float positiony=posy[index];
float positionz=posz[index];
float mingravdistsqr=mingravdist[index]*mingravdist[index];
float accelerationx=0;
float accelerationy=0;
float accelerationz=0;
float thismass=mass[index];
for(int a=0; a<get_local_size(0); a++) {
if (a!=index) {
dx=posx[a]-positionx;
dy=posy[a]-positiony;
dz=posz[a]-positionz;
distance=sqrt(dx*dx+dy*dy+dz*dz);
dx=dx/distance;
dy=dy/distance;
dz=dz/distance;
//old method - all objects are assumed to have the same mass
//force=1/(distance*distance+mingravdistsqr);
//new method - allows objects to have different masses
force=(thismass*mass[a])/(distance*distance+mingravdistsqr);
accelerationx+=dx*force;
accelerationy+=dy*force;
accelerationz+=dz*force;
}
}
velx[index]+=accelerationx;
vely[index]+=accelerationy;
velz[index]+=accelerationz;
accx[index]=accelerationx;
accy[index]=accelerationy;
accz[index]=accelerationz;
}

The kernel code loops through every object and calculates the forces against every other object. This is the naive unoptimized O(n2) version of the algorithm. Once all the loops are finished the new object velocity and acceleration values are read back from the GPU memory into local memory and then the CPU can access the results. All of the object positions are then updated using the new velocity and acceleration values and then displayed. For displaying the objects I am using the old software only OpenGL billboard quads. A billboard quad is a texture on a quad (rectangle) that always faces the “camera” in OpenGL. If you put a nicely shaded and transparent “blob” as the texture it looks like a simple star and blends in with other stars.

Calculation Times

Using the above code allows me to process millions of objects in a reasonable time frame. For these simulations the slowest part is always the display and CPU calculations. The OpenCL is always the fastest part of the simulations.

Here are a few stats for how quick (slow) these simulations are on a PC with an i7 CPU and GTX 1080 GPU rendering 4K resolution frames;

250,000 objects – 856 ms per frame – 31 ms OpenCL time
500,000 objects – 1606 ms per frame – 57 ms OpenCL time
1,000,000 objects – 3124 ms per frame – 91 ms OpenCL time
5,000,000 objects – 15336 ms per frame – 430 ms OpenCL time
10,000,000 objects – 30638 ms per frame – 864 ms OpenCL time

The reason the per frame time is so much longer than the OpenCL is that I am still using software OpenGL to render all the particles. Software OpenGL falls back to the older v1.1 OpenGL DLLs provided by Microsoft in Windows and has no benefits of hardware acceleration. On my to do list is getting the OpenGL code up to date to use hardware acceleration, then the above times should come way down. Rendering the objects as pixels rather than billboard quads is roughly 1/3 the time, but the billboard quads give a much nicer blended display.

Results

Here is a new sample 4K resolution 3D Gravity movie.

After the last movie I went back and improved the color shading code and added the option for a “black hole”, which in this case is only a single object with a larger mass than the others. The black hole has a mass of around 100 to 500 times the other stars. Any higher and all the stars are flung out of the simulation area too quickly. Here are are some of the latest results.

The spiral galaxy like results are mostly a fluke. I started the simulation with a disk or oblate spheroid (squished sphere) of particles rotating around the origin (Y axis) with a central black hole and let it run.

After a bit more tweaking of the code and coloring algorithms the next movie was ready.

Try It Yourself

The latest 3D Gravity code is now updated and included in Visions of Chaos. For now I am finally happy that at least the basic Newtonian gravity is working and using all objects for the calculations. Future changes will be converting the “old” OpenGL v1.1 display into “modern” OpenGL v3.3+ so the display time can be reduced and get closer to realtime.

A new variety of cellular automata from Jeff Heaton. His original paper describing MergeLife is here, but he also made the following video that clearly explains how the rules work.

Jeff’s MergLife page also has more info and examples you can run online.

MergeLife is now included with Visions of Chaos. I haven’t added the genetic mutations yet, but you can repeatedly click the random and mutate buttons and see what new patterns emerge.