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

Aperta
  • Qualcosa non va?
#5,266 1 commento 0 reazioni 1 assegnatario Vedi su GitHub

I maintainer di solito rispondono entro 1 giorno

@connorjward ci sta già lavorando.

Dal 21/7/2026.

Valutazione

Questa issue non è ancora stata valutata.

Descrizione

bug

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}"
Lingua principale
Python
Stelle
684
Fork
204
Merge medio
4g 14h
PR unite (30g)
46

Preparare l'ambiente

  • Nessun Dockerfile né file Docker Compose
  • Ha un modello di pull request
  • Nessuna guida per i contributori

Come iniziare

  1. Leggi tutta la issue e poi la guida ai contributi del progetto.
  2. Commenta sulla issue per dire che te ne occupi tu — evita che due persone facciano lo stesso lavoro.
  3. Fai un fork del repository e lavora su un branch.
  4. Apri una pull request che faccia riferimento al numero della issue.

Altre issue di firedrakeproject/firedrake

Tutte le issue di firedrakeproject/firedrake

Issue simili

Altre issue su Python

Ricevi le nuove issue nella tua casella

Un breve riepilogo di issue GitHub adatte ai principianti.