Add ABC methods (PML, Higdon, HABC), EM chapter, and fix Devito solvers

Major additions:
- Implement split-field PML (Grote-Sim), second-order Higdon ABC (P=2),
  and Hybrid ABC with weighted absorption layer in abc_methods.py
- Add electromagnetics chapter (Maxwell equations, Yee scheme,
  verification, GPR and waveguide applications)
- Add tested book snippets for all ABC methods

Fixes and improvements:
- Fix sigma_max default to theory-derived 3c/W instead of hardcoded 50
- Replace linear damping ramp with polynomial d^3 in wave1D_features
- Add dtype parameter to all Devito solvers for FP64 support
- Standardize LaTeX notation (\hbox→\text, font-size removal)
- Tighten test tolerances to match O(dx^2) theory
- Add elliptic spatial convergence test and Pint unit checks
- Update .gitignore for TeX artifacts and review documents

All 579 tests pass.

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>

Author: ggorman

.gitignore

@@ -53,6 +53,24 @@ doc/pub/
 /_book/
 **/*.quarto_ipynb
 
+# Root-level TeX artifacts (Quarto/Pandoc intermediates)
+/*.tex
+
+# Quarto HTML render artifacts
+**/*_files/
+
+# Review and planning documents
+book-review.md
+codex-review.md
+codex-reviewer.md
+review.md
+review-2.md
+review-3.md
+review-4.md
+ianmgdev_review.txt
+*.docx
+Proposed Extension Roadmap*.md
+
 # Coverage
 .coverage
 coverage.xml

.typos.toml

@@ -8,3 +8,5 @@ strang = "strang"
 Strang = "Strang"
 # Variable name for "p at iteration n" in Jacobi iteration
 pn = "pn"
+# Journal abbreviation: "J. Comput. Phys."
+Comput = "Comput"

Makefile

@@ -3,7 +3,7 @@
 .PHONY: pdf html all preview clean test test-devito test-no-devito lint format check help
 
 # Default target
-all: pdf
+all: pdf html
 
 # Build targets
 pdf:

chapters/advec/advec.qmd

@@ -456,7 +456,7 @@ u(x,0) = Ae^{-\half\left(\frac{x-L/10}{\sigma}\right)^2},
 $$
 $$
 u(x,0) = A\cos\left(\frac{5\pi}{L}\left( x - \frac{L}{10}\right)\right),\quad
