Viscous Diffusion (and Poisson Equations)

In this chapter:

We now move to the second step of calculating our intermediate velocity field, simulating viscous diffusion. Viscous diffusion is often just called “diffusion”, the reason we say it is viscous diffusion is because it is primarily governed by the viscosity of the fluid. Here’s a review from the start of this chapter to jog your memory:

Diffusion is the movement of some quantity in the fluid from a region of higher concentration to a region of lower concentration. Diffusion is primarily controlled by viscosity.

Like advection (and all of our future steps) diffusion is also carried out as a slab operation, and picks up after our advection calculation is done.

Deriving the Diffuse Update Equation

At this point I’m sure you’re probably tired of me saying this , but as always, to figure out what we need to do we start from our Navier-Stokes equations:

∂u⃗∂t=−(u⃗⋅∇)u⃗−1ρ∇p+ν∇2u⃗+F⃗\frac{\partial \vec{u}}{\partial t} = -(\vec{u}\cdot\nabla)\vec{u} - \frac{1}{\rho}\nabla p + \nu \nabla^2 \vec{u} + \vec{F}

In this step we are only concerned about applying viscous diffusion, so ignore the other terms and let’s just look at how the diffuse term affects the change in velocity:

∂u⃗∂t=−(u⃗⋅∇)u⃗−1ρ∇p+ν∇2u⃗+F⃗\frac{\partial \vec{u}}{\partial t} = \xcancel{\textcolor{red}{-(\vec{u}\cdot\nabla)\vec{u} - \frac{1}{\rho}\nabla p}} + \nu \nabla^2 \vec{u} + \xcancel{\textcolor{red}{\vec{F}}} ∂u⃗∂t=ν∇2u⃗\frac{\partial \vec{u}}{\partial t} = \nu \nabla^2 \vec{u}

Recall that $\partial \vec{u}/\partial t$ is telling us about the “rate of change of velocity” or “acceleration” of the velocity vector field. Of course, our goal during this step is to change the velocity, so this is exactly what we’re looking for.

The only issue is that this equation is in acceleration terms, and in order to update our velocity, we want to know what the new velocity is after a known time-step. The key to this of course, is time.

What we’re going to do is approximate this solution using an explicit Euler step (which sounds complicated but is actually really simple).

From a Car to a Fluid

We create an analogy of our fluid particle as a car, moving only in one dimension for simplicity:

A car is driving on the road. Suppose that at the current moment, it’s velocity is $u$, and the driver has his foot on the gas pedal creating a constant acceleration of $a$. We want to know how much the velocity changed after $\Delta t$ seconds.

Since we know that the definition of acceleration is “the rate of change of velocity with respect to time”, or in other words, “how much the velocity changed after a certain period of time has passed”, then we can infer that the change in velocity is simply the acceleration multiplied by the time passed:

Δu=a⋅Δt(1)\tag{1} \Delta u = a \cdot \Delta t

If we know how much the velocity changed, and we know the old velocity, then the new velocity is:

unew=uold+Δuu_{\text{new}} = u_{\text{old}} + \Delta u

Equation 1 already tells us what $\Delta u$ is, so we can substitute it to find:

unew=uold+(a⋅Δt)u_{\text{new}} = u_{\text{old}} + (a \cdot \Delta t)

Our last step is to figure out what $a$ is. Just a little while ago we said:

Recall that $\partial \vec{u}/\partial t$ is telling us about the “rate of change of velocity” or “acceleration” of the velocity vector field.

And that’s it really. Extend the equation to a vector equation and substitute $a$ as our viscous term $\nu \nabla^2 \vec{u}$ from the Navier-Stokes equations:

u⃗new=u⃗old+(ν⋅(∇2u⃗))⋅Δt\vec{u}_{\text{new}} = \vec{u}_{\text{old}} + (\nu \cdot (\nabla^2 \vec{u})) \cdot \Delta t

We can write this a little more neatly as:

u⃗new=u⃗old+νΔt⋅(∇2u⃗)      (Diffuse Update Equation)\boxed{\vec{u}_{\text{new}} = \vec{u}_{\text{old}} + \nu \Delta t \cdot (\nabla^2 \vec{u})} \;\;\; \text{(Diffuse Update Equation)}

Congratulations, you have just learned how to numerically approximate a differential equation using Euler’s method.

Diffusion Method 1: Direct Euler Step

What? Did you really think that was the end of the math?

It is possible to apply the explicit Euler step we just derived in our shader. It requires computing the Laplacian of the velocity field, which we can do using finite difference methods that we talked about here.

