Boo-AI — Master Artificial Intelligence by Building from Scratch

Learning Objectives

By the end of this section you will be able to:

State the continuity equation in its general form $\frac{\partial \rho}{\partial t} + \nabla\!\cdot(\rho\,\mathbf{v}) = 0$ and explain every symbol.
Derive it from scratch by doing mass bookkeeping on a tiny control volume.
Read the incompressible form $\nabla\!\cdot \mathbf{v} = 0$ as a geometric statement about the divergence of a vector field.
Verify continuity numerically with the divergence theorem, in plain Python and in PyTorch.
Diagnose whether a given flow is compressible, incompressible, a source, or a sink — by inspection.

The Bookkeeping Question

One question, no equations yet: if I draw a tiny box in the middle of a flowing river, and I count every drop of water that crosses the wall of the box per second, what must the total come to?

For an ordinary river — incompressible water, no taps, no drains — the answer is obviously zero. Every drop that enters one wall must leave through another, because the water inside the box has nowhere to hide. It cannot be compressed. It cannot be destroyed. So the books must balance.

That single sentence — the books must balance — is the continuity equation. Everything that follows is just translating it into the language of partial derivatives so we can apply it at every point of the fluid simultaneously, not just on one box.

The big idea, before any math

The continuity equation is a conservation law for mass. It says: if the density inside a region is changing, that change must be accounted for by mass crossing the region's boundary. Mass cannot appear or disappear without a flux to explain it.

Intuition: A Tiny Glass Box

Imagine an infinitesimally small transparent box sitting in a stream. Call it a control volume. Its sides are $dx$ , $dy$ , $dz$ . It does not move. Water flows freely through its faces.

Three things can happen inside this box, and only three:

1. Mass flows in

Water crosses one face heading into the box. We count this as positive inflow.

2. Mass flows out

Water crosses some face heading away. We count this as positive outflow.

3. Density inside changes

The fluid already in the box gets denser or thinner: density $\rho$ ticks up or down.

Conservation of mass says these three are not independent. The change in what is inside must equal inflow minus outflow:

Rate of mass accumulating inside the box = Mass flowing in − Mass flowing out

This is the same sentence a shopkeeper says when balancing a cash register at the end of the day. It is also the continuity equation. The only thing we are about to do is translate "mass flowing in minus out" into the calculus expression $-\nabla\!\cdot(\rho\,\mathbf{v})$ .

Flux, Divergence, and the Theorem That Links Them

Before we can write the continuity equation we need two ideas: flux (how much mass crosses a surface per second), and divergence (the local source rate of a vector field).

Flux

For a fluid moving with velocity $\mathbf{v}$ and density $\rho$ , the mass flux is the vector $\rho\,\mathbf{v}$ (units: kg/m²/s). It points in the direction mass is moving, and its magnitude says how much mass per second is crossing a unit area perpendicular to it.

The total mass crossing a surface $S$ per second is the surface integral:

\dot{m} \;=\; \iint_S \rho\,\mathbf{v}\cdot \hat{\mathbf{n}} \; dS

where $\hat{\mathbf{n}}$ is the outward unit normal. The dot product picks out only the component of velocity that actually crosses the surface — fluid moving along $S$ doesn't count.

Divergence

The divergence of a vector field at a point measures the net outflow rate per unit volume at that point:

\nabla\!\cdot\mathbf{F} \;=\; \frac{\partial F_x}{\partial x} + \frac{\partial F_y}{\partial y} + \frac{\partial F_z}{\partial z}

Picture a tiny sphere around the point. If $\nabla\!\cdot\mathbf{F} > 0$ , more is flowing out than in — the point acts like a source. If $\nabla\!\cdot\mathbf{F} < 0$ , more is flowing in — a sink. If zero, the books balance locally.

The Divergence Theorem

The bridge between the two is the divergence theorem of Gauss:

\iint_{\partial V} \mathbf{F}\cdot \hat{\mathbf{n}} \; dS \;=\; \iiint_V \nabla\!\cdot \mathbf{F} \; dV

The total flux through the boundary of any region equals the integral of the divergence over the region. This is the multivariable cousin of the Fundamental Theorem of Calculus, and it is the engine that turns "mass bookkeeping on a box" into a differential equation at every point.

Why the divergence theorem is the key. The bookkeeping statement "net inflow = accumulation" is naturally written as an integral over a region. The divergence theorem lets us shrink the region to a point — turning an integral law into a pointwise differential equation. Without it, we'd only ever have the continuity equation in bulk, not at every individual point of the fluid.

Deriving the Continuity Equation

We now translate "rate of mass accumulating = mass flowing in minus out" into symbols. Take a fixed region $V$ of space, with boundary $\partial V$ .

Left-hand side: rate of mass accumulating

The total mass inside $V$ at time $t$ is the volume integral of density:

M(t) \;=\; \iiint_V \rho(\mathbf{x},t) \; dV

The rate at which that mass changes is just its time derivative:

\frac{dM}{dt} \;=\; \iiint_V \frac{\partial \rho}{\partial t}\, dV

We pulled the derivative inside the integral because $V$ is fixed in space — it doesn't move with the fluid, so its limits of integration don't depend on $t$ .

Right-hand side: net mass entering

The mass crossing the boundary $\partial V$ per second is the surface integral of mass flux. Outflow is positive when the velocity points along the outward normal. So the net inflow is the negative of the outflow:

