search

firedrakeproject / firedrake

Firedrake is an automated system for the portable solution of partial differential equations using the finite element method (FEM)
https://firedrakeproject.org
Other
517 stars 160 forks source link

BUG: When making a mesh hierarchy from a base mesh loaded from file, the fine mesh does not match the coarse mesh #3739

Open pefarrell opened 3 months ago

pefarrell commented 3 months ago

Describe the bug I am loading a mesh from a file and making a mesh hierarchy. I expect the refined mesh to coincide with the coarse mesh. It does not.

Steps to Reproduce

Save and gunzip the following mesh.

mesh_0.h5.gz

Run

from firedrake import *

with CheckpointFile('mesh_0.h5', 'r') as afile:
    mesh_import = afile.load_mesh('naca0012')

VTKFile(f"mesh_initial.pvd").write(mesh_import.coordinates)

mh = MeshHierarchy(mesh_import, 1)

for i in range(2):
    VTKFile(f"mesh_after_{i}.pvd").write(mh[i].coordinates)

I find that mesh_initial.pvd matches perfectly with mesh_after_0.pvd (as expected, but that mesh_after_0.pvd and mesh_after_1.pvd do not coincide geometrically.

mesh-0

Zoom in at the leading edge of the coarse mesh, coloured in red (mesh_after_0.pvd)

mesh-1

The same view of the fine mesh, coloured in blue (mesh_after_1.pvd)

Expected behavior I expected the two meshes to coincide geometrically.

Environment:

Additional Info Add any other context about the problem here.

ksagiyam commented 3 months ago

I think the coordinates of the mesh saved in mesh_0.h5 had been modified, and Firedrake coordinates and plex coordinates of the stored mesh are different. One solution would be to overwrite plex coordinates (from which coordinates of the refined meshes are made) after loading to match with the Firedrake coordinates by doing the reverse of reordered_coords.

ksagiyam commented 2 months ago

Details:

A Firedrake mesh object wraps a PETSc DMPlex object, which contains the topology information and the coordinates. As Firedrake uses different DoF orderings from PETSc, we convert DMPlex coordinates to Firedrake coordinates when the mesh is initialised using reorderd_coords function:https://github.com/firedrakeproject/firedrake/blob/d041208aade32f19904d177e824d160fffb5fb28/firedrake/mesh.py#L1948.

When Firedrake's coordinates are modified, coordinates of the underlying DMPlex are not automatically updated accordingly. (If we then save/load the mesh, both Firedrake coordinates and the unsynchronised DMPlex coordinates are saved/loaded, but checkpointing is not directly relevant here.) This causes a problem as MeshHierarchy uses a DMPlex function to do refinement, and this DMPlex function will not know that coordinates have been modified. So, when you modify Firedrake coordinates, you might also want to modify DMPlex coordinates accordingly. This can be done (if you do not have high order coordinates) by doing the reverse of reordered_coords as mentioned in the above.

connorjward commented 2 months ago

What do we need to do to start natively using the DMPlex coordinates? Is it not simply that we need to build a new cell-node map that has an extra permutation? In the same way that we do hex meshes?

I'm guessing that there might be issues to do with higher-order coordinates - does DMPlex support them?

ksagiyam commented 2 months ago

I think DMPlex supports higher-order coordinates.

We first need to have plex.coordSection to match up with our CG section; we need to push ghost DoFs to the end. Then DoF layout on each entity must match (given plex cone of that entity, we lay out DoFs using FIAT definition). I think we could work around the latter by introducing additional permutation in our current cell_node_map, but note that we use the same cell_node_map for non-coordinate CG function spaces.

connorjward commented 2 months ago

We first need to have plex.coordSection to match up with our CG section; we need to push ghost DoFs to the end.

This doesn't strike me as being too difficult. We already build Sections like this.

Then DoF layout on each entity must match (given plex cone of that entity, we lay out DoFs using FIAT definition). I think we could work around the latter by introducing additional permutation in our current cell_node_map, but note that we use the same cell_node_map for non-coordinate CG function spaces.

Indeed, I was thinking we could represent this as an alternative variant of a normal FunctionSpace in some way, such that the maps are permuted.