Refactor book code to use Quarto includes with tested source files

- Add callout notes to chapters pointing to tested source implementations
- Extract inline Python code to src/book_snippets/ with RESULT variables for testing
- Create snippet wrapper .qmd files in chapters/*/snippets/ directories
- Add comprehensive test suite in tests/test_book_snippets.py (21 tests)
- Create src/nonlin/split_logistic.py and split_diffu_react.py modules

Chapters updated with callout notes:
- elliptic/elliptic.qmd -> src/elliptic/
- vib/vib_undamped.qmd -> src/vib/
- wave/wave1D_fd1.qmd, wave1D_fd2.qmd -> src/wave/
- diffu/diffu_fd1.qmd, diffu_fd3.qmd, diffu_rw.qmd -> src/diffu/
- appendices/softeng2/softeng2.qmd -> src/softeng2/

Chapters refactored to use includes:
- devito_intro/ (what_is_devito, first_pde, boundary_conditions, verification)
- nonlin/ (burgers, nonlin_ode, nonlin_split)
- advec/advec1D_devito.qmd

Author: ggorman

chapters/advec/advec1D_devito.qmd

@@ -70,46 +70,11 @@ $$
 u^{n+1}_i = u^n_i - C(u^n_i - u^n_{i-1})
 $$ {#eq-advec-upwind-update}
 
-In Devito, we express this using shifted indexing:
+In Devito, we express this using shifted indexing. The key technique is
+using `u.subs(x_dim, x_dim - x_dim.spacing)` to create a reference to
+$u^n_{i-1}$:
 
-```python
-from devito import Grid, TimeFunction, Eq, Operator, Constant
-import numpy as np
-
-def solve_advection_upwind(L, c, Nx, T, C, I):
-    """Upwind scheme for 1D advection."""
-    # Grid setup
-    dx = L / Nx
-    dt = C * dx / c
-
-    grid = Grid(shape=(Nx + 1,), extent=(L,))
-    x_dim, = grid.dimensions
-
-    u = TimeFunction(name='u', grid=grid, time_order=1, space_order=1)
-
-    # Set initial condition
-    x_coords = np.linspace(0, L, Nx + 1)
-    u.data[0, :] = I(x_coords)
-
-    # Courant number as constant
-    courant = Constant(name='C', value=C)
-
-    # Upwind stencil: u^{n+1} = u - C*(u - u[x-dx])
-    u_minus = u.subs(x_dim, x_dim - x_dim.spacing)
-    stencil = u - courant * (u - u_minus)
-    update = Eq(u.forward, stencil)
-
-    op = Operator([update])
-    # ... time stepping loop
-```
-
-The key line is:
-```python
-u_minus = u.subs(x_dim, x_dim - x_dim.spacing)
-```
-
-This creates a reference to $u^n_{i-1}$ by substituting `x_dim - x_dim.spacing`
-for `x_dim` in the `TimeFunction` `u`.
+{{< include snippets/advec_upwind.qmd >}}
 
 ### Lax-Wendroff Scheme Implementation
 
@@ -120,36 +85,11 @@ $$
 u^{n+1}_i = u^n_i - \frac{C}{2}(u^n_{i+1} - u^n_{i-1}) + \frac{C^2}{2}(u^n_{i+1} - 2u^n_i + u^n_{i-1})
 $$
 
-This can be written using Devito's derivative operators:
-
-```python
-def solve_advection_lax_wendroff(L, c, Nx, T, C, I):
-    """Lax-Wendroff scheme for 1D advection."""
-    dx = L / Nx
-    dt = C * dx / c
-
-    grid = Grid(shape=(Nx + 1,), extent=(L,))
-    u = TimeFunction(name='u', grid=grid, time_order=1, space_order=2)
-
-    x_coords = np.linspace(0, L, Nx + 1)
-    u.data[0, :] = I(x_coords)
-
-    courant = Constant(name='C', value=C)
-
-    # Lax-Wendroff: u - (C/2)*dx*u.dx + (C²/2)*dx²*u.dx2
-    # u.dx  = centered first derivative
-    # u.dx2 = centered second derivative
-    stencil = u - 0.5*courant*dx*u.dx + 0.5*courant**2*dx**2*u.dx2
-    update = Eq(u.forward, stencil)
-
-    op = Operator([update])
-    # ... time stepping loop
-```
-
-Here we use Devito's built-in derivative operators:
+This can be written using Devito's built-in derivative operators where
+`u.dx` computes the centered first derivative and `u.dx2` computes the
+centered second derivative:
 
-- `u.dx` computes the centered first derivative $(u_{i+1} - u_{i-1})/(2\Delta x)$
-- `u.dx2` computes the centered second derivative $(u_{i+1} - 2u_i + u_{i-1})/\Delta x^2$
+{{< include snippets/advec_lax_wendroff.qmd >}}
 
 ### Lax-Friedrichs Scheme Implementation
 

chapters/advec/snippets/advec_lax_wendroff.qmd

@@ -0,0 +1,7 @@
+::: {.callout-note title="Snippet (tested)"}
+This code is included from `src/book_snippets/advec_lax_wendroff.py` and exercised by `pytest` (see `tests/test_book_snippets.py`).
+:::
+
+```python
+{{< include ../../src/book_snippets/advec_lax_wendroff.py >}}
+```

chapters/advec/snippets/advec_upwind.qmd

@@ -0,0 +1,7 @@
+::: {.callout-note title="Snippet (tested)"}
+This code is included from `src/book_snippets/advec_upwind.py` and exercised by `pytest` (see `tests/test_book_snippets.py`).
+:::
+
+```python
+{{< include ../../src/book_snippets/advec_upwind.py >}}
+```

chapters/appendices/softeng2/softeng2.qmd

@@ -41,6 +41,20 @@ and @sec-wave-pde2-var-c.
 
 ## A solver function
 
+::: {.callout-note}
+## Source Files
+
+The software engineering patterns presented in this appendix (classes like
+`Storage`, `Parameters`, `Problem`, `Mesh`, `Function`, `Solver`) are
+implemented in `src/softeng2/`. Key files include:
+
+- `src/softeng2/Storage.py` - Data persistence with joblib
+- `src/softeng2/wave1D_oo.py` - Object-oriented wave solver with `Parameters` class
+- `src/softeng2/wave2D_u0.py` - 2D wave implementations
+
+These serve as reference implementations for the patterns discussed.
+:::
+
 The general initial-boundary value problem
 solved by finite difference methods. For modern implementations using
 Devito, see @sec-wave-devito. This section covers software engineering

chapters/devito_intro/boundary_conditions.qmd

@@ -15,41 +15,11 @@ $$
 
 The most direct approach adds equations that set boundary values:
 
-```python
-from devito import Grid, TimeFunction, Eq, Operator
-
-grid = Grid(shape=(101,), extent=(1.0,))
-u = TimeFunction(name='u', grid=grid, time_order=2, space_order=2)
-
-# Get the time dimension for indexing
-t = grid.stepping_dim
-
-# Interior update (wave equation)
-update = Eq(u.forward, 2*u - u.backward + dt**2 * c**2 * u.dx2)
-
-# Boundary conditions: u = 0 at both ends
-bc_left = Eq(u[t+1, 0], 0)
-bc_right = Eq(u[t+1, 100], 0)
-
-# Include all equations in the operator
-op = Operator([update, bc_left, bc_right])
-```
+{{< include snippets/boundary_dirichlet_wave.qmd >}}
 
 **Method 2: Using subdomain**
 
-For interior-only updates, use `subdomain=grid.interior`:
-
-```python
-# Update only interior points (automatically excludes boundaries)
-update = Eq(u.forward, 2*u - u.backward + dt**2 * c**2 * u.dx2,
-            subdomain=grid.interior)
-
-# Set boundaries explicitly
-bc_left = Eq(u[t+1, 0], 0)
-bc_right = Eq(u[t+1, 100], 0)
-
-op = Operator([update, bc_left, bc_right])
-```
+The snippet above already uses `subdomain=grid.interior` to keep the interior PDE update separate from boundary treatment.
 
 The `subdomain=grid.interior` approach is often cleaner because it
 explicitly separates the physics (interior PDE) from the boundary treatment.
@@ -71,87 +41,25 @@ $$
 
 This gives $u_{-1} = u_1$, which we substitute into the interior equation:
 
-```python
-grid = Grid(shape=(101,), extent=(1.0,))
-u = TimeFunction(name='u', grid=grid, time_order=1, space_order=2)
-x = grid.dimensions[0]
-t = grid.stepping_dim
-
-# Interior update (diffusion equation)
-update = Eq(u.forward, u + alpha * dt * u.dx2, subdomain=grid.interior)
-
-# Neumann BC at left (du/dx = 0): use one-sided update
-# u_new[0] = u[0] + alpha*dt * 2*(u[1] - u[0])/dx^2
-dx = grid.spacing[0]
-bc_left = Eq(u[t+1, 0], u[t, 0] + alpha * dt * 2 * (u[t, 1] - u[t, 0]) / dx**2)
-
-# Neumann BC at right (du/dx = 0)
-bc_right = Eq(u[t+1, 100], u[t, 100] + alpha * dt * 2 * (u[t, 99] - u[t, 100]) / dx**2)
-
-op = Operator([update, bc_left, bc_right])
-```
+{{< include snippets/neumann_bc_diffusion_1d.qmd >}}
 
 ### Mixed Boundary Conditions
 
 Often we have different conditions on different boundaries:
 
-```python
-# Dirichlet on left, Neumann on right
-bc_left = Eq(u[t+1, 0], 0)  # u(0,t) = 0
-bc_right = Eq(u[t+1, 100], u[t+1, 99])  # du/dx(L,t) = 0 (copy from interior)
-
-op = Operator([update, bc_left, bc_right])
-```
+{{< include snippets/mixed_bc_diffusion_1d.qmd >}}
 
 ### 2D Boundary Conditions
 
 For 2D problems, boundary conditions apply to all four edges:
 
-```python
-grid = Grid(shape=(101, 101), extent=(1.0, 1.0))
-u = TimeFunction(name='u', grid=grid, time_order=2, space_order=2)
-
-x, y = grid.dimensions
-t = grid.stepping_dim
-Nx, Ny = 100, 100
-
-# Interior update
-update = Eq(u.forward, 2*u - u.backward + dt**2 * c**2 * u.laplace,
-            subdomain=grid.interior)
-
-# Dirichlet BCs on all four edges
-bc_left = Eq(u[t+1, 0, y], 0)
-bc_right = Eq(u[t+1, Nx, y], 0)
-bc_bottom = Eq(u[t+1, x, 0], 0)
-bc_top = Eq(u[t+1, x, Ny], 0)
-
-op = Operator([update, bc_left, bc_right, bc_bottom, bc_top])
-```
+{{< include snippets/bc_2d_dirichlet_wave.qmd >}}
 
 ### Time-Dependent Boundary Conditions
 
 For boundaries that vary in time, use the time index:
 
-```python
-from devito import Constant
-
-# Time-varying amplitude
-A = Constant(name='A')
-
-# Sinusoidal forcing at left boundary
-# u(0, t) = A * sin(omega * t)
-import sympy as sp
-omega = 2 * sp.pi  # Angular frequency
-
-# The time value at step n
-t_val = t * dt  # Symbolic time value
-
-bc_left = Eq(u[t+1, 0], A * sp.sin(omega * t_val))
-
-# Set the amplitude before running
-op = Operator([update, bc_left, bc_right])
-op(time=Nt, dt=dt, A=1.0)  # Pass A as keyword argument
-```
+{{< include snippets/time_dependent_bc_sine.qmd >}}
 
 ### Absorbing Boundary Conditions
 
@@ -163,14 +71,7 @@ $$
 
 This can be discretized as:
 
-```python
-# Absorbing BC at right boundary (waves traveling right)
-dx = grid.spacing[0]
-bc_right_absorbing = Eq(
-    u[t+1, Nx],
-    u[t, Nx] - c * dt / dx * (u[t, Nx] - u[t, Nx-1])
-)
-```
+{{< include snippets/absorbing_bc_right_wave.qmd >}}
 
 More sophisticated absorbing conditions use damping layers (sponges)
 near the boundaries. This is covered in detail in @sec-wave-1d-absorbing.
@@ -185,11 +86,7 @@ $$
 Devito doesn't directly support periodic BCs, but they can be implemented
 by copying values:
 
-```python
-# Periodic BCs: u[0] = u[Nx-1], u[Nx] = u[1]
-bc_periodic_left = Eq(u[t+1, 0], u[t+1, Nx-1])
-bc_periodic_right = Eq(u[t+1, Nx], u[t+1, 1])
-```
+{{< include snippets/periodic_bc_advection_1d.qmd >}}
 
 Note: The order of equations matters. Update the interior first, then
 copy for periodicity.
@@ -212,42 +109,7 @@ copy for periodicity.
 
 Here's a complete example combining interior updates with boundary conditions:
 
-```python
-from devito import Grid, TimeFunction, Eq, Operator
-import numpy as np
-
-# Setup
-L, c, T = 1.0, 1.0, 2.0
-Nx = 100
-C = 0.9  # Courant number
-dx = L / Nx
-dt = C * dx / c
-Nt = int(T / dt)
-
-# Grid and field
-grid = Grid(shape=(Nx + 1,), extent=(L,))
-u = TimeFunction(name='u', grid=grid, time_order=2, space_order=2)
-t = grid.stepping_dim
-
-# Initial condition: plucked string
-x_vals = np.linspace(0, L, Nx + 1)
-u.data[0, :] = np.sin(np.pi * x_vals)
-u.data[1, :] = u.data[0, :]  # Zero initial velocity
-
-# Equations
-update = Eq(u.forward, 2*u - u.backward + (c*dt)**2 * u.dx2,
-            subdomain=grid.interior)
-bc_left = Eq(u[t+1, 0], 0)
-bc_right = Eq(u[t+1, Nx], 0)
-
-# Solve
-op = Operator([update, bc_left, bc_right])
-op(time=Nt, dt=dt)
-
-# Verify: solution should return to initial shape at t = 2L/c
-print(f"Initial max: {np.max(u.data[1, :]):.6f}")
-print(f"Final max: {np.max(u.data[0, :]):.6f}")
-```
+{{< include snippets/boundary_dirichlet_wave.qmd >}}
 
 For a string with fixed ends and initial shape $\sin(\pi x)$, the solution
 oscillates with period $2L/c$. After one period, it should return to the

chapters/devito_intro/first_pde.qmd

@@ -43,52 +43,7 @@ for $C \le 1$.
 
 Let's implement this step by step:
 
-```python
-from devito import Grid, TimeFunction, Eq, Operator
-import numpy as np
-
-# Problem parameters
-L = 1.0       # Domain length
-c = 1.0       # Wave speed
-T = 1.0       # Final time
-Nx = 100      # Number of grid points
-C = 0.5       # Courant number (for stability)
-
-# Derived parameters
-dx = L / Nx
-dt = C * dx / c
-Nt = int(T / dt)
-
-# Create the computational grid
-grid = Grid(shape=(Nx + 1,), extent=(L,))
-
-# Create a time-varying field
-# time_order=2 because we have second derivative in time
-# space_order=2 for standard second-order accuracy
-u = TimeFunction(name='u', grid=grid, time_order=2, space_order=2)
-
-# Set initial condition: a Gaussian pulse
-x = grid.dimensions[0]
-x_coord = 0.5 * L  # Center of domain
-sigma = 0.1        # Width of pulse
-u.data[0, :] = np.exp(-((np.linspace(0, L, Nx+1) - x_coord)**2) / (2*sigma**2))
-u.data[1, :] = u.data[0, :]  # Zero initial velocity
-
-# Define the update equation
-# u.forward is u at time n+1, u is at time n, u.backward is at time n-1
-# u.dx2 is the second spatial derivative
-eq = Eq(u.forward, 2*u - u.backward + (c*dt)**2 * u.dx2)
-
-# Create the operator
-op = Operator([eq])
-
-# Run the simulation
-op(time=Nt, dt=dt)
-
-# The solution is now in u.data
-print(f"Simulation complete: {Nt} time steps")
-print(f"Max amplitude at t={T}: {np.max(np.abs(u.data[0, :])):.6f}")
-```
+{{< include snippets/first_pde_wave1d.qmd >}}
 
 ### Understanding the Code
 

chapters/devito_intro/snippets/absorbing_bc_right_wave.qmd

@@ -0,0 +1,7 @@
+::: {.callout-note title="Snippet (tested)"}
+This code is included from `src/book_snippets/absorbing_bc_right_wave.py` and exercised by `pytest` (see `tests/test_book_snippets.py`).
+:::
+
+```python
+{{< include ../../src/book_snippets/absorbing_bc_right_wave.py >}}
+```