Stencil Operator Examples

Stencil routines are useful both as direct updates and as matrix-free linear operators.

Apply a 2D Laplacian

num::Vector u(N * N, 0.0);
num::Vector Lu(N * N, 0.0);
 
num::pde::laplacian_stencil_2d(u, Lu, N);

Matrix-Free Operator

auto A = num::operators::make_op(
    [N](const num::Vector& x, num::Vector& y) {
        num::pde::laplacian_stencil_2d(x, y, N);
        num::scale(y, -1.0);
    },
    N * N);
 
num::SolverResult info = num::cg(A, rhs, sol, 1e-8, 1000);

The operator represents

\[ y = -L_h x . \]

Periodic Diffusion Step

double coeff = kappa * dt / (h * h);
num::pde::diffusion_step_2d(u, N, coeff, num::best_backend);

Fourth-Order Dirichlet Laplacian

num::Vector Lu4(N * N, 0.0);
num::pde::laplacian_stencil_2d_4th(u, Lu4, N);

Use matrix-free operators when the sparse matrix would be expensive to assemble or when the stencil is applied many times inside a Krylov method.