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

Add preparation-free surface carrier trace diagnostics

Browse Source
This commit is contained in:
Emily Boudreaux
2026-09-17 08:34:49 -04:00
parent 9e1b65dbde
commit b1dc2acf35
1 changed files with 123 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

123
experiments/plateau_pressure_obstruction_diagnostics.hpp
View File

@@ -277,4 +277,127 @@ void ProbePressureObstruction(State &s, const mean_field::fem::FEM &f, double n,
std::rethrow_exception(error);
}
}
// Supplemental, preparation-free inspection of the ACTUAL algebraic surface
// carrier rows. Boundary-only interpolation is a contribution to the interior
// polynomial, not its trace evaluated on the physical stellar surface.
template <class State>
void ProbeSurfaceTrace(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;
using EnthalpyField = field::Field<field::Enthalpy>;
if (!f.mesh || !f.enthalpyFes || !f.quadratureFactory)
throw std::invalid_argument("Surface trace inspection requires an enthalpy FEM.");
if (normalizedBase.Size() != s.Problem().StateSize() || normalizedDirection.Size() != normalizedBase.Size() ||
!detail::Finite(normalizedBase) || !detail::Finite(normalizedDirection))
throw std::invalid_argument("Surface trace inspection requires finite compatible vectors.");
if (!std::filesystem::is_directory(output))
throw std::invalid_argument("Surface trace inspection requires an existing output directory.");
const mfem::Vector savedCorrection(s.normalizedCorrection), direction(normalizedDirection);
mfem::Vector physicalBase, physicalDirection;
s.normalizedOperator->DenormalizeState(normalizedBase, physicalBase);
s.normalizedOperator->DenormalizeState(direction, physicalDirection);
// No PrepareTrial, Refresh, BuildResidual, or accepted-state mutation here.
try {
s.normalizedCorrection = direction;
const auto measures = MeasureEnthalpyStep(s, f, n, &physicalBase);
s.normalizedCorrection = savedCorrection;
int offset = -1, size = -1;
for (const auto &block : s.Problem().GetManifest().valueBlocks())
if (block.stableId == "specific_enthalpy") { offset = block.offset; size = block.size; }
const auto &surface = s.Problem().GetPressureSurfaceRows();
if (offset < 0 || size <= 0 || surface.field_size() != size || surface.size() <= 0)
throw std::runtime_error("Surface trace inspection found an incompatible or empty carrier map.");
mfem::Vector boundaryH(size), boundaryDh(size);
boundaryH = 0.0;
boundaryDh = 0.0;
auto rows = detail::Open(output / "surface_trace_rows.csv");
rows << "surface_index,reduced_h_index,root_index,base_h,direction_dh,h_plus_dh\n";
std::array<double, 3> minimum, maximum, maximumAbs{};
minimum.fill(std::numeric_limits<double>::infinity());
maximum.fill(-std::numeric_limits<double>::infinity());
std::array<long double, 3> squared{};
std::array<std::uint64_t, 3> negative{}, nonzero{};
double alphaSurface = std::numeric_limits<double>::infinity();
bool baseNonnegative = true;
for (int i = 0; i < surface.size(); ++i) {
const int reduced = surface.reduced_dofs()[i];
if (reduced < 0 || reduced >= size)
throw std::runtime_error("Surface carrier row is outside the enthalpy block.");
const double h = physicalBase(offset + reduced), dh = physicalDirection(offset + reduced);
if (!std::isfinite(h) || !std::isfinite(dh) || !std::isfinite(h + dh))
throw std::runtime_error("Surface trace coefficients are nonfinite.");
boundaryH(reduced) = h;
boundaryDh(reduced) = dh;
rows << i << ',' << reduced << ',' << offset + reduced << ',' << h << ',' << dh << ',' << h + dh << '\n';
const std::array values{h, dh, h + dh};
for (int j = 0; j < 3; ++j) {
minimum[j] = std::min(minimum[j], values[j]);
maximum[j] = std::max(maximum[j], values[j]);
maximumAbs[j] = std::max(maximumAbs[j], std::abs(values[j]));
squared[j] += static_cast<long double>(values[j]) * values[j];
negative[j] += values[j] < 0.0;
nonzero[j] += values[j] != 0.0;
}
baseNonnegative = baseNonnegative && h >= 0.0;
if (dh < 0.0) alphaSurface = std::min(alphaSurface, -h / dh);
}
auto summary = detail::Open(output / "surface_trace_summary.csv");
summary << "metric,value\nphysical_state_size," << physicalBase.Size() << "\nh_root_offset," << offset
<< "\nh_block_size," << size << "\nsurface_row_count," << surface.size() << '\n';
const std::array<const char *, 3> names{"base_h", "direction_dh", "h_plus_dh"};
for (int j = 0; j < 3; ++j) {
summary << names[j] << "_minimum," << minimum[j] << '\n'
<< names[j] << "_maximum," << maximum[j] << '\n'
<< names[j] << "_maximum_absolute," << maximumAbs[j] << '\n'
<< names[j] << "_norm," << std::sqrt(squared[j]) << '\n'
<< names[j] << "_negative_count," << negative[j] << '\n'
<< names[j] << "_nonzero_count," << nonzero[j] << '\n';
}
summary << "surface_coefficient_base_nonnegative," << baseNonnegative << '\n'
<< "surface_coefficient_alpha_h," << (baseNonnegative ? alphaSurface : std::numeric_limits<double>::quiet_NaN()) << '\n'
<< "closure_alpha_h," << measures.at("closure_alpha_h") << '\n'
<< "pressure_alpha_h," << measures.at("pressure_alpha_h") << '\n';
const auto adapter = field::make_field_dof_grid_function_adapter<field::Enthalpy, DomainSchema>(*f.enthalpyFes);
mfem::ParGridFunction hBoundary(f.enthalpyFes.get()), dhBoundary(f.enthalpyFes.get());
adapter.scatter(boundaryH, hBoundary);
adapter.scatter(boundaryDh, dhBoundary);
auto limits = detail::Open(output / "surface_trace_limits.csv");
limits << "family,element,point,attribute,geometry_type,transformation_order_w,eos_extra_order,rule_points,"
"ip_x,ip_y,ip_z,alpha_h,full_h,full_dh,boundary_h,boundary_dh,interior_h,interior_dh,"
"full_h_plus_dh,boundary_h_plus_dh,interior_h_plus_dh\n";
for (bool pressure : {false, true}) {
const std::string prefix = pressure ? "pressure_" : "";
const int element = static_cast<int>(measures.at(prefix + "limiting_element"));
const int point = static_cast<int>(measures.at(prefix + "limiting_point"));
if (element < 0) continue; // No decreasing sample: the summary retains the infinite bound.
auto *transformation = f.mesh->GetElementTransformation(element);
const int extraOrder = static_cast<int>(measures.at(prefix + "eos_extra_order"));
const auto query = pressure ? EnthalpyField::template make_query<field::Enthalpy::Form::PressureForce>(
quadrature::QuadratureRole::discretization, transformation->OrderW(), std::array<int, 1>{extraOrder},
utils::DOMAINS::STELLAR, quadrature::MappingKind::general)
: EnthalpyField::template make_query<field::Enthalpy::Form::EosClosureSource>(
quadrature::QuadratureRole::discretization, transformation->OrderW(), std::array<int, 1>{extraOrder},
utils::DOMAINS::STELLAR, quadrature::MappingKind::general);
const auto rule = f.quadratureFactory->get(query, transformation->GetGeometryType());
if (!rule.integration_rule || point < 0 || point >= rule.integration_rule->GetNPoints())
throw std::runtime_error("Surface trace inspection could not resolve its production limiting point.");
const auto &ip = rule.integration_rule->IntPoint(point);
const double h = measures.at(prefix + "limiting_h"), dh = measures.at(prefix + "limiting_dh");
const double bh = hBoundary.GetValue(element, ip), bdh = dhBoundary.GetValue(element, ip);
limits << (pressure ? "pressure" : "closure") << ',' << element << ',' << point << ','
<< f.mesh->GetAttribute(element) << ',' << static_cast<int>(transformation->GetGeometryType()) << ','
<< transformation->OrderW() << ',' << extraOrder << ',' << rule.integration_rule->GetNPoints() << ','
<< ip.x << ',' << ip.y << ',' << ip.z << ',' << measures.at(prefix + "alpha_h") << ','
<< h << ',' << dh << ',' << bh << ',' << bdh << ',' << h - bh << ',' << dh - bdh << ','
<< h + dh << ',' << bh + bdh << ',' << (h - bh) + (dh - bdh) << '\n';
}
} catch (...) {
s.normalizedCorrection = savedCorrection;
throw;
}
}
} // namespace experiment::plateau_diagnostics