Parallel Compatible Cross-Mesh Interpolation

[ReubenHill]

Description

This implements cross-mesh interpolation for interpolating into function spaces which have point evaluation nodes. This is an attempt to integrate this gist into Firedrake and makes use of the recent merging of #2936.

We use VertexOnlyMesh as an intermediary for the global locations of the point evaluation nodes of the target function space:

1. We get the point evaluation node locations using the method described in the interpolation from external data section of the manual. This will have the parallel domain decomposition of the target mesh.
2. Next we create a VertexOnlyMesh (A) at those locations within the source mesh such that we inherit the source mesh's parallel domain decomposition.
3. We interpolate our expression in our source function space onto a P0DG function space on VertexOnlyMesh (A), which has the effect of point evaluating at the target function space node locations.
4. This VertexOnlyMesh (A) has an input_ordering VertexOnlyMesh (B) whose vertices have the ordering and parallel domain decomposition of the target function space global node locations. We interpolate from P0DG on (A) onto P0DG on (B). Under the hood, this is an SF reduce operation which moves the point evaluations from (A) to (B).
5. We now have a Function on the input_ordering VertexOnlyMesh (B) which has point evaluations from our source mesh function space at the target mesh function space node locations. These are in the correct order and have the correct domain decomposition. We can therefore safely set the dat of a function in our target function space to the values of this function.

For this to work for the general case, we would need to create a VertexOnlyMesh at the global quadrature points of the target function space, which is rather more complicated than the work I've done here.

Some important notes:

• This does not require one mesh to be a structured refinement of the other. This, for example, should allow you to solve a PDE on two different unstructured meshes, one of which is finer than the other, and directly check the difference in solutions by interpolating from one mesh to the other within Firedrake.
• The target mesh merely needs to share some of its domain with the source mesh and we should get the results we expect (see to do note below).
• Crucially, this is entirely parallel compatible!

Other notes:

• The VertexOnlyMesh required is stored in the Interpolator, as are the underlying Interpolators
• For interpolation between mixed spaces, I create sub_interpolators for each space and evaluate them as necessary when calling interpolate
  • Interpolation into a mixed space therefore requires the function space being interpolated from to be another mixed space. I throw a NotImplementedError if not.

To do:

• Test with more meshes
  • This ought to work with meshes of different geometric dimensions as long as the target mesh is fully within the source mesh. I should check this, because this would be really cool! You need the geometric dimensions to match or you would need to manually specify where the lower dimensional mesh was within the higher dimensional one! This is the whole point of immersed meshes and there's no reason why an immersed mesh wouldn't work as expected methinks.
  • Check that meshes with different domains work as expected.
    • Decide whether to optionally allow interpolation into a mesh with extent beyond the target mesh - I think this is possible with how I've implemented VertexOnlyMesh's input_ordering property, so long as I switch off giving an error when not all the points are found, and would result in a Function with some values that are unchanged by the interpolation operation. Implemented!
• Work out how to ensure we only try to interpolate into point evaluation function spaces
• Test the adjoint - since this uses pure Firedrake, I think it ought to work fine. But perhaps setting the dat values will cause an issue? EDIT: Testing adjoint looks to be too tricky
• Implement tests for exact interpolation and transpose interpolation by using meshes which are refinements of each other
• Document in the manual

Associated Pull Requests:

None

Fixes Issues:

None, this is a new feature.

Checklist for author:

• I have linted the codebase (make lint in the firedrake source directory).
• My changes generate no new warnings.
• All of my functions and classes have appropriate docstrings.
• I have commented my code where its purpose may be unclear.
• I have included and updated any relevant documentation.
• Documentation builds locally (make linkcheck; make html; make latexpdf in firedrake/docs directory)
• I have added tests specific to the issues fixed in this PR.
• I have added tests that exercise the new functionality I have introduced
• Tests pass locally (pytest tests in the firedrake source directory) (useful, but not essential if you don't have suitable hardware).
• I have performed a self-review of my own code using the below guidelines.

Checklist for reviewer:

• Docstrings present
• New tests present
• Code correctly commented
• No bad "code smells"
• No issues in parallel
• No CI issues (excessive parallelism/memory usage/time/warnings generated)
• Upstream/dependent branches and PRs are ready

[connorjward]

Exciting!

[stephankramer]

Very exciting! Quick question out of interest: will this be adjointable? I'm guessing step 3 already is, step 5 should be trivial to adjoint - but step 4 requires reverse communication, is that also easy?

[dham]

Very exciting! Quick question out of interest: will this be adjointable? I'm guessing step 3 already is, step 5 should be trivial to adjoint - but step 4 requires reverse communication, is that also easy?

Happily PETSc SF gives you the reverse communication so this stage is adjointed (already). You are correct that this means that the mesh to mesh interpolation is adjoinable.

[ReubenHill]

Very exciting! Quick question out of interest: will this be adjointable? I'm guessing step 3 already is, step 5 should be trivial to adjoint - but step 4 requires reverse communication, is that also easy?

Once I have correctly implemented the transpose interpolation operation, this will be adjointable out of the box since interpolation is already annotated by making use of the transpose (i.e. adjoint) interpolation operation. What use cases do you have in mind?

[stephankramer]

That's fantastic! For use cases: any time-dependent models with dynamic mesh adaptivity where I want an adjoint, say for sensitivity analysis. For optimisation and goal-based adaptivity it gets a bit more involved as you have to manage the sequence of meshes in the outer loop being adapted which requires re-taping (e.g. https://github.com/pyroteus/pyroteus)

Stale review

[JDBetteridge]

Previously approved