From that page, the finite difference approximation for $\nabla^2 \vec{u}$ is

∇2u⃗=[ui+1,j−2ui,j+ui−1,j(δx)2+ui,j+1−2ui,j+ui,j−1(δy)2vi+1,j−2vi,j+vi−1,j(δx)2+vi,j+1−2vi,j+vi,j−1(δy)2]\nabla^2 \vec{u} = \LARGE \begin{bmatrix} \frac{u_{i+1,j} - 2u_{i,j} + u_{i-1,j}}{(\delta x)^2} + \frac{u_{i,j+1} - 2u_{i,j} + u_{i,j-1}}{(\delta y)^2} \\[20pt] \frac{v_{i+1,j} - 2v_{i,j} + v_{i-1,j}}{(\delta x)^2} + \frac{v_{i,j+1} - 2v_{i,j} + v_{i,j-1}}{(\delta y)^2} \end{bmatrix}

Since we’re only looking one pixel over to the left or right or up or down, we can set $\delta x = \delta y = 1$. We can rewrite this much simpler as:

∇2u⃗=[ui+1,j+ui−1,j+ui,j+1+ui,j−1−4ui,jvi+1,j+vi−1,j+vi,j+1+vi,j−1−4vi,j]\nabla^2 \vec{u} = \Large \begin{bmatrix} u_{i+1,j} + u_{i-1,j} + u_{i,j+1} + u_{i,j-1} - 4u_{i,j} \\ v_{i+1,j} + v_{i-1,j} + v_{i,j+1} + v_{i,j-1} - 4v_{i,j} \end{bmatrix}

Which can be further reduced to vector form as:

∇2u⃗=u⃗i+1,j+u⃗i−1,j+u⃗i,j+1+u⃗i,j−1−4u⃗i,j\nabla^2 \vec{u} = \vec{u}_{i+1,j} + \vec{u}_{i-1,j} + \vec{u}_{i,j+1} + \vec{u}_{i,j-1} - 4\vec{u}_{i,j}

For our own understanding, we can rewrite the subscripts using their relative position to the current velocity vector:

∇2u⃗=u⃗right+u⃗left+u⃗up+u⃗down−4u⃗\nabla^2 \vec{u} = \vec{u}_{\text{right}} + \vec{u}_{\text{left}} + \vec{u}_{\text{up}} + \vec{u}_{\text{down}} - 4\vec{u}

In code we can write this as:

vec2 c     = Tv(p).xy;
vec2 up    = Tv(p + vec2(0,1)).xy;
vec2 right = Tv(p + vec2(1,0)).xy;
vec2 down  = Tv(p - vec2(0,1)).xy;
vec2 left  = Tv(p - vec2(1,0)).xy;

// Viscous Diffusion (explicit Euler step using Laplacian calculation)
vec2 laplacian = (right + left + up + down - 4.0 * c);
vec2 new_vel = c + dt * nu * laplacian;

Which utilizes the same Tv() function we defined on the previous page. nu is the value of the viscosity $\nu$ which is user-defined.

Compute Shader Implementation

As a full compute pass, method 1 would look like:

#version 450

layout(set = 0, binding = 0, rg32f) uniform image2D VELOCITY_INPUT;
[...]
layout(push_constant, std430) uniform MouseData {
	vec2 mouse_pos;
};

vec4 Tv(vec2 _texel)
{
	[...] // Copy from before
}

layout(local_size_x = 8, local_size_y = 8, local_size_z = 1) in;
void main()
{
	ivec2 texel = ivec2(gl_GlobalInvocationID.xy);
	vec2 ftexel = vec2(texel); // Convert from int to float
	vec2 p = ftexel;

	vec2 c     = Tv(p).xy;
	vec2 up    = Tv(p + vec2(0,1)).xy;
	vec2 right = Tv(p + vec2(1,0)).xy;
	vec2 down  = Tv(p - vec2(0,1)).xy;
	vec2 left  = Tv(p - vec2(1,0)).xy;

	// Viscous Diffusion (explicit Euler step using Laplacian calculation)
	float nu = parameters.viscosity;
	float dt = parameters.dt;
	vec2 laplacian = (right + left + up + down - 4.0 * c);
	vec2 new_vel = dt * nu * laplacian;

	// Write out data
	vec4 color = vec4(new_vel, 0.0, 0.0);
	imageStore(VELOCITY_OUTPUT, texel, color);
}

Keep in mind that after this pass is run you still need to copy the data from VELOCITY_OUTPUT to VELOCITY_OUTPUT in order to complete the full slab operation like we discussed on the last page.

