prolong gives incorrect results for higher-order (k>1) KMV elements

[Thingwayy]

I am attempting to perform a convergence test for a problem using KMV elements. Since the exact solution is unknown, I am calculating the convergence rate using the formula

log(E_{h}-E_{h/2} / E_{h/2}-E(h/4) ) / log 2.

This requires transferring a solution from a coarse grid to a finer one. Normally, I would use interpolate or project, but as these are not fully supported for N1curl elements (another variable in my context), I had to resort to using the prolong operator.

The issue arises when I use prolong with KMV elements of degree greater than 1 (i.e., k > 1). The operation executes, but the generated values on the refined mesh appear to be incorrect, which breaks my convergence test.

Has anyone encountered this problem with prolong and high-order KMV elements? Is there a known workaround, a different operator I should use for this transfer, or a configuration I might be missing?

Here is a simple code that demonstrates the problem:

from firedrake import *
from firedrake import MeshHierarchy, TransferManager
for order in [1,2]:
    for n in [4, 8]:
        # tm.prolong(m_n_1, m_n1)
        mesh = UnitSquareMesh(n, n)
        hierarchy = MeshHierarchy(mesh, 1)
        mesh_fined = hierarchy[1]
        V = FunctionSpace(mesh, "KMV", order)
        V_fined = FunctionSpace(mesh_fined, "KMV", order)

        x = SpatialCoordinate(mesh)
        f_expr = sin(x[0])*cos(x[1])  # 符号表达式

        u = Function(V); u_fined = Function(V_fined)
        u.interpolate(f_expr); u_fined.interpolate(f_expr)

        u_1 = Function(V_fined); u_1.interpolate(u)
        tm = TransferManager()
        u_2 = Function(V_fined); tm.prolong(u, u_2)

        dx_fined = dx(domain=mesh_fined)
        # 直接计算与符号表达式的误差
        error_interpolate = sqrt(assemble((u_1 - u_fined) ** 2 * dx_fined))
        error_prolong = sqrt(assemble((u_2 - u_fined) ** 2 * dx_fined))
        print(f"KMV_order={order}\n"
              f"{n=} error_interpolate: {error_interpolate}\n"
              f"{n=} error_prolong: {error_prolong}\n") 

mpirun -n 1 python3 test.py            
/home/thingway/.local/lib/python3.10/site-packages/numpy/_core/getlimits.py:551: UserWarning: Signature b'\x00\xd0\xcc\xcc\xcc\xcc\xcc\xcc\xfb\xbf\x00\x00\x00\x00\x00\x00' for <class 'numpy.longdouble'> does not match any known type: falling back to type probe function.
This warnings indicates broken support for the dtype!
  machar = _get_machar(dtype)
firedrake:WARNING OMP_NUM_THREADS is not set or is set to a value greater than 1, we suggest setting OMP_NUM_THREADS=1 to improve performance
KMV_order=1
n=4 error_interpolate: 0.0030940381005530398
n=4 error_prolong: 0.0030940381005531

KMV_order=1
n=8 error_interpolate: 0.000785567405199398
n=8 error_prolong: 0.0007855674051994479

KMV_order=2
n=4 error_interpolate: 0.00018095600904842547
n=4 error_prolong: 0.044837447325300414

KMV_order=2
n=8 error_interpolate: 2.2609439801423873e-05
n=8 error_prolong: 0.035130566642550536

[pbrubeck]

You mentioned two issues, the first one not being able to prolong an N1curl Function, and the second one getting unexpected results for KMV elements. I will need to have a closer look to the second issue, but for the meantime, you could try to get prolong properly working with N1curl elements .

We support prolong for N1curl elements (or more generally non-Lagrange elements). There are two mechanisms for multigrid prolongation, one is firedrake.prolong() (which is limited to Lagrange elements), and the general one is TransferManager.prolong().

Here's an example on how to prolong an N1curl Function

from firedrake import *

mesh = UnitSquareMesh(5, 5)
mh = MeshHierarchy(mesh, 2)
V0 = FunctionSpace(mh[0], "N1curl", 1)
V1 = V0.reconstruct(mesh=mh[1])

u0 = Function(V0)
u1 = Function(V1)

tm = TransferManager()
tm.prolong(u0, u1)

[Thingwayy]

Thank you for your kind reply. Perhaps I did not describe my problem clearly. The prolong operator does work for the N1curl element. Indeed, I'm using the prolongation operator for the N1curl element because interpolate and project both failed to transfer my coarse-grid solution to the refined space.

