From 9c22934d52b369f96364e781d6a4843939ec9523 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Sat, 19 Sep 2026 09:46:52 -0400 Subject: [PATCH] Add frozen negative density localization and independent EOS projection diagnostics --- experiments/audit_negative_density.py | 230 ++++++++++ .../negative_density_cases_2026-09-19.json | 18 + experiments/negative_density_diagnostics.hpp | 362 ++++++++++++++++ experiments/negative_density_projection.hpp | 396 ++++++++++++++++++ .../negative_density_reference_model.py | 347 +++++++++++++++ experiments/plateau_diagnostic.cpp | 120 ++++++ experiments/run_plateau_campaign.py | 23 + 7 files changed, 1496 insertions(+) create mode 100644 experiments/audit_negative_density.py create mode 100644 experiments/negative_density_cases_2026-09-19.json create mode 100644 experiments/negative_density_diagnostics.hpp create mode 100644 experiments/negative_density_projection.hpp create mode 100644 experiments/negative_density_reference_model.py diff --git a/experiments/audit_negative_density.py b/experiments/audit_negative_density.py new file mode 100644 index 0000000..4190a28 --- /dev/null +++ b/experiments/audit_negative_density.py @@ -0,0 +1,230 @@ +#!/usr/bin/env python3 +"""Audit frozen density observations and summarize their physical location. + +This checks saved arithmetic and provenance, not a second FEM implementation. +The projection and production closure are independently assembled in C++. +""" +import argparse +import csv +import hashlib +import json +import math +from pathlib import Path + +import numpy as np + + +def digest(path): + return hashlib.sha256(path.read_bytes()).hexdigest() + + +def rows(path): + with path.open() as stream: + return list(csv.DictReader(stream)) + + +def vector(path): + value = np.loadtxt(path) + if value.ndim != 1 or int(value[0]) != len(value) - 1 or not np.all(np.isfinite(value)): + raise ValueError(f"Invalid saved vector: {path}") + return value[1:] + + +def require(condition, message): + if not condition: + raise AssertionError(message) + + +def close(left, right, label, rtol=2e-10, atol=1e-15): + require(math.isfinite(left) and math.isfinite(right) and math.isclose(left, right, rel_tol=rtol, abs_tol=atol), + f"{label}: {left} != {right}") + + +def clean(value): + if isinstance(value, dict): + return {str(k): clean(v) for k, v in value.items()} + if isinstance(value, list): + return [clean(v) for v in value] + if isinstance(value, float) and not math.isfinite(value): + return None + return value + + +def numeric(row, omit=()): + return {key: float(value) for key, value in row.items() if key not in omit} + + +def audit_case(folder): + meta = dict(line.split('=', 1) for line in (folder / 'metadata.txt').read_text().splitlines() if '=' in line) + n, K = float(meta['n']), float(meta['K']) + summary = {(int(r['quadrature']), r['field']): numeric(r, ('quadrature', 'field')) + for r in rows(folder / 'density_summary.csv')} + element_rows = rows(folder / 'density_elements.csv') + for (order, field), totals in summary.items(): + matching = [numeric(r, ('quadrature', 'field', 'element', 'attribute', 'surface_layer')) + for r in element_rows if int(r['quadrature']) == order and r['field'] == field] + for key in ('samples', 'negative_samples', 'negative_h_samples', 'joint_negative_samples', + 'volume', 'mass', 'negative_mass', 'negative_volume', 'negative_h_volume', 'closure_squared_integral'): + close(math.fsum(r[key] for r in matching), totals[key], f"element sum {order}/{field}/{key}") + close(min(r['minimum_density'] for r in matching), totals['minimum_density'], 'element minimum') + close(max(r['maximum_density'] for r in matching), totals['maximum_density'], 'element maximum') + radial = [numeric(r, ('quadrature', 'field', 'radial_bin_left', 'radial_bin_right')) + for r in rows(folder / 'density_radial_bins.csv') + if int(r['quadrature']) == order and r['field'] == field] + for key in ('volume', 'mass', 'negative_mass', 'negative_volume'): + close(math.fsum(r[key] for r in radial), totals[key], f"radial sum {order}/{field}/{key}") + + physical = {r['metric']: float(r['value']) for r in rows(folder / 'density_physical.csv')} + observed18, observed26 = summary[(18, 'state')], summary[(26, 'state')] + for key, physical_key in [('minimum_density', 'minimum_density'), ('minimum_h', 'minimum_enthalpy'), + ('volume', 'volume'), ('mass', 'mass'), ('negative_mass', 'negative_density_mass_magnitude'), + ('negative_samples', 'negative_density_samples'), ('negative_h_samples', 'negative_enthalpy_samples')]: + close(observed18[key], physical[physical_key], f"physical observer {key}") + close(observed18['negative_volume'] / observed18['volume'], physical['negative_density_volume_fraction'], + 'physical negative volume fraction') + negative = rows(folder / 'negative_density_samples.csv') + require(len(negative) == observed26['negative_samples'], 'Sparse negative sample count does not bind to complete observation.') + for r in negative: + rho, h = float(r['rho']), float(r['h']) + require(rho < 0 and int(r['quadrature']) == 26, 'Unexpected negative sample record.') + close(float(r['radius']), math.sqrt(sum(float(r[d]) ** 2 for d in ('x', 'y', 'z'))), 'negative sample radius') + close(float(r['rho_from_h']), (h / ((n + 1) * K)) ** n if h > 0 else 0., 'negative sample EOS') + close(math.fsum(float(r['weight']) for r in negative), observed26['negative_volume'], 'negative export volume') + close(math.fsum(-float(r['weight']) * float(r['rho']) for r in negative), observed26['negative_mass'], 'negative export mass') + + projection_rows = rows(folder / 'negative_density_projection_elements.csv') + production = {int(r['element']): r for r in projection_rows if r['control'] == 'production'} + actual_rules = {int(r['element']): r for r in rows(folder / 'production_closure_rules.csv')} + require(production.keys() == actual_rules.keys(), 'Projection support differs from the actual closure.') + for element, r in production.items(): + rule = actual_rules[element] + # MFEM may report the achieved odd Gauss degree for an even requested + # degree. Equal points and this exact one-degree alias are intentional. + requested, achieved = int(r['quadrature_order']), int(rule['quadrature_order']) + require(achieved in (requested, requested + 1) and r['quadrature_points'] == rule['quadrature_points'], + f"Production quadrature mismatch in element {element}.") + actual = vector(folder / 'production_closure_residual_full.txt') + reconstructed = vector(folder / 'projection_production_closure_residual_full.txt') + require(actual.shape == reconstructed.shape, 'Full closure residual layout mismatch.') + difference = actual - reconstructed + rhs_norm = math.sqrt(math.fsum(float(r['rhs_norm']) ** 2 for r in production.values())) + allowance = 2048 * np.finfo(float).eps * rhs_norm + defect_norm = float(np.linalg.norm(difference)) + require(defect_norm <= allowance, f"Independent closure residual mismatch {defect_norm}, allowance {allowance}.") + + saved = vector(folder / 'observed_density_coefficients.txt') + coefficient_comparison = {} + for control in ('production', 'q18', 'q26'): + projected = vector(folder / f'negative_density_projection_{control}_coefficients.txt') + require(projected.shape == saved.shape, 'Projected coefficient layout mismatch.') + delta = projected - saved + matching = [r for r in projection_rows if r['control'] == control] + maximum = float(np.max(np.abs(delta))) + close(maximum, max(float(r['coefficient_maximum_difference']) for r in matching), 'coefficient error summary') + for row in matching: + require(float(row['constant_coefficient_error_maximum']) < 1e-10 and + float(row['constant_sample_error_maximum']) < 1e-10, + 'Independent constant reproduction failed.') + coefficient_comparison[control] = dict(maximum_absolute_difference=maximum, + l2_coefficient_difference=float(np.linalg.norm(delta)), + relative_l2_coefficient_difference=float(np.linalg.norm(delta) / np.linalg.norm(saved))) + invariants = {r['metric']: float(r['value']) for r in rows(folder / 'density_invariants.csv')} + require(all(invariants[k] == 1 for k in ('accepted_physical_unchanged', 'accepted_normalized_unchanged', + 'accepted_residual_unchanged')), 'Accepted state changed during observation.') + topology = {int(r['element']): r for r in rows(folder / 'density_element_topology.csv')} + affected = {int(r['element']) for r in negative} + layers, attributes = {}, {} + for r in negative: + element, weight = int(r['element']), float(r['weight']) + require(r['surface_layer'] == topology[element]['surface_layer'], 'Negative sample topology mismatch.') + for table, key in ((layers, int(r['surface_layer'])), (attributes, int(r['attribute']))): + item = table.setdefault(key, {'volume': 0., 'negative_mass': 0.}) + item['volume'] += weight + item['negative_mass'] -= weight * float(r['rho']) + traces = rows(folder / 'density_surface_traces.csv') + surface = [r for r in traces if float(r['inward_coordinate']) == 0.] + for r in traces: + require(all(math.isfinite(float(r[key])) for key in r), 'Nonfinite surface trace observation.') + h = float(r['h']) + close(float(r['rho_from_h']), (h / ((n + 1) * K)) ** n if h > 0 else 0., 'surface trace EOS') + result = dict(n=n, K=K, source=meta['density_source'], normalized_residual=invariants['prepared_residual_norm'], + state={str(q): summary[(q, 'state')] for q in (18, 26)}, + controls={field: {str(q): summary[(q, field)] for q in (18, 26)} for field in ('eos_pointwise', 'production', 'q18', 'q26')}, + coefficient_comparison=coefficient_comparison, + closure_residual_check=dict(actual_norm=float(np.linalg.norm(actual)), independent_difference_norm=defect_norm, + arithmetic_allowance=allowance, rhs_norm=rhs_norm), + affected_elements=len(affected), stellar_elements=len(topology), + surface_elements=sum(int(r['surface_layer']) == 0 for r in topology.values()), + negative_volume_by_surface_layer=layers, negative_volume_by_attribute=attributes, + density_dof_minimum=min(float(r['density_dof_minimum']) for r in topology.values()), + surface_trace=dict(samples=len(surface), negative_samples=sum(float(r['rho']) < 0 for r in surface), + minimum_density=min(float(r['rho']) for r in surface), + maximum_abs_enthalpy=max(abs(float(r['h'])) for r in surface)), + full_trace_minimum_density=min(float(r['rho']) for r in traces), + physical_q18=physical, artifact_passed=True) + return clean(result) + + +def plots(campaign, output): + import matplotlib + matplotlib.use('Agg') + import matplotlib.pyplot as plt + output.mkdir(parents=True, exist_ok=True) + fig, axes = plt.subplots(1, 2, figsize=(11, 4.4), constrained_layout=True) + for ax, name, title in zip(axes, ('n1p5_equilibrium', 'n3_seed'), ('n = 1.5 converged state', 'n = 3 initial state')): + table = rows(campaign / name / 'density_surface_traces.csv') + centers = [r for r in table if float(r['tangent_u']) == .5 and float(r['tangent_v']) == .5] + witness = min(centers, key=lambda r: float(r['rho'])) + line = [r for r in centers if r['face'] == witness['face']] + line.sort(key=lambda r: float(r['distance_from_face'])) + x = [float(r['distance_from_face']) for r in line] + for key, label, style in [('rho_from_h', 'Positive EOS density from h', 'k-'), + ('rho', 'Stored quadratic density', 'o-'), + ('rho_projected_production', 'Independent closure projection', '--')]: + ax.plot(x, [float(r[key]) for r in line], style, label=label, markersize=3) + ax.axhline(0., color='0.55', linewidth=.8) + 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 + 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) + ax.grid(alpha=.18) + axes[0].legend(fontsize=8) + fig.savefig(output / 'surface_density_traces.png', dpi=180) + plt.close(fig) + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument('campaign', type=Path) + parser.add_argument('--output', type=Path, required=True) + parser.add_argument('--plots', type=Path) + args = parser.parse_args() + campaign = args.campaign.resolve() + provenance = json.loads((campaign / 'provenance.json').read_text()) + require(digest(campaign / 'plateau_diagnostic') == provenance['binary_sha256'], 'Producer binary changed.') + require(digest(Path(provenance['mesh'])) == provenance['mesh_sha256'], 'Input mesh changed.') + results = {} + for case in provenance['cases']: + name = case['id'] + run = json.loads((campaign / (name + '_run.json')).read_text()) + require(run['exit_code'] == 0 and run['inputs_unchanged'], f'Incomplete frozen diagnostic: {name}') + for relative, value in run['output_hashes'].items(): + require(digest(campaign / name / relative) == value, f'Changed output: {name}/{relative}') + for path, value in run['input_state_hashes'].items(): + require(digest(Path(path)) == value, f'Changed input: {path}') + results[name] = audit_case(campaign / name) + payload = dict(source_commit=provenance['commit'], artifact_passed=True, auditor_sha256=digest(Path(__file__)), + scope='Saved algebra, projection/production agreement and provenance; sampled spatial localization is not a sign certificate.', + cases=results) + args.output.write_text(json.dumps(clean(payload), indent=2, allow_nan=False) + '\n') + if args.plots: + plots(campaign, args.plots) + print('PASS', args.output) + + +if __name__ == '__main__': + main() diff --git a/experiments/negative_density_cases_2026-09-19.json b/experiments/negative_density_cases_2026-09-19.json new file mode 100644 index 0000000..fa3cb0b --- /dev/null +++ b/experiments/negative_density_cases_2026-09-19.json @@ -0,0 +1,18 @@ +[ + { + "id": "n1_equilibrium", + "arguments": ["--index", "1", "--angular-momentum", "0", "--surface-prescription", "harmonics", "--surface-degree", "4", "--no-central", "--max-newton", "0", "--density-diagnostic", "--density-snapshot", "../../diagnostics/spherical_harmonic_surface_2026-09-19/campaign/n1_J0_L4"] + }, + { + "id": "n1p5_seed", + "arguments": ["--index", "1.5", "--angular-momentum", "0", "--surface-prescription", "harmonics", "--surface-degree", "4", "--no-central", "--max-newton", "0", "--density-diagnostic"] + }, + { + "id": "n1p5_equilibrium", + "arguments": ["--index", "1.5", "--angular-momentum", "0", "--surface-prescription", "harmonics", "--surface-degree", "4", "--no-central", "--max-newton", "0", "--density-diagnostic", "--density-snapshot", "../../diagnostics/spherical_harmonic_surface_2026-09-19/campaign/n1p5_J0_L4"] + }, + { + "id": "n3_seed", + "arguments": ["--index", "3", "--angular-momentum", "0", "--surface-prescription", "harmonics", "--surface-degree", "4", "--max-newton", "0", "--density-diagnostic"] + } +] diff --git a/experiments/negative_density_diagnostics.hpp b/experiments/negative_density_diagnostics.hpp new file mode 100644 index 0000000..e6d8e61 --- /dev/null +++ b/experiments/negative_density_diagnostics.hpp @@ -0,0 +1,362 @@ +#pragma once + +// Experiment-only observations of a frozen physical state. No sampled field, +// projection, or positivity condition is installed in the equilibrium problem. +#include "negative_density_projection.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace experiment::plateau_diagnostics { + namespace negative_density_detail { + inline std::ofstream Open(const std::filesystem::path &path) { + if (std::filesystem::exists(path)) + throw std::runtime_error("Diagnostic output already exists: " + path.string()); + std::ofstream out(path); + out.exceptions(std::ios::badbit | std::ios::failbit); + out << std::setprecision(17); + return out; + } + + struct SurfaceFace { + int face, element, axis; + double coordinate; + }; + + struct Topology { + std::vector surface; + std::vector layer; + std::vector> faces; + }; + + inline Topology BuildTopology(polytrope_validation::PhysicalState &state) { + auto &mesh = *state.finiteElements.mesh; + Topology result; + result.layer.assign(mesh.GetNE(), -1); + result.faces.resize(mesh.GetNE()); + std::vector> neighbors(mesh.GetNE()); + std::queue queue; + for (int face = 0; face < mesh.GetNumFaces(); ++face) { + int left = -1, right = -1; + mesh.GetFaceElements(face, &left, &right); + const bool leftStar = left >= 0 && state.isStellar(left); + const bool rightStar = right >= 0 && state.isStellar(right); + if (leftStar && rightStar) { + neighbors[left].push_back(right); + neighbors[right].push_back(left); + } + if (leftStar == rightStar) + continue; + const int element = leftStar ? left : right; + auto *transformation = mesh.GetFaceElementTransformations(face); + if (transformation == nullptr) + throw std::runtime_error("Missing stellar interface transformation."); + mfem::IntegrationPoint facePoint, point; + facePoint.Set2(0.5, 0.5); + if (leftStar) + transformation->Loc1.Transform(facePoint, point); + else + transformation->Loc2.Transform(facePoint, point); + const std::array coordinate{point.x, point.y, point.z}; + int axis = -1; + double value = 0.; + for (int d = 0; d < 3; ++d) { + if (std::abs(coordinate[d]) < 1e-12 || std::abs(coordinate[d] - 1.) < 1e-12) { + if (axis >= 0) + throw std::runtime_error("Ambiguous reference surface face."); + axis = d; + value = std::round(coordinate[d]); + } + } + if (axis < 0) + throw std::runtime_error("Density localization currently requires tensor-product cells."); + result.faces[element].push_back(static_cast(result.surface.size())); + result.surface.push_back({face, element, axis, value}); + if (result.layer[element] < 0) { + result.layer[element] = 0; + queue.push(element); + } + } + while (!queue.empty()) { + const int element = queue.front(); + queue.pop(); + for (const int neighbor : neighbors[element]) + if (result.layer[neighbor] < 0) { + result.layer[neighbor] = result.layer[element] + 1; + queue.push(neighbor); + } + } + for (int element = 0; element < mesh.GetNE(); ++element) + if (state.isStellar(element) && result.layer[element] < 0) + throw std::runtime_error("Stellar component is disconnected from its surface."); + return result; + } + + struct Aggregate { + std::size_t samples = 0, negative = 0, negativeH = 0, jointNegative = 0; + long double volume = 0., mass = 0., negativeMass = 0., negativeVolume = 0., negativeHVolume = 0., + closure2 = 0.; + double minimum = INFINITY, maximum = -INFINITY, minimumH = INFINITY; + double negativeHMinimum = INFINITY, negativeHMaximum = -INFINITY; + double negativeRadiusMinimum = INFINITY, negativeRadiusMaximum = -INFINITY; + double negativeReferenceDepthMinimum = INFINITY, negativeReferenceDepthMaximum = -INFINITY; + void + Add(double rho, + double h, + double rhoH, + double radius, + double depth, + double weight) { + ++samples; + volume += weight; + mass += weight * rho; + closure2 += weight * (rho - rhoH) * (rho - rhoH); + minimum = std::min(minimum, rho); + maximum = std::max(maximum, rho); + minimumH = std::min(minimumH, h); + if (h < 0.) { + ++negativeH; + negativeHVolume += weight; + } + if (rho >= 0.) + return; + ++negative; + if (h < 0.) + ++jointNegative; + negativeMass -= weight * rho; + negativeVolume += weight; + negativeHMinimum = std::min(negativeHMinimum, h); + negativeHMaximum = std::max(negativeHMaximum, h); + negativeRadiusMinimum = std::min(negativeRadiusMinimum, radius); + negativeRadiusMaximum = std::max(negativeRadiusMaximum, radius); + if (std::isfinite(depth)) { + negativeReferenceDepthMinimum = std::min(negativeReferenceDepthMinimum, depth); + negativeReferenceDepthMaximum = std::max(negativeReferenceDepthMaximum, depth); + } + } + }; + + inline double Observed(double v) { + return std::isfinite(v) ? v : std::numeric_limits::quiet_NaN(); + } + inline void Header(std::ostream &out) { + out << "samples,negative_samples,negative_h_samples,joint_negative_samples,volume,mass,negative_mass," + "negative_volume,negative_h_volume,closure_squared_integral,minimum_density,maximum_density,minimum_" + "h,minimum_h_at_negative_rho,maximum_h_at_negative_rho,negative_radius_minimum,negative_radius_" + "maximum,negative_reference_depth_minimum,negative_reference_depth_maximum\n"; + } + inline void Write( + std::ostream &out, + const Aggregate &a + ) { + out << a.samples << ',' << a.negative << ',' << a.negativeH << ',' << a.jointNegative << ',' << a.volume + << ',' << a.mass << ',' << a.negativeMass << ',' << a.negativeVolume << ',' << a.negativeHVolume << ',' + << a.closure2 << ',' << Observed(a.minimum) << ',' << Observed(a.maximum) << ',' << Observed(a.minimumH) + << ',' << Observed(a.negativeHMinimum) << ',' << Observed(a.negativeHMaximum) << ',' + << Observed(a.negativeRadiusMinimum) << ',' << Observed(a.negativeRadiusMaximum) << ',' + << Observed(a.negativeReferenceDepthMinimum) << ',' << Observed(a.negativeReferenceDepthMaximum) + << '\n'; + } + + inline double ReferenceDepth( + const Topology &topology, + int element, + const mfem::IntegrationPoint &point + ) { + const std::array coordinates{point.x, point.y, point.z}; + double depth = INFINITY; + for (const int index : topology.faces[element]) { + const auto &face = topology.surface[index]; + depth = std::min(depth, std::abs(coordinates[face.axis] - face.coordinate)); + } + return Observed(depth); + } + } // namespace negative_density_detail + + inline void RecordNegativeDensity( + polytrope_validation::PhysicalState &state, + double n, + double K, + const std::filesystem::path &output + ) { + namespace d = negative_density_detail; + const auto topology = d::BuildTopology(state); + const auto projections = MeasureNegativeDensityProjection(state, n, K, output); + auto &mesh = *state.finiteElements.mesh; + auto topologyFile = d::Open(output / "density_element_topology.csv"); + topologyFile << "element,attribute,surface_layer,surface_face_count,density_order,enthalpy_order,density_dof_" + "minimum,density_dof_maximum\n"; + mfem::Array dofs; + mfem::Vector coefficients; + for (int element = 0; element < mesh.GetNE(); ++element) + if (state.isStellar(element)) { + auto *transform = state.finiteElements.densityFes->GetElementVDofs(element, dofs); + state.density.GetSubVector(dofs, coefficients); + if (transform) + transform->InvTransformPrimal(coefficients); + topologyFile << element << ',' << mesh.GetAttribute(element) << ',' << topology.layer[element] << ',' + << topology.faces[element].size() << ',' + << state.finiteElements.densityFes->GetFE(element)->GetOrder() << ',' + << state.finiteElements.enthalpyFes->GetFE(element)->GetOrder() << ',' + << coefficients.Min() << ',' << coefficients.Max() << '\n'; + } + auto globalFile = d::Open(output / "density_summary.csv"); + globalFile << "quadrature,field,"; + d::Header(globalFile); + auto elementFile = d::Open(output / "density_elements.csv"); + elementFile << "quadrature,field,element,attribute,surface_layer,"; + d::Header(elementFile); + auto radialFile = d::Open(output / "density_radial_bins.csv"); + radialFile << "quadrature,field,radial_bin_left,radial_bin_right,"; + d::Header(radialFile); + auto negativeFile = d::Open(output / "negative_density_samples.csv"); + negativeFile << "quadrature,element,attribute,surface_layer,xi,eta,zeta,x,y,z,radius,reference_inward_depth," + "weight,rho,h,rho_from_h"; + for (const auto &control : projections.controls) + negativeFile << ",rho_projected_" << control.label; + negativeFile << '\n'; + std::vector labels{"state", "eos_pointwise"}; + for (const auto &control : projections.controls) + labels.push_back(control.label); + const double factor = (n + 1.) * K; + mean_field::mapping::MappingPointContext mapped; + for (const int order : {18, 26}) { + std::cout << "Density localization quadrature " << order << std::endl; + std::vector global(labels.size()); + std::vector> radial(labels.size()); + for (int element = 0; element < mesh.GetNE(); ++element) { + if (!state.isStellar(element)) + continue; + auto *transformation = mesh.GetElementTransformation(element); + if (transformation->GetGeometryType() != mfem::Geometry::CUBE) + throw std::runtime_error("Density localization requires hexahedral elements."); + const auto &rule = mfem::IntRules.Get(transformation->GetGeometryType(), order); + std::vector local(labels.size()); + for (int q = 0; q < rule.GetNPoints(); ++q) { + const auto &point = rule.IntPoint(q); + if (state.mapping.EvaluatePoint(*transformation, point, mapped) != + mean_field::mapping::MappingStatus::valid) + throw std::runtime_error("Invalid mapping during density localization."); + transformation->SetIntPoint(&point); + const double weight = point.weight * transformation->Jacobian().Det() * mapped.mapping_determinant; + if (!std::isfinite(weight) || weight <= 0.) + throw std::runtime_error("Invalid density observation weight."); + const double rho = state.density.GetValue(element, point); + const double h = state.enthalpy.GetValue(element, point); + if (!std::isfinite(rho) || !std::isfinite(h)) + throw std::runtime_error("Nonfinite observed density or enthalpy."); + const double rhoH = h > 0. ? std::pow(h / factor, n) : 0.; + const double radius = mapped.physical_position.Norml2(); + const double depth = d::ReferenceDepth(topology, element, point); + std::vector densities{rho, rhoH}; + for (const auto &control : projections.controls) + densities.push_back(control.density->GetValue(element, point)); + const int bin = static_cast(std::floor(radius * 200.)); + for (std::size_t field = 0; field < labels.size(); ++field) { + if (!std::isfinite(densities[field])) + throw std::runtime_error("Nonfinite projected density."); + local[field].Add(densities[field], h, rhoH, radius, depth, weight); + global[field].Add(densities[field], h, rhoH, radius, depth, weight); + radial[field][bin].Add(densities[field], h, rhoH, radius, depth, weight); + } + // Retain every negative-state point at the finer rule. Volume + // statistics above include all points, so this sparse export + // cannot silently redefine the denominator of a fraction. + if (rho < 0. && order == 26) { + negativeFile << order << ',' << element << ',' << mesh.GetAttribute(element) << ',' + << topology.layer[element] << ',' << point.x << ',' << point.y << ',' << point.z; + for (int axis = 0; axis < 3; ++axis) + negativeFile << ',' << mapped.physical_position(axis); + negativeFile << ',' << radius << ',' << depth << ',' << weight << ',' << rho << ',' << h << ',' + << rhoH; + for (std::size_t field = 2; field < densities.size(); ++field) + negativeFile << ',' << densities[field]; + negativeFile << '\n'; + } + } + for (std::size_t field = 0; field < labels.size(); ++field) { + elementFile << order << ',' << labels[field] << ',' << element << ',' << mesh.GetAttribute(element) + << ',' << topology.layer[element] << ','; + d::Write(elementFile, local[field]); + } + } + for (std::size_t field = 0; field < labels.size(); ++field) { + globalFile << order << ',' << labels[field] << ','; + d::Write(globalFile, global[field]); + for (const auto &[bin, value] : radial[field]) { + radialFile << order << ',' << labels[field] << ',' << bin / 200. << ',' << (bin + 1) / 200. << ','; + d::Write(radialFile, value); + } + } + globalFile.flush(); + elementFile.flush(); + radialFile.flush(); + negativeFile.flush(); + } + + // Surface-inward traces include the face itself, which volume Gauss rules + // never sample. The direction is a reference-cell coordinate line, not a + // physical normal or an assumed spherical radial line. + auto traces = d::Open(output / "density_surface_traces.csv"); + traces << "face,element,tangent_u,tangent_v,inward_coordinate,x,y,z,radius,distance_from_face,rho,h,rho_from_h"; + for (const auto &control : projections.controls) + traces << ",rho_projected_" << control.label; + traces << '\n'; + constexpr std::array tangents{0., .25, .5, .75, 1.}; + constexpr std::array depths{0., .001, .003, .007, .01, .02, .03, .04, .06, + .08, .1, .15, .2, .3, .5, .75, 1.}; + for (const auto &face : topology.surface) { + auto *transformation = mesh.GetElementTransformation(face.element); + const int axis1 = (face.axis + 1) % 3, axis2 = (face.axis + 2) % 3; + for (const double u : tangents) + for (const double v : tangents) { + mfem::Vector boundaryPosition(3); + for (const double depth : depths) { + std::array coordinate{}; + coordinate[face.axis] = face.coordinate == 0. ? depth : 1. - depth; + coordinate[axis1] = u; + coordinate[axis2] = v; + mfem::IntegrationPoint point; + point.Set3(coordinate[0], coordinate[1], coordinate[2]); + if (state.mapping.EvaluatePoint(*transformation, point, mapped) != + mean_field::mapping::MappingStatus::valid) + throw std::runtime_error("Invalid surface-inward trace mapping."); + if (depth == 0.) + boundaryPosition = mapped.physical_position; + mfem::Vector distance(mapped.physical_position); + distance -= boundaryPosition; + const double rho = state.density.GetValue(face.element, point), + h = state.enthalpy.GetValue(face.element, point); + traces << face.face << ',' << face.element << ',' << u << ',' << v << ',' << depth; + for (int axis = 0; axis < 3; ++axis) + traces << ',' << mapped.physical_position(axis); + traces << ',' << mapped.physical_position.Norml2() << ',' << distance.Norml2() << ',' << rho + << ',' << h << ',' << (h > 0. ? std::pow(h / factor, n) : 0.); + for (const auto &control : projections.controls) + traces << ',' << control.density->GetValue(face.element, point); + traces << '\n'; + } + } + } + auto scope = d::Open(output / "density_scope.txt"); + scope << "Frozen observations only; no equilibrium state, operator or field space is modified.\n" + << "density_basis=" << state.finiteElements.densityFes->FEColl()->Name() << '\n' + << "enthalpy_basis=" << state.finiteElements.enthalpyFes->FEColl()->Name() << '\n' + << "surface_faces=" << topology.surface.size() << '\n' + << "Surface layer is graph distance through stellar face neighbors to a stellar/nonstellar interface.\n" + << "Reference inward depth is available only in cells touching that interface.\n" + << "Radial bins are in current physical coordinates with the diagnostic's zero center and unit reference " + "radius.\n" + << "Quadrature extrema and trace extrema are observations, not polynomial sign certificates.\n"; + } +} // namespace experiment::plateau_diagnostics diff --git a/experiments/negative_density_projection.hpp b/experiments/negative_density_projection.hpp new file mode 100644 index 0000000..b6bc646 --- /dev/null +++ b/experiments/negative_density_projection.hpp @@ -0,0 +1,396 @@ +#pragma once + +// Include after import mean_field. These element-local projections diagnose the +// weak closure at a frozen geometry and enthalpy; they are never state updates. +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "polytrope_physical_state.hpp" + +namespace experiment::plateau_diagnostics { + +struct DensityProjectionControl final { + std::string label; + // -1 denotes the production query, whose resolved order is recorded for + // each element. Explicit controls use ordinary MFEM rules of this order. + int quadratureOrder; + bool productionQuadrature; + std::unique_ptr density; +}; + +struct NegativeDensityProjectionReport final { + // The spaces remain borrowed from PhysicalState::finiteElements. Values on + // vacuum elements are zero storage, not a physical extension of the star. + std::vector controls; +}; + +namespace negative_density_projection_detail { + +inline double Norm(const mfem::Vector& vector) { + long double squared = 0.0L; + for (int i = 0; i < vector.Size(); ++i) { + if (!std::isfinite(vector(i))) throw std::domain_error("Nonfinite projection vector."); + squared += static_cast(vector(i)) * vector(i); + } + const double norm = static_cast(std::sqrt(squared)); + if (!std::isfinite(norm)) throw std::domain_error("Projection norm overflow."); + return norm; +} + +inline double Relative(const double numerator, const double denominator) { + return denominator > 0.0 ? numerator / denominator + : (numerator == 0.0 ? 0.0 : std::numeric_limits::infinity()); +} + +inline void WriteVector(const std::filesystem::path& path, const mfem::Vector& vector) { + std::ofstream stream(path); + stream.exceptions(std::ios::failbit | std::ios::badbit); + stream << std::setprecision(17) << vector.Size() << '\n'; + for (int i = 0; i < vector.Size(); ++i) stream << vector(i) << '\n'; +} + +// The independent, small Cholesky solve keeps this diagnostic separate from +// the production density inverse. A bad mass matrix must not be regularized: +// that would change the equation whose negative solution we are investigating. +inline mfem::DenseMatrix Factor(const mfem::DenseMatrix& mass) { + const int size = mass.Height(); + mfem::DenseMatrix lower(size); + lower = 0.0; + double scale = 0.0; + for (int i = 0; i < size; ++i) scale = std::max(scale, mass(i, i)); + if (!(scale > 0.0) || !std::isfinite(scale)) + throw std::domain_error("Invalid density projection mass diagonal."); + for (int i = 0; i < size; ++i) { + for (int j = 0; j <= i; ++j) { + long double value = mass(i, j); + for (int k = 0; k < j; ++k) value -= static_cast(lower(i, k)) * lower(j, k); + if (i == j) { + if (!std::isfinite(value) || value <= 64.0L * std::numeric_limits::epsilon() * scale) + throw std::domain_error("Density projection mass is singular or numerically unresolved."); + lower(i, j) = static_cast(std::sqrt(value)); + } else { + lower(i, j) = static_cast(value / lower(j, j)); + if (!std::isfinite(lower(i, j))) throw std::domain_error("Nonfinite mass factor."); + } + } + } + return lower; +} + +inline mfem::Vector Solve(const mfem::DenseMatrix& lower, const mfem::Vector& rhs) { + mfem::Vector result(rhs); + for (int i = 0; i < result.Size(); ++i) { + long double value = result(i); + for (int j = 0; j < i; ++j) value -= static_cast(lower(i, j)) * result(j); + result(i) = static_cast(value / lower(i, i)); + } + for (int i = result.Size() - 1; i >= 0; --i) { + long double value = result(i); + for (int j = i + 1; j < result.Size(); ++j) value -= static_cast(lower(j, i)) * result(j); + result(i) = static_cast(value / lower(i, i)); + } + (void)Norm(result); + return result; +} + +inline mfem::Vector Defect(const mfem::DenseMatrix& mass, const mfem::Vector& coefficients, + const mfem::Vector& rhs) { + mfem::Vector result(rhs.Size()); + for (int i = 0; i < rhs.Size(); ++i) { + long double value = -static_cast(rhs(i)); + for (int j = 0; j < coefficients.Size(); ++j) + value += static_cast(mass(i, j)) * coefficients(j); + result(i) = static_cast(value); + } + return result; +} + +inline void VerifySolve(const mfem::DenseMatrix& mass, const mfem::Vector& coefficients, + const mfem::Vector& rhs, const mfem::Vector& defect) { + long double massSquared = 0.0L; + for (int i = 0; i < mass.Height(); ++i) + for (int j = 0; j < mass.Width(); ++j) + massSquared += static_cast(mass(i, j)) * mass(i, j); + const long double scale = Norm(rhs) + std::sqrt(massSquared) * Norm(coefficients); + if (!std::isfinite(scale) || Norm(defect) > 1.0e-11L * scale) + throw std::runtime_error("Local density projection failed its independent normal-equation check."); +} + +inline bool Identical(const mfem::Vector& first, const mfem::Vector& second) { + return first.Size() == second.Size() && + (first.Size() == 0 || std::memcmp(first.GetData(), second.GetData(), + sizeof(double) * static_cast(first.Size())) == 0); +} + +struct Totals final { + int elements{0}; + long double volume{0.0L}, projectedSavedSquared{0.0L}, projectedTargetSquared{0.0L}, savedTargetSquared{0.0L}; + long double projectedDefectSquared{0.0L}, savedDefectSquared{0.0L}; + double minimumSavedCoefficient{std::numeric_limits::infinity()}; + double minimumProjectedCoefficient{std::numeric_limits::infinity()}; + double maximumCoefficientDifference{0.0}, maximumConstantCoefficientError{0.0}, maximumConstantSampleError{0.0}; +}; + +} // namespace negative_density_projection_detail + +inline NegativeDensityProjectionReport MeasureNegativeDensityProjection( + experiment::polytrope_validation::PhysicalState& state, const double n, const double K, + const std::filesystem::path& output) { + namespace detail = negative_density_projection_detail; + using namespace mean_field; + const eos::Polytrope equationOfState(n, K); + const auto& fem = state.finiteElements; + if (!fem.mesh || !fem.densityFes || !fem.enthalpyFes || !fem.compactificationCoordinate || + !fem.quadratureFactory || fem.mesh->Dimension() != 3 || !fem.mesh->Conforming() || + fem.densityFes->GetVDim() != 1 || + dynamic_cast(fem.densityFes->FEColl()) == nullptr || + state.density.ParFESpace() != fem.densityFes.get() || state.enthalpy.ParFESpace() != fem.enthalpyFes.get()) + throw std::invalid_argument("Density projection diagnostics require the current scalar conforming L2 density space."); + int ranks = 0; + MPI_Comm_size(fem.mesh->GetComm(), &ranks); + if (ranks != 1) throw std::invalid_argument("Density projection diagnostics require one MPI rank."); + + const double extra = (n - 1.0) * static_cast(field::Enthalpy::Scalar::familyOrder); + if (!std::isfinite(extra) || extra < 0.0 || extra > std::numeric_limits::max()) + throw std::invalid_argument("Invalid production EOS quadrature order."); + const int extraOrder = static_cast(std::ceil(extra)); + std::filesystem::create_directories(output); + const auto elementsPath = output / "negative_density_projection_elements.csv"; + const auto summaryPath = output / "negative_density_projection_summary.csv"; + const auto productionResidualPath = output / "projection_production_closure_residual_full.txt"; + if (std::filesystem::exists(elementsPath) || std::filesystem::exists(summaryPath) || + std::filesystem::exists(productionResidualPath)) + throw std::invalid_argument("Density projection diagnostics require fresh output files."); + + NegativeDensityProjectionReport report; + for (const auto& [label, order] : std::array, 3>{ + std::pair{"production", -1}, std::pair{"q18", 18}, std::pair{"q26", 26}}) { + const auto vectorPath = output / (std::string("negative_density_projection_") + label + "_coefficients.txt"); + if (std::filesystem::exists(vectorPath)) + throw std::invalid_argument("Density projection diagnostics would overwrite saved coefficients."); + auto projected = std::make_unique(fem.densityFes.get()); + *projected = 0.0; + report.controls.push_back({label, order, order < 0, std::move(projected)}); + } + const mfem::Vector savedDensity(state.density), savedEnthalpy(state.enthalpy), savedPotential(state.potential), + savedGravity(state.gravityGradientReference), savedDisplacement(state.displacement); + // Even the observer's mapping cache is private here. No trial preparation, + // normalization, accepted buffer, or prepared operator is involved. + mapping::GridFunctionMappingEvaluator mapper(state.domainMapper, state.displacement, *fem.compactificationCoordinate); + mfem::ParGridFunction constant(fem.densityFes.get()); + mfem::ConstantCoefficient one(1.0); + constant.ProjectCoefficient(one); + mfem::Vector fullClosureResidual(fem.densityFes->GetVSize()); + fullClosureResidual = 0.0; + std::array totals; + std::vector visited(static_cast(fem.densityFes->GetVSize()), false); + + std::ofstream elements(elementsPath), summary(summaryPath); + elements.exceptions(std::ios::failbit | std::ios::badbit); + summary.exceptions(std::ios::failbit | std::ios::badbit); + elements << std::setprecision(17); + summary << std::setprecision(17) << "control,metric,value\n"; + elements << "control,requested_order,quadrature_order,quadrature_points,element,attribute,density_order,enthalpy_order," + "basis_family,dofs,volume,minimum_h,minimum_target_density,saved_coefficient_minimum,saved_coefficient_maximum," + "projected_coefficient_minimum,projected_coefficient_maximum,coefficient_difference_norm,coefficient_relative_difference," + "coefficient_maximum_difference,rhs_norm,projected_normal_defect_norm,projected_normal_relative_defect," + "saved_normal_defect_norm,saved_normal_relative_defect,projection_vs_saved_l2,projection_vs_target_l2,saved_vs_target_l2," + "minimum_projected_sampled_density,minimum_saved_sampled_density,constant_normal_defect_norm," + "constant_coefficient_error_maximum,constant_sample_error_maximum\n"; + + for (int element = 0; element < fem.mesh->GetNE(); ++element) { + if (!state.isStellar(element)) continue; + const auto* densityElement = fem.densityFes->GetFE(element); + const auto* enthalpyElement = fem.enthalpyFes->GetFE(element); + if (!densityElement || !enthalpyElement || densityElement->GetMapType() != mfem::FiniteElement::VALUE || + densityElement->GetOrder() != field::Density::Scalar::familyOrder || + enthalpyElement->GetOrder() != field::Enthalpy::Scalar::familyOrder) + throw std::invalid_argument("Density projection FE elements do not match the production scalar fields."); + mfem::Array dofs; + if (fem.densityFes->GetElementDofs(element, dofs) != nullptr) + throw std::invalid_argument("Density projection does not support transformed element DOFs."); + if (dofs.Size() != densityElement->GetDof()) throw std::logic_error("Density element layout mismatch."); + for (int i = 0; i < dofs.Size(); ++i) { + if (dofs[i] < 0 || dofs[i] >= fullClosureResidual.Size() || visited[static_cast(dofs[i])]) + throw std::invalid_argument("Density projection requires independent, unsigned element-local DOFs."); + visited[static_cast(dofs[i])] = true; + } + mfem::Vector saved, constantExpected; + state.density.GetSubVector(dofs, saved); + constant.GetSubVector(dofs, constantExpected); + (void)detail::Norm(saved); + auto* transformation = fem.mesh->GetElementTransformation(element); + + for (std::size_t controlIndex = 0; controlIndex < report.controls.size(); ++controlIndex) { + auto& control = report.controls[controlIndex]; + auto& total = totals[controlIndex]; + const mfem::IntegrationRule* rule = nullptr; + int resolvedOrder = control.quadratureOrder; + if (control.productionQuadrature) { + // PreparedBarotropicClosureOperator uses this source query for + // both rho*q and rho(h)*q, not the separate legacy mass query. + const auto query = field::Field::make_query( + quadrature::QuadratureRole::discretization, transformation->OrderW(), std::array{extraOrder}, + utils::DOMAINS::STELLAR, quadrature::MappingKind::general); + const auto resolution = fem.quadratureFactory->get(query, transformation->GetGeometryType()); + rule = resolution.integration_rule; + resolvedOrder = resolution.resolution.order; + } else { + rule = &mfem::IntRules.Get(transformation->GetGeometryType(), resolvedOrder); + } + if (!rule || rule->GetNPoints() == 0) throw std::runtime_error("Missing density projection quadrature."); + const int size = dofs.Size(), points = rule->GetNPoints(); + mfem::DenseMatrix mass(size), samples(points, size); + mfem::Vector rhs(size), constantRhs(size), shape(size), weights(points), targets(points); + mass = 0.0; + rhs = 0.0; + constantRhs = 0.0; + double minimumH = std::numeric_limits::infinity(); + double minimumTarget = minimumH; + long double volume = 0.0L; + for (int q = 0; q < points; ++q) { + const auto& ip = rule->IntPoint(q); + mapping::MappingPointContext mapped; + if (mapper.EvaluatePoint(*transformation, ip, mapped) != mapping::MappingStatus::valid || mapped.compactified) + throw std::domain_error("Invalid stellar mapping during density projection."); + const double h = state.enthalpy.GetValue(element, ip); + if (!std::isfinite(h) || h < 0.0) + throw std::domain_error("Density projection encountered negative or nonfinite enthalpy; EOS was not clipped."); + const double target = eos::evaluate( + equationOfState, dimensions::SpecificEnthalpyValue{h}).value(); + transformation->SetIntPoint(&ip); + const double referenceDeterminant = transformation->Jacobian().Det(); + const double weight = ip.weight * referenceDeterminant * mapped.mapping_determinant; + if (!(referenceDeterminant > 0.0) || !std::isfinite(referenceDeterminant) || + !(mapped.mapping_determinant > 0.0) || !std::isfinite(mapped.mapping_determinant) || + !(weight > 0.0) || !std::isfinite(weight) || target < 0.0 || !std::isfinite(target)) + throw std::domain_error("Invalid density projection weight or EOS value."); + densityElement->CalcShape(ip, shape); + (void)detail::Norm(shape); + weights(q) = weight; + targets(q) = target; + volume += weight; + minimumH = std::min(minimumH, h); + minimumTarget = std::min(minimumTarget, target); + for (int i = 0; i < size; ++i) { + samples(q, i) = shape(i); + rhs(i) += weight * shape(i) * target; + constantRhs(i) += weight * shape(i); + for (int j = 0; j < size; ++j) mass(i, j) += weight * shape(i) * shape(j); + } + } + const auto lower = detail::Factor(mass); + const auto coefficients = detail::Solve(lower, rhs); + const auto constantCoefficients = detail::Solve(lower, constantRhs); + const auto projectedDefect = detail::Defect(mass, coefficients, rhs); + const auto savedDefect = detail::Defect(mass, saved, rhs); + const auto constantDefect = detail::Defect(mass, constantCoefficients, constantRhs); + detail::VerifySolve(mass, coefficients, rhs, projectedDefect); + detail::VerifySolve(mass, constantCoefficients, constantRhs, constantDefect); + control.density->SetSubVector(dofs, coefficients); + if (control.productionQuadrature) fullClosureResidual.SetSubVector(dofs, savedDefect); + + mfem::Vector difference(coefficients); + difference -= saved; + double minimumSaved = std::numeric_limits::infinity(), maximumSaved = -minimumSaved; + double minimumProjected = minimumSaved, maximumProjected = -minimumSaved; + double maximumDifference = 0.0, constantCoefficientError = 0.0, constantSampleError = 0.0; + for (int i = 0; i < size; ++i) { + minimumSaved = std::min(minimumSaved, saved(i)); + maximumSaved = std::max(maximumSaved, saved(i)); + minimumProjected = std::min(minimumProjected, coefficients(i)); + maximumProjected = std::max(maximumProjected, coefficients(i)); + maximumDifference = std::max(maximumDifference, std::abs(difference(i))); + constantCoefficientError = std::max(constantCoefficientError, std::abs(constantCoefficients(i) - constantExpected(i))); + } + long double projectedSavedSquared = 0.0L, projectedTargetSquared = 0.0L, savedTargetSquared = 0.0L; + double minimumProjectedSample = std::numeric_limits::infinity(), minimumSavedSample = minimumProjectedSample; + for (int q = 0; q < points; ++q) { + long double projectedValue = 0.0L, savedValue = 0.0L, constantValue = 0.0L; + for (int i = 0; i < size; ++i) { + projectedValue += static_cast(samples(q, i)) * coefficients(i); + savedValue += static_cast(samples(q, i)) * saved(i); + constantValue += static_cast(samples(q, i)) * constantCoefficients(i); + } + const long double differenceValue = projectedValue - savedValue; + const long double projectedError = projectedValue - targets(q), savedError = savedValue - targets(q); + projectedSavedSquared += weights(q) * differenceValue * differenceValue; + projectedTargetSquared += weights(q) * projectedError * projectedError; + savedTargetSquared += weights(q) * savedError * savedError; + minimumProjectedSample = std::min(minimumProjectedSample, static_cast(projectedValue)); + minimumSavedSample = std::min(minimumSavedSample, static_cast(savedValue)); + constantSampleError = std::max(constantSampleError, static_cast(std::abs(constantValue - 1.0L))); + } + if (!std::isfinite(constantCoefficientError) || !std::isfinite(constantSampleError) || + constantCoefficientError > 1.0e-10 || constantSampleError > 1.0e-10) + throw std::runtime_error("Density projection failed independent constant reconstruction."); + const double rhsNorm = detail::Norm(rhs), projectedDefectNorm = detail::Norm(projectedDefect); + const double savedDefectNorm = detail::Norm(savedDefect); + elements << control.label << ',' << control.quadratureOrder << ',' << resolvedOrder << ',' << points << ',' + << element << ',' << fem.mesh->GetAttribute(element) << ',' << densityElement->GetOrder() << ',' + << enthalpyElement->GetOrder() << ',' << std::quoted(fem.densityFes->FEColl()->Name()) << ',' << size << ',' + << volume << ',' << minimumH << ',' << minimumTarget << ',' << minimumSaved << ',' << maximumSaved << ',' + << minimumProjected << ',' << maximumProjected << ',' << detail::Norm(difference) << ',' + << detail::Relative(detail::Norm(difference), detail::Norm(saved)) << ',' << maximumDifference << ',' + << rhsNorm << ',' << projectedDefectNorm << ',' << detail::Relative(projectedDefectNorm, rhsNorm) << ',' + << savedDefectNorm << ',' << detail::Relative(savedDefectNorm, rhsNorm) << ',' + << std::sqrt(projectedSavedSquared) << ',' << std::sqrt(projectedTargetSquared) << ',' + << std::sqrt(savedTargetSquared) << ',' << minimumProjectedSample << ',' << minimumSavedSample << ',' + << detail::Norm(constantDefect) << ',' << constantCoefficientError << ',' << constantSampleError << '\n'; + ++total.elements; + total.volume += volume; + total.projectedSavedSquared += projectedSavedSquared; + total.projectedTargetSquared += projectedTargetSquared; + total.savedTargetSquared += savedTargetSquared; + total.projectedDefectSquared += static_cast(projectedDefectNorm) * projectedDefectNorm; + total.savedDefectSquared += static_cast(savedDefectNorm) * savedDefectNorm; + total.minimumSavedCoefficient = std::min(total.minimumSavedCoefficient, minimumSaved); + total.minimumProjectedCoefficient = std::min(total.minimumProjectedCoefficient, minimumProjected); + total.maximumCoefficientDifference = std::max(total.maximumCoefficientDifference, maximumDifference); + total.maximumConstantCoefficientError = std::max(total.maximumConstantCoefficientError, constantCoefficientError); + total.maximumConstantSampleError = std::max(total.maximumConstantSampleError, constantSampleError); + } + } + const bool fieldsIdentical = detail::Identical(savedDensity, state.density) && detail::Identical(savedEnthalpy, state.enthalpy) && + detail::Identical(savedPotential, state.potential) && detail::Identical(savedGravity, state.gravityGradientReference) && + detail::Identical(savedDisplacement, state.displacement); + if (!fieldsIdentical) throw std::logic_error("Density projection diagnostic changed a source field."); + summary << "all,source_fields_bitwise_identical,1\nall,trial_preparations,0\nall,mpi_ranks,1\n" + "all,unsupported_elements_are_zero_storage,1\nall,negative_enthalpy_is_rejected,1\n" + "all,sampled_extrema_are_not_polynomial_certificates,1\n"; + for (std::size_t index = 0; index < report.controls.size(); ++index) { + const auto& control = report.controls[index]; + const auto& total = totals[index]; + if (total.elements == 0) throw std::invalid_argument("Density projection found no stellar elements."); + const auto metric = [&](const char* name, const auto value) { summary << control.label << ',' << name << ',' << value << '\n'; }; + metric("elements", total.elements); + metric("local_mass_solves", 2 * total.elements); + metric("volume", total.volume); + metric("projection_vs_saved_l2", std::sqrt(total.projectedSavedSquared)); + metric("projection_vs_target_l2", std::sqrt(total.projectedTargetSquared)); + metric("saved_vs_target_l2", std::sqrt(total.savedTargetSquared)); + metric("projected_normal_defect_norm", std::sqrt(total.projectedDefectSquared)); + metric("saved_normal_defect_norm", std::sqrt(total.savedDefectSquared)); + metric("saved_coefficient_minimum", total.minimumSavedCoefficient); + metric("projected_coefficient_minimum", total.minimumProjectedCoefficient); + metric("coefficient_maximum_difference", total.maximumCoefficientDifference); + metric("constant_coefficient_error_maximum", total.maximumConstantCoefficientError); + metric("constant_sample_error_maximum", total.maximumConstantSampleError); + detail::WriteVector(output / ("negative_density_projection_" + control.label + "_coefficients.txt"), *control.density); + } + detail::WriteVector(productionResidualPath, fullClosureResidual); + return report; +} + +} // namespace experiment::plateau_diagnostics diff --git a/experiments/negative_density_reference_model.py b/experiments/negative_density_reference_model.py new file mode 100644 index 0000000..316a2ad --- /dev/null +++ b/experiments/negative_density_reference_model.py @@ -0,0 +1,347 @@ +#!/usr/bin/env python3 +"""Independent surface-cell controls; no FEM state, solver, or limiter is changed. + +For inward distance s in [0, width], assume h=kappa*s and rho=A*s**n. +Compare exact, uniformly weighted Q2 L2 projection with interpolation at MFEM's +default L2 Q2 Gauss-Legendre nodes. In one dimension Q2 means degree <= 2. +Extrema, zeros, and negative integrals use quadratic algebra, not plot sampling. + +The normalized model is t**n, t=s/width. Physical density is +A*width**n*p(t). Curved metrics, tangential variation, higher terms in h, +and changes in kappa under refinement are deliberately outside this control. + +Example: + python negative_density_reference_model.py --self-test + python negative_density_reference_model.py --output /tmp/surface-control --png + +Only the optional PNG writer requires matplotlib; all arithmetic is stdlib. +""" + +from __future__ import annotations + +import argparse +import csv +import hashlib +import json +import math +from decimal import Decimal, localcontext +from fractions import Fraction +from pathlib import Path +import unittest + + +INDICES = (Fraction(1), Fraction(3, 2), Fraction(3)) +INVERSE_MOMENTS = ((9, -36, 30), (-36, 192, -180), (30, -180, 180)) +METHODS = ("exact_l2_projection", "gauss_legendre_interpolation") + + +def decimal(value): + if isinstance(value, Fraction): + return Decimal(value.numerator) / Decimal(value.denominator) + return Decimal(str(value)) + + +def nodes(): + r = (Decimal(3) / 5).sqrt() + return ((1 - r) / 2, Decimal("0.5"), (1 + r) / 2) + + +def exact_value(t, n): + return t * t.sqrt() if n == Fraction(3, 2) else t ** int(n) + + +def projection(n): + """Invert the degree-two Hilbert moment matrix with rational arithmetic.""" + moments = [1 / (n + j + 1) for j in range(3)] + return tuple(sum(a * b for a, b in zip(row, moments)) + for row in INVERSE_MOMENTS) + + +def interpolation(n): + x = nodes() + values = [exact_value(t, n) for t in x] + slope = (values[2] - values[0]) / (x[2] - x[0]) + curvature = (values[2] + values[0] - 2 * values[1]) * 10 / 3 + return (values[1] - slope / 2 + curvature / 4, + slope - curvature, curvature) + + +def coefficients(n, method): + if method == METHODS[0] or n in (1, 3): + # For n=3 the cubic nodal polynomial is the shifted Legendre P3; + # hence its remainder is also the exact L2 projection. Using that + # identity avoids classifying a rounded zero curvature as an n=1 root. + return tuple(decimal(c) for c in projection(n)) + return interpolation(n) + + +def value(c, t): + return c[0] + t * (c[1] + t * c[2]) + + +def integral(c, left=Decimal(0), right=Decimal(1)): + return sum(a * (right ** (j + 1) - left ** (j + 1)) / (j + 1) + for j, a in enumerate(c)) + + +def analyze(c): + a, b, q = c + if q == 0: + roots = [] if b == 0 else [-a / b] + else: + discriminant = b * b - 4 * a * q + roots = ([] if discriminant < 0 else + sorted({(-b - discriminant.sqrt()) / (2 * q), + (-b + discriminant.sqrt()) / (2 * q)})) + unit_roots = [t for t in roots if 0 <= t <= 1] + extrema = [Decimal(0), Decimal(1)] + if q != 0 and 0 < -b / (2 * q) < 1: + extrema.append(-b / (2 * q)) + minimum_at = min(extrema, key=lambda t: value(c, t)) + bounds = sorted({Decimal(0), Decimal(1), *unit_roots}) + negative = [(left, right) for left, right in zip(bounds, bounds[1:]) + if value(c, (left + right) / 2) < 0] + return { + "roots_in_cell": unit_roots, + "minimum_at": minimum_at, + "minimum": value(c, minimum_at), + "negative_intervals": negative, + "negative_width_fraction": sum((b - a for a, b in negative), Decimal(0)), + "negative_mass": -sum((integral(c, a, b) for a, b in negative), Decimal(0)), + "cell_mean": integral(c), + "boundary_value": a, + "interior_endpoint_value": value(c, Decimal(1)), + "nonnegative_on_closed_cell": value(c, minimum_at) >= 0, + } + + +def as_numbers(item): + if isinstance(item, (Decimal, Fraction)): + result = float(item) + if not math.isfinite(result): + raise ValueError("Nonfinite derived output") + return result + if isinstance(item, dict): + return {k: as_numbers(v) for k, v in item.items()} + if isinstance(item, (tuple, list)): + return [as_numbers(v) for v in item] + return item + + +def case(n, method): + c = coefficients(n, method) + result = analyze(c) + exact_mean = 1 / (decimal(n) + 1) + if n == 1: + exact_roots, exact_minimum = "0", "0" + elif n == 3: + exact_roots, exact_minimum = "1/5-sqrt(6)/30, 1/5+sqrt(6)/30", "-1/100" + elif method == METHODS[0]: + exact_roots, exact_minimum = "(sqrt(174)-12)/30", "-2/105" + else: + exact_roots = "(-c1+sqrt(c1^2-4*c0*c2))/(2*c2); c_j from the algebraic interpolation formula" + exact_minimum = "c0 from the algebraic interpolation formula" + result.update({ + "n": n, "method": method, "coefficients_c0_c1_c2": c, + "coefficients_decimal": [str(v) for v in c], + "exact_projection_coefficients": [str(v) for v in projection(n)], + "roots_in_cell_exact": exact_roots, "minimum_exact": exact_minimum, + "exact_cell_mean": exact_mean, + "mean_error": result["cell_mean"] - exact_mean, + "negative_mass_over_exact_cell_mass": result["negative_mass"] / exact_mean, + "values_at_gl_nodes": [value(c, t) for t in nodes()], + "exact_values_at_gl_nodes": [exact_value(t, n) for t in nodes()], + }) + return result + + +def refinement_row(result, width, amplitude): + n = result["n"] + scale = amplitude * exact_value(width, n) + return { + "n": n, "method": result["method"], "width": width, + "amplitude_A": amplitude, "minimum_density": scale * result["minimum"], + "minimum_distance": width * result["minimum_at"], + "negative_width": width * result["negative_width_fraction"], + "negative_mass_per_unit_area": scale * width * result["negative_mass"], + "exact_mass_per_unit_area": scale * width * result["exact_cell_mean"], + "negative_mass_over_exact_cell_mass": result["negative_mass_over_exact_cell_mass"], + } + + +class ReferenceChecks(unittest.TestCase): + def test_exact_moments(self): + for n in INDICES: + c = projection(n) + for j in range(3): + self.assertEqual(sum(a / (i + j + 1) for i, a in enumerate(c)), + 1 / (n + j + 1)) + + def test_known_coefficients(self): + self.assertEqual(projection(Fraction(1)), (0, 1, 0)) + self.assertEqual(projection(Fraction(3, 2)), + (Fraction(-2, 105), Fraction(16, 35), Fraction(4, 7))) + self.assertEqual(projection(Fraction(3)), + (Fraction(1, 20), Fraction(-3, 5), Fraction(3, 2))) + + def test_interpolation_and_integer_identity(self): + for n in INDICES: + c = interpolation(n) + for t in nodes(): + self.assertLess(abs(value(c, t) - exact_value(t, n)), Decimal("1e-55")) + if n in (1, 3): + for actual, expected in zip(c, coefficients(n, METHODS[1])): + self.assertLess(abs(actual - expected), Decimal("1e-55")) + + def test_positive_dofs_do_not_certify_polynomial(self): + for n in (Fraction(3, 2), Fraction(3)): + result = case(n, METHODS[1]) + self.assertTrue(all(v > 0 for v in result["values_at_gl_nodes"])) + self.assertLess(result["minimum"], 0) + + def test_roots_and_extrema(self): + for method in METHODS: + self.assertTrue(case(Fraction(1), method)["nonnegative_on_closed_cell"]) + for n in (Fraction(3, 2), Fraction(3)): + result = case(n, method) + c = coefficients(n, method) + for t in result["roots_in_cell"]: + self.assertLess(abs(value(c, t)), Decimal("1e-55")) + for left, right in result["negative_intervals"]: + self.assertLess(value(c, (left + right) / 2), 0) + self.assertGreater(result["negative_mass"], 0) + cubic = case(Fraction(3), METHODS[0]) + self.assertEqual(cubic["minimum_at"], Decimal("0.2")) + self.assertEqual(cubic["minimum"], Decimal("-0.01")) + self.assertGreater(cubic["boundary_value"], 0) + self.assertEqual(cubic["roots_in_cell"], + [Decimal(".2") - Decimal(6).sqrt() / 30, + Decimal(".2") + Decimal(6).sqrt() / 30]) + + def test_gauss_mean_and_nonconservation_of_interpolation(self): + weights = (Decimal(5) / 18, Decimal(4) / 9, Decimal(5) / 18) + for n in INDICES: + mean = sum(w * exact_value(t, n) for w, t in zip(weights, nodes())) + self.assertLess(abs(mean - integral(interpolation(n))), Decimal("1e-55")) + self.assertLess(case(Fraction(3, 2), METHODS[1])["mean_error"], Decimal("-1e-4")) + + def test_refinement_scales(self): + for n in (Fraction(3, 2), Fraction(3)): + for method in METHODS: + result = case(n, method) + coarse = refinement_row(result, Decimal(1), Decimal(7)) + fine = refinement_row(result, Decimal(".5"), Decimal(7)) + factor = exact_value(Decimal(".5"), n) + for key, expected in (("minimum_density", factor), + ("negative_width", Decimal(".5")), + ("negative_mass_per_unit_area", factor / 2)): + self.assertLess(abs(fine[key] / coarse[key] - expected), Decimal("1e-55")) + + +def write_csv(path, rows): + rows = [as_numbers(row) for row in rows] + with path.open("w", newline="") as stream: + writer = csv.DictWriter(stream, fieldnames=list(rows[0])) + writer.writeheader() + for row in rows: + writer.writerow({k: json.dumps(v) if isinstance(v, (list, dict)) else v + for k, v in row.items()}) + + +def plot(output, results, samples): + import matplotlib + matplotlib.use("Agg") + import matplotlib.pyplot as plt + + fig, axes = plt.subplots(2, 3, figsize=(12, 6), constrained_layout=True) + ts = [i / (samples - 1) for i in range(samples)] + for column, n in enumerate(INDICES): + for row in range(2): + ax = axes[row, column] + ax.plot(ts, [t ** float(n) for t in ts], "k", label="Exact density") + for result in results: + if result["n"] != n: + continue + c = [float(v) for v in result["coefficients_c0_c1_c2"]] + ax.plot(ts, [value(c, t) for t in ts], + "--" if result["method"] == METHODS[1] else "-", + label="GL interpolation" if result["method"] == METHODS[1] else "L2 projection") + ax.axhline(0, color="0.7", linewidth=0.8) + ax.scatter([float(t) for t in nodes()], + [float(exact_value(t, n)) for t in nodes()], s=18, color="k") + ax.set_title(f"n = {float(n):g}" + ("; near-surface detail" if row else "")) + ax.set_xlabel("Inward distance / cell width") + if row: + ax.set_xlim(0, 0.35) + ax.set_ylim(-0.027, 0.075 if n == 3 else 0.23) + if column == 0: + ax.set_ylabel("Density / (A width^n)") + axes[0, 0].legend(fontsize=8) + fig.suptitle("Uniform 1D surface-cell model: positive nodal data can yield a negative polynomial") + fig.savefig(output / "profiles.png", dpi=170) + plt.close(fig) + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--output", type=Path, help="Fresh directory for JSON and CSV artifacts") + parser.add_argument("--self-test", action="store_true") + parser.add_argument("--png", action="store_true", help="Also plot profiles (requires matplotlib)") + parser.add_argument("--widths", nargs="+", default=["1", ".5", ".25", ".125"]) + parser.add_argument("--amplitude", default="1", help="A in rho=A*s^n") + parser.add_argument("--samples", type=int, default=501, help="Plot/CSV samples, not extrema checks") + args = parser.parse_args() + widths, amplitude = [Decimal(x) for x in args.widths], Decimal(args.amplitude) + if any(not x.is_finite() or x <= 0 for x in (*widths, amplitude)): + parser.error("Widths and amplitude must be finite and positive") + if args.samples < 2 or (args.png and args.output is None): + parser.error("At least two samples are required; --png requires --output") + with localcontext() as context: + context.prec = 70 + if args.self_test: + test_result = unittest.TextTestRunner(verbosity=2).run( + unittest.defaultTestLoader.loadTestsFromTestCase(ReferenceChecks)) + if not test_result.wasSuccessful(): + raise SystemExit(1) + if args.output is None: + return + results = [case(n, method) for n in INDICES for method in METHODS] + report = { + "schema": "negative_density_reference_model_v1", + "producer_sha256": hashlib.sha256(Path(__file__).read_bytes()).hexdigest(), + "model": "h=kappa*s; rho=A*s^n; t=s/width in [0,1]; uniform physical weight", + "scope": "Analytic mechanism control, not a fit or accuracy certificate for the 3D star", + "arithmetic": "Exact rational projection; 70-digit decimal algebraic interpolation/roots", + "gl_nodes": nodes(), "gl_nodes_exact": "(1-sqrt(3/5))/2, 1/2, (1+sqrt(3/5))/2", + "node_source": "MFEM L2_FECollection defaults to BasisType::GaussLegendre; Q2 tensor element", + "projection_derivation": "c=H^-1 b; H_ij=1/(i+j+1); b_i=1/(n+i+1), i,j=0,1,2", + "interpolation_derivation": "v_j=t_j^n; b=(v_+-v_-)/sqrt(3/5); c=(10/3)*(v_++v_--2v_0); p=v_0+b*(t-.5)+c*(t-.5)^2", + "negative_interval_convention": "Interior points are strictly negative; endpoints are roots except a negative physical-cell endpoint", + "refinement": "rho=A*width^n*p(t); negative mass/area scales width^(n+1); relative first-cell negativity need not disappear", + "cases": results, + } + if args.output is None: + print(json.dumps(as_numbers(report), indent=2, allow_nan=False)) + return + args.output.mkdir(parents=True, exist_ok=False) + write_csv(args.output / "summary.csv", results) + write_csv(args.output / "refinement.csv", + [refinement_row(result, width, amplitude) for result in results for width in widths]) + profiles = [] + for result in results: + for i in range(args.samples): + t = Decimal(i) / (args.samples - 1) + profiles.append({"n": result["n"], "method": result["method"], "t": t, + "exact": exact_value(t, result["n"]), + "approximation": value(result["coefficients_c0_c1_c2"], t)}) + write_csv(args.output / "profiles.csv", profiles) + if args.png: + plot(args.output, results, args.samples) + report["artifacts_sha256"] = {p.name: hashlib.sha256(p.read_bytes()).hexdigest() + for p in sorted(args.output.iterdir())} + (args.output / "summary.json").write_text(json.dumps(as_numbers(report), indent=2, allow_nan=False) + "\n") + print(f"Wrote six analytic controls to {args.output}") + + +if __name__ == "__main__": + main() diff --git a/experiments/plateau_diagnostic.cpp b/experiments/plateau_diagnostic.cpp index eff07a9..cfd62a6 100644 --- a/experiments/plateau_diagnostic.cpp +++ b/experiments/plateau_diagnostic.cpp @@ -33,6 +33,7 @@ import mean_field; #include "shape_field_weak_diagnostics.hpp" #include "reduced_surface_diagnostics.hpp" #include "spherical_harmonic_surface_diagnostics.hpp" +#include "negative_density_diagnostics.hpp" using namespace mean_field; namespace { @@ -45,6 +46,7 @@ struct Options { bool preconditionerProbe = false, nonlinearStepProbe = false; bool weakShapeProbe = false, weakShapeFinal = false; bool reducedSurfaceProbe = false, reducedSurfaceAlgebraCheck = false; + bool densityDiagnostic = false; int surfaceDegree = 4, surfaceMaximumOrder = -1; std::string surfacePrescription = "nodal"; bool coupledTrialProbe = false, coupledSeedControl = false; @@ -64,6 +66,7 @@ struct Options { std::vector captureIterations; std::string mesh = "sandbox.smesh"; std::filesystem::path output, probeState, probeDirection, replayCandidate; + std::filesystem::path densitySnapshot; }; Options Parse(int argc, char **argv) { Options o; @@ -124,6 +127,10 @@ Options Parse(int argc, char **argv) { o.reducedSurfaceProbe = true; else if (a == "--reduced-surface-algebra-check") o.reducedSurfaceAlgebraCheck = true; + else if (a == "--density-diagnostic") + o.densityDiagnostic = true; + else if (a == "--density-snapshot") + o.densitySnapshot = v(); else if (a == "--surface-prescription") o.surfacePrescription = v(); else if (a == "--surface-maximum-order") @@ -194,6 +201,15 @@ Options Parse(int argc, char **argv) { throw std::invalid_argument("Invalid options"); if (o.surfacePrescription != "nodal" && o.surfacePrescription != "harmonics") throw std::invalid_argument("Surface prescriptions are nodal or harmonics"); + if (!o.densitySnapshot.empty() && !o.densityDiagnostic) + throw std::invalid_argument("--density-snapshot requires --density-diagnostic"); + if (o.densityDiagnostic && (o.maxNewton != 0 || o.strategy != "baseline" || o.probes || + o.pressureProbe || o.shapeCurvature || o.surfaceTrace || o.strategyProbe || o.positivityProbe || + o.positivityCertificate || o.preconditionerProbe || o.nonlinearStepProbe || o.reducedSurfaceProbe || + o.reducedSurfaceAlgebraCheck || o.weakShapeProbe || o.weakShapeFinal || o.coupledTrialProbe || + o.dampedCoupledProbe || !o.probeState.empty() || !o.probeDirection.empty() || + !o.replayCandidate.empty() || !o.captureIterations.empty())) + throw std::invalid_argument("Density diagnostics require --max-newton 0 and no solver/probe strategies"); if (o.surfacePrescription == "harmonics") { deformation::SphericalHarmonicSurfaceOptions surfaceOptions{.maximumDegree = o.surfaceDegree}; if (o.surfaceMaximumOrder != -1) surfaceOptions.maximumAzimuthalOrder = o.surfaceMaximumOrder; @@ -372,6 +388,56 @@ mfem::Vector LoadVector(const std::filesystem::path &path, int expectedSize) { return result; } +template +mfem::Vector LoadDensitySnapshot(const Problem &problem, const Options &o, const Scales &scales) { + std::ifstream metadata(o.densitySnapshot / "metadata.txt"); + if (!metadata) throw std::invalid_argument("Density replay requires snapshot metadata."); + std::map saved; + std::string line; + while (std::getline(metadata, line)) { + const auto separator = line.find('='); + if (separator != std::string::npos) saved.emplace(line.substr(0, separator), line.substr(separator + 1)); + } + const auto number = [&](const std::string &key, double expected) { + const double actual = std::stod(saved.at(key)); + if (!std::isfinite(actual) || std::abs(actual - expected) > 1e-13 * std::max(1., std::abs(expected))) + throw std::invalid_argument("Density snapshot model mismatch: " + key); + }; + number("n", o.n); number("K", scales.K); number("J", o.angularMomentum); + number("M", 1.); number("R_reference", 1.); number("central_phase", !o.noCentral); + if (saved.at("surface_prescription") != o.surfacePrescription) + throw std::invalid_argument("Density snapshot surface prescription mismatch."); + if (o.surfacePrescription == "harmonics") { + number("surface_degree", o.surfaceDegree); + number("surface_maximum_order", o.surfaceMaximumOrder); + } + // A matching vector size does not establish matching coordinates. Require + // every named state and residual block to match the current manifest. + std::map> expected; + for (const auto &block : problem.GetManifest().valueBlocks()) + expected.emplace("value," + std::string(block.stableId), std::pair{block.offset, block.size}); + for (const auto &block : problem.GetManifest().residualBlocks()) + expected.emplace("residual," + std::string(block.stableId), std::pair{block.offset, block.size}); + std::ifstream layout(o.densitySnapshot / "layout.csv"); + if (!std::getline(layout, line) || line != "kind,block,offset,size,scale") + throw std::invalid_argument("Density replay requires the saved block layout."); + while (std::getline(layout, line)) { + std::istringstream row(line); + std::array fields; + for (auto &field : fields) if (!std::getline(row, field, ',')) + throw std::invalid_argument("Malformed saved density layout."); + const auto key = fields[0] + "," + fields[1]; + const auto found = expected.find(key); + if (found == expected.end() || found->second != std::pair{std::stoi(fields[2]), std::stoi(fields[3])}) + throw std::invalid_argument("Density snapshot coordinate layout mismatch: " + key); + expected.erase(found); + } + if (!expected.empty() || layout.bad() || metadata.bad()) + throw std::invalid_argument("Incomplete density snapshot provenance."); + return experiment::plateau_diagnostics::coupled_physical_detail::Load( + o.densitySnapshot / "accepted_state.txt", problem.StateSize()); +} + template void FocusedProbes(State &s, const fem::FEM &f, const Options &o, const mfem::Vector &physicalBase) { mfem::Vector base; @@ -615,6 +681,60 @@ int RunContext(Context &context, const Options &o, const Scales &scales) { << "\ncontext_gravity_schur_relaxation=" << (contextUsesGravitySchur ? o.schurRelaxation : 0.0) << "\nschur_correction_sweeps=" << o.schurCorrectionSweeps << "\nschur_relaxation=" << o.schurRelaxation << "\ncompiler=" << __VERSION__ << "\ncompiled=" << __DATE__ << ' ' << __TIME__ << '\n'; + if (o.densityDiagnostic) { + metadata << "density_diagnostic=1\ndensity_source=" + << (o.densitySnapshot.empty() ? "lane_emden_seed" : o.densitySnapshot.string()) << '\n'; + metadata.flush(); + solver::detail::StellarEquilibriumContextDiagnostics::WithState(context, [&](auto &s, const fem::FEM &f) { + namespace d = experiment::plateau_diagnostics::coupled_physical_detail; + const mfem::Vector accepted(s.AcceptedPhysicalState()), normalized(s.acceptedNormalizedState), + residual(s.acceptedNormalizedResidual); + const auto snapshot = o.densitySnapshot.empty() ? mfem::Vector(accepted) : LoadDensitySnapshot(s.Problem(), o, scales); + s.normalizedOperator->NormalizeState(snapshot, s.trialNormalizedState); + if (!s.PrepareTrial()) throw std::runtime_error("Density snapshot cannot be prepared by the physical operator."); + d::PreparedSnapshotView view(s); + experiment::polytrope_validation::PhysicalState physical(view, f); + std::cout << "Frozen density diagnostic F=" << s.trialNormalizedResidual.Norml2() << std::endl; + SaveVector(o.output / "observed_state.txt", s.normalizedOperator->GetPhysicalState()); + auto layout = File(o.output / "layout.csv"); + layout << "kind,block,offset,size,scale\n"; + for (const auto &block : s.Problem().GetManifest().valueBlocks()) + layout << "value," << block.stableId << ',' << block.offset << ',' << block.size << ',' << block.scale << '\n'; + for (const auto &block : s.Problem().GetManifest().residualBlocks()) + layout << "residual," << block.stableId << ',' << block.offset << ',' << block.size << ',' << block.scale << '\n'; + auto physicalOutput = File(o.output / "density_physical.csv"); + physicalOutput << "quadrature,metric,value\n"; + for (const auto &[key, value] : experiment::plateau_diagnostics::MeasurePhysical(physical, o.n, scales.K, 18)) + physicalOutput << "18," << key << ',' << value << '\n'; + physicalOutput.flush(); + const auto &closure = s.Problem().GetPhysicalOperator().GetBarotropicClosureOperator(); + auto ruleOutput = File(o.output / "production_closure_rules.csv"); + ruleOutput << "element,quadrature_order,quadrature_points\n"; + closure.VisitMappedGeometryRules([&](int element, const mfem::IntegrationRule &rule) { + ruleOutput << element << ',' << rule.GetOrder() << ',' << rule.GetNPoints() << '\n'; + }); + mfem::Vector closureResidual; + closure.BuildResidual(closureResidual); + mfem::ParGridFunction closureFull(f.densityFes.get()); + const auto adapter = field::make_field_dof_grid_function_adapter(*f.densityFes); + adapter.scatter(closureResidual, closureFull); + SaveVector(o.output / "production_closure_residual_full.txt", closureFull); + SaveVector(o.output / "observed_density_coefficients.txt", physical.density); + experiment::plateau_diagnostics::RecordNegativeDensity(physical, o.n, scales.K, o.output); + if (!d::Identical(accepted, s.AcceptedPhysicalState()) || !d::Identical(normalized, s.acceptedNormalizedState) || + !d::Identical(residual, s.acceptedNormalizedResidual)) + throw std::runtime_error("Density diagnostic changed the accepted state."); + s.RestoreAccepted(); + auto invariants = File(o.output / "density_invariants.csv"); + invariants << "metric,value\naccepted_physical_unchanged," << d::Identical(accepted, s.AcceptedPhysicalState()) + << "\naccepted_normalized_unchanged," << d::Identical(normalized, s.acceptedNormalizedState) + << "\naccepted_residual_unchanged," << d::Identical(residual, s.acceptedNormalizedResidual) + << "\nprepared_residual_norm," << physical.normalizedBorderedResidualNorm << '\n'; + if (!d::Identical(accepted, s.AcceptedPhysicalState())) throw std::runtime_error("Density restoration failed."); + }); + return 0; + } if (o.weakShapeProbe || o.weakShapeFinal) metadata << "weak_shape_direction_space=normalized_current_factors\n" << "weak_shape_operator=full_monolithic\n" diff --git a/experiments/run_plateau_campaign.py b/experiments/run_plateau_campaign.py index 9152926..6b580e2 100644 --- a/experiments/run_plateau_campaign.py +++ b/experiments/run_plateau_campaign.py @@ -132,6 +132,26 @@ def main(): path = source / path path = path.resolve() state_files[str(path)] = digest(path) + elif word == '--density-snapshot': + path = Path(case['arguments'][i + 1]) + path = (path if path.is_absolute() else source / path).resolve() + # Modal coordinates also depend on the mesh and named block + # layout. Bind replay to its original recorded campaign rather + # than trusting a coincidentally matching vector length. + prior_manifest = path.parent / 'provenance.json' + prior_run = path.parent / (path.name + '_run.json') + prior = json.loads(prior_manifest.read_text()) + saved_run = json.loads(prior_run.read_text()) + if prior['mesh_sha256'] != provenance['mesh_sha256']: + raise RuntimeError('Density replay mesh differs from the snapshot campaign.') + for filename in ('accepted_state.txt', 'metadata.txt', 'layout.csv'): + snapshot_file = path / filename + value = digest(snapshot_file) + if value != saved_run['output_hashes'][filename]: + raise RuntimeError('Density replay input no longer matches its original run: ' + filename) + state_files[str(snapshot_file)] = value + state_files[str(prior_manifest)] = digest(prior_manifest) + state_files[str(prior_run)] = digest(prior_run) print('START', name, flush=True) start = time.monotonic() cpu_start = resource.getrusage(resource.RUSAGE_CHILDREN) @@ -161,12 +181,15 @@ def main(): for path in sorted(folder.rglob('*')) if path.is_file() and path.suffix in ('.txt', '.csv', '.json') } + record['inputs_unchanged'] = all(digest(Path(path)) == value for path, value in state_files.items()) write_json(output / (name + '_run.json'), record) summaries.append({**record, **summarize(folder)}) write_json(output / 'summary.json', summaries) print('FINISH', name, record['exit_code'], round(record['seconds'], 2), flush=True) if record['exit_code'] not in (0, 3): raise RuntimeError('Case failed to execute; see its retained log and run record.') + if not record['inputs_unchanged']: + raise RuntimeError('A frozen input changed during the diagnostic run.') # Exit 0 denotes completed observations; exit 3 from a model remains explicitly recorded as nonconvergence. return 0