BUG: Are facet normals correct for periodic rectangle meshes of width 1 or 2?

[tommbendall]

One of the Gusto "case study" tests seems to have been failing since we picked up the Firedrake 2026.04.1 release. The Skamarock-Klemp hydrostatic gravity wave uses a PeriodicRectangleMesh which is 1 cell wide in the y-direction.

The pressure gradient is calculated via a term which includes cp*avg(exner)*jump(theta*w, n)*dS_v, where:
cp is the specific heat capacity at constant pressure, exner is the "Exner" pressure, theta is the potential temperature, w is the velocity test function and n is the FacetNormal. In our test, the pressure gradient in the y-direction should be essentially zero, but this term now seems to be enormous and causes the model to blown up in the first timestep.

I think the issue is that there has been a change in behaviour for FacetNormal on the PeriodicRectangleMesh, specifically when the mesh is only 1 or 2 cells wide.

To investigate this, I checked whether n takes equal but opposite values on the '+' and '-' sides of each facet. It seems that for the PeriodicRectangleMesh it now doesn't when the mesh is 1 or 2 cells wide.

I wonder if this was a side effect of #4836?

If I change the Gusto test to use a width of 3, it runs successfully.

Here is an MFE:

"""
Minimum failing example demonstrating that facet normals are not computed
correctly on periodic meshes of widths 1 or 2.

Facet normals on the '+' and '-' sides of a facet should sum to zero, as they
should be of equal size but opposite direction.

At the time of writing, this test passes for a width of 3 or more but fails
for a width
of 1 or 2.
"""

from firedrake import (
    PeriodicRectangleMesh, FunctionSpace, Function, FacetNormal,
    TestFunction, TrialFunction, inner, dS,
    LinearVariationalProblem, LinearVariationalSolver, VTKFile
)
import numpy as np

to_output = False
mesh_widths = [4, 3, 2, 1]  # Number of cells in y-direction

for width in mesh_widths:
    # Periodic rectangle mesh
    mesh = PeriodicRectangleMesh(5, width, 10.0, 10.0/width, quadrilateral=True)

    # RTCF1 has one degree of freedom per facet
    HDiv = FunctionSpace(mesh, "RTCF", 1)

    n = FacetNormal(mesh)
    psi = TestFunction(HDiv)
    p = TrialFunction(HDiv)
    n_plus = Function(HDiv, name='n_plus')    # Takes value of n on + side
    n_minus = Function(HDiv, name='n_minus')  # Takes value of n on - side
    n_sum = Function(HDiv, name='n_sum')      # Sum of n+ and n-

    # Project the facet normals onto the HDiv space
    # As psi and p are continuous, their label on the facet doesn't matter
    a = inner(p('+'), psi('+'))*dS
    L_plus = inner(n('+'), psi('+'))*dS
    L_minus = inner(n('-'), psi('+'))*dS
    problem_plus = LinearVariationalProblem(a, L_plus, n_plus)
    problem_minus = LinearVariationalProblem(a, L_minus, n_minus)
    solver_plus = LinearVariationalSolver(problem_plus)
    solver_minus = LinearVariationalSolver(problem_minus)
    solver_plus.solve()
    solver_minus.solve()

    n_sum.project(n_plus + n_minus)

    if to_output:
        output = VTKFile(f"facet_normals_{width}.pvd")
        output.write(n_plus, n_minus, n_sum)

    assert np.allclose(n_sum.dat.data_ro, 0.0, 1e-16), \
        f"Facet normals do not sum to zero for mesh width {width}"

[connorjward]

I did not consider facet normals when I made those changes. I can easily imagine there's a bug there. Happy to investigate.