440 lines
27 KiB
C++
440 lines
27 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);
|
|
}
|
|
}
|
|
// 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);
|
|
mfem::Vector fullH(size), fullDh(size);
|
|
for (int i = 0; i < size; ++i) {
|
|
fullH(i) = physicalBase(offset + i);
|
|
fullDh(i) = physicalDirection(offset + i);
|
|
}
|
|
mfem::ParGridFunction hFull(f.enthalpyFes.get()), dhFull(f.enthalpyFes.get());
|
|
adapter.scatter(fullH, hFull);
|
|
adapter.scatter(fullDh, dhFull);
|
|
auto nodes = detail::Open(output / "surface_trace_nodes.csv");
|
|
nodes << "family,element,attribute,geometry_type,fe_order,local_node,node_count,ip_x,ip_y,ip_z,"
|
|
"h,dh,alpha_union,h_plus_alpha_union_dh,h_plus_dh\n";
|
|
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 auto *finiteElement = f.enthalpyFes->GetFE(element);
|
|
const auto &elementNodes = finiteElement->GetNodes();
|
|
const double alphaUnion = measures.at("alpha_thermodynamic");
|
|
// Evaluate the scattered physical fields at FE nodes; do not infer
|
|
// nodal values by indexing potentially signed/oriented local DOFs.
|
|
// A shared limiting element is exported once per family for identity.
|
|
for (int node = 0; node < elementNodes.GetNPoints(); ++node) {
|
|
const auto &ip = elementNodes.IntPoint(node);
|
|
const double h = hFull.GetValue(element, ip), dh = dhFull.GetValue(element, ip);
|
|
nodes << (pressure ? "pressure" : "closure") << ',' << element << ','
|
|
<< f.mesh->GetAttribute(element) << ',' << static_cast<int>(transformation->GetGeometryType()) << ','
|
|
<< finiteElement->GetOrder() << ',' << node << ',' << elementNodes.GetNPoints() << ','
|
|
<< ip.x << ',' << ip.y << ',' << ip.z << ',' << h << ',' << dh << ',' << alphaUnion << ','
|
|
<< (std::isfinite(alphaUnion) ? h + alphaUnion * dh : std::numeric_limits<double>::quiet_NaN()) << ','
|
|
<< h + dh << '\n';
|
|
}
|
|
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);
|
|
auto rule = f.quadratureFactory->get(query, transformation->GetGeometryType());
|
|
if (!pressure) {
|
|
const auto massQuery = field::Field<field::Density>::template make_query<
|
|
field::Density::Form::EosClosureMass>(
|
|
quadrature::QuadratureRole::discretization, transformation->OrderW(), {},
|
|
utils::DOMAINS::STELLAR, quadrature::MappingKind::general);
|
|
const auto massResolution = f.quadratureFactory->get(massQuery, transformation->GetGeometryType());
|
|
if (massResolution.resolution.order > rule.resolution.order)
|
|
rule = massResolution;
|
|
}
|
|
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
|