Added support for 3D conservative resampling to regridding.weights(). - #24
Added support for 3D conservative resampling to regridding.weights().#24roytsmart wants to merge 26 commits into
regridding.weights().#24Conversation
8673d30 to
b60b60d
Compare
Codecov Report❌ Patch coverage is 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
Flags with carried forward coverage won't be shown. Click here to find out more. ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
…ce correct results.
7d85b52 to
7298d2c
Compare
Investigation into the sheared-grid failuresChasing the "little dipoles all over the image" symptom. Findings below, including several dead ends recorded so nobody repeats them. First: this branch predates
|
| 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:
- Coalesced weights must be non-negative — they are overlap volumes. Negative entries are an unambiguous bug signal.
- 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()vscell_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_rightis 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_statictakes the first face crossing, not the nearest. Not_mintracking, unlike the 2D_step_outside_static, and no last-face exclusion — only a1e-8threshold. 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_staticandnormal_sweepare 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 feeddirection = 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
No description provided.