# ADI for a Mixed-Derivative (Correlated) 2D Parabolic PDE — Task Specification You are working on a Linux coding task about alternating-direction-implicit finite-difference methods for parabolic PDEs with a mixed second derivative. ## Task directory - `base` ## Visible input files - Problem specification: `base/input/problem_spec.md` - Python runtime manifest: `base/input/runtime_env/pyproject.toml` - Python lockfile: `base/input/runtime_env/uv.lock` - Canonical Python wrapper: `base/software/python.sh` ## Recommended workflow 1. Change into `base`. 2. Read `base/input/problem_spec.md` in full before coding. 3. Use `"base/software/python.sh" your_script.py` or `"base/software/python.sh" -c "..."` to run Python with the staged NumPy/SciPy environment. 4. Implement your solver and save it to `base/output/adi_solver.py`. 5. Run the solver and write every required artifact under `base/output`. ## Required outputs under `base/output` - `adi_solver.py` - `reference_arrays.npz` - `results.json` --- ## Background You will build and validate an **Alternating Direction Implicit (ADI)** solver for a two-dimensional parabolic PDE that contains a **mixed second derivative** `u_xy`. This is the model operator behind correlated-diffusion problems (it arises, for example, whenever a diffusion process has correlated components). The mixed term is what makes the problem hard: a directional ADI sweep cannot make `u_xy` implicit, and a naive treatment either destroys the convergence order or goes unstable. ## The problem Solve $$ u_t = a\,u_{xx} + 2c\,u_{xy} + b\,u_{yy}, \qquad (x, y) \in D = (0,1)^2, \quad t \in [0, T], \quad T = 1, $$ with **diffusion matrix** $$ M = \begin{pmatrix} a & c \\ c & b \end{pmatrix} = \begin{pmatrix} 1.0 & 0.5 \\ 0.5 & 0.8 \end{pmatrix}, $$ which is symmetric positive definite (`a>0`, `b>0`, `ab − c² = 0.55 > 0`), so the PDE is uniformly parabolic. **Boundary and initial data** are taken from the manufactured exact solution $$ u_e(x, y, t) = \sum_{m=1}^{3} A_m \, e^{-\lambda_m t} \sin(\alpha_m x + \beta_m y + \varphi_m), \qquad \lambda_m = a\,\alpha_m^2 + 2c\,\alpha_m\beta_m + b\,\beta_m^2, $$ with parameters | m | A_m | α_m | β_m | φ_m | |---|-----|-----|-----|-----| | 1 | 1.0 | 0.7 | 1.3 | 0.0 | | 2 | 0.7 | 1.4 | -0.6 | 0.3 | | 3 | 0.5 | 1.0 | 0.8 | 0.7 | Each rotated wave solves the PDE exactly (`λ_m ≥ 0` by SPD), so `u_e` does too. Dirichlet data is `g = u_e|_{∂D}` and the initial condition is `u_0 = u_e(·,·,0)`. > **Note.** A single separable product `sin(px)·sin(qy)` is **not** a solution: > the mixed term turns `sin·sin` into `cos·cos`. The cross term genuinely > couples the directions — you cannot ignore it. ## Discretization - Uniform grid `h = 1/m`, points `(x_i, y_j) = (ih, jh)`, `0 ≤ i, j ≤ m`. - Time step tied to the mesh: `k = h`. - Use **second-order central differences** for all three operators, including the standard nine-point mixed stencil `δ_xy u_{i,j} = (u_{i+1,j+1} − u_{i+1,j−1} − u_{i−1,j+1} + u_{i−1,j−1}) / (4h²)`. ## Required scheme properties Your time integrator must be an **ADI-type scheme whose only linear solves are tridiagonal one-dimensional sweeps** (one per spatial direction), and it must: 1. **Handle the mixed derivative so that the global error is second order** in `h` (with `k = h`). Because a directional sweep cannot make `u_xy` implicit, the mixed term has to be incorporated by other means while preserving second-order accuracy — this is the heart of the task. Dropping the term or lagging it naively will not converge to this PDE. 2. Be **stable** on this mesh sequence (no growth / blow-up). Selecting a time-integration scheme — and any parameter it carries — that is *simultaneously* second order *and* stable in the presence of a mixed derivative is part of what is being tested. A scheme (or parameter value) that is fine for a pure `u_xx + u_yy` problem is not automatically second order or stable once the cross term is present; verify both properties yourself via the convergence study. **Any** ADI scheme is acceptable provided it **demonstrates** the two properties above on this mesh sequence, and that demonstration is defined precisely by the published criteria below — nothing else. You are not being asked to reproduce a particular published scheme, and the grader never compares your errors against those of one privileged method. Note the emphasis on *demonstrates*. A scheme that is second order only in the asymptotic limit, whose asymptotic regime lies beyond `h = 0.00625`, will not show order → 2 on the five prescribed meshes and is therefore not acceptable here, however sound it is in theory. The task asks you to exhibit second-order convergence on `h = 0.1·2^-p`, `p = 0..4`; if your scheme carries a free parameter, choose a value that does so. (The criteria are not narrow: several standard treatments of the mixed derivative qualify, each across a wide range of its own parameter. Identifying one is part of the task.) ### Published acceptance criteria So that nothing is hidden, these are the exact numerical criteria your convergence study must satisfy. They contain no reference to any particular scheme: | criterion | requirement | |---|---| | stability | every `e_h` finite and positive | | monotonicity | `e_h` strictly decreasing under refinement | | total reduction | `e(0.1) / e(0.00625) ≥ 100` | | fitted order | least-squares slope of `log2(e_h)` vs `p` in `[1.70, 2.30]` | | pointwise order | both of the two finest `order_h` in `[1.60, 2.60]` | | sanity bounds | `e(0.1) < 1.0` and `e(0.00625) < 1e-2` | | **hidden instance** | `solve()` called on an undisclosed problem must land within `3e-2` of an independently computed reference at `h = 0.05`, and its error must fall by at least `2×` from `h = 0.1` | The two sanity bounds are deliberately loose — they exclude nonsense, not schemes. Every scheme that converges at second order on this mesh sequence clears them by a wide margin. The real requirement is the convergence behaviour, which is exactly what the task asks you to demonstrate. ## What to compute ### Tier 1 — Single-grid validation (`h = 0.1`) Solve to `T = 1` and report the max-norm error `e_h = max_{(x_i,y_j)} |u(x_i,y_j,1) − u_e(x_i,y_j,1)|`. The result must be finite (a stable run). ### Tier 2 — Convergence study (the core check) For `h = 0.1·2^{−p}`, `p = 0,1,2,3,4`, compute `e_h` and the experimental order `order_h = log2(e_{2h} / e_h)`. Build the `(h, e_h, order_h)` table; the order must approach **2**. (A solver that mishandles the mixed term will instead show the error plateauing and the order collapsing toward 0.) ### Tier 3 — Finest-mesh artifacts On the finest mesh `h = 0.00625` (161×161 grid), save the numerical solution, the exact solution, and the error surface `u − u_e` at `T = 1`. ## Required solver interface `adi_solver.py` must be importable **without side effects** (guard any demo code behind `if __name__ == "__main__":`) and must expose exactly this entry point: ```python def solve(h, T=1.0, problem=None): """Advance from t=0 to t=T on a uniform grid of size h, with k = h. problem : None -> the benchmark problem defined above. Otherwise a dict {"a": float, "b": float, "c": float, # SPD diffusion coefficients "u0": callable(X, Y) -> initial field, "g": callable(X, Y, t) -> Dirichlet data (boundary values are used)} Returns ------- (X, Y, U) : three (m+1, m+1) float arrays, m = round(1/h), built with numpy.meshgrid(..., indexing='ij'), U being the numerical solution at time T (Dirichlet boundary values included). """ ``` **Your solver must be general.** The grader calls `solve()` on a problem it chooses — different diffusion coefficients, and an initial state with **no closed-form solution** — and checks that your result converges to an independently computed reference. Write a solver for the PDE, not for the manufactured solution: anything specialised to `u_e` fails this check, and so does any field produced without actually integrating the equation. Reproducing the benchmark's own numbers is necessary but not sufficient. The grader imports your module in an isolated subprocess and calls `solve()` to regenerate the convergence study. It checks that all five regenerated errors match the errors you submitted and that the regenerated finest numerical field matches `U_finest` **pointwise** (`rtol=1e-10`, `atol=1e-12`). Arrays or surfaces that your code cannot reproduce will not be accepted, so do not hand-edit the outputs. ## Required outputs (write to `base/output`) | File | Format | Contents | |------|--------|----------| | `adi_solver.py` | Python | Your solver, exposing `solve(h, T=1.0, problem=None)` as above. | | `reference_arrays.npz` | NPZ | `h_values`, `max_norm_errors`, `eoc`, and finest-mesh `X_finest`, `Y_finest`, `U_finest`, `U_exact_finest`, `error_finest`. | | `results.json` | JSON | `convergence_table` (5 rows of `p`, `h`, `max_norm_error`, `eoc_order`), `observed_order_of_convergence`, `finest_mesh`, `diffusion_matrix`, `final_time`, `time_step_rule`, and `scheme`. Must be consistent with the arrays. | ### Array schema | key | shape | contents | |---|---|---| | `h_values` | `(5,)` | `0.1·2^-p`, `p = 0..4`, descending | | `max_norm_errors` | `(5,)` | `e_h` for each mesh, same order as `h_values` | | `eoc` | `(5,)` | `log2(e_{2h}/e_h)`; entry 0 has no predecessor — use `NaN` | | `X_finest`, `Y_finest` | `(161,161)` | finest-mesh coordinates, `indexing='ij'` | | `U_finest` | `(161,161)` | numerical solution at `T = 1` | | `U_exact_finest` | `(161,161)` | `u_e` at `T = 1` on the same grid | | `error_finest` | `(161,161)` | `U_finest − U_exact_finest` | `max_norm_errors[4]` must equal `max|U_finest − U_exact_finest|`: the last row of the table and the Tier 3 surfaces must come from the same run. ### `results.json` schema The grader checks these keys and value formats exactly, so use them verbatim: ```json { "scheme": "", "diffusion_matrix": {"a": 1.0, "b": 0.8, "c": 0.5}, "final_time": 1.0, "time_step_rule": "k = h", "convergence_table": [ {"p": 0, "h": 0.1, "max_norm_error": , "eoc_order": null}, {"p": 1, "h": 0.05, "max_norm_error": , "eoc_order": }, {"p": 2, "h": 0.025, "max_norm_error": , "eoc_order": }, {"p": 3, "h": 0.0125, "max_norm_error": , "eoc_order": }, {"p": 4, "h": 0.00625, "max_norm_error": , "eoc_order": } ], "observed_order_of_convergence": , "finest_mesh": {"h": 0.00625, "grid_points_per_side": 161, "max_norm_error": } } ``` *(angle brackets are placeholders — this block shows the required **schema**; no reference magnitudes are disclosed in the solver-visible specification)* | field | requirement | |---|---| | `scheme` | non-empty string | | `diffusion_matrix` | an **object** with keys `a`, `b`, `c` — not a nested list | | `final_time` | `1.0` | | `time_step_rule` | must literally contain `k = h` | | `convergence_table` | exactly 5 rows; each carries `p` (0..4), `h`, `max_norm_error`, `eoc_order`; `eoc_order` is `null` for `p = 0` | | `observed_order_of_convergence` | a number consistent with your errors (within 0.25 of the least-squares fitted order) | | `finest_mesh` | object with `h` = 0.00625, `grid_points_per_side` = 161, `max_norm_error` | Every value must agree with `reference_arrays.npz`. If you cannot complete Tier 3, **omit the five finest-mesh arrays entirely** rather than fabricating them. The three Tier 1 / Tier 2 arrays (`h_values`, `max_norm_errors`, `eoc`) are always required; omitting the Tier 3 arrays costs you Tier 3 credit only. ## Rules - **Python + NumPy + SciPy only.** Do **not** use a PDE/FEM/FDM package (FEniCS, FiPy, py-pde, scikit-fdiff, Dedalus, findiff, etc.). Implement the scheme and the tridiagonal solves yourself (`scipy.linalg.solve_banded`, `scipy.linalg.lu_factor`, or a hand-written Thomas algorithm). - Use the diffusion matrix, mesh sequence, and `k = h` exactly as specified. - Only the one-dimensional sweeps may be solved implicitly; the full 2D operator (including the mixed term) must not be assembled and inverted directly as a single sparse system — that defeats the purpose of ADI. **The grader instruments your solver at run time and records both the dimension and the bandwidth of every linear system it solves**, so a monolithic 2-D solve is detected, not merely discouraged. Concretely, on an `m×m` grid every implicit solve must satisfy: - **dimension ≤ `2(m+1)`** — a genuine sweep is `m−1`; the forbidden monolithic system is `(m−1)²` (at `h = 0.1`: 9 versus 81, limit 22); - **half-bandwidth ≤ 1**, i.e. actually tridiagonal. A *dense* solve fails this even when it is small enough to pass the dimension limit. A hand-written Thomas algorithm makes no library call and so records nothing; that is permitted, and the grader reports the audit as inconclusive for human review rather than failing you. - On an evaluator-selected general-data case, the grader also reruns the solver after replacing every intercepted linear-solve result with invalid data. The returned field must change pointwise by more than `1e-12`, become non-finite, or fail: an unchanged finite field proves that the recorded solves were discarded dummy calls. This causal audit is skipped only when no library solve is observable, which triggers the manual-review path described above. - All saved values must be truthful and internally consistent — the grader recomputes the exact solution and the order from your arrays, re-runs your solver on all five meshes, and compares its finest numerical field pointwise with `U_finest`. If you cannot reach Tier 3, still write honest Tier 1 / Tier 2 outputs rather than fabricating the finest-mesh arrays. - Do not modify files outside `base/output`. (Running the staged `python.sh` wrapper materialises its virtual environment in a temporary location **outside** the task directory — that is expected and permitted, and is not a modification of task files. The staged inputs are left untouched, so they may be mounted read-only.)