# Rebuilding ESP32-fluid-simulation: an outline of the sim task (Part 3)

Okay, I wondered if I should have led this three-part series with the physics, but I think saving it for last was the right call. As I was writing about the FreeRTOS tasks involved and their communication and the touch and render tasks specifically, I started to think about how I could write about this with the detail and approachability it deserves. In fact, this final part will have to be two parts.

To start, I’ll be honest: I’m not presenting anything groundbreaking here. In 1999, Jos Stam introduced a simple and fast form of fluid simulation in his conference paper called “Stable Fluids”, and in 2003, he published a straightforward version of it in “Realtime Fluid Simulation for Games”. Many people have written guides to “fluid simulation” that have been based on these two papers since. Two key examples to me: a chapter of NVIDIA’s *GPU Gems* and a blog post by Jamie Wong (that was, in turn, based on that very chapter!). These were the guides I followed when I first wrote ESP32-fluid-simulation. In both was Stam’s technique, and between everything I just linked to, you could probably write your own implementation of it eventually.

That said, between then and when I rewrote it recently, I picked up some background knowledge that proved incredibly useful. I’m not saying here that I became an expert on fluid sims—I can’t advise you on designing a new technique from scratch. Rather, if I has known it back then, I wouldn’t have made nearly as many wrong turns. It turns out that implementing Stam’s technique gets easier when you understand the *whats* of the operations he would have you write, if not the whys.

#### Review: Vector fields and scalar fields

If you recall any vector calculus, then “vector fields” and “scalar fields” may be an obvious concept to you already, but if not, we can start with the fact that they’re a part of the foundation of fluid dynamics. For now, I’ll review what they are. However, I highly recommend picking up a total understanding of vector calculus somewhere else before looking at any fluid sim techniques besides Stam’s. In fact, perhaps fluid simulations make just the right concrete example to keep in mind while learning!

Anyway, let’s sketch out what vector fields and scalar fields are, and hopefully, the picture is filled in as you keep reading this article. The ordinary idea of a mathematical function is a thing that outputs a number when given a number input. Vector fields and scalar fields are functions—though of different kinds.

Consider a flat, two-dimensional space, and then consider a function that outputs a number when given a *location in this space* as the input. Furthermore, this location can be expressed as a pair of numbers if we used a coordinate system (three if we worked in three dimensions, and we could, but we won’t here). A concrete example of this would be a function of a location that gives the temperature there—the location being written as the latitude and longitude on the map. It’s 48 degrees Fahrenheit in Arkhangelsk and 84 degrees in Singapore. Considering that Arkangelsk can be found at 64.5°N, 40.5°E and Singapore at 1.2°N, 103.8°E, we can define a temperature function that gives $T(64.5, 40.5) = 48$ and $T(1.2, 103.8) = 89$. We can call it a temperature field, but more generally, it’s a “scalar field”. It’s a scalar-valued function of the location, possibly written as $f(x, y)$.

Now on the other hand, a “vector-valued function” is any function that outputs a vector, and a vector-valued function of a location is a “vector field”! A concrete example of this would be a wind velocity field. For any location, such a field would give how fast the wind there blows and the direction in which it goes, and it would be given as the magnitude and direction of a single vector. Like for scalar fields, we could possibly write them as $\bold{f}(x, y)$, the boldface font meaning that we have a vector output.

That said, though functions of location they are, written like one they are really not. Rather, the dependence on location is assumed, and then $f(x, y)$ and $\bold{f}(x, y)$ are just written as $f$ and $\bold{f}$ instead. Another thing to keep in mind: coordinates are just a pair of numbers, but we can also think of them as a single coordinate *vector*. Though we may never actually draw that arrow, the interchangeability is relevant. For example, I briefly talked in the previous post about the similarity between a velocity vector and a *change* in the coordinate vector over a finite period of time.

First, it would be helpful to picture what we want to simulate. The input and output are the *touch* and *screen* of a touchscreen, and the user dragging around the stylus on it should stir around the fluid on display. The physical scenario this should match is if someone stuck their arm into a bed of dyed water and then stirred it around. In such a scenario, the color would be determined by the concentration of the dye, but the dye itself moves! To capture this physical behavior with a computer simulation, we can start by describing it with a mathematical model.

In Stam’s “Real-time Fluid Dynamics for Games”, he wanted to capture smoke rising from a cigarette and being blown around by air currents. To do so, he ascribed a velocity field (a vector field) and a smoke density field (a scalar field) to the air. But that was it for his model: everything else about it he threw out. In the same way, we can reduce the bed of water to just a velocity field and a dye concentration field.

