1. Intro Paragraph
The Marker-and-Cell (MAC) method is one of the most widely used fluid simulation methods in computer graphics. Compared with the finite volume method (FVM), MAC may feel less familiar to readers outside the graphics community, but the underlying idea is similar: discretize the Navier-Stokes equations and solve them on a grid. MAC can be used to simulate many fluids, including water and smoke. One benefit of the MAC formulation is that the pressure projection and velocity update are naturally aligned with the staggered grid layout, which makes the method easier to implement and debug.
This tutorial is structured as follows. First, we briefly derive the incompressible Navier-Stokes equations. Then we introduce the MAC grid layout. After that, we walk through each substep of one simulation time step. Finally, an interactive browser demo lets you experiment with the method directly.
2. Governing Equations
Everything starts from Newton's second law:
\[ \mathbf{F} = m\mathbf{a}. \]Here, $\mathbf{F}$ is force, $\mathbf{a}$ is acceleration, and $m$ is mass. Bold symbols denote vectors.
Suppose the velocity field $\mathbf{u}$ depends on position $\mathbf{x}$ and time $t$:
\[ \mathbf{u} = \mathbf{u}(\mathbf{x}, t). \]Taking the derivative of velocity with respect to time gives acceleration. Since $\mathbf{u}$ depends on both position and time, we need the material derivative:
\[ \frac{D\mathbf{u}}{Dt} = \frac{\partial \mathbf{u}}{\partial t} + \frac{\partial \mathbf{u}}{\partial \mathbf{x}} \cdot \frac{\partial \mathbf{x}}{\partial t}. \]At a moving material point, the derivative of position with respect to time is the local velocity:
\[ \frac{\partial \mathbf{x}}{\partial t} = \mathbf{u}. \]The spatial derivative of velocity can be written using the gradient operator:
\[ \nabla = \left( \frac{\partial}{\partial x}, \frac{\partial}{\partial y}, \frac{\partial}{\partial z} \right). \]Therefore, the material acceleration is:
\[ \frac{D\mathbf{u}}{Dt} = \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u}\cdot\nabla)\mathbf{u}. \]For a small fluid volume with density $\rho$, the force terms usually include pressure, viscosity, and external body forces. Pressure contributes $-\frac{1}{\rho}\nabla p$. Viscosity is represented by $\nu\nabla^2\mathbf{u}$, where $\nu = \mu / \rho$ is the kinematic viscosity. External body forces such as gravity or buoyancy are denoted by $\mathbf{f}$.
The incompressible Navier-Stokes equations are:
\[ \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u}\cdot\nabla)\mathbf{u} = -\frac{1}{\rho}\nabla p + \nu\nabla^2\mathbf{u} + \mathbf{f}, \]with the incompressibility constraint:
\[ \nabla\cdot\mathbf{u} = 0. \]The constraint says that the velocity field should have zero divergence. In other words, the fluid should not locally expand or compress.
3. MAC Grid Layout
MAC stands for Marker-and-Cell. The markers are particles that indicate where the fluid is. The cells form the grid on which we discretize the fluid quantities. A MAC grid is a staggered grid: pressure is stored at cell centers, while velocity components are stored at face centers. The velocity component is stored on the face normal to that direction. For example, $u_x$ is stored on x-faces, $u_y$ is stored on y-faces, and $u_z$ is stored on z-faces.
Each marker particle $p$ stores its position $\mathbf{x}_p^n$, velocity $\mathbf{u}_p^n$, and mass $m_p$.
Specifically, the velocity components are stored as:
\[ \begin{aligned} u_{i+\frac12,j,k} &: \text{x-face velocity},\\ v_{i,j+\frac12,k} &: \text{y-face velocity},\\ w_{i,j,k+\frac12} &: \text{z-face velocity}. \end{aligned} \]Pressure is stored at:
\[ p_{i,j,k}, \]where $i,j,k$ are cell indices.
3.1 Divergence at the Cell Center
The divergence at cell center $(i,j,k)$ is:
\[ \begin{aligned} (\nabla\cdot\mathbf{u})_{i,j,k} =& \frac{u_{i+\frac12,j,k}-u_{i-\frac12,j,k}}{h}\\ &+ \frac{v_{i,j+\frac12,k}-v_{i,j-\frac12,k}}{h}\\ &+ \frac{w_{i,j,k+\frac12}-w_{i,j,k-\frac12}}{h}, \end{aligned} \]where $h$ is the grid spacing.
3.2 Discrete Pressure Gradient
At an x-face, the pressure gradient can be approximated as:
\[ \left(\frac{\partial p}{\partial x}\right)_{i+\frac12,j,k} \approx \frac{p_{i+1,j,k}-p_{i,j,k}}{h}. \]The y- and z-direction gradients are defined similarly.
3.3 Interpolation
Information is transferred between marker particles and grid faces. The transfer weight can be defined using linear, quadratic, or cubic interpolation functions.
For an x-face, one possible weight is:
\[ \begin{aligned} N^u_{i+\frac12,j,k}(\mathbf{x}_p) =& w\left(\frac{x_p^x}{h}-\left(i+\frac12\right)\right)\\ &\times w\left(\frac{x_p^y}{h}-j\right)\\ &\times w\left(\frac{x_p^z}{h}-k\right). \end{aligned} \]The weights for y- and z-faces are defined in the same way, but shifted to match their face-centered locations.
4. Practical Implementation Guideline
At each time step, particle information is first transferred to the grid faces. After the grid velocities are updated, information is transferred back to the particles.
4.1 Particle-to-Grid Transfer
For x-faces, the mass accumulated at a face is:
\[ m^u_{i+\frac12,j,k} = \sum_p m_p N^u_{i+\frac12,j,k}(\mathbf{x}_p). \]The x-momentum at the same face is:
\[ (mu)_{i+\frac12,j,k} = \sum_p m_p v_p^x N^u_{i+\frac12,j,k}(\mathbf{x}_p). \]Then the face velocity is:
\[ u_{i+\frac12,j,k} = \frac{(mu)_{i+\frac12,j,k}}{m^u_{i+\frac12,j,k}}. \]The y- and z-face velocities are computed similarly.
4.2 Classify Cells
Cells are classified as air, fluid, or solid. A simple rule is: if a cell contains particles, it is a fluid cell. Solid cells can be assigned from scene geometry or boundary conditions.
4.3 Apply Body Forces
Body forces are applied to face-centered velocities explicitly. For example, gravity is applied as:
\[ \mathbf{u}^{(g)} = \mathbf{u}^n + \Delta t\,\mathbf{g}. \]If gravity acts in the negative y-direction:
\[ v^{(g)}_{i,j+\frac12,k} = v^n_{i,j+\frac12,k} - \Delta t\,g. \]4.4 Solve Viscosity
For thick fluids, viscosity caused by internal friction may be important. The viscosity equation is:
\[ \frac{\mathbf{u}^*-\mathbf{u}^{(g)}}{\Delta t} = \nu\nabla^2\mathbf{u}^*, \]where:
\[ \nu = \frac{\mu}{\rho}. \]In implicit form:
\[ \left(I-\Delta t\,\nu\nabla^2\right)\mathbf{u}^* = \mathbf{u}^{(g)}. \]This gives a linear system for the face-centered velocity components. For example, the Laplacian of an x-face velocity can be discretized as:
\[ \begin{aligned} (\nabla^2 u)_{i+\frac12,j,k} =& \frac{ u_{i+\frac32,j,k} + u_{i-\frac12,j,k} + u_{i+\frac12,j+1,k} + u_{i+\frac12,j-1,k} }{h^2}\\ &+ \frac{ u_{i+\frac12,j,k+1} + u_{i+\frac12,j,k-1} - 6u_{i+\frac12,j,k} }{h^2}. \end{aligned} \]For many simple graphics demos, viscosity can be omitted at first. Pressure projection is usually the more important step to implement correctly.
4.5 Pressure Projection
The purpose of pressure projection is to project the velocity field into a divergence-free space so the fluid volume remains approximately constant. Assume the pressure field is known. The projected velocity is:
\[ \mathbf{u}^{n+1} = \mathbf{u}^* - \frac{\Delta t}{\rho}\nabla p. \]Impose incompressibility:
\[ \nabla\cdot\mathbf{u}^{n+1}=0. \]Substitute the projected velocity:
\[ \nabla\cdot \left( \mathbf{u}^* - \frac{\Delta t}{\rho}\nabla p \right) = 0. \]Therefore:
\[ \nabla\cdot \left( \frac{1}{\rho}\nabla p \right) = \frac{1}{\Delta t} \nabla\cdot\mathbf{u}^*. \]For constant density:
\[ \nabla^2 p = \frac{\rho}{\Delta t} \nabla\cdot\mathbf{u}^*. \]For fluid cell $(i,j,k)$:
\[ \begin{aligned} & \frac{p_{i+1,j,k}-2p_{i,j,k}+p_{i-1,j,k}}{h^2}\\ &+ \frac{p_{i,j+1,k}-2p_{i,j,k}+p_{i,j-1,k}}{h^2}\\ &+ \frac{p_{i,j,k+1}-2p_{i,j,k}+p_{i,j,k-1}}{h^2}\\ &= \frac{\rho}{\Delta t} (\nabla\cdot\mathbf{u}^*)_{i,j,k}. \end{aligned} \]The right-hand side is:
\[ \begin{aligned} (\nabla\cdot\mathbf{u}^*)_{i,j,k} =& \frac{u^*_{i+\frac12,j,k}-u^*_{i-\frac12,j,k}}{h}\\ &+ \frac{v^*_{i,j+\frac12,k}-v^*_{i,j-\frac12,k}}{h}\\ &+ \frac{w^*_{i,j,k+\frac12}-w^*_{i,j,k-\frac12}}{h}. \end{aligned} \]The left-hand side involves pressures from neighboring cells, so the pressure equation is a global linear system:
\[ \mathbf{A}\mathbf{x}=\mathbf{b}. \]Solving this system gives the pressure field that makes the final velocity field divergence-free. For a simple free-surface treatment, the air pressure can be assumed to be zero.
4.6 Enforce Boundary Conditions
For a static solid wall, the normal velocity should be zero:
\[ \mathbf{u}\cdot\mathbf{n}=0, \]where $\mathbf{n}$ is the outward normal of the wall. In practice, this means setting the velocity component normal to the wall to zero on solid faces.
4.7 Grid-to-Particle Transfer
There are several ways to transfer velocity back to marker particles. Particle-in-Cell (PIC) transfer is:
\[ v_p^{PIC} = \sum_f N_f(\mathbf{x}_p)u_f^{n+1}. \]Component-wise, the x-velocity is:
\[ (v_p^{PIC})^x = \sum_{i,j,k} N^u_{i+\frac12,j,k}(\mathbf{x}_p) u^{n+1}_{i+\frac12,j,k}. \]PIC is stable but highly dissipative.
Fluid-Implicit-Particle (FLIP) transfer is less dissipative, but it can be less stable:
\[ v_p^{FLIP} = v_p^n + \sum_f N_f(\mathbf{x}_p) \left( u_f^{n+1} - u_f^n \right). \]A common practical choice is to blend PIC and FLIP:
\[ v_p^{n+1} = (1-\alpha)v_p^{PIC} + \alpha v_p^{FLIP}. \]A typical value is:
\[ \alpha \approx 0.95. \]4.8 Advect Particles
Finally, marker particle positions are updated using the interpolated velocity. A simple explicit Euler update is:
\[ \mathbf{x}_p^{n+1} = \mathbf{x}_p^n + \Delta t\,\mathbf{v}_p^{n+1}. \]5. Interactive Preview
Here is an interactive demo running a MAC fluid simulation in the browser. You can adjust the FLIP ratio, time step, particle sampling resolution, and total simulation time. Press Run to compute the simulation, then use the progress slider to inspect the generated frames.