BUG: Checkpointing Hdiv elements on higher order meshes is not correct.

[JaroslavHron]

When trying to restore checkpointed function, it seems not to be identical for some elements.

Steps to Reproduce

from firedrake import *
from firedrake.__future__ import interpolate

def report(z):
    l2_norm=sqrt(assemble(inner(z, z)*dx))
    h1_norm=sqrt(assemble(inner(grad(z), grad(z))*dx))
    div_norm=sqrt(assemble(inner(div(z), div(z))*dx))
    print(f"{l2_norm=} {h1_norm=} {div_norm=}")

mesh = UnitDiskMesh(refinement_level=2, name='disk')
k=1
V = FunctionSpace(mesh, "RT", k)
#V = VectorFunctionSpace(mesh, "CG", k)
x=SpatialCoordinate(mesh)

z=assemble(interpolate(as_vector((x[1],x[0])), V))
report(z)

with CheckpointFile(f"temp_test.h5", 'w') as chp:
    chp.save_function(z, name='z')
report(z)

with CheckpointFile(f"temp_test.h5", 'r') as chp:
    mesh=chp.load_mesh(name='disk')
    z=chp.load_function(mesh, name='z')
report(z)

Expected behavior When run for CG element the restored function seems to be identical

l2_norm=1.2443484447999211 h1_norm=2.4976354854818696 div_norm=5.739510232624121e-17
l2_norm=1.2443484447999211 h1_norm=2.4976354854818696 div_norm=5.739510232624121e-17
l2_norm=1.2443484447999211 h1_norm=2.4976354854818696 div_norm=5.739510232624121e-17

While for RT/BDM element, there is a small difference:

l2_norm=1.2468831137910228 h1_norm=1.582540644473559e-15 div_norm=2.1481119431386426e-15
l2_norm=1.2468831137910228 h1_norm=1.582540644473559e-15 div_norm=2.1481119431386426e-15
l2_norm=1.2468831137910226 h1_norm=2.8526732776154944e-15 div_norm=4.014406586513413e-15

Is this to be expected?

[ksagiyam]

That is indeed expected. We currently project functions on HDiv/HCurl spaces to ones on appropriate DG spaces before saving and project them back after loading. Those are projection errors.

[JaroslavHron]

Hi, thanks - I see - I was investigating this since we see much bigger difference with higher order mesh. If I do the same with quadratic mesh, the CG vector field is restored from the checkpoint file correctly, while the HDiv is restored with much higer value of the divergence:

from firedrake import *
from firedrake.__future__ import interpolate

def report(z):
    l2_norm=sqrt(assemble(inner(z, z)*dx))
    h1_norm=sqrt(assemble(inner(grad(z), grad(z))*dx))
    div_norm=sqrt(assemble(inner(div(z), div(z))*dx))
    print(f"{l2_norm=} {h1_norm=} {div_norm=}")

mesh = UnitDiskMesh(refinement_level=2, name='disk')

l=2   # mesh order
R=1.0
V1 = mesh.coordinates.function_space()
V2 = VectorFunctionSpace(mesh, "CG", l)
x = SpatialCoordinate(mesh)
f = Function(V2).interpolate(mesh.coordinates)
ur=R*as_vector([x[0]/sqrt(x[0]**2+x[1]**2), x[1]/sqrt(x[0]**2+x[1]**2)])
bc=DirichletBC(V2, ur, 1)
bc.apply(f)  # correct the boundary to conform to circle
mesh = Mesh(f, name='disk')
mesh.init()
print(mesh.coordinates.ufl_element())

k=1
V = FunctionSpace(mesh, "BDM", k)
#V = VectorFunctionSpace(mesh, "CG", k)
x=SpatialCoordinate(mesh)

z=assemble(interpolate(as_vector((x[1],x[0])), V))
report(z)

with CheckpointFile(f"temp_test.h5", 'w') as chp:
    chp.save_function(z, name='z')

report(z)

with CheckpointFile(f"temp_test.h5", 'r') as chp:
    mesh=chp.load_mesh(name='disk')
    z=chp.load_function(mesh, name='z')

report(z)

then i see this:

l2_norm=1.2533075678366274 h1_norm=2.5050750966777287 div_norm=1.2924204487211691e-14
l2_norm=1.2533075678366274 h1_norm=2.5050750966777287 div_norm=1.2924204487211691e-14
l2_norm=1.2538698595522295 h1_norm=2.5088860745287382 div_norm=0.010262263358943974

But maybe I am constructing the higher order mesh in wrong way?

[ksagiyam]

That indeed looks like a bug. Probably, it is caused by that projection onto the DG spaces are not exact if the mappings are non-affine. I will have a closer look when I am back from holiday. Are you only interested in simplex meshes for now?

[JaroslavHron]

Yes, we use simplex meshes for now. I will modified the bug title to reflect the problem. Enjoy your holiday.