Now, what was the relationship between these two fields? Stam wrote that the density field undergoes “advection” by the velocity field. That’s the process of fluid carrying around (smoke particles, dye, or anything in general), and this happens everywhere. He also wrote that it undergoes “diffusion”, which is the spontaneous spreading of a thing in a fluid from areas of higher density *without* being carried by the velocity. He provided an “advection-diffusion” equation that captures both, and it’s a “partial differential equation”.

#### Review: Partial derivatives and the differential operators

Just like how we can take the derivative of your ordinary function, we can take a differential operator of a field. However, these differential operators don’t just mean the slope of a tangent line, but rather they each represent a different way the field changes over a change in location. The critical ones to understand here are the “divergence” and the “gradient”, but the “Laplacian” is also worth touching on. (A formal vector calculus course would also cover the “curl”, the identities, and the associated theorems.)

First of all, differential operators are constructed from the “partial derivatives”. These are the derivatives you already know, but we strictly take them with respect to *one* of the components while holding the others constant. The reason? Formally, your ordinary derivative is the limit of the change in your ordinary function $f(x)$ over the change in the input $x$ as that change in the input approaches zero.

However, in the case of fields, by doing this to only *one* of the components of the location coordinate, the partial derivative just formally means the change in the field $f(x, y)$ over the change in *that component*. Keeping the other components constant is naturally a part of measuring this change. In two dimensions, fields can have a partial derivative with respect to $x$ or one with respect to $y$. Then, $y$ or $x$ respectively is held constant.

A good example would actually be to perform a derivation. Given the function $f(x, y) = x^2 + 2xy + y^2$ as a field, let’s find the partial derivative with respect to $x$.

\[\begin{align*} \frac{\partial}{\partial x}(x^2 + 2xy + y^2) & = 2x + 2y + 0 \\ & \boxed{ = 2x + 2y } \end{align*}\]Notice that—because $y$ is taken as a constant—$y^2$ drops out and $2xy$ is treated as an $x$-term with a coefficient of $2y$. And finally, to expand on this a bit with a geometric picture, we know that the derivative is the slope of the tangent line, but to be exact, it’s the line tangent to the curve of $f(x)$ at the point $(x, f(x))$. The partial derivative is still the slope of *a* line that *is tangent* to the surface of the field at the point $(x, y, f(x, y))$, but it is also strictly running in the $x$-direction for $\partial/\partial x$ or in the $y$-direction for $\partial/\partial y$. Technically, infinitely many lines satisfy the conditions of being tangent to the surface at that point, and these lines form a tangent plane, but we only concern ourselves with the two.

That aside, taking a partial derivative with respect to some single component is not as useful as taking *every* partial derivative with respect to *each* component. This set is written like a vector of sorts (though a vector it is not) called the “del operator”. For two dimensions, that is

The constructions out of this set that we call the differential operators can absolutely be written without using the del operator, but you’d usually see that they are.

The “gradient” is the simplest construction: line up each and every partial derivative of a scalar field into a vector. Keeping in mind here that the partial derivative of a field (like $x^2+2xy+y^2$) is actually yet another function of the location (like $x^2 + 2y$), a vector composed of these will itself vary by the location. The gradient of a scalar field is a vector field! We can get to exactly how the gradient gets applied to our fluid sim later, but one useful fact to picture here is that it can be shown that the gradient always happens to point in the direction of steepest ascent in the scalar field. Walking in the direction of the gradient of the temperature field, for example, would warm you up the fastest!

\[\nabla f = \begin{bmatrix} \displaystyle \frac{\partial f}{\partial x} \\[1em] \displaystyle \frac{\partial f}{\partial y} \end{bmatrix}\]Using the del operator, it looks kind of like scalar multiplication from the right.

And remember, the gradient is just one shockingly meaningful operator that we can construct from the partial derivatives, which were just slopes of tangent lines! The “divergence” is a slightly more complicated construction: if we write out a vector field using its components

\[\bold{f}(x, y) = \begin{bmatrix} f_x(x, y) \\ f_y(x, y) \end{bmatrix}\]then we can take the partial derivative of each component with respect to its *associated component* of the coordinates (that’s $f_x$ to $\partial/\partial x$ and $f_y$ to $\partial/\partial y$) and then add them up. We should be able to recognize here that the divergence of a vector field is a scalar field. And what is the meaning of this scalar field? For now, it can be imagined as the degree to which the vectors surrounding an input location are pointing away from it, though Gauss’s theorem expresses this more formally (a bit out-of-scope for now).

Using the del operator, it looks kind of like a dot product.

Finally, the “Laplacian” is actually the divergence of the gradient of a scalar field, and this ultimately means that it’s also a scalar field! It is also the sum of the second-order partial derivatives (besides the mixed ones, but that’s totally out-of-scope).

\[\nabla^2 f \equiv \nabla \cdot (\nabla f) = \frac{\partial^2 f}{\partial x^2} + \frac{\partial^2 f}{\partial y^2}\]Using the del operator, some take the liberty of expressing this composition as a single $\nabla^2$ operator.

