281 lines
16 KiB
C++
281 lines
16 KiB
C++
#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
|