\dot{M}_{\text{in}} \;=\; -\iint_{\partial V} \rho\,\mathbf{v} \cdot \hat{\mathbf{n}} \; dS

Apply the divergence theorem

The right-hand side is begging for the divergence theorem applied to the vector field $\mathbf{F} = \rho\,\mathbf{v}$ :

\iint_{\partial V} \rho\,\mathbf{v} \cdot \hat{\mathbf{n}} \; dS \;=\; \iiint_V \nabla\!\cdot(\rho\,\mathbf{v}) \; dV

Set accumulation equal to inflow

Combining the last three statements, $dM/dt = \dot{M}_{\text{in}}$ becomes:

\iiint_V \frac{\partial \rho}{\partial t}\, dV \;=\; -\iiint_V \nabla\!\cdot(\rho\,\mathbf{v})\, dV

or, moving both integrands to one side,

\iiint_V \left[\, \frac{\partial \rho}{\partial t} + \nabla\!\cdot(\rho\,\mathbf{v}) \,\right]\, dV \;=\; 0.

Shrink the region to a point

The crucial step. This integral identity must hold for any choice of $V$ — any shape, any size, anywhere in the fluid. The only way that's possible is if the integrand is itself zero everywhere. (If it were positive at a point, we'd pick a tiny ball around that point and break the equation.) That gives us the differential form:

The Continuity Equation

\boxed{\;\dfrac{\partial \rho}{\partial t} \;+\; \nabla\!\cdot(\rho\,\mathbf{v}) \;=\; 0\;}

Local rate of density change + Net outflow of mass per unit volume = 0

Read it left to right: if the density at a point is going up, the divergence of mass flux at that point must be negative — mass must be piling in faster than it leaves. If density is going down, mass is escaping. If density is steady, the books balance pointwise.

The Incompressible Form

For liquids like water, density barely changes — call it constant. Plugging in $\rho = \text{const}$ :

The time-derivative term vanishes: $\partial\rho/\partial t = 0$ .
The constant slides out of the divergence: $\nabla\!\cdot(\rho \mathbf{v}) = \rho\,\nabla\!\cdot\mathbf{v}$ .
Dividing through by the (non-zero) constant $\rho$ leaves a beautifully clean equation.

The Incompressible Continuity Equation

\boxed{\;\nabla\!\cdot \mathbf{v} \;=\; 0\;}

For water, blood, oil at modest speeds — the velocity field is divergence-free.

Geometrically: at every point in the fluid, the velocity field has no sources and no sinks. Whatever flows in must flow out — at every point, all at once. This is one of the cleanest equations in physics, and it constrains the velocity field of a Navier-Stokes flow in a profound way: of the three velocity components $(u, v, w)$ , only two are actually independent, because $\nabla\!\cdot\mathbf{v} = 0$ ties them together.

When you cannot assume incompressibility

For gases at high speed (Mach > 0.3), shock waves, sound waves, atmospheric science, plasma physics, and combustion, density changes matter and you must keep the full $\partial\rho/\partial t$ term. The incompressible shortcut is a wonderful, almost-universally-correct approximation for liquids and slow gas flows — but it is a shortcut, not a law of physics.

Interactive: Watch Mass Flow Through a Box

Drag the dashed amber box anywhere in the field. Move the divergence slider $\alpha$ and the rotation slider $\omega$ . The four numbers on the box show the flux through each face — red is outward, blue is inward — and the panel on the right keeps a running ledger.

Loading interactive flow explorer…

Try these experiments — they make the math concrete:

Click Pure Rotation. Watch every face have non-zero flux, yet the net out always pinches back to zero. That is the continuity equation for $\nabla\!\cdot\mathbf{v} = 0$ .
Click Source. Notice every face now has positive outflux (red). The net is positive, equal to $2\alpha \cdot \text{Area}$ — exactly what the divergence theorem predicts. Compare the "NET OUT" row to the "2α · Area" row in the panel.
Drag the box around in the source case. Net outflux changes with where you put the box? No — only with the box's area. That's because divergence is uniform here.
Set Spiral (Mixed). The flow looks complicated but the bookkeeping is the same: net out = 2α · Area, independent of ω. Rotation does not violate continuity.

Worked Example: Spinning vs Exploding Fluids

Let's compute the divergence (and net outflux through a unit square) of two specific fields by hand. They look similar, but one is incompressible and the other isn't.

Field	Components	What it looks like
Field A (Rotation)	$\mathbf{v}_A = (-y,\; x)$	Counter-clockwise solid-body rotation about the origin
Field B (Outflow)	$\mathbf{v}_B = (x,\; y)$	Radial expansion away from the origin

We'll compute $\nabla\!\cdot\mathbf{v}$ analytically, then verify with a surface integral on the unit square $[0,1]\times[0,1]$ . Try it yourself first before expanding the solution.

Expand: full hand calculation (try it yourself first!)

Field A: rotation

Step 1 — divergence. With $u_A = -y$ and $v_A = x$ :

\nabla\!\cdot\mathbf{v}_A = \frac{\partial(-y)}{\partial x} + \frac{\partial(x)}{\partial y} = 0 + 0 = 0.

The field is everywhere divergence-free. The continuity equation is satisfied at every single point.

Step 2 — flux through the unit square. The outward normal on the right edge (x = 1) is $\hat{\mathbf{n}} = (1,0)$ , so $\mathbf{v}_A \cdot \hat{\mathbf{n}} = -y$ . Integrate over $0 \le y \le 1$ :

