Poisson's equation on cells

Download notebookView on GitHubRun when the site was built.

solve_poisson solves $-\nabla^2 \phi = f$ with one value per cell, on all patches at once: a conservative finite-volume scheme in which a coarse cell next to a patch exchanges flux with the fine cells along their shared face. The boundary value is one constant on the whole boundary.

import numpy as np

from patchwork import AMRHierarchy, plot_patches, solve_poisson

A problem with a known answer

On the unit square, $\phi = \sin(\pi x)\sin(\pi y)$ vanishes on the boundary and has $f = 2\pi^2 \phi$. The right-hand side is a cell-centered field, set from the cell centres.

def exact(x, y):
    return np.sin(np.pi * x) * np.sin(np.pi * y)


def solve(base, patches=()):
    mesh = AMRHierarchy(base_shape=(base, base))
    for level, box in patches:
        mesh.add_patch(level=level, box=box)
    mesh.allocate_field("rhs")
    for patch in mesh:
        x, y = patch.local_cell_centers(mesh)
        patch.field("rhs").local[...] = 2 * np.pi**2 * exact(x, y)
    solve_poisson(mesh, "rhs", potential="phi")
    error = max(
        np.abs(patch.field("phi").local - exact(*patch.local_cell_centers(mesh))).max()
        for patch in mesh
    )
    return mesh, error


mesh, error = solve(32, [(1, ((16, 48), (16, 48)))])
print(f"largest error: {error:.2e}")
fig, ax = plot_patches(mesh, field="phi")

largest error: 5.00e-03

solve_poisson writes the solution into the field named by potential, allocating it if needed. The cells of a coarse patch that a finer patch covers hold the mean of the fine cells over them.

Convergence

Halving the cells divides the error by about four: the scheme is second order.

previous = None
for base in (16, 32, 64):
    _, error = solve(base)
    ratio = "" if previous is None else f", {previous / error:.1f} times smaller"
    print(f"{base} x {base} cells: error {error:.2e}{ratio}")
    previous = error

16 x 16 cells: error 3.19e-03
32 x 32 cells: error 8.02e-04, 4.0 times smaller
64 x 64 cells: error 2.01e-04, 4.0 times smaller

For faces with different conditions, periodic axes or electrodes inside the domain, use the node-centered solver of the next tutorial.