Can't label facets on periodic meshes after construction

[pefarrell]

I have a code that uses interior facet integrals on a PeriodicUnitSquareMesh. It works when I use ... * dS. I then colour all interior facets with colour 1. The code does not work if I use ... * dS(1) instead. On investigation, this issue appears to be specific to periodic domains; on a UnitSquareMesh the issue does not arise.

Minimal reproducing code:

from firedrake import *
from firedrake.cython import dmplex
from firedrake.petsc import PETSc
import numpy

# Colour all interior facets to be colour 1
def mark_mesh(mesh):
    dm = mesh._plex
    sec = dm.getCoordinateSection()
    coords = dm.getCoordinatesLocal()
    dm.removeLabel(dmplex.FACE_SETS_LABEL)
    dm.createLabel(dmplex.FACE_SETS_LABEL)

    faces = dm.getStratumIS("interior_facets", 1).indices
    for face in faces:
        dm.setLabelValue(dmplex.FACE_SETS_LABEL, face, 1)
    return mesh

def run_cgh(N, r, periodic=False, use_markers=False):
    if periodic:
        mesh = PeriodicUnitSquareMesh(N, N, direction="x")
    else:
        mesh = UnitSquareMesh(N, N)

    mesh = mark_mesh(mesh)
    x = SpatialCoordinate(mesh)

    # symbolic expressions for exact solution and f
    u_e = as_vector([sin(pi*x[0])*sin(pi*x[1]), sin(pi*x[0])*sin(pi*x[1])])
    f_e = -div(grad(u_e))
    V_a = VectorFunctionSpace(mesh, 'DG', r + 3)
    f = interpolate(f_e, V_a)

    # set up variational problem
    V_element = FiniteElement('CG', mesh.ufl_cell(), r)
    V = FunctionSpace(mesh, VectorElement(BrokenElement(V_element)))
    Qhat = FunctionSpace(mesh, VectorElement(BrokenElement(V_element[facet])))
    Vhat = FunctionSpace(mesh, VectorElement(V_element[facet]))
    W = V * Qhat * Vhat

    z = Function(W)
    (u, qhat, uhat) = split(z)

    J = 0.5 * inner(grad(u), grad(u))*dx - inner(f, u)*dx
    L_match = inner(qhat, uhat - u)

    if use_markers:
        L = J + (L_match('+') + L_match('-'))*dS(1) + L_match*ds
    else:
        L = J + (L_match('+') + L_match('-'))*dS    + L_match*ds
    F = derivative(L, z, TestFunction(W))
    bcs = [DirichletBC(W.sub(2), 0, 'on_boundary')]

    params = {
        'snes_type': 'ksponly',
        'ksp_type': 'preonly',
        'pc_type': 'svd',
        'pc_svd_monitor': None,
        'pc_factor_mat_solver_type': 'mumps',
        'mat_type': 'aij'
    }

    solve(F == 0, z, bcs, solver_parameters=params)

# It works except when periodic = True and use_markers = True
for periodic in [False, True]:
    for use_markers in [False, True]:
        print("periodic = %s, use_markers = %s" % (periodic, use_markers))
        run_cgh(N=4, r=2, periodic=periodic, use_markers=use_markers)

which prints

periodic = False, use_markers = False
      SVD: condition number 4.483440553745e+03, 0 of 930 singular values are (nearly) zero
      SVD: smallest singular values: 1.058848554785e-03 1.058848554785e-03 1.063620151543e-03 1.063620151543e-03 1.066708024458e-03
      SVD: largest singular values : 4.747284006177e+00 4.747284162579e+00 4.747284162579e+00 4.747284550799e+00 4.747284550799e+00
periodic = False, use_markers = True
      SVD: condition number 4.483440553745e+03, 0 of 930 singular values are (nearly) zero
      SVD: smallest singular values: 1.058848554785e-03 1.058848554785e-03 1.063620151543e-03 1.063620151543e-03 1.066708024458e-03
      SVD: largest singular values : 4.747284006177e+00 4.747284162579e+00 4.747284162579e+00 4.747284550799e+00 4.747284550799e+00
periodic = True, use_markers = False
      SVD: condition number 4.483440725012e+03, 0 of 912 singular values are (nearly) zero
      SVD: smallest singular values: 1.058848554785e-03 1.058848554785e-03 1.074053173472e-03 1.074053173472e-03 1.074053173472e-03
      SVD: largest singular values : 4.747284166340e+00 4.747284229622e+00 4.747284229622e+00 4.747284732145e+00 4.747284732145e+00
periodic = True, use_markers = True
      SVD: condition number            inf, 496 of 912 singular values are (nearly) zero
      SVD: smallest singular values: 0.000000000000e+00 0.000000000000e+00 0.000000000000e+00 0.000000000000e+00 0.000000000000e+00
      SVD: largest singular values : 4.742944178293e+00 4.742944178293e+00 4.742944178293e+00 4.742944178293e+00 4.742944178293e+00

i.e. the use_markers True vs False makes no difference on UnitSquareMesh, but does on PeriodicUnitSquareMesh.

[pefarrell]

Fixed in pull request #1596 .

[wence-]

I think this is basically a cache invalidation problem. Periodic...Mesh gets further through mesh construction than Unit...Mesh and so the facet labels are baked in (your mesh marking happens too late).

We could fix this by making the subset and boundary node creation that goes along with mesh labels "dynamic" and not cached.

I guess this is still a problem because #1596 was not merged.