.
First off, congrats to all the WoA participants. I didn't have the time to do it this year, but I'll try again for it the next go around ;)
.
So, in my last little let-my-figure-everything-out-by-writing-stuff-down, I got some stuff wrong, and I got some stuff misleading. Minor things, mostly:
- The last little bit I said about density. Scratch all of that, I'm just an idiot.
- Some of the stuff I said about the actual equation was misleading. No big deal items were incorrect though.
So fluid simulation is a large topic, even though it's governed by only a handful of really intuitive equations. So I'm going to break everything up into pieces. This time, I'm going to be talking about the simplest piece of the equation: advection.
.
Two ways to represent a fluid
.
There are many ways to represent a fluid, but they all fall into one of two categories, or are a hybrid of the two. I discussed them last time, so here they are:
- Lagrangian: fluids are made up of a ton of particles/molecules that interact with one another. This is probably the most intuitive way to visualize a fluid
- Eulerian: fluids are partitioned up into a bunch of little cells that have properties like temperature, velocity, etc. This is probably the most GPU efficient way to visualize a fluid.
There's a really nice connection between the two as well:
.
.
The term on the left side is called the Lagrangian derivative (along with a bunch of other things), and the right side looks very much like part of the original Navier Stokes equation that was reviewed last entry.
.
So what does it mean? Well, imagine a 1D fluid, with a temperature T(x) at each particle, with a constant speed at each particle, and a thermometer placed in the middle at some point:
.
.
The big thing to notice here is that, even IF the individual particles don't change their own temperatures, the thermometer will still feel a change in temperature of the fluid over time. This is because, as different parts of the fluid move, they carry the quantities such as temperature, velocity, etc with them. You can write the situation that was just described as such:
.
.
The first term is 0 because the individual particles don't change temperature (hence, why it's called the Lagrangian derivative). Rearranging terms gets you:
.
.
, which tells you precisely how the thermometer is being effected by the flow of a fluid.
.
We'll take this idea and apply it to the Navier Stokes equation, and drop the viscosity term while we're at it (it's rather annoying):
.
.
This is a much cleaner and simpler form, and it also encompasses 2 different ways to view a fluid: as a group of particles, and as a bunch of thermometers (or really any tool for taking data, like speedometers or barometers) monitoring the current state of the fluid.
.
Advection
.
The problem I'm going to talk about here is an even simpler form of the above equation:
.
.
This is a pure advection equation. It says that some quantity in the fluid is carried along by the fluid itself, and isn't affected by anything else (like pressure or gravity or stuff). We'll look at some of the simplest ways to solve this relationship.
.
Do demonstrate the different techniques, I'll be squirting ink into the velocity field and letting the fluid advect that around, because that's the simplest way to show the differences between the different methods.
.
[a comparison of different techniques. My bet is that Eulerian is going places]
.
It should be abundantly clear that Eulerian is extremely underperformant compared to the semi-Lagrangian versions. This is due to the nature of the problem itself. The Eulerian viewpoint talks about a gradient operator, which is problematic, as we will see. The Semi-Lagrangian techniques are based on interpreting the field q as a bunch of particles that we can then advect simply by moving them.
.
Eulerian Advection
.
Eulerian advection involves directly solving the equation:
.
.
Here's the code snippet for naively solving this thing:
void advectEulerian(T q[cellW][cellH], T q_next[cellW][cellH]){ for (int y = 0; y < cellH; ++y) { for (int x = 0; x < cellW; ++x) { if (x - 1 < 0 || x + 1 > cellW - 1 || y < 0 || y + 1 > cellH - 1) continue; q_next[x][y] = q[x][y]; T dq_dx = (q[x + 1][y] - q[x - 1][y]) / 2.f; T dq_dy = (q[x][y + 1] - q[x][y - 1]) / 2.f; T u_dot_del_q = vel[x][y].x * dq_dx + vel[x][y].y * dq_dy; T dq_dt = -u_dot_del_q; q_next[x][y] += dq_dt * dt; } }}There's nothing fancy here. It's literally evaluating the gradient and then storing the result in a new buffer. There ARE better ways to do this (or so I've been told) by, instead of aligning the velocity components to grid centers, to instead align them to the edges of the cells. This avoids the issue that you see where everything looks all ripply because there's no more null space.
.
Here's an animation for q being ink squirted into the liquid. I don't have a video for q being the velocity buffer though, because it blows up in about 3 seconds
[media]
[/media]
.
Forward Semi-Lagrangian
.
This technique and the next one are, again, exactly what it sounds like. You treat the cell like a particle in the Lagrangian sense, and then advect it forward by the current velocity, not caring at all about how everything is changing around it. Here's the code snippet for that:
template void advectForwardLangrangian(T q[cellW][cellH], T q_next[cellW][cellH]){ for (int y = 0; y < cellH; ++y) { for (int x = 0; x < cellW; ++x) { vec2 advectTarget = vec2(x, y) + vel[x][y] / 60.f; float tu = advectTarget.x - floor(advectTarget.x); float tv = advectTarget.y - floor(advectTarget.y); int xt = floor(advectTarget.x); int yt = floor(advectTarget.y); float ll = (1 - tu) * (1 - tv); float lh = (1 - tu) * tv; float hl = tu * (1 - tv); float hh = tu * tv; if (xt < 0 || xt + 1 > cellW - 1 || yt < 0 || yt + 1 > cellH - 1) { continue; } q_next[xt][yt] += ll * q[x][y]; q_next[xt + 1][yt] += hl * q[x][y]; q_next[xt][yt + 1] += lh * q[x][y]; [xt + 1][yt + 1] += hh * q[x][y]; } }}Again, really simple. The contents of each q in the cell are pushed forward by the cell's velocity and interpolates the difference. As you can see, this looks a LOT smoother. Unfortunately, this method suffers from something called numerical diffusion, which I might talk about at a later date. Here is an animation of that forward semi-Lagrangian action:
.
[media]
[/media]
.
Backwards Semi-Lagrangian
.
So it turns out that, in theory, forward semi-Lagrangian is unconditionally unstable if ?t > ?x / max(vel). I personally haven't had any issues, but whatever, I'm not going to start doubting math now. Fortunately there's a simple solution: just find what particles WILL BE at each grid position after interpolation, and advect backwards that way. This also has the added bonus of being ridiculously simple to do on the GPU, so that's a plus as well. Here's a code sample:
template void advectBackwardsLangrangian(T q[cellW][cellH], T q_next[cellW][cellH]){ for (int y = 0; y < cellH; ++y) { for (int x = 0; x < cellW; ++x) { vec2 advectTarget = vec2(x, y) - vel[x][y] / 60.f; float tu = advectTarget.x - floor(advectTarget.x); float tv = advectTarget.y - floor(advectTarget.y); int xt = floor(advectTarget.x); int yt = floor(advectTarget.y); float ll = (1 - tu) * (1 - tv); float lh = (1 - tu) * tv; float hl = tu * (1 - tv); float hh = tu * tv; if (xt < 0 || xt + 1 > cellW - 1 || yt < 0 || yt + 1 > cellH - 1) { q_next[x][y] = T(); continue; } q_next[x][y] += ll * q[xt][yt]; q_next[x][y] += hl * q[xt + 1][yt]; q_next[x][y] += lh * q[xt][yt + 1]; q_next[x][y] += hh * q[xt + 1][yt + 1]; } }}This looks very much like the forward-Lagrangian. The only difference is that there's a minus where a plus used to be and that you're now writing to the current cell instead of 4 different cells. Upon looking at the animations, there seems to be just about no difference, in this case, to forward semi-Lagrangian. The only strange things that happen are near the edges where backward semi-Lagrangian predicts incorrectly. Here's an animation:
.
[media]
[/media]
.
