The Star That Started Moving — Well-Balanced FVM and Angular Momentum in Polar Coordinates

Same equations, same slope limiter, same Riemann solver. Only the grid changed from Cartesian to polar. Gas began bleeding out of the center of a quiescent star. Pressure was uniform. Every velocity was zero.

Drag $\Delta r$ up and the textbook source $2P/r \cdot V$ pulls away from the face-area jump $P(S_{\text{out}} - S_{\text{in}})$.

The gap shrinks as $O(\Delta r^{3})$ and still is not zero. A coarser mesh opens it first.

The grid is the culprit. Curvilinear coordinates add geometric terms to the momentum equations, and one bad line in that discretization kills a static equilibrium. The swap from φ-momentum to angular momentum sits next to that line.

What polar coordinates add to the momentum equation#

In polar coordinates $(r, \theta, \varphi)$ the radial momentum equation reads

where $u, v, w$ are the $r, \theta, \varphi$ velocity components and $P$ is pressure. Two differences from the Cartesian form: ① the divergence is wrapped by $1/r^{2}$, and ② the right side carries a pseudo-force $-\rho(v^{2}+w^{2})/r$.

These pseudo-forces are not real accelerations. The local basis (the vector frame attached to each point) rotates with position, so a momentum vector that points the same way globally has different components in the local frame. But numerically they enter the cell update as if they were forces, and the way we discretize them matters.

Where 2P/r really comes from — face-area mismatch#

If we move $\partial P/\partial r$ inside the conservative divergence, the discrete radial momentum equation becomes

The flux $F_{rr} = \rho u^{2} + P$ is exactly what your Riemann solver returns. The $2P/r \cdot V$ source is today's story.

We get it by rewriting $\partial P / \partial r$ as $\frac{1}{r^{2}}\partial (P r^{2})/\partial r - 2P/r$. The first piece slides into the flux divergence; the second becomes the source. Discretized cell-by-cell, the source reads

That looks innocent. Now drop in a static solution — pressure constant, velocities zero — and check the books. The flux side gives $\Delta_{r}(F_{rr} S_{r}) = P (S_{i+1/2} - S_{i-1/2})$, which is not zero because the inner and outer face areas of a polar cell differ. The source side, the Taylor-derived $2P/r \cdot V$, only matches that difference approximately. The mismatch leaks into the cell as a spurious acceleration every step.

Back to the printed numbers. $N=40$, radius from 1 to 4, 200 steps. Naive $\max|u|$ is $1.834\times 10^{-3}$. Balanced prints $0.000\times 10^{0}$. The initial state is $P=1$, $\rho=1$, $u=0$, $\gamma=1.4$. $\Delta r=(4-1)/40=0.075$. The code steps with $\Delta t=0.4\,\Delta r/c_0$ and $c_0=\sqrt{\gamma P/\rho}$. Two hundred of those steps pile a spurious velocity up to $1.834\times 10^{-3}$. The balanced form does not move inside the same loop.

Refining to 200 cells cuts the velocity residual as $(\Delta r)^{2}$, not to zero. The figure quotes an $O(\Delta r^{3})$ gap. One order apart. The cell source error is the Taylor remainder of the area mismatch, so it is smaller relative to the volume. Divide by $\rho V$ and it is an acceleration. Repeat $\Delta t$ and it becomes a velocity. That is why the velocity residual looks one order softer. Swap in $P(S_{i+1/2}-S_{i-1/2})$ and the order talk ends, because the flux difference and the source are the same expression.

Volume is $V=(r_R^{3}-r_L^{3})/3$ and area is $S=r^{2}$. The factor $\sin\theta\,\Delta\theta\,\Delta\varphi$ is omitted on both sides, as the code comment says. Equal factors keep the cancellation. Change only the area definition and $P\cdot\Delta S$ no longer matches the flux difference. The same hole shows up when a curved-cell area comes from a diagonal cross product: $\partial x/\partial x$ leaves 1. Discrete divergence and face area have to be a pair, or a static solution will not hold.

One-line fix that restores equilibrium#

Replace the source term with the same expression as the flux difference:

Now the flux contribution and the source cancel to machine precision for a static solution. The star stops leaking.

The same trick applies to the $\cos\theta P / (r \sin\theta)$ source in the $\theta$-momentum equation:

The $\varphi$-momentum equation has no such geometric pressure source — the two azimuthal faces share the same area.

This is the simplest example of a well-balanced scheme: the discrete divergence's fake acceleration is subtracted by the same discrete expression on the other side.

Sixty lines of NumPy — does the equilibrium hold?#

On a 1D radial grid we start with $P = 1$, $\rho = 1$, $u = 0$ and watch the momentum residual under each formulation.

import numpy as np
 
