Skip to content

Commit 1612506

Browse files
committed
Add Devito-focused appendix content (Phase 6)
- Add "Devito and Truncation Errors" section to truncation appendix covering space_order parameter, stencil viewing, and accuracy tradeoffs - Add "Software Engineering with Devito" section to softeng appendix covering project structure, pytest fixtures, convergence testing, and performance profiling - Fix Strang splitting typo in tests (strange -> strang) - Update devito-plan.md to mark Phase 6 complete
1 parent fea7a24 commit 1612506

4 files changed

Lines changed: 457 additions & 10 deletions

File tree

chapters/appendices/softeng2/softeng2.qmd

Lines changed: 277 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2065,6 +2065,283 @@ classes are used. The arrays coming from Python, and looking like
20652065
plain C/C++ arrays, can be efficiently wrapped in more user-friendly
20662066
C++ array classes in the C++ code, if desired.
20672067
2068+
## Software Engineering with Devito {#sec-softeng2-devito}
2069+
2070+
The previous sections described traditional approaches to migrating
2071+
Python loops to compiled languages. Devito provides an alternative
2072+
paradigm: write the mathematics symbolically in Python, and let the
2073+
framework generate optimized C code automatically.
2074+
2075+
### The Devito Approach
2076+
2077+
Instead of manually writing C, Cython, or Fortran code, Devito:
2078+
2079+
1. Accepts symbolic PDE specifications in Python
2080+
2. Automatically generates optimized C/C++ code
2081+
3. Compiles and caches the generated code
2082+
4. Provides OpenMP parallelization and optional GPU support
2083+
2084+
This eliminates the need for manual low-level coding while achieving
2085+
competitive performance with hand-tuned implementations.
2086+
2087+
### Project Structure for Devito Applications
2088+
2089+
A well-organized Devito project follows standard Python package conventions:
2090+
2091+
```
2092+
my_pde_solver/
2093+
├── src/
2094+
│ ├── __init__.py
2095+
│ ├── solvers/
2096+
│ │ ├── __init__.py
2097+
│ │ ├── wave.py # Wave equation solvers
2098+
│ │ └── diffusion.py # Diffusion equation solvers
2099+
│ └── utils/
2100+
│ ├── __init__.py
2101+
│ ├── visualization.py # Plotting utilities
2102+
│ └── convergence.py # Convergence testing
2103+
├── tests/
2104+
│ ├── conftest.py # Pytest fixtures
2105+
│ ├── test_wave.py
2106+
│ └── test_diffusion.py
2107+
├── examples/
2108+
│ └── run_simulation.py
2109+
├── pyproject.toml
2110+
└── README.md
2111+
```
2112+
2113+
### Pytest Fixtures for Devito Testing
2114+
2115+
Devito's Grid and Function objects can be reused across tests using
2116+
pytest fixtures:
2117+
2118+
```python
2119+
# tests/conftest.py
2120+
import pytest
2121+
import numpy as np
2122+
from devito import Grid, TimeFunction, Function
2123+
2124+
2125+
@pytest.fixture
2126+
def grid_1d():
2127+
"""Create a standard 1D grid for testing."""
2128+
return Grid(shape=(101,), extent=(1.0,))
2129+
2130+
2131+
@pytest.fixture
2132+
def grid_2d():
2133+
"""Create a standard 2D grid for testing."""
2134+
return Grid(shape=(101, 101), extent=(1.0, 1.0))
2135+
2136+
2137+
@pytest.fixture
2138+
def wave_field(grid_2d):
2139+
"""Create a TimeFunction for wave equation testing."""
2140+
return TimeFunction(name='u', grid=grid_2d,
2141+
time_order=2, space_order=4)
2142+
2143+
2144+
@pytest.fixture
2145+
def velocity_model(grid_2d):
2146+
"""Create a velocity model with constant value."""
2147+
c = Function(name='c', grid=grid_2d)
2148+
c.data[:] = 1500.0 # m/s
2149+
return c
2150+
```
2151+
2152+
Usage in tests:
2153+
2154+
```python
2155+
# tests/test_wave.py
2156+
def test_wave_propagation(grid_2d, wave_field, velocity_model):
2157+
"""Test that wave equation solver runs without error."""
2158+
from src.solvers.wave import solve_acoustic_wave
2159+
2160+
result = solve_acoustic_wave(
2161+
grid=grid_2d,
2162+
u=wave_field,
2163+
c=velocity_model,
2164+
T=0.1,
2165+
)
2166+
2167+
assert result is not None
2168+
assert not np.isnan(result.u.data).any()
2169+
```
2170+
2171+
### Convergence Testing Pattern
2172+
2173+
Verifying numerical schemes against manufactured solutions is essential.
2174+
Here's a reusable pattern:
2175+
2176+
```python
2177+
def convergence_test(solver_func, exact_solution, grid_sizes, **solver_kwargs):
2178+
"""
2179+
Run a convergence test for a Devito solver.
2180+
2181+
Parameters
2182+
----------
2183+
solver_func : callable
2184+
Solver function that returns a result with .u attribute
2185+
exact_solution : callable
2186+
Function(x, t) returning exact solution
2187+
grid_sizes : list
2188+
List of N values to test
2189+
**solver_kwargs : dict
2190+
Additional arguments passed to solver
2191+
2192+
Returns
2193+
-------
2194+
rates : list
2195+
Computed convergence rates between successive refinements
2196+
"""
2197+
import numpy as np
2198+
2199+
errors = []
2200+
dx_values = []
2201+
2202+
for N in grid_sizes:
2203+
result = solver_func(Nx=N, **solver_kwargs)
2204+
2205+
# Compute error at final time
2206+
x = np.linspace(0, result.L, N + 1)
2207+
u_exact = exact_solution(x, result.t)
2208+
error = np.max(np.abs(result.u - u_exact))
2209+
2210+
errors.append(error)
2211+
dx_values.append(result.L / N)
2212+
2213+
# Compute convergence rates
2214+
rates = []
2215+
for i in range(len(errors) - 1):
2216+
rate = np.log(errors[i] / errors[i + 1]) / np.log(2)
2217+
rates.append(rate)
2218+
2219+
return rates
2220+
2221+
2222+
# Usage in test
2223+
def test_diffusion_convergence():
2224+
from src.solvers.diffusion import solve_diffusion
2225+
2226+
rates = convergence_test(
2227+
solver_func=solve_diffusion,
2228+
exact_solution=lambda x, t: np.exp(-np.pi**2 * t) * np.sin(np.pi * x),
2229+
grid_sizes=[20, 40, 80, 160],
2230+
T=0.01,
2231+
a=1.0,
2232+
)
2233+
2234+
# Expect second-order convergence
2235+
assert all(r > 1.9 for r in rates), f"Convergence rates {rates} < 2"
2236+
```
2237+
2238+
### Performance Profiling with Devito
2239+
2240+
Devito provides built-in profiling through environment variables:
2241+
2242+
```python
2243+
import os
2244+
2245+
# Enable performance logging
2246+
os.environ['DEVITO_LOGGING'] = 'PERF'
2247+
2248+
# Run your simulation
2249+
from src.solvers.wave import solve_acoustic_wave
2250+
result = solve_acoustic_wave(...)
2251+
2252+
# Output will include timing information for each operator
2253+
```
2254+
2255+
For more detailed analysis:
2256+
2257+
```python
2258+
from devito import configuration
2259+
2260+
# Enable detailed profiling
2261+
configuration['profiling'] = 'advanced'
2262+
2263+
# Create and run operator
2264+
op = Operator([update_eq])
2265+
summary = op.apply(time=nt, dt=dt)
2266+
2267+
# Access timing information
2268+
print(f"Total time: {summary.globals['fdlike'].time:.3f} s")
2269+
print(f"GFLOPS: {summary.globals['fdlike'].gflopss:.2f}")
2270+
```
2271+
2272+
### Caching and Compilation
2273+
2274+
Devito caches compiled operators to avoid recompilation:
2275+
2276+
```python
2277+
from devito import configuration
2278+
2279+
# View cache location
2280+
print(configuration['cachedir'])
2281+
2282+
# Force recompilation (useful during development)
2283+
configuration['jit-backdoor'] = True
2284+
```
2285+
2286+
For production runs, ensure the cache is preserved between runs to
2287+
avoid recompilation overhead.
2288+
2289+
### Result Classes for Solver Output
2290+
2291+
Using dataclasses provides clean interfaces for solver results:
2292+
2293+
```python
2294+
from dataclasses import dataclass, field
2295+
import numpy as np
2296+
2297+
2298+
@dataclass
2299+
class SolverResult:
2300+
"""Container for solver output."""
2301+
u: np.ndarray # Solution at final time
2302+
x: np.ndarray # Spatial grid
2303+
t: float # Final time
2304+
L: float # Domain length
2305+
dx: float # Grid spacing
2306+
dt: float # Time step used
2307+
nsteps: int # Number of time steps
2308+
u_history: list = field(default_factory=list) # Optional history
2309+
t_history: list = field(default_factory=list) # Time points
2310+
2311+
2312+
def solve_with_result(...):
2313+
"""Solver that returns a SolverResult."""
2314+
# ... solver code ...
2315+
2316+
return SolverResult(
2317+
u=np.array(u.data[0, :]),
2318+
x=x_values,
2319+
t=t_final,
2320+
L=L,
2321+
dx=dx,
2322+
dt=dt,
2323+
nsteps=nt,
2324+
)
2325+
```
2326+
2327+
### Comparison with Manual Optimization
2328+
2329+
The following table compares Devito with manual optimization approaches:
2330+
2331+
| Approach | Development Time | Performance | Portability | Maintainability |
2332+
|----------|-----------------|-------------|-------------|-----------------|
2333+
| Pure Python | Low | Poor | High | High |
2334+
| NumPy vectorized | Medium | Medium | High | Medium |
2335+
| Cython | High | Good | Medium | Low |
2336+
| Fortran/f2py | High | Excellent | Low | Low |
2337+
| C/C++ | Very High | Excellent | Low | Low |
2338+
| **Devito** | Low | Excellent | High | High |
2339+
2340+
Devito achieves performance comparable to hand-tuned code while
2341+
maintaining the simplicity and portability of Python. This makes it
2342+
an excellent choice for scientific computing projects where both
2343+
productivity and performance matter.
2344+
20682345
## Exercise: Explore computational efficiency of numpy.sum versus built-in sum {#sec-softeng2-exer-sum}
20692346
20702347
Using the task of computing the sum of the first `n` integers, we want to

0 commit comments

Comments
 (0)