diff --git a/12-steps-to-Navier-Stokes_filled.ipynb b/12-steps-to-Navier-Stokes_filled.ipynb
new file mode 100644
index 0000000..f9b1faa
--- /dev/null
+++ b/12-steps-to-Navier-Stokes_filled.ipynb
@@ -0,0 +1,260945 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Navier–Stokes in 12 Easy Steps (Computational Workflow)\n",
+ "\n",
+ "This notebook is a **hands-on CFD mini-course** inspired by the classic *CFDPython* “12 steps” series.\n",
+ "It walks from simple 1D advection all the way to **2D incompressible Navier–Stokes** (cavity flow) and **channel flow**.\n",
+ "\n",
+ "## What you’ll learn\n",
+ "- How to discretize PDEs with **finite differences**\n",
+ "- How to reason about **stability (CFL)** and **numerical diffusion**\n",
+ "- How the **pressure Poisson equation** enforces incompressibility\n",
+ "- A reusable workflow: *model → discretize → implement → verify/visualize*\n",
+ "\n",
+ "> **Conventions**\n",
+ "> - Spatial grid: uniform, Cartesian\n",
+ "> - Time stepping: explicit (Forward Euler) unless noted\n",
+ "> - Boundary conditions: stated per step\n",
+ "> - Arrays are indexed as `u[j, i]` with **`j` = y-index** and **`i` = x-index**\n",
+ "\n",
+ "---\n",
+ "\n",
+ "### Setup\n",
+ "Run the next cell to import dependencies and define helper plotting functions used throughout.\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "id": "ef724e97",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Requirement already satisfied: numpy in /Users/asheshghosh/.julia/conda/3/aarch64/lib/python3.10/site-packages (1.26.4)\n",
+ "Collecting numpy\n",
+ " Downloading numpy-2.2.6-cp310-cp310-macosx_14_0_arm64.whl.metadata (62 kB)\n",
+ "\u001b[2K \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m62.0/62.0 kB\u001b[0m \u001b[31m1.7 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
+ "\u001b[?25hCollecting matplotlib\n",
+ " Downloading matplotlib-3.10.8-cp310-cp310-macosx_11_0_arm64.whl.metadata (52 kB)\n",
+ "\u001b[2K \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m52.8/52.8 kB\u001b[0m \u001b[31m2.7 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
+ "\u001b[?25hCollecting contourpy>=1.0.1 (from matplotlib)\n",
+ " Downloading contourpy-1.3.2-cp310-cp310-macosx_11_0_arm64.whl.metadata (5.5 kB)\n",
+ "Collecting cycler>=0.10 (from matplotlib)\n",
+ " Using cached cycler-0.12.1-py3-none-any.whl.metadata (3.8 kB)\n",
+ "Collecting fonttools>=4.22.0 (from matplotlib)\n",
+ " Downloading fonttools-4.61.1-cp310-cp310-macosx_10_9_universal2.whl.metadata (114 kB)\n",
+ "\u001b[2K \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m114.2/114.2 kB\u001b[0m \u001b[31m7.1 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
+ "\u001b[?25hCollecting kiwisolver>=1.3.1 (from matplotlib)\n",
+ " Downloading kiwisolver-1.4.9-cp310-cp310-macosx_11_0_arm64.whl.metadata (6.3 kB)\n",
+ "Requirement already satisfied: packaging>=20.0 in /Users/asheshghosh/.julia/conda/3/aarch64/lib/python3.10/site-packages (from matplotlib) (23.2)\n",
+ "Collecting pillow>=8 (from matplotlib)\n",
+ " Downloading pillow-12.1.1-cp310-cp310-macosx_11_0_arm64.whl.metadata (8.8 kB)\n",
+ "Collecting pyparsing>=3 (from matplotlib)\n",
+ " Downloading pyparsing-3.3.2-py3-none-any.whl.metadata (5.8 kB)\n",
+ "Requirement already satisfied: python-dateutil>=2.7 in /Users/asheshghosh/.julia/conda/3/aarch64/lib/python3.10/site-packages (from matplotlib) (2.9.0.post0)\n",
+ "Requirement already satisfied: six>=1.5 in /Users/asheshghosh/.julia/conda/3/aarch64/lib/python3.10/site-packages (from python-dateutil>=2.7->matplotlib) (1.17.0)\n",
+ "Downloading numpy-2.2.6-cp310-cp310-macosx_14_0_arm64.whl (5.3 MB)\n",
+ "\u001b[2K \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m5.3/5.3 MB\u001b[0m \u001b[31m46.8 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m00:01\u001b[0m00:01\u001b[0m\n",
+ "\u001b[?25hDownloading matplotlib-3.10.8-cp310-cp310-macosx_11_0_arm64.whl (8.1 MB)\n",
+ "\u001b[2K \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m8.1/8.1 MB\u001b[0m \u001b[31m45.9 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0ma \u001b[36m0:00:01\u001b[0m\n",
+ "\u001b[?25hDownloading contourpy-1.3.2-cp310-cp310-macosx_11_0_arm64.whl (253 kB)\n",
+ "\u001b[2K \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m253.4/253.4 kB\u001b[0m \u001b[31m31.7 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
+ "\u001b[?25hUsing cached cycler-0.12.1-py3-none-any.whl (8.3 kB)\n",
+ "Downloading fonttools-4.61.1-cp310-cp310-macosx_10_9_universal2.whl (2.9 MB)\n",
+ "\u001b[2K \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m2.9/2.9 MB\u001b[0m \u001b[31m49.9 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0ma \u001b[36m0:00:01\u001b[0m\n",
+ "\u001b[?25hDownloading kiwisolver-1.4.9-cp310-cp310-macosx_11_0_arm64.whl (65 kB)\n",
+ "\u001b[2K \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m65.3/65.3 kB\u001b[0m \u001b[31m7.9 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
+ "\u001b[?25hDownloading pillow-12.1.1-cp310-cp310-macosx_11_0_arm64.whl (4.7 MB)\n",
+ "\u001b[2K \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m4.7/4.7 MB\u001b[0m \u001b[31m51.0 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0ma \u001b[36m0:00:01\u001b[0m\n",
+ "\u001b[?25hDownloading pyparsing-3.3.2-py3-none-any.whl (122 kB)\n",
+ "\u001b[2K \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m122.8/122.8 kB\u001b[0m \u001b[31m18.1 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
+ "\u001b[?25hInstalling collected packages: pyparsing, pillow, numpy, kiwisolver, fonttools, cycler, contourpy, matplotlib\n",
+ " Attempting uninstall: numpy\n",
+ " Found existing installation: numpy 1.26.4\n",
+ " Uninstalling numpy-1.26.4:\n",
+ " Successfully uninstalled numpy-1.26.4\n",
+ "Successfully installed contourpy-1.3.2 cycler-0.12.1 fonttools-4.61.1 kiwisolver-1.4.9 matplotlib-3.10.8 numpy-2.2.6 pillow-12.1.1 pyparsing-3.3.2\n"
+ ]
+ }
+ ],
+ "source": [
+ "import sys\n",
+ "!{sys.executable} -m pip install -U numpy matplotlib"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 15,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Core scientific stack\n",
+ "import numpy as np\n",
+ "import matplotlib.pyplot as plt\n",
+ "from matplotlib import cm\n",
+ "\n",
+ "# Nicer inline figures (works in Jupyter/Quarto)\n",
+ "# If you run this as a .py script, you can remove the next line.\n",
+ "%config InlineBackend.figure_formats = ['svg']\n",
+ "\n",
+ "plt.rcParams.update({\n",
+ " \"figure.figsize\": (6, 5),\n",
+ " \"figure.dpi\": 300,\n",
+ " \"axes.grid\": True,\n",
+ " \"axes.spines.top\": False,\n",
+ " \"axes.spines.right\": False,\n",
+ "})\n",
+ "\n",
+ "def plot1d(x, u, title=\"\", xlabel=\"x\", ylabel=\"u\"):\n",
+ " \"\"\"Quick 1D line plot.\"\"\"\n",
+ " fig, ax = plt.subplots()\n",
+ " ax.plot(x, u, lw=2)\n",
+ " ax.set_title(title)\n",
+ " ax.set_xlabel(xlabel)\n",
+ " ax.set_ylabel(ylabel)\n",
+ " fig.tight_layout()\n",
+ " plt.show()\n",
+ "\n",
+ "def plot2d_contour(X, Y, Z, title=\"\", xlabel=\"x\", ylabel=\"y\", cmap=cm.viridis):\n",
+ " \"\"\"2D filled contour plot with colorbar.\"\"\"\n",
+ " fig, ax = plt.subplots()\n",
+ " c = ax.contourf(X, Y, Z, levels=50, cmap=cmap)\n",
+ " fig.colorbar(c, ax=ax, shrink=0.85)\n",
+ " ax.set_title(title)\n",
+ " ax.set_xlabel(xlabel)\n",
+ " ax.set_ylabel(ylabel)\n",
+ " ax.set_aspect(\"equal\", adjustable=\"box\")\n",
+ " fig.tight_layout()\n",
+ " plt.show()\n",
+ "\n",
+ "def plot2d_surface(X, Y, Z, title=\"\", xlabel=\"x\", ylabel=\"y\", zlabel=\"\"):\n",
+ " \"\"\"3D surface plot (robust across Matplotlib versions).\"\"\"\n",
+ " fig = plt.figure(figsize=(11, 7), dpi=110)\n",
+ " ax = fig.add_subplot(111, projection=\"3d\") # reliable vs fig.gca(projection='3d')\n",
+ " surf = ax.plot_surface(X, Y, Z, cmap=cm.viridis, rstride=1, cstride=1,\n",
+ " linewidth=0, antialiased=True)\n",
+ " fig.colorbar(surf, ax=ax, shrink=0.6, pad=0.08)\n",
+ " ax.set_title(title)\n",
+ " ax.set_xlabel(xlabel)\n",
+ " ax.set_ylabel(ylabel)\n",
+ " if zlabel:\n",
+ " ax.set_zlabel(zlabel)\n",
+ " fig.tight_layout()\n",
+ " plt.show()\n",
+ "\n",
+ "def cfl_number(c, dt, dx):\n",
+ " \"\"\"CFL number for 1D advection-like terms.\"\"\"\n",
+ " return c * dt / dx\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Step 1 — 1D Linear Convection\n",
+ "\n",
+ "We start with the simplest hyperbolic PDE:\n",
+ "\n",
+ "\n",
+ "$$\\frac{\\partial u}{\\partial t} + c\\,\\frac{\\partial u}{\\partial x} = 0$$\n",
+ "\n",
+ "where **`c`** is a constant wave speed. This equation transports the profile to the right (for `c>0`) without changing its shape (in the continuous limit).\n",
+ "\n",
+ "### Discretization (Forward Euler + backward difference)\n",
+ "For stability with `c>0`, we use an **upwind** (backward) difference:\n",
+ "\n",
+ "$$\n",
+ "u_i^{n+1} = u_i^n - c\\,\\frac{\\Delta t}{\\Delta x}\\left(u_i^n-u_{i-1}^n\\right).\n",
+ "$$\n",
+ "\n",
+ "### Boundary conditions\n",
+ "We hold the left boundary fixed (`u[0]=1`), and the right boundary is also kept at `1` for this demo.\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 16,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "CFL = 0.20000000000000004\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "# --- Parameters ---\n",
+ "nx = 81\n",
+ "nt = 100\n",
+ "c = 1.0\n",
+ "dx = 2.0 / (nx - 1)\n",
+ "sigma = 0.2 # CFL-like safety factor\n",
+ "dt = sigma * dx / c # choose dt from dx and c\n",
+ "\n",
+ "print(\"CFL =\", cfl_number(c, dt, dx))\n",
+ "\n",
+ "# --- Grid ---\n",
+ "x = np.linspace(0, 2, nx)\n",
+ "\n",
+ "# --- Initial condition: \"hat\" function ---\n",
+ "u = np.ones(nx)\n",
+ "u[int(0.5 / dx):int(1.0 / dx) + 1] = 2.0\n",
+ "\n",
+ "plot1d(x, u, title=\"Step 1: Initial condition (1D linear convection)\")\n",
+ "\n",
+ "# --- Time integration ---\n",
+ "for n in range(nt):\n",
+ " un = u.copy()\n",
+ " # Update interior points; i=0 is boundary\n",
+ " u[1:] = un[1:] - c * dt / dx * (un[1:] - un[:-1])\n",
+ "\n",
+ " # Boundary conditions\n",
+ " u[0] = 1.0\n",
+ " u[-1] = 1.0\n",
+ "\n",
+ "plot1d(x, u, title=f\"Step 1: After nt={nt} steps\")\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Step 2 — 1D Nonlinear Convection\n",
+ "\n",
+ "Now the wave speed depends on the solution:\n",
+ "\n",
+ "$$\n",
+ "\\frac{\\partial u}{\\partial t} + u\\,\\frac{\\partial u}{\\partial x} = 0.\n",
+ "$$\n",
+ "\n",
+ "This creates **steepening**: higher values move faster, potentially forming shocks.\n",
+ "\n",
+ "### Discretization (upwind)\n",
+ "$$\n",
+ "u_i^{n+1} = u_i^n - u_i^n\\,\\frac{\\Delta t}{\\Delta x}\\left(u_i^n-u_{i-1}^n\\right).\n",
+ "$$\n",
+ "\n",
+ "We again use a backward difference (upwind) for stability.\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 17,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "nx = 101\n",
+ "nt = 200\n",
+ "dx = 2.0 / (nx - 1)\n",
+ "dt = 0.2 * dx # simple choice; for robustness you can use dt ~ dx / max(u)\n",
+ "\n",
+ "x = np.linspace(0, 2, nx)\n",
+ "u = np.ones(nx)\n",
+ "u[int(0.5 / dx):int(1.0 / dx) + 1] = 2.0\n",
+ "\n",
+ "plot1d(x, u, title=\"Step 2: Initial condition (1D nonlinear convection)\")\n",
+ "\n",
+ "for n in range(nt):\n",
+ " un = u.copy()\n",
+ " u[1:] = un[1:] - un[1:] * dt / dx * (un[1:] - un[:-1])\n",
+ " u[0] = 1.0\n",
+ " u[-1] = 1.0\n",
+ "\n",
+ "plot1d(x, u, title=f\"Step 2: After nt={nt} steps\")\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Step 3 — 1D Diffusion\n",
+ "\n",
+ "Diffusion smooths gradients:\n",
+ "\n",
+ "$$\n",
+ "\\frac{\\partial u}{\\partial t} = \\nu\\,\\frac{\\partial^2 u}{\\partial x^2},\n",
+ "$$\n",
+ "where **`ν`** is the diffusion coefficient (viscosity in fluid contexts).\n",
+ "\n",
+ "### Discretization (Forward Euler + centered second derivative)\n",
+ "$$\n",
+ "u_i^{n+1} = u_i^n + \\nu\\,\\frac{\\Delta t}{\\Delta x^2}\\left(u_{i+1}^n - 2u_i^n + u_{i-1}^n\\right).\n",
+ "$$\n",
+ "\n",
+ "### Stability (important!)\n",
+ "For explicit diffusion in 1D, a typical condition is:\n",
+ "\n",
+ "$$\n",
+ "\\nu\\,\\frac{\\Delta t}{\\Delta x^2} \\le \\frac{1}{2}.\n",
+ "$$\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 18,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "nx = 101\n",
+ "nt = 200\n",
+ "nu = 0.3\n",
+ "dx = 2.0 / (nx - 1)\n",
+ "sigma = 0.2\n",
+ "dt = sigma * dx**2 / nu # stable explicit diffusion step\n",
+ "\n",
+ "x = np.linspace(0, 2, nx)\n",
+ "u = np.ones(nx)\n",
+ "u[int(0.5 / dx):int(1.0 / dx) + 1] = 2.0\n",
+ "\n",
+ "plot1d(x, u, title=\"Step 3: Initial condition (1D diffusion)\")\n",
+ "\n",
+ "for n in range(nt):\n",
+ " un = u.copy()\n",
+ " u[1:-1] = un[1:-1] + nu * dt / dx**2 * (un[2:] - 2*un[1:-1] + un[:-2])\n",
+ " u[0] = 1.0\n",
+ " u[-1] = 1.0\n",
+ "\n",
+ "plot1d(x, u, title=f\"Step 3: After nt={nt} steps\")\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Step 4 — 1D Burgers’ Equation (nonlinear convection + diffusion)\n",
+ "\n",
+ "Burgers combines the previous two effects:\n",
+ "\n",
+ "$$\n",
+ "\\frac{\\partial u}{\\partial t} + u\\,\\frac{\\partial u}{\\partial x} = \\nu\\,\\frac{\\partial^2 u}{\\partial x^2}.\n",
+ "$$\n",
+ "\n",
+ "It is a standard playground for shock formation (from nonlinearity) regularized by viscosity (diffusion).\n",
+ "\n",
+ "### Discretization\n",
+ "- Nonlinear convection term: upwind\n",
+ "- Diffusion term: centered second derivative\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 19,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "/var/folders/l_/_p2_9s4x6sgbl7sh8tjclkyc0000gn/T/ipykernel_15673/1776893760.py:16: RuntimeWarning: overflow encountered in multiply\n",
+ " - un[1:-1] * dt/dx * (un[1:-1] - un[:-2])\n",
+ "/var/folders/l_/_p2_9s4x6sgbl7sh8tjclkyc0000gn/T/ipykernel_15673/1776893760.py:16: RuntimeWarning: invalid value encountered in subtract\n",
+ " - un[1:-1] * dt/dx * (un[1:-1] - un[:-2])\n",
+ "/var/folders/l_/_p2_9s4x6sgbl7sh8tjclkyc0000gn/T/ipykernel_15673/1776893760.py:15: RuntimeWarning: invalid value encountered in subtract\n",
+ " u[1:-1] = (un[1:-1]\n",
+ "/var/folders/l_/_p2_9s4x6sgbl7sh8tjclkyc0000gn/T/ipykernel_15673/1776893760.py:17: RuntimeWarning: invalid value encountered in subtract\n",
+ " + nu * dt/dx**2 * (un[2:] - 2*un[1:-1] + un[:-2]))\n",
+ "/var/folders/l_/_p2_9s4x6sgbl7sh8tjclkyc0000gn/T/ipykernel_15673/1776893760.py:17: RuntimeWarning: invalid value encountered in add\n",
+ " + nu * dt/dx**2 * (un[2:] - 2*un[1:-1] + un[:-2]))\n",
+ "/var/folders/l_/_p2_9s4x6sgbl7sh8tjclkyc0000gn/T/ipykernel_15673/1776893760.py:15: RuntimeWarning: invalid value encountered in add\n",
+ " u[1:-1] = (un[1:-1]\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "nx = 101\n",
+ "nt = 200\n",
+ "nu = 0.07\n",
+ "dx = 2.0 / (nx - 1)\n",
+ "dt = 0.2 * dx # keep small; for shocks you may need smaller dt\n",
+ "\n",
+ "x = np.linspace(0, 2, nx)\n",
+ "u = np.ones(nx)\n",
+ "u[int(0.5 / dx):int(1.0 / dx) + 1] = 2.0\n",
+ "\n",
+ "plot1d(x, u, title=\"Step 4: Initial condition (1D Burgers)\")\n",
+ "\n",
+ "for n in range(nt):\n",
+ " un = u.copy()\n",
+ " u[1:-1] = (un[1:-1]\n",
+ " - un[1:-1] * dt/dx * (un[1:-1] - un[:-2])\n",
+ " + nu * dt/dx**2 * (un[2:] - 2*un[1:-1] + un[:-2]))\n",
+ " u[0] = 1.0\n",
+ " u[-1] = 1.0\n",
+ "\n",
+ "plot1d(x, u, title=f\"Step 4: After nt={nt} steps\")\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Step 5 — 2D Linear Convection\n",
+ "\n",
+ "We move to 2D transport:\n",
+ "\n",
+ "$$\n",
+ "\\frac{\\partial u}{\\partial t} + c\\,\\frac{\\partial u}{\\partial x} + c\\,\\frac{\\partial u}{\\partial y} = 0.\n",
+ "$$\n",
+ "\n",
+ "### Discretization (upwind in both x and y)\n",
+ "$$\n",
+ "u_{i,j}^{n+1} = u_{i,j}^n\n",
+ "- c\\,\\frac{\\Delta t}{\\Delta x}(u_{i,j}^n-u_{i-1,j}^n)\n",
+ "- c\\,\\frac{\\Delta t}{\\Delta y}(u_{i,j}^n-u_{i,j-1}^n).\n",
+ "$$\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 20,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "nx = 81\n",
+ "ny = 81\n",
+ "nt = 100\n",
+ "c = 1.0\n",
+ "\n",
+ "dx = 2.0 / (nx - 1)\n",
+ "dy = 2.0 / (ny - 1)\n",
+ "sigma = 0.2\n",
+ "dt = sigma * min(dx, dy) / c\n",
+ "\n",
+ "x = np.linspace(0, 2, nx)\n",
+ "y = np.linspace(0, 2, ny)\n",
+ "X, Y = np.meshgrid(x, y)\n",
+ "\n",
+ "u = np.ones((ny, nx))\n",
+ "u[int(0.5/dy):int(1.0/dy)+1, int(0.5/dx):int(1.0/dx)+1] = 2.0\n",
+ "\n",
+ "plot2d_surface(X, Y, u, title=\"Step 5: Initial condition (2D linear convection)\", zlabel=\"u\")\n",
+ "\n",
+ "for n in range(nt):\n",
+ " un = u.copy()\n",
+ " u[1:-1, 1:-1] = (un[1:-1, 1:-1]\n",
+ " - c * dt/dx * (un[1:-1, 1:-1] - un[1:-1, 0:-2])\n",
+ " - c * dt/dy * (un[1:-1, 1:-1] - un[0:-2, 1:-1]))\n",
+ "\n",
+ " # BCs\n",
+ " u[0, :] = 1\n",
+ " u[-1, :] = 1\n",
+ " u[:, 0] = 1\n",
+ " u[:, -1] = 1\n",
+ "\n",
+ "plot2d_surface(X, Y, u, title=f\"Step 5: After nt={nt} steps\", zlabel=\"u\")\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Step 6 — 2D Nonlinear Convection\n",
+ "\n",
+ "Now we advect **two velocity components** $u(x,y,t)$ and $v(x,y,t)$:\n",
+ "\n",
+ "\n",
+ "\\begin{aligned}\n",
+ "\\frac{\\partial u}{\\partial t} + u\\frac{\\partial u}{\\partial x} + v\\frac{\\partial u}{\\partial y} &= 0,\\\\\n",
+ "\\frac{\\partial v}{\\partial t} + u\\frac{\\partial v}{\\partial x} + v\\frac{\\partial v}{\\partial y} &= 0.\n",
+ "\\end{aligned}\n",
+ "\n",
+ "\n",
+ "This is essentially “nonlinear transport” in 2D. We use upwind differences.\n",
+ "\n",
+ "> **Note:** In this step there is **no pressure** and **no viscosity** yet.\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 21,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "nx = 81\n",
+ "ny = 81\n",
+ "nt = 100\n",
+ "\n",
+ "dx = 2.0 / (nx - 1)\n",
+ "dy = 2.0 / (ny - 1)\n",
+ "sigma = 0.2\n",
+ "dt = sigma * min(dx, dy)\n",
+ "\n",
+ "x = np.linspace(0, 2, nx)\n",
+ "y = np.linspace(0, 2, ny)\n",
+ "X, Y = np.meshgrid(x, y)\n",
+ "\n",
+ "u = np.ones((ny, nx))\n",
+ "v = np.ones((ny, nx))\n",
+ "\n",
+ "u[int(0.5/dy):int(1.0/dy)+1, int(0.5/dx):int(1.0/dx)+1] = 2.0\n",
+ "v[int(0.5/dy):int(1.0/dy)+1, int(0.5/dx):int(1.0/dx)+1] = 2.0\n",
+ "\n",
+ "plot2d_contour(X, Y, u, title=\"Step 6: u initial (2D nonlinear convection)\")\n",
+ "plot2d_contour(X, Y, v, title=\"Step 6: v initial (2D nonlinear convection)\")\n",
+ "\n",
+ "for n in range(nt):\n",
+ " un = u.copy()\n",
+ " vn = v.copy()\n",
+ "\n",
+ " u[1:-1, 1:-1] = (un[1:-1, 1:-1]\n",
+ " - (un[1:-1, 1:-1] * dt/dx) * (un[1:-1, 1:-1] - un[1:-1, 0:-2])\n",
+ " - (vn[1:-1, 1:-1] * dt/dy) * (un[1:-1, 1:-1] - un[0:-2, 1:-1]))\n",
+ "\n",
+ " v[1:-1, 1:-1] = (vn[1:-1, 1:-1]\n",
+ " - (un[1:-1, 1:-1] * dt/dx) * (vn[1:-1, 1:-1] - vn[1:-1, 0:-2])\n",
+ " - (vn[1:-1, 1:-1] * dt/dy) * (vn[1:-1, 1:-1] - vn[0:-2, 1:-1]))\n",
+ "\n",
+ " # BCs\n",
+ " u[0, :] = 1; u[-1, :] = 1; u[:, 0] = 1; u[:, -1] = 1\n",
+ " v[0, :] = 1; v[-1, :] = 1; v[:, 0] = 1; v[:, -1] = 1\n",
+ "\n",
+ "plot2d_contour(X, Y, u, title=f\"Step 6: u after nt={nt} steps\")\n",
+ "plot2d_contour(X, Y, v, title=f\"Step 6: v after nt={nt} steps\")\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Step 7 — 2D Diffusion\n",
+ "\n",
+ "Diffusion in 2D:\n",
+ "\n",
+ "$$\n",
+ "\\frac{\\partial u}{\\partial t} = \\nu\\left(\\frac{\\partial^2 u}{\\partial x^2}+\\frac{\\partial^2 u}{\\partial y^2}\\right).\n",
+ "$$\n",
+ "\n",
+ "### Discretization (explicit)\n",
+ "$$\n",
+ "u_{i,j}^{n+1} = u_{i,j}^{n} + \\nu\\Delta t\\left(\n",
+ "\\frac{u_{i+1,j}^n - 2u_{i,j}^n + u_{i-1,j}^n}{\\Delta x^2} +\n",
+ "\\frac{u_{i,j+1}^n - 2u_{i,j}^n + u_{i,j-1}^n}{\\Delta y^2}\n",
+ "\\right).\n",
+ "$$\n",
+ "\n",
+ "### Stability\n",
+ "A typical sufficient condition is:\n",
+ "\n",
+ "$$\n",
+ "\\nu\\Delta t\\left(\\frac{1}{\\Delta x^2}+\\frac{1}{\\Delta y^2}\\right) \\lesssim \\frac{1}{2}.\n",
+ "$$\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 22,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "nx = 81\n",
+ "ny = 81\n",
+ "nt = 100\n",
+ "nu = 0.1\n",
+ "\n",
+ "dx = 2.0 / (nx - 1)\n",
+ "dy = 2.0 / (ny - 1)\n",
+ "sigma = 0.2\n",
+ "dt = sigma * min(dx, dy)**2 / nu\n",
+ "\n",
+ "x = np.linspace(0, 2, nx)\n",
+ "y = np.linspace(0, 2, ny)\n",
+ "X, Y = np.meshgrid(x, y)\n",
+ "\n",
+ "u = np.ones((ny, nx))\n",
+ "u[int(0.5/dy):int(1.0/dy)+1, int(0.5/dx):int(1.0/dx)+1] = 2.0\n",
+ "\n",
+ "plot2d_surface(X, Y, u, title=\"Step 7: Initial condition (2D diffusion)\", zlabel=\"u\")\n",
+ "\n",
+ "for n in range(nt):\n",
+ " un = u.copy()\n",
+ " u[1:-1, 1:-1] = (un[1:-1, 1:-1]\n",
+ " + nu * dt/dx**2 * (un[1:-1, 2:] - 2*un[1:-1, 1:-1] + un[1:-1, 0:-2])\n",
+ " + nu * dt/dy**2 * (un[2:, 1:-1] - 2*un[1:-1, 1:-1] + un[0:-2, 1:-1]))\n",
+ " u[0, :] = 1; u[-1, :] = 1; u[:, 0] = 1; u[:, -1] = 1\n",
+ "\n",
+ "plot2d_surface(X, Y, u, title=f\"Step 7: After nt={nt} steps\", zlabel=\"u\")\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Step 8 — 2D Burgers’ Equations\n",
+ "\n",
+ "Add diffusion to the 2D nonlinear convection system:\n",
+ "\n",
+ "\n",
+ "\\begin{aligned}\n",
+ "\\frac{\\partial u}{\\partial t} + u\\frac{\\partial u}{\\partial x} + v\\frac{\\partial u}{\\partial y} &= \\nu\\left(\\frac{\\partial^2 u}{\\partial x^2}+\\frac{\\partial^2 u}{\\partial y^2}\\right),\\\\\n",
+ "\\frac{\\partial v}{\\partial t} + u\\frac{\\partial v}{\\partial x} + v\\frac{\\partial v}{\\partial y} &= \\nu\\left(\\frac{\\partial^2 v}{\\partial x^2}+\\frac{\\partial^2 v}{\\partial y^2}\\right).\n",
+ "\\end{aligned}\n",
+ "\n",
+ "\n",
+ "This is a classic nonlinear + dissipative system and a stepping stone toward Navier–Stokes.\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 23,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "nx = 61\n",
+ "ny = 61\n",
+ "nt = 100\n",
+ "nu = 0.01\n",
+ "\n",
+ "dx = 2.0 / (nx - 1)\n",
+ "dy = 2.0 / (ny - 1)\n",
+ "sigma = 0.2\n",
+ "dt = sigma * min(dx, dy)\n",
+ "\n",
+ "x = np.linspace(0, 2, nx)\n",
+ "y = np.linspace(0, 2, ny)\n",
+ "X, Y = np.meshgrid(x, y)\n",
+ "\n",
+ "u = np.ones((ny, nx))\n",
+ "v = np.ones((ny, nx))\n",
+ "u[int(0.5/dy):int(1.0/dy)+1, int(0.5/dx):int(1.0/dx)+1] = 2.0\n",
+ "v[int(0.5/dy):int(1.0/dy)+1, int(0.5/dx):int(1.0/dx)+1] = 2.0\n",
+ "\n",
+ "for n in range(nt):\n",
+ " un = u.copy()\n",
+ " vn = v.copy()\n",
+ "\n",
+ " u[1:-1, 1:-1] = (un[1:-1, 1:-1]\n",
+ " - (un[1:-1, 1:-1] * dt/dx) * (un[1:-1, 1:-1] - un[1:-1, 0:-2])\n",
+ " - (vn[1:-1, 1:-1] * dt/dy) * (un[1:-1, 1:-1] - un[0:-2, 1:-1])\n",
+ " + nu * dt/dx**2 * (un[1:-1, 2:] - 2*un[1:-1, 1:-1] + un[1:-1, 0:-2])\n",
+ " + nu * dt/dy**2 * (un[2:, 1:-1] - 2*un[1:-1, 1:-1] + un[0:-2, 1:-1]))\n",
+ "\n",
+ " v[1:-1, 1:-1] = (vn[1:-1, 1:-1]\n",
+ " - (un[1:-1, 1:-1] * dt/dx) * (vn[1:-1, 1:-1] - vn[1:-1, 0:-2])\n",
+ " - (vn[1:-1, 1:-1] * dt/dy) * (vn[1:-1, 1:-1] - vn[0:-2, 1:-1])\n",
+ " + nu * dt/dx**2 * (vn[1:-1, 2:] - 2*vn[1:-1, 1:-1] + vn[1:-1, 0:-2])\n",
+ " + nu * dt/dy**2 * (vn[2:, 1:-1] - 2*vn[1:-1, 1:-1] + vn[0:-2, 1:-1]))\n",
+ "\n",
+ " u[0, :] = 1; u[-1, :] = 1; u[:, 0] = 1; u[:, -1] = 1\n",
+ " v[0, :] = 1; v[-1, :] = 1; v[:, 0] = 1; v[:, -1] = 1\n",
+ "\n",
+ "plot2d_contour(X, Y, u, title=\"Step 8: u (2D Burgers) after evolution\")\n",
+ "plot2d_contour(X, Y, v, title=\"Step 8: v (2D Burgers) after evolution\")\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Step 9 — Laplace Equation\n",
+ "\n",
+ "Laplace’s equation appears in steady-state diffusion, potential flow, electrostatics, etc.:\n",
+ "\n",
+ "$$\n",
+ "\\frac{\\partial^2 p}{\\partial x^2} + \\frac{\\partial^2 p}{\\partial y^2} = 0.\n",
+ "$$\n",
+ "\n",
+ "We solve it iteratively using a simple Jacobi-like update.\n",
+ "\n",
+ "### Discretization\n",
+ "$$\n",
+ "p_{i,j} = \\frac{1}{2(\\Delta x^2+\\Delta y^2)}\\Big[\n",
+ "\\Delta y^2(p_{i+1,j}+p_{i-1,j}) + \\Delta x^2(p_{i,j+1}+p_{i,j-1})\n",
+ "\\Big].\n",
+ "$$\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 24,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "nx = 51\n",
+ "ny = 51\n",
+ "\n",
+ "dx = 2.0 / (nx - 1)\n",
+ "dy = 1.0 / (ny - 1)\n",
+ "\n",
+ "x = np.linspace(0, 2, nx)\n",
+ "y = np.linspace(0, 1, ny)\n",
+ "X, Y = np.meshgrid(x, y)\n",
+ "\n",
+ "p = np.zeros((ny, nx))\n",
+ "\n",
+ "# Boundary conditions (example)\n",
+ "# p = 0 at x=0, p = y at x=2, dp/dy = 0 at y=0 and y=1\n",
+ "p[:, 0] = 0\n",
+ "p[:, -1] = y\n",
+ "\n",
+ "def laplace2d(p, nit=500):\n",
+ " pn = np.empty_like(p)\n",
+ " for it in range(nit):\n",
+ " pn[:] = p\n",
+ "\n",
+ " p[1:-1, 1:-1] = ((dy**2 * (pn[1:-1, 2:] + pn[1:-1, 0:-2]) +\n",
+ " dx**2 * (pn[2:, 1:-1] + pn[0:-2, 1:-1]))\n",
+ " / (2*(dx**2 + dy**2)))\n",
+ "\n",
+ " # Re-apply BCs every iteration\n",
+ " p[:, 0] = 0\n",
+ " p[:, -1] = y\n",
+ " p[0, :] = p[1, :] # dp/dy=0 at y=0\n",
+ " p[-1, :] = p[-2, :] # dp/dy=0 at y=1\n",
+ "\n",
+ " return p\n",
+ "\n",
+ "p = laplace2d(p, nit=800)\n",
+ "plot2d_contour(X, Y, p, title=\"Step 9: Laplace solution p(x,y)\", cmap=cm.viridis)\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Step 10 — Poisson Equation\n",
+ "\n",
+ "Poisson’s equation adds a source term \\(b(x,y)\\):\n",
+ "\n",
+ "$$\n",
+ "\\frac{\\partial^2 p}{\\partial x^2} + \\frac{\\partial^2 p}{\\partial y^2} = b.\n",
+ "$$\n",
+ "\n",
+ "In incompressible CFD, a Poisson equation for pressure is central to enforcing $\\nabla\\cdot\\mathbf{u}=0$.\n",
+ "\n",
+ "We again use an iterative solver (Jacobi-like).\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 25,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "nx = 51\n",
+ "ny = 51\n",
+ "dx = 2.0 / (nx - 1)\n",
+ "dy = 1.0 / (ny - 1)\n",
+ "\n",
+ "x = np.linspace(0, 2, nx)\n",
+ "y = np.linspace(0, 1, ny)\n",
+ "X, Y = np.meshgrid(x, y)\n",
+ "\n",
+ "p = np.zeros((ny, nx))\n",
+ "b = np.zeros((ny, nx))\n",
+ "\n",
+ "# Example point sources\n",
+ "b[int(ny/4), int(nx/4)] = 100\n",
+ "b[int(3*ny/4), int(3*nx/4)] = -100\n",
+ "\n",
+ "def poisson2d(p, b, nit=800):\n",
+ " pn = np.empty_like(p)\n",
+ " for it in range(nit):\n",
+ " pn[:] = p\n",
+ " p[1:-1, 1:-1] = (((pn[1:-1, 2:] + pn[1:-1, 0:-2]) * dy**2 +\n",
+ " (pn[2:, 1:-1] + pn[0:-2, 1:-1]) * dx**2 -\n",
+ " b[1:-1, 1:-1] * dx**2 * dy**2)\n",
+ " / (2*(dx**2 + dy**2)))\n",
+ "\n",
+ " # Dirichlet BCs (example)\n",
+ " p[0, :] = 0\n",
+ " p[-1, :] = 0\n",
+ " p[:, 0] = 0\n",
+ " p[:, -1] = 0\n",
+ " return p\n",
+ "\n",
+ "p = poisson2d(p, b, nit=1200)\n",
+ "plot2d_contour(X, Y, p, title=\"Step 10: Poisson solution p(x,y)\", cmap=cm.viridis)\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Step 11 — Incompressible Navier–Stokes: Lid-Driven Cavity\n",
+ "\n",
+ "Now we solve the **2D incompressible Navier–Stokes equations**:\n",
+ "\n",
+ "\n",
+ "\\begin{aligned}\n",
+ "\\frac{\\partial u}{\\partial t} + u\\frac{\\partial u}{\\partial x} + v\\frac{\\partial u}{\\partial y} &= -\\frac{1}{\\rho}\\frac{\\partial p}{\\partial x} + \\nu\\left(\\frac{\\partial^2 u}{\\partial x^2}+\\frac{\\partial^2 u}{\\partial y^2}\\right),\\\\\n",
+ "\\frac{\\partial v}{\\partial t} + u\\frac{\\partial v}{\\partial x} + v\\frac{\\partial v}{\\partial y} &= -\\frac{1}{\\rho}\\frac{\\partial p}{\\partial y} + \\nu\\left(\\frac{\\partial^2 v}{\\partial x^2}+\\frac{\\partial^2 v}{\\partial y^2}\\right),\\\\\n",
+ "\\frac{\\partial u}{\\partial x} + \\frac{\\partial v}{\\partial y} &= 0.\n",
+ "\\end{aligned}\n",
+ "\n",
+ "\n",
+ "### Projection idea\n",
+ "1. Predict a tentative velocity from advection + diffusion (ignoring pressure).\n",
+ "2. Solve a **pressure Poisson equation** derived from enforcing \\(\\nabla\\cdot\\mathbf{u}=0\\).\n",
+ "3. Correct velocity using the pressure gradient.\n",
+ "\n",
+ "### Boundary conditions (lid-driven cavity)\n",
+ "- Velocity is zero on all walls **except** the top lid:\n",
+ " - `u = 1` at `y = 2` (top)\n",
+ " - `u = v = 0` elsewhere\n",
+ "\n",
+ "This produces a circulating vortex and is a standard benchmark.\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 26,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "nx = 41\n",
+ "ny = 41\n",
+ "nt = 300 # time steps\n",
+ "nit = 50 # pressure iterations each step\n",
+ "\n",
+ "dx = 2.0 / (nx - 1)\n",
+ "dy = 2.0 / (ny - 1)\n",
+ "\n",
+ "rho = 1.0\n",
+ "nu = 0.1\n",
+ "dt = 0.001\n",
+ "\n",
+ "x = np.linspace(0, 2, nx)\n",
+ "y = np.linspace(0, 2, ny)\n",
+ "X, Y = np.meshgrid(x, y)\n",
+ "\n",
+ "u = np.zeros((ny, nx))\n",
+ "v = np.zeros((ny, nx))\n",
+ "p = np.zeros((ny, nx))\n",
+ "b = np.zeros((ny, nx))\n",
+ "\n",
+ "def build_up_b(b, rho, dt, u, v, dx, dy):\n",
+ " \"\"\"Build the RHS of the pressure Poisson equation.\"\"\"\n",
+ " b[1:-1, 1:-1] = (rho * (1/dt * (\n",
+ " (u[1:-1, 2:] - u[1:-1, 0:-2]) / (2*dx) +\n",
+ " (v[2:, 1:-1] - v[0:-2, 1:-1]) / (2*dy)\n",
+ " ) -\n",
+ " ((u[1:-1, 2:] - u[1:-1, 0:-2]) / (2*dx))**2 -\n",
+ " 2 * ((u[2:, 1:-1] - u[0:-2, 1:-1]) / (2*dy) *\n",
+ " (v[1:-1, 2:] - v[1:-1, 0:-2]) / (2*dx)) -\n",
+ " ((v[2:, 1:-1] - v[0:-2, 1:-1]) / (2*dy))**2))\n",
+ " return b\n",
+ "\n",
+ "def pressure_poisson(p, dx, dy, b, nit):\n",
+ " pn = np.empty_like(p)\n",
+ " for it in range(nit):\n",
+ " pn[:] = p\n",
+ " p[1:-1, 1:-1] = (((pn[1:-1, 2:] + pn[1:-1, 0:-2]) * dy**2 +\n",
+ " (pn[2:, 1:-1] + pn[0:-2, 1:-1]) * dx**2) /\n",
+ " (2*(dx**2 + dy**2)) -\n",
+ " (dx**2 * dy**2)/(2*(dx**2 + dy**2)) * b[1:-1, 1:-1])\n",
+ "\n",
+ " # Pressure BCs for cavity:\n",
+ " p[:, -1] = p[:, -2] # dp/dx = 0 at x=2\n",
+ " p[:, 0] = p[:, 1] # dp/dx = 0 at x=0\n",
+ " p[0, :] = p[1, :] # dp/dy = 0 at y=0\n",
+ " p[-1, :] = 0 # p = 0 at y=2 (reference)\n",
+ " return p\n",
+ "\n",
+ "for n in range(nt):\n",
+ " un = u.copy()\n",
+ " vn = v.copy()\n",
+ "\n",
+ " b = build_up_b(b, rho, dt, u, v, dx, dy)\n",
+ " p = pressure_poisson(p, dx, dy, b, nit)\n",
+ "\n",
+ " # Velocity update (explicit)\n",
+ " u[1:-1, 1:-1] = (un[1:-1, 1:-1]\n",
+ " - un[1:-1, 1:-1] * dt/dx * (un[1:-1, 1:-1] - un[1:-1, 0:-2])\n",
+ " - vn[1:-1, 1:-1] * dt/dy * (un[1:-1, 1:-1] - un[0:-2, 1:-1])\n",
+ " - dt/(2*rho*dx) * (p[1:-1, 2:] - p[1:-1, 0:-2])\n",
+ " + nu * (dt/dx**2 * (un[1:-1, 2:] - 2*un[1:-1, 1:-1] + un[1:-1, 0:-2]) +\n",
+ " dt/dy**2 * (un[2:, 1:-1] - 2*un[1:-1, 1:-1] + un[0:-2, 1:-1])))\n",
+ "\n",
+ " v[1:-1, 1:-1] = (vn[1:-1, 1:-1]\n",
+ " - un[1:-1, 1:-1] * dt/dx * (vn[1:-1, 1:-1] - vn[1:-1, 0:-2])\n",
+ " - vn[1:-1, 1:-1] * dt/dy * (vn[1:-1, 1:-1] - vn[0:-2, 1:-1])\n",
+ " - dt/(2*rho*dy) * (p[2:, 1:-1] - p[0:-2, 1:-1])\n",
+ " + nu * (dt/dx**2 * (vn[1:-1, 2:] - 2*vn[1:-1, 1:-1] + vn[1:-1, 0:-2]) +\n",
+ " dt/dy**2 * (vn[2:, 1:-1] - 2*vn[1:-1, 1:-1] + vn[0:-2, 1:-1])))\n",
+ "\n",
+ " # Velocity BCs\n",
+ " u[0, :] = 0\n",
+ " u[:, 0] = 0\n",
+ " u[:, -1] = 0\n",
+ " u[-1, :] = 1 # moving lid\n",
+ "\n",
+ " v[0, :] = 0\n",
+ " v[-1, :] = 0\n",
+ " v[:, 0] = 0\n",
+ " v[:, -1] = 0\n",
+ "\n",
+ "# Visualize: pressure + velocity field\n",
+ "plot2d_contour(X, Y, p, title=\"Step 11: Pressure field p(x,y) in cavity\")\n",
+ "\n",
+ "fig, ax = plt.subplots(figsize=(8, 6), dpi=110)\n",
+ "ax.contourf(X, Y, p, levels=50, cmap=cm.viridis, alpha=0.8)\n",
+ "ax.streamplot(X, Y, u, v, density=1.6, linewidth=1)\n",
+ "ax.set_title(\"Step 11: Lid-driven cavity — pressure contours + streamlines\")\n",
+ "ax.set_xlabel(\"x\")\n",
+ "ax.set_ylabel(\"y\")\n",
+ "ax.set_aspect(\"equal\")\n",
+ "fig.tight_layout()\n",
+ "plt.show()\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Step 12 — Incompressible Navier–Stokes: Channel Flow (Pressure-Driven)\n",
+ "\n",
+ "We solve the same incompressible Navier–Stokes system, but now in a channel with a constant pressure gradient\n",
+ "(or equivalently a body force) driving flow.\n",
+ "\n",
+ "A common pedagogical setup:\n",
+ "- Periodic-like condition in x can be approximated by enforcing pressure BCs with a fixed $\\Delta p$ across the domain.\n",
+ "- No-slip on the top and bottom walls.\n",
+ "\n",
+ "Below is a compact explicit solver with pressure Poisson iterations, following the same projection idea as Step 11.\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 27,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "nx = 41\n",
+ "ny = 41\n",
+ "nt = 500\n",
+ "nit = 50\n",
+ "\n",
+ "dx = 2.0 / (nx - 1)\n",
+ "dy = 2.0 / (ny - 1)\n",
+ "\n",
+ "rho = 1.0\n",
+ "nu = 0.1\n",
+ "dt = 0.001\n",
+ "\n",
+ "F = 1.0 # body force driving the flow in x-direction (acts like a pressure gradient)\n",
+ "\n",
+ "x = np.linspace(0, 2, nx)\n",
+ "y = np.linspace(0, 2, ny)\n",
+ "X, Y = np.meshgrid(x, y)\n",
+ "\n",
+ "u = np.zeros((ny, nx))\n",
+ "v = np.zeros((ny, nx))\n",
+ "p = np.zeros((ny, nx))\n",
+ "b = np.zeros((ny, nx))\n",
+ "\n",
+ "def build_up_b_channel(b, rho, dt, u, v, dx, dy):\n",
+ " b[1:-1, 1:-1] = (rho * (1/dt * (\n",
+ " (u[1:-1, 2:] - u[1:-1, 0:-2]) / (2*dx) +\n",
+ " (v[2:, 1:-1] - v[0:-2, 1:-1]) / (2*dy)\n",
+ " ) -\n",
+ " ((u[1:-1, 2:] - u[1:-1, 0:-2]) / (2*dx))**2 -\n",
+ " 2 * ((u[2:, 1:-1] - u[0:-2, 1:-1]) / (2*dy) *\n",
+ " (v[1:-1, 2:] - v[1:-1, 0:-2]) / (2*dx)) -\n",
+ " ((v[2:, 1:-1] - v[0:-2, 1:-1]) / (2*dy))**2))\n",
+ " return b\n",
+ "\n",
+ "def pressure_poisson_channel(p, dx, dy, b, nit):\n",
+ " pn = np.empty_like(p)\n",
+ " for it in range(nit):\n",
+ " pn[:] = p\n",
+ " p[1:-1, 1:-1] = (((pn[1:-1, 2:] + pn[1:-1, 0:-2]) * dy**2 +\n",
+ " (pn[2:, 1:-1] + pn[0:-2, 1:-1]) * dx**2) /\n",
+ " (2*(dx**2 + dy**2)) -\n",
+ " (dx**2 * dy**2)/(2*(dx**2 + dy**2)) * b[1:-1, 1:-1])\n",
+ "\n",
+ " # Pressure BCs for channel (approximate periodic in x)\n",
+ " p[:, -1] = p[:, -2] # dp/dx = 0 at x=2\n",
+ " p[:, 0] = p[:, 1] # dp/dx = 0 at x=0\n",
+ " p[0, :] = p[1, :] # dp/dy = 0 at y=0\n",
+ " p[-1, :] = p[-2, :] # dp/dy = 0 at y=2\n",
+ " return p\n",
+ "\n",
+ "for n in range(nt):\n",
+ " un = u.copy()\n",
+ " vn = v.copy()\n",
+ "\n",
+ " b = build_up_b_channel(b, rho, dt, u, v, dx, dy)\n",
+ " p = pressure_poisson_channel(p, dx, dy, b, nit)\n",
+ "\n",
+ " u[1:-1, 1:-1] = (un[1:-1, 1:-1]\n",
+ " - un[1:-1, 1:-1] * dt/dx * (un[1:-1, 1:-1] - un[1:-1, 0:-2])\n",
+ " - vn[1:-1, 1:-1] * dt/dy * (un[1:-1, 1:-1] - un[0:-2, 1:-1])\n",
+ " - dt/(2*rho*dx) * (p[1:-1, 2:] - p[1:-1, 0:-2])\n",
+ " + nu * (dt/dx**2 * (un[1:-1, 2:] - 2*un[1:-1, 1:-1] + un[1:-1, 0:-2]) +\n",
+ " dt/dy**2 * (un[2:, 1:-1] - 2*un[1:-1, 1:-1] + un[0:-2, 1:-1]))\n",
+ " + F * dt)\n",
+ "\n",
+ " v[1:-1, 1:-1] = (vn[1:-1, 1:-1]\n",
+ " - un[1:-1, 1:-1] * dt/dx * (vn[1:-1, 1:-1] - vn[1:-1, 0:-2])\n",
+ " - vn[1:-1, 1:-1] * dt/dy * (vn[1:-1, 1:-1] - vn[0:-2, 1:-1])\n",
+ " - dt/(2*rho*dy) * (p[2:, 1:-1] - p[0:-2, 1:-1])\n",
+ " + nu * (dt/dx**2 * (vn[1:-1, 2:] - 2*vn[1:-1, 1:-1] + vn[1:-1, 0:-2]) +\n",
+ " dt/dy**2 * (vn[2:, 1:-1] - 2*vn[1:-1, 1:-1] + vn[0:-2, 1:-1])))\n",
+ "\n",
+ " # Velocity BCs: no-slip at y=0 and y=2\n",
+ " u[0, :] = 0\n",
+ " u[-1, :] = 0\n",
+ " v[0, :] = 0\n",
+ " v[-1, :] = 0\n",
+ "\n",
+ " # \"Periodic-ish\" in x: copy from near-boundary (simple approximation)\n",
+ " u[:, 0] = u[:, 1]\n",
+ " u[:, -1] = u[:, -2]\n",
+ " v[:, 0] = v[:, 1]\n",
+ " v[:, -1] = v[:, -2]\n",
+ "\n",
+ "plot2d_contour(X, Y, u, title=\"Step 12: Channel flow — u(x,y)\")\n",
+ "fig, ax = plt.subplots(figsize=(8, 6), dpi=110)\n",
+ "ax.contourf(X, Y, u, levels=50, cmap=cm.viridis)\n",
+ "ax.streamplot(X, Y, u, v, density=1.5, linewidth=1)\n",
+ "ax.set_title(\"Step 12: Channel flow — u contours + streamlines\")\n",
+ "ax.set_xlabel(\"x\"); ax.set_ylabel(\"y\")\n",
+ "ax.set_aspect(\"equal\")\n",
+ "fig.tight_layout()\n",
+ "plt.show()\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "---\n",
+ "\n",
+ "## Notes, extensions, and “next steps”\n",
+ "- Replace Forward Euler with a higher-order time integrator (RK2/RK3).\n",
+ "- Use higher-order spatial schemes (QUICK/WENO) for sharper convection.\n",
+ "- Swap the simple pressure solver for a faster method (SOR, multigrid).\n",
+ "- Validate with known solutions (Taylor–Green vortex, Poiseuille profile, etc.).\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "3a4f365b",
+ "metadata": {},
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "base",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.10.13"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}