Add frozen negative density localization and independent EOS projection diagnostics
This commit is contained in:
230
experiments/audit_negative_density.py
Normal file
230
experiments/audit_negative_density.py
Normal file
@@ -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()
|
||||
18
experiments/negative_density_cases_2026-09-19.json
Normal file
18
experiments/negative_density_cases_2026-09-19.json
Normal file
@@ -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"]
|
||||
}
|
||||
]
|
||||
362
experiments/negative_density_diagnostics.hpp
Normal file
362
experiments/negative_density_diagnostics.hpp
Normal file
@@ -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 <array>
|
||||
#include <cmath>
|
||||
#include <filesystem>
|
||||
#include <fstream>
|
||||
#include <iomanip>
|
||||
#include <iostream>
|
||||
#include <limits>
|
||||
#include <map>
|
||||
#include <queue>
|
||||
#include <stdexcept>
|
||||
#include <string>
|
||||
#include <vector>
|
||||
|
||||
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<SurfaceFace> surface;
|
||||
std::vector<int> layer;
|
||||
std::vector<std::vector<int>> 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<std::vector<int>> neighbors(mesh.GetNE());
|
||||
std::queue<int> 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<double, 3> 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<int>(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<double>::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<double, 3> 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<int> 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<std::string> 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<d::Aggregate> global(labels.size());
|
||||
std::vector<std::map<int, d::Aggregate>> 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<d::Aggregate> 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<double> densities{rho, rhoH};
|
||||
for (const auto &control : projections.controls)
|
||||
densities.push_back(control.density->GetValue(element, point));
|
||||
const int bin = static_cast<int>(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<double, 5> tangents{0., .25, .5, .75, 1.};
|
||||
constexpr std::array<double, 17> 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<double, 3> 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
|
||||
396
experiments/negative_density_projection.hpp
Normal file
396
experiments/negative_density_projection.hpp
Normal file
@@ -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 <algorithm>
|
||||
#include <array>
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
#include <filesystem>
|
||||
#include <fstream>
|
||||
#include <iomanip>
|
||||
#include <limits>
|
||||
#include <memory>
|
||||
#include <stdexcept>
|
||||
#include <string>
|
||||
#include <utility>
|
||||
#include <vector>
|
||||
|
||||
#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<mfem::ParGridFunction> 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<DensityProjectionControl> 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<long double>(vector(i)) * vector(i);
|
||||
}
|
||||
const double norm = static_cast<double>(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<double>::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<long double>(lower(i, k)) * lower(j, k);
|
||||
if (i == j) {
|
||||
if (!std::isfinite(value) || value <= 64.0L * std::numeric_limits<double>::epsilon() * scale)
|
||||
throw std::domain_error("Density projection mass is singular or numerically unresolved.");
|
||||
lower(i, j) = static_cast<double>(std::sqrt(value));
|
||||
} else {
|
||||
lower(i, j) = static_cast<double>(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<long double>(lower(i, j)) * result(j);
|
||||
result(i) = static_cast<double>(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<long double>(lower(j, i)) * result(j);
|
||||
result(i) = static_cast<double>(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<long double>(rhs(i));
|
||||
for (int j = 0; j < coefficients.Size(); ++j)
|
||||
value += static_cast<long double>(mass(i, j)) * coefficients(j);
|
||||
result(i) = static_cast<double>(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<long double>(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<std::size_t>(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<double>::infinity()};
|
||||
double minimumProjectedCoefficient{std::numeric_limits<double>::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<const mfem::L2_FECollection*>(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<double>(field::Enthalpy::Scalar::familyOrder);
|
||||
if (!std::isfinite(extra) || extra < 0.0 || extra > std::numeric_limits<int>::max())
|
||||
throw std::invalid_argument("Invalid production EOS quadrature order.");
|
||||
const int extraOrder = static_cast<int>(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<std::pair<const char*, int>, 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<mfem::ParGridFunction>(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<detail::Totals, 3> totals;
|
||||
std::vector<bool> visited(static_cast<std::size_t>(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<int> 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<std::size_t>(dofs[i])])
|
||||
throw std::invalid_argument("Density projection requires independent, unsigned element-local DOFs.");
|
||||
visited[static_cast<std::size_t>(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<field::Enthalpy>::make_query<field::Enthalpy::Form::EosClosureSource>(
|
||||
quadrature::QuadratureRole::discretization, transformation->OrderW(), std::array<int, 1>{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<double>::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<eos::quantity::Density>(
|
||||
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<double>::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<double>::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<long double>(samples(q, i)) * coefficients(i);
|
||||
savedValue += static_cast<long double>(samples(q, i)) * saved(i);
|
||||
constantValue += static_cast<long double>(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<double>(projectedValue));
|
||||
minimumSavedSample = std::min(minimumSavedSample, static_cast<double>(savedValue));
|
||||
constantSampleError = std::max(constantSampleError, static_cast<double>(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<long double>(projectedDefectNorm) * projectedDefectNorm;
|
||||
total.savedDefectSquared += static_cast<long double>(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
|
||||
347
experiments/negative_density_reference_model.py
Normal file
347
experiments/negative_density_reference_model.py
Normal file
@@ -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()
|
||||
@@ -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<int> 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 <class Problem>
|
||||
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<std::string, std::string> 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<std::string, std::pair<int, int>> 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<std::string, 5> 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 <class State>
|
||||
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<field::Density,
|
||||
utils::domain::CoreEnvelopeVacuumDomainSchema>(*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"
|
||||
|
||||
@@ -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
|
||||
|
||||
|
||||
Reference in New Issue
Block a user