I will look into the two prolong mechanisms you mentioned. When performing convergence tests, transferring the KMV element can be done directly using interpolate, it makes correct convergence rate.

The issue with prolonging the KMV element does seem to exist. I look forward to any further insights.

[pbrubeck]

Seems like Lagrange-only approach firedrake.prolong does the right thing for KMV(2) elements.
If I replace

tm = TransferManager()
u_2 = Function(V_fined); tm.prolong(u, u_2)

in your code with

u_2 = Function(V_fined); prolong(u, u_2)

I get

KMV_order=1
n=4 error_interpolate: 0.003094038100553039
n=4 error_prolong: 0.0030940381005530415

KMV_order=1
n=8 error_interpolate: 0.0007855674051994033
n=8 error_prolong: 0.0007855674051994119

KMV_order=2
n=4 error_interpolate: 0.00018095600904842745
n=4 error_prolong: 0.00018095600904843067

KMV_order=2
n=8 error_interpolate: 2.2609439801428745e-05
n=8 error_prolong: 2.2609439801424334e-05

[pbrubeck]

We can change the logic in TransferManager.prolong() to call firedrake.prolong() on KMV elements.

[pbrubeck]

#4747

[Thingwayy]

This seems reasonable. I have one more question: In the context of computing the convergence rate with the aforementioned formula, is the use of the prolong function to transfer a solution from a coarse mesh to a finer mesh appropriate for Nédélec (N1curl) elements?

As I noted earlier, both the interpolate and project operators fail in Firedrake for this purpose.

Is it possible to implement interpolation between different mesh spaces for N1curl (or similar) elements? In FEniCSx, this can be achieved using the interpolate_nonmatching function. Does Firedrake offer a comparable capability, or is there a recommended workaround?

[pbrubeck]

We were not labelling the KMV elements with correct embedding degree, this is fixed here firedrakeproject/fiat#199

The issue is that KMV(p) is a subspace of CG(p+bump), where bump depends on p.

[pbrubeck]

project should be fine on N1curl. There are plans to support cross-mesh interpolation for N1curl elements in the near future.

[Thingwayy]

project should be fine on N1curl.

In my code, the transfer part is

    m_n_1, B_n_1 = mllg(mesh); m_n2, B_n2 = mllg(mesh_fined)
    m_n1 = Function(V_fined); B_n1 = Function(ND_fined)
    m_n1.interpolate(m_n_1)
    tm = TransferManager()
    # tm.prolong(m_n_1, m_n1) 
    # tm.prolong(B_n_1, B_n1)
    B_n1.project(B_n_1, ND_fined) 

where B_n_1 is numerical solution in a coarse Nedelec space, when i run the code with

B_n1.project(B_n_1, ND_fined) 

i get error

 mpiexec -n 1 python3 convergence_without_known_sol.py 
Traceback (most recent call last):
  File "/mnt/c/Users/Thingway/Desktop/code_LLGp2/tri_tet_KMV_Nedelec/convergence_without_known_sol.py", line 178, in <module>
    B_n1.project(B_n_1, ND_fined) 
  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/thingway/.local/lib/python3.10/site-packages/firedrake/adjoint_utils/function.py", line 47, in wrapper
    output = project(self, b, *args, **kwargs)
  File "/home/thingway/.local/lib/python3.10/site-packages/firedrake/function.py", line 354, in project
    return projection.project(b, self, *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/thingway/.local/lib/python3.10/site-packages/firedrake/adjoint_utils/projection.py", line 34, in wrapper
    output = project(*args, **kwargs)
  File "/home/thingway/.local/lib/python3.10/site-packages/firedrake/projection.py", line 105, in project
    val = Projector(
  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/thingway/.local/lib/python3.10/site-packages/firedrake/projection.py", line 359, in Projector
    raise ValueError("Haven't implemented supermesh projection with boundary conditions yet, sorry!")
ValueError: Haven't implemented supermesh projection with boundary conditions yet, sorry!
--------------------------------------------------------------------------
Primary job  terminated normally, but 1 process returned
a non-zero exit code. Per user-direction, the job has been aborted.
--------------------------------------------------------------------------
--------------------------------------------------------------------------
mpiexec detected that one or more processes exited with non-zero status, thus causing
the job to be terminated. The first process to do so was:

  Process name: [[6430,1],0]
  Exit code:    1
--------------------------------------------------------------------------

[pbrubeck]

Ah, that's right. We don't have cross-mesh projection working either. Hopefully we will start supporting it soon. For a MeshHierarchy, the recommended approach is prolong.