421 lines
25 KiB
C++
421 lines
25 KiB
C++
#pragma once
|
|
|
|
// Include after import mean_field. Root splitting removes the discontinuous
|
|
// sign indicator from the integration coordinate; tangential and mapped-volume
|
|
// integration still require the quadrature convergence study recorded here.
|
|
#include "polytrope_physical_state.hpp"
|
|
|
|
#include <algorithm>
|
|
#include <array>
|
|
#include <cmath>
|
|
#include <cstring>
|
|
#include <filesystem>
|
|
#include <fstream>
|
|
#include <iomanip>
|
|
#include <limits>
|
|
#include <string>
|
|
#include <vector>
|
|
|
|
namespace experiment::plateau_diagnostics {
|
|
namespace negative_density_regions_detail {
|
|
|
|
using Real = long double;
|
|
|
|
struct Polynomial final {
|
|
Real c, b, a;
|
|
[[nodiscard]] Real Value(Real t) const { return std::fma(std::fma(a, t, b), t, c); }
|
|
[[nodiscard]] Real Integral(Real left, Real right) const {
|
|
const Real width = right - left;
|
|
// Subtracting antiderivatives loses a thin pocket's entire mass when
|
|
// both antiderivatives are large. This centered quadratic identity does
|
|
// not subtract nearly equal endpoint values.
|
|
return width * (Value(left + width / 2) + a * width * width / 12);
|
|
}
|
|
[[nodiscard]] Real Scale() const { return std::max({std::abs(a), std::abs(b), std::abs(c)}); }
|
|
};
|
|
|
|
struct Interval final { Real left, right; bool negative; };
|
|
|
|
inline Polynomial Interpolate(Real left, Real middle, Real right) {
|
|
return {left, 4 * middle - 3 * left - right, 2 * (left + right - 2 * middle)};
|
|
}
|
|
|
|
inline std::vector<Interval> Partition(const Polynomial& p) {
|
|
if (!std::isfinite(p.a) || !std::isfinite(p.b) || !std::isfinite(p.c))
|
|
throw std::domain_error("Nonfinite density slice polynomial.");
|
|
std::vector<Real> points{0, 1};
|
|
const auto add = [&](Real root) {
|
|
if (std::isfinite(root) && root > 0 && root < 1) points.push_back(root);
|
|
};
|
|
const Real scale = p.Scale();
|
|
if (scale > 0) {
|
|
const Real a = p.a / scale, b = p.b / scale, c = p.c / scale;
|
|
// No relative cutoff turns a small quadratic into a line: a thin
|
|
// negative interval can be the very effect this diagnostic must retain.
|
|
if (a == 0) {
|
|
if (b != 0) add(-c / b);
|
|
} else {
|
|
const Real discriminant = std::fma(-4 * a, c, b * b);
|
|
if (discriminant == 0) {
|
|
add(-b / (2 * a));
|
|
} else if (discriminant > 0) {
|
|
const Real q = -Real{0.5} * (b + std::copysign(std::sqrt(discriminant), b));
|
|
// The product-of-roots form avoids cancellation at a root near
|
|
// the face. q cannot vanish for a positive discriminant.
|
|
add(q / a);
|
|
add(c / q);
|
|
}
|
|
}
|
|
}
|
|
std::sort(points.begin(), points.end());
|
|
points.erase(std::unique(points.begin(), points.end()), points.end());
|
|
std::vector<Interval> result;
|
|
for (std::size_t i = 1; i < points.size(); ++i) {
|
|
const Real left = points[i - 1], right = points[i];
|
|
if (!(right > left)) throw std::logic_error("Invalid density root partition.");
|
|
const Real value = p.Value(left + (right - left) / 2);
|
|
const bool negative = value < 0;
|
|
// Two extra interior signs catch missed roots without assigning a
|
|
// finite-width pocket to an exactly tangent double root.
|
|
for (const Real fraction : {Real{0.25}, Real{0.75}}) {
|
|
const Real check = p.Value(left + fraction * (right - left));
|
|
if ((negative && check > 0) || (!negative && check < 0))
|
|
throw std::runtime_error("Density slice root partition has inconsistent interval signs.");
|
|
}
|
|
result.push_back({left, right, negative});
|
|
}
|
|
return result;
|
|
}
|
|
|
|
inline std::ofstream Open(const std::filesystem::path& path) {
|
|
if (std::filesystem::exists(path)) throw std::invalid_argument("Density region output already exists: " + path.string());
|
|
std::ofstream out(path);
|
|
out.exceptions(std::ios::failbit | std::ios::badbit);
|
|
out << std::setprecision(20);
|
|
return out;
|
|
}
|
|
|
|
// These are analytic, uniform-measure one-dimensional controls, independent of
|
|
// FEM and geometry. The fractional/cubic data are interpolated at three Gauss
|
|
// nodes, reproducing the underresolved positive-profile mechanism of interest.
|
|
inline void SelfChecks(const std::filesystem::path& path) {
|
|
auto out = Open(path);
|
|
out << "case,negative_length,expected_negative_length,negative_mass,expected_negative_mass,passed\n";
|
|
const auto check = [&](const char* name, const Polynomial& p, Real expectedLength, Real expectedMass) {
|
|
Real length = 0, mass = 0;
|
|
for (const auto& interval : Partition(p)) if (interval.negative) {
|
|
length += interval.right - interval.left;
|
|
mass -= p.Integral(interval.left, interval.right);
|
|
}
|
|
const bool passed = std::abs(length - expectedLength) <= 2e-14L * std::max(Real{1e-12}, expectedLength) &&
|
|
std::abs(mass - expectedMass) <= 2e-13L * std::max(Real{1e-12}, expectedMass);
|
|
out << name << ',' << length << ',' << expectedLength << ',' << mass << ',' << expectedMass << ',' << passed << '\n';
|
|
if (!passed) throw std::runtime_error(std::string("Density region analytic self-check failed: ") + name);
|
|
};
|
|
check("positive_quadratic", {1, 0, 1}, 0, 0);
|
|
check("negative_constant", {-2, 0, 0}, 1, 2);
|
|
check("linear_crossing", {-0.25L, 1, 0}, 0.25L, 0.03125L);
|
|
check("double_root_nonnegative", {0.25L, -1, 1}, 0, 0);
|
|
check("near_linear_positive", {1, 1, 1e-28L}, 0, 0);
|
|
// Exact binary coefficients isolate the root algorithm from the separate
|
|
// conditioning of subtracting a tiny decimal offset from a large value.
|
|
check("near_double_negative", {0.25L - std::ldexp(Real{1}, -40), -1, 1},
|
|
std::ldexp(Real{1}, -19), 4 * std::ldexp(Real{1}, -60) / 3);
|
|
const Real d = std::sqrt(Real{15}) / 10;
|
|
const auto fractional = [](Real t) { return t * std::sqrt(t); };
|
|
const Real y0 = fractional(Real{0.5} - d), y1 = fractional(Real{0.5}), y2 = fractional(Real{0.5} + d);
|
|
const Real a = (y0 + y2 - 2 * y1) / (2 * d * d);
|
|
const Real b = (y2 - y0) / (2 * d) - a;
|
|
const Polynomial fractionalInterpolant{y1 - b / 2 - a / 4, b, a};
|
|
check("gauss_q2_t_power_1p5", fractionalInterpolant,
|
|
0.04548403175047501373888L, 0.0005262295634284655L);
|
|
check("gauss_q2_t_power_3", {0.05L, -0.6L, 1.5L}, std::sqrt(Real{6}) / 15, std::sqrt(Real{6}) / 2250);
|
|
}
|
|
|
|
struct Axis final {
|
|
int coordinate{0};
|
|
int orientation{1};
|
|
int face{-1};
|
|
int faceCount{0};
|
|
};
|
|
|
|
inline std::vector<Axis> ChooseAxes(polytrope_validation::PhysicalState& state) {
|
|
auto& mesh = *state.finiteElements.mesh;
|
|
std::vector<Axis> axes(static_cast<std::size_t>(mesh.GetNE()));
|
|
for (int face = 0; face < mesh.GetNumFaces(); ++face) {
|
|
int left = -1, right = -1;
|
|
mesh.GetFaceElements(face, &left, &right);
|
|
const bool leftStar = left >= 0 && state.isStellar(left);
|
|
const bool rightStar = right >= 0 && state.isStellar(right);
|
|
if (leftStar == rightStar) continue;
|
|
const int element = leftStar ? left : right;
|
|
auto& axis = axes[static_cast<std::size_t>(element)];
|
|
++axis.faceCount;
|
|
if (axis.face >= 0) continue;
|
|
auto* transformation = mesh.GetFaceElementTransformations(face);
|
|
if (!transformation) throw std::runtime_error("Missing density region face transformation.");
|
|
mfem::IntegrationPoint facePoint, point;
|
|
facePoint.Set2(0.5, 0.5);
|
|
if (leftStar) transformation->Loc1.Transform(facePoint, point);
|
|
else transformation->Loc2.Transform(facePoint, point);
|
|
const std::array<double, 3> coordinates{point.x, point.y, point.z};
|
|
int selected = -1;
|
|
for (int d = 0; d < 3; ++d) {
|
|
if (std::abs(coordinates[d]) < 1e-12 || std::abs(coordinates[d] - 1) < 1e-12) {
|
|
if (selected >= 0) throw std::runtime_error("Ambiguous density region reference face.");
|
|
selected = d;
|
|
}
|
|
}
|
|
if (selected < 0) throw std::runtime_error("Density region face is not a tensor-product cell face.");
|
|
axis.coordinate = selected;
|
|
axis.orientation = coordinates[selected] < 0.5 ? 1 : -1;
|
|
axis.face = face;
|
|
}
|
|
return axes;
|
|
}
|
|
|
|
inline mfem::IntegrationPoint Point(const Axis& axis, double u, double v, Real t) {
|
|
std::array<double, 3> coordinate{};
|
|
coordinate[axis.coordinate] = static_cast<double>(axis.orientation > 0 ? t : 1 - t);
|
|
coordinate[(axis.coordinate + 1) % 3] = u;
|
|
coordinate[(axis.coordinate + 2) % 3] = v;
|
|
mfem::IntegrationPoint result;
|
|
result.Set3(coordinate[0], coordinate[1], coordinate[2]);
|
|
return result;
|
|
}
|
|
|
|
struct Measure final {
|
|
Real volume{0}, mass{0}, negativeVolume{0}, negativeMass{0};
|
|
std::size_t lines{0}, negativeIntervals{0}, points{0};
|
|
double minimumRho{std::numeric_limits<double>::infinity()};
|
|
double witnessH{std::numeric_limits<double>::quiet_NaN()};
|
|
double witnessT{std::numeric_limits<double>::quiet_NaN()}, witnessU{witnessT}, witnessV{witnessT};
|
|
std::array<double, 3> witnessPosition{};
|
|
int witnessElement{-1};
|
|
double negativeRadiusMin{std::numeric_limits<double>::infinity()}, negativeRadiusMax{-negativeRadiusMin};
|
|
Real negativeTMin{std::numeric_limits<Real>::infinity()}, negativeTMax{-negativeTMin};
|
|
double maximumReconstructionError{0}, maximumReconstructionTolerance{0};
|
|
|
|
void Add(const Measure& other) {
|
|
volume += other.volume; mass += other.mass;
|
|
negativeVolume += other.negativeVolume; negativeMass += other.negativeMass;
|
|
lines += other.lines; negativeIntervals += other.negativeIntervals; points += other.points;
|
|
if (other.minimumRho < minimumRho) {
|
|
minimumRho = other.minimumRho; witnessH = other.witnessH;
|
|
witnessT = other.witnessT; witnessU = other.witnessU; witnessV = other.witnessV;
|
|
witnessPosition = other.witnessPosition; witnessElement = other.witnessElement;
|
|
}
|
|
negativeRadiusMin = std::min(negativeRadiusMin, other.negativeRadiusMin);
|
|
negativeRadiusMax = std::max(negativeRadiusMax, other.negativeRadiusMax);
|
|
negativeTMin = std::min(negativeTMin, other.negativeTMin);
|
|
negativeTMax = std::max(negativeTMax, other.negativeTMax);
|
|
maximumReconstructionError = std::max(maximumReconstructionError, other.maximumReconstructionError);
|
|
maximumReconstructionTolerance = std::max(maximumReconstructionTolerance, other.maximumReconstructionTolerance);
|
|
}
|
|
};
|
|
|
|
inline Real Observed(Real value) { return std::isfinite(value) ? value : std::numeric_limits<Real>::quiet_NaN(); }
|
|
|
|
inline void Header(std::ostream& out) {
|
|
out << "tangent_lines,negative_intervals,integration_points,volume,mass,negative_volume,negative_mass,"
|
|
"negative_volume_fraction,negative_mass_fraction,minimum_density,witness_h,witness_element,witness_t,witness_u,witness_v,"
|
|
"witness_x,witness_y,witness_z,negative_radius_minimum,negative_radius_maximum,negative_interval_t_minimum,"
|
|
"negative_interval_t_maximum,maximum_reconstruction_error,maximum_reconstruction_tolerance\n";
|
|
}
|
|
|
|
inline void Write(std::ostream& out, const Measure& m) {
|
|
out << m.lines << ',' << m.negativeIntervals << ',' << m.points << ',' << m.volume << ',' << m.mass << ','
|
|
<< m.negativeVolume << ',' << m.negativeMass << ',' << m.negativeVolume / m.volume << ','
|
|
<< (m.mass != 0 ? m.negativeMass / std::abs(m.mass) : std::numeric_limits<Real>::quiet_NaN()) << ','
|
|
<< m.minimumRho << ',' << m.witnessH << ',' << m.witnessElement << ',' << m.witnessT << ',' << m.witnessU << ',' << m.witnessV;
|
|
for (double x : m.witnessPosition) out << ',' << x;
|
|
out << ',' << Observed(m.negativeRadiusMin) << ',' << Observed(m.negativeRadiusMax) << ','
|
|
<< Observed(m.negativeTMin) << ',' << Observed(m.negativeTMax) << ','
|
|
<< m.maximumReconstructionError << ',' << m.maximumReconstructionTolerance << '\n';
|
|
}
|
|
|
|
inline bool Identical(const mfem::Vector& a, const mfem::Vector& b) {
|
|
return a.Size() == b.Size() && (a.Size() == 0 || std::memcmp(a.GetData(), b.GetData(), sizeof(mfem::real_t) * a.Size()) == 0);
|
|
}
|
|
|
|
} // namespace negative_density_regions_detail
|
|
|
|
inline void RecordDensityNegativeRegions(polytrope_validation::PhysicalState& state, double n, double K,
|
|
const std::filesystem::path& output) {
|
|
namespace d = negative_density_regions_detail;
|
|
using d::Real;
|
|
if (!std::isfinite(n) || n < 1 || !std::isfinite(K) || K <= 0)
|
|
throw std::invalid_argument("Density regions require finite n>=1 and K>0.");
|
|
std::filesystem::create_directories(output);
|
|
// Fail the one-dimensional arithmetic controls before any FEM evaluation.
|
|
d::SelfChecks(output / "density_negative_regions_selfchecks.csv");
|
|
const auto& fem = state.finiteElements;
|
|
if (!fem.mesh || !fem.densityFes || !fem.enthalpyFes || !fem.compactificationCoordinate ||
|
|
fem.mesh->Dimension() != 3 || !fem.mesh->Conforming() || fem.densityFes->GetVDim() != 1 ||
|
|
dynamic_cast<const mfem::L2_FECollection*>(fem.densityFes->FEColl()) == nullptr ||
|
|
state.density.ParFESpace() != fem.densityFes.get() || state.enthalpy.ParFESpace() != fem.enthalpyFes.get())
|
|
throw std::invalid_argument("Density regions require the current scalar L2 density space on a conforming 3D mesh.");
|
|
int ranks = 0;
|
|
MPI_Comm_size(fem.mesh->GetComm(), &ranks);
|
|
if (ranks != 1) throw std::invalid_argument("Density regions currently require one MPI rank.");
|
|
for (int element = 0; element < fem.mesh->GetNE(); ++element) if (state.isStellar(element)) {
|
|
const auto* fe = fem.densityFes->GetFE(element);
|
|
if (!fe || fe->GetOrder() != 2 || fe->GetMapType() != mfem::FiniteElement::VALUE ||
|
|
fe->GetGeomType() != mfem::Geometry::CUBE || fe->GetDof() != 27)
|
|
throw std::invalid_argument("Density root splitting requires scalar VALUE Q2 hexahedral density elements.");
|
|
}
|
|
const auto axes = d::ChooseAxes(state);
|
|
const mfem::Vector savedDensity(state.density), savedEnthalpy(state.enthalpy), savedPotential(state.potential),
|
|
savedGravity(state.gravityGradientReference), savedDisplacement(state.displacement);
|
|
mean_field::mapping::GridFunctionMappingEvaluator mapper(state.domainMapper, state.displacement, *fem.compactificationCoordinate);
|
|
mean_field::mapping::MappingPointContext mapped;
|
|
auto elements = d::Open(output / "density_negative_regions_elements.csv");
|
|
elements << "tangential_order,interval_order,element,attribute,axis,orientation,surface_face,surface_face_count,";
|
|
d::Header(elements);
|
|
auto summary = d::Open(output / "density_negative_regions_summary.csv");
|
|
summary << "tangential_order,interval_order,";
|
|
d::Header(summary);
|
|
auto intervals = d::Open(output / "density_negative_regions_intervals.csv");
|
|
intervals << "tangential_order,element,attribute,axis,orientation,surface_face,u,v,tangential_weight,t_left,t_right,"
|
|
"coefficient_c,coefficient_b,coefficient_a,negative_volume,negative_mass,minimum_density,witness_h,"
|
|
"radius_minimum,radius_maximum,left_x,left_y,left_z,right_x,right_y,right_z\n";
|
|
constexpr int intervalOrder = 14;
|
|
const auto& intervalRule = mfem::IntRules.Get(mfem::Geometry::SEGMENT, intervalOrder);
|
|
if (intervalRule.GetNPoints() != 8) throw std::logic_error("Unexpected density interval Gauss rule.");
|
|
for (int tangentialOrder : {14, 22, 30}) {
|
|
const auto& tangentRule = mfem::IntRules.Get(mfem::Geometry::SEGMENT, tangentialOrder);
|
|
if (tangentRule.GetNPoints() != tangentialOrder / 2 + 1)
|
|
throw std::logic_error("Unexpected density tangential Gauss rule.");
|
|
d::Measure global;
|
|
for (int element = 0; element < fem.mesh->GetNE(); ++element) {
|
|
if (!state.isStellar(element)) continue;
|
|
const auto& axis = axes[static_cast<std::size_t>(element)];
|
|
auto* transformation = fem.mesh->GetElementTransformation(element);
|
|
d::Measure local;
|
|
for (int iu = 0; iu < tangentRule.GetNPoints(); ++iu) for (int iv = 0; iv < tangentRule.GetNPoints(); ++iv) {
|
|
const auto& up = tangentRule.IntPoint(iu);
|
|
const auto& vp = tangentRule.IntPoint(iv);
|
|
const double u = up.x, v = vp.x;
|
|
const Real tangentWeight = static_cast<Real>(up.weight) * vp.weight;
|
|
const auto value = [&](Real t) {
|
|
const auto point = d::Point(axis, u, v, t);
|
|
const double rho = state.density.GetValue(element, point);
|
|
if (!std::isfinite(rho)) throw std::domain_error("Nonfinite density slice value.");
|
|
return rho;
|
|
};
|
|
const auto polynomial = d::Interpolate(value(0), value(Real{0.5}), value(1));
|
|
const Real tolerance = 512 * std::numeric_limits<double>::epsilon() * polynomial.Scale();
|
|
const auto verify = [&](Real t, double actual) {
|
|
const Real error = std::abs(polynomial.Value(t) - actual);
|
|
local.maximumReconstructionError = std::max(local.maximumReconstructionError, static_cast<double>(error));
|
|
local.maximumReconstructionTolerance = std::max(local.maximumReconstructionTolerance, static_cast<double>(tolerance));
|
|
if (!std::isfinite(error) || error > tolerance)
|
|
throw std::runtime_error("Frozen density slice fails independent Q2 reconstruction.");
|
|
};
|
|
for (Real t : {Real{0.125}, Real{0.25}, Real{0.375}, Real{0.625}, Real{0.75}, Real{0.875}})
|
|
verify(t, value(t));
|
|
const auto parts = d::Partition(polynomial);
|
|
++local.lines;
|
|
const auto evaluate = [&](Real t, d::Measure& measure, bool negativeRegion) {
|
|
const auto point = d::Point(axis, u, v, t);
|
|
if (mapper.EvaluatePoint(*transformation, point, mapped) != mean_field::mapping::MappingStatus::valid || mapped.compactified)
|
|
throw std::domain_error("Invalid stellar density-region mapping.");
|
|
const double rho = state.density.GetValue(element, point), h = state.enthalpy.GetValue(element, point);
|
|
if (!std::isfinite(rho) || !std::isfinite(h)) throw std::domain_error("Nonfinite density-region witness.");
|
|
verify(t, rho);
|
|
if (rho < measure.minimumRho) {
|
|
measure.minimumRho = rho; measure.witnessH = h; measure.witnessElement = element;
|
|
measure.witnessT = static_cast<double>(t); measure.witnessU = u; measure.witnessV = v;
|
|
for (int component = 0; component < 3; ++component) measure.witnessPosition[component] = mapped.physical_position(component);
|
|
}
|
|
if (negativeRegion) {
|
|
const double radius = mapped.physical_position.Norml2();
|
|
if (!std::isfinite(radius)) throw std::domain_error("Nonfinite density-region radius.");
|
|
measure.negativeRadiusMin = std::min(measure.negativeRadiusMin, radius);
|
|
measure.negativeRadiusMax = std::max(measure.negativeRadiusMax, radius);
|
|
}
|
|
transformation->SetIntPoint(&point);
|
|
const double referenceDeterminant = transformation->Jacobian().Det();
|
|
const Real determinant = static_cast<Real>(referenceDeterminant) * mapped.mapping_determinant;
|
|
if (!(referenceDeterminant > 0) || !std::isfinite(referenceDeterminant) ||
|
|
!(mapped.mapping_determinant > 0) || !std::isfinite(mapped.mapping_determinant) ||
|
|
!(determinant > 0) || !std::isfinite(determinant))
|
|
throw std::domain_error("Nonpositive or nonfinite density-region physical measure.");
|
|
return determinant;
|
|
};
|
|
for (const auto& part : parts) {
|
|
d::Measure piece;
|
|
std::array<double, 3> leftPosition{}, rightPosition{};
|
|
(void)evaluate(part.left, piece, part.negative);
|
|
for (int component = 0; component < 3; ++component) leftPosition[component] = mapped.physical_position(component);
|
|
(void)evaluate(part.right, piece, part.negative);
|
|
for (int component = 0; component < 3; ++component) rightPosition[component] = mapped.physical_position(component);
|
|
if (polynomial.a != 0) {
|
|
const Real stationary = -polynomial.b / (2 * polynomial.a);
|
|
if (stationary > part.left && stationary < part.right) (void)evaluate(stationary, piece, part.negative);
|
|
}
|
|
for (int it = 0; it < intervalRule.GetNPoints(); ++it) {
|
|
const auto& qp = intervalRule.IntPoint(it);
|
|
const Real t = part.left + (part.right - part.left) * qp.x;
|
|
const Real determinant = evaluate(t, piece, part.negative);
|
|
const Real weight = tangentWeight * (part.right - part.left) * qp.weight * determinant;
|
|
const Real rho = polynomial.Value(t);
|
|
if (!(weight > 0) || !std::isfinite(weight) || !std::isfinite(rho) ||
|
|
(part.negative && rho > 0) || (!part.negative && rho < 0))
|
|
throw std::runtime_error("Density region quadrature disagrees with its root partition.");
|
|
piece.volume += weight;
|
|
piece.mass += weight * rho;
|
|
++piece.points;
|
|
}
|
|
if (part.negative) {
|
|
piece.negativeIntervals = 1;
|
|
piece.negativeVolume = piece.volume;
|
|
piece.negativeMass = -piece.mass;
|
|
piece.negativeTMin = part.left;
|
|
piece.negativeTMax = part.right;
|
|
if (tangentialOrder == 30) {
|
|
intervals << tangentialOrder << ',' << element << ',' << fem.mesh->GetAttribute(element) << ','
|
|
<< axis.coordinate << ',' << axis.orientation << ',' << axis.face << ',' << u << ',' << v << ','
|
|
<< tangentWeight << ',' << part.left << ',' << part.right << ',' << polynomial.c << ',' << polynomial.b << ','
|
|
<< polynomial.a << ',' << piece.negativeVolume << ',' << piece.negativeMass << ',' << piece.minimumRho << ','
|
|
<< piece.witnessH << ',' << piece.negativeRadiusMin << ',' << piece.negativeRadiusMax;
|
|
for (double x : leftPosition) intervals << ',' << x;
|
|
for (double x : rightPosition) intervals << ',' << x;
|
|
intervals << '\n';
|
|
}
|
|
}
|
|
local.Add(piece);
|
|
}
|
|
}
|
|
if (!(local.volume > 0) || !std::isfinite(local.volume) || !std::isfinite(local.mass) ||
|
|
!std::isfinite(local.negativeVolume) || !std::isfinite(local.negativeMass) || local.negativeMass < 0)
|
|
throw std::runtime_error("Invalid integrated density-region measure.");
|
|
elements << tangentialOrder << ',' << intervalOrder << ',' << element << ',' << fem.mesh->GetAttribute(element) << ','
|
|
<< axis.coordinate << ',' << axis.orientation << ',' << axis.face << ',' << axis.faceCount << ',';
|
|
d::Write(elements, local);
|
|
global.Add(local);
|
|
}
|
|
if (!(global.volume > 0)) throw std::invalid_argument("Density region diagnostic found no stellar volume.");
|
|
summary << tangentialOrder << ',' << intervalOrder << ',';
|
|
d::Write(summary, global);
|
|
summary.flush(); elements.flush(); intervals.flush();
|
|
}
|
|
const bool unchanged = d::Identical(savedDensity, state.density) && d::Identical(savedEnthalpy, state.enthalpy) &&
|
|
d::Identical(savedPotential, state.potential) && d::Identical(savedGravity, state.gravityGradientReference) &&
|
|
d::Identical(savedDisplacement, state.displacement);
|
|
if (!unchanged) throw std::logic_error("Density region observation changed a source field.");
|
|
auto scope = d::Open(output / "density_negative_regions_scope.txt");
|
|
scope << "source_fields_bitwise_identical=1\ntrial_preparations=0\nmpi_ranks=1\n"
|
|
<< "n=" << n << "\nK=" << K << "\ndensity_basis=" << fem.densityFes->FEColl()->Name() << '\n'
|
|
<< "Mass uses the reconstructed Q2 polynomial, verified against direct FE values; no density or EOS clipping.\n"
|
|
<< "Volume weights contain one reference determinant and one deformation determinant.\n"
|
|
<< "All stellar cells are included; surface axes point inward from the lowest-index stellar/nonstellar face.\n"
|
|
<< "Other cells use reference axis0. Tangential orders14/22/30 use8/12/16 points per axis; interval order14 uses8.\n"
|
|
<< "Whole-cell denominators integrate every root-partition interval with the same sliced quadrature.\n"
|
|
<< "Negative intervals are exact for each reconstructed Q2 line up to floating arithmetic; tangential integration is sampled.\n"
|
|
<< "Witnesses and radial bounds are samples, not global extrema or polynomial sign certificates.\n"
|
|
<< "Near-double-root regions below reconstruction uncertainty are not certified; compare quadrature levels and errors.\n";
|
|
}
|
|
|
|
} // namespace experiment::plateau_diagnostics
|