PDE Module

Header: #include "pde/pde.hpp" Namespace: num, num::pde

Finite-difference PDE infrastructure: stencil operators, backward-Euler matrix assembly, Crank-Nicolson ADI, and explicit diffusion steps. Grid geometry and field storage live in fields/; time integration lives in ode/. The pde/ module is spatial operators only.

---

Contents

Header	What it provides
pde/stencil.hpp	2D/3D Laplacians, gradient, divergence, curl, fiber sweeps
pde/diffusion.hpp	backward_euler_matrix, make_cg_solver; implicit backward-Euler setup
pde/fields.hpp	ScalarField3D, VectorField3D, FieldSolver, MagneticSolver
pde/adi.hpp	CrankNicolsonADI – Strang-split CN solver for parabolic equations
pde/poisson.hpp	poisson2d_fd, poisson2d – DST-based 2D Poisson solver, O(N^2 log N)

---

Stencil operators (<tt>pde/stencil.hpp</tt>)

2D operators

// Dirichlet BCs (zero on boundary), NxN grid row-major (idx = i*N + j)
void num::laplacian_stencil_2d(const Vector& x, Vector& y, int N);
 
// Periodic BCs
void num::laplacian_stencil_2d_periodic(const Vector& x, Vector& y, int N);

All three stencils return h² ∇²u — divide by h² to recover the Laplacian:

Function	Points	Order
laplacian_stencil_2d	5 (cross)	O(h²), Dirichlet
laplacian_stencil_2d_periodic	5 (cross)	O(h²), periodic
laplacian_stencil_2d_4th	13 (extended cross)	O(h⁴), Dirichlet

The standard 5-point stencil:

\[y_{i,j} = x_{i+1,j} + x_{i-1,j} + x_{i,j+1} + x_{i,j-1} - 4\,x_{i,j}\]

The 4th-order 13-point stencil:

\[ y_{i,j} = \frac{1}{12}\bigl( -x_{i-2,j} + 16\,x_{i-1,j} - 30\,x_{i,j} + 16\,x_{i+1,j} - x_{i+2,j} -x_{i,j-2} + 16\,x_{i,j-1} + 16\,x_{i,j+1} - x_{i,j+2} \bigr) \]

2D grid utilities

// Initialise an NxN field from a callable f(x,y) -> real
template<typename F>
void num::fill_grid(Vector& u, int N, double h, F&& f);
 
// Extract a row/column as a plottable Series (node k at (k+1)*h)
num::Series num::row_slice(const Vector& u, int N, double h, int row);
num::Series num::col_slice(const Vector& u, int N, double h, int col);

with Dirichlet or periodic treatment at the boundaries.

2D fiber sweeps

Used internally by CrankNicolsonADI and available directly:

// For each column j: extract fiber psi[*,j], call f(fiber), write back
template<typename F>
void num::col_fiber_sweep(CVector& psi, int N, F&& f);
 
// For each row i: extract fiber psi[i,*], call f(fiber), write back
template<typename F>
void num::row_fiber_sweep(CVector& psi, int N, F&& f);

3D operators

All operate on Grid3D objects with central differences and one-sided stencils at boundaries:

void num::neg_laplacian_3d(const Vector& v, Vector& Av, int nx, int ny, int nz, double inv_dx2);
void num::gradient_3d(const Grid3D& phi, Grid3D& gx, Grid3D& gy, Grid3D& gz);
void num::divergence_3d(const Grid3D& fx, const Grid3D& fy, const Grid3D& fz, Grid3D& div);
void num::curl_3d(const Grid3D& Ax, const Grid3D& Ay, const Grid3D& Az,
                  Grid3D& Bx, Grid3D& By, Grid3D& Bz);

neg_laplacian_3d computes \((-\nabla^2 x)\) with Dirichlet BCs (boundary rows = identity), giving an SPD system suitable for CG.

---

3D fields (<tt>pde/fields.hpp</tt>)

See fields.md for full documentation. The classes live in num:::