\Phi_{\text{right}} = \int_0^1 (-y)\, dy = -\tfrac{1}{2}.

By similar work on the other three faces:

Left face (x = 0, normal = (−1, 0)): $\Phi_{\text{left}} = \int_0^1 -(-y)\, dy = +\tfrac{1}{2}$ — but wait, on x = 0 the field component u = −y still, so this is $\int_0^1 -(-y)\, dy = \tfrac{1}{2}$ .
Top face (y = 1, normal = (0, 1)): $\Phi_{\text{top}} = \int_0^1 x\, dx = +\tfrac{1}{2}$ .
Bottom face (y = 0, normal = (0, −1)): $\Phi_{\text{bot}} = \int_0^1 -x\, dx = -\tfrac{1}{2}$ .

Net: $-\tfrac{1}{2} + \tfrac{1}{2} + \tfrac{1}{2} - \tfrac{1}{2} = 0$ . The fluid spins around but the books are balanced — consistent with the pointwise divergence being zero.

Field B: outflow

Step 1 — divergence. Now $u_B = x$ , $v_B = y$ :

\nabla\!\cdot\mathbf{v}_B = \frac{\partial x}{\partial x} + \frac{\partial y}{\partial y} = 1 + 1 = 2.

Divergence is the constant 2 everywhere. Every point of space is acting as a uniform mass source.

Step 2 — flux through the unit square.

Right face (x = 1): $\Phi_{\text{right}} = \int_0^1 1\, dy = 1$ .
Left face (x = 0): $\Phi_{\text{left}} = \int_0^1 -0\, dy = 0$ .
Top face (y = 1): $\Phi_{\text{top}} = \int_0^1 1\, dx = 1$ .
Bottom face (y = 0): $\Phi_{\text{bot}} = \int_0^1 -0\, dx = 0$ .

Net: $1 + 0 + 1 + 0 = 2$ .

Check against divergence theorem: $\iint_{V} \nabla\!\cdot\mathbf{v}_B \, dA = \iint_{[0,1]^2} 2\, dA = 2\cdot 1 = 2.$ The two answers match exactly.

The Punchline

Field A passes the continuity test (∇·v = 0). Field B fails it uniformly — every point is a source of mass at rate 2. If you saw Field B in a real flow, you'd know mass is being created from nothing somewhere, which is physically impossible: in a real fluid you'd have to add a source term (a tap, an injector, a chemical reaction) to the right-hand side of the continuity equation.

Plain Python: Checking Continuity Numerically

Algebra is reassuring. Numbers are more reassuring. We'll now compute both sides of the divergence theorem — the pointwise divergence and the net flux through the boundary — in plain Python with finite differences, and confirm they match for two test fields.

What we're about to do. Build the velocity field

\mathbf{v} = (\alpha x - \omega y,\; \alpha y + \omega x)

, compute its divergence numerically at a point, compute the net outflux through a 1.8 × 1.8 box numerically by walking around the four sides, and watch the two numbers agree to many decimals.

Continuity Equation Verified by Finite Differences

🐍python

Explanation(43)

Code(76)

1Import NumPy

We will use NumPy only for cleaner numerical types. The actual logic is plain Python loops — that is on purpose, because seeing the four sides of the box as four separate loops is exactly what teaches continuity.

9Define the velocity field

We pick a two-parameter family v(x,y) = (αx − ωy, αy + ωx) because we can hand-compute its divergence and curl: ∂u/∂x = α, ∂v/∂y = α, so div v = 2α; and ∂v/∂x − ∂u/∂y = ω + ω = 2ω. Setting α = 0 turns the field into pure rotation — exactly the test case the continuity equation should pass.

EXAMPLE

u = αx − ωy
v = αy + ωx

10Compute u (x-component)

u is the x-component of the velocity. At the point (2, 1) with α = 0.3 and ω = 0.7 we get u = 0.3·2 − 0.7·1 = 0.6 − 0.7 = −0.1. The first term grows the field outward from the origin; the second term twists it.

11Compute v (y-component)

v is the y-component. At the same point: v = 0.3·1 + 0.7·2 = 0.3 + 1.4 = 1.7. Notice u and v share the same α (so divergence is symmetric in both directions) and the same ω (so rotation is rigid).

12Return the velocity as a 2-tuple

Returning a tuple keeps the code language-agnostic (no NumPy arrays needed). When we sample the four sides of the box below, we will only ever look at one component at a time — u on vertical sides, v on horizontal sides.

17Central-difference divergence

This is a direct numerical check of ∇·v = ∂u/∂x + ∂v/∂y at a single point. The continuity equation says this should be zero everywhere for incompressible flow — so we will use this to verify our analytical claim (div = 2α) without trusting symbolic math.

18Sample u to the right

We evaluate the x-component slightly to the right of the target point: u(x + h, y). For h = 10⁻⁴ at (1.5, 0.5) with α = 0.3, ω = 0, this gives u_xp = 0.3 · 1.5001 = 0.450030. The trailing _ discards the v-component because we only need u here.

19Sample u to the left

Same x-component, but slightly to the left: u(x − h, y). For the same parameters: u_xm = 0.3 · 1.4999 = 0.449970. The pair (u_xp, u_xm) brackets the true x-derivative at the centre.

20Sample v above

Now we need the y-derivative of v, so we move in y by +h and read off the y-component: v(x, y + h). The leading _ throws away the u-component because central differences in y only need v.

21Sample v below

Mirror of the previous line: evaluate v at y − h. Combined with line 20 these two values let us approximate ∂v/∂y.

22Approximate ∂u/∂x

Central difference: (u(x+h) − u(x−h)) / (2h). For our incompressible case (α = 0) both u_xp and u_xm equal −ω·y so du_dx ≈ 0 — exactly what calculus predicts. The error of a central difference is O(h²), much better than the O(h) error of a one-sided difference.

EXAMPLE

du_dx ≈ (u(x+h) - u(x-h)) / (2h)

23Approximate ∂v/∂y

Same idea in the y-direction. We deliberately do not mix axes — du_dx uses only u, dv_dy uses only v. The divergence is then the sum of these two slopes.

24Return the divergence sum

∇·v = ∂u/∂x + ∂v/∂y. For α = 0 we expect 0; for α = 0.3 we expect 0.6. Comparing this single number with the volume integral computed below is the heart of the divergence theorem.

32Net outward flux through a square

This is the surface-integral side of the divergence theorem. Instead of evaluating ∇·v at a point, we add up how much mass actually crosses each face of a square. The continuity equation says these two quantities — the local divergence and the net outflux — must match (up to factor of area).

33x-range of the box

We centre the box at (xc, yc) and give it half-width `half`. Then a = xc − half is the left edge and b = xc + half is the right edge. For xc = 1.5, half = 0.9 we get a = 0.6, b = 2.4.

34y-range of the box

Same construction in the y-direction. Together a, b, c, d are the four corners of the control volume — the region we are doing the mass bookkeeping over.

35x step size for the midpoint rule

We will integrate along each side using N rectangles of equal width. dx = (b − a) / N is the step size in x. With N = 200 and a length of 1.8, dx = 0.009. The error of the midpoint rule is O(dx²), so doubling N quarters the error.

36y step size

Same idea for the y-direction. dx and dy are independent here — for a square box they happen to be equal, but the code does not assume that, which makes it reusable for rectangles too.

38Running total of flux

We accumulate the contribution of each face into one number. We use the convention: positive = mass leaving the box, negative = mass entering. The continuity equation says this total should equal ∇·v · Area for any choice of box.

40Right face: outward normal is +x̂

On the right edge the outward-pointing unit normal is (+1, 0). The dot product v · n is therefore just u, the x-component of velocity. We sample u at the midpoint of each strip along the right edge.

41Midpoint of each strip on the right face

For the midpoint rule we sample at the centre of each strip: y_mid = c + (i + 0.5)·dy. With c = −0.4, dy = 0.009, i = 0 gives y_mid = −0.3955. The midpoint rule converges twice as fast as left- or right-Riemann.

42Read u on the right edge

We evaluate the velocity at (b, y_mid) — i.e., on the right edge at vertical position y_mid. We keep only u because the outward normal is along x̂. Concretely u = α·b − ω·y_mid; for α = 0, ω = 0.7, b = 2.4, y_mid = 0 we get u = −0.7·0 = 0 — mass moves vertically there, not horizontally.

43Add u·dy to total

Each strip contributes u(b, y_mid)·dy to the right-face flux. Summing all N strips approximates ∫ u(b, y) dy from c to d. For our incompressible α = 0 case this whole sum will cancel against the left face — but the cancellation is non-obvious until you watch it happen.

45Left face: outward normal is −x̂

On the left edge the outward normal points in the −x direction, so v · n = −u. That single minus sign is the reason the right and left faces can cancel each other out for incompressible flow.

46Midpoint of each strip on the left face

Same midpoint rule, same y-grid. We re-use the (i + 0.5)·dy construction so the right and left faces are sampled at the exact same y values — that is what guarantees clean cancellation when the field is divergence-free.

47Read u on the left edge

Now we evaluate u at x = a (the left edge), not b. For α = 0 we get u(a, y_mid) = −ω·y_mid, exactly mirroring u(b, y_mid) — so when we add −u(a) to the total it cancels +u(b) strip by strip.

48Add −u·dy to total

Notice the minus sign. We are explicitly applying v · n with n = −x̂. After the right and left loops finish, we will have computed ∫₍right₎ u dy − ∫₍left₎ u dy = ∫∫ ∂u/∂x dA — which is half of the divergence theorem in 2D.

50Top face: outward normal is +ŷ

On the top edge the outward normal is (0, +1), so v · n = v (the y-component). We now sample along the x-direction at fixed y = d.

51Midpoint of each strip on the top face

x_mid = a + (i + 0.5)·dx. With a = 0.6, dx = 0.009, i = 0 we get x_mid = 0.6045. We are now sweeping horizontally across the top of the box.

52Read v on the top edge

Velocity at (x_mid, d). We keep only v because the outward normal is along ŷ. The pattern mirrors the right face exactly — that symmetry is why the divergence theorem looks so clean.

53Add v·dx to total

Contribution of one strip on the top face. Summing all N strips approximates ∫ v(x, d) dx from a to b.

55Bottom face: outward normal is −ŷ

Last of the four faces. Outward normal points in −y direction, so v · n = −v. Together with the top face this will compute ∫∫ ∂v/∂y dA — the other half of the divergence theorem.

56Midpoint along the bottom face

Same x-grid as the top face. Identical sampling guarantees that for incompressible flow the top and bottom contributions cancel cleanly.

57Read v on the bottom edge

