From b1dc2acf35c4a2b339f61a3dd7d6ebc84318d22f Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Thu, 17 Sep 2026 08:34:49 -0400 Subject: [PATCH] Add preparation-free surface carrier trace diagnostics --- ...ateau_pressure_obstruction_diagnostics.hpp | 123 ++++++++++++++++++ 1 file changed, 123 insertions(+) diff --git a/experiments/plateau_pressure_obstruction_diagnostics.hpp b/experiments/plateau_pressure_obstruction_diagnostics.hpp index 5c1197e..7e19854 100644 --- a/experiments/plateau_pressure_obstruction_diagnostics.hpp +++ b/experiments/plateau_pressure_obstruction_diagnostics.hpp @@ -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 +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; + 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 minimum, maximum, maximumAbs{}; + minimum.fill(std::numeric_limits::infinity()); + maximum.fill(-std::numeric_limits::infinity()); + std::array squared{}; + std::array negative{}, nonzero{}; + double alphaSurface = std::numeric_limits::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(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 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::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(*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(measures.at(prefix + "limiting_element")); + const int point = static_cast(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(measures.at(prefix + "eos_extra_order")); + const auto query = pressure ? EnthalpyField::template make_query( + quadrature::QuadratureRole::discretization, transformation->OrderW(), std::array{extraOrder}, + utils::DOMAINS::STELLAR, quadrature::MappingKind::general) + : EnthalpyField::template make_query( + quadrature::QuadratureRole::discretization, transformation->OrderW(), std::array{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(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