Skip to content

Conserved quantities and budgets

Questions: does the code keep what the physics keeps? And does it keep exact budgets exactly?

Energy drift that shrinks at the right order

A 2nd-order integrator may let the energy of an oscillator drift, but only by O(h^2). invariant_study measures the drift at several step sizes and checks its order.

"""Energy of a harmonic oscillator: velocity Verlet keeps it to O(h^2); a wrong velocity update does not."""
import numpy as np

from physverify.studies import invariant_study


def make(buggy):
    def integrate(source, y0, t_end, n):  # velocity Verlet for x'' = -x
        h = t_end / n
        x, v = np.array(y0, dtype=float)
        for _ in range(n):
            a = -x
            x = x + h * v + 0.5 * h * h * a
            a_new = -x
            v = v + (h * a if buggy else 0.5 * h * (a + a_new))  # bug: only the old acceleration
        return np.array([x, v])
    return integrate


def energy(y):
    return 0.5 * (y[0] ** 2 + y[1] ** 2)


for buggy in (False, True):
    r = invariant_study(make(buggy), energy, np.array([1.0, 0.0]), 10.0, [50, 100, 200, 400], 2)
    print("buggy:  " if buggy else "correct:", r.status, "-", r.message)

Output:

correct: PASS - observed order 2.02 (expected 2)
buggy:   FAIL - observed order 1.21 is below the expected 2

The buggy velocity update uses only the old acceleration. The positions still look like an oscillator, but the energy drifts at about first order.

Budgets kept exactly on every grid

Some schemes conserve a quantity exactly: a finite-volume heat solver with insulated walls keeps the total heat to round-off, on every grid and for every step. A convergence study is the wrong test there, because a leak that shrinks with the grid would pass it. balance_check asks for the exact budget at every resolution.

"""Heat budget of a 1D explicit heat solver with insulated walls: the total heat must not change."""
import numpy as np

from physverify.studies import balance_check


def initial(n):
    x = (np.arange(n) + 0.5) / n  # cell centres on [0, 1]
    return np.exp(-50 * (x - 0.3) ** 2)


def heat_change(n, buggy, steps=200):
    dx = 1.0 / n
    u = initial(n)
    dt = 0.4 * dx * dx
    for _ in range(steps):
        if buggy:  # "insulated" by copying the neighbour value after the update
            lap = np.zeros_like(u)
            lap[1:-1] = u[2:] - 2 * u[1:-1] + u[:-2]
            u = u + dt / dx**2 * lap
            u[0], u[-1] = u[1], u[-2]
        else:  # zero flux through both walls, written as fluxes between cells
            flux = np.zeros(n + 1)
            flux[1:-1] = -(u[1:] - u[:-1]) / dx
            u = u - dt / dx * (flux[1:] - flux[:-1])
    return (u.sum() - initial(n).sum()) * dx


for buggy in (False, True):
    r = balance_check(lambda n: heat_change(n, buggy), 0.0, [40, 80, 160], tol=1e-12)
    print("buggy:  " if buggy else "correct:", r.status, "-", r.message)

Output:

correct: PASS - balance 0 kept to 1e-12 at every resolution (n=40: -4.44089e-17, n=80: 4.44089e-17, n=160: 4.44089e-17)
buggy:   FAIL - expected 0, got n=40: 0.0119969, n=80: 0.00389703, n=160: 0.000583943; the mismatch shrinks with resolution, so the scheme converges but does not conserve

The buggy wall copies the neighbour value after each update, a common way to write an "insulated" boundary. It converges (the mismatch shrinks) but does not conserve heat. A real code written this way lost 8.5% of its heat.

Other budgets worth checking

  • Heat gained through a wall with a known flux, or from a known source, over a run.
  • Total mass in an advection code with periodic boundaries.
  • Momentum in an N-body code (exact for pairwise forces that obey Newton's third law).
  • Angular momentum in a central-force orbit.