It’s pointless to say. Everyone lives in a soundscape and it’s also true that I work in a field related to acoustics. And yes, this animation is also about the acoustic pressure field. I made the RGB image as the initial hydrodynamic pressure and velocity to see how the waves behave. This is how I made it (btw this is just for fun.)
The implementation conceptually follows a standard hydrodynamic/acoustic split: an incompressible base flow solves a Navier–Stokes–type system, while a linear acoustic field is driven by the time derivative of the incompressible (hydrodynamic) pressure.
For a low‑Mach, incompressible fluid in 2D with constant density \(\rho\) and kinematic viscosity \(\nu\), the governing equations for the hydrodynamic velocity \(\mathbf{u} = (u, v)\) and hydrodynamic pressure \(p_{\mathrm{fluid}}\) are
\[ \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} = -\frac{1}{\rho}\nabla p_{\mathrm{fluid}} + \nu \nabla^2 \mathbf{u} \]
\[ \nabla \cdot \mathbf{u} = 0. \]
In the code, the incompressible constraint is enforced using a projection method that solves a pressure Poisson equation for \(p_{\mathrm{fluid}}\).
The acoustic pressure \(p_{\mathrm{acoustic}}\) satisfies a linear wave equation with a source based on the time derivative of the incompressible hydrodynamic pressure, which is a common form in hydrodynamic/acoustic splitting formulations. The model equation used is
\[ \frac{\partial^2 p_{\mathrm{acoustic}}}{\partial t^2} - c_0^2 \nabla^2 p_{\mathrm{acoustic}} = -\frac{\partial}{\partial t} p_{\mathrm{fluid}}, \]
where \(c_0\) is the speed of sound. This couples the flow to the acoustic field while keeping the acoustic dynamics linear and defined on the same grid as the flow.
The implementation uses a uniform grid with spacing \(\Delta x = dx\), \(\Delta y = dy\) and explicit time-stepping with step size \(\Delta t = dt\).
All spatial derivatives are approximated by centered second‑order finite differences on a collocated grid:
First derivative in \(x\): \[ \left.\frac{\partial \phi}{\partial x}\right|_{i,j} \approx \frac{\phi_{i+1,j} - \phi_{i-1,j}}{2\Delta x} \]
First derivative in \(y\): \[ \left.\frac{\partial \phi}{\partial y}\right|_{i,j} \approx \frac{\phi_{i,j+1} - \phi_{i,j-1}}{2\Delta y} \]
Laplacian: \[ \left.\nabla^2 \phi\right|_{i,j} \approx \frac{\phi_{i+1,j} - 2\phi_{i,j} + \phi_{i-1,j}}{\Delta x^2} + \frac{\phi_{i,j+1} - 2\phi_{i,j} + \phi_{i,j-1}}{\Delta y^2}. \]
After advection and diffusion, an intermediate velocity \(\mathbf{u}^*\) is computed. Its discrete divergence at grid node \((i,j)\) is
\[ (\nabla \cdot \mathbf{u}^*)_{i,j} \approx \frac{u^*_{i+1,j} - u^*_{i-1,j}}{2\Delta x} + \frac{v^*_{i,j+1} - v^*_{i,j-1}}{2\Delta y}. \]
The projection method enforces incompressibility by solving a Poisson equation for a scalar pressure \(p_{\mathrm{prj}}\):
\[ -\nabla^2 p_{\mathrm{prj}} = \frac{\rho}{\Delta t} (\nabla \cdot \mathbf{u}^*). \]
In discrete form, using the centered Laplacian,
\[ -(\nabla^2 p_{\mathrm{prj}})_{i,j} \approx -\left[ \frac{p_{\mathrm{prj},i+1,j} - 2p_{\mathrm{prj},i,j} + p_{\mathrm{prj},i-1,j}}{\Delta x^2} + \frac{p_{\mathrm{prj},i,j+1} - 2p_{\mathrm{prj},i,j} + p_{\mathrm{prj},i,j-1}}{\Delta y^2} \right] = \frac{\rho}{\Delta t} (\nabla \cdot \mathbf{u}^*)_{i,j}. \]
The code solves this Poisson equation with a Jacobi iteration:
\[ p_{\mathrm{prj},i,j}^{(n+1)} = \frac{1}{4} \left( p_{\mathrm{prj},i+1,j}^{(n)} + p_{\mathrm{prj},i-1,j}^{(n)} + p_{\mathrm{prj},i,j+1}^{(n)} + p_{\mathrm{prj},i,j-1}^{(n)} - \Delta x^2 \frac{\rho}{\Delta t} (\nabla \cdot \mathbf{u}^*)_{i,j} \right). \]
Once \(p_{\mathrm{prj}}\) is obtained, the velocity is projected to a divergence‑free field:
\[ \mathbf{u}^{n+1} = \mathbf{u}^* - \frac{\Delta t}{\rho} \nabla p_{\mathrm{prj}}. \]
Componentwise in discrete form,
\[ u_{i,j}^{n+1} = u^*_{i,j} - \frac{\Delta t}{\rho} \frac{p_{\mathrm{prj},i+1,j} - p_{\mathrm{prj},i-1,j}}{2\Delta x}, \]
\[ v_{i,j}^{n+1} = v^*_{i,j} - \frac{\Delta t}{\rho} \frac{p_{\mathrm{prj},i,j+1} - p_{\mathrm{prj},i,j-1}}{2\Delta y}. \]
The hydrodynamic pressure \(p_{\mathrm{fluid}}\) is then updated from this projection pressure:
\[ p_{\mathrm{fluid}}^{n+1} \gets p_{\mathrm{prj}}. \]
The code keeps both \(p_{\mathrm{fluid}}^n\) and \(p_{\mathrm{fluid}}^{n-1}\) to construct a discrete time derivative for the acoustic source.
The acoustic equation
\[ \frac{\partial^2 p_{\mathrm{acoustic}}}{\partial t^2} - c_0^2 \nabla^2 p_{\mathrm{acoustic}} = -\frac{\partial}{\partial t} p_{\mathrm{fluid}} \]
is discretized with a second‑order central difference in time (leapfrog) and the same spatial Laplacian. Let \(p_{\mathrm{acoustic}}^n\) be the acoustic pressure at time \(t_n\), similarly \(p_{\mathrm{fluid}}^n\). The second time derivative is approximated as
\[ \left.\frac{\partial^2 p_{\mathrm{acoustic}}}{\partial t^2}\right|^n \approx \frac{p_{\mathrm{acoustic}}^{n+1} - 2 p_{\mathrm{acoustic}}^{n} + p_{\mathrm{acoustic}}^{n-1}}{\Delta t^2}, \]
and the source term is approximated by a backward difference in time,
\[ \left.\frac{\partial p_{\mathrm{fluid}}}{\partial t}\right|^n \approx \frac{p_{\mathrm{fluid}}^{n} - p_{\mathrm{fluid}}^{n-1}}{\Delta t}. \]
Substituting these into the PDE and solving for \(p_{\mathrm{acoustic}}^{n+1}\) yields
\[ p_{\mathrm{acoustic}}^{n+1} = 2 p_{\mathrm{acoustic}}^{n} - p_{\mathrm{acoustic}}^{n-1} + c_0^2 \Delta t^2 \nabla^2 p_{\mathrm{acoustic}}^{n} - \Delta t^2 \frac{p_{\mathrm{fluid}}^{n} - p_{\mathrm{fluid}}^{n-1}}{\Delta t}, \]
which simplifies to
\[ p_{\mathrm{acoustic}}^{n+1} = 2 p_{\mathrm{acoustic}}^{n} - p_{\mathrm{acoustic}}^{n-1} + c_0^2 \Delta t^2 \nabla^2 p_{\mathrm{acoustic}}^{n} - \Delta t (p_{\mathrm{fluid}}^{n} - p_{\mathrm{fluid}}^{n-1}). \]
The Laplacian is computed via the centered stencil and the time history \((p_{\mathrm{acoustic}}^n, p_{\mathrm{acoustic}}^{n-1}, p_{\mathrm{fluid}}^n, p_{\mathrm{fluid}}^{n-1})\) is stored in fields.
The overall algorithm per time step is:
Advection of velocity A semi‑Lagrangian–style backtracking step traces grid points backward along the current velocity field and samples the velocity to build an intermediate field \(\mathbf{u}^*\). In the code this is done using a simple nearest‑neighbor sampling of the backtracked position.
Viscous diffusion (optional) If viscosity \(\nu \neq 0\), a Laplacian term \(\nu \nabla^2 \mathbf{u}\) is applied explicitly to \(\mathbf{u}^*\).
Compute divergence of \(\mathbf{u}^*\) Evaluate the centered finite‑difference divergence of the intermediate velocity field.
Jacobi iterations for the pressure Poisson equation Perform Jacobi iterations to approximate the solution of the Poisson equation for the projection pressure \(p_{\mathrm{prj}}\). Boundary behavior is enforced via clamping indices at the domain edges, effectively approximating a zero‑normal‑gradient condition for pressure.
Projection and hydrodynamic pressure update Compute the gradient of \(p_{\mathrm{prj}}\) to correct \(\mathbf{u}^*\) into a divergence‑free velocity \(\mathbf{u}^{n+1}\) and store \(p_{\mathrm{prj}}\) into \(p_{\mathrm{fluid}}\). The previous hydrodynamic pressure \(p_{\mathrm{fluid}}^{n-1}\) is also shifted to keep a one‑step history.
Acoustic wave equation step Compute the Laplacian of \(p_{\mathrm{acoustic}}^n\), form the discrete source \(-\partial_t p_{\mathrm{fluid}}\), and advance the acoustic pressure to \(p_{\mathrm{acoustic}}^{n+1}\) using the second‑order leapfrog formula. Then store the previous step values to build the second‑order time discretization.
The solver uses an RGB image as a compact representation of the initial hydrodynamic state. Given the image of shape \((H, W, 3)\), the simulation grid is set to \((nx, ny) = (W, H)\). In the initialization kernel, each pixel at image coordinates \((i,j)\) (with \(i\) as \(x\), \(j\) as \(y\)) is mapped to the fluid fields as
The acoustic pressure \(p_{\mathrm{acoustic}}\) is initially set to zero, assuming that the hydrodynamic field defines the base flow and the acoustic part starts from rest. For output, the mapping is mirrored: \(R \leftarrow p_{\mathrm{fluid}}(i,j)\), \(G \leftarrow u(i,j)\), and \(B \leftarrow v(i,j)\).