Issues

Applying the viscous diffusion using the direct Euler step we derived has one big issue. Technically, the fluid particles can move through multiple cells with different velocities during a single time-step, but we only calculate the Laplacian from the nearby pixels, which means different cells can apply different parts of diffuse resistance (Fernando, 2004).

The result is that this method is not stable at large time-steps or large viscosities and can numerically explode. However, it does have one HUGE benefit over the next method we are going to discuss, and that is that it only requires a single pass to calculate (which makes it extremely fast).

Diffusion Method 2: Poisson Solver

The diffusion update equation we derived earlier uses Euler’s method, specifically, this is an explicit solution derived using Forward-Euler approximation. In order to mitigate the issues with method 1 described above, we make use of an implicit method described by Stam (1999).

Stam (1999) proposes the following implicit equation:

(I−ν⋅Δt⋅∇2)u⃗new=u⃗old(\textbf{I} - \nu \cdot \Delta t \cdot \nabla^2) \vec{u}_{\text{new}} = \vec{u}_{\text{old}}

where $\textbf{I}$ is the identity matrix.

From Explicit to Implicit

In my research, I wasn’t able to find anyone who directly explained how you get to this implicit equation so I am going to show you here.

First, we need to understand the difference between an explicit method and implicit method and why the implicit method mitigates the issues explained earlier.

The car analogy I used earlier is just a nicer physical interpretation for how to come up with a Forward-Euler approximation, but at it’s core, what Euler’s method really does is figure out a tangent line or gradient to the function, and determines the next solution step by nudging the current solution a small amount along the tangent line / gradient. The “small amount” specifically is the time-step.

We already see this in our equation from earlier:

unew=uold+(a⋅Δt)u_{\text{new}} = u_{\text{old}} + (a \cdot \Delta t)

Which in simple terms says “take $u_{\text{old}}$ and nudge it by the gradient $a$ over a small amount $\Delta t$ and you will find an approximation for $u_{\text{new}}$”.

The difference between Forward-Euler and Backward-Euler is simply where you evaluate the gradient (Zeltkevic, 1998).

If you evaluate the gradient at the current time-step, it is an explicit method, if you evaluate the gradient at the future time-step, that’s an implicit method.

Forward-Euler method evaluates the gradient at the current time-step, making it an explicit method. Backward-Euler method evaluates the gradient at the future time-step making it an implicit method.

Typically, explicit methods are only conditionally stable. That means that they are only stable if certain conditions are met. For the Forward-Euler method this is “the existence of a critical time step size beyond which numerical instabilities manifest” (Zeltkevic, 1998).

Implicit methods can combat this but have their own issue: implicit methods are usually a lot more expensive to calculate.

For the Forward-Euler method derived earlier, the exact gradient calculation is done as follows:

u⃗new−u⃗oldΔt=u⃗old\frac{\vec{u}_{\text{new}} - \vec{u}_{\text{old}}}{\Delta t} = \vec{u}_{\text{old}}

The left hand side of this equation is just the straight-line gradient calculation, and the right hand side is the known gradient (a.k.a the viscous acceleration term). Notice how the viscous acceleration term contains $\vec{u}_{\text{old}}$, which means that we are evaluating the gradient at the current time-step.

So to convert this into an implicit method we simply change the gradient evaluation to $\vec{u}_{\text{new}}$ giving us:

u⃗new−u⃗oldΔt=u⃗new(2)\tag{2} \frac{\vec{u}_{\text{new}} - \vec{u}_{\text{old}}}{\Delta t} = \vec{u}_{\text{new}}

Rearranging equation 2 gives us:

u⃗new−νΔt∇2u⃗new=u⃗old(3)\tag{3} \vec{u}_{\text{new}} - \nu \Delta t \nabla^2 \vec{u}_{\text{new}} = \vec{u}_{\text{old}}

and if you factor out the $\vec{u}_{\text{new}}$ term you get the same equation derived by Stam (1999):

(I−ν⋅Δt⋅∇2)u⃗new=u⃗old(\textbf{I} - \nu \cdot \Delta t \cdot \nabla^2) \vec{u}_{\text{new}} = \vec{u}_{\text{old}}

Although I personally don’t like this representation because it makes explaining the next step more difficult. So for simplicity we’ll just stick to what’s written in equation 3, but now at least you know where Stam’s equation comes from.

Introduction to Poisson Equations

Equation 3 is a type of equation with a special name, it is called a Poisson equation. This is actually a good thing, because Poisson equations pop up in many different fields of science for many different problems. So, a lot of different methods to solve these types of equations already exist.

