Boo-AI — Master Artificial Intelligence by Building from Scratch

Learning Objectives

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

Explain why the Navier-Stokes equations are just Newton's second law applied to a chunk of fluid you can shrink to a point.
Derive the continuity equation from conservation of mass using the Reynolds Transport Theorem.
Construct the material derivative $\frac{D}{Dt} = \frac{\partial}{\partial t} + (\mathbf{v}\cdot\nabla)$ and explain its two pieces physically.
Derive Cauchy's momentum equation from a force balance on an infinitesimal cube.
Apply the Newtonian constitutive law to close the system and obtain the full Navier-Stokes equations.
Implement the viscous diffusion term in both plain Python and PyTorch, and predict its output by hand.

The Big Picture: Newton's 2nd Law on a Tiny Box

Every equation in this section is one sentence: for any small chunk of fluid, the rate of change of momentum equals the net force on it. The work is just turning that one sentence into mathematics that survives shrinking the chunk to a single point.

Picture a swimming pool. Mentally draw a small cube anywhere inside it. The water in that cube has some momentum. A moment later that water has moved a little, swirled a little, and its momentum has changed. By how much? Newton's second law tells us exactly:

\text{mass} \times \text{acceleration} \;=\; \text{net force}

That equation is true for a baseball, a planet — and a cube of water. The problem is that "the cube of water" is being constantly reshuffled by neighbouring water flowing in and out. We need a clean way to talk about how a property of the fluid changes when the fluid is itself moving. That tool is the Reynolds Transport Theorem, and it is the trick that unlocks all of fluid dynamics.

Goal A

Conserve mass — get the continuity equation.

Goal B

Conserve momentum — get Cauchy's equation.

Goal C

Close Cauchy with a Newtonian stress law — get Navier-Stokes.

Two Views of a Fluid: Lagrangian vs Eulerian

Before writing a single derivative, we need to agree on what we are differentiating with respect to. There are two natural choices, and they look at the same fluid from opposite vantage points.

🚣 Lagrangian view — ride along with a particle

Pick one fluid particle. Follow it. Track its position $\mathbf{x}(t)$ , velocity $\mathbf{v}(t)$ , temperature $T(t)$ . This is what Newton's second law is literally about — "m a = F" for one particle.

Analogy: dye dropped into the pool; you watch the dye blob.

📷 Eulerian view — fix a window in space

Pick a fixed point $\mathbf{x}$ in space. Track what is passing through it: $\mathbf{v}(\mathbf{x},t)$ , $p(\mathbf{x},t)$ , $T(\mathbf{x},t)$ . This is what experiments and numerical grids actually measure.

Analogy: a thermometer bolted to the lane line; you watch whatever water happens to be touching it.

Why both views matter

Newton's law lives in the Lagrangian world. Measurements and computer grids live in the Eulerian world. The bridge between them is the material derivative, which we'll build in Step 3. The Reynolds Transport Theorem is the same bridge applied to volume integrals.

Step 1 — Reynolds Transport Theorem

Let $B(t) = \int_{V(t)} \phi(\mathbf{x},t)\, dV$ be the total amount of some property (mass, momentum, energy) inside a material volume $V(t)$ — meaning a blob of fluid that we follow as it moves and deforms. The Reynolds Transport Theorem translates the rate of change of $B(t)$ from the moving Lagrangian volume to a fixed Eulerian volume $V$ with the same instantaneous shape:

\displaystyle \frac{d}{dt}\int_{V(t)} \phi\, dV \;=\; \int_{V} \frac{\partial \phi}{\partial t}\, dV \;+\; \oint_{\partial V} \phi\,(\mathbf{v}\cdot\mathbf{n})\, dA

In words: total change inside the moving blob equals local change at every point of the (fixed-shape) box plus net flux of $\phi$ being carried across its boundary. The second term is a bookkeeping correction: stuff that flowed into our snapshot box wasn't actually "created" — it was just transported in by the velocity field.

LHS

\frac{d}{dt}\int_{V(t)} \phi\, dV