def polar_cell_volume(rL: float, rR: float) -> float:
    """spherical-shell volume between rL and rR (sin θ Δθ Δφ omitted)."""
    return (rR**3 - rL**3) / 3.0
 
def polar_r_surface(r: float) -> float:
    """r²-weighted face area in spherical coords."""
    return r * r
 
def naive_pressure_source(P: float, r_mid: float, V: float) -> float:
    """The textbook 2P/r · V form — wrong at machine precision for P=const."""
    return 2.0 * P / r_mid * V
 
def balanced_pressure_source(P: float, S_in: float, S_out: float) -> float:
    """Same expression as the flux difference, so they cancel exactly."""
    return P * (S_out - S_in)
 
def radial_static_run(N: int = 40, scheme: str = "balanced", steps: int = 200):
    rIn, rOut = 1.0, 4.0
    dr = (rOut - rIn) / N
    rho = np.ones(N)
    u = np.zeros(N)
    P = np.ones(N)
    gamma = 1.4
    c0 = np.sqrt(gamma * P[0] / rho[0])
    dt = 0.4 * dr / c0
 
    r_face = np.linspace(rIn, rOut, N + 1)
    S = polar_r_surface(r_face)
    V = np.array([polar_cell_volume(r_face[i], r_face[i + 1]) for i in range(N)])
    r_mid = 0.5 * (r_face[:-1] + r_face[1:])
 
    for _ in range(steps):
        # pressure flux through faces (P constant → exact value)
        flux_diff = P * (S[1:] - S[:-1])
        # the source-term line we want to scrutinise
        if scheme == "naive":
            src = naive_pressure_source(P, r_mid, V)
        else:
            src = balanced_pressure_source(P, S[:-1], S[1:])
        # momentum update: d(ρ u V)/dt = -flux_diff + src
        du = dt / (rho * V) * (src - flux_diff)
        u += du
    return r_mid, u
 
r1, u_naive    = radial_static_run(N=40, scheme="naive",    steps=200)
r2, u_balanced = radial_static_run(N=40, scheme="balanced", steps=200)
print(f"naive    : max |u| = {np.abs(u_naive).max():.3e}")
print(f"balanced : max |u| = {np.abs(u_balanced).max():.3e}")

After 200 steps on a 40-cell grid:

naive    : max |u| = 1.834e-03
balanced : max |u| = 0.000e+00

Leave the toggle on naive and the curve walks off zero as steps accumulate. Switch to well-balanced and it stays on the axis. Drop cells and only the naive residual grows. On 40 cells and 200 steps that walk is $\max|u|=1.834\times 10^{-3}$; balanced prints $0.000\times 10^{0}$.

Replace φ-momentum with angular momentum#

Pressure is not the only source. The $\varphi$-momentum equation carries $\rho u w / r + \cos\theta \rho v w / r$. In strongly rotating problems — stars, accretion disks, atmospheres — these terms erode accuracy and conservation.

The cure is a change of variable. Define $l = w r \sin\theta$ — plain angular momentum — and the source terms in the $\varphi$-momentum equation vanish entirely:

Discretely we take the Riemann fluxes $F_{\varphi r}$, $F_{\varphi\theta}$, $F_{\varphi\varphi}$, multiply by $r \sin\theta$ to get $\tilde F$, and difference them. No corrections, no source-term juggling. Even on a coarse mesh angular momentum is conserved by construction (Kley 1998, A&A 338, L37).

Differentially rotating systems (Keplerian disks) need an additional rotating-frame trick — FARGO — to relax the CFL bound. That is a separate post.

A decision guide#

Standing in front of the star again#

A Cartesian grid hides this trap. Curvilinear coordinates add terms to the momentum equation, and unless their discrete form exactly matches the discrete flux difference, a quiet gas will not stay quiet. The one-line edit of this post is

One line, one commit. The star is quiet again.

Keep 2P/r, or replace it with ΔS?#

If only the pressure balance wobbles, replace $2P/r\,V$ with $P(S_{i+1/2}-S_{i-1/2})$. That is the edit that turns $\max|u|=1.834\times 10^{-3}$ on 40 cells into 0. The $\theta$ source $\cos\theta\,P/(r\sin\theta)$ takes the same swap, $P\cdot\Delta S_{\theta}$. The $\varphi$ direction does not need it: the two azimuthal faces share one area.

Strong rotation is a different decision. $\rho u w/r + \cos\theta\,\rho v w/r$ still sits on the $\varphi$-momentum equation. Defining $l = w r \sin\theta$ removes that source from the equation. That is the form in Kley 1998 (A&A 338, L37). A radially varying rotation still wants FARGO on top. A static interface that starts to flow on its own is the same disease as parasitic currents. Unless the fake force from the discrete divergence is subtracted with the same discrete expression, the equilibrium does not come back.