MeanField/experiments/GEOMETRY_QUALITY_FINDINGS_2026-09-08.md

# Geometry failure diagnosis — 2026-09-08

## Conclusion

The first production Newton correction is small in surface amplitude but produces
large, oscillatory radial displacement gradients near the eight core corners.
The tensor-P3 representation of the prescribed extension distorts the intended
ray-wise profile; the poorly conditioned reference mesh amplifies its gradients.

There is also a verified geometry-safety defect: the production quadrature union
misses inversion between its samples. At the first accepted step length
`alpha = 0.086127771767489855`, the relative mapping determinant is **-0.670218 at
the core corner and -0.428915 at an interior point close to it**, while the
reported limiting quadrature point has determinant **+0.100006**.

Thus this is not merely an approach toward a future singularity: the first
accepted step already corresponds to a folded represented geometry. The
production operators can continue evaluating because their own sampled points
remain admissible.

No physical equations, mapping/extension implementations, Newton policies, or
production tolerances were changed. This investigation implemented an opt-in
experiment and an internal, guarded diagnostic-access hook only.

## Reproduction and artifacts

Build:

```sh
cmake --build cmake-build-release-homebrew --target geometry_quality_experiment -j 6
```

Runs performed, all single-rank Release:

```sh
./cmake-build-release-homebrew/geometry_quality_experiment --output geometry_quality_results_2026-09-08_baseline --tolerances 0.03,0.003 --max-linear-iterations 100 --save-vectors
./cmake-build-release-homebrew/geometry_quality_experiment --output geometry_quality_results_2026-09-08_profile --replay-vectors geometry_quality_results_2026-09-08_baseline/newton_0_vectors.csv --no-fd --no-block-actions
./cmake-build-release-homebrew/geometry_quality_experiment --output geometry_quality_results_2026-09-08_vertices --replay-vectors geometry_quality_results_2026-09-08_baseline/newton_0_vectors.csv --diagonal-only
python3 experiments/summarize_geometry_quality.py geometry_quality_results_2026-09-08_baseline
python3 experiments/summarize_geometry_quality.py geometry_quality_results_2026-09-08_vertices
```

The runs took 510.721, 152.037, and 92.3526 seconds, respectively. Replay checks
the saved accepted state against a fresh context and verifies `Jp+F`, rather than
re-solving for a new correction. Output directories are refused if they already
exist. Help, invalid-tolerance rejection, existing-directory rejection, and normal
context restoration were exercised. `git diff --check` passed. No test suite was
built or run.

The existing build cache contained a nonexistent Eigen 5.0.1 include path. CMake
dependency discovery was refreshed with `-U '*eigen3*'`, locating installed Eigen
3.4.0. See the run provenance and metadata files. The fresh baseline reproduces
the earlier sandbox's seed residual, 23 Krylov iterations, and first safe step.

## 1. Element 0 is the representative of a core-corner pattern

The eight limiting elements are **0, 9, 18, 27, 36, 45, 54, 63**. They are core
elements (attribute 1), not compactified exterior elements. For the Newton
direction their quadrature boundaries differ by only about 0.013%; element 0
has the smallest value. They are near-ties, not exact ties. For uniform
contraction they agree within one part per million.

The actual Newton limiting sample is on element 0's body diagonal:

- Integration coordinates: `(0.010885670927, 0.010885670927, 0.010885670927)`.
- Physical reference position: approximately `(-0.144337126, -0.144337126, -0.144337126)`.
- Physical radius: `0.249999235`, immediately inside the core boundary.
- Reference element Jacobian singular values: minimum `5.51738e-5`, maximum
  `0.108818`, condition number approximately **1,972**.
- At the exact corner the corresponding condition number is approximately **3,151**.

The relative displacement mapping at the seed is the identity. Its determinant
of one therefore conceals the conditioning of the underlying reference element.

Evidence: baseline `newton_0_geometry_elements.csv`,
`uniform_contraction_reference_corner_probes.csv`.

## 2. The dangerous correction is small and nonuniform

Unweighted surface-parameter statistics for the first correction:

| Quantity | Fraction of reference radius |
| --- | ---: |
| Mean | -0.000201759 |
| RMS | 0.000360659 |
| Nonmean RMS | 0.000298945 |
| Most negative | -0.000751271 |
| Most positive | +0.000661834 |

The largest displacement is only **0.0751% of the radius**. The negative extrema
occur at the eight cube-corner surface directions; positive extrema occur near
the twelve edge-center directions. The pattern is nearly invariant under cube
symmetries: the largest spread among symmetry-equivalent nodal directions is
`1.3544e-6`, compared with the full amplitude range `0.0014131`.

The mean-only correction has no sampled boundary through alpha=1. Removing the
mean changes the quadrature boundary from `0.0956975242` to `0.0957208926`, only
0.0244%. The nonmean component therefore accounts for almost the entire local
compression. It is not an oversized uniform radius change or a warm-start-only
artifact; the first solve starts cold.

## 3. Direct diagonal inspection identifies the interpolation mechanism

On the negative core-corner ray, the target surface point is fixed. For the
default radial-power-2 prescription and logical stellar radius one, the intended
continuous radial displacement is exactly

`u_r(s) = a_corner * r_logical(s)^2`, with `a_corner = -0.0007504485104433`.

At the Newton limiting sample:

| Quantity | Intended continuous ray profile | Represented P3 FE field |
| --- | ---: | ---: |
| Radial displacement | -4.6393850e-5 | -3.5188044e-5 |
| Physical radial derivative | -0.4881315 | **-10.4495911** |

