Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

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

Open
JaroslavHron opened this issue Sep 25, 2024 · 4 comments · May be fixed by #3838
Open

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

JaroslavHron opened this issue Sep 25, 2024 · 4 comments · May be fixed by #3838
Labels
bug

Comments

@JaroslavHron
Copy link
Contributor

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
Copy link
Contributor

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
Copy link
Contributor Author

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
Copy link
Contributor

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
Copy link
Contributor Author

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

@JaroslavHron JaroslavHron changed the title BUG: Checkpointing some elements gives slightly different results after restoring the checkpoint. BUG: Checkpointing Hdiv elements on higher order meshes is not correct. Sep 27, 2024
@ksagiyam ksagiyam linked a pull request Oct 24, 2024 that will close this issue
Open
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
bug
Projects
None yet
Milestone
No milestone
Development

Successfully merging a pull request may close this issue.

io: fix saving/loading of HDiv/HCurl functions on a high-order mesh firedrakeproject/firedrake
2 participants
@JaroslavHron @ksagiyam