There is another reason I bring up the topic of Poisson equations. We are going to need to solve one more Poisson equation later in our fluid simulator, so it is a good idea to outline a Poisson equation solver now that we can re-use later.

Poisson equations have the general form:

∇2ϕ=f\nabla^2 \phi = f

where $\phi$ is some unknown function we wish to solve for, and $f$ is a known function.

Immediately, you might notice a problem. Equation 3 that we derived earlier is not in this form at all:

u⃗new−νΔt∇2u⃗new=u⃗old\vec{u}_{\text{new}} - \nu \Delta t \nabla^2 \vec{u}_{\text{new}} = \vec{u}_{\text{old}}

Converting equation 3 to Poisson form is simply a matter of algebraic manipulation. This is why I wanted to avoid doing the extra step to get to Stam’s equation:

∇2u⃗new=u⃗new−u⃗oldνΔt(Poisson Velocity Diffusion Equation)(4)\tag{4} \boxed{\nabla^2 \vec{u}_{\text{new}} = \frac{\vec{u}_{\text{new}} - \vec{u}_{\text{old}}}{\nu\Delta t}} \\[6pt] \text{(Poisson Velocity Diffusion Equation)}

Discretizing Poisson Equations

Set aside equation 4 for now.

We want to write an efficient, GPU-friendly Poisson equation solver that we can use for this equation and the one other Poisson equation we will see in later in our simulator. In order to do this, we have to first generalize these equations.

The most common method to solve these equations on the GPU is by employing a discretization of the Laplacian term and then applying an iterative solver to this discretization. According to Fernando (2004) the Poisson equations can be discretized and rewritten into the form:

xnew=xleftold+xrightold+xupold+xdownold+α⋅bβ(5)\tag{5} x^{\text{new}} = \frac{x^{\text{old}}_{\text{left}} + x^{\text{old}}_{\text{right}} + x^{\text{old}}_{\text{up}} + x^{\text{old}}_{\text{down}} + \alpha \cdot b}{\beta}

where $x$ is whatever quantity you are looking for.

Let’s see how we can reconfigure our Poisson diffusion equation into this format. Recall that method 1 taught us how to calculate Laplacians using finite differences. Let’s just go ahead and write out the finite difference approximation for $\nabla^2 \vec{u}_{\text{new}}$.

Note: I am going to write $\vec{u}_{\text{new}}$ as $\vec{u}^{\text{new}}$ instead so we can use the subscript space to denote which cell we are reading

∇2u⃗new=u⃗rightnew+u⃗leftnew+u⃗upnew+u⃗downnew−4u⃗new(6)\tag{6} \nabla^2 \vec{u}^{\text{new}} = \vec{u}_{\text{right}}^{\text{new}} + \vec{u}_{\text{left}}^{\text{new}} + \vec{u}_{\text{up}}^{\text{new}} + \vec{u}_{\text{down}}^{\text{new}} - 4\vec{u}^{\text{new}}

Plug equation 6 into equation 4:

u⃗rightnew+u⃗leftnew+u⃗upnew+u⃗downnew−4u⃗new=u⃗new−u⃗oldνΔt\vec{u}_{\text{right}}^{\text{new}} + \vec{u}_{\text{left}}^{\text{new}} + \vec{u}_{\text{up}}^{\text{new}} + \vec{u}_{\text{down}}^{\text{new}} - 4\vec{u}^{\text{new}} = \frac{\vec{u}^{\text{new}} - \vec{u}^{\text{old}}}{\nu\Delta t}

Rearrange this and make $\vec{u}^{\text{new}}$ the subject:

u⃗new=u⃗rightnew+u⃗leftnew+u⃗upnew+u⃗downnew+1νΔtu⃗old(1νΔt+4)\vec{u}^{\text{new}} = \frac{\vec{u}_{\text{right}}^{\text{new}} + \vec{u}_{\text{left}}^{\text{new}} + \vec{u}_{\text{up}}^{\text{new}} + \vec{u}_{\text{down}}^{\text{new}} + \frac{1}{\nu \Delta t}\vec{u}^{\text{old}}}{(\frac{1}{\nu \Delta t} + 4)}

This is slightly different from the format we see in equation 5. We cannot really compute the Laplacian of the new field until we have calculated the values of the new field. Really all we can do is approximate the Laplacian of the new field as the Laplacian of the old field. As you will see later, with enough iterations this won’t matter as the values of the old field and new field will become really close. So our equation changes to:

u⃗new=u⃗rightold+u⃗leftold+u⃗upold+u⃗downold+1νΔtu⃗old(4+1νΔt)(Discrete Poisson Diffusion Equation)(7)\tag{7} \boxed{\vec{u}^{\text{new}} = \frac{\vec{u}_{\text{right}}^{\text{old}} + \vec{u}_{\text{left}}^{\text{old}} + \vec{u}_{\text{up}}^{\text{old}} + \vec{u}_{\text{down}}^{\text{old}} + \frac{1}{\nu \Delta t}\vec{u}^{\text{old}}}{(4 + \frac{1}{\nu \Delta t})}} \\[6pt] \text{(Discrete Poisson Diffusion Equation)}

Comparing equation 7 to equation 5 we can determine what all of the values should be set to in our Poisson equation solver that we are going to develop next:

{b=u⃗oldx=u⃗oldα=1νΔtβ=4+α\Large \begin{cases} b = \vec{u}^{\text{old}} \\ x = \vec{u}^{\text{old}} \\ \alpha = \frac{1}{\nu \Delta t} \\ \beta = 4 + \alpha \end{cases}

Parallel Computation of Numerical Solvers for Poisson Equations

Method 1, although only conditionally stable, doesn’t require any solver routines and can be directly applied using the code previously shown. In this section we focus only method 2, which requires the attachment of a numerical solver to compute the solution to the discrete Poisson equation we just derived.

These methods are applicable to any discrete Poisson equation. Here we are applying it for the viscous diffusion step, but later, we can use any of these methods also to solve the pressure Poisson equation we are going to develop later.

While it is possible for you to skip developing a Poisson solver now and just stick to using method 1. When we get to pressure solving to remove the velocity field divergence, we encounter an unavoidable Poisson equation that we need to solve.

So, it is a good idea to understand the Poisson solver types available to you now so you can choose the one best suited to your needs.

Jacobi’s Iterative Method

The main benefit of this method is how easy it is to understand and implement in any GPU program. Other than that, in my opinion, the Jacobi method is terrible in almost all other aspects.

That doesn’t mean that it shouldn’t be used however. Many simple incompressible fluid simulators you will find on Shadertoy or various other parts of the internet employ the Jacobi method to solve their Poisson equations. Sometimes, you might be even forced to use this method because of framework limitations. Just because something isn’t as good as other methods, it doesn’t mean that it is useless.

If you have studied a Gauss-Seidel routine at any point in your life, this method is very similar but essentially just a step down.

Jacobi’s iterative method is done by direct application of equation 5. It is also carried out as a slab operation, that means:

  1. Read the values of each pixel and its 4 neighbors from old texture
  2. Plug these values into equation 5 to calculate new values
  3. Write new values to new texture
  4. Copy data from new texture to old texture
  5. Repeat steps 1-4 many times

Steps 1 to 4 is composes 1 Jacobi iteration. The main issue with Jacobi’s method is how many iterations it requires to converge to a good answer. For our fluid sim it is recommended to run 50 Jacobi iterations! (Fernando, 2004)

In code it looks something like this (Fernando, 2004):

vec4 jacobi(ivec2 texel_coords, float alpha, float beta, sampler2D x, sampler2D b)
{
	vec4 xL = texelFetch(x, texel_coords + ivec2(-1, 0), 0);
	vec4 xR = texelFetch(x, texel_coords + ivec2( 1, 0), 0);
	vec4 xD = texelFetch(x, texel_coords + ivec2( 0,-1), 0);
	vec4 xU = texelFetch(x, texel_coords + ivec2( 0, 1), 0);

	vec4 bC = texelFetch(b, texel_coords, 0);

	return (xL + xR + xU + xD + alpha * bC) / beta;
}

Poisson Filters

The first paper on the application of this topic to solving real-time fluid sims was developed by Rabbani & Khiat (2020), and therefore this is considered a fairly new method of solving Poisson equations on the GPU.

This method actually also uses Jacobi’s method to solve the Poisson equation. The difference is that it precomputes a fixed number of Jacobi iterations as a convolution filter.

Implementing convolution filters on the GPU is pretty easy. The only challenge with this method is precomputing the convolution filters so that you can use them in your fluid simulator. Calculating these filters requires some information about the fluid including its viscosity, which means you have to recalculate the filters any time one of these parameters is changed.

However, from a simulation standpoint, it is possibly the best balance between accuracy and speed compared to all the other methods we are going to discuss.

Next:
References
Previous: Gpu fluids › Chapter 9.
Jos Stam Advection Scheme