The derivative is amplified **21.4 times relative to the intended profile**.
The values agree at the P3 interpolation nodes, but disagree between them. At
`s=0.1`, the FE radial displacement even becomes positive although the intended
profile remains negative.

The code constructs the extension from logical-radius weights and surface-trace
interpolation at volume nodes, then represents it in the tensor-P3 volume space.
The corner element spans three logical max-coordinate sectors, projecting toward
three surface faces. The resulting angular/radial composition is not a single
low-degree tensor polynomial inside that element. Its nodal interpolant need not
preserve the prescribed ray-wise radial behavior.

The physical gradient further multiplies by the inverse reference Jacobian.
Near the core corner, the physical radial coordinate changes extremely slowly
along the logical ray. At the limiting sample, `dr_logical/ds=-0.125` but
`dr_physical/ds=-9.55639e-5`.

This explanation is supported by the measured saved mesh, not just by mesher
source. Sibling STROID source suggests why the core map flattens there, but its
working tree is dirty and was not treated as proof of mesh provenance.

Evidence: profile/vertices `newton_0_core_diagonal.csv`; production implementation
`libmeanfield/impl/deformation/radial_extensions.cpp` and field registry P3
displacement versus saved mesh P4 geometry.

## 4. Quadrature safety is not element safety

The unchanged production estimator was run with a second sampling plan containing
the eight vertices of each of the eight core-corner elements. Independently,
the diagnostic formed the polynomial

`det(J_reference + alpha * dU/dxi) / det(J_reference)`

directly from element matrices, without inverse Jacobians. The methods agree.

| Direction | Production quadrature boundary | Vertex boundary |
| --- | ---: | ---: |
| Original Newton correction, radial power 2 | 0.0956975 | **0.0515682** |
| Same surface correction, radial power 3 | 0.382625 | **0.206273** |
| Same surface correction, radial power 4 | No boundary through 1 | **0.825091** |

For the original correction at alpha=0.0861277718:

| Diagonal coordinate s | Relative determinant |
| --- | ---: |
| 0, exact corner | **-0.670218** |
| 0.005, interior point | **-0.428915** |
| 0.01088567, production limiting quadrature point | +0.100006 |
| 0.02 | +0.609367 |

The negative interior value establishes an actual fold, not merely a singular
boundary vertex. Adding vertices would catch this case, but vertex sampling alone
is not a general positivity certificate: the uniform-contraction control has a
much stricter interior quadrature boundary than vertex boundary.

## 5. Controls and ruled-down explanations

### Linear accuracy

| Requested linear tolerance | Krylov iterations | Verified relative residual | Geometry-safe step |
| --- | ---: | ---: | ---: |
| 0.03 | 23 | 0.0244757 | 0.0861278 |
| 0.003 | 37 | 0.00299622 | 0.0865524 |

Tenfold tightening changes the normalized surface correction by 1.30% and the
safe step by only 0.49%. Linear inaccuracy is not the main explanation for this
first-step failure, although this does not establish behavior at every later
state or arbitrary tolerance.

### Mapping and derivative consistency

- Fresh mapped geometry agrees with the preflight affine prediction to about
  `6e-16` for the actual Newton correction.
- Finite differences of generated volume displacement agree with the production
  extension direction to about `1.6e-16`.
- Material/shape/hydrostatic directional residual differences improve roughly
  tenfold with the first tenfold perturbation reduction, reaching relative errors
  around `1e-7`, then show roundoff/cancellation.
- Gravity checks improve initially but have smaller net Jacobian actions and
  noisier relative errors. Mass finite differences are cancellation-sensitive;
  the smallest perturbation is worse. These checks do not prove the entire
  Jacobian correct, but reveal no gross inconsistency along this correction.
- The central-density relative error of one has absolute magnitude `1.22e-23`;
  it is not evidence of an important phase-border defect.

### Extension and representation controls

Uniform contraction through the default extension has a quadrature boundary of
0.04819 (4.82% surface contraction), independent of Newton/EOS balances.
Core-only P3 projections of physical affine and physical-radius-squared
displacements have boundaries 0.2014 and 0.7693. These are diagnostic bypass
controls, not production solutions. Even physical affine scaling is not exactly
representable by P3 displacement on P4 geometry.

Higher radial powers reduce core sensitivity, but **power 4 is not a demonstrated
fix**: despite all production quadrature points permitting alpha=1, its corner
determinant at alpha=1 is **-0.211997**, and a nearby interior determinant is also
negative. Increasing order or power alone should not be assumed to solve the
problem. In particular, faithfully interpolating the original logical-radial
profile can make the uniform control more singular near a poorly conditioned
reference corner; interpolation sometimes suppresses as well as amplifies it.

## Recommended next work

1. Add this replay direction and the negative interior/corner determinant as a
   targeted geometry regression. Extend admissibility sampling to vertices and
   near-corner/edge regions, then consider adaptive or bounded determinant
   certification for high-order elements. Do not simply reduce the minimum step.
2. Design an extension that controls gradients in physical reference coordinates
   and preserves simple radial/affine modes in the represented FE space. Treat
   reference-mesh quality and displacement/geometry order compatibility together.
3. Rebuild the coupled residual/Jacobian and solve anew under any changed
   extension. Do not accept a post-processed old Newton direction as though it
   solved the new Newton equation.
4. Recheck the mesh-symmetric nonuniform surface pattern and block residuals after
   repairing the geometry representation. The precise origin of that discrete
   pattern, and full nonlinear convergence after a repair, remain unestablished.

The immediate geometric failure mechanism and the sampling defect are now
reproduced. A production remedy and its nonlinear convergence are deliberately
not claimed by this diagnostic experiment.