-x < \frac{L}{5} \hbox{ else } 0\tp
+x < \frac{L}{5} \text{ else } 0\tp
 $$ {#eq-advec-1D-case_gaussian}
 The parameter $A$ is the maximum value of the initial condition.
 
@@ -723,16 +723,16 @@ A general solution may be viewed as a collection of long and
 short waves with different amplitudes. Algebraically, the work
 simplifies if we introduce the complex Fourier component
 $$
-u(x,t)=\Aex e^{ikx},
+u(x,t)=A_{\text{e}} e^{ikx},
 $$
 with
 $$
-\Aex=Be^{-ikv\Delta t} = Be^{-iCk\Delta x}\tp
+A_{\text{e}}=Be^{-ikv\Delta t} = Be^{-iCk\Delta x}\tp
 $$
-Note that $|\Aex| \leq 1$.
+Note that $|A_{\text{e}}| \leq 1$.
 
 It turns out that many schemes also allow a Fourier wave component as
-solution, and we can use the numerically computed values of $\Aex$
+solution, and we can use the numerically computed values of $A_{\text{e}}$
 (denoted $A$) to learn about the
 quality of the scheme. Hence, to analyze the difference scheme we have just
 implemented, we look at how it treats the Fourier component
@@ -862,7 +862,7 @@ A = -iC\sin p \pm \sqrt{1-C^2\sin^2 p},
 $$
 and is to be compared to the exact amplification factor
 $$
-\Aex = e^{-ikv\Delta t} = e^{-ikC\Delta x} = e^{-iCp}\tp
+A_{\text{e}} = e^{-ikv\Delta t} = e^{-ikC\Delta x} = e^{-iCp}\tp
 $$
 Section @sec-advec-1D-disprel compares numerical amplification factors
 of many schemes with the exact expression.
@@ -995,7 +995,7 @@ is constant:
 \end{align*}
 as long as $u(0)=u(L)=0$. We can therefore use the property
 $$
-\int_0^L u(x,t)dx = \hbox{const}
+\int_0^L u(x,t)dx = \text{const}
 $$
 as a partial verification during the simulation. Now, any numerical method
 with $C\neq 1$ will deviate from the constant, expected value, so
@@ -1469,7 +1469,7 @@ reason why this model problem has been so successful in designing and
 investigating numerical methods for mixed convection/advection and
 diffusion.  The exact solution reads
 $$
-\uex (x) = \frac{e^{x/\epsilon} - 1}{e^{1/\epsilon} - 1}\tp
+u_{\text{e}} (x) = \frac{e^{x/\epsilon} - 1}{e^{1/\epsilon} - 1}\tp
 $$
 The forthcoming plots illustrate this function for various values of
 $\epsilon$.

chapters/advec/advec_devito_exercises.qmd

@@ -170,7 +170,7 @@ plt.loglog(sizes_lf, err_lf, 'g-^', label=f'Lax-Friedrichs (rate={rate_lf:.2f})'
 # Reference slopes
 h = np.array(sizes_up)
 plt.loglog(h, err_up[0]*(h[0]/h), 'k--', alpha=0.5, label='O(h)')
-plt.loglog(h, err_lw[0]*(h[0]/h)**2, 'k:', alpha=0.5, label='O(h²)')
+plt.loglog(h, err_lw[0]*(h[0]/h)**2, 'k:', alpha=0.5, label='O(h^2)')
 
 plt.xlabel('Grid points')
 plt.ylabel('L2 Error')

chapters/appendices/formulas/formulas.qmd

@@ -57,69 +57,69 @@ derivative at the point $t_{n+\theta}=\theta t_{n+1} + (1-\theta) t_{n}$.
 
 $$
 \begin{split}
-\uex'(t_n) &=
-[D_t\uex]^n + R^n = \frac{\uex^{n+\half} - \uex^{n-\half}}{\Delta t} +R^n,\\
-R^n &= -\frac{1}{24}\uex'''(t_n)\Delta t^2 + {\cal O}(\Delta t^4)
+u_{\text{e}}'(t_n) &=
+[D_tu_{\text{e}}]^n + R^n = \frac{u_{\text{e}}^{n+\half} - u_{\text{e}}^{n-\half}}{\Delta t} +R^n,\\
+R^n &= -\frac{1}{24}u_{\text{e}}'''(t_n)\Delta t^2 + \Oof{\Delta t^4}
 \end{split}
 $$
 $$
 \begin{split}
-\uex'(t_n) &=
-[D_{2t}\uex]^n +R^n = \frac{\uex^{n+1} - \uex^{n-1}}{2\Delta t} +
+u_{\text{e}}'(t_n) &=
+[D_{2t}u_{\text{e}}]^n +R^n = \frac{u_{\text{e}}^{n+1} - u_{\text{e}}^{n-1}}{2\Delta t} +
 R^n,\\
-R^n &= -\frac{1}{6}\uex'''(t_n)\Delta t^2 + {\cal O}(\Delta t^4)
+R^n &= -\frac{1}{6}u_{\text{e}}'''(t_n)\Delta t^2 + \Oof{\Delta t^4}
 \end{split}
 $$
 $$
 \begin{split}
-\uex'(t_n) &=
-[D_t^-\uex]^n +R^n = \frac{\uex^{n} - \uex^{n-1}}{\Delta t}
+u_{\text{e}}'(t_n) &=
+[D_t^-u_{\text{e}}]^n +R^n = \frac{u_{\text{e}}^{n} - u_{\text{e}}^{n-1}}{\Delta t}
 +R^n,\\
-R^n &= -\half\uex''(t_n)\Delta t + {\cal O}(\Delta t^2)
+R^n &= -\half u_{\text{e}}''(t_n)\Delta t + \Oof{\Delta t^2}
 \end{split}
 $$
 $$
 \begin{split}
-\uex'(t_n) &=
-[D_t^+\uex]^n +R^n = \frac{\uex^{n+1} - \uex^{n}}{\Delta t}
+u_{\text{e}}'(t_n) &=
+[D_t^+u_{\text{e}}]^n +R^n = \frac{u_{\text{e}}^{n+1} - u_{\text{e}}^{n}}{\Delta t}
 +R^n,\\
-R^n &= \half\uex''(t_n)\Delta t + {\cal O}(\Delta t^2)
+R^n &= \half u_{\text{e}}''(t_n)\Delta t + \Oof{\Delta t^2}
 \end{split}
 $$
 $$
 \begin{split}
-\uex'(t_{n+\theta}) &=
-[\bar D_t\uex]^{n+\theta} +R^{n+\theta} = \frac{\uex^{n+1} - \uex^{n}}{\Delta t}
+u_{\text{e}}'(t_{n+\theta}) &=
+[\bar D_tu_{\text{e}}]^{n+\theta} +R^{n+\theta} = \frac{u_{\text{e}}^{n+1} - u_{\text{e}}^{n}}{\Delta t}
 +R^{n+\theta},\\
-R^{n+\theta} &= -\half(1-2\theta)\uex''(t_{n+\theta})\Delta t +
-\frac{1}{6}((1 - \theta)^3 - \theta^3)\uex'''(t_{n+\theta})\Delta t^2 +
+R^{n+\theta} &= -\half(1-2\theta)u_{\text{e}}''(t_{n+\theta})\Delta t +
+\frac{1}{6}((1 - \theta)^3 - \theta^3)u_{\text{e}}'''(t_{n+\theta})\Delta t^2 +
 \\
-&\quad {\cal O}(\Delta t^3)
+&\quad \Oof{\Delta t^3}
 \end{split}
 $$
 $$
 \begin{split}
-\uex'(t_n) &=
-[D_t^{2-}\uex]^n +R^n = \frac{3\uex^{n} - 4\uex^{n-1} + \uex^{n-2}}{2\Delta t}
+u_{\text{e}}'(t_n) &=
+[D_t^{2-}u_{\text{e}}]^n +R^n = \frac{3u_{\text{e}}^{n} - 4u_{\text{e}}^{n-1} + u_{\text{e}}^{n-2}}{2\Delta t}
 +R^n,\\
-R^n &= \frac{1}{3}\uex'''(t_n)\Delta t^2 + {\cal O}(\Delta t^3)
+R^n &= \frac{1}{3}u_{\text{e}}'''(t_n)\Delta t^2 + \Oof{\Delta t^3}
 \end{split}
 $$
 $$
 \begin{split}
-\uex''(t_n) &=
-[D_tD_t \uex]^n +R^n = \frac{\uex^{n+1} - 2\uex^{n} + \uex^{n-1}}{\Delta t^2}
+u_{\text{e}}''(t_n) &=
+[D_tD_t u_{\text{e}}]^n +R^n = \frac{u_{\text{e}}^{n+1} - 2u_{\text{e}}^{n} + u_{\text{e}}^{n-1}}{\Delta t^2}
 +R^n,\\
-R^n &= -\frac{1}{12}\uex''''(t_n)\Delta t^2 + {\cal O}(\Delta t^4)
+R^n &= -\frac{1}{12}u_{\text{e}}''''(t_n)\Delta t^2 + \Oof{\Delta t^4}
 \end{split}
 $$ {#eq-form-trunc-fd1-center}
 
 $$
 \begin{split}
-\uex(t_{n+\theta}) &= [\overline{\uex}^{t,\theta}]^{n+\theta} +R^{n+\theta}
-= \theta \uex^{n+1} + (1-\theta)\uex^n +R^{n+\theta},\\
-R^{n+\theta} &= -\half\uex''(t_{n+\theta})\Delta t^2\theta (1-\theta) +
-{\cal O}(\Delta t^3)\tp
+u_{\text{e}}(t_{n+\theta}) &= [\overline{u_{\text{e}}}^{t,\theta}]^{n+\theta} +R^{n+\theta}
+= \theta u_{\text{e}}^{n+1} + (1-\theta)u_{\text{e}}^n +R^{n+\theta},\\
+R^{n+\theta} &= -\half u_{\text{e}}''(t_{n+\theta})\Delta t^2\theta (1-\theta) +
+\Oof{\Delta t^3}\tp
 \end{split}
 $$ {#eq-form-trunc-avg-theta}
 

chapters/appendices/softeng2/softeng2.qmd

@@ -16,10 +16,10 @@ $$
 u_t(x,0) = V(t),\quad x\in [0,L]
 $$
 $$
-u(0,t) = U_0(t)\hbox{ or } u_x(0,t)=0,\quad t\in (0,T]
+u(0,t) = U_0(t)\text{ or } u_x(0,t)=0,\quad t\in (0,T]
 $$
 $$
-u(L,t) = U_L(t)\hbox{ or } u_x(L,t)=0,\quad t\in (0,T]
+u(L,t) = U_L(t)\text{ or } u_x(L,t)=0,\quad t\in (0,T]
 $$ {#eq-wave-pde2-software-ueq2}
 We allow variable wave velocity $c^2(x)=q(x)$, and Dirichlet or homogeneous
 Neumann conditions at the boundaries.