[5/n][Adjoint Module] Adjoint gradients with respect to block geometries - #3285
Draft
smartalecH wants to merge 10 commits into
Draft
[5/n][Adjoint Module] Adjoint gradients with respect to block geometries#3285smartalecH wants to merge 10 commits into
smartalecH wants to merge 10 commits into
Conversation
geometry_addgradient finite-differences eff_chi1inv_row with respect to an object's centre or size, the same way the design gradient does with respect to a design weight, and contracts against the adjoint field. No extra timestepping. Two things differ from the design gradient. A design weight is local so it is perturbed inside the point loop; a centre is global, so the perturbation hoists out -- twelve passes for six parameters rather than twelve per point. And the geometry cannot simply be mutated, because geom_box_tree caches bounding boxes per node and per object; stale bounds would misreport which object owns a point near the moved boundary, which is exactly where a shape derivative lives. The tree is rebuilt around each perturbation. The step has units here. FD_DEFAULT is a dimensionless perturbation of a weight; a centre is a length, so the step is set as a fraction of a pixel. Subpixel smoothing is required and refused if absent, since without it the permittivity is a step function of position and a finite difference of that is not a derivative. State: correct in 2D TM, 6e-6 to 6e-4 against finite differences for centre and size, and exactly zero for a block translated through a uniform medium. TE is wrong by 1-38% and does not converge with resolution, so the off-diagonal contraction is defective -- see the next commit message. Not yet fit to land.
Three hypotheses eliminated, none of them the cause. The off-diagonal tensor terms are dead code: disabling them is bit-identical in both polarizations. Meep's own design gradient is accurate in TE (4.8e-5) without ever taking that branch either, so second-order smoothing's off-diagonal entries are not what TE needs here. The contraction is kept because it is correct, not because it is load-bearing. The design gradient itself is not at fault: same geometry, same polarization, differentiating the density instead of the boundary gives 4.8e-5. Monitor padding is not the cause. The support of d(epsilon)/d(parameter) is the shell of voxels a boundary sweeps through, which extends outside the object, so padding by two pixels is right in principle -- it grows the monitor from 437 to 621 points -- but the added points contribute exactly zero and TE is unchanged. Kept as correct-in-principle, with no evidence it is needed. TM remains accurate (6e-6 to 6e-4, and exactly zero for a block translated through a uniform medium). TE is 1-38% and diverges with resolution.
The step sweep from the design doc, which I had not run, explains the 2D failures. There is no plateau when a block face sits on a voxel edge: the reported gradient sweeps monotonically from -0.19 to +0.12 as the step goes from 0.2 px to 0.001 px, passing through the true value without settling. Moving the same face half a pixel makes the sweep converge -- 10.1, 9.6, 6.5, 3.4, 1.3, 0.72 percent over that same range. A voxel's filling fraction is piecewise linear in the position of the boundary crossing it, kinked wherever the boundary reaches a voxel edge, so epsilon and the objective are C0 but not C1 in an object's position. On an edge the left and right derivatives differ and no derivative exists; a central difference returns a step-dependent mixture of the two. This is easy to hit by accident, since round sizes at round resolutions land on edges -- a 0.8 block at resolution 20 has faces exactly 8 pixels from its centre. Now warned about. It also means the one case that looked right was a false positive. 2D TM had its faces on edges too, but with a single field component and a symmetric geometry the two neighbouring voxels carried near-identical fields and the mixture cancelled. The default step was also far too large: 0.02 px costs 3.4% on its own. Now 0.002 px. 3D is not explained by either. A smaller step moves centre.x from 38% to 17% but it plateaus near 15%, so there is a second cause there.
Measured rather than assumed. With smoothing off, the objective is piecewise constant in an object's position, swinging 77% across one pixel. With it on the objective is continuous and the swing falls to 14%, so smoothing does deliver the smooth position dependence it promises. The residual oscillates with a period of exactly one pixel, which identifies it as a function of the sub-pixel phase rather than anything physical. Writing J(p) = J_exact(p) + Delta^2 E(p/Delta) and differentiating turns the Delta^2 into Delta, because d/dp of E(p/Delta) carries a 1/Delta. So second-order smoothing gives a first-order-accurate shape derivative. Confirmed: the per-pixel swing is 9.9, 4.0 and 1.9 percent at resolutions 20, 40 and 80 -- halving per doubling, not quartering. This corrects two earlier claims of mine. The kink story does not survive arithmetic: a central difference across a kink returns the average of the one-sided slopes and is step-independent, which is not what was measured. And the 2D TM agreement at 6e-6 was not verification of anything -- it was a symmetric configuration in which the artifact cancelled. The implementation should therefore be tested against a finite difference of the discrete objective, which is what it actually computes. Agreement with the physically intended derivative is separately limited to O(Delta).
…rence
The old approach perturbed the geometry and re-differenced eff_chi1inv_row.
That cannot work: the fill fraction comes from box_overlap_with_object, an
adaptive quadrature, so differencing it amplifies its tolerance by 1/step.
Measured, there was no stable step regime at all -- the reported gradient swept
monotonically across two decades of step without settling.
The smoothed tensor depends on geometry only through the fill fraction and the
interface normal, so
d(chi1inv)/dp = d(chi1inv)/d(fill) * d(fill)/dp
d(fill)/dp is analytic and separable, because a block is an intersection of
three slabs and the pixel overlap factorizes into clamped one-dimensional
overlaps whose derivatives are -1, 0, +1 for a centre and 0, +/-1/2 for a size.
d(chi1inv)/d(fill) comes from varying fill through eff_chi1inv_matrix itself,
via a new fill_override argument, rather than re-deriving Kottke's algebra
here; is linear in fill by construction so that variation is pure local
algebra with no quadrature in the loop.
The result is step-independent to the bit across a 200x range of the internal
fill step, which is the behaviour an analytic derivative should have and the
old code never had.
The loop now visits only pixels with 0 < fill < 1, which is the discrete form
of a shape derivative being a surface integral.
Magnitudes are still wrong, by a factor that varies with resolution (2.56 at
res 20 against 0.54 at res 40) and with index contrast (2.56 at n=2.5 against
0.74 at n=1.5), so something in the contraction is still structurally wrong
rather than merely scaled. Guarding that get_front_object returned the object
being differentiated, and matching its shiftby in the analytic overlap, changes
nothing -- both were already right.
Found with a sum rule that needs no fields and no timestepping: summing d(fill)/d(parameter) times the pixel volume over all pixels must reproduce the change in the object's own volume -- zero for a translation, and the product of the other extents for a size. Two bugs, both of which silently dropped contributions rather than erring. Strict inequalities in axis_overlap. A face lying exactly on a pixel boundary satisfied neither `blo > lo` nor `bhi < hi`, so both neighbouring pixels reported no face inside them and the derivative for that axis vanished entirely. Round geometry on a round grid hits this constantly: a block of width 1.0 at resolution 20 has faces exactly 10 pixels from its centre. Now half-open, which assigns such a face to exactly one pixel -- nothing dropped, nothing double counted, and the derivative one-sided there, which is the truth. The support test was the wrong set. I argued the support of the shape derivative is where d(fill)/d(parameter) is non-zero, then implemented `0 < fill < 1`. Those differ exactly where a face is grid-aligned: no pixel straddles it, fill is 0 or 1 everywhere along it, and meep applies no smoothing there at all -- so the pixels carrying that face's derivative were all excluded. The sum rule showed it plainly: 40 boundary pixels found, which is exactly the two y faces, with none from the two x faces. The sum rule itself belongs in the test suite. It is exact, costs no simulation, and would have caught both of these immediately.
MEEP_CHECK_DFILL compares the analytic d(fill)/d(parameter) against a direct measurement, per pixel, with no fields and no chi1inv in the path. The measurement calls box_overlap_with_object on the object itself, so it needs no geometry tree and cannot be tripped by the stale bounding boxes that moving an object leaves behind. Validating the diagnostic before trusting it: its reference scale comes out at exactly 0.5/Delta -- 10.0 at resolution 20 and 20.0 at resolution 40 -- which is what a size parameter must give when one face sweeps a pixel. Two earlier versions of this same check were wrong and I did trust them. One used a quarter-pixel step, large enough that fill saturated and the reference merely measured the clamp. The other read fill back through interface_fill after moving the object without rebuilding the tree, so get_front_object returned garbage. With a sound reference, the analytic derivative is off by up to 50% on some pixels. So the defect is in d(fill)/d(parameter) specifically, not in d(chi1inv)/d(fill), not in the field contraction, and not in the off-diagonal terms -- those are exactly zero for an axis-aligned block anyway, since normal_to_fixed_object returns the nearest face's normal and the Kottke tensor stays diagonal in Cartesian axes. Next: dump the per-pixel comparison and find which pixels disagree. The error is a constant 5.0 in absolute terms at both resolutions, which does not fit a corner-pixel explanation and needs the pattern to identify.
The field-free per-pixel check now agrees with a direct measurement to 1.2e-12 against a scale of 10, at resolutions 20 and 40. So the analytic fill derivative is right, and whatever remains is elsewhere. Getting there needed one correction and one piece of structure. The voxel-edge warning was inverted. Pixel centres are at integer multiples of dx, so pixel edges are at half-integers: a face at a pixel *centre* is straddled and smooths normally, and a face at a pixel *edge* is straddled by nothing. The warning fired on integers -- the good case -- and the half-pixel offset I had introduced to "fix" an earlier sweep moved faces from centres onto edges. I had been measuring the degenerate configuration and calling it well conditioned. The structure is that this is component-dependent. Yee components sit half a pixel apart, so one component's grid puts a given face mid-pixel while another's puts it exactly on an edge, simultaneously. There is no placement that is well conditioned for every component at once. Faces have to avoid both integer and half-integer multiples of dx, so the safe offset is a quarter pixel, not a half. With that, every component straddles every face and the derivative is two-sided everywhere. Where a face does land on some component's pixel edge, the disagreement is not an error in either quantity: the overlap is clamped at zero or full, so the analytic value is the correct one-sided derivative while a central difference reports the average of the two sides. Both are right; they answer different questions. Still wrong: the assembled gradient, by a factor that varies (0.77, 1.15, 2.18 across resolution and index). With d(fill)/d(parameter) exact, that isolates to d(chi1inv)/d(fill) or the field contraction. The same field-free treatment applies to the former and is the next step.
Verified at two resolutions and two index contrasts, with the residual being
the finite difference's own truncation -- a step sweep shows it still
converging toward the adjoint at h = 0.02 pixels.
The last several rounds of apparent failure were a broken test, not broken
code. t_ratio.py perturbed the block's centre to form the finite difference
while reading the adjoint's *size* entry, left over from an earlier edit that
switched it to size for a sum-rule check. It had been comparing
d/d(center.y) against d/d(size.y). A test that compares centre to centre
agrees to 0.19%.
What actually needed fixing, in order of discovery:
- the geometry finite difference, which differenced an adaptive quadrature
and so amplified its tolerance by 1/step, replaced by the analytic chain
rule d(chi1inv)/dp = d(chi1inv)/d(fill) * d(fill)/dp;
- strict inequalities in axis_overlap, which zeroed a whole axis whenever a
face lay on a pixel boundary;
- the support test, which used 0 < fill < 1 rather than d(fill)/dp != 0;
- an inverted voxel-edge warning: pixel centres are the well conditioned
case and pixel edges the degenerate one, and Yee components sit half a
pixel apart, so faces must avoid both integer and half-integer multiples
of dx -- a quarter-pixel offset is safe for every component at once.
Both operator factors are independently verified rather than inferred:
d(fill)/dp agrees with a direct measurement to 1.2e-12, and d(chi1inv)/d(fill)
matches the closed form implied by Kottke's algebra -- -(eps1-eps2)/eps_avg^2
for the tangential component -- to six digits at two different fill values.
An object opts in with `differentiable=['center', 'size']` and its gradient
comes back under its `name`, alongside the design gradient. Costs no extra
simulation.
Subpixel smoothing makes the permittivity depend on geometry only through each
pixel's filling fraction and the interface normal, so
d(chi1inv)/dp = d(chi1inv)/d(fill) * d(fill)/dp
The first factor is obtained by varying `fill` through eff_chi1inv_matrix
itself, via a new fill_override argument, rather than re-deriving Kottke's
algebra -- `delta` is linear in fill by construction, so that variation is pure
local algebra. The second is analytic: a block is an intersection of slabs, so
the pixel overlap factorizes into clamped one-dimensional overlaps.
Differencing the geometry instead would mean differencing
box_overlap_with_object, an adaptive quadrature, which amplifies its tolerance
by 1/step and has no stable step regime at all.
Only pixels the boundary passes through contribute, since d(fill)/dp vanishes
wherever a pixel is wholly inside or outside. That is the discrete form of a
shape derivative being a surface integral.
Two conventions the caller has to know, both documented and one enforced.
Subpixel smoothing is required, and the gradient refuses without it rather than
returning a number, because an unsmoothed permittivity is a step function of
position. And object faces want to sit a quarter pixel clear of pixel edges:
Yee components are half a pixel apart, so a face at a pixel centre for one
component is on a pixel edge for another, where nothing straddles it and the
derivative is one-sided. A warning fires when that is detected.
Verified against finite differences at 0.1-0.9% across resolution and index
contrast, with the residual being the finite difference's own truncation -- a
step sweep shows it still converging toward the adjoint. Both operator factors
are also checked independently: d(fill)/dp against a direct measurement to
1.2e-12, and d(chi1inv)/d(fill) against the closed form Kottke's algebra
implies, to six digits.
smartalecH
changed the title
smartalecH
changed the title
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
1 participant
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.
No description provided.