Velocity at (x_mid, c). We are now on y = c, the bottom edge. For pure rotation v(x, c) = ω·x — and v on the top edge equals ω·x at the same x_mid, so the difference (top − bottom) is exactly zero strip by strip.

58Add −v·dx to total

Closes the loop. After all four loops, `total` holds the line integral of v · n over the boundary — i.e., the net mass leaving the box per unit time. The continuity equation predicts this equals 2α · Area, independent of where the box is.

60Return the net flux

One number summarising every side. If positive, the box is a net source; if negative, a net sink; if (approximately) zero, the flow is incompressible inside the box.

67Test 1: incompressible case

We set α = 0, ω = 0.7. Analytically ∇·v = 2·0 = 0 everywhere. We expect both the pointwise divergence and the net outflux to be zero (up to numerical noise of order 10⁻¹⁰).

68Direct numerical flux

We integrate around the box. For pure rotation, every drop of fluid that enters through the right re-emerges on the left, every drop that goes up through the top comes back down through the bottom — net is zero. This is the continuity equation in action.

69Print incompressible result

Both numbers should print as essentially 0 (perhaps 1e-14 from rounding). That single agreement is more convincing than any algebraic proof, because we computed the two sides of the divergence theorem independently and they matched.

EXAMPLE

Incompressible:  div = 0.0   net out ≈ 0.0

72Test 2: source case

Now α = 0.3, ω = 0. Analytically ∇·v = 2·0.3 = 0.6 at every point, meaning the fluid is being created uniformly. The continuity equation predicts net outflux = div · Area.

73Numerical net flux for the source

We expect this to equal 0.6 · 3.24 = 1.944. Run the code and watch it land within a few parts in a million of that value — the leftover is the O(dx²) error of the midpoint rule.

74Hand-computed area

Area of the box = (2 · half)² = 1.8 · 1.8 = 3.24. We separate it out so the printed line directly shows that div × area matches the surface integral — that is the divergence theorem in one numerical line.

75Print source result

When you run this you see three numbers in a row: div = 0.6, net out ≈ 1.944, div × area = 1.944. The fact that the second and third agree to many decimals is the numerical witness for the continuity equation.

EXAMPLE

Source:  div = 0.6   net out = 1.94399...   div * area = 1.944

33 lines without explanation

1import numpy as np
2
3# ---------------------------------------------------------
4# A 2D velocity field. We will parameterize a family:
5#   v(x, y) = (alpha*x - omega*y,  alpha*y + omega*x)
6# Divergence: div v = 2*alpha     (creation / destruction)
7# Curl:       curl v = 2*omega    (rotation, no div effect)
8# ---------------------------------------------------------
9def velocity(x, y, alpha, omega):
10    u = alpha * x - omega * y
11    v = alpha * y + omega * x
12    return u, v
13
14# ---------------------------------------------------------
15# 1. Numerical divergence at a point, by central differences
16# ---------------------------------------------------------
17def divergence_at(x, y, alpha, omega, h=1e-4):
18    u_xp, _ = velocity(x + h, y, alpha, omega)
19    u_xm, _ = velocity(x - h, y, alpha, omega)
20    _, v_yp = velocity(x, y + h, alpha, omega)
21    _, v_ym = velocity(x, y - h, alpha, omega)
22    du_dx = (u_xp - u_xm) / (2 * h)
23    dv_dy = (v_yp - v_ym) / (2 * h)
24    return du_dx + dv_dy
25
26# ---------------------------------------------------------
27# 2. Net outward flux through a square, by direct integration
28#    along the four sides. This is the LHS of the divergence
29#    theorem:   \oint  v . n  dS
30# ---------------------------------------------------------
31def net_outflux(xc, yc, half, alpha, omega, N=200):
32    a, b = xc - half, xc + half   # x-range
33    c, d = yc - half, yc + half   # y-range
34    dx = (b - a) / N
35    dy = (d - c) / N
36
37    total = 0.0
38    # Right side: outward normal = +x_hat, integrand = u
39    for i in range(N):
40        y_mid = c + (i + 0.5) * dy
41        u, _ = velocity(b, y_mid, alpha, omega)
42        total += u * dy
43    # Left side: outward normal = -x_hat, integrand = -u
44    for i in range(N):
45        y_mid = c + (i + 0.5) * dy
46        u, _ = velocity(a, y_mid, alpha, omega)
47        total += -u * dy
48    # Top side: outward normal = +y_hat, integrand = v
49    for i in range(N):
50        x_mid = a + (i + 0.5) * dx
51        _, v = velocity(x_mid, d, alpha, omega)
52        total += v * dx
53    # Bottom side: outward normal = -y_hat, integrand = -v
54    for i in range(N):
55        x_mid = a + (i + 0.5) * dx
56        _, v = velocity(x_mid, c, alpha, omega)
57        total += -v * dx
58
59    return total
60
61# ---------------------------------------------------------
62# 3. Run the experiment and confirm the divergence theorem
63# ---------------------------------------------------------
64if __name__ == "__main__":
65    # Incompressible test: alpha = 0 -> div = 0 -> net outflux must be 0
66    inc_div = divergence_at(1.5, 0.5, alpha=0.0, omega=0.7)
67    inc_out = net_outflux(1.5, 0.5, 0.9, alpha=0.0, omega=0.7)
68    print("Incompressible:  div =", inc_div, "  net out =", inc_out)
69
70    # Source test: alpha = 0.3 -> div = 0.6 -> net outflux = 0.6 * area
71    src_div  = divergence_at(1.5, 0.5, alpha=0.3, omega=0.0)
72    src_out  = net_outflux(1.5, 0.5, 0.9, alpha=0.3, omega=0.0)
73    area = (2 * 0.9) ** 2          # = 3.24
74    print("Source:          div =", src_div,
75          "  net out =", src_out,
76          "  div * area =", src_div * area)

