Refactor Wave and Diffusion chapters to Devito-only

- Fix broken LaTeX subscripts (**{ → _{) in wave and diffu chapters
- Fix Python syntax errors (missing * operators) in wave1D_fd1.qmd
- Fix broken cross-references (@sec- → @eq-) in wave_app_exer.qmd
- Remove NumPy implementation files from chapters and src directories
- Update all file references to point to Devito implementations
- Add Neumann BC, verification, and convergence sections to wave1D_devito.qmd
- Update exercise solutions to use Devito solvers

Deleted files:
- chapters/wave/wave1D_prog.qmd, wave2D_prog.qmd
- chapters/wave/exer-wave/ (exercise solutions)
- chapters/diffu/exer-diffu/ (exercise solutions)
- src/wave/wave1D/, wave2D/, wave2D_u0/ (NumPy solvers)
- src/diffu/diffu*_u0.py, diffu1D_vc.py, etc. (NumPy solvers)

All 247 tests pass, PDF builds without errors.

Author: ggorman

chapters/appendices/softeng2/softeng2.qmd

@@ -42,22 +42,9 @@ and @sec-wave-pde2-var-c.
 ## A solver function
 
 The general initial-boundary value problem
-solved by finite difference methods can be implemented as shown in
-the following `solver` function (taken from the
-file [`wave1D_dn_vc.py`](https://github.com/devitocodes/devito_book/tree/main/src/wave/wave1D/wave1D_dn_vc.py)).
-This function builds on
-simpler versions described in
-@sec-wave-pde1-impl, @sec-wave-pde1-impl-vec,
-@sec-wave-pde2-Neumann, and @sec-wave-pde2-var-c.
-There are several quite advanced
-constructs that will be commented upon later.
-The code is lengthy, but that is because we provide a lot of
-flexibility with respect to input arguments,
-boundary conditions, and optimization
-(scalar versus vectorized loops).
-
-```{.python include="../../../src/wave/wave1D/wave1D_dn_vc.py" start-after="def solver" end-before="def test_quadratic"}
-```
+solved by finite difference methods. For modern implementations using
+Devito, see @sec-wave-devito. This section covers software engineering
+principles that apply broadly to scientific computing.
 
 ## Storing simulation data in files {#sec-softeng2-wave1D-filestorage}
 
@@ -442,9 +429,8 @@ The returned sequence `r` should converge to 2 since the error
 analysis in @sec-wave-pde1-analysis predicts various error measures to behave
 like $\Oof{\Delta t^2} + \Oof{\Delta x^2}$. We can easily run
 the case with standing waves and the analytical solution
-$u(x,t) = \cos(\frac{2\pi}{L}t)\sin(\frac{2\pi}{L}x)$. The call will
-be very similar to the one provided in the `test_convrate_sincos` function
-in @sec-wave-pde1-impl-verify-rate, see the file `src/wave/wave1D/wave1D_dn_vc.py` for details.
+$u(x,t) = \cos(\frac{2\pi}{L}t)\sin(\frac{2\pi}{L}x)$. See @sec-wave-devito-convergence
+for details on convergence rate testing with Devito.
 
 Many who know about class programming prefer to organize their software
 in terms of classes. This gives a richer application programming interface
@@ -1221,9 +1207,7 @@ The `wave2D_u0.py` file contains a `solver` function, which calls an
 in time.  The function `advance_scalar` applies standard Python loops
 to implement the scheme, while `advance_vectorized` performs
 corresponding vectorized arithmetics with array slices. The statements
-of this solver are explained in @sec-wave-2D3D-impl, in
-particular @sec-wave2D3D-impl-scalar and
-@sec-wave2D3D-impl-vectorized.
+of this solver are explained in @sec-wave-2D3D-models.
 
 Although vectorization can bring down the CPU time dramatically
 compared with scalar code, there is still some factor 5-10 to win in
@@ -1506,9 +1490,8 @@ to a single index.
 
 We write a Fortran subroutine `advance` in a file
 [`wave2D_u0_loop_f77.f`](https://github.com/devitocodes/devito_book/tree/main/src/softeng2/wave2D_u0_loop_f77.f)
-for implementing the updating formula
-(@eq-wave-2D3D-impl1-2Du0-ueq-discrete) and setting the solution to zero
-at the boundaries:
+for implementing the updating formula (@eq-wave-2D3D-models-unp1) and
+setting the solution to zero at the boundaries:
 
 ```fortran
 subroutine advance(u, u_1, u_2, f, Cx2, Cy2, dt2, Nx, Ny)

chapters/diffu/diffu_devito_exercises.qmd

@@ -521,63 +521,51 @@ The 2D solver should also achieve second-order spatial convergence
 when the Fourier number is held fixed.
 :::
 
-### Exercise 10: Comparison with Legacy Code {#exer-diffu-legacy}
+### Exercise 10: Performance Scaling {#exer-diffu-scaling}
 
-Compare the Devito solver with the legacy NumPy implementation.
+Investigate how Devito's performance scales with problem size.
 
-a) Run both solvers with the same parameters.
-b) Verify they produce the same results.
-c) Compare execution times.
+a) Run the 1D solver with increasing grid sizes (Nx = 100, 500, 1000, 5000).
+b) Measure and plot the execution time vs grid size.
+c) Determine if the scaling is linear in Nx.
 
 ::: {.callout-note collapse="true" title="Solution"}
 ```python
 from src.diffu import solve_diffusion_1d
-from src.diffu.diffu1D_u0 import solver_FE_simple
 import numpy as np
 import time
+import matplotlib.pyplot as plt
 
 # Parameters
 L = 1.0
 a = 1.0
-Nx = 200
 F = 0.5
-T = 0.1
+T = 0.01  # Short time for timing tests
 
-dx = L / Nx
-dt = F * dx**2 / a
-
-# Devito solver
-t0 = time.perf_counter()
-result_devito = solve_diffusion_1d(
-    L=L, a=a, Nx=Nx, T=T, F=F,
-    I=lambda x: np.sin(np.pi * x),
-)
-t_devito = time.perf_counter() - t0
-
-# Legacy NumPy solver
-t0 = time.perf_counter()
-u_legacy, x_legacy, t_legacy, cpu_legacy = solver_FE_simple(
-    I=lambda x: np.sin(np.pi * x),
-    a=a,
-    f=lambda x, t: 0,
-    L=L,
-    dt=dt,
-    F=F,
-    T=T,
-)
-t_numpy = time.perf_counter() - t0
+grid_sizes = [100, 500, 1000, 5000]
+times = []
 
-# Compare results
-diff = np.max(np.abs(result_devito.u - u_legacy))
-print(f"Maximum difference: {diff:.2e}")
-print(f"Devito time: {t_devito:.4f} s")
-print(f"NumPy time: {t_numpy:.4f} s")
+for Nx in grid_sizes:
+    t0 = time.perf_counter()
+    result = solve_diffusion_1d(
+        L=L, a=a, Nx=Nx, T=T, F=F,
+        I=lambda x: np.sin(np.pi * x),
+    )
+    times.append(time.perf_counter() - t0)
+    print(f"Nx={Nx}: {times[-1]:.4f} s")
 
-# Note: For small problems, NumPy may be faster due to compilation
-# overhead. For large problems, Devito's optimized C code wins.
+# Plot scaling
+plt.figure(figsize=(8, 6))
+plt.loglog(grid_sizes, times, 'bo-', label='Measured')
+plt.loglog(grid_sizes, times[0]*(np.array(grid_sizes)/grid_sizes[0]),
+           'r--', label='O(N)')
+plt.xlabel('Grid size (Nx)')
+plt.ylabel('Time (s)')
+plt.legend()
+plt.title('Devito 1D Diffusion Solver Scaling')
+plt.grid(True)
 ```
 
-For large grids, Devito's automatically generated and optimized C code
-typically outperforms pure Python/NumPy implementations. The advantage
-grows with problem size.
+The Devito solver typically shows linear scaling in Nx for 1D problems,
+as expected for an explicit scheme where each time step is O(Nx).
 :::