There is also the extension of the Laplacian onto vector fields, but it really is just the Laplacian on each component.

\[\nabla^2 \bold{f} = \begin{bmatrix} \nabla^2 f_x \\ \nabla^2 f_y \end{bmatrix}\]The gradient, divergence, and Laplacian are all the differential operators that are relevant here, and hopefully these will become more concrete to you as we use them from here on to describe Stam’s fluid sim technique. However, I’d again recommend formally learning vector calculus if you’d like to look at other techniques.

A “partial differential equation” is kind of like a system of linear equations in this context. Here, they still relate known and unknown variables, and they still have a solution which is the value of the unknowns. However, these “values” are entire fields, not just numbers! Given this, partial differential equations also involve the differential operators of these field-valued variables.

The advection-diffusion equation that Stam provides is a simple example of one: advection and diffusion are *independent terms*, and their sum is exactly how the density field evolves over time. It is

where $\rho$ is the density field and $\bold{v}$ is the velocity field. $-(\bold{v} \cdot \nabla) \rho$ is the advection term, and $\kappa \nabla^2 \rho$ is the diffusion term—$\kappa$ being a constant for us to control the strength of the diffusion. $S$ is just a term that lets us add density (of smoke, or concentration of dye in our case) to the scene. Notice how this equation is a definition of $\partial \rho / \partial t$. It’s the partial derivative of the density field with respect to time, and it means that $\rho$ is a variable whose value is a function of location and time. However, it’s more useful for us to think of it equivalently as a field that evolves over time. An evolving density field is exactly what we want to show on the screen!

You may also notice that $(\bold{v} \cdot \nabla)$ is clearly some kind of construction from the partial derivatives because it uses the del operator $\nabla$. That is the “advection” operator. I’ve only seen it in fluid dynamics papers and yet still don’t totally understand it. Still, we’ll see how Stam treats it, but that’ll have to be in the next post.

That aside, I’m going to adjust this equation while we’re here. Because the overall project was about simulating dye in water on an ESP32 and not smoke in air on a GPU, I didn’t use the entire equation. Anyway, this can be thought of as an exercise in finding what part of the physics can be ignored while still looking sorta-physical, I suppose. I really have to wash my hands of any assertions I’m making at this moment, for I am no expert in this field. With that said, I can confirm that deleting the diffusion term by letting $\kappa = 0$ doesn’t look so egregious. We also don’t have to add more dye, so we can delete $S$ while we’re at it. That actually just leaves the advection alone.

\[\frac{\partial \rho}{\partial t} = -(\bold{v} \cdot \nabla) \rho\]All said though, where is the room in this model for the user’s input? Is the velocity field just a thing we get to set? (Right now, we have two unknowns, $\rho$ and $\bold{v}$, but one equation!) No, it’s more complicated than that: the way water and air move continues to change even after we stop stirring it. That leads to the missing piece to stirring digital water: we need a physical way to define $\partial \bold{v} / \partial t$ (a.k.a. the acceleration!) just like how $\partial \rho / \partial t$ was defined. That missing piece is the famous “Navier-Stokes equations”.

The “Navier-Stokes equations” are also partial differential equations. A definition of Navier-Stokes can be found in any fluid dynamics article, but the one Stam provided in “Stable Fluids” is the most directly relevant one.

\[\frac{\partial \bold{v}}{\partial t} = -(\bold{v} \cdot \nabla) \bold{v} - \frac{1}{\rho} \nabla p + \nu \nabla^2 \bold{v} + \bold{f}\] \[\nabla \cdot \bold{v} = 0\]The first one is a definition of $\partial \bold{v} / \partial t$. Here, $-(\bold{v} \cdot \nabla) \bold{v}$ and $\nu \nabla^2 \bold{v}$ represent advection and diffusion again, though these are also known as “convection” and acceleration due to “viscosity”, respectively. That is to say, the velocity is carried around and diffused just like how the smoke density was. The only difference is that the constant $\nu$ here is the “kinematic viscosity”, and it’s higher for fluids like honey and lower for fluids like water. A quick note here: I can also confirm that letting $\nu = 0$ has no severe consequences. That aside, $\bold{f}$ represents the acceleration due to external forces, and there is the place in our mathematical model where the user input would go!

On the other hand, $-\frac{1}{\rho} \nabla p$ is an interesting term—it’s *not* independent here. Let me try to explain. Nominally, it’s the acceleration due to a difference in pressure, and the negative of the gradient represents the tendency for fluids to flow from regions of high pressure to regions of low pressure. Since the gradient points in the direction of steepest ascent, then going in the opposite direction gives the steepest *descent*. But this is only in name! As Stam had put it, “[t]he pressure and the velocity fields which appear in the Navier-Stokes equations are in fact related”. Ultimately, he used it as a correction term in order to guarantee that the second equation holds. It’s a constraint on the velocity field that is critical to looking like water, called “incompressibility”.