What you should see when you run this. For the incompressible case (α = 0), both numbers are essentially zero (≈ 10⁻¹⁴). For the source case (α = 0.3), the pointwise divergence is 0.6, the net outflux is ≈ 1.944, and 0.6 × 3.24 = 1.944 — three independent numbers locking together to confirm the divergence theorem.

PyTorch: Continuity via Autograd

Finite differences are honest but noisy. Autograd is exact (to floating-point precision) and scales effortlessly. Better still, in 2D we can construct a divergence-free field automatically using a stream function: pick any scalar ψ(x, y), set $\mathbf{v} = (\partial \psi/\partial y,\; -\partial \psi/\partial x)$ , and the divergence is identically zero — guaranteed by the equality of mixed partials.

The plan. Define a stream function. Get the velocity by differentiating it. Then differentiate that to compute the divergence. PyTorch threads the whole chain rule automatically. Compare side by side with a deliberately non-divergence-free radial outflow field.

Divergence-free Flow from a Stream Function

🐍python

Explanation(33)

Code(61)

1Import PyTorch

Plain Python made the divergence concrete with finite differences. PyTorch lets us go further: we will construct a field whose divergence is zero by construction, and let autograd verify it for us at every grid point.

9Stream function ψ(x, y)

In 2D incompressible flow there is a remarkable trick: any smooth scalar field ψ(x, y) generates a divergence-free velocity field via v = (∂ψ/∂y, −∂ψ/∂x). This is because mixed partials commute: ∂²ψ/∂x∂y = ∂²ψ/∂y∂x, so the two terms in ∇·v cancel identically. The continuity equation is built into the math, not enforced at runtime.

EXAMPLE

v = ( ∂ψ/∂y , -∂ψ/∂x )    ⇒   ∇·v ≡ 0

11Choose a stream function

We pick two Gaussian bumps of opposite sign, centred at x = 1 and x = -1. This creates two counter-rotating vortices side by side — a classic 2D ideal-fluid example. The exact shape doesn't matter for the divergence test; what matters is that any twice-differentiable ψ gives a divergence-free velocity.

14Define the velocity v(x, y)

We take x and y as tensors with `requires_grad=True` upstream, build ψ from them, and ask autograd to differentiate. This is the same machinery that trains neural networks — but here we're using it to do calculus on a fluid.

16Compute ψ

We need a scalar to call .grad on. Because x and y are grids, ψ is a tensor of shape (21, 21). We will reduce to a scalar with .sum() inside the grad call. Autograd is smart enough to keep track of the chain rule even through the sum.

17Differentiate ψ with respect to x and y

`torch.autograd.grad(psi.sum(), (x, y), create_graph=True)` returns the pair (∂(Σψ)/∂x, ∂(Σψ)/∂y). Because Σ is just a sum of independent contributions at each grid point, this reduces to the elementwise partial derivatives we want. `create_graph=True` is essential — without it we couldn't differentiate again later to compute the divergence.

18Unpack the two derivatives

`grad` is a tuple of two tensors, each shape (21, 21). We give them the meaningful names ∂ψ/∂x and ∂ψ/∂y. From here we just plug into the stream-function definition of velocity.

19u = ∂ψ/∂y

The x-component of velocity is the y-derivative of the stream function. Geometrically, where ψ has steep contours in the y direction, the fluid moves fast in the x direction — streamlines are level sets of ψ.

20v = −∂ψ/∂x

The minus sign is what makes the field divergence-free. If we had v = +∂ψ/∂x instead, the mixed partials would add up instead of cancelling, and we'd get back to a field with non-zero divergence in general.

23Compute the divergence using autograd

This is the test. We re-differentiate u with respect to x and v with respect to y, then add. If our construction is correct, every entry of the result is zero to floating-point precision — and PyTorch will confirm that for us at 441 grid points simultaneously.

24Velocity at the grid

We get u and v as tensors of shape (21, 21). Each entry is the velocity component at one grid node.

25∂u/∂x

We index `[0]` because `torch.autograd.grad` always returns a tuple, even with a single output. `create_graph=True` keeps the computational graph alive in case we want to take a third derivative — useful for vorticity dynamics, not used here.

26∂v/∂y

The other half of the divergence. Note that we differentiate u with respect to x and v with respect to y — never crossing the axes. That ‘matching’ is exactly what the divergence operator does.

27Return the divergence sum

∇·v = ∂u/∂x + ∂v/∂y. For our stream-function field this should be (numerically) zero. For a different field it gives us a direct read-out of how much mass is being created or destroyed at each point.

33A field that is NOT incompressible

For contrast, we build a radial outflow: v = α·(x, y). Every fluid particle moves directly away from the origin. ∂u/∂x = α and ∂v/∂y = α, so the divergence is the constant 2α. The continuity equation is not satisfied — mass is being created uniformly.

34u = α·x

Pure x-stretching with rate α. At α = 0.5 the velocity at x = 3 is 1.5 units/sec outward.

35v = α·y

Symmetric stretching in y. Together these two lines define a uniformly expanding flow — every point moves away from the origin at speed proportional to its distance.

