`Meshes.containing_element` disagrees with `Fields.coordinate_field` with MPI
- Dominant language
- Julia
- Stars
- 117
- Forks
- 19
- Avg merge
- 3d 8h
- Merged PRs (30d)
- 38
Description
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.
```julia
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:
```julia
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:
```julia
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?
Contributor guide
No contributing guide indexed for this repository
Research direction
Run the provided two-process MPI reproduction first, then inspect Fields.coordinate_field alongside Meshes.containing_element and the containing_pid helper. Trace containing_panel, _inv_coordinates, refindex, and the topology PID lookup, focusing on boundary coordinates and roundoff. Done means the PIDs returned for extracted coordinates agree without the reported mismatches.
Written by the indexing model from the issue text.
Assessment
- Tech stack
- julia
- Domain
- distributed-systems
- Issue type
- Bug
- Difficulty
- 4/5
- Estimated time
- 3-5 days
- Activity status
- Stale
- Clarity
- Mostly clear
- Newbie friendliness
- 38/100