diff --git a/experiments/PRESSURE_OBSTRUCTION_METHODOLOGY_2026-09-17.md b/experiments/PRESSURE_OBSTRUCTION_METHODOLOGY_2026-09-17.md new file mode 100644 index 0000000..8acef6e --- /dev/null +++ b/experiments/PRESSURE_OBSTRUCTION_METHODOLOGY_2026-09-17.md @@ -0,0 +1,82 @@ +# Frozen pressure-obstruction diagnostic + +`plateau_pressure_obstruction_diagnostics.hpp` exports the experiment-only helper +`experiment::plateau_diagnostics::ProbePressureObstruction(state, fem, n, normalizedBase, normalizedDirection, output)`. +Include it after `import mean_field` and call it inside the existing context +diagnostics guard. The output directory must exist and its four +`pressure_obstruction_*.csv` files must not already exist. The helper is limited +to one MPI rank and the current registered polytropic field orders. No public +configuration, model equation, solver acceptance rule, or accepted state changes. + +The caller supplies an owned or borrowed frozen normalized state and correction; +the helper copies both before mutating workspaces. It restores the frozen state +before each candidate, and restores that state, inverse freshness, and the +caller's correction on return. It never calls `CommitPreparedTrial` or changes +accepted coefficients. The outer `WithState` guard remains responsible for +restoring the context's actual accepted state when a loaded frozen state differs +from it. Fatal restoration errors are propagated instead of silently claiming a +clean context. + +## Bound and location + +The existing `MeasureEnthalpyStep` scans the exact production pressure-force and +barotropic-closure quadrature rules. For scalar pullbacks, sampled enthalpy along +a fixed direction obeys `h_q(alpha) = h_q(0) + alpha * dh_q`. Each family reports +its smallest `-h/dh` for negative `dh`, and their minimum is the union bound. +Zero means that a sampled zero-enthalpy point is decreasing; infinity means no +sampled upper bound. Neither means that all other preparation checks pass. + +`pressure_obstruction_summary.csv` retains the complete existing scan, including +point counts, negative samples, both limiting points and both bounds. +`pressure_obstruction_limits.csv` adds the limiting element's Cartesian reference +and frozen physical positions. Element-local integration coordinates are recorded +separately. The exact topological flag `element_touches_stellar_surface` locates +the element relative to the tagged stellar/vacuum interface. It is not a physical +distance or a certificate that the failure occurs at the boundary. Radii are +measured about the coordinate origin and are not normalized by a critical radius. + +The scanner uses a pointwise polynomial contraction; optimized production may +use tensor contractions. These sample the same polynomial but need not agree +bitwise. Trials avoid alpha equal to the bound. Full-direction factors are +0.5, 0.9, 0.99, 1.01, 1.1, and 2 times the union bound, plus distinct fixed +steps 1, 1/2, 1/4, 1/8, and 1/16. Fixed steps coincident with the boundary are +omitted. The helper also tests enthalpy-only and shape-only components at the +same factors 0.9, 0.99, 1.01, and 1.1, retaining their original amplitudes. +This is at most 19 candidate preparations plus frozen-state restorations. For a +zero or infinite bound, a smaller fixed grid replaces the relative-bound grid. + +## Exact rejection and interpretation + +`pressure_obstruction_trials.csv` records the full root preparation reason and +stage enums, an explicit reason name, the EOS error enum, and whether that EOS +code is applicable. Nonthermodynamic failures can carry a default EOS code; +do not interpret it without the applicability flag. Missing determinant data +remain NaN. Successful preparation records normalized and physical residual +norms, with normalized residual blocks in `pressure_obstruction_blocks.csv`. +The physical whole-vector norm combines units and is for reproduction, not a +physical acceptance metric. The helper does not run the line-search merit test. + +For this production pressure operator, `thermodynamic_domain` at the pressure +stage identifies a finite sampled negative enthalpy; `inverted_geometry` +identifies the mapping's nonpositive-determinant/weight branch. Other mapping +statuses and nonfinite input/results have separate reasons. Thus the exact +failure payload resolves the ambiguity of the earlier observer's generic +`inadmissible_state` plus `pressure force` label. The recorded reason identifies +the selected preparation failure, not necessarily every defect in the candidate. + +Both predicted and independently re-evaluated enthalpy at the two frozen limiting +points are saved for every candidate. This is not a rescan for the most negative +trial point: the most negative point need not be the earliest-crossing point. +If full and enthalpy-only trials cross the same pressure thermodynamic boundary +while shape-only remains admissible, that isolates the enthalpy obstruction from +shape geometry along this correction. It does not establish the spatial cause +of the Newton direction or suggest permitting negative enthalpy. + +If the caller reconstructs the correction from a recorded trial as +`(candidate - base) / alpha`, cancellation at tiny alpha can perturb the direction +and bound. Preserve that provenance and replay the original candidate directly +as an independent control. This helper diagnoses only its supplied direction. + +Implementation review only: no compilation or numerical runs have been performed +for this helper at the time of this note. Root integration and validation remain +pending. diff --git a/experiments/plateau_pressure_obstruction_diagnostics.hpp b/experiments/plateau_pressure_obstruction_diagnostics.hpp new file mode 100644 index 0000000..5c1197e --- /dev/null +++ b/experiments/plateau_pressure_obstruction_diagnostics.hpp @@ -0,0 +1,280 @@ +#pragma once + +// Diagnostic only. Include after import mean_field; invoke inside WithState. +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "plateau_enthalpy_step_diagnostics.hpp" + +namespace experiment::plateau_diagnostics { +namespace pressure_obstruction_detail { +inline std::ofstream Open(const std::filesystem::path &path) { + if (std::filesystem::exists(path)) + throw std::runtime_error("Refusing to overwrite pressure-obstruction output: " + path.string()); + std::ofstream stream(path); + stream.exceptions(std::ios::failbit | std::ios::badbit); + stream << std::setprecision(17); + return stream; +} + +inline const char *ReasonName(mean_field::operators::StellarEquilibriumPreparationRejectionReason reason) { + using R = mean_field::operators::StellarEquilibriumPreparationRejectionReason; + switch (reason) { + case R::inverted_geometry: return "inverted_geometry"; + case R::non_finite_geometry: return "non_finite_geometry"; + case R::thermodynamic_domain: return "thermodynamic_domain"; + case R::non_finite_thermodynamics: return "non_finite_thermodynamics"; + case R::inadmissible_physics: return "inadmissible_physics"; + case R::non_finite_physics: return "non_finite_physics"; + } + return "unknown"; +} + +inline bool Finite(const mfem::Vector &v) { + for (int i = 0; i < v.Size(); ++i) + if (!std::isfinite(v(i))) return false; + return true; +} + +struct Limit { + std::string family; + int element{-1}, point{-1}, attribute{-1}; + mfem::IntegrationPoint ip; + double alpha{std::numeric_limits::infinity()}; + double h{std::numeric_limits::quiet_NaN()}; + double dh{std::numeric_limits::quiet_NaN()}; +}; + +inline void AddAlpha(std::vector &values, double alpha, double boundary) { + if (!std::isfinite(alpha) || alpha <= 0.0) return; + // The boundary itself is deliberately excluded: algebraically equivalent + // state addition / polynomial contractions can differ at roundoff there. + if (std::isfinite(boundary) && boundary > 0.0 && + std::abs(alpha / boundary - 1.0) < 1.e-10) return; + for (double previous : values) + if (std::abs(alpha - previous) <= 1.e-12 * std::max(alpha, previous)) return; + values.push_back(alpha); +} +} // namespace pressure_obstruction_detail + +template +void ProbePressureObstruction(State &s, const mean_field::fem::FEM &f, double n, + const mfem::Vector &normalizedBase, const mfem::Vector &normalizedDirection, + const std::filesystem::path &output) { + using namespace mean_field; + namespace detail = pressure_obstruction_detail; + using DomainSchema = utils::domain::CoreEnvelopeVacuumDomainSchema; + if (!f.mesh || !f.enthalpyFes || !f.displacementFes || !f.compactificationCoordinate) + throw std::invalid_argument("Pressure obstruction requires complete field/mapping data."); + int ranks = 0; + if (MPI_Comm_size(f.mesh->GetComm(), &ranks) != MPI_SUCCESS || ranks != 1) + throw std::invalid_argument("Pressure obstruction requires exactly one MPI rank."); + if (normalizedBase.Size() != s.Problem().StateSize() || + normalizedDirection.Size() != normalizedBase.Size() || + !detail::Finite(normalizedBase) || !detail::Finite(normalizedDirection)) + throw std::invalid_argument("Pressure obstruction requires finite compatible base/direction vectors."); + if (!std::filesystem::is_directory(output)) + throw std::invalid_argument("Pressure obstruction requires an existing output directory."); + + // Own these snapshots: callers may pass one of the live workspaces. + const mfem::Vector frozen(normalizedBase), direction(normalizedDirection), savedCorrection(s.normalizedCorrection); + auto restore = [&]() { + s.trialNormalizedState = frozen; + if (!s.PrepareTrial()) throw std::runtime_error("Pressure-obstruction frozen base is inadmissible."); + }; + auto finish = [&]() { + restore(); + s.physicalInverse->Refresh(); + s.normalizedCorrection = savedCorrection; + }; + + try { + restore(); + s.normalizedCorrection = direction; + mfem::Vector physicalBase, physicalDirection; + s.normalizedOperator->DenormalizeState(frozen, physicalBase); + s.normalizedOperator->DenormalizeState(direction, physicalDirection); + const auto measures = MeasureEnthalpyStep(s, f, n, &physicalBase); + const double boundary = measures.at("alpha_thermodynamic"); + if (measures.at("thermodynamic_base_admissible") != 1.0 || std::isnan(boundary)) + throw std::runtime_error("Production-rule scan found an invalid frozen enthalpy base."); + + auto summary = detail::Open(output / "pressure_obstruction_summary.csv"); + summary << "metric,value\n"; + for (const auto &[key, value] : measures) summary << key << ',' << value << '\n'; + summary << "base_normalized_residual," << s.trialNormalizedResidual.Norml2() << '\n'; + summary << "normalized_direction_norm," << direction.Norml2() << '\n'; + + mfem::Vector enthalpyDirection(direction.Size()), shapeDirection(direction.Size()); + enthalpyDirection = 0.0; + shapeDirection = 0.0; + int hOffset = -1, hSize = -1, shapeSize = -1; + for (const auto &block : s.Problem().GetManifest().valueBlocks()) { + if (block.stableId == "specific_enthalpy") { + hOffset = block.offset; + hSize = block.size; + for (int i = 0; i < block.size; ++i) enthalpyDirection(block.offset + i) = direction(block.offset + i); + } + if (block.stableId == "surface_deformation") { + shapeSize = block.size; + for (int i = 0; i < block.size; ++i) shapeDirection(block.offset + i) = direction(block.offset + i); + } + } + if (hOffset < 0 || hSize <= 0 || shapeSize <= 0) + throw std::runtime_error("Pressure obstruction could not locate enthalpy and shape blocks."); + summary << "normalized_enthalpy_direction_norm," << enthalpyDirection.Norml2() << '\n'; + summary << "normalized_shape_direction_norm," << shapeDirection.Norml2() << '\n'; + + std::array limits; + for (int i = 0; i < 2; ++i) { + auto &limit = limits[i]; + const std::string prefix = i == 0 ? "" : "pressure_"; + limit.family = i == 0 ? "closure" : "pressure"; + limit.alpha = measures.at(prefix + "alpha_h"); + limit.element = static_cast(measures.at(prefix + "limiting_element")); + limit.point = static_cast(measures.at(prefix + "limiting_point")); + limit.attribute = static_cast(measures.at(prefix + "limiting_attribute")); + limit.h = measures.at(prefix + "limiting_h"); + limit.dh = measures.at(prefix + "limiting_dh"); + limit.ip.Set3(measures.at(prefix + "limiting_ip_x"), measures.at(prefix + "limiting_ip_y"), + measures.at(prefix + "limiting_ip_z")); + } + + // Surface is a material interface, not necessarily the outer mesh boundary. + std::vector touchesSurface(static_cast(f.mesh->GetNE()), false); + const int surfaceAttribute = DomainSchema::template boundary_attribute(); + for (int b = 0; b < f.mesh->GetNBE(); ++b) { + if (f.mesh->GetBdrAttribute(b) != surfaceAttribute) continue; + const int face = f.mesh->GetBdrElementFaceIndex(b); + if (face < 0) throw std::runtime_error("Missing stellar-surface mesh face."); + auto *faces = f.mesh->GetFaceElementTransformations(face); + if (!faces) throw std::runtime_error("Missing stellar-surface face transformation."); + for (int element : {faces->Elem1No, faces->Elem2No}) + if (element >= 0 && element < f.mesh->GetNE() && + DomainSchema::template attribute_belongs_to(f.mesh->GetAttribute(element))) + touchesSurface[static_cast(element)] = true; + } + + mfem::ParGridFunction displacement(f.displacementFes.get()), trialH(f.enthalpyFes.get()); + displacement.SetFromTrueDofs(s.Problem().GetPhysicalOperator().GetGeneratedVolumeDisplacement()); + mapping::GridFunctionMappingEvaluator mapper(s.Problem().GetDiscretization().domainMapper(), displacement, + *f.compactificationCoordinate); + auto location = detail::Open(output / "pressure_obstruction_limits.csv"); + location << "family,element,point,attribute,alpha_h,h,dh,ip_x,ip_y,ip_z,element_touches_stellar_surface," + "mapping_status,reference_x,reference_y,reference_z,reference_radius,physical_x,physical_y," + "physical_z,physical_radius\n"; + const double nan = std::numeric_limits::quiet_NaN(); + for (const auto &limit : limits) { + location << limit.family << ',' << limit.element << ',' << limit.point << ',' << limit.attribute << ',' + << limit.alpha << ',' << limit.h << ',' << limit.dh << ',' << limit.ip.x << ',' << limit.ip.y + << ',' << limit.ip.z << ','; + if (limit.element < 0) { + location << "-1,-1,nan,nan,nan,nan,nan,nan,nan,nan\n"; + continue; + } + mapping::MappingPointContext point; + auto *transformation = f.mesh->GetElementTransformation(limit.element); + const auto status = mapper.EvaluatePoint(*transformation, limit.ip, point); + location << touchesSurface[static_cast(limit.element)] << ',' << static_cast(status); + if (status == mapping::MappingStatus::valid) { + for (int d = 0; d < 3; ++d) location << ',' << point.reference_position(d); + location << ',' << point.reference_position.Norml2(); + for (int d = 0; d < 3; ++d) location << ',' << point.physical_position(d); + location << ',' << point.physical_position.Norml2() << '\n'; + } else location << ",nan,nan,nan,nan,nan,nan,nan,nan\n"; + } + summary.flush(); + location.flush(); + + auto trials = detail::Open(output / "pressure_obstruction_trials.csv"); + trials << "direction,alpha,alpha_over_union_bound,admissible,reason,reason_name,stage,eos_code," + "eos_code_applicable,min_jacobian,normalized_residual_norm,physical_residual_norm," + "closure_limiting_h_prediction,pressure_limiting_h_prediction," + "closure_limiting_h_evaluated,pressure_limiting_h_evaluated\n"; + auto blocks = detail::Open(output / "pressure_obstruction_blocks.csv"); + blocks << "direction,alpha,block,normalized_residual_norm\n"; + const auto adapter = field::make_field_dof_grid_function_adapter(*f.enthalpyFes); + mfem::Vector candidatePhysical, reducedH(hSize); + std::vector fullAlphas, splitAlphas; + if (std::isfinite(boundary) && boundary > 0.0) { + for (double factor : {.5, .9, .99, 1.01, 1.1, 2.0}) detail::AddAlpha(fullAlphas, factor * boundary, boundary); + for (double factor : {.9, .99, 1.01, 1.1}) detail::AddAlpha(splitAlphas, factor * boundary, boundary); + } else { + // Zero: every positive step may violate h>=0. Infinity: no sampled bound. + for (double a : (boundary == 0.0 ? std::array{1.e-6, 1.e-4} : std::array{.5, 1.0})) + detail::AddAlpha(splitAlphas, a, boundary); + for (double a : splitAlphas) detail::AddAlpha(fullAlphas, a, boundary); + } + for (double a : {1.0, .5, .25, .125, .0625}) detail::AddAlpha(fullAlphas, a, boundary); + summary << "candidate_trial_count," << fullAlphas.size() + 2 * splitAlphas.size() << '\n'; + summary.flush(); + + auto probe = [&](const char *name, const mfem::Vector &d, const std::vector &alphas, bool changesH) { + for (double alpha : alphas) { + restore(); + s.normalizedCorrection = d; + s.trialNormalizedState = frozen; + s.trialNormalizedState.Add(alpha, d); + s.normalizedOperator->DenormalizeState(s.trialNormalizedState, candidatePhysical); + if (!detail::Finite(candidatePhysical) || !detail::Finite(s.trialNormalizedState)) + throw std::runtime_error("Pressure-obstruction trial coefficients overflowed; no trial was committed."); + for (int i = 0; i < hSize; ++i) reducedH(i) = candidatePhysical(hOffset + i); + adapter.scatter(reducedH, trialH); + std::array evaluated{nan, nan}; + for (int i = 0; i < 2; ++i) + if (limits[i].element >= 0) evaluated[i] = trialH.GetValue(limits[i].element, limits[i].ip); + const auto result = s.PrepareTrial(); + trials << name << ',' << alpha << ',' + << (std::isfinite(boundary) && boundary > 0.0 ? alpha / boundary : nan) << ',' + << result.has_value() << ','; + if (result) { + trials << "-1,none,-1,-1,0," << result->minimumJacobianDeterminant << ',' + << s.trialNormalizedResidual.Norml2() << ',' + << s.normalizedOperator->GetPhysicalResidual().Norml2(); + for (const auto &row : s.Problem().GetManifest().residualBlocks()) { + long double squared = 0.0; + for (int i = 0; i < row.size; ++i) { + const long double value = s.trialNormalizedResidual(row.offset + i); + squared += value * value; + } + blocks << name << ',' << alpha << ',' << row.stableId << ',' << std::sqrt(squared) << '\n'; + } + } else { + const auto &error = result.error(); + const bool thermodynamic = error.reason == operators::StellarEquilibriumPreparationRejectionReason::thermodynamic_domain || + error.reason == operators::StellarEquilibriumPreparationRejectionReason::non_finite_thermodynamics; + trials << static_cast(error.reason) << ',' << detail::ReasonName(error.reason) << ',' + << static_cast(error.stage) << ',' << static_cast(error.thermodynamicErrorCode) + << ',' << thermodynamic << ',' << error.minimumJacobianDeterminant << ",nan,nan"; + } + for (const auto &limit : limits) trials << ',' << (limit.h + (changesH ? alpha * limit.dh : 0.0)); + for (double value : evaluated) trials << ',' << value; + trials << '\n'; + trials.flush(); + blocks.flush(); + } + }; + probe("full", direction, fullAlphas, true); + probe("enthalpy_only", enthalpyDirection, splitAlphas, true); + probe("shape_only", shapeDirection, splitAlphas, false); + finish(); + } catch (...) { + const auto error = std::current_exception(); + try { finish(); } + catch (...) { + s.normalizedCorrection = savedCorrection; + std::throw_with_nested(std::runtime_error("Pressure-obstruction diagnostics could not restore the frozen base.")); + } + std::rethrow_exception(error); + } +} +} // namespace experiment::plateau_diagnostics