Devito phase 1

.typos.toml

@@ -6,3 +6,5 @@ ue = "ue"
 # Strang splitting (named after mathematician Gilbert Strang)
 strang = "strang"
 Strang = "Strang"
+# Variable name for "p at iteration n" in Jacobi iteration
+pn = "pn"

Makefile

@@ -0,0 +1,70 @@
+# Makefile for Finite Difference Computing with PDEs book
+
+.PHONY: pdf html all preview clean test test-devito test-no-devito lint format check help
+
+# Default target
+all: pdf
+
+# Build targets
+pdf:
+	quarto render --to pdf
+
+html:
+	quarto render --to html
+
+# Build both PDF and HTML
+book:
+	quarto render
+
+# Live preview with hot reload
+preview:
+	quarto preview
+
+# Clean build artifacts
+clean:
+	rm -rf _book/
+	rm -rf .quarto/
+	find . -name "*.aux" -delete
+	find . -name "*.log" -delete
+	find . -name "*.out" -delete
+
+# Test targets
+test:
+	pytest tests/ -v
+
+test-devito:
+	pytest tests/ -v -m devito
+
+test-no-devito:
+	pytest tests/ -v -m "not devito"
+
+test-phase1:
+	pytest tests/test_elliptic_devito.py tests/test_burgers_devito.py tests/test_swe_devito.py -v
+
+# Linting and formatting
+lint:
+	ruff check src/
+
+format:
+	ruff check --fix src/
+	isort src/
+
+check:
+	pre-commit run --all-files
+
+# Help
+help:
+	@echo "Available targets:"
+	@echo "  pdf          - Build PDF (default)"
+	@echo "  html         - Build HTML"
+	@echo "  book         - Build all formats (PDF + HTML)"
+	@echo "  preview      - Live preview with hot reload"
+	@echo "  clean        - Remove build artifacts"
+	@echo "  test         - Run all tests"
+	@echo "  test-devito  - Run only Devito tests"
+	@echo "  test-no-devito - Run tests without Devito"
+	@echo "  test-phase1  - Run Phase 1 tests (elliptic, burgers, swe)"
+	@echo "  lint         - Check code with ruff"
+	@echo "  format       - Auto-format code with ruff and isort"
+	@echo "  check        - Run all pre-commit hooks"
+	@echo "  help         - Show this help message"

_quarto.yml

@@ -6,8 +6,8 @@ book:
   title: "Finite Difference Computing with PDEs"
   subtitle: "A Devito Approach"
   author:
-    - name: Hans Petter Langtangen
-    - name: Svein Linge
+    - name: Gerard J. Gorman
+      affiliation: Imperial College London
   date: today
   chapters:
     - index.qmd
@@ -19,6 +19,8 @@ book:
         - chapters/diffu/index.qmd
         - chapters/advec/index.qmd
         - chapters/nonlin/index.qmd
+        - chapters/elliptic/index.qmd
+        - chapters/systems/index.qmd
     - part: "Appendices"
       chapters:
         - chapters/appendices/formulas/index.qmd
@@ -169,6 +171,8 @@ src_diffu: "https://github.com/devitocodes/devito_book/tree/devito/src/diffu"
 src_nonlin: "https://github.com/devitocodes/devito_book/tree/devito/src/nonlin"
 src_trunc: "https://github.com/devitocodes/devito_book/tree/devito/src/trunc"
 src_advec: "https://github.com/devitocodes/devito_book/tree/devito/src/advec"
+src_elliptic: "https://github.com/devitocodes/devito_book/tree/devito/src/elliptic"
+src_systems: "https://github.com/devitocodes/devito_book/tree/devito/src/systems"
 src_formulas: "https://github.com/devitocodes/devito_book/tree/devito/src/formulas"
 src_softeng2: "https://github.com/devitocodes/devito_book/tree/devito/src/softeng2"
 

chapters/advec/advec.qmd

@@ -72,7 +72,7 @@ u(i\Delta x, (n+1)\Delta t) &= I(i\Delta x - v(n+1)\Delta t) \nonumber \\
 provided $v = \Delta x/\Delta t$. So, whenever we see a scheme that
 collapses to
 $$
