Skip to content

Added support for 3D conservative resampling to regridding.weights(). - #24

Open
roytsmart wants to merge 26 commits into
mainfrom
feature/weights-conservative-3d
Open

Added support for 3D conservative resampling to regridding.weights().#24
roytsmart wants to merge 26 commits into
mainfrom
feature/weights-conservative-3d

Conversation

@roytsmart

Copy link
Copy Markdown
Collaborator

No description provided.

@roytsmart
roytsmart force-pushed the feature/weights-conservative-3d branch 2 times, most recently from 8673d30 to b60b60d Compare May 5, 2025 22:49
@codecov

codecov Bot commented Feb 3, 2026

Copy link
Copy Markdown

Codecov Report

❌ Patch coverage is 91.45861% with 65 lines in your changes missing coverage. Please review.
✅ Project coverage is 91.33%. Comparing base (6c339fd) to head (4c600e9).
⚠️ Report is 4 commits behind head on main.

Files with missing lines Patch % Lines
...eights_conservative_3d/_weights_conservative_3d.py 80.40% 58 Missing ⚠️
regridding/_weights/_weights_conservative.py 20.00% 4 Missing ⚠️
...idding/_weights/_weights_conservative_3d/_grids.py 99.50% 1 Missing ⚠️
...g/_weights/_weights_conservative_3d/_intercepts.py 99.22% 1 Missing ⚠️
...ights/_weights_conservative_3d/_intercepts_test.py 97.05% 1 Missing ⚠️
Additional details and impacted files
@@            Coverage Diff             @@
##             main      #24      +/-   ##
==========================================
+ Coverage   91.27%   91.33%   +0.05%     
==========================================
  Files          33       42       +9     
  Lines        1662     2423     +761     
==========================================
+ Hits         1517     2213     +696     
- Misses        145      210      +65     
Flag Coverage Δ
unittests 91.33% <91.45%> (+0.05%) ⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

@roytsmart
roytsmart force-pushed the feature/weights-conservative-3d branch from 7d85b52 to 7298d2c Compare May 7, 2026 19:43
@roytsmart

Copy link
Copy Markdown
Collaborator Author

Investigation into the sheared-grid failures

Chasing the "little dipoles all over the image" symptom. Findings below, including several dead ends recorded so nobody repeats them.

First: this branch predates perturb

weights() here takes six parameters — no perturb, bounds, or seed. So the degeneracy protection that main applies by default is absent, and a good fraction of what looks broken on this branch is that. Applying main's jitter by hand to the output grid:

wrong entries max abs error negative weights
no perturbation 30 / 729 8.85e-02 5
perturbed output 13 / 729 1.95e-09 0