Simply, it states the divergence is zero. Unfortunately, knowing *what* the incompressibility constraint is turns out to not be the same as knowing *why* it is nor why incompressibility matters. Those are things beyond what I can comfortably explain in the first place (which I’d think is worse than it being out-of-scope). On the other hand, I definitely have my understanding of how the so-called pressure term is used to ensure zero divergence! However, that’ll again have to be another matter to cover in the next post.

There is one final thing I want to end with. What we’ve covered so far is the mathematical model that Stam used to capture how fluids move (air in his case, but water being captured equivalently). It’s partial differential equations, and the solution is a density field that evolves over time. But computers can’t think in terms of fields. Instead, continuous space can be approximated by a mesh of discrete points, each point taking on the *value* of the field there. Remember that fields are functions of the location, and so the value of a field at a point is a single scalar or vector. This step is generally called “discretization”, and its approximate solutions are generally called “numerical solutions”. In our case, Stam used a grid of points, and the incredibly convenient thing is that fields discretized with a grid become an *array of values*. Each point on the grid $(x_j, y_i)$ corresponds with a value in that array `f[i, j]`

. (If you’re wondering why I write $(x_j, y_i)$ and not $(x_i, y_j)$, see my review of the AdafruitGFX coordinate system in my last post. I’m just keeping it consistent, though the math still holds if you work in ordinary Cartesian coordinates.)

That is not to say that this array is a *matrix*. The array is only two-dimensional because the space is two-dimensional. If the space was three-dimensional, then so would be the array. And you can forget about arrays if the mesh wasn’t a grid! Most matrix operations wouldn’t mean anything either. It’d be more correct to think of it as a very long vector, but that is another matter for another day.

Anyway, computers can handle these arrays. The other thing discretization does is turn the differential operators into discrete sums, another thing computers can handle! In our case, the grid does another wonderful thing: it makes these discrete sums incredibly simple. To begin with, in a grid, a point has four neighbors. From there $\partial / \partial x$ can be approximated by subtracting the left value from the right value, and the same can be done for $\partial / \partial y$ by subtracting the top value from the bottom value (again with AdafruitGFX: “next row *down*” and not “next row *up*”, and positive $y$-velocity means downward movement). The morphing of the partial derivative into just a difference is why discretizations using a grid are called “finite difference” methods.

From there, we can repeat the construction of the divergence and the gradient, just as if we were actually using partial derivatives.

\[\nabla f \approx \begin{bmatrix} \displaystyle \frac{f[i, j+1]-f[i, j-1]}{2 \Delta x} \\[1em] \displaystyle \frac{f[i+1, j]-f[i-1, j]}{2 \Delta y} \end{bmatrix}\] \[\nabla \cdot \bold{f} \approx \frac{f_x[i, j+1]-f_x[i, j-1]}{2 \Delta x} + \frac{f_y[i+1, j]-f_y[i-1, j]}{2 \Delta y}\]That said, though the Laplacian is just the divergence of the gradient, constructions of the discrete Laplacian are *not* done this way. Instead, replacing $\partial^2 / \partial x^2$ and $\partial^2 / \partial y^2$ with a different finite difference is the starting point.

And remember, the Laplacian is defined for vector fields *and* scalar fields, where it’s vector definition is just the scalar Laplacian on each component. So, we can define the discrete vector Laplacian as the discrete scalar Laplacian on each component. But incidentally, we can also just write $\bold{f}$ for $f$ here, and that definition would be exactly the same!

With one exception (another next post matter), Stam used finite differences to solve the Navier-Stokes equations and also the advection-diffusion equation from there. That said, there are many, many possible discretizations of the differential operators, each with its own advantages and disadvantages. (I’ve written a Smoothed Particle Hydrodynamics (SPH) solver before, which used discretization with *moving* points!) Though I can’t rightly compare one discretization against another on my own authority, I can certainly comment that they don’t get much simpler than finite differences. In fact, I find the advection-diffusion equation and the Navier-Stokes equations easier to think about if I look at them in terms of finite differences, instead of just looking at them in terms of the pure differential operators.

So ends this post. With the governing equations (advection-diffusion and Navier-Stokes) and their discretization, we’ve laid out the fundamental outline of Stam’s technique. We’ve also reviewed the relevant vector calculus, though no more than that. Though I didn’t have all the authority I needed to get the whys, we should be equipped to understand the whats. In the second part, we’ll fill in the outline to get an end-to-end fluid simulation. In particular, that “exception” will become the star of the show. If you’re still here before that post comes out though, there’s always the ESP32-fluid-simulation source code on GitHub.