Total in the moving blob. Newton's view.

Unsteady term

\int_V \frac{\partial \phi}{\partial t}\, dV

How the field is changing point-by-point with the box held fixed.

Flux term

\oint_{\partial V} \phi\,(\mathbf{v}\cdot\mathbf{n})\, dA

Net amount of $\phi$ leaving through the box surface.

Where RTT comes from

It's a multivariable generalisation of the Leibniz rule for differentiating an integral with moving limits. The boundary moves with velocity $\mathbf{v}$ , so the surface integral picks up the familiar "evaluate at the moving endpoint" correction — just spread out over a closed surface instead of two endpoints.

Step 2 — Conservation of Mass

We pick $\phi = \rho$ in RTT. Conservation of mass says: the total mass inside a material volume cannot change (no nuclear reactions, please).

\displaystyle \frac{d}{dt}\int_{V(t)} \rho\, dV \;=\; 0

Apply RTT to the left-hand side and use the divergence theorem on the flux term (turning the surface integral into a volume integral):

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

Since this holds for every control volume — and only the zero function integrates to zero on every region — the integrand itself must vanish. We get the continuity equation:

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

For an incompressible fluid $\rho$ is constant, and the continuity equation collapses to the famously elegant

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

— "the velocity field is divergence-free." Visually: whatever enters a tiny box, the same amount leaves. The visualizer below makes this physical.

Interactive: Flux Through a Control Volume

Loading control-volume visualizer…

Try this: pick the Source field. Net outflow ≈ 2a · (box area). Now pick Vortex — net outflow is essentially zero no matter where you drag the box. That's ∇·v = 0 in action: the fluid is just spinning, not accumulating.

Step 3 — The Material Derivative

Newton's law wants the acceleration of a fluid particle — its Lagrangian acceleration. But the velocity field we work with is Eulerian: $\mathbf{v}(\mathbf{x},t)$ tells you what the velocity is at a fixed point. How are the two related?

Let $\mathbf{x}(t)$ be the trajectory of a particle. Its velocity is $\mathbf{v}(\mathbf{x}(t), t)$ . Take the time derivative by the chain rule:

\displaystyle \frac{d\mathbf{v}}{dt} \;=\; \frac{\partial \mathbf{v}}{\partial t} \;+\; \frac{\partial \mathbf{v}}{\partial x}\frac{dx}{dt} + \frac{\partial \mathbf{v}}{\partial y}\frac{dy}{dt} + \frac{\partial \mathbf{v}}{\partial z}\frac{dz}{dt}

But $(dx/dt, dy/dt, dz/dt) = \mathbf{v}$ . So

\displaystyle \boxed{\;\frac{D\mathbf{v}}{Dt} \;=\; \underbrace{\frac{\partial \mathbf{v}}{\partial t}}_{\text{local}} \;+\; \underbrace{(\mathbf{v}\!\cdot\!\nabla)\mathbf{v}}_{\text{convective}}\;}

This is the material derivative. It packages "the rate of change a particle actually experiences" into one symbol $D/Dt$ . The two terms answer two different questions:

Term	What it answers	Zero when…
∂v/∂t (local)	How is the field changing right here over time?	the flow is steady — the field doesn't depend on t
(v·∇)v (convective)	How does v differ from one point to the next, weighted by where the particle is heading?	the particle moves perpendicular to all gradients, or the field is uniform

Both terms can vanish independently

A garden hose with the spigot fully open produces a steady jet — the velocity at the nozzle doesn't change in time, so the local term is zero. But a fluid particle accelerates from rest in the pipe to high speed at the nozzle — the convective term $(\mathbf{v}\cdot\nabla)\mathbf{v}$ is doing all the work.

Interactive: Local vs Convective Rate of Change

We replace $\mathbf{v}$ with a scalar field $T(\mathbf{x})$ for clarity. The same formula applies:

\displaystyle \frac{DT}{Dt} \;=\; \frac{\partial T}{\partial t} \;+\; \mathbf{v}\!\cdot\!\nabla T

Loading material-derivative visualizer…

