tboudreaux/MeanField
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Record negative density mechanism, interval audits, and resolution options

Browse Source
This commit is contained in:
Emily Boudreaux
2026-09-19 10:16:03 -04:00
parent d1a594808f
commit 7ed4d232d0
3 changed files with 668 additions and 2 deletions
Download Patch File Download Diff File Expand all files Collapse all files

184
experiments/NEGATIVE_DENSITY_INVESTIGATION_2026-09-19.md Normal file
View File

@@ -0,0 +1,184 @@
# Negative density: location, mechanism, and implications
Branch: `research/negative-density-diagnosis-20260919`.
Parent: `1148c09` (native spherical harmonic surface validation).
Scope: frozen-state observations and independent reconstructions only. No equilibrium equation, field order, surface parameterization, preconditioner, normalization, or Newton strategy was changed. No new equilibrium was solved in this investigation.
## Main finding
At the measured converged n=1.5 endpoint, the negative density is an undershoot of the **quadratic discontinuous density representation required by the weak EOS closure**. An independently assembled elementwise projection of the nonnegative EOS target sampled by the closure reproduces the stored coefficients within **1.40e-14**, while the negative density is of order 1e-4. This identifies a property of the measured discrete root, rather than an unfinished density-equation solve.
The negative region is in the outermost stellar cells. All 96 surface-adjacent cells contain it; all stored density DOF values remain positive. The initial n=1.5 seed is already negative there before Newton, through a related but distinct interpolation mechanism. A nonrotating n=3 seed has a negative pocket inside the same outer cell layer. A thinner negative layer also exists in the converged n=1 control, despite being missed by q18 and q26 volume sampling.
The practical implication is that a better preconditioner or tighter solve cannot remove the negative polynomial at this same root. We need to assess the density representation/closure before treating convergence as physical validation. This does not prove that the previous Newton plateaus, or the separate hydrostatic-gradient error, have the same cause.
## What the code actually solves
Density and enthalpy are **separate unknowns**. The current defaults are discontinuous L2/Q2 density and continuous H1/Q3 enthalpy; the uniform order increment is zero in these runs. The surface is the existing L=4 all-m candidate, solved simultaneously with the interior in the saved equilibrium runs.
On each stellar element, the density equation is
\[
R_i=\sum_q w_q\,\varphi_i(q)\left[\rho_h(q)-\rho_{\rm EOS}(h_h(q))\right],\qquad
\rho_{\rm EOS}(h)=\left[\frac{h}{(n+1)K}\right]^n\quad(h\ge0).
\]
The weights include the physical deformation. Consequently, a zero density residual gives
\[
M\boldsymbol\rho=\boldsymbol b(h),\qquad
M_{ij}=\sum_q w_q\varphi_i\varphi_j,\quad
b_i=\sum_q w_q\varphi_i\rho_{\rm EOS}(h).
\]
This is a quadrature-defined weighted L2 projection. It is not pointwise EOS equality, and its unique polynomial need not be nonnegative. The local mass matrices were independently assembled and solved with Cholesky; no regularization or clipping was used.
In particular, the implemented dependency is not “negative density produces negative enthalpy.” Preparation requires finite density and rejects sampled negative enthalpy. Pressure is evaluated from enthalpy. Hydrostatic balance constrains `h + phi - centrifugal_potential - C`; density affects it indirectly through the coupled problem. At every negative-density q18/q26 volume sample in this study, enthalpy was positive. The surface traces have exactly zero enthalpy even where the density trace is negative.
Source locations:
- [Independent density and enthalpy declarations](../libmeanfield/interface/field/field_registry.cppm), `Density` and `Enthalpy`.
- [Prepared closure](../libmeanfield/impl/operators/prepared_barotropic_closure.cpp), state checks and `weightedResidual` assembly around lines 603–652; `BuildResidual` around line 681.
- [Pressure preparation](../libmeanfield/impl/operators/prepared_pressure_force.cpp), pressure from enthalpy around line 787.
- [Hydrostatic preparation](../libmeanfield/impl/operators/prepared_hydrostatic_equilibrium.cpp), Bernoulli imbalance around line 863.
- [Seed initialization](../libmeanfield/impl/seed/stellar_equilibrium_projection.cpp), independent `ProjectCoefficient` calls around lines 164–178.
These relative links are repository references. The implementation paths and line numbers refer to the recorded diagnostic branch; the numerical findings are tied to immutable producer commits below.
## Cases and controls
Four frozen states were examined on the existing modest mesh:
| State | Origin | Central-density constraint | Purpose |
|---|---|---|---|
| n=1, J=0 equilibrium | Saved converged harmonic campaign | Absent | Lower-index control |
| n=1.5, J=0 seed | Original Lane–Emden initialization | Absent | Determine whether Newton creates the defect |
| n=1.5, J=0 equilibrium | Saved converged harmonic campaign | Absent | Identify the endpoint closure mechanism |
| n=3, J=0 seed | Original Lane–Emden initialization | Present | Assess relevance to the target index |
The n=3 state is **not a converged n=3 model**, and no rotating n=3 model was tested. Its normalized seed residual is about 3.10e-3. The two saved equilibria retain their previous normalized residuals, approximately 1.53e-12 and 1.30e-12.
The independent projection used the production closure query and separate q18/q26 rules. The production operator's actual rules and assembled residual were exported separately. Independent versus production closure residual differences have L2 norms of roughly 1.2e-17 to 2.1e-17 across the four states. Constant reconstruction errors were also checked independently. These checks distinguish a correct projection reconstruction from an accidental reuse of the production inverse.
For the n=1.5 equilibrium, the production-rule projected density differs from the saved field by only 2.25e-15 in the physical L2 norm. Raising projection quadrature to q26 changes its coefficients by at most 2.36e-7; the negative density persists. Quadrature error in the EOS projection is therefore far too small to account for the measured undershoot.
### Location and sampling limits
The q26 volume samples locate n=1.5 equilibrium negativity at physical radii approximately 0.99839–0.99993 in unit-reference-radius coordinates. All negative samples belong to the 96 cells touching the stellar/nonstellar interface; none occur in the other 736 stellar cells at these quadrature points. The surface is slightly deformed, so a global radius is not the same as distance from the local surface.
For the n=3 seed, q26 samples locate a negative pocket at radii approximately 0.93615–0.97070, with positive enthalpy about 0.0302–0.0682 there. Its surface density trace is positive even though the density becomes negative farther into the same cell. Thus checking only the surface or density DOFs would miss this defect.
| State | Minimum density at q26 samples | Minimum density DOF value | Minimum density on sampled surface traces |
|---|---:|---:|---:|
| n=1 equilibrium | +9.14281e-4 | +1.70519e-2 | -1.14780e-4 |
| n=1.5 seed | -7.46708e-4 | +2.90301e-3 | -9.26358e-4 |
| n=1.5 equilibrium | -5.59609e-4 | +2.83809e-3 | -7.21193e-4 |
| n=3 seed | -8.48632e-5 | +3.26982e-6 | +2.49230e-4 |
All 2,400 sampled surface trace points are negative in the n=1 and n=1.5 cases; none are negative on the n=3 surface itself. These counts include coincident cell-edge points and are not surface-area fractions or global minimum certificates.
Ordinary quadrature of the sign indicator is not reliable for the volume of a thin negative layer. For the n=1.5 equilibrium, the negative mass estimate changes from 2.89e-5 at q18 to 2.18e-5 at q26; its estimated negative volume nearly halves. The field is unchanged. We therefore added a separate diagnostic that solves the density polynomial's zero crossings along reference-cell lines and integrates its negative intervals directly.
### Integrated negative regions
Define negative-mass magnitude as `integral_{rho<0} |rho| dV`, divided below by the absolute total stellar mass. Volume fractions use the entire stellar volume. The finest interval results are:
| Frozen state | Negative mass / stellar mass | Negative volume / stellar volume | Physical radius range of measured negative intervals / reference radius | Inward reference-cell depth range |
|---|---:|---:|---:|---:|
| n=1 equilibrium | 6.09585e-8 | 0.0327346% | 0.999643–0.999886 | 0–0.0007655 |
| n=1.5 seed | 3.66771e-5 | 1.83795% | 0.993813–1.000019 | 0–0.0329556 |
| n=1.5 equilibrium | 2.24535e-5 | 1.49252% | 0.994917–1.001222 | 0–0.0286204 |
| n=3 seed | 2.88372e-5 | 12.0835% | 0.934409–0.978469 | 0.114925–0.349767 |
The inward coordinate is zero on each surface face and one on the opposite face of that cell; it is not a fraction of the whole stellar radius. Physical radii are measured from the center and do not identify the local surface location. Small deviations above unit radius arise in the represented surface, including the curved reference mesh.
Each state has measured negative intervals in all 96 surface-adjacent cells and none on the sampled lines through the other 736 stellar cells. The n=1/n=1.5 intervals reach the surface; the n=3 intervals form an interior pocket. All 24,576 exported finest negative intervals per state lie above the independent auditor's coefficient-uncertainty threshold.
The minimum-density witnesses in this integration are approximately -1.14756e-4, -9.26332e-4, -7.21171e-4, and -8.54870e-5 respectively. These are different tangential samples from the earlier surface-trace table, not certified global minima. Enthalpy at the first three witnesses is zero on the surface; at the n=3 pocket minimum it is +0.0455581.
Tangential rules with 8, 12, and 16 points per direction agree on negative integrals to roundoff in these cases; the two finest negative-mass integrals differ by at most 1.0e-14 relative. Full mass and volume agree with the independently sampled q18/q26 observations within 1.77e-13 relative. Each root-partitioned interval uses eight integration points. Eight analytic C++ controls and a separate decimal-arithmetic audit check roots, signs, minima, weights, sums, source identity, and provenance. The auditor also has four synthetic test methods, including deliberate corruptions.
This is strong convergence evidence for the measured integrals, not a rigorous three-dimensional sign certificate or a bound on every integration error. Tangential locations remain sampled; a sufficiently narrow angular pocket could be missed. Near-double roots below reconstruction uncertainty are not certified. The auditor verifies saved mapped integrals and their agreement with separate whole-volume observations; it is not a second FEM mapping implementation.
![Representative frozen density traces](../../../diagnostics/negative_density_2026-09-19/figures/surface_density_traces.png)
The lines connect saved samples along representative surface-cell traces. The n=1.5 endpoint and independent closure projection overlap; the n=3 panel is an initial seed. The displayed samples illustrate location and are not global extrema.
## Why positive data can produce a negative polynomial
The density DOFs are values at interior Gauss–Legendre points. Positive values at those points do not imply that the interpolating polynomial is positive everywhere in the element. This explains why simply constraining the existing density coefficients would not address the observed defect.
The initial seed and the converged closure use different approximation operations:
1. The seed interpolates an independently tabulated, nonnegative Lane–Emden density into its nodal density space.
2. At equilibrium, the weak EOS closure determines the weighted L2 projection of the separately represented enthalpy's EOS density.
Both operations can undershoot. The nonrotating seeds establish that neither rotational distortion nor an accepted Newton update is necessary for the defect. Its occurrence throughout the outer cell layer points to a radial approximation mechanism, rather than an isolated bad angular mode. This does not establish independence from mesh quality or deformation in more difficult runs. In a simple one-dimensional surface cell, take positive enthalpy proportional to inward distance and normalize the density to `t^n`, with `0 <= t <= 1`. The exact quadratic L2 projections are
\[
p_{1.5}(t)=-\frac{2}{105}+\frac{16}{35}t+\frac47t^2,
\qquad
p_3(t)=\frac1{20}-\frac35t+\frac32t^2.
\]
The first is negative for `0 <= t < 0.0396969`. The cubic-density example has an **interior negative pocket**, `0.118350 < t < 0.281650`, and minimum -0.01 at `t=0.2`. Both retain the exact cell mean. The n=3 example is a smooth cubic, so the problem is not restricted to fractional powers or a nonsmooth EOS derivative. Three-point Gauss–Legendre interpolation also undershoots: it agrees with the cubic example exactly, and produces a slightly deeper surface undershoot for n=1.5. All interpolated DOF values are positive.
This independently explains the two observed kinds of location, without Newton, gravity, surface deformation, rotation, or exterior-domain approximation. It is a mechanism control, not a quantitative prediction for the curved three-dimensional mesh. The n=1 leading-order example is linear and exact. In the actual discretization, Q3 enthalpy need not lie in the Q2 density space; its projection uses the curved physical-volume weights. The toy model omits that higher-order reference-cell variation. Geometry weights alone would not spoil exact reproduction of an EOS density already in Q2. The n=1 endpoint is nevertheless reproduced by the same independent projection to roundoff.
For this fixed-slope toy model, refining the surface-cell width by a factor of two reduces undershoot amplitude by `2^-n` and negative mass per unit area by `2^-(n+1)`. It does not remove the negative fraction of that cell. Seven exact/algebraic self-tests verify the projections, roots, positive DOFs, mean preservation, and refinement laws.
## Why virial balance does not exclude the defect
The constant test function in the density space preserves the EOS density's cell integral under the projection quadrature. Positive and negative pointwise errors can therefore cancel in mass even when they remain locally substantial. Pressure is computed from `h`, while gravitational and other density-dependent terms use the independent density representation. A small enthalpy-based virial error is compatible with this discrepancy.
The n=1.5 endpoint still has enthalpy-based energy and force virial errors approximately 6.55e-7 and 6.01e-7, respectively. Its density-based pressure and virial are undefined wherever density is negative. Neither successful algebraic convergence nor the enthalpy virial establishes pointwise EOS consistency. Negative-mass fractions describe the defect; they are not virial-error estimates and must not be compared directly with the 1e-5 virial tolerance.
The newly detected n=1 face layer also means that the earlier statement “no negative density at the sampled volume points” must remain qualified. Positive volume-quadrature observations are not a positivity certificate, and density-based virials computed from those samples do not prove that the complete field admits a physical density EOS everywhere.
## Candidate resolutions, not implemented here
**A small controlled order/refinement study is the next inexpensive discriminator.** Change density order alone while retaining the surface, enthalpy, gravity, and displacement choices. Q3 density can represent Q3 enthalpy's linear n=1 EOS exactly on compatible reference cells, making the n=1 layer a clean control. For n=3, Q3 density exactly represents the leading linear-enthalpy surface model, but a generic Q3 enthalpy cubed is Q9; Q3 density is not a positivity guarantee. For n=1.5 the composition remains nonpolynomial. Local refinement should reduce the defect but does not generally enforce its sign.
There is no declared stability contract requiring density order to equal gravity-potential order: the existing RT/L2 pairing constraint belongs to gravity flux and potential. The available uniform-order build option would raise multiple fields together and confound a density-only study. Any eventual independent choice should remain in the field registry and preserve its compile-time checks.
**Pointwise EOS-derived density is the preferred architectural candidate for the barotropic model.** Evaluate physical `rho(h)` consistently wherever density is required, as pressure is already evaluated from enthalpy. This removes the particular independent-density projection undershoot, conditional on admissible enthalpy, and gives zero density wherever the surface enthalpy is zero. It also removes an independent field/closure pair. It is a formulation change: gravity sources and forces, centrifugal forces, mass, angular momentum/inertia, the central-density condition, diagnostics, and output must all use the same constitutive density and correct derivatives.
The compiler's state blocks and dependency graph, normalization and preconditioner blocks must change together through a typed formulation choice. Replacing density values in a cache while leaving density-column dependencies intact would be incorrect. This can remain an internal numerical choice without adding complexity to the user's EOS/mass/rotation specification. Enthalpy admissibility still needs to hold; composing the EOS does not by itself prevent an enthalpy polynomial from undershooting. Its convergence, cost, and virial quality require measurement; none is established by this investigation.
Some weak quantities could remain unchanged at fixed enthalpy and geometry: with matched quadrature, the existing projection already preserves the constant mass moment and any gravity-source moments whose test functions lie in the density space. The benefit is a consistent constitutive density throughout the model, not a guarantee that every source term changes or that virial balance automatically improves.
**A positivity-constrained conservative projection is an alternative if an independent polynomial density is needed.** A local weighted projection could enforce nonnegativity and preserve the constant EOS moment. It cannot retain every original Q2 closure moment when the unique unconstrained solution is negative. Nor does preserving cell mass automatically preserve inertia, forces or virial balance. The altered residual and derivatives through its active constraints would need consistent coupling; sampling positivity constraints only guarantees the chosen samples unless a stronger polynomial certificate is supplied.
Mean-preserving scaling constructions are well established for time-dependent high-order DG methods, but their time-step and admissibility hypotheses do not supply a convergence theorem for this steady stellar system. [Zhang and Shu's review](https://www.math.purdue.edu/~zhan1966/research/paper/review.pdf) describes that construction. [Xu and Shu's stationary analysis](https://bpb-us-w2.wpmucdn.com/sites.brown.edu/dist/0/348/files/2022/10/On-the-conservation-of-Positivity-preservin-discontinuous-Galerkind-methods-for-stationary-hyperbolic-equations.pdf) further illustrates why preserving cell averages alone can fail to preserve a stationary discrete balance. Those are methodological references, not ready-made stellar solvers.
### Changes that would not resolve the identified mechanism
- A stronger preconditioner or tighter Newton tolerance would reproduce the same negative density at this root.
- Exact algebraic elimination `rho = M(q)^-1 b(h,q)` also reproduces that polynomial. Its correct derivative is `delta rho = M^-1(delta b - delta M rho)`; eliminating a block does not make this a pointwise EOS formulation.
- A basis change alone leaves the represented polynomial unchanged. Nonnegative coefficients in the current nodal basis are already present.
- Clipping density in diagnostics would conceal the error. Clipping or rescaling it after solving would alter the solved equations and density-dependent physical quantities.
## Implication for the n=3 rotating target
This is useful progress: the investigation isolates a concrete density-representation defect that exists before Newton and remains at a converged root, including a directly relevant n=3 control. It also exposes a validation blind spot in ordinary quadrature sampling. It does not yet establish a valid n=3 equilibrium, explain every nonlinear plateau, or solve the hydrostatic-gradient discrepancy.
Any chosen resolution still needs an accurate n=3 nonrotating anchor, the intended central-density normalization for its homology freedom, then a rotating sequence to the agreed moderate-rotation target. Global virial error <=1e-5 remains required, together with EOS-consistent internal fields and resolution checks suitable for comparison with ESTER.
## Reproducibility and validation
Initial spatial/projection producer: `9c22934d52b369f96364e781d6a4843939ec9523`.
Root-interval producer: `d1a594808f960bbaeeb6df8a8ce6a95458765b19`.
Evidence directory: `../../diagnostics/negative_density_2026-09-19` relative to this checkout.
- `campaign/` contains four serial observations, immutable binary, source/build/mesh provenance, exact commands, input/output hashes, and retained CSVs.
- `campaign_audit.json` verifies the saved algebra, independent production/projection residual agreement, physical-observer agreement, sums, coefficient errors, and hash bindings.
- `regions_campaign/` contains the same four states, observed serially with the root-interval diagnostic. All four completed successfully; no Newton steps were taken.
- `regions_campaign_audit_v1.json` and its log record the independent interval audit. Ninety-six consumed files are hash-bound. Observed state vectors, full density coefficients, and block layouts are byte-identical between the two campaigns.
- `audit_self_checks.json` verifies acceptance of a valid case and rejection of missing hash coverage and a nonmaximal coefficient corruption that does not change the reported maximum error.
- `cli_scope_checks.json` records rejection of accidental Newton execution and incomplete snapshot options.
- `reference_model/` and `reference_model.log` record the independent one-dimensional controls and seven passing self-tests.
- `figures/surface_density_traces.png` displays representative stored-field and independent-projection traces; the n=3 panel is an initial seed, not a solution.
The exact case arguments are committed in `negative_density_cases_2026-09-19.json` and `negative_density_region_cases_2026-09-19.json`. Each campaign contains its private executable, exact command records, build settings and source commit. Both use `sandbox.smesh`, SHA-256 `ce11ae99e21e6c3bbfbee402acd2e191c1da0d8261d2227b4203f73f2f337a74`, with uniform field-order increment zero. Run `audit_negative_density.py` on `campaign/` and `audit_negative_density_regions.py` on `regions_campaign/`; the latter binds automatically to its sibling `campaign/`. The analytic control and region auditor each expose `--self-test`.
The numerical observations are single-rank. No production test suite was rerun because production code was unchanged; the compiled diagnostic, independent mathematical controls, actual-operator residual comparisons, frozen-state invariants and artifact audits are the validation specific to this investigation.

5
experiments/audit_negative_density.py
View File

@@ -198,8 +198,9 @@ def plots(campaign, output):
ax.set_title(title)
ax.set_xlabel('Distance inward from surface / reference radius')
ax.set_ylabel('Density (model units)')
ax.set_xlim(0., x[-1] * .32)
shown = [float(r[k]) for r in line if float(r['inward_coordinate']) <= .3
depth_limit = .08 if name == 'n1p5_equilibrium' else .5
ax.set_xlim(0., x[-1] * depth_limit)
shown = [float(r[k]) for r in line if float(r['inward_coordinate']) <= depth_limit
for k in ('rho', 'rho_from_h', 'rho_projected_production')]
span = max(shown) - min(shown)
ax.set_ylim(min(shown) - .08 * span, max(shown) + .08 * span)

481
experiments/audit_negative_density_regions.py Normal file
View File

@@ -0,0 +1,481 @@
#!/usr/bin/env python3
"""Audit root-split density artifacts against the original frozen campaign.
Saved polynomial algebra and aggregation are checked independently. Mapped
quadrature weights are not exported, so this is not an independent FEM integral
or a global positivity certificate. Tangential convergence remains an observation.
Run after the serial numerical campaign, never alongside a timed FEM run.
"""
from __future__ import annotations
import argparse
from collections import defaultdict
import csv
from decimal import Decimal, localcontext
import hashlib
import json
import math
from pathlib import Path
import tempfile
import unittest
EPS = 2.220446049250313e-16
PREFIX = "density_negative_regions_"
ORDERS = (14, 22, 30)
SUM_FIELDS = ("tangent_lines", "negative_intervals", "integration_points", "volume", "mass",
"negative_volume", "negative_mass")
NEGATIVE_BOUNDS = ("negative_radius_minimum", "negative_radius_maximum",
"negative_interval_t_minimum", "negative_interval_t_maximum")
CONSUMED = ("metadata.txt", "layout.csv", "observed_state.txt", "observed_density_coefficients.txt",
"density_invariants.csv", "density_physical.csv")
def require(condition, message):
if not condition:
raise AssertionError(message)
def close(a, b, label, rtol=3e-11, atol=1e-25):
require(math.isfinite(a) and math.isfinite(b) and math.isclose(a, b, rel_tol=rtol, abs_tol=atol),
f"{label}: {a:.17g} != {b:.17g}")
def digest(path):
with path.open("rb") as stream:
return hashlib.file_digest(stream, "sha256").hexdigest()
def rows(path):
with path.open(newline="") as stream:
reader = csv.DictReader(stream)
require(reader.fieldnames is not None and len(set(reader.fieldnames)) == len(reader.fieldnames),
f"Missing or duplicate CSV columns: {path}")
result = list(reader)
require(all(None not in row and None not in row.values() for row in result), f"Malformed CSV: {path}")
return result
def metadata(path):
entries = [line.split("=", 1) for line in path.read_text().splitlines() if "=" in line]
require(len({key for key, _ in entries}) == len(entries), f"Duplicate metadata: {path}")
return dict(entries)
def keyed(table, key):
result = {key(row): row for row in table}
require(len(result) == len(table), "Duplicate artifact rows")
return result
def numeric(row):
return {key: float(value) for key, value in row.items()}
def finite(row, allowed_nan=()):
for key, value in row.items():
require(math.isfinite(value) or (key in allowed_nan and math.isnan(value)),
f"Unexpected nonfinite value in {key}")
class Campaign:
def __init__(self, path):
self.path = path.resolve()
self.ledger = {}
self.provenance = json.loads(self.record(self.path / "provenance.json").read_text())
p = self.provenance
require(digest(self.record(self.path / "plateau_diagnostic")) == p["binary_sha256"], "Producer binary changed")
require(digest(self.record(Path(p["mesh"]))) == p["mesh_sha256"], "Producer mesh changed")
require(json.loads(self.record(self.path / "cases.json").read_text()) == p["cases"], "Case manifest changed")
self.cases = keyed(p["cases"], lambda case: case["id"])
self.runs = {}
def record(self, path):
path = path.resolve()
self.ledger[str(path)] = digest(path)
return path
def run(self, name):
if name not in self.runs:
r = json.loads(self.record(self.path / f"{name}_run.json").read_text())
require(r["id"] == name and r["exit_code"] == 0 and r["inputs_unchanged"] is True,
f"Incomplete or changed case: {name}")
expected = [str(self.path / "plateau_diagnostic"), *self.cases[name]["arguments"],
"--mesh", self.provenance["mesh"], "--output", str(self.path / name)]
require(r["command"] == expected, f"Recorded command differs from manifest: {name}")
for path, expected_hash in r["input_state_hashes"].items():
require(self.ledger.get(str(self.record(Path(path)))) == expected_hash, f"Changed source input: {path}")
self.runs[name] = r
return self.runs[name]
def artifact(self, name, relative):
r = self.run(name)
require(relative in r["output_hashes"], f"Consumed artifact lacks producer hash: {name}/{relative}")
path = self.record(self.path / name / relative)
require(self.ledger[str(path)] == r["output_hashes"][relative], f"Changed output: {name}/{relative}")
return path
def vector_size(path):
with path.open() as stream:
size = int(next(stream).strip())
require(size > 0, "Empty saved vector")
count = 0
for line in stream:
require(math.isfinite(float(line)), "Nonfinite saved vector")
count += 1
require(size == count, "Saved vector size differs from its header")
return size
def gauss_rule(count):
"""Independent Legendre recurrence on [0,1], used only for node binding."""
result = []
for j in range(1, count + 1):
z = math.cos(math.pi * (j - .25) / (count + .5))
for _ in range(30):
previous, current = 1., z
for degree in range(2, count + 1):
previous, current = current, ((2 * degree - 1) * z * current - (degree - 1) * previous) / degree
derivative = count * (z * current - previous) / (z * z - 1)
update = current / derivative
z -= update
if abs(update) < 2e-16:
break
result.append(((1 + z) / 2, 1 / ((1 - z * z) * derivative * derivative)))
return sorted(result)
def polynomial_intervals(c, b, a):
"""Direct high-precision quadratic formula differs from producer q/a,c/q."""
with localcontext() as context:
context.prec = 70
c, b, a = [Decimal(str(x)) for x in (c, b, a)]
points = [Decimal(0), Decimal(1)]
if a == 0:
roots = [] if b == 0 else [-c / b]
else:
discriminant = b * b - 4 * a * c
roots = [] if discriminant < 0 else [(-b - discriminant.sqrt()) / (2 * a),
(-b + discriminant.sqrt()) / (2 * a)]
points = sorted(set(points + [t for t in roots if 0 < t < 1]))
value = lambda t: c + t * (b + t * a)
return [(float(left), float(right)) for left, right in zip(points, points[1:])
if value((left + right) / 2) < 0]
def interval_check(row, rule):
r = numeric(row)
finite(r)
left, right = r["t_left"], r["t_right"]
require(r["tangential_order"] == 30 and 0 <= left < right <= 1, "Invalid finest negative interval")
require(r["axis"] in (0, 1, 2) and r["orientation"] in (-1, 1), "Invalid interval axis")
a, b, c = (r[key] for key in ("coefficient_a", "coefficient_b", "coefficient_c"))
scale = max(abs(a), abs(b), abs(c))
require(scale > 0, "Zero polynomial cannot have a negative interval")
allowance = 4096 * EPS * scale
evaluate = lambda t: math.fsum((c, b * t, a * t * t))
for endpoint in (left, right):
if 0 < endpoint < 1:
require(abs(evaluate(endpoint)) <= allowance, "Interior interval endpoint is not a root")
for fraction in (.125, .25, .5, .75, .875):
require(evaluate(left + fraction * (right - left)) <= allowance, "Negative interval contains a positive value")
candidates = [left, right]
if a and left < -b / (2 * a) < right:
candidates.append(-b / (2 * a))
minimum = min(map(evaluate, candidates))
require(minimum <= allowance, "Claimed negative interval has positive minimum")
close(r["minimum_density"], minimum, "Direct FE and polynomial interval minimum", rtol=0, atol=allowance)
require(r["negative_volume"] > 0 and r["negative_mass"] > 0, "Nonpositive negative-region integral")
require(r["negative_mass"] / r["negative_volume"] <= -minimum + allowance,
"Negative mean density exceeds the polynomial's maximum magnitude")
weights = []
for coordinate in ("u", "v"):
nearest = min(rule, key=lambda pair: abs(pair[0] - r[coordinate]))
close(nearest[0], r[coordinate], "Tangential Gauss node", rtol=0, atol=4e-15)
weights.append(nearest[1])
close(weights[0] * weights[1], r["tangential_weight"], "Tangential Gauss weight", rtol=2e-12)
radii = [math.sqrt(math.fsum(r[f"{side}_{component}"] ** 2 for component in "xyz"))
for side in ("left", "right")]
require(0 <= r["radius_minimum"] <= r["radius_maximum"], "Invalid radius bounds")
tolerance = 2e-13 * max(1., r["radius_maximum"])
require(r["radius_minimum"] <= min(radii) + tolerance and r["radius_maximum"] >= max(radii) - tolerance,
"Sampled radius bounds exclude a saved endpoint")
return r, minimum < -allowance
def measure_check(r):
empty = r["negative_intervals"] == 0
finite(r, NEGATIVE_BOUNDS if empty else ())
require(r["volume"] > 0 and 0 <= r["negative_volume"] <= r["volume"], "Invalid measured volume")
require(r["negative_mass"] >= 0, "Invalid negative mass")
for key in SUM_FIELDS[:3]:
require(r[key] >= 0 and r[key].is_integer(), f"Noninteger count: {key}")
close(r["negative_volume_fraction"], r["negative_volume"] / r["volume"], "Negative volume fraction")
require(r["mass"] != 0, "Zero total mass requires an explicitly unsupported fraction convention")
close(r["negative_mass_fraction"], r["negative_mass"] / abs(r["mass"]), "Negative mass fraction")
require(0 <= r["maximum_reconstruction_error"] <= r["maximum_reconstruction_tolerance"], "Q2 reconstruction failed")
require(all(0 <= r[key] <= 1 for key in ("witness_t", "witness_u", "witness_v")), "Witness outside reference cell")
if empty:
require(r["negative_mass"] == 0 and r["negative_volume"] == 0 and
all(math.isnan(r[key]) for key in NEGATIVE_BOUNDS), "Empty negative region has nonempty observations")
else:
require(r["negative_mass"] > 0 and r["negative_volume"] > 0 and
0 <= r["negative_interval_t_minimum"] < r["negative_interval_t_maximum"] <= 1 and
0 <= r["negative_radius_minimum"] <= r["negative_radius_maximum"], "Invalid nonempty negative bounds")
def audit_selfchecks(table):
expected = {"positive_quadratic": (0., 0.), "negative_constant": (1., 2.),
"linear_crossing": (.25, .03125), "double_root_nonnegative": (0., 0.),
"near_linear_positive": (0., 0.), "near_double_negative": (2. ** -19, 4 * 2. ** -60 / 3),
"gauss_q2_t_power_1p5": (.04548403175047501373888, .0005262295634284655),
"gauss_q2_t_power_3": (math.sqrt(6) / 15, math.sqrt(6) / 2250)}
indexed = keyed(table, lambda r: r["case"])
require(indexed.keys() == expected.keys(), "Analytic selfcheck coverage changed")
for name, (length, mass) in expected.items():
r = indexed[name]
require(r["passed"] == "1", "Producer arithmetic selfcheck failed")
for key, target in (("negative_length", length), ("negative_mass", mass)):
close(float(r[key]), target, f"Independent selfcheck {name}/{key}", rtol=3e-12, atol=1e-27)
close(float(r["expected_" + key]), target, "Producer selfcheck expected value", rtol=3e-12, atol=1e-27)
def audit_case(current, original, name):
new = lambda filename: current.artifact(name, filename)
old = lambda filename: original.artifact(name, filename)
meta, prior_meta = metadata(new("metadata.txt")), metadata(old("metadata.txt"))
for key in ("n", "K", "J", "rho_c", "M", "R_reference", "density_source", "surface_prescription",
"surface_degree", "surface_maximum_order", "central_phase"):
require(key in meta and key in prior_meta and meta[key] == prior_meta[key], f"Changed model metadata: {key}")
require(meta["density_diagnostic"] == "1" and meta["density_region_integral"] == "1", "Wrong diagnostic mode")
require(current.run(name)["input_state_hashes"] == original.run(name)["input_state_hashes"], "Changed frozen inputs")
for filename in ("observed_state.txt", "observed_density_coefficients.txt", "layout.csv"):
require(digest(new(filename)) == digest(old(filename)), f"Source is not byte-identical: {filename}")
state_size = vector_size(new("observed_state.txt"))
density_size = vector_size(new("observed_density_coefficients.txt"))
layout = rows(new("layout.csv"))
for kind in ("value", "residual"):
blocks = sorted((int(r["offset"]), int(r["size"])) for r in layout if r["kind"] == kind)
end = 0
for offset, size in blocks:
require(offset == end and size > 0, "Invalid saved block layout")
end += size
require(end == state_size, "Layout does not cover the observed state")
require(density_size % 27 == 0, "Full density storage is not Q2 tensor sized")
invariants = {r["metric"]: float(r["value"]) for r in rows(new("density_invariants.csv"))}
previous_invariants = {r["metric"]: float(r["value"]) for r in rows(old("density_invariants.csv"))}
finite(invariants)
for key in ("accepted_physical_unchanged", "accepted_normalized_unchanged", "accepted_residual_unchanged"):
require(invariants[key] == previous_invariants[key] == 1, "Accepted-state observation/restoration failed")
close(invariants["prepared_residual_norm"], previous_invariants["prepared_residual_norm"], "Frozen residual", atol=1e-14)
scope_file = new(PREFIX + "scope.txt")
scope = metadata(scope_file)
require(scope["source_fields_bitwise_identical"] == "1" and scope["trial_preparations"] == "0" and
scope["mpi_ranks"] == "1", "Region observer invariants failed")
for key in ("n", "K"):
close(float(scope[key]), float(meta[key]), "Observer model binding")
text = scope_file.read_text()
for phrase in ("no density or EOS clipping", "All stellar cells are included", "tangential integration is sampled",
"not global extrema or polynomial sign certificates", "Near-double-root"):
require(phrase in text, "Missing observer scope limitation")
audit_selfchecks(rows(new(PREFIX + "selfchecks.csv")))
totals = keyed([numeric(r) for r in rows(new(PREFIX + "summary.csv"))], lambda r: int(r["tangential_order"]))
elements = keyed([numeric(r) for r in rows(new(PREFIX + "elements.csv"))],
lambda r: (int(r["tangential_order"]), int(r["element"])))
topology = keyed(rows(old("density_element_topology.csv")), lambda r: int(r["element"]))
require(totals.keys() == set(ORDERS), "Missing quadrature levels")
require(elements.keys() == {(order, element) for order in ORDERS for element in topology}, "Stellar support differs")
original_summary = keyed(rows(old("density_summary.csv")), lambda r: (int(r["quadrature"]), r["field"]))
physical = {r["metric"]: float(r["value"]) for r in rows(new("density_physical.csv")) if r["quadrature"] == "18"}
denominator_differences = {}
for order, total in totals.items():
require(total["interval_order"] == 14, "Changed interval quadrature")
measure_check(total)
subset = [elements[(order, element)] for element in topology]
for r in subset:
measure_check(r)
top = topology[int(r["element"])]
require(r["interval_order"] == 14 and r["tangent_lines"] == (order // 2 + 1) ** 2 and
r["attribute"] == int(top["attribute"]) and r["surface_face_count"] == int(top["surface_face_count"]),
"Element integration layout or topology differs")
require(r["witness_element"] == r["element"] and r["axis"] in (0, 1, 2) and
r["orientation"] in (-1, 1) and r["integration_points"] % 8 == 0 and
8 * r["tangent_lines"] <= r["integration_points"] <= 24 * r["tangent_lines"], "Element count/axis contradiction")
require((r["surface_face"] >= 0) == (r["surface_face_count"] > 0), "Surface axis metadata inconsistent")
if r["surface_face_count"] == 0:
require(r["axis"] == 0 and r["orientation"] == 1, "Non-surface slicing convention changed")
for key in SUM_FIELDS:
close(math.fsum(r[key] for r in subset), total[key], f"Global element sum {order}/{key}")
close(min(r["minimum_density"] for r in subset), total["minimum_density"], "Global minimum")
witness = elements[(order, int(total["witness_element"]))]
for key in ("minimum_density", "witness_h", "witness_t", "witness_u", "witness_v", "witness_x", "witness_y", "witness_z"):
close(witness[key], total[key], "Global witness binding", atol=1e-15)
nonempty = [r for r in subset if r["negative_intervals"]]
for key, function in (("negative_radius_minimum", min), ("negative_radius_maximum", max),
("negative_interval_t_minimum", min), ("negative_interval_t_maximum", max)):
if nonempty:
close(function(r[key] for r in nonempty), total[key], "Negative-region bounds")
for key in ("maximum_reconstruction_error", "maximum_reconstruction_tolerance"):
close(max(r[key] for r in subset), total[key], "Reconstruction maximum")
denominator_differences[str(order)] = {}
for key in ("mass", "volume"):
for old_order in (18, 26):
reference = float(original_summary[(old_order, "state")][key])
close(total[key], reference, f"Independent full {key}, q{order}/original q{old_order}", rtol=2e-9, atol=1e-13)
denominator_differences[str(order)][f"{key}_relative_to_q{old_order}"] = (total[key] - reference) / abs(reference)
close(physical[key], float(original_summary[(18, "state")][key]), "Repeated q18 physical measure")
interval_rows = rows(new(PREFIX + "intervals.csv"))
by_element, by_line = defaultdict(list), defaultdict(list)
unresolved = 0
rule = gauss_rule(16)
for row in interval_rows:
r, resolved = interval_check(row, rule)
unresolved += not resolved
element = int(r["element"])
require((30, element) in elements, "Interval on an unsupported element")
for key in ("attribute", "axis", "orientation", "surface_face"):
require(r[key] == elements[(30, element)][key], "Interval geometry metadata mismatch")
by_element[element].append(r)
by_line[(element, r["u"], r["v"])].append((row, r, resolved))
for line in by_line.values():
first = line[0][0]
for row, _, _ in line:
require(all(row[key] == first[key] for key in ("coefficient_c", "coefficient_b", "coefficient_a")), "One line has multiple polynomials")
expected = polynomial_intervals(first["coefficient_c"], first["coefficient_b"], first["coefficient_a"])
actual = sorted((r["t_left"], r["t_right"]) for _, r, _ in line)
require(all(actual[i][1] <= actual[i + 1][0] for i in range(len(actual) - 1)), "Overlapping negative intervals")
if all(resolved for _, _, resolved in line):
require(len(expected) == len(actual), "Missing or extra resolved negative interval")
a, b, c = [float(first[k]) for k in ("coefficient_a", "coefficient_b", "coefficient_c")]
uncertainty = 4096 * EPS * max(abs(a), abs(b), abs(c))
for got, wanted in zip(actual, expected):
for x, y in zip(got, wanted):
if y in (0., 1.):
require(x == y, "Negative interval misses the physical cell endpoint")
continue
# Near a double root, position uncertainty scales as the
# square root of coefficient uncertainty, not machine eps.
allowance = max(64 * EPS, uncertainty / max(abs(b + 2 * a * y), math.sqrt(abs(a) * uncertainty)))
close(x, y, "Independent Decimal root", rtol=0, atol=allowance)
for element in topology:
r = elements[(30, element)]
pieces = by_element[element]
require(len(pieces) == r["negative_intervals"], "Finest interval export coverage mismatch")
for key in ("negative_mass", "negative_volume"):
close(math.fsum(piece[key] for piece in pieces), r[key], "Finest interval-to-element sum")
if pieces:
allowance = r["maximum_reconstruction_tolerance"] + 128 * EPS * max(abs(r["minimum_density"]), 1e-30)
close(min(piece["minimum_density"] for piece in pieces), r["minimum_density"], "Finest negative minimum", rtol=0, atol=allowance)
for interval_key, element_key, function in (("t_left", "negative_interval_t_minimum", min),
("t_right", "negative_interval_t_maximum", max), ("radius_minimum", "negative_radius_minimum", min),
("radius_maximum", "negative_radius_maximum", max)):
close(function(p[interval_key] for p in pieces), r[element_key], "Finest interval bounds")
comparison = {}
for key in ("negative_mass", "negative_volume", "negative_mass_fraction", "negative_volume_fraction"):
comparison[key] = {"q14": totals[14][key], "q22": totals[22][key], "q30": totals[30][key],
"q30_minus_q22": totals[30][key] - totals[22][key],
"relative_q30_minus_q22": ((totals[30][key] - totals[22][key]) / totals[30][key]
if totals[30][key] else None)}
return {"artifact_passed": True, "n": float(meta["n"]), "state_size": state_size,
"density_storage_size": density_size, "stellar_elements": len(topology),
"finest_exported_negative_intervals": len(interval_rows), "finest_negative_lines": len(by_line),
"intervals_below_coefficient_uncertainty": unresolved,
"finest_affected_elements": len([v for v in by_element.values() if v]),
"finest_negative_elements_not_touching_surface": sum(bool(by_element[e]) and int(topology[e]["surface_face_count"]) == 0 for e in topology),
"denominator_comparisons": denominator_differences, "quadrature_comparison": comparison,
"summary": {str(k): v for k, v in totals.items()}}
class SelfTests(unittest.TestCase):
def test_decimal_roots(self):
self.assertEqual(polynomial_intervals("-.25", "1", "0"), [(0., .25)])
self.assertEqual(polynomial_intervals(".25", "-1", "1"), [])
roots = polynomial_intervals(".05", "-.6", "1.5")
close(roots[0][0], .2 - math.sqrt(6) / 30, "Cubic root")
close(roots[0][1], .2 + math.sqrt(6) / 30, "Cubic root")
def fixture(self):
rule = gauss_rule(16)
u, weight = rule[3]
r = dict(tangential_order=30, element=1, attribute=1, axis=0, orientation=1, surface_face=2,
u=u, v=u, tangential_weight=weight * weight, t_left=0, t_right=.25,
coefficient_c=-.25, coefficient_b=1, coefficient_a=0,
negative_volume=.25 * weight * weight, negative_mass=.03125 * weight * weight,
minimum_density=-.25, witness_h=1, radius_minimum=math.sqrt(2 * u * u),
radius_maximum=math.sqrt(.25 ** 2 + 2 * u * u), left_x=0, left_y=u, left_z=u,
right_x=.25, right_y=u, right_z=u)
return {k: str(v) for k, v in r.items()}, rule
def test_valid_interval(self):
r, rule = self.fixture()
self.assertTrue(interval_check(r, rule)[1])
def test_corrupt_interval(self):
for key, value in (("t_right", ".3"), ("negative_mass", "1"), ("u", ".5"),
("minimum_density", "-.5"), ("witness_h", "nan")):
with self.subTest(key=key):
r, rule = self.fixture()
r[key] = value
with self.assertRaises(AssertionError):
interval_check(r, rule)
def test_hash_coverage(self):
with tempfile.TemporaryDirectory() as directory:
folder = Path(directory)
(folder / "case").mkdir()
(folder / "case/file.txt").write_text("original")
campaign = object.__new__(Campaign)
campaign.path, campaign.ledger = folder, {}
campaign.runs = {"case": {"output_hashes": {"file.txt": digest(folder / "case/file.txt")}}}
campaign.artifact("case", "file.txt")
(folder / "case/file.txt").write_text("changed")
with self.assertRaises(AssertionError):
campaign.artifact("case", "file.txt")
with self.assertRaises(AssertionError):
campaign.artifact("case", "unrecorded.txt")
def clean(item):
if isinstance(item, dict):
return {str(k): clean(v) for k, v in item.items()}
if isinstance(item, list):
return [clean(v) for v in item]
return None if isinstance(item, float) and not math.isfinite(item) else item
def main():
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument("campaign", type=Path, nargs="?")
parser.add_argument("--original", type=Path, help="Original negative-density campaign; defaults to sibling campaign")
parser.add_argument("--output", type=Path)
parser.add_argument("--self-test", action="store_true")
args = parser.parse_args()
if args.self_test:
result = unittest.TextTestRunner(verbosity=2).run(unittest.defaultTestLoader.loadTestsFromTestCase(SelfTests))
if not result.wasSuccessful():
raise SystemExit(1)
if args.campaign is None:
return
if args.campaign is None or args.output is None:
parser.error("campaign and --output are required for artifact auditing")
current = Campaign(args.campaign)
original = Campaign(args.original or args.campaign.resolve().parent / "campaign")
require(current.path != original.path, "A distinct source campaign is required")
require(current.provenance["mesh_sha256"] == original.provenance["mesh_sha256"], "Campaign meshes differ")
require(current.cases.keys() <= original.cases.keys(), "Region case has no original source case")
report = {"schema": "negative_density_regions_audit_v1", "artifact_passed": True,
"auditor_sha256": digest(Path(__file__)), "producer_commit": current.provenance["commit"],
"original_commit": original.provenance["commit"],
"scope": "Independent saved quadratic algebra, aggregation and source binding; not independent mapped FEM integration or global positivity certification.",
"tolerances": {"aggregation_relative": 3e-11, "aggregation_absolute": 1e-25,
"root_coefficient_roundoff_factor": 4096, "whole_measure_relative": 2e-9,
"whole_measure_absolute": 1e-13},
"cases": {name: audit_case(current, original, name) for name in current.cases}}
report["consumed_sha256"] = {**original.ledger, **current.ledger}
args.output.write_text(json.dumps(clean(report), indent=2, allow_nan=False) + "\n")
print("PASS", args.output)
if __name__ == "__main__":
main()