Normalize the gradient by the dual cell and add the spherical metric terms to curl and divergence - #1663
Open
rajeeja wants to merge 4 commits into
Open
Normalize the gradient by the dual cell and add the spherical metric terms to curl and divergence#1663rajeeja wants to merge 4 commits into
rajeeja wants to merge 4 commits into
Conversation
…terms Green-Gauss on a contour must divide by the area that contour encloses. The face gradient integrates around the dual cell but divided by the primal face area, inflating results by A_dual/A_primal: about 4 on quadrilateral meshes and 3 on hexagonal ones. Refining healpix z3 to z6 converged to 3.9278, 3.9811, 3.9952, 3.9988, so this was a normalization error rather than truncation. Curl and divergence used the planar formulas and dropped the u*tan(lat)/a and -v*tan(lat)/a metric terms. On solid-body rotation that costs exactly a factor of two, which partly cancelled the area inflation and made curl look like a clean 2x on quads but 1.5x on hexagons. Both corrections have to land together: fixing the area alone drives curl to 0.5. Faces whose contour is left open by a missing edge neighbor now return NaN, since no enclosed area exists to divide by. Interior faces of SCRIP-derived grids stay finite, so #1452 remains fixed. Adds manufactured-solution tests with non-zero closed-form answers. The previous suite relied on null tests and ordering comparisons, both of which are satisfied by an operator that is off by a constant factor. Fixes #1662
cmdupuis3
reviewed
Aug 10, 2026
…closed form _compute_gradients_on_faces spent half its time re-deriving the dual-cell area with 4th-order Gaussian quadrature (16 points per fan triangle, ~580 heap allocations per face) when the area of a polygon bounded by great-circle arcs has an exact closed form (spherical excess). The other half was lost to per-edge np.cross/np.linalg.norm temporaries and array-based fill-value checks; both are now scalar, following the zero-allocation convention used elsewhere in the grid kernels. Gradient values shift by up to 5e-6 relative to the old quadrature-based result -- that gap is the retired quadrature's own truncation error, so the new values are the more accurate ones. All existing test assertions hold to 6 significant digits. Added a regression test pinning scale invariance, since the dyamond-30km test grid stores face_x/y/z in meters rather than as unit vectors, which is the one input that distinguishes a correct closed-form implementation from a plausible wrong one. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
The closed-form area reproduced the retired quadrature's unsigned accumulation for bit-compatibility. That convention is wrong for a non-convex cell: a fan from vertex 0 sweeps some triangles backwards and they have to cancel, but taking abs per triangle double-counts them. About 46% of HealPix z4 dual cells are non-convex, differing by up to 20% in area. Measured against d/dlat sin(lat) = cos(lat), the unsigned sum leaves a worst-case gradient error of ~19% that does not improve with resolution (0.186 / 0.192 / 0.196 at z4 / z5 / z6). Signed converges as expected (0.0127 / 0.0088 / 0.0051). Only the tail is affected, which is why the median looked healthy either way; the new test bounds the worst face so a regression to unsigned fails.
erogluorhan
reviewed
Aug 11, 2026
erogluorhan
left a comment
erogluorhan
left a comment
Member
There was a problem hiding this comment.
This looks good to me as long as @cmdupuis3 's concerns are addressed. Thanks @rajeeja !
Collaborator
|
Fine with me, I just wanted a second opinion on the new code before I signed off on it |
cmdupuis3
approved these changes
Aug 11, 2026
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
3 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Fixes #1662.
Green-Gauss on a contour has to divide by the area that same contour encloses.
_compute_gradients_on_faceswalks the dual cell but_normalize_and_project_gradientdivided by the primal face area, so gradients came out inflated byA_dual/A_primal— about 4 on quads, 3 on hexagons. Separately,curlanddivergenceused the planar formulas and dropped theu*tan(lat)/aand-v*tan(lat)/ametric terms, which costs a factor of two on solid-body rotation and partly cancelled the area inflation. Both corrections have to land together, since fixing the area alone drives curl to 0.5.Verified against closed-form answers on the unit sphere, median ratio to exact over faces below 60 degrees:
All were 4.0 / 3.0 / 2.0 / 1.5 before, depending on operator and topology. The quad and hex columns now agree, which they did not previously.
Two behavior changes worth flagging for review. First,
curlanddivergenceof a constant field are no longer zero — on a sphere they reduce to the metric term, and the new values matchtan(lat)/ato machine zero, so I updated those two tests rather than the code. Second, a face whose contour is left open by a missing edge neighbor now returns NaN, because there is no enclosed area to divide by; the all-boundary quad-hexagon fixture is therefore all-NaN now. I checked that #1452 stays fixed — ne30pg2 is 1.2% NaN, not 100% — and added a regression guard for it.There is a cost: gradient on healpix z7 goes from 0.408 s to 1.037 s, about 2.5x, since the dual-cell area is rebuilt per face inside the njit loop. It is not cached the way
Grid.boundsis. Happy to add caching here or as a follow-up if you would rather not take that hit now.Also added manufactured-solution tests with non-zero exact answers. The existing suite was all null tests and ordering comparisons, which a constant-factor error passes — that is how this shipped.
test/coreis green, 137 passed, and pre-commit is clean.