39Divergence of the source field

Same routine as before, but now applied to a field with built-in non-zero divergence. Comparing the two outputs side by side is the punchline: one is zero everywhere, the other is uniformly 2α.

40Get u, v

Evaluate the radial outflow at every grid point. u and v each have shape (21, 21).

41∂u/∂x = α (analytically)

Analytically ∂u/∂x = ∂(αx)/∂x = α. Autograd will return a tensor of shape (21, 21) where every entry equals α. We will confirm this in the print statement.

42∂v/∂y = α (analytically)

Symmetric. So the total divergence is 2α everywhere.

43Return the (non-zero) divergence

Tensor of shape (21, 21), every entry approximately 2α. This is what an incompressibility check would flag.

49Build a 21×21 grid in x

`torch.linspace(-2, 2, 21)` creates 21 evenly spaced points from −2 to +2. `requires_grad=True` is what tells autograd to track gradients with respect to these inputs.

50Same in y

Identical construction in y. The grid has 21·21 = 441 nodes, plenty for a smooth visualisation and enough to test divergence-freeness convincingly.

51Make a 2D meshgrid

`indexing="xy"` makes X[i, j] = xs[j] and Y[i, j] = ys[i]. This is the Cartesian convention used by NumPy and matplotlib — the alternative `"ij"` swaps the two and is the source of countless silent bugs in scientific code.

53Divergence of the stream-function field

Call the function we built. PyTorch silently runs the whole chain — ψ → ∇ψ → v → ∇·v — and gives us back a (21, 21) tensor.

54Divergence of the radial source

Same call shape, different field. After this line we have two tensors of identical shape and very different content — one effectively zero, one effectively 2α.

56Print stream-function results

Two summary statistics: the maximum absolute divergence anywhere on the grid, and the mean. Both should be ≈ 1e-7 or smaller, which is the floating-point noise of single-precision differentiation. Mathematically zero; computationally negligible.

57Max |div v| for the stream field