The rectilinear case becomes essentially exact. This wants a rebase onto main before anything else is judged. In particular, the unguarded det == 0 in line_triangle_intersection_parameters (axis-aligned segments are parallel to four of every cell's six faces) needs no guard — perturbation already covers it.

A sharper test than conservation

Row sums are a weak probe here, and so is regridding a constant field: the failure moves flux between output cells in cancelling pairs, so both come back clean. Two probes that do detect it:

  1. Coalesced weights must be non-negative — they are overlap volumes. Negative entries are an unambiguous bug signal.
  2. Ground truth. For axis-aligned rectilinear grids the exact weights are the product of 1D interval overlaps. For an affine shear (x += s*y) cells stay parallelepipeds with constant Jacobian, so uniform Monte Carlo sampling is exact. Both are cheap and give objective pass/fail.

With perturbation on and an affine shear of 0.4, the sheared bug is real and survives: 26 wrong entries, 3 negatives, max error 3.3e-01.

Verified correct (ruled out)

  • Face triangulation consistency. Adjacent cells agree on the diagonal for every shared face along all three axes, so non-planarity is handled correctly.
  • grid_volume() vs cell_boundary(). These agree to 5e-17 even for bilinear (non-planar) shear. The weights are normalized by a volume consistent with the surface all the geometry uses. This is the property that most often sinks curvilinear 3D schemes, and it holds.
  • cell_axes / cell_normals. All 12 triangles verified against geometric outward normals.
  • _index_input_output. align_axis_right is a cyclic roll and (k,i,j) / (j,k,i) are its exact inverses.

Real defects found, none of which explain the symptom

  • _step_inside_static takes the first face crossing, not the nearest. No t_min tracking, unlike the 2D _step_outside_static, and no last-face exclusion — only a 1e-8 threshold. Sheared cells are slightly non-convex (measured: chords cross 4 triangles rather than 2 in ~0.4% of cases at shear 1.0), so the first hit in index order need not be the exit face. Fixing it changed nothing measurable.
  • The "normals" are edge vectors. normal_static and normal_sweep are differences of adjacent vertices, which equal the face normal only on an orthogonal grid. Measured divergence: 4.0 deg at shear 0.05, 12.0 deg at 0.15, 29.5 deg at 0.4, 54.7 deg at 1.0. Since they feed direction = cross(normal_input, normal_output) and then a sign test that selects which cell a contribution is filed against, this is a discrete failure, not a small perturbation.

What I tried that made things worse

Recorded so they aren't retried. Baseline is 26 wrong / 3 negatives:

attempt wrong negatives
baseline 26 3
t_min nearest crossing 26 3
true face normal for normal_static only 33 5
ditto, as outward normal of the next cell 138 51
true face normals for both 29 5
emit the telescoping pairs together or not at all 97 54

Two lessons. The normals are consumed as a pair whose cross product is the direction of the line where an input face meets an output face, so replacing one and not the other is meaningless. And the three terms in _calc_and_save_weights summing to zero on interior edges is an artifact of uniform cell volumes in my test grids, not an invariant — each cell's share is meant to be emitted independently, as in 2D.

The sharpest remaining lead

The largest errors are exact simple fractions — 1/3 and 1/6 — deposited into cell pairs whose true overlap is exactly zero, and max|error| = 3.3345e-01 was invariant across all six variants above. One is traceable to a single contribution:

weights_00 (lower_left)  index_sweep=(2, 3, 2)  index_static=(2, 0, 2)  +0.333333

index_sweep = (2,3,2) is a boundary vertex. An exact simple fraction is a whole untelescoped tetrahedron, not accumulated round-off, and all six worst entries sit at the last index along the first axis.

High-level: is the approach sound?

The decomposition is valid. ∂(C_i ∩ D_j) consists of patches of ∂C_i inside D_j and patches of ∂D_j inside C_i, and each patch is bounded by exactly two kinds of curve — grid-edge segments inside the other cell, and face-face intersection curves — which are precisely the two things this computes. Same hierarchy (divergence theorem, then Green's theorem per face) that exact polyhedron-intersection methods use.

The difficulty is a qualitative jump from 2D. In 2D each boundary element is a segment with exactly two adjacent cells, so the +/- assignment is purely local: look left, look right. In 3D each boundary element is a patch whose boundary is a cycle of segments, assembled from pieces discovered in two separate phases and arriving in arbitrary order. Integrating over it requires sorting that cycle and orienting it consistently, so the correctness condition becomes global and combinatorial rather than local. A single missing or misordered segment silently corrupts one patch, and no local invariant detects it.

That fits the evidence: every ingredient checks out while the assembled result is wrong, with errors that look like whole tetrahedra filed against the wrong cell. It also explains the rectilinear/sheared split — on an axis-aligned grid the cycles are trivially rectangular and orientation follows from axis alignment, so the fragile machinery is never exercised.

For contrast, direct polyhedron clipping (hexahedra to tets, clip each tet against each candidate cell, sum signed volumes) has a local correctness burden: each (input cell, output cell) pair is independent and self-contained, with no cycle assembly, no ordering, no orientation sign tests. That is why ESMF, YAC, and Omega_h/R3D all take that route for unstructured 3D conservative remapping. #59 is the 2D version of that argument, and it needs no perturbation because it handles degeneracies exactly.

None of which is to say this can't be finished — the decomposition is correct and the ingredients are sound. But the remaining work is in the assembly, which is the part with no local invariant to test against.

Reusable test harnesses

Happy to push these somewhere if useful:

  • exact analytic diff for axis-aligned rectilinear grids (product of 1D interval overlaps)
  • exact Monte Carlo diff for affine-sheared grids, with perturbation applied
  • negative-weight check, non-convexity probe, and a per-edge telescoping diagnostic

🤖 Generated with Claude Code

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant