Loading of periodic meshes generated by gmsh by miguelcoolchips · Pull Request #4935 · firedrakeproject/firedrake

miguelcoolchips · 2026-03-01T17:09:29Z

Summary

This PR adds support for loading periodic meshes created by Gmsh (.msh files with $Periodic sections) into Firedrake. This addresses #1852 and #1246.

Note

This PR was generated by AI. The tests seem to be working, but frankly I am not sure if the implementation is correct. Any feedback is welcome. I will do my best to address the comments.

Problem

When PETSc loads a Gmsh file with periodic boundaries via DMPlexCreateFromFile, it calls DMLocalizeCoordinates which creates cell-local (DG) coordinates only for cells that straddle the periodic boundary. However, Firedrake's reordered_coords function expects every cell to have an entry in the cell coordinate vector when DG coordinates are present, causing a ValueError reshape failure.

Solution

A new function _fully_localize_coordinates(dm) is added to firedrake/mesh.py and called in make_mesh_from_mesh_topology before the coordinate element type is determined. It expands PETSc's partially localized coordinates to cover all cells by filling in missing entries from CG vertex coordinates via vecGetClosure.

The function includes early-return guards so it is a no-op for:

Non-periodic meshes (getCoordinatesLocalized() returns False)
Already fully-localized periodic meshes (e.g. Firedrake's built-in PeriodicRectangleMesh)

Tests

Three Gmsh geometry files and corresponding pre-generated .msh files (MSH 2 format) are included as test fixtures:

p2d.geo / p2d.msh — 2D rectangle, periodic in x
p2d_xy.geo / p2d_xy.msh — 2D rectangle, periodic in x and y
p3d.geo / p3d.msh — 3D box, periodic in x

Each test solves a Poisson problem and verifies solution periodicity using PointEvaluator at paired probe points across the periodic boundaries. All tests run in both serial and parallel (2 MPI ranks).

Test plan

test_periodic_2d_coordinates — DG coordinate element is used
test_periodic_2d_x_solve / _parallel — Poisson on x-periodic 2D mesh
test_periodic_2d_xy_solve / _parallel — Poisson on doubly-periodic 2D mesh
test_periodic_3d_solve / _parallel — Poisson on x-periodic 3D mesh
Existing regression tests (Helmholtz, non-periodic Gmsh, utility periodic meshes) still pass

PETSc's DMLocalizeCoordinates only creates cell-local (DG) coordinates for cells that straddle the periodic boundary, but Firedrake's reordered_coords expects every cell to have an entry. Add _fully_localize_coordinates() which expands the partial cell coordinate vector to cover all cells, filling non-periodic cells from the CG vertex data via vecGetClosure. The function is called from make_mesh_from_mesh_topology before the existing DG coordinate path runs, so periodic Gmsh meshes now work transparently with Mesh("file.msh") -- no special parameters needed. Closes firedrakeproject#1246 Closes firedrakeproject#1852 Made-with: Cursor

Every solve test now asserts that the solution matches across the periodic boundaries using PointEvaluator probe pairs, rather than only checking that the solution is nonzero. Made-with: Cursor

pbrubeck · 2026-03-01T17:15:55Z

Each test solves a Poisson problem and verifies solution periodicity using PointEvaluator at paired probe points across the periodic boundaries

If we are solving Poisson, it'd be preferable to remove PointEvaluator from these tests and use errornorm with a manufactured (periodic and polynomial) solution

Replace PointEvaluator-based periodicity checks with manufactured solution approach using errornorm, following test_periodic_2d.py. - x-periodic cases (2D, 3D): polynomial manufactured solutions (quadratic in y/z, constant in x) with tight tolerances - doubly-periodic case: trigonometric manufactured solution with refined mesh (h=0.05) for adequate resolution per wavelength - All tests solve Helmholtz (-div(grad(u)) + u = f) to avoid null space issues Made-with: Cursor

tests/periodic/test_periodic.py

Co-authored-by: Pablo Brubeck <brubeck@protonmail.com>

pbrubeck · 2026-03-02T07:15:49Z

tests/periodic/test_periodic.py

+    x and y.  No boundary conditions needed.
+
+    Uses a wider tolerance than the polynomial tests because a
+    non-trivial doubly-periodic solution must be trigonometric,


The function at the boundary is just zero, and for the other polynomial solutions the boundary data is still just constant. If you use cosine instead of sine you could get non-constant boundary conditions.

I would suggest something like

Function(V).interpolate(cos(pi*x)*cos(pi*y))

Good point. I changed it to cosine.

connorjward · 2026-03-02T09:23:07Z

firedrake/mesh.py

+    PETSc's ``DMLocalizeCoordinates`` only creates cell-local (DG)
+    coordinates for cells that straddle the periodic boundary.  The


This isn't quite true. If 'sparse localize' is false then the whole coordinate field is localised. I think DMLocalizeCoordinates should be able to do the right thing for us.

Is it then just a matter of exposing DMSetSparseLocalize with petsc4py and setting sparseLocalize = False? Or maybe just set the option dm_sparse_localize 0?

I also realized that this implementation iterates over the elements. Definitely not the most efficient.

@connorjward are we happy with the current implementation minus the fact that could be done in cython to avoid iterating over cells in python? I can do that next

I've asked Matt for advice. Doing it in Cython might be good but honestly I have found that the difference is performance isn't that great.

connorjward · 2026-03-03T10:07:06Z

firedrake/mesh.py

+        else:
+            needs_expansion = True
+
+    if not needs_expansion or dofs_per_cell is None:


I don't really understand this bit. What case is this capturing?

connorjward · 2026-03-03T10:08:06Z

firedrake/mesh.py

+    For file-based periodic meshes (e.g. Gmsh) the localization
+    happens inside ``setFromOptions`` with ``sparseLocalize = True``


This is wrong. setFromOptions is not called in that routine: https://petsc.org/release/src/dm/impls/plex/plexgmsh.c.html#DMPlexCreateGmshFromFile

connorjward · 2026-03-03T10:09:29Z

firedrake/mesh.py

+        Remove this function once petsc4py exposes
+        ``DMSetSparseLocalize``; at that point we can set
+        ``sparseLocalize = False`` before ``setFromOptions`` and let
+        PETSc handle full localization natively.


I have asked Matt for input on this. Adding bindings to petsc4py is easy so I could just do that once we know the right approach.

pbrubeck · 2026-03-03T23:14:30Z

tests/periodic/test_periodic.py

+def test_periodic_2d_coordinates(periodic_2d_mesh):
+    """Mesh uses a DG coordinate element after loading."""
+    elem = periodic_2d_mesh.ufl_coordinate_element()
+    assert "DG" in str(elem)


Suggested change

assert "DG" in str(elem)

assert elem.family() == "Discontinuous Galerkin"

pbrubeck · 2026-03-03T23:18:35Z

tests/periodic/test_periodic.py

+    f_coeff = (2 * pi / Lx) ** 2 + (2 * pi / Ly) ** 2 + 1.0
+    f = Function(V).assign(f_coeff * u_exact)


Suggested change

f_coeff = (2 * pi / Lx) ** 2 + (2 * pi / Ly) ** 2 + 1.0

f = Function(V).assign(f_coeff * u_exact)

pbrubeck · 2026-03-03T23:19:01Z

tests/periodic/test_periodic.py

+    u = TrialFunction(V)
+    v = TestFunction(V)
+    a = (inner(grad(u), grad(v)) + inner(u, v)) * dx
+    L = inner(f, v) * dx


Suggested change

L = inner(f, v) * dx

L = a(v, u_exact)

pbrubeck · 2026-03-03T23:19:42Z

tests/periodic/test_periodic.py

+    L = inner(f, v) * dx
+
+    uh = Function(V)
+    solve(a == L, uh, solver_parameters={"ksp_type": "cg"})


Suggested change

solve(a == L, uh, solver_parameters={"ksp_type": "cg"})

solve(a == L, uh)

pbrubeck · 2026-03-03T23:21:35Z

tests/periodic/test_periodic.py

+
+    l2err = sqrt(assemble(inner(uh - u_exact, uh - u_exact) * dx))
+    l2norm = sqrt(assemble(inner(u_exact, u_exact) * dx))
+    assert l2err / l2norm < 0.15


This should give zero if u_exact is a Function in V. Additionally, it is more convenient to use errornorm and norm in the lines above

miguelcoolchips added 3 commits March 1, 2026 12:13

Verify periodicity in all periodic Gmsh mesh tests

08a0413

Every solve test now asserts that the solution matches across the periodic boundaries using PointEvaluator probe pairs, rather than only checking that the solution is nonzero. Made-with: Cursor

parallel tests

59ef529

miguelcoolchips added 2 commits March 1, 2026 18:15

fix docstring

1ba5374

pbrubeck reviewed Mar 1, 2026

View reviewed changes

tests/periodic/test_periodic.py Outdated Show resolved Hide resolved

tests/periodic/test_periodic.py Outdated Show resolved Hide resolved

tests/periodic/test_periodic.py Outdated Show resolved Hide resolved

tests/periodic/test_periodic.py Outdated Show resolved Hide resolved

miguelcoolchips and others added 6 commits March 2, 2026 05:54

Update tests/periodic/test_periodic.py

3c3653f

Co-authored-by: Pablo Brubeck <brubeck@protonmail.com>

Update tests/periodic/test_periodic.py

7027ce1

Co-authored-by: Pablo Brubeck <brubeck@protonmail.com>

Update tests/periodic/test_periodic.py

0dd3b99

Co-authored-by: Pablo Brubeck <brubeck@protonmail.com>

Update tests/periodic/test_periodic.py

9b58211

Co-authored-by: Pablo Brubeck <brubeck@protonmail.com>

revert

458553a

fix this test

e9f14f1

pbrubeck reviewed Mar 2, 2026

View reviewed changes

non-constant boundary data

fd0d1c7

connorjward reviewed Mar 2, 2026

View reviewed changes

edit comment

57508e2

connorjward reviewed Mar 3, 2026

View reviewed changes

pbrubeck reviewed Mar 3, 2026

View reviewed changes

		PETSc's ``DMLocalizeCoordinates`` only creates cell-local (DG)
		coordinates for cells that straddle the periodic boundary. The

		For file-based periodic meshes (e.g. Gmsh) the localization
		happens inside ``setFromOptions`` with ``sparseLocalize = True``

	assert "DG" in str(elem)
	assert elem.family() == "Discontinuous Galerkin"

		f_coeff = (2 * pi / Lx) ** 2 + (2 * pi / Ly) ** 2 + 1.0
		f = Function(V).assign(f_coeff * u_exact)

	solve(a == L, uh, solver_parameters={"ksp_type": "cg"})
	solve(a == L, uh)

Conversation

miguelcoolchips commented Mar 1, 2026

Summary

Note

Problem

Solution

Tests

Test plan

Uh oh!

pbrubeck commented Mar 1, 2026

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

miguelcoolchips Mar 3, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

miguelcoolchips Mar 3, 2026 •

edited

Loading