What is the recommended PETSc way to project a field from one DMPlex to another (same bounding box)?

[gthyagi]

Hi @knepley,

I'm trying to project a field from one DMPlex to another using PETSc within Underworld3. Both meshes cover the same bounding box, but have different resolutions.

Here's a minimal example:

import underworld3 as uw

# Mesh bounding box
minX, minY, minZ = 0., 0., -40.
maxX, maxY, maxZ = 150., 150., 0.

# Create first mesh (coarser)
mesh_1 = uw.meshing.UnstructuredSimplexBox(minCoords=(minX, minY, minZ), 
                                           maxCoords=(maxX, maxY, maxZ),
                                           cellSize=16)

v_soln_1 = uw.discretisation.MeshVariable('V', mesh_1, mesh_1.data.shape[1], degree=2)

# Assign values
with mesh_1.access(v_soln_1):
    v_soln_1.data[...] = 1.

# Create second mesh (finer)
mesh_2 = uw.meshing.UnstructuredSimplexBox(minCoords=(minX, minY, minZ), 
                                           maxCoords=(maxX, maxY, maxZ),
                                           cellSize=10)

v_soln_2 = uw.discretisation.MeshVariable('V', mesh_2, mesh_2.data.shape[1], degree=2)

# Question: How do I project v_soln_1 to v_soln_2?

What is the recommended PETSc approach to project v_soln_1 onto mesh_2's v_soln_2, assuming the meshes have the same bounding box?

Thanks!
Thyagi

[knepley]

Beyond the bounding box, the two really important things are whether the meshes are nested and what you want to preserve. If the meshes are not nested, we have to do large scale point location. This can be done by locating dual basis vectors on the other mesh (Plex does this), using a particle cloud (Underworld does this, and DMSwarm can too now, and Firedrake), or using a common refinement of the mesh (supermesh) (there is a package and Firedrake used to do it this way). The supermesh way can do accurate integrals, so it can preserve moments. The point cloud method can also preserve moments with solves (this is what DMSwarm does). If you do not care about preserving moments, just interpolation using the dual basis or particles will work.

I wrote the above since there are a lot of choices. The Plex default is to create an interpolation matrix from one mesh to the other using DMPlexCreateInterpolator(), but as you see, that is not the only thing you can do.

[gthyagi]

@knepley Thanks for the detailed explanation.
Is there a PETSc example that demonstrates using DMSwarm with a particle cloud for interpolation?

[gthyagi]

@lmoresi Currently, we are using rbf method to interpolate. Does this same as DMSwarm with a particle cloud?

[lmoresi]

Not really - the function.evaluate call will use the DMInterpolate routines. In the updated code, I added an rbf fallback for points that fall outside the mesh. This is important for curved meshes and semi-Lagrange advection where it is easy to specify points lying just outside the boundary. This is now handled behind the scenes.

There is a pure rbf option that can be used for visualisation etc when exact values are less important than speed.

The new code also has a global evaluation that handles non-local points and this also uses interpolation paired with rbf.

L