• ScalarField3D – uniform-grid scalar potential with trilinear sampling
• VectorField3D – three ScalarField3D components
• FieldSolver – static methods: solve_poisson, gradient, divergence, curl
• MagneticSolver – static methods: current_density, solve_magnetic_field

---

Crank-Nicolson ADI (<tt>pde/adi.hpp</tt>)

Prefactored Crank-Nicolson solver for 2D parabolic equations via Strang splitting. Suitable for TDSE, heat equations with complex coefficients.

num::CrankNicolsonADI adi(N, dt, h);

Factorises two ComplexTriDiag systems on construction ( \(O(N)\)):

• td_half_ – for half-step sweeps ( \(\tau = dt/2\))
• td_full_ – for full-step sweep ( \(\tau = dt\))

Sweep

adi.sweep(psi, x_axis, tau);

Applies one CN sub-step along the chosen axis:

• x_axis = true – processes each column fiber (x-direction)
• x_axis = false – processes each row fiber (y-direction)
• Selects td_half_ when tau < 0.75*dt, td_full_ otherwise

Strang splitting (TDSE)

// Full time step:
adi.sweep(psi, true,  dt * 0.5);  // x, half-step
adi.sweep(psi, false, dt);         // y, full-step
adi.sweep(psi, true,  dt * 0.5);  // x, half-step

Each sweep solves the 1D tridiagonal system per fiber:

\[\left(I - i\alpha\,\nabla^2_{1D}\right)\psi^{n+1} = \left(I + i\alpha\,\nabla^2_{1D}\right)\psi^n, \qquad \alpha = \frac{\tau}{4h^2}\]

The LHS tridiagonal has sub/super-diagonal \(-i\alpha\) and main diagonal \(1 + 2i\alpha\); it is factored once on construction and reused for every fiber.

---

Implicit backward Euler (<tt>pde/diffusion.hpp</tt>)

Assembly for the implicit system A = I - coeff*L where L is the 2D Dirichlet Laplacian. A is sparse SPD — conjugate gradient converges in O(N) iterations.

// Build once; reuse A and solver across all time steps
num::Grid2D       grid{N, h};
num::SparseMatrix A      = num::pde::backward_euler_matrix(grid, coeff);
num::LinearSolver solver = num::make_cg_solver(A);
 
// Integrate: solve A u^{n+1} = u^n  for nstep steps
num::solve(u, num::BackwardEuler{.solver=solver, .dt=dt, .nstep=nstep});

coeff = kappa * dt / h². Implicit stepping removes the dt < 0.25 h² Von Neumann stability constraint.

Explicit diffusion (<tt>pde/stencil.hpp</tt>)

Forward Euler diffusion steps for 2D uniform grids. Von Neumann stability requires coeff <= 0.25.

// Periodic BCs:  u += coeff * Lap_periodic(u)
num::pde::diffusion_step_2d(u, N, coeff);
num::pde::diffusion_step_2d(u, N, coeff, num::best_backend);
 
// Dirichlet BCs: u += coeff * Lap_dirichlet(u)
num::pde::diffusion_step_2d_dirichlet(u, N, coeff);

Typical usage (Navier-Stokes viscosity)

const double coeff = dt * nu / (h * h);
num::pde::diffusion_step_2d(u_star, N, coeff, num::best_backend);
num::pde::diffusion_step_2d(v_star, N, coeff, num::best_backend);

---

2D Poisson solver (<tt>pde/poisson.hpp</tt>)

Solves \(-\nabla^2 u = f\) on \((0,1)^2\) with homogeneous Dirichlet BCs via the Discrete Sine Transform. Two variants:

• pde::poisson2d_fd – FD eigenvalues \(D_k = 2(1-\cos(k\pi/(N+1)))\), error \(O(h^2)\).
• pde::poisson2d – exact eigenvalues \((k\pi)^2\), machine-precision error for \(f\) in the DST eigenbasis.

Both run in \(O(N^2 \log N)\) and require \(N = 2^p - 1\). See poisson.md for the full derivation and convergence table.

Matrix u_fd   = pde::poisson2d_fd(f, N);   // O(h^2)
Matrix u_spec = pde::poisson2d(f, N);      // machine precision