`Meshes.containing_element` disagrees with `Fields.coordinate_field` with MPI

[Sbozzolo]

I run the example below with 2 processes (export CLIMACOMMS_CONTEXT="MPI"; mpirun -np 2 julia --project myexample.lj). The code does the following:

• First, I extract the coordinate of a distributed 2D Field into the hcoords array. hcoords will be different for the
• Second, I loop over each of them and extract those that Meshes/Topologies claim to have a PID that is the same as the PID for the given process. Given that we extracted the points from coordinate_field, we should expect that every point has the same PID as the current process
• Third, I take the difference between the points obtained in the first and the second step. I show that the difference is not empty.

containing_pid is a simple function that just reads of the underlying topology.

function containing_pid(
    target_point::P,
    topology::T,
) where {
    P <: Union{Geometry.LatLongPoint, Geometry.XYPoint},
    T <: Topologies.Topology2D,
}
    containing_elem = Meshes.containing_element(topology.mesh, target_point)
    gidx = topology.orderindex[containing_elem]
    return topology.elempid[gidx]
end

containing_element, however, is not a very simple function:

function containing_element(
    mesh::AbstractCubedSphere,
    coord::Geometry.Cartesian123Point,
)
    ne = mesh.ne
    panel = containing_panel(coord)
    ϕx, ϕy = _inv_coordinates(mesh, coord, panel)

    x = refindex(ϕx, ne)
    y = refindex(ϕy, ne)
    return CartesianIndex(x, y, panel)
end
containing_element(mesh::AbstractCubedSphere, coord::Geometry.LatLongPoint) =
    containing_element(
        mesh,
        Geometry.Cartesian123Point(
            coord,
            Geometry.SphericalGlobalGeometry(mesh.domain.radius),
        ),
    )

So my first instintct is that there's is something incorrect here (maybe the boundary points are not handled correctly?).

Code to reproduce:

import ClimaComms
ClimaComms.@import_required_backends
using ClimaCore.CommonSpaces
import ClimaCore
import ClimaCore.Remapping: containing_pid

comms_ctx = ClimaComms.context()
ClimaComms.init(comms_ctx)

space = CubedSphereSpace(;
    radius = 10,
    n_quad_points = 4,
    h_elem = 2,
)

# These are coordinates owned by this process according to `coordinate_field`
coords = ClimaCore.to_cpu(ClimaCore.Fields.coordinate_field(space))
lats = ClimaCore.Fields.field2array(coords.lat)
lons = ClimaCore.Fields.field2array(coords.long)
hcoords = ClimaCore.Geometry.LatLongPoint.(lats, lons)

# Let's see if `containing_pid`, which uses `Meshes.containing_element`, agrees
hcoords_containing_element = filter(x -> containing_pid(x, space.grid.topology) == ClimaComms.mypid(comms_ctx), hcoords)

@show ClimaComms.mypid(comms_ctx), setdiff(Set(hcoords), Set(hcoords_containing_element))

This returns

(ClimaComms.mypid(comms_ctx), setdiff(Set(hcoords), Set(hcoords_containing_element))) = (1, Set(ClimaCore.Geometry.LatLongPoint{Float64}[LatLongPoint(-35.264389682754654, 135.0), LatLongPoint(0.0, 135.0), LatLongPoint(-9.542300130845922, 135.0), LatLongPoint(-25.72208955190873, 135.0), LatLongPoint(44.21261139871027, -166.62721893004286), LatLongPoint(35.264389682754654, 135.0), LatLongPoint(44.21261139871027, 166.62721893004286), LatLongPoint(-44.21261139871027, 103.37278106995717), LatLongPoint(25.72208955190873, 135.0), LatLongPoint(44.21261139871027, -76.62721893004284), LatLongPoint(9.542300130845922, 135.0), LatLongPoint(-0.0, 135.0)]))
(ClimaComms.mypid(comms_ctx), setdiff(Set(hcoords), Set(hcoords_containing_element))) = (2, Set(ClimaCore.Geometry.LatLongPoint{Float64}[LatLongPoint(-44.21261139871027, 13.372781069957151), LatLongPoint(-44.21261139871027, 76.62721893004284), LatLongPoint(-9.542300130845922, -45.0), LatLongPoint(44.21261139871027, -103.37278106995717), LatLongPoint(9.542300130845922, -45.0), LatLongPoint(-44.21261139871027, -13.372781069957151)]))

Note that just a handful of points are incorrect (length(hcoords) = 4800), so maybe this is a roundoff error?