-u^{n+1}**i = u**{i-1}^n,
+u^{n+1}_i = u_{i-1}^n,
 $$ {#eq-advec-1D-pde1-uprop2}
 for the PDE in question, we have in fact a scheme that reproduces the
 analytical solution, and many of the schemes to be presented possess
@@ -771,11 +771,12 @@ $$
 $$
 which results in the updating formula
 $$
-u^{n+1}_i = u^{n-1}**i - C(u**{i+1}^n-u_{i-1}^n)\tp
+u^{n+1}_i = u^{n-1}_i - C(u_{i+1}^n-u_{i-1}^n)\tp
 $$
 A special scheme is needed to compute $u^1$, but we leave that problem for
-now. Anyway, this special scheme can be found in
-[`advec1D.py`](https://github.com/devitocodes/devito_book/tree/main/src/advec/advec1D.py).
+now. The Devito implementation handles this automatically; see
+[`advec1D_devito.py`](https://github.com/devitocodes/devito_book/tree/main/src/advec/advec1D_devito.py)
+and @sec-advec-devito for details.
 
 ### Implementation
 We now need to work with three time levels and must modify our solver a bit:
@@ -888,7 +889,7 @@ $$
 $$ {#eq-advec-1D-upwind}
 written out as
 $$
-u^{n+1}_i = u^n_i - C(u^{n}**{i}-u^{n}**{i-1}),
+u^{n+1}_i = u^n_i - C(u^{n}_{i}-u^{n}_{i-1}),
 $$
 gives a generally popular and robust scheme that is stable if $C\leq 1$.
 As with the Leapfrog scheme, it becomes exact if $C=1$, exactly as shown in
@@ -929,7 +930,7 @@ $$
 $$
 by a forward difference in time and centered differences in space,
 $$
-D^+**t u + vD**{2x} u = \nu D_xD_x u]^n_i,
+D^+_t u + vD_{2x} u = \nu D_xD_x u]^n_i,
 $$
 actually gives the upwind scheme (@eq-advec-1D-upwind) if
 $\nu = v\Delta x/2$. That is, solving the PDE $u_t + vu_x=0$
@@ -1170,30 +1171,31 @@ def run(scheme='UP', case='gaussian', C=1, dt=0.01):
         os.system(cmd)
     print 'Integral of u:', integral.max(), integral.min()
 ```
-The complete code is found in the file
-[`advec1D.py`](https://github.com/devitocodes/devito_book/tree/main/src/advec/advec1D.py).
+The Devito implementation is found in
+[`advec1D_devito.py`](https://github.com/devitocodes/devito_book/tree/main/src/advec/advec1D_devito.py).
+See @sec-advec-devito for the complete implementation.
 
 ## A Crank-Nicolson discretization in time and centered differences in space {#sec-advec-1D-CN}
 
 Another obvious candidate for time discretization is the Crank-Nicolson
 method combined with centered differences in space:
 $$
-[D_t u]^n_i + v\half([D_{2x} u]^{n+1}**i + [D**{2x} u]^{n}_i) = 0\tp
+[D_t u]^n_i + v\half([D_{2x} u]^{n+1}_i + [D_{2x} u]^{n}_i) = 0\tp
 $$
 It can be nice to include the Backward Euler scheme too, via the
 $\theta$-rule,
 $$
-[D_t u]^n_i + v\theta [D_{2x} u]^{n+1}**i + v(1-\theta)[D**{2x} u]^{n}_i = 0\tp
+[D_t u]^n_i + v\theta [D_{2x} u]^{n+1}_i + v(1-\theta)[D_{2x} u]^{n}_i = 0\tp
 $$
 When $\theta$ is different from zero, this gives rise to an *implicit* scheme,
 $$
-u^{n+1}**i + \frac{\theta}{2} C (u^{n+1}**{i+1} - u^{n+1}_{i-1})
-= u^n_i - \frac{1-\theta}{2} C (u^{n}**{i+1} - u^{n}**{i-1})
+u^{n+1}_i + \frac{\theta}{2} C (u^{n+1}_{i+1} - u^{n+1}_{i-1})
+= u^n_i - \frac{1-\theta}{2} C (u^{n}_{i+1} - u^{n}_{i-1})
 $$
 for $i=1,\ldots,N_x-1$. At the boundaries we set $u=0$ and simulate just to
 the point of time when the signal hits the boundary (and gets reflected).
 $$
-u^{n+1}**0 = u^{n+1}**{N_x} = 0\tp
+u^{n+1}_0 = u^{n+1}_{N_x} = 0\tp
 $$
 The elements on the diagonal in the matrix become:
 $$
@@ -1208,7 +1210,7 @@ And finally, the right-hand side becomes
 
 \begin{align*}
 b_0 &= u^n_{N_x}\\
-b_i &= u^n_i - \frac{1-\theta}{2} C (u^{n}**{i+1} - u^{n}**{i-1}),\quad i=1,\ldots,N_x-1\\
+b_i &= u^n_i - \frac{1-\theta}{2} C (u^{n}_{i+1} - u^{n}_{i-1}),\quad i=1,\ldots,N_x-1\\
 b_{N_x} &= u^n_0
 \end{align*}
 
@@ -1273,7 +1275,7 @@ u^{n+1}_i = u^n_i -v \Delta t [D_{2x} u]^n_i
 $$
 or written out,
 $$
-u^{n+1}_i = u^n_i - \frac{1}{2} C (u^{n}**{i+1} - u^{n}**{i-1})
+u^{n+1}_i = u^n_i - \frac{1}{2} C (u^{n}_{i+1} - u^{n}_{i-1})
 + \frac{1}{2} C^2 (u^{n}_{i+1}-2u^n_i+u^n_{i-1})\tp
 $$
 This is the explicit Lax-Wendroff scheme.
@@ -1576,15 +1578,15 @@ is the dominating term, collect its information in the flow direction, i.e.,
 upstream or upwind of the point in question. So, instead of using a
 centered difference
 $$
-\frac{du}{dx}**i\approx \frac{u**{i+1}-u_{i-1}}{2\Delta x},
+\frac{du}{dx}_i\approx \frac{u_{i+1}-u_{i-1}}{2\Delta x},
 $$
 we use the one-sided *upwind* difference
 $$
-\frac{du}{dx}**i\approx \frac{u**{i}-u_{i-1}}{\Delta x},
+\frac{du}{dx}_i\approx \frac{u_{i}-u_{i-1}}{\Delta x},
 $$
 in case $v>0$. For $v<0$ we set
 $$
-\frac{du}{dx}**i\approx \frac{u**{i+1}-u_{i}}{\Delta x},
+\frac{du}{dx}_i\approx \frac{u_{i+1}-u_{i}}{\Delta x},
 $$
 On compact operator notation form, our upwind scheme can be expressed
 as

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/appendices/trunc/trunc.qmd

@@ -1644,12 +1644,12 @@ $[D_xD_x (D_tD_t u)]^n_i$ gives
 
 \begin{align*}
 \frac{1}{\Delta t^2}\biggl(
-&\frac{u^{n+1}**{i+1} - 2u^{n}**{i+1} + u^{n-1}_{i+1}}{\Delta x^2} -2\\
-&\frac{u^{n+1}**{i} - 2u^{n}**{i} + u^{n-1}_{i}}{\Delta x^2} +
-&\frac{u^{n+1}**{i-1} - 2u^{n}**{i-1} + u^{n-1}_{i-1}}{\Delta x^2}
+&\frac{u^{n+1}_{i+1} - 2u^{n}_{i+1} + u^{n-1}_{i+1}}{\Delta x^2} -2\\
+&\frac{u^{n+1}_{i} - 2u^{n}_{i} + u^{n-1}_{i}}{\Delta x^2} +
+&\frac{u^{n+1}_{i-1} - 2u^{n}_{i-1} + u^{n-1}_{i-1}}{\Delta x^2}
 \biggr)
 \end{align*}
-Now the unknown values $u^{n+1}**{i+1}$, $u^{n+1}**{i}$,
+Now the unknown values $u^{n+1}_{i+1}$, $u^{n+1}_{i}$,
 and $u^{n+1}_{i-1}$ are *coupled*, and we must solve a tridiagonal
 system to find them. This is in principle straightforward, but it
 results in an implicit finite difference scheme, while we had

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).
 :::