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:

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

Runs performed, all single-rank Release:

./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.