From d1a594808f960bbaeeb6df8a8ce6a95458765b19 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Sat, 19 Sep 2026 10:06:10 -0400 Subject: [PATCH] Integrate frozen negative density over polynomial root intervals --- experiments/audit_negative_density.py | 13 + ...ative_density_region_cases_2026-09-19.json | 73 +++ experiments/negative_density_regions.hpp | 420 ++++++++++++++++++ experiments/plateau_diagnostic.cpp | 13 +- 4 files changed, 517 insertions(+), 2 deletions(-) create mode 100644 experiments/negative_density_region_cases_2026-09-19.json create mode 100644 experiments/negative_density_regions.hpp diff --git a/experiments/audit_negative_density.py b/experiments/audit_negative_density.py index 4190a28..1abaec8 100644 --- a/experiments/audit_negative_density.py +++ b/experiments/audit_negative_density.py @@ -40,6 +40,16 @@ def close(left, right, label, rtol=2e-10, atol=1e-15): f"{label}: {left} != {right}") +def require_complete_hashes(record): + required = {'metadata.txt', 'density_summary.csv', 'density_elements.csv', 'density_radial_bins.csv', + 'density_physical.csv', 'negative_density_samples.csv', 'negative_density_projection_elements.csv', + 'production_closure_rules.csv', 'production_closure_residual_full.txt', + 'projection_production_closure_residual_full.txt', 'observed_density_coefficients.txt', + 'density_invariants.csv', 'density_element_topology.csv', 'density_surface_traces.csv'} + required.update(f'negative_density_projection_{control}_coefficients.txt' for control in ('production', 'q18', 'q26')) + require(required <= record['output_hashes'].keys(), 'Run record omits a consumed artifact hash.') + + def clean(value): if isinstance(value, dict): return {str(k): clean(v) for k, v in value.items()} @@ -121,6 +131,8 @@ def audit_case(folder): 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') + close(float(np.linalg.norm(delta)), math.sqrt(math.fsum(float(r['coefficient_difference_norm']) ** 2 for r in matching)), + 'coefficient difference norm assembled from disjoint elements') for row in matching: require(float(row['constant_coefficient_error_maximum']) < 1e-10 and float(row['constant_sample_error_maximum']) < 1e-10, @@ -212,6 +224,7 @@ def main(): 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}') + require_complete_hashes(run) 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(): diff --git a/experiments/negative_density_region_cases_2026-09-19.json b/experiments/negative_density_region_cases_2026-09-19.json new file mode 100644 index 0000000..35c5be2 --- /dev/null +++ b/experiments/negative_density_region_cases_2026-09-19.json @@ -0,0 +1,73 @@ +[ + { + "id": "n1_equilibrium", + "arguments": [ + "--index", + "1", + "--angular-momentum", + "0", + "--surface-prescription", + "harmonics", + "--surface-degree", + "4", + "--no-central", + "--max-newton", + "0", + "--density-regions", + "--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-regions" + ] + }, + { + "id": "n1p5_equilibrium", + "arguments": [ + "--index", + "1.5", + "--angular-momentum", + "0", + "--surface-prescription", + "harmonics", + "--surface-degree", + "4", + "--no-central", + "--max-newton", + "0", + "--density-regions", + "--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-regions" + ] + } +] diff --git a/experiments/negative_density_regions.hpp b/experiments/negative_density_regions.hpp new file mode 100644 index 0000000..3bcd0c2 --- /dev/null +++ b/experiments/negative_density_regions.hpp @@ -0,0 +1,420 @@ +#pragma once + +// Include after import mean_field. Root splitting removes the discontinuous +// sign indicator from the integration coordinate; tangential and mapped-volume +// integration still require the quadrature convergence study recorded here. +#include "polytrope_physical_state.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace experiment::plateau_diagnostics { +namespace negative_density_regions_detail { + +using Real = long double; + +struct Polynomial final { + Real c, b, a; + [[nodiscard]] Real Value(Real t) const { return std::fma(std::fma(a, t, b), t, c); } + [[nodiscard]] Real Integral(Real left, Real right) const { + const Real width = right - left; + // Subtracting antiderivatives loses a thin pocket's entire mass when + // both antiderivatives are large. This centered quadratic identity does + // not subtract nearly equal endpoint values. + return width * (Value(left + width / 2) + a * width * width / 12); + } + [[nodiscard]] Real Scale() const { return std::max({std::abs(a), std::abs(b), std::abs(c)}); } +}; + +struct Interval final { Real left, right; bool negative; }; + +inline Polynomial Interpolate(Real left, Real middle, Real right) { + return {left, 4 * middle - 3 * left - right, 2 * (left + right - 2 * middle)}; +} + +inline std::vector Partition(const Polynomial& p) { + if (!std::isfinite(p.a) || !std::isfinite(p.b) || !std::isfinite(p.c)) + throw std::domain_error("Nonfinite density slice polynomial."); + std::vector points{0, 1}; + const auto add = [&](Real root) { + if (std::isfinite(root) && root > 0 && root < 1) points.push_back(root); + }; + const Real scale = p.Scale(); + if (scale > 0) { + const Real a = p.a / scale, b = p.b / scale, c = p.c / scale; + // No relative cutoff turns a small quadratic into a line: a thin + // negative interval can be the very effect this diagnostic must retain. + if (a == 0) { + if (b != 0) add(-c / b); + } else { + const Real discriminant = std::fma(-4 * a, c, b * b); + if (discriminant == 0) { + add(-b / (2 * a)); + } else if (discriminant > 0) { + const Real q = -Real{0.5} * (b + std::copysign(std::sqrt(discriminant), b)); + // The product-of-roots form avoids cancellation at a root near + // the face. q cannot vanish for a positive discriminant. + add(q / a); + add(c / q); + } + } + } + std::sort(points.begin(), points.end()); + points.erase(std::unique(points.begin(), points.end()), points.end()); + std::vector result; + for (std::size_t i = 1; i < points.size(); ++i) { + const Real left = points[i - 1], right = points[i]; + if (!(right > left)) throw std::logic_error("Invalid density root partition."); + const Real value = p.Value(left + (right - left) / 2); + const bool negative = value < 0; + // Two extra interior signs catch missed roots without assigning a + // finite-width pocket to an exactly tangent double root. + for (const Real fraction : {Real{0.25}, Real{0.75}}) { + const Real check = p.Value(left + fraction * (right - left)); + if ((negative && check > 0) || (!negative && check < 0)) + throw std::runtime_error("Density slice root partition has inconsistent interval signs."); + } + result.push_back({left, right, negative}); + } + return result; +} + +inline std::ofstream Open(const std::filesystem::path& path) { + if (std::filesystem::exists(path)) throw std::invalid_argument("Density region output already exists: " + path.string()); + std::ofstream out(path); + out.exceptions(std::ios::failbit | std::ios::badbit); + out << std::setprecision(20); + return out; +} + +// These are analytic, uniform-measure one-dimensional controls, independent of +// FEM and geometry. The fractional/cubic data are interpolated at three Gauss +// nodes, reproducing the underresolved positive-profile mechanism of interest. +inline void SelfChecks(const std::filesystem::path& path) { + auto out = Open(path); + out << "case,negative_length,expected_negative_length,negative_mass,expected_negative_mass,passed\n"; + const auto check = [&](const char* name, const Polynomial& p, Real expectedLength, Real expectedMass) { + Real length = 0, mass = 0; + for (const auto& interval : Partition(p)) if (interval.negative) { + length += interval.right - interval.left; + mass -= p.Integral(interval.left, interval.right); + } + const bool passed = std::abs(length - expectedLength) <= 2e-14L * std::max(Real{1e-12}, expectedLength) && + std::abs(mass - expectedMass) <= 2e-13L * std::max(Real{1e-12}, expectedMass); + out << name << ',' << length << ',' << expectedLength << ',' << mass << ',' << expectedMass << ',' << passed << '\n'; + if (!passed) throw std::runtime_error(std::string("Density region analytic self-check failed: ") + name); + }; + check("positive_quadratic", {1, 0, 1}, 0, 0); + check("negative_constant", {-2, 0, 0}, 1, 2); + check("linear_crossing", {-0.25L, 1, 0}, 0.25L, 0.03125L); + check("double_root_nonnegative", {0.25L, -1, 1}, 0, 0); + check("near_linear_positive", {1, 1, 1e-28L}, 0, 0); + // Exact binary coefficients isolate the root algorithm from the separate + // conditioning of subtracting a tiny decimal offset from a large value. + check("near_double_negative", {0.25L - std::ldexp(Real{1}, -40), -1, 1}, + std::ldexp(Real{1}, -19), 4 * std::ldexp(Real{1}, -60) / 3); + const Real d = std::sqrt(Real{15}) / 10; + const auto fractional = [](Real t) { return t * std::sqrt(t); }; + const Real y0 = fractional(Real{0.5} - d), y1 = fractional(Real{0.5}), y2 = fractional(Real{0.5} + d); + const Real a = (y0 + y2 - 2 * y1) / (2 * d * d); + const Real b = (y2 - y0) / (2 * d) - a; + const Polynomial fractionalInterpolant{y1 - b / 2 - a / 4, b, a}; + check("gauss_q2_t_power_1p5", fractionalInterpolant, + 0.04548403175047501373888L, 0.0005262295634284655L); + check("gauss_q2_t_power_3", {0.05L, -0.6L, 1.5L}, std::sqrt(Real{6}) / 15, std::sqrt(Real{6}) / 2250); +} + +struct Axis final { + int coordinate{0}; + int orientation{1}; + int face{-1}; + int faceCount{0}; +}; + +inline std::vector ChooseAxes(polytrope_validation::PhysicalState& state) { + auto& mesh = *state.finiteElements.mesh; + std::vector axes(static_cast(mesh.GetNE())); + 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) continue; + const int element = leftStar ? left : right; + auto& axis = axes[static_cast(element)]; + ++axis.faceCount; + if (axis.face >= 0) continue; + auto* transformation = mesh.GetFaceElementTransformations(face); + if (!transformation) throw std::runtime_error("Missing density region face 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 coordinates{point.x, point.y, point.z}; + int selected = -1; + for (int d = 0; d < 3; ++d) { + if (std::abs(coordinates[d]) < 1e-12 || std::abs(coordinates[d] - 1) < 1e-12) { + if (selected >= 0) throw std::runtime_error("Ambiguous density region reference face."); + selected = d; + } + } + if (selected < 0) throw std::runtime_error("Density region face is not a tensor-product cell face."); + axis.coordinate = selected; + axis.orientation = coordinates[selected] < 0.5 ? 1 : -1; + axis.face = face; + } + return axes; +} + +inline mfem::IntegrationPoint Point(const Axis& axis, double u, double v, Real t) { + std::array coordinate{}; + coordinate[axis.coordinate] = static_cast(axis.orientation > 0 ? t : 1 - t); + coordinate[(axis.coordinate + 1) % 3] = u; + coordinate[(axis.coordinate + 2) % 3] = v; + mfem::IntegrationPoint result; + result.Set3(coordinate[0], coordinate[1], coordinate[2]); + return result; +} + +struct Measure final { + Real volume{0}, mass{0}, negativeVolume{0}, negativeMass{0}; + std::size_t lines{0}, negativeIntervals{0}, points{0}; + double minimumRho{std::numeric_limits::infinity()}; + double witnessH{std::numeric_limits::quiet_NaN()}; + double witnessT{std::numeric_limits::quiet_NaN()}, witnessU{witnessT}, witnessV{witnessT}; + std::array witnessPosition{}; + int witnessElement{-1}; + double negativeRadiusMin{std::numeric_limits::infinity()}, negativeRadiusMax{-negativeRadiusMin}; + Real negativeTMin{std::numeric_limits::infinity()}, negativeTMax{-negativeTMin}; + double maximumReconstructionError{0}, maximumReconstructionTolerance{0}; + + void Add(const Measure& other) { + volume += other.volume; mass += other.mass; + negativeVolume += other.negativeVolume; negativeMass += other.negativeMass; + lines += other.lines; negativeIntervals += other.negativeIntervals; points += other.points; + if (other.minimumRho < minimumRho) { + minimumRho = other.minimumRho; witnessH = other.witnessH; + witnessT = other.witnessT; witnessU = other.witnessU; witnessV = other.witnessV; + witnessPosition = other.witnessPosition; witnessElement = other.witnessElement; + } + negativeRadiusMin = std::min(negativeRadiusMin, other.negativeRadiusMin); + negativeRadiusMax = std::max(negativeRadiusMax, other.negativeRadiusMax); + negativeTMin = std::min(negativeTMin, other.negativeTMin); + negativeTMax = std::max(negativeTMax, other.negativeTMax); + maximumReconstructionError = std::max(maximumReconstructionError, other.maximumReconstructionError); + maximumReconstructionTolerance = std::max(maximumReconstructionTolerance, other.maximumReconstructionTolerance); + } +}; + +inline Real Observed(Real value) { return std::isfinite(value) ? value : std::numeric_limits::quiet_NaN(); } + +inline void Header(std::ostream& out) { + out << "tangent_lines,negative_intervals,integration_points,volume,mass,negative_volume,negative_mass," + "negative_volume_fraction,negative_mass_fraction,minimum_density,witness_h,witness_element,witness_t,witness_u,witness_v," + "witness_x,witness_y,witness_z,negative_radius_minimum,negative_radius_maximum,negative_interval_t_minimum," + "negative_interval_t_maximum,maximum_reconstruction_error,maximum_reconstruction_tolerance\n"; +} + +inline void Write(std::ostream& out, const Measure& m) { + out << m.lines << ',' << m.negativeIntervals << ',' << m.points << ',' << m.volume << ',' << m.mass << ',' + << m.negativeVolume << ',' << m.negativeMass << ',' << m.negativeVolume / m.volume << ',' + << (m.mass != 0 ? m.negativeMass / std::abs(m.mass) : std::numeric_limits::quiet_NaN()) << ',' + << m.minimumRho << ',' << m.witnessH << ',' << m.witnessElement << ',' << m.witnessT << ',' << m.witnessU << ',' << m.witnessV; + for (double x : m.witnessPosition) out << ',' << x; + out << ',' << Observed(m.negativeRadiusMin) << ',' << Observed(m.negativeRadiusMax) << ',' + << Observed(m.negativeTMin) << ',' << Observed(m.negativeTMax) << ',' + << m.maximumReconstructionError << ',' << m.maximumReconstructionTolerance << '\n'; +} + +inline bool Identical(const mfem::Vector& a, const mfem::Vector& b) { + return a.Size() == b.Size() && (a.Size() == 0 || std::memcmp(a.GetData(), b.GetData(), sizeof(mfem::real_t) * a.Size()) == 0); +} + +} // namespace negative_density_regions_detail + +inline void RecordDensityNegativeRegions(polytrope_validation::PhysicalState& state, double n, double K, + const std::filesystem::path& output) { + namespace d = negative_density_regions_detail; + using d::Real; + if (!std::isfinite(n) || n < 1 || !std::isfinite(K) || K <= 0) + throw std::invalid_argument("Density regions require finite n>=1 and K>0."); + std::filesystem::create_directories(output); + // Fail the one-dimensional arithmetic controls before any FEM evaluation. + d::SelfChecks(output / "density_negative_regions_selfchecks.csv"); + const auto& fem = state.finiteElements; + if (!fem.mesh || !fem.densityFes || !fem.enthalpyFes || !fem.compactificationCoordinate || + 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 regions require the current scalar L2 density space on a conforming 3D mesh."); + int ranks = 0; + MPI_Comm_size(fem.mesh->GetComm(), &ranks); + if (ranks != 1) throw std::invalid_argument("Density regions currently require one MPI rank."); + for (int element = 0; element < fem.mesh->GetNE(); ++element) if (state.isStellar(element)) { + const auto* fe = fem.densityFes->GetFE(element); + if (!fe || fe->GetOrder() != 2 || fe->GetMapType() != mfem::FiniteElement::VALUE || + fe->GetGeomType() != mfem::Geometry::CUBE || fe->GetDof() != 27) + throw std::invalid_argument("Density root splitting requires scalar VALUE Q2 hexahedral density elements."); + } + const auto axes = d::ChooseAxes(state); + const mfem::Vector savedDensity(state.density), savedEnthalpy(state.enthalpy), savedPotential(state.potential), + savedGravity(state.gravityGradientReference), savedDisplacement(state.displacement); + mean_field::mapping::GridFunctionMappingEvaluator mapper(state.domainMapper, state.displacement, *fem.compactificationCoordinate); + mean_field::mapping::MappingPointContext mapped; + auto elements = d::Open(output / "density_negative_regions_elements.csv"); + elements << "tangential_order,interval_order,element,attribute,axis,orientation,surface_face,surface_face_count,"; + d::Header(elements); + auto summary = d::Open(output / "density_negative_regions_summary.csv"); + summary << "tangential_order,interval_order,"; + d::Header(summary); + auto intervals = d::Open(output / "density_negative_regions_intervals.csv"); + intervals << "tangential_order,element,attribute,axis,orientation,surface_face,u,v,tangential_weight,t_left,t_right," + "coefficient_c,coefficient_b,coefficient_a,negative_volume,negative_mass,minimum_density,witness_h," + "radius_minimum,radius_maximum,left_x,left_y,left_z,right_x,right_y,right_z\n"; + constexpr int intervalOrder = 14; + const auto& intervalRule = mfem::IntRules.Get(mfem::Geometry::SEGMENT, intervalOrder); + if (intervalRule.GetNPoints() != 8) throw std::logic_error("Unexpected density interval Gauss rule."); + for (int tangentialOrder : {14, 22, 30}) { + const auto& tangentRule = mfem::IntRules.Get(mfem::Geometry::SEGMENT, tangentialOrder); + if (tangentRule.GetNPoints() != tangentialOrder / 2 + 1) + throw std::logic_error("Unexpected density tangential Gauss rule."); + d::Measure global; + for (int element = 0; element < fem.mesh->GetNE(); ++element) { + if (!state.isStellar(element)) continue; + const auto& axis = axes[static_cast(element)]; + auto* transformation = fem.mesh->GetElementTransformation(element); + d::Measure local; + for (int iu = 0; iu < tangentRule.GetNPoints(); ++iu) for (int iv = 0; iv < tangentRule.GetNPoints(); ++iv) { + const auto& up = tangentRule.IntPoint(iu); + const auto& vp = tangentRule.IntPoint(iv); + const double u = up.x, v = vp.x; + const Real tangentWeight = static_cast(up.weight) * vp.weight; + const auto value = [&](Real t) { + const auto point = d::Point(axis, u, v, t); + const double rho = state.density.GetValue(element, point); + if (!std::isfinite(rho)) throw std::domain_error("Nonfinite density slice value."); + return rho; + }; + const auto polynomial = d::Interpolate(value(0), value(Real{0.5}), value(1)); + const Real tolerance = 512 * std::numeric_limits::epsilon() * polynomial.Scale(); + const auto verify = [&](Real t, double actual) { + const Real error = std::abs(polynomial.Value(t) - actual); + local.maximumReconstructionError = std::max(local.maximumReconstructionError, static_cast(error)); + local.maximumReconstructionTolerance = std::max(local.maximumReconstructionTolerance, static_cast(tolerance)); + if (!std::isfinite(error) || error > tolerance) + throw std::runtime_error("Frozen density slice fails independent Q2 reconstruction."); + }; + for (Real t : {Real{0.125}, Real{0.25}, Real{0.375}, Real{0.625}, Real{0.75}, Real{0.875}}) + verify(t, value(t)); + const auto parts = d::Partition(polynomial); + ++local.lines; + const auto evaluate = [&](Real t, d::Measure& measure, bool negativeRegion) { + const auto point = d::Point(axis, u, v, t); + if (mapper.EvaluatePoint(*transformation, point, mapped) != mean_field::mapping::MappingStatus::valid || mapped.compactified) + throw std::domain_error("Invalid stellar density-region mapping."); + const double rho = state.density.GetValue(element, point), h = state.enthalpy.GetValue(element, point); + if (!std::isfinite(rho) || !std::isfinite(h)) throw std::domain_error("Nonfinite density-region witness."); + verify(t, rho); + if (rho < measure.minimumRho) { + measure.minimumRho = rho; measure.witnessH = h; measure.witnessElement = element; + measure.witnessT = static_cast(t); measure.witnessU = u; measure.witnessV = v; + for (int component = 0; component < 3; ++component) measure.witnessPosition[component] = mapped.physical_position(component); + } + if (negativeRegion) { + const double radius = mapped.physical_position.Norml2(); + if (!std::isfinite(radius)) throw std::domain_error("Nonfinite density-region radius."); + measure.negativeRadiusMin = std::min(measure.negativeRadiusMin, radius); + measure.negativeRadiusMax = std::max(measure.negativeRadiusMax, radius); + } + transformation->SetIntPoint(&point); + const double referenceDeterminant = transformation->Jacobian().Det(); + const Real determinant = static_cast(referenceDeterminant) * mapped.mapping_determinant; + if (!(referenceDeterminant > 0) || !std::isfinite(referenceDeterminant) || + !(mapped.mapping_determinant > 0) || !std::isfinite(mapped.mapping_determinant) || + !(determinant > 0) || !std::isfinite(determinant)) + throw std::domain_error("Nonpositive or nonfinite density-region physical measure."); + return determinant; + }; + for (const auto& part : parts) { + d::Measure piece; + std::array leftPosition{}, rightPosition{}; + (void)evaluate(part.left, piece, part.negative); + for (int component = 0; component < 3; ++component) leftPosition[component] = mapped.physical_position(component); + (void)evaluate(part.right, piece, part.negative); + for (int component = 0; component < 3; ++component) rightPosition[component] = mapped.physical_position(component); + if (polynomial.a != 0) { + const Real stationary = -polynomial.b / (2 * polynomial.a); + if (stationary > part.left && stationary < part.right) (void)evaluate(stationary, piece, part.negative); + } + for (int it = 0; it < intervalRule.GetNPoints(); ++it) { + const auto& qp = intervalRule.IntPoint(it); + const Real t = part.left + (part.right - part.left) * qp.x; + const Real determinant = evaluate(t, piece, part.negative); + const Real weight = tangentWeight * (part.right - part.left) * qp.weight * determinant; + const Real rho = polynomial.Value(t); + if (!(weight > 0) || !std::isfinite(weight) || !std::isfinite(rho) || + (part.negative && rho > 0) || (!part.negative && rho < 0)) + throw std::runtime_error("Density region quadrature disagrees with its root partition."); + piece.volume += weight; + piece.mass += weight * rho; + ++piece.points; + } + if (part.negative) { + piece.negativeIntervals = 1; + piece.negativeVolume = piece.volume; + piece.negativeMass = -piece.mass; + piece.negativeTMin = part.left; + piece.negativeTMax = part.right; + if (tangentialOrder == 30) { + intervals << tangentialOrder << ',' << element << ',' << fem.mesh->GetAttribute(element) << ',' + << axis.coordinate << ',' << axis.orientation << ',' << axis.face << ',' << u << ',' << v << ',' + << tangentWeight << ',' << part.left << ',' << part.right << ',' << polynomial.c << ',' << polynomial.b << ',' + << polynomial.a << ',' << piece.negativeVolume << ',' << piece.negativeMass << ',' << piece.minimumRho << ',' + << piece.witnessH << ',' << piece.negativeRadiusMin << ',' << piece.negativeRadiusMax; + for (double x : leftPosition) intervals << ',' << x; + for (double x : rightPosition) intervals << ',' << x; + intervals << '\n'; + } + } + local.Add(piece); + } + } + if (!(local.volume > 0) || !std::isfinite(local.volume) || !std::isfinite(local.mass) || + !std::isfinite(local.negativeVolume) || !std::isfinite(local.negativeMass) || local.negativeMass < 0) + throw std::runtime_error("Invalid integrated density-region measure."); + elements << tangentialOrder << ',' << intervalOrder << ',' << element << ',' << fem.mesh->GetAttribute(element) << ',' + << axis.coordinate << ',' << axis.orientation << ',' << axis.face << ',' << axis.faceCount << ','; + d::Write(elements, local); + global.Add(local); + } + if (!(global.volume > 0)) throw std::invalid_argument("Density region diagnostic found no stellar volume."); + summary << tangentialOrder << ',' << intervalOrder << ','; + d::Write(summary, global); + summary.flush(); elements.flush(); intervals.flush(); + } + const bool unchanged = d::Identical(savedDensity, state.density) && d::Identical(savedEnthalpy, state.enthalpy) && + d::Identical(savedPotential, state.potential) && d::Identical(savedGravity, state.gravityGradientReference) && + d::Identical(savedDisplacement, state.displacement); + if (!unchanged) throw std::logic_error("Density region observation changed a source field."); + auto scope = d::Open(output / "density_negative_regions_scope.txt"); + scope << "source_fields_bitwise_identical=1\ntrial_preparations=0\nmpi_ranks=1\n" + << "n=" << n << "\nK=" << K << "\ndensity_basis=" << fem.densityFes->FEColl()->Name() << '\n' + << "Mass uses the reconstructed Q2 polynomial, verified against direct FE values; no density or EOS clipping.\n" + << "Volume weights contain one reference determinant and one deformation determinant.\n" + << "All stellar cells are included; surface axes point inward from the lowest-index stellar/nonstellar face.\n" + << "Other cells use reference axis0. Tangential orders14/22/30 use8/12/16 points per axis; interval order14 uses8.\n" + << "Whole-cell denominators integrate every root-partition interval with the same sliced quadrature.\n" + << "Negative intervals are exact for each reconstructed Q2 line up to floating arithmetic; tangential integration is sampled.\n" + << "Witnesses and radial bounds are samples, not global extrema or polynomial sign certificates.\n" + << "Near-double-root regions below reconstruction uncertainty are not certified; compare quadrature levels and errors.\n"; +} + +} // namespace experiment::plateau_diagnostics diff --git a/experiments/plateau_diagnostic.cpp b/experiments/plateau_diagnostic.cpp index cfd62a6..0a1a256 100644 --- a/experiments/plateau_diagnostic.cpp +++ b/experiments/plateau_diagnostic.cpp @@ -34,6 +34,7 @@ import mean_field; #include "reduced_surface_diagnostics.hpp" #include "spherical_harmonic_surface_diagnostics.hpp" #include "negative_density_diagnostics.hpp" +#include "negative_density_regions.hpp" using namespace mean_field; namespace { @@ -47,6 +48,7 @@ struct Options { bool weakShapeProbe = false, weakShapeFinal = false; bool reducedSurfaceProbe = false, reducedSurfaceAlgebraCheck = false; bool densityDiagnostic = false; + bool densityRegions = false; int surfaceDegree = 4, surfaceMaximumOrder = -1; std::string surfacePrescription = "nodal"; bool coupledTrialProbe = false, coupledSeedControl = false; @@ -129,6 +131,10 @@ Options Parse(int argc, char **argv) { o.reducedSurfaceAlgebraCheck = true; else if (a == "--density-diagnostic") o.densityDiagnostic = true; + else if (a == "--density-regions") { + o.densityDiagnostic = true; + o.densityRegions = true; + } else if (a == "--density-snapshot") o.densitySnapshot = v(); else if (a == "--surface-prescription") @@ -682,7 +688,7 @@ int RunContext(Context &context, const Options &o, const Scales &scales) { << "\nschur_correction_sweeps=" << o.schurCorrectionSweeps << "\nschur_relaxation=" << o.schurRelaxation << "\ncompiler=" << __VERSION__ << "\ncompiled=" << __DATE__ << ' ' << __TIME__ << '\n'; if (o.densityDiagnostic) { - metadata << "density_diagnostic=1\ndensity_source=" + metadata << "density_diagnostic=1\ndensity_region_integral=" << o.densityRegions << "\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) { @@ -721,7 +727,10 @@ int RunContext(Context &context, const Options &o, const Scales &scales) { 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 (o.densityRegions) + experiment::plateau_diagnostics::RecordDensityNegativeRegions(physical, o.n, scales.K, o.output); + else + 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.");