`.detach()` removes the tensor from the autograd graph (we don't want to differentiate the print statement). `.abs().max()` takes the largest absolute value over all 441 points. `.item()` extracts a Python float. The whole point is one number.

EXAMPLE

max |div v| ≈ 6e-8   (numerical zero)

58Mean div v for the stream field

Mean is even smaller because positive and negative noise cancels. If our construction were wrong this would be the first sign — a non-zero mean would mean the field has a net source or sink.

60Print source results

Now the same diagnostics for the radial outflow. The max and mean should both be very close to 2α = 1.0 (and equal to each other, because the divergence is the same constant at every point).

61Max |div v| for the source

Equals 2α = 1.0 (up to floating-point noise). This is the loudest possible statement that the field is not incompressible — every point of the grid registers the same positive divergence.

62Mean div v matches 2α

Mean = max here because the divergence is uniform. Side by side with the stream-function field, the contrast is striking: one field has divergence floating around 1e-8 everywhere, the other has divergence locked at exactly 1.0 everywhere. That is what the continuity equation tests.

EXAMPLE

Radial:  max |div v| = 1.0   mean div v = 1.0   (expected 2α = 1.0)

28 lines without explanation

1import torch
2
3# ---------------------------------------------------------
4# A divergence-free field, built from a scalar stream function.
5# In 2D, if  v = ( d psi/d y, -d psi/d x ),  then
6#   div v = d^2 psi/(dx dy) - d^2 psi/(dy dx) = 0    automatically.
7# We let autograd handle every derivative.
8# ---------------------------------------------------------
9def stream_function(x, y):
10    # A vortex pair: two opposite swirls offset along x.
11    return torch.exp(-((x - 1)**2 + y**2)) - torch.exp(-((x + 1)**2 + y**2))
12
13def velocity(x, y):
14    # v = (d psi/d y, -d psi/d x)
15    psi = stream_function(x, y)
16    grad = torch.autograd.grad(psi.sum(), (x, y), create_graph=True)
17    dpsi_dx, dpsi_dy = grad
18    u =  dpsi_dy
19    v = -dpsi_dx
20    return u, v
21
22def divergence(x, y):
23    u, v = velocity(x, y)
24    du_dx = torch.autograd.grad(u.sum(), x, create_graph=True)[0]
25    dv_dy = torch.autograd.grad(v.sum(), y, create_graph=True)[0]
26    return du_dx + dv_dy
27
28# ---------------------------------------------------------
29# Compare to a NOT-divergence-free field, for contrast
30# ---------------------------------------------------------
31def velocity_source(x, y, alpha=0.5):
32    # Radial outflow: v = alpha * (x, y). div = 2*alpha != 0.
33    u = alpha * x
34    v = alpha * y
35    return u, v
36
37def divergence_source(x, y, alpha=0.5):
38    u, v = velocity_source(x, y, alpha)
39    du_dx = torch.autograd.grad(u.sum(), x, create_graph=True)[0]
40    dv_dy = torch.autograd.grad(v.sum(), y, create_graph=True)[0]
41    return du_dx + dv_dy
42
43# ---------------------------------------------------------
44# Evaluate both fields on a grid and confirm continuity
45# ---------------------------------------------------------
46if __name__ == "__main__":
47    xs = torch.linspace(-2, 2, 21, requires_grad=True)
48    ys = torch.linspace(-2, 2, 21, requires_grad=True)
49    X, Y = torch.meshgrid(xs, ys, indexing="xy")
50
51    div_psi    = divergence(X, Y)
52    div_source = divergence_source(X, Y, alpha=0.5)
53
54    print("Stream-function field:")
55    print("  max |div v| =", div_psi.detach().abs().max().item())
56    print("  mean div v  =", div_psi.detach().mean().item())
57
58    print("Radial source field:")
59    print("  max |div v| =", div_source.detach().abs().max().item())
60    print("  mean div v  =", div_source.detach().mean().item(),
61          "  (expected 2*alpha = 1.0)")

The deep connection. The stream-function trick is the same idea behind the magnetic vector potential in electromagnetism, the Hamiltonian formulation in classical mechanics, and the "divergence-free normalizing flows" in modern generative models — wherever you need a constraint like

\nabla\!\cdot\mathbf{v} = 0

to hold by construction, you parameterize by a scalar potential and let differentiation do the work.

What the Equation Really Says

Putting it all together: the continuity equation is the local form of mass conservation. It bridges two scales:

Integral (bulk) form

\frac{d}{dt}\!\!\iiint_V \rho\, dV = -\!\!\iint_{\partial V} \rho\,\mathbf{v}\cdot \hat{\mathbf{n}}\, dS

"The rate of mass accumulating inside a region equals the rate mass flows in across its boundary."

Differential (pointwise) form

\frac{\partial \rho}{\partial t} + \nabla\!\cdot(\rho\,\mathbf{v}) = 0

"At every point, the rate density changes is balanced by the divergence of mass flux."

These are the same statement, written at two different zoom levels. The divergence theorem is the optical instrument that switches between them.

Three ways to read the incompressible form

Reading	Statement	Useful when...
Volumetric	Fluid parcels conserve their volume as they move.	Tracking individual fluid elements (Lagrangian view).
Geometric	The velocity field has no sources or sinks.	Visualizing streamlines or doing CFD post-processing.
Algebraic	u_x + v_y + w_z = 0 ties the three velocity components.	Numerically integrating the Navier-Stokes equations.

Applications

🚿 Pipe flow

For incompressible flow in a pipe, continuity says A₁v₁ = A₂v₂. A pipe that narrows must speed the fluid up — the basis for nozzles, venturis, and your garden hose.

🌬️ Atmospheric science

Weather models cannot use the incompressible shortcut. The full continuity equation links wind, temperature, and pressure changes — it's what propagates pressure systems across continents.

❤️ Hemodynamics

Blood flow through arteries is treated as incompressible. When an artery narrows (stenosis), continuity demands the velocity rises — and ultrasound exploits this to see the blockage.

⚡ Electromagnetism

Maxwell's ∇·B = 0 is literally a continuity equation for magnetic field lines — there are no magnetic monopoles, so field lines never start or end. Same math, different physics.

🎮 Fluid simulation

Realistic water and smoke in games and films are made divergence-free at each frame by solving a Poisson equation — that projection step is the continuity equation imposed numerically.

🤖 ML & normalizing flows

Continuous normalizing flows model probability density as it transports through a velocity field. The same continuity equation governs how a probability cloud reshapes during sampling.

Common Pitfalls

Pitfall	What goes wrong	Fix
Forgetting the density factor	Writing $\nabla\!\cdot\mathbf{v} = 0$ for a compressible flow.	Always keep the full $\partial\rho/\partial t + \nabla\!\cdot(\rho\mathbf{v}) = 0$ unless ρ is genuinely constant.
Confusing divergence and curl	Treating a rotating flow as compressible because it ‘looks’ like it’s moving sideways.	Rotation gives non-zero curl but zero divergence. Solid-body rotation passes continuity.
Sign errors on the outward normal	Adding +u on the left face instead of −u — flux comes out double or zero.	Always evaluate v · n with n pointing outward. The minus on left/bottom faces is the whole point.
Mixing Eulerian and Lagrangian frames	Computing ‘rate of mass change’ inside a moving fluid parcel using ∂ρ/∂t.	For a moving region, use the material derivative Dρ/Dt and the Reynolds transport theorem.
Assuming a source term doesn't exist	In combustion, evaporation, or reactive flows, mass really is created or destroyed locally.	Add a source term S on the RHS: $\partial\rho/\partial t + \nabla\!\cdot(\rho\mathbf{v}) = S$ .

Summary

The continuity equation is the mathematical statement that mass is conserved at every point of a fluid. It takes the simplest possible accounting principle — "what flows in must either accumulate or flow out" — and translates it into a partial differential equation that holds at every point simultaneously.

Key Takeaways

General form: $\partial\rho/\partial t + \nabla\!\cdot(\rho\,\mathbf{v}) = 0$ — local mass conservation.
Incompressible form: $\nabla\!\cdot\mathbf{v} = 0$ — when density is constant, the velocity field is divergence-free.
The divergence theorem is the bridge between the bulk (surface-integral) statement and the pointwise (differential) statement.
Rotation ≠ compressibility: a fluid can swirl as violently as it likes and still satisfy continuity, as long as no mass appears or disappears.
In 2D, the stream function automatically generates divergence-free velocity fields — a trick used across physics and deep learning whenever an incompressibility constraint matters.

The Continuity Equation:

"What flows in must either pile up or flow out — at every point, at every instant."

Coming Next: Section 4 examines viscosity and the stress tensor — the physical origin of the $\mu\nabla^2\mathbf{v}$ term in the momentum equation, and the reason real fluids resist being deformed.