Account for spherical face bulging during hashgrid construction and point in cell checks

[wyatt-fluidnumerics]

Description

Issue #2878 raised awareness of two separate bugs:

1. During point in cell checking on spherical unstructured or cuvilinear grids the old code did a projection of the particle position by calculating the direction normal to the face and moving the particle in that direction. This meant that for particles close to the edge of a face the normal projection could project them outside of the cell face, which was incorrect. This could result in a GridSearchingError.
2. The spatialhash grid bounds were previously calculated only by using the positions of a cell's corners. This does not account for the bulged component of the cell that the use a spherical mesh creates. As a result, the spatialhash grid does not fully enclose spherical grids.

Both of these issues are magnified for grids that contain large cells, which is why for most of the small celled grids that Parcels is tested on, they did not show up earlier.

The fixes are the following:

1. Change from a projection normal to the face to a gnomonic projection which moves the particle on a ray that extends through the face and the center of the unit sphere. This way particles very close to the edge of their cell remain in the cell after projection.
2. Add a much more extensive bounding box calculation for the spatialhash grid. For any give cell their are three possibilities for its minima/maxima. The first is that is on one of the nodes (this is all that was used for the calculation previously). The second is that the face overlaps a global minama/maxima of the unit sphere, think a face that extends over the North Pole, in this case the minima/maxima can be anywhere on the inside of the face, for this reason we introduce and explicit check to see whether a face contains one of the 6 global minima/maxima of the unit sphere, and if it does set the respective xyz bound accordingly. Third, the minima maxima could sits on the face edge. In this case picture a face that has a bottom edge extending from lon -5 to 5 at a lat of 10, the final node forming the face could sit at lon 0, lat 20. In this case the point lon 0, lat 10 serves as an extrema for the face due to its bulge between longitudes of -5 and 5. However this extrema lies only on an edge rather than on a node, or on a global extrema. These types of points can be calculated analytically by parameterizing the curve between to vertices as a cos function and then determining whether that function has a local extrema along the curve. If it does, then that point is one of the faces extrema, if it doesn't, then the node at the end of the curve is the extrema. I am happy to help put this calculation into the docs that @fluidnumericsJoe is writing up.

The PR in this state covers the case for unstrucutred grids. Structured curvilinear grids still suffer the same bug and I have yet to implement these fixes there. I'll try and do that tomorrow so this is resolved before the v4 release!

Checklist

• Closes Fix tutorial_nested_grids.ipynb - grid searching #2878
• Tests added
• This PR targets the correct branch (main for normal development, v3-support for v3 support)

AI Disclosure

• This PR contains AI-generated content.
  • I have tested any AI-generated content in my PR.
  • I take responsibility for any AI-generated content in my PR.
  • Describe how you used it (e.g., by pasting your prompt): I guided claude code to generate parts of the new _spherical_triangle_bounds function. I made significant modifications to the output for clarity, structure, and scope. I also guided claude code to modify and improve the new testing associated with the PR. I have thoroughly reviewed and understand all code associated with the PR.

Currently a draft PR until the spherical curvilinear fixes are in.

[erikvansebille]

Looks good; one comment and one question below

[erikvansebille]

Is this the most efficient implementation? Does using np.stack here (and later in this PR) not lead to a similar performance loss as #2846? Or are the arrays so small here that that doesn't matter?

[fluidnumericsJoe]

Hmm, yeah, maybe preallocating face_matrix and concatenating into this array would be better? I wonder if concatenating in Python concats by reference or by value

[fluidnumericsJoe]

Looks like we'd have to use the out= argument with concatenate (or stack) to avoid additional allocation costs, see https://numpy.org/devdocs/reference/generated/numpy.concatenate.html and https://numpy.org/doc/stable/reference/generated/numpy.stack.html

Source code for stack shows it just calls concatenate under the hood .

Neither avoids the copy cost here though. Buffers are copied . Only numpy only route I see possible here is to use an array view to avoid copies.

[erikvansebille]

Thanks for digging into this, @fluidnumericsJoe. I don't think we need to go too deep into optimisation of this statement specifically before we know whether it is a bottle-neck or not. So I'd avoid creating our own views. I just raised the question in case there was a low-hanging-fruit solution to make this faster

[fluidnumericsJoe]

Pre-allocating and using np.stack with the out parameter could be helpful - but that pre-allocation might need to have broader scope than just a single PIC call to be of any value now that I think of it. IMO, after thinking on this a bit, I wouldn't worry about this too much here since its a single call per vectorized PIC check.

[erikvansebille]

I updated the nested grids tutorial so that gets run again - hope that's ok?

[wyatt-fluidnumerics]

That's great, thank you!

[wyatt-fluidnumerics]

@erikvansebille @fluidnumericsJoe I've added the changes for spherical curvilinear grids. I generalized the _spherical_face_bounds function to work for both cases, and then adjusted _spherical_project_cell_and_query to use a gnomonic projection instead of orthonormal. I also changed some of the new tests to use xgrid/uxgrid parameterizations. I think things should be good to go now, but would be great to get another review on the changes!