BUG: Compilation error for MG transfer of HDiv element on immersed mesh

[JHopeCollins]

Describe the bug
I'm getting a compilation error from pyop2 when trying use MG transfer operators with HDiv functions on immersed meshes with firedrake-2026.4.0

Steps to Reproduce

from firedrake import *
basemesh = UnitIcosahedralSphereMesh(
    refinement_level=0, degree=1
)
mc, mf = MeshHierarchy(
    basemesh, refinement_levels=1)
mc.init_cell_orientations(SpatialCoordinate(mc))
mf.init_cell_orientations(SpatialCoordinate(mf))
family = "BDM"
degree = 1
Vc = FunctionSpace(mc, family, degree)
Vf = FunctionSpace(mf, family, degree)
uc = Function(Vc)
uf = Function(Vf)
transfers = fd.TransferManager()
transfers.inject(uf, uc)

Expected behavior
This worked fine in firedrake-2025.10.3

Error message
Running with PYOP2_DEBUG=1 python mfe.py gives the following error, and similar errors if you try transfers.restrict or transfers.prolong with Cofunctions.

Traceback (most recent call last):
  File "/home/jhc/codes/fd/fd-main/src/asQ/utils/mg.py", line 119, in <module>
    transfers.inject(uf, uc)
  File "/home/jhc/codes/fd/fd-main/src/firedrake/firedrake/mg/embedded.py", line 294, in inject
    self.op(uf, uc, transfer_op=Op.INJECT)
  File "petsc4py/PETSc/Log.pyx", line 250, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "petsc4py/PETSc/Log.pyx", line 251, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "/home/jhc/codes/fd/fd-main/src/firedrake/firedrake/mg/embedded.py", line 245, in op
    self._native_transfer(target_element, transfer_op)(source, target)
  File "petsc4py/PETSc/Log.pyx", line 250, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "petsc4py/PETSc/Log.pyx", line 251, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "/home/jhc/codes/fd/fd-main/src/firedrake/firedrake/mg/interface.py", line 241, in inject
    op2.par_loop(kernel, coarse.node_set,
  File "/home/jhc/codes/fd/fd-main/src/firedrake/pyop2/parloop.py", line 747, in par_loop
    parloop(*args, **kwargs)
  File "/home/jhc/codes/fd/fd-main/src/firedrake/pyop2/parloop.py", line 760, in parloop
    LegacyParloop(knl, *args, **kwargs)()
  File "petsc4py/PETSc/Log.pyx", line 250, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "petsc4py/PETSc/Log.pyx", line 251, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "/home/jhc/codes/fd/fd-main/src/firedrake/pyop2/parloop.py", line 251, in __call__
    self._compute(self.iterset.core_part)
  File "/home/jhc/codes/fd/fd-main/src/firedrake/pyop2/parloop.py", line 232, in _compute
    self.global_kernel(self.comm, part.offset, part.offset+part.size, *self.arglist)
  File "/home/jhc/codes/fd/fd-main/src/firedrake/pyop2/global_kernel.py", line 334, in __call__
    func = compile_global_kernel(self, comm)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "petsc4py/PETSc/Log.pyx", line 250, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "petsc4py/PETSc/Log.pyx", line 251, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "/home/jhc/codes/fd/fd-main/src/firedrake/pyop2/caching.py", line 583, in wrapper
    value = func(*args, **kwargs)
            ^^^^^^^^^^^^^^^^^^^^^
  File "/home/jhc/codes/fd/fd-main/src/firedrake/pyop2/global_kernel.py", line 445, in compile_global_kernel
    dll = load(
          ^^^^^
  File "petsc4py/PETSc/Log.pyx", line 250, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "petsc4py/PETSc/Log.pyx", line 251, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "/home/jhc/codes/fd/fd-main/src/firedrake/pyop2/caching.py", line 583, in wrapper
    value = func(*args, **kwargs)
            ^^^^^^^^^^^^^^^^^^^^^
  File "petsc4py/PETSc/Log.pyx", line 250, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "petsc4py/PETSc/Log.pyx", line 251, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "/home/jhc/codes/fd/fd-main/src/firedrake/pyop2/compilation.py", line 464, in load
    so_name = make_so(compiler_instance, code, extension, comm)
              ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "petsc4py/PETSc/Log.pyx", line 250, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "petsc4py/PETSc/Log.pyx", line 251, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "/home/jhc/codes/fd/fd-main/src/firedrake/pyop2/caching.py", line 583, in wrapper
    value = func(*args, **kwargs)
            ^^^^^^^^^^^^^^^^^^^^^
  File "petsc4py/PETSc/Log.pyx", line 250, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "petsc4py/PETSc/Log.pyx", line 251, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "/home/jhc/codes/fd/fd-main/src/firedrake/pyop2/compilation.py", line 622, in make_so
    return mpi.safe_noncollective(ccomm, compile_single_rank, root=0)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/jhc/codes/fd/fd-main/src/firedrake/pyop2/mpi.py", line 602, in safe_noncollective
    raise result
pyop2.exceptions.CompilationError: 
Command "('mpicc', '-fPIC', '-Wall', '-std=gnu11', '-I/home/jhc/codes/libs/petsc/include', '-I/home/jhc/codes/libs/petsc/arch-fd-main/include', '-I/home/jhc/codes/fd/fd-main/src/firedrake/pyop2', '-march=native', '-O3', '-ffast-math', '-o', '/tmp/pyop2-tempcache-uid1000/0c3e55/b5/tmp52cuhkcx.so', '/tmp/pyop2-tempcache-uid1000/0c3e55/b5/tmp52cuhkcx.c', '-shared', '-L/home/jhc/codes/libs/petsc/lib', '-L/home/jhc/codes/libs/petsc/arch-fd-main/lib', '-Wl,-rpath,/home/jhc/codes/libs/petsc/lib', '-Wl,-rpath,/home/jhc/codes/libs/petsc/arch-fd-main/lib', '-lpetsc', '-lm')" return error status 1.
Unable to compile code

                    Code is:
                    #include <complex.h>
#include <math.h>
#include <petsc.h>
#include <petsc.h>
        #include <math.h>
#include <stdio.h>
#include <petsc.h>


void to_reference_coords_newton_step(double const *__restrict__ C, double const *__restrict__ x0, double const *__restrict__ X, double *__restrict__ dX)
{
  double t0;
  double t1;
  double t10;
  double t11;
  double t12;
  double t13;
  double t14;
  double t15;
  double t16;
  double t17;
  double t18;
  double t19;
  double t2;
  double t20;
  double t21;
  double t22;
  double t23;
  double t24;
  double t25;
  double t26;
  double t27;
  double t3;
  double t4;
  double t5;
  double t6;
  double t7;
  double t8;
  double t9;

  t0 = -1.0 * C[0];
  t1 = t0 + C[3] + 0.0 * C[6];
  t2 = -1.0 * C[1];
  t3 = t2 + C[4] + 0.0 * C[7];
  t4 = -1.0 * C[2];
  t5 = t4 + C[5] + 0.0 * C[8];
  t6 = t1 * t1 + t3 * t3 + t5 * t5;
  t7 = t0 + 0.0 * C[3] + C[6];
  t8 = t2 + 0.0 * C[4] + C[7];
  t9 = t4 + 0.0 * C[5] + C[8];
  t10 = t7 * t7 + t8 * t8 + t9 * t9;
  t11 = t1 * t7 + t3 * t8 + t5 * t9;
  t12 = t7 * t1 + t8 * t3 + t9 * t5;
  t13 = t6 * t10 + -1.0 * t11 * t12;
  t14 = t10 / t13;
  t15 = (-1.0 * t11) / t13;
  t16 = X[1] * 2.0 + -1.0;
  t17 = 0.5 * (t16 + -1.0);
  t18 = (-0.5 * (X[0] * 2.0 + -1.0 + t17 + 1.0) + -1.0 * -0.5 * t17) * -1.0 * 1.224744871391589 * 1.118033988749895;
  t19 = (-0.5 * t16 + -0.5) * -1.224744871391589 * 1.118033988749895;
  t20 = 0.7302967433402214 * (1.3693063937629153 + -1.0 * t18 + -1.0 * t19);
  t21 = 0.7302967433402214 * t18;
  t22 = 0.7302967433402214 * t19;
  t23 = t20 * C[0] + t21 * C[3] + t22 * C[6] + -1.0 * x0[0];
  t24 = t20 * C[1] + t21 * C[4] + t22 * C[7] + -1.0 * x0[1];
  t25 = t20 * C[2] + t21 * C[5] + t22 * C[8] + -1.0 * x0[2];
  dX[0] = dX[0] + (t14 * t1 + t15 * t7) * t23 + (t14 * t3 + t15 * t8) * t24 + (t14 * t5 + t15 * t9) * t25;
  t26 = (-1.0 * t12) / t13;
  t27 = t6 / t13;
  dX[1] = dX[1] + (t26 * t1 + t27 * t7) * t23 + (t26 * t3 + t27 * t8) * t24 + (t26 * t5 + t27 * t9) * t25;

}

static inline void to_reference_coords_kernel(PetscScalar *X, const PetscScalar *x0, const PetscScalar *C)
{
    const int space_dim = 3;

    /*
     * Mapping coordinates from physical to reference space
     */

X[0] = 0.3333333333333333;
X[1] = 0.3333333333333333;

    int converged = 0;
    for (int it = 0; !converged && it < 1; it++) {
        double dX[2] = { 0.0 };
        to_reference_coords_newton_step(C, x0, X, dX);

        if (PetscRealPart(dX[0])*PetscRealPart(dX[0]) + PetscRealPart(dX[1])*PetscRealPart(dX[1]) < 1e-12 * 1e-12) {
            converged = 1;
        }

	X[0] -= dX[0];
	X[1] -= dX[1];
    }
}
        #include <stdbool.h>
#include <stdint.h>
#include <stdbool.h>
#include <math.h>

void pyop2_kernel_evaluate(double *__restrict__ A, int32_t const *__restrict__ cell_orientations_0, double const *__restrict__ w_0, double const *__restrict__ w_1, double const *__restrict__ rt_X)
{
  double t0;
  double t1;
  double t10;
  double t11;
  double t12;
  double t13;
  double t14;
  double t15;
  double t16;
  double t17;
  double t18;
  double t19;
  double t2;
  double t20;
  double t21;
  double t22[1] = { 1.0 };
  double t23;
  double t24[3];
  double t25;
  double t3;
  double t4[6];
  double t5;
  double t6[6];
  double t7;
  double t8;
  double t9;

  t8 = (double) (0.0);
  t0 = rt_X[1] * 2.0 + -1.0;
  t1 = (1.5 * t0 + 0.5) * 0.7071067811865476 * 1.4142135623730951;
  t2 = 0.4714045207910316 * t1;
  t3 = (rt_X[0] * 2.0 + -1.0 + 0.5 * (t0 + -1.0) + 1.0) * 0.7071067811865476 * 1.7320508075688772 * 1.4142135623730951;
  t4[0] = 0.4714045207910317 + t2;
  t4[1] = 0.8164965809277261 + 0.8164965809277261 * t1;
  t4[2] = -0.4714045207910317 + -0.4714045207910316 * t1;
  t4[3] = -0.8164965809277261 + -0.8164965809277261 * t1;
  t4[4] = -0.9428090415820636 + t2;
  t4[5] = -1.4142135623730951 * t3;
  t7 = (double) (0.0);
  t5 = 0.7071067811865476 * t3;
  t6[0] = 0.4714045207910318 + 0.408248290463863 * t3 + -0.2357022603955158 * t1;
  t6[1] = -0.8164965809277261 + -0.7071067811865476 * t3 + 0.40824829046386313 * t1;
  t6[2] = 0.9428090415820634 + -0.408248290463863 * t3 + 0.2357022603955158 * t1;
  t6[3] = t5 + 1.2247448713915892 * t1;
  t6[4] = 0.4714045207910318 + 0.40824829046386296 * t3 + -0.23570226039551584 * t1;
  t6[5] = 0.8164965809277261 + t5 + -0.4082482904638631 * t1;
  for (int32_t i = 0; i <= 5; ++i)
  {
    t7 = t7 + t6[i] * w_1[i];
    t8 = t8 + t4[i] * w_1[i];
  }
  t9 = -1.0 * w_0[1];
  t10 = t9 + w_0[4];
  t11 = -1.0 * w_0[2];
  t12 = t11 + w_0[8];
  t13 = t9 + w_0[7];
  t14 = t11 + w_0[5];
  t15 = t10 * t12 + -1.0 * t13 * t14;
  t16 = -1.0 * w_0[0];
  t17 = t16 + w_0[6];
  t18 = t16 + w_0[3];
  t19 = t17 * t14 + -1.0 * t18 * t12;
  t20 = t18 * t13 + -1.0 * t17 * t10;
  t21 = 1.0 / (((cell_orientations_0[0] == 1.0) ? -1.0 : ((cell_orientations_0[0] == 0.0f) ? 1.0 : (double) (NAN))) * sqrt(t15 * t15 + t19 * t19 + t20 * t20));
  for (int32_t p1 = 0; p1 <= 2; ++p1)
  {
    t23 = -1.0 * w_0[p1];
    t24[p1] = (t23 + w_0[3 + p1]) * t7 + (t23 + w_0[6 + p1]) * t8;
  }
  {
    int32_t const p0 = 0;

    t25 = t22[0] * t21;
    for (int32_t p1_0 = 0; p1_0 <= 2; ++p1_0)
      A[p1_0] = A[p1_0] + t24[p1_0] * t25;
  }

}
        __attribute__((noinline)) /* Clang bug */
        static void pyop2_kernel_prolong(PetscScalar *R, PetscScalar *f, const PetscScalar *X, const PetscScalar *Xc)
        {
            PetscScalar Xref[2];
            int cell = -1;
            int bestcell = -1;
            double bestdist = 1e10;
            for (int i = 0; i < 4; i++) {
                const PetscScalar *Xci = Xc + i*9;
                double celldist = 2*bestdist;
                to_reference_coords_kernel(Xref, X, Xci);
                if (0.5*fabs(1.0*fabs(PetscRealPart(Xref[0])) + 1.0*fabs(PetscRealPart(Xref[1])) + fabs(1.0*PetscRealPart(Xref[0]) + 1.0*PetscRealPart(Xref[1]) - 1.0) - 1.0) <= 1.0e-8) {
                    cell = i;
                    break;
                }

                celldist = 0.5*fabs(1.0*fabs(PetscRealPart(Xref[0])) + 1.0*fabs(PetscRealPart(Xref[1])) + fabs(1.0*PetscRealPart(Xref[0]) + 1.0*PetscRealPart(Xref[1]) - 1.0) - 1.0);
                if (celldist < bestdist) {
                    bestdist = celldist;
                    bestcell = i;
                }

            }
            if (cell == -1) {
                /* We didn't find a cell that contained this point exactly.
                   Did we find one that was close enough? */
                if (bestdist < 10) {
                    cell = bestcell;
                } else {
                    fprintf(stderr, "Could not identify cell in transfer operator. Point: ");
                    for (int coord = 0; coord < 2; coord++) {
                      fprintf(stderr, "%.14e ", X[coord]);
                    }
                    fprintf(stderr, "\n");
                    fprintf(stderr, "Number of candidates: %d. Best distance located: %14e", 4, bestdist);
                    abort();
                }
            }
            const PetscScalar *fi = f + cell*6;
            const PetscScalar *Xci = Xc + cell*9;
            for ( int i = 0; i < 3; i++ ) {
                R[i] = 0;
            }
            pyop2_kernel_evaluate(R, Xci, fi, Xref);
        }

#include <stdint.h>
#include <stdbool.h>

void wrap_pyop2_kernel_prolong(int32_t const start, int32_t const end, double *__restrict__ dat0, double const *__restrict__ dat2, double const *__restrict__ dat1, double const *__restrict__ dat3, int32_t const *__restrict__ map0, int32_t const *__restrict__ map1)
{
  double t0[24];
  double t1[12 * 3];

  for (int32_t n = start; n <= -1 + end; ++n)
  {
    {
      int32_t const i10 = 0;

      {
        int32_t const i8 = 0;

        for (int32_t i9 = 0; i9 <= 23; ++i9)
          t0[i9] = dat2[map0[24 * n + i9]];
      }
    }
    {
      int32_t const i11 = 0;

      for (int32_t i12 = 0; i12 <= 11; ++i12)
        for (int32_t i13 = 0; i13 <= 2; ++i13)
          t1[3 * i12 + i13] = dat3[3 * map1[12 * n + i12] + i13];
    }
    pyop2_kernel_prolong(&(dat0[3 * n]), &(t0[0]), &(dat1[3 * n]), &(t1[0]));
  }
}

                        Compiler output is:
                        /tmp/pyop2-tempcache-uid1000/0c3e55/b5/tmp52cuhkcx.c: In function ‘to_reference_coords_kernel’:
/tmp/pyop2-tempcache-uid1000/0c3e55/b5/tmp52cuhkcx.c:76:15: warning: unused variable ‘space_dim’ [-Wunused-variable]
   76 |     const int space_dim = 3;
      |               ^~~~~~~~~
/tmp/pyop2-tempcache-uid1000/0c3e55/b5/tmp52cuhkcx.c: In function ‘pyop2_kernel_evaluate’:
/tmp/pyop2-tempcache-uid1000/0c3e55/b5/tmp52cuhkcx.c:175:19: warning: unused variable ‘p0’ [-Wunused-variable]
  175 |     int32_t const p0 = 0;
      |                   ^~
/tmp/pyop2-tempcache-uid1000/0c3e55/b5/tmp52cuhkcx.c: In function ‘pyop2_kernel_prolong’:
/tmp/pyop2-tempcache-uid1000/0c3e55/b5/tmp52cuhkcx.c:226:38: warning: passing argument 2 of ‘pyop2_kernel_evaluate’ from incompatible pointer type [-Wincompatible-pointer-types]
  226 |             pyop2_kernel_evaluate(R, Xci, fi, Xref);
      |                                      ^~~
      |                                      |
      |                                      const PetscScalar * {aka const double *}
/tmp/pyop2-tempcache-uid1000/0c3e55/b5/tmp52cuhkcx.c:103:80: note: expected ‘const int32_t * restrict’ {aka ‘const int * restrict’} but argument is of type ‘const PetscScalar *’ {aka ‘const double *’}
  103 | void pyop2_kernel_evaluate(double *__restrict__ A, int32_t const *__restrict__ cell_orientations_0, double const *__restrict__ w_0, double const *__restrict__ w_1, double const *__restrict__ rt_X)
      |                                                    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~
/tmp/pyop2-tempcache-uid1000/0c3e55/b5/tmp52cuhkcx.c:226:13: error: too few arguments to function ‘pyop2_kernel_evaluate’
  226 |             pyop2_kernel_evaluate(R, Xci, fi, Xref);
      |             ^~~~~~~~~~~~~~~~~~~~~
/tmp/pyop2-tempcache-uid1000/0c3e55/b5/tmp52cuhkcx.c:103:6: note: declared here
  103 | void pyop2_kernel_evaluate(double *__restrict__ A, int32_t const *__restrict__ cell_orientations_0, double const *__restrict__ w_0, double const *__restrict__ w_1, double const *__restrict__ rt_X)
      |      ^~~~~~~~~~~~~~~~~~~~~
/tmp/pyop2-tempcache-uid1000/0c3e55/b5/tmp52cuhkcx.c: In function ‘wrap_pyop2_kernel_prolong’:
/tmp/pyop2-tempcache-uid1000/0c3e55/b5/tmp52cuhkcx.c:243:23: warning: unused variable ‘i8’ [-Wunused-variable]
  243 |         int32_t const i8 = 0;
      |                       ^~
/tmp/pyop2-tempcache-uid1000/0c3e55/b5/tmp52cuhkcx.c:240:21: warning: unused variable ‘i10’ [-Wunused-variable]
  240 |       int32_t const i10 = 0;
      |                     ^~~
/tmp/pyop2-tempcache-uid1000/0c3e55/b5/tmp52cuhkcx.c:250:21: warning: unused variable ‘i11’ [-Wunused-variable]
  250 |       int32_t const i11 = 0;
      |                     ^~~

[JHopeCollins]

@pbrubeck any idea?

[connorjward]

I think it's a failure to pass cell orientations through. Very very similar to #5043

[JHopeCollins]

Very very similar to #5043

I thought it might be related, but I still get the error using that branch

[connorjward]

Yeah it's a different bonus kernel parameter that we're not passing