Step 4 — Momentum Balance (Cauchy's Equation)

Now we apply RTT with $\phi = \rho\mathbf{v}$ (momentum per unit volume) and equate the rate of change of momentum to the sum of forces. There are two kinds of forces on a fluid element:

Body forces per unit volume $\rho\mathbf{f}$ — gravity, electromagnetic, etc. They act on every molecule of the cube.
Surface forces per unit area — pressure pushing on each face, and viscous stresses tangential to each face. Both are captured in the stress tensor $\boldsymbol{\sigma}$ .

The momentum balance reads

\displaystyle \frac{d}{dt}\int_{V(t)} \rho \mathbf{v}\, dV \;=\; \int_V \rho \mathbf{f}\, dV \;+\; \oint_{\partial V} \boldsymbol{\sigma}\cdot\mathbf{n}\, dA

Apply RTT to the left, the divergence theorem to the surface integral on the right, and use the continuity equation to simplify. The volume integrals must agree for every blob, so the integrands match. After the dust settles you get Cauchy's equation of motion:

\displaystyle \rho\,\frac{D\mathbf{v}}{Dt} \;=\; \rho\mathbf{f} \;+\; \nabla\!\cdot\!\boldsymbol{\sigma}

Read it like Newton's law: "density × acceleration of the particle = body force + divergence of stress". We're done with calculus. Everything left is physics — what form does $\boldsymbol{\sigma}$ take?

Step 5 — The Stress Tensor

The stress tensor $\boldsymbol{\sigma}$ is a 3×3 matrix whose entry $\sigma_{ij}$ is the force per unit area in direction $i$ acting on a face whose outward normal points in direction $j$ . Diagonals push or pull perpendicular to the face (normal stress); off-diagonals shear it (tangential stress).

\boldsymbol{\sigma} \;=\; \begin{pmatrix} \sigma_{xx} & \sigma_{xy} & \sigma_{xz} \\ \sigma_{yx} & \sigma_{yy} & \sigma_{yz} \\ \sigma_{zx} & \sigma_{zy} & \sigma_{zz} \end{pmatrix}

For a fluid at rest, only pressure acts, and it pushes inward equally in every direction:

\boldsymbol{\sigma}_{\text{static}} \;=\; -p\,\mathbf{I}

For a fluid in motion we add a deviatoric piece $\boldsymbol{\tau}$ that captures how molecules drag past their neighbours:

\boldsymbol{\sigma} \;=\; -p\,\mathbf{I} \;+\; \boldsymbol{\tau}

Cauchy's equation now becomes:

\displaystyle \rho\frac{D\mathbf{v}}{Dt} \;=\; -\nabla p \;+\; \nabla\!\cdot\!\boldsymbol{\tau} \;+\; \rho\mathbf{f}

We have one equation and an unknown 3×3 tensor field $\boldsymbol{\tau}$ . We need a constitutive law — an extra physical assumption that tells us how $\boldsymbol{\tau}$ depends on the flow.

Step 6 — The Newtonian Closure

Newton's hypothesis (1687, restated rigorously by Stokes in 1845): viscous stress is linear in the rate of strain. For an incompressible Newtonian fluid:

\tau_{ij} \;=\; \mu\!\left( \frac{\partial v_i}{\partial x_j} + \frac{\partial v_j}{\partial x_i} \right)

The combination in parentheses is twice the symmetric part of the velocity gradient — the rate-of-strain tensor. The constant $\mu$ is the dynamic viscosity; honey has a big $\mu$ , air has a tiny one.

Taking the divergence (and using $\nabla\cdot\mathbf{v}=0$ for incompressibility) gives the marvellous simplification:

\nabla\!\cdot\!\boldsymbol{\tau} \;=\; \mu\,\nabla^2 \mathbf{v}

That is the famous Laplacian-of-velocity viscous term. Physically: viscosity tries to diffuse velocity — sharp gradients get smeared out, just like heat diffusing through a metal rod.

Why linear is reasonable

Linearity is exact when the deformation is slow compared to molecular relaxation times — true for water, air, oil under normal conditions. Non-Newtonian fluids (ketchup, blood, polymer melts) need a fancier $\boldsymbol{\tau}$ law, but the rest of the derivation is unchanged.

Step 7 — Assembling the Navier-Stokes Equations

Substitute the Newtonian closure into Cauchy's equation:

Incompressible Navier-Stokes

\displaystyle \rho\!\left( \frac{\partial \mathbf{v}}{\partial t} + (\mathbf{v}\!\cdot\!\nabla)\mathbf{v} \right) \;=\; -\nabla p \;+\; \mu\,\nabla^2 \mathbf{v} \;+\; \rho\mathbf{f}

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

Four scalar equations (three components of momentum + continuity) in four scalar unknowns ( $u, v, w, p$ ). Let's name every piece one last time so the equation never looks foreign again:

Piece	Physical meaning	Source
ρ ∂v/∂t	Mass × local acceleration (unsteady inertia)	Reynolds Transport Theorem
ρ (v·∇)v	Mass × convective acceleration — being swept through gradients	Material derivative chain rule
−∇p	Net pressure push (high → low)	Pressure part of σ
μ ∇²v	Viscous diffusion of momentum	Newtonian closure of τ
ρf	Body forces (gravity, etc.)	Newton's second law
∇·v = 0	Mass conservation for incompressible flow	Continuity equation

The nonlinearity that breaks everything

$(\mathbf{v}\!\cdot\!\nabla)\mathbf{v}$ is the only nonlinear term — and it is the source of turbulence, vortex shedding, weather chaos, and the open Millennium Prize problem. Without it the equations would collapse to a heat equation for velocity and we could solve them on paper.

Worked Example: 1D Viscous Diffusion of a Velocity Bump

We strip the full equation to its bones to isolate the viscous term. Take $p$ , $\mathbf{f}$ , and the convective term all to zero. Only one component $u(y,t)$ of velocity survives, depending on one coordinate $y$ between two parallel walls. Navier-Stokes collapses to:

\displaystyle \frac{\partial u}{\partial t} \;=\; \nu\,\frac{\partial^2 u}{\partial y^2}, \qquad u(0,t)=u(L,t)=0

— pure 1D diffusion of momentum, with no-slip walls at $y=0,L$ . This is the same equation as heat in a rod. Below we work it out by hand on a tiny grid, then verify with code.

📐 Click to expand: hand calculation on a 5-node grid

Use 5 grid points $y_0,\dots,y_4$ with spacing $\Delta y = 1$ , and a time step $\Delta t = 0.1$ , kinematic viscosity $\nu = 1$ . The discrete second derivative at interior node $j$ is:

\displaystyle \left.\frac{\partial^2 u}{\partial y^2}\right|_j \;\approx\; \frac{u_{j+1} - 2u_j + u_{j-1}}{\Delta y^{2}}

The explicit Euler update is then

u_j^{n+1} \;=\; u_j^n + r\,(u_{j+1}^n - 2 u_j^n + u_{j-1}^n), \quad r = \frac{\nu \Delta t}{\Delta y^{2}}

For our numbers $r = 1 \cdot 0.1 / 1^2 = 0.1$ . Initial profile: a unit spike at the middle node $u^0 = (0,\,0,\,1,\,0,\,0)$ . Walls pinned to zero. One step:

j	u^0[j]	u^0[j+1] − 2u^0[j] + u^0[j−1]	u^1[j] = u^0[j] + 0.1 × …
0	0	— (wall, pinned to 0)	0
1	0	1 − 0 + 0 = 1	0 + 0.1 × 1 = 0.1
2	1	0 − 2 + 0 = −2	1 + 0.1 × (−2) = 0.8
3	0	0 − 0 + 1 = 1	0 + 0.1 × 1 = 0.1
4	0	— (wall, pinned to 0)	0

After one step: $u^1 = (0,\,0.1,\,0.8,\,0.1,\,0)$ . The spike lost $2r = 0.2$ from the middle and donated $r = 0.1$ to each neighbour. That redistribution is literally what $\nu\,\nabla^2 u$ means in this problem.

Total "mass of velocity" $\sum_j u_j$ is conserved at the interior level (0.1 + 0.8 + 0.1 = 1.0), but the walls drain it. As $n\to\infty$ the profile flattens to zero everywhere — momentum has fully dissipated through the walls.

Interactive: Viscosity Smears Momentum

Loading viscous-diffusion visualizer…

Python Implementation

The whole worked example is below as 20 lines of NumPy. The output of the toy run matches our table by hand exactly:

1D viscous diffusion — plain Python / NumPy

🐍viscous_diffusion.py

Explanation(15)

Code(29)

1Import NumPy

We use NumPy for array slicing. The vectorised expression u[2:] - 2*u[1:-1] + u[:-2] computes the second difference at every interior point in one shot — without writing a Python loop over j.

EXAMPLE

import numpy as np

3Function signature — what the inputs mean

u0 is the initial velocity profile sampled at every grid node along y. nu is the kinematic viscosity ν = μ/ρ. dt is the time step. dy is the spacing between grid nodes. n_steps is how many time steps to advance. We return the final state and a history list so you can plot the evolution.

EXAMPLE

viscous_diffusion_1d(u0=[0,0,1,0,0], nu=1.0, dt=0.1, dy=1.0, n_steps=1)

10Copy so we don't mutate the caller's array

u0.copy() makes an independent array. If you wrote u = u0 instead, every change inside the loop would also change the array the caller passed in — a classic NumPy footgun.

EXAMPLE

a = np.array([1,2,3]); b = a; b[0] = 9   # a is now [9,2,3] — bad

11Diffusion number r = ν dt / dy²

This is the dimensionless ratio that controls how aggressively each step diffuses. It appears because the central second difference has units 1/dy², the time step contributes dt, and ν provides the physical scale. For our example: r = 1.0 × 0.1 / 1² = 0.1.

EXAMPLE

r = 1.0 * 0.1 / 1.0**2  →  0.1

12Stability check — von Neumann criterion

Explicit FTCS for the heat/viscous equation is only stable when r ≤ 1/2. Cross that line and tiny round-off errors get amplified every step until the numerical solution explodes. The assert turns that mathematical fact into a runtime guard.

EXAMPLE

r = 0.6  → assertion fails because  u[j+1] - 2 u[j] + u[j-1]  amplifies, not damps

14Record history

We snapshot the profile at every step so we can plot it, validate against an analytical solution, or animate. Storing copies (not references) is essential — otherwise every slot in history would point at the same array that keeps changing.

15Time loop — march forward in time

Each iteration advances the solution by one dt. We don't care about the loop variable, only that it runs n_steps times, so we use _ as the throwaway name.

EXAMPLE

for _ in range(1):  # only one step in our toy example

16Buffer the new state

We compute u_new from u — never overwriting u while we still need its old values. If we wrote u[1:-1] = ... directly, the j=2 update would already see the freshly-computed u[1], silently using mixed time levels.

18Vectorised second difference (the viscous term)

This single line is the heart of ∂²u/∂y² ≈ (u[j+1] − 2 u[j] + u[j−1])/dy². For our example with u = [0, 0, 1, 0, 0]: at j=1 we get (1 − 0 + 0)/1 = 1, so u_new[1] = 0 + 0.1·1 = 0.1. At j=2: (0 − 2 + 0)/1 = −2, so u_new[2] = 1 + 0.1·(−2) = 0.8. At j=3: (0 − 0 + 1)/1 = 1, so u_new[3] = 0 + 0.1·1 = 0.1.

EXAMPLE

j=2: r*(u[3] - 2*u[2] + u[1]) = 0.1*(0 - 2*1 + 0) = -0.2  →  u_new[2] = 1 + (-0.2) = 0.8

19No-slip wall: u(y=0) = 0

Real fluid molecules in contact with a stationary wall don't slip past it. So the velocity at the wall is pinned to zero. We hard-set it after the interior update, overriding whatever the stencil produced. This is a Dirichlet boundary condition.

EXAMPLE

Imagine syrup against a glass plate — the layer touching the glass is motionless.

20No-slip wall: u(y=L) = 0

Same condition on the other wall. Together these two pins are what makes momentum eventually decay to zero — the walls keep sucking energy out of the column.

21Roll forward: u ← u_new

We swap names. Next iteration u refers to the just-computed state. This is the time-marching pattern every explicit PDE solver uses.

23Return final state + history

u is the final profile; history is every intermediate snapshot. Pair (final, trajectory) is the standard interface for simulation drivers — give the caller the answer they probably want plus everything they might want to inspect.

26Initial condition — a spike at the middle node

Five grid points along y. All velocity concentrated at j=2 (the middle). This isolates the question: 'what does ν ∇²u do to a sharp feature?' The answer (next line) is that it immediately spreads.

27Run one step and print

We march exactly one step so the values are easy to verify by hand. r = 0.1, the spike of height 1 loses 2r = 0.2 and gives r = 0.1 to each neighbour. After one step the spike has flattened from [0,0,1,0,0] to [0, 0.1, 0.8, 0.1, 0]. Confirmed by the printed output.

EXAMPLE

after one step: [0.   0.1  0.8  0.1  0. ]

14 lines without explanation

1import numpy as np
2
3def viscous_diffusion_1d(u0, nu, dt, dy, n_steps):
4    """Solve  ∂u/∂t = ν ∂²u/∂y²  on a 1D grid with no-slip walls.
5
6    This is the *viscous term* of Navier-Stokes in isolation —
7    no pressure, no advection, no body force. Just internal friction
8    redistributing momentum.
9    """
10    u = u0.copy()
11    r = nu * dt / (dy ** 2)
12    assert r <= 0.5, "Explicit scheme unstable: pick smaller dt or larger dy"
13
14    history = [u.copy()]
15    for _ in range(n_steps):
16        u_new = u.copy()
17        # Central second difference: (u[j+1] - 2 u[j] + u[j-1]) / dy²
18        u_new[1:-1] = u[1:-1] + r * (u[2:] - 2 * u[1:-1] + u[:-2])
19        u_new[0]  = 0.0      # no-slip at bottom wall
20        u_new[-1] = 0.0      # no-slip at top wall
21        u = u_new
22        history.append(u.copy())
23    return u, history
24
25# ── Run a tiny example you can verify by hand ───────────────────────────────
26u0 = np.array([0.0, 0.0, 1.0, 0.0, 0.0])   # a momentum spike at the middle node
27u_final, hist = viscous_diffusion_1d(u0, nu=1.0, dt=0.1, dy=1.0, n_steps=1)
28print("after one step:", u_final)
29# → [0.   0.1  0.8  0.1  0. ]

What you just simulated: the right-hand side of NS with everything but the viscous term zeroed out. Add back the pressure gradient and you get the classic Stokes problem; add back the convective term and you get the full nonlinearity — and the open Millennium Prize question of whether it can ever blow up in 3D.

PyTorch Tensor Implementation

Same physics, same arithmetic, swap NumPy for PyTorch. For 1D it's overkill — but the moment you move to 2D / 3D Navier-Stokes on a real grid the GPU speedup is enormous, and you also get autograd for free (handy for inverse-problem and ML-flavoured CFD).

1D viscous diffusion — PyTorch tensor version

🐍viscous_diffusion_torch.py

Explanation(14)

Code(29)

1Import torch

torch.Tensor looks and acts like a NumPy array but lives in a compute graph and can be moved to GPU with .cuda(). For this 1D toy it makes no speed difference; for 2D NS on a 512×512 grid it makes a 50–200× difference.

3Type-annotated signature

Same arguments as the NumPy version — but u0 is now declared torch.Tensor. That single annotation tells future-you and your IDE that this function expects PyTorch, will use torch operations, and can participate in autograd or GPU dispatch.

EXAMPLE

u0 = torch.zeros(101)

10u.clone() — the torch equivalent of np.copy

Tensors share storage by default the same way NumPy arrays do. .clone() makes an independent copy. Always clone before mutating an input tensor inside a function — otherwise you'll silently corrupt the caller's data and (worse) confuse autograd.

EXAMPLE

a = torch.tensor([1.,2.,3.]);  b = a;  b[0] = 9   # a is now [9,2,3] too — bad

11Diffusion number r

Identical to the NumPy version. r is a plain Python float, computed once outside the loop. We don't make it a tensor — it would add unnecessary autograd bookkeeping.

EXAMPLE

r = 1.0 * 0.1 / 1.0**2  →  0.1

12Same stability assert

Stability of the FTCS scheme is a property of the discretisation, not the backend. r ≤ 1/2 still applies whether you compute in NumPy, PyTorch, or hand-calculate on graph paper.

14Snapshot history

u.clone() again, for the same reason as in NumPy. If we just appended u, every slot in history would point at the same tensor, and the final list would show only the last state repeated n+1 times.

15Time loop

Plain Python loop. Each step is a vectorised tensor update — the for-loop overhead is one iteration per dt, which is fine. The expensive part is inside the brackets, not the loop itself.

17u_new = u.clone() — buffer

We need u_new to start equal to u so that the walls (positions 0 and -1) keep their boundary values, and so any nodes we don't touch (none here, but in 2D there are corners) don't get random uninitialised data.

19Tensor slicing identical to NumPy

u[2:] - 2*u[1:-1] + u[:-2] is the second-difference stencil. PyTorch broadcasts the scalar r the same way NumPy does. For u = [0, 0, 1, 0, 0]: u[2:]=[1,0,0], u[1:-1]=[0,1,0], u[:-2]=[0,0,1]. The stencil gives [1-0+0, 0-2+0, 0-0+1] = [1, -2, 1]. Multiply by r=0.1 and add to u[1:-1]: [0+0.1, 1-0.2, 0+0.1] = [0.1, 0.8, 0.1].

EXAMPLE

stencil at interior:  [1, -2, 1] × 0.1  →  [0.1, -0.2, 0.1]

20No-slip wall at bottom

torch.Tensor assignment with a scalar broadcasts it. u_new[0] = 0.0 pins the velocity to zero at the first grid node. Note that 0.0 (Python float) is implicitly cast to the tensor's dtype.

21No-slip wall at top

Same for the last node. Together these two assignments encode the Dirichlet condition u(y=0,t) = u(y=L,t) = 0.

22Roll forward in time

u now refers to the freshly-computed state. PyTorch handles the underlying storage; the old tensor will be garbage-collected once nothing references it.

28Initial spike (torch.tensor)

torch.tensor(...) builds a 1D tensor on CPU with dtype torch.float32 by default. The values mirror the NumPy version exactly, so we can sanity-check that the two implementations agree.

EXAMPLE

u0.shape  → torch.Size([5])

29Run and print

The printed tensor matches NumPy's output to the bit — the two stencils computed the same arithmetic, just routed through different libraries. That's the parity check.

EXAMPLE

tensor([0.0000, 0.1000, 0.8000, 0.1000, 0.0000])

15 lines without explanation

1import torch
2
3def viscous_diffusion_1d_torch(u0: torch.Tensor, nu: float, dt: float, dy: float, n_steps: int):
4    """Same physics, tensor-flavoured.
5
6    Why bother? Because the moment we move to 2D / 3D Navier-Stokes,
7    torch.Tensor + GPU + autograd become the difference between
8    a 4-hour simulation and a 4-second one.
9    """
10    u = u0.clone()
11    r = nu * dt / (dy ** 2)
12    assert r <= 0.5
13
14    history = [u.clone()]
15    for _ in range(n_steps):
16        u_new = u.clone()
17        # Same stencil — but torch.Tensor slicing instead of NumPy
18        u_new[1:-1] = u[1:-1] + r * (u[2:] - 2 * u[1:-1] + u[:-2])
19        u_new[0]  = 0.0
20        u_new[-1] = 0.0
21        u = u_new
22        history.append(u.clone())
23    return u, history
24
25# ── Tiny example, mirror of the NumPy version ───────────────────────────────
26u0 = torch.tensor([0.0, 0.0, 1.0, 0.0, 0.0])
27u_final, hist = viscous_diffusion_1d_torch(u0, nu=1.0, dt=0.1, dy=1.0, n_steps=1)
28print("after one step:", u_final)
29# → tensor([0.0000, 0.1000, 0.8000, 0.1000, 0.0000])

Parity check

Both the NumPy and PyTorch versions print [0.0, 0.1, 0.8, 0.1, 0.0] after one step. That matches our hand calculation in the collapsible example above. Whenever you port a numerical routine across libraries, this kind of bit-for-bit equivalence is the test you want to keep around — it catches indexing bugs, off-by-one stencil shifts, and dtype mismatches.

Common Mistakes to Avoid

Confusing d/dt with ∂/∂t.

$\partial/\partial t$ holds $\mathbf{x}$ fixed; $d/dt$ (and the material derivative $D/Dt$ ) follows the particle. A steady flow has $\partial \mathbf{v}/\partial t = 0$ but the particle can still accelerate.

Forgetting that (v·∇)v is nonlinear.

It looks like a product of derivatives, but $\mathbf{v}$ itself is the unknown — that's what makes NS unsolved analytically in 3D. Treating it as a small perturbation is the Stokes-flow approximation, valid only for very low Reynolds numbers.

Dropping no-slip boundary conditions.

Without no-slip the viscous term has nothing to push against and energy is never dissipated. Every numerical NS solver pins the wall velocity to the wall's velocity — usually zero for stationary walls.

Choosing dt too large in explicit schemes.

For the FTCS viscous scheme, $r = \nu\Delta t/\Delta y^2 \le 1/2$ . Violate it and round-off blows up exponentially. Implicit schemes (Crank-Nicolson, backward Euler) have no such restriction but cost a linear solve per step.

Assuming μ ∇²v is the "friction" force.

It looks like friction, but it's really diffusion of momentum. Friction usually dissipates energy locally; viscosity transports momentum from fast layers to slow layers, then ultimately dumps it into the walls. The first sentence of every fluids paper that mishandles this is wrong.

Summary

We turned Newton's second law into the Navier-Stokes equations in seven steps:

Picked an arbitrary blob of fluid and decided to use Reynolds Transport Theorem as the bridge between Lagrangian (follow-the-particle) and Eulerian (fixed grid) descriptions.
Set $\phi = \rho$ and got the continuity equation $\partial_t \rho + \nabla\!\cdot\!(\rho\mathbf{v}) = 0$ , which becomes $\nabla\!\cdot\!\mathbf{v}=0$ for incompressible flow.
Built the material derivative $D/Dt = \partial_t + (\mathbf{v}\!\cdot\!\nabla)$ from the chain rule — local rate plus convective rate.
Set $\phi = \rho\mathbf{v}$ and balanced momentum against body + surface forces to get Cauchy's equation $\rho\,D\mathbf{v}/Dt = \rho\mathbf{f} + \nabla\!\cdot\!\boldsymbol{\sigma}$ .
Decomposed the stress tensor into pressure + deviatoric: $\boldsymbol{\sigma} = -p\mathbf{I} + \boldsymbol{\tau}$ .
Closed the system with the Newtonian linear law $\tau_{ij} = \mu(\partial_j v_i + \partial_i v_j)$ , which for incompressible flow gives $\nabla\!\cdot\!\boldsymbol{\tau} = \mu\nabla^2\mathbf{v}$ .
Assembled everything into the Navier-Stokes equations and verified the viscous term numerically in both NumPy and PyTorch.

The Navier-Stokes Equations:

"Newton's second law for a fluid: density × material derivative of velocity = pressure push + viscous diffusion + body forces."

Coming Next: we'll spend a whole section on the continuity equation alone — what it implies geometrically, how compressibility breaks it, and why divergence-free vector fields are central to fluid simulation algorithms like projection methods.