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

experiment: localize pressure preparation admissibility boundaries

Browse Source
This commit is contained in:
Emily Boudreaux
2026-09-17 08:16:37 -04:00
parent fd08f7b0e9
commit 6753cc5ddc
2 changed files with 362 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

82
experiments/PRESSURE_OBSTRUCTION_METHODOLOGY_2026-09-17.md Normal file
View File

@@ -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.

280
experiments/plateau_pressure_obstruction_diagnostics.hpp Normal file
View File

@@ -0,0 +1,280 @@
#pragma once
// Diagnostic only. Include after import mean_field; invoke inside WithState.
#include <algorithm>
#include <array>
#include <cmath>
#include <exception>
#include <filesystem>
#include <fstream>
#include <iomanip>
#include <limits>
#include <stdexcept>
#include <string>
#include <vector>
#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<double>::infinity()};
double h{std::numeric_limits<double>::quiet_NaN()};
double dh{std::numeric_limits<double>::quiet_NaN()};
};
inline void AddAlpha(std::vector<double> &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 <class State>
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<detail::Limit, 2> 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<int>(measures.at(prefix + "limiting_element"));
limit.point = static_cast<int>(measures.at(prefix + "limiting_point"));
limit.attribute = static_cast<int>(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<bool> touchesSurface(static_cast<std::size_t>(f.mesh->GetNE()), false);
const int surfaceAttribute = DomainSchema::template boundary_attribute<utils::domain::StellarSurface>();
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<utils::domain::Stellar>(f.mesh->GetAttribute(element)))
touchesSurface[static_cast<std::size_t>(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<double>::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<std::size_t>(limit.element)] << ',' << static_cast<int>(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<field::Enthalpy, DomainSchema>(*f.enthalpyFes);
mfem::Vector candidatePhysical, reducedH(hSize);
std::vector<double> 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<double> &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<double, 2> 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<int>(error.reason) << ',' << detail::ReasonName(error.reason) << ','
<< static_cast<int>(error.stage) << ',' << static_cast<int>(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