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

Add controlled density-order projection and coupled-solve observations

Browse Source
This commit is contained in:
Emily Boudreaux
2026-09-19 10:48:14 -04:00
parent dd955abca5
commit 2ca6890a18
5 changed files with 970 additions and 1 deletions
Download Patch File Download Diff File Expand all files Collapse all files

117
experiments/density_order_coupled_cases_2026-09-19.json Normal file
View File

@@ -0,0 +1,117 @@
[
{
"id": "n1_J0",
"arguments": [
"--index",
"1",
"--angular-momentum",
"0",
"--surface-prescription",
"harmonics",
"--surface-degree",
"4",
"--max-newton",
"8",
"--restart",
"160",
"--max-linear",
"200",
"--linear-tol",
"0.03",
"--rtol",
"1e-8",
"--atol",
"0",
"--strategy",
"baseline",
"--density-order-final",
"--no-central"
]
},
{
"id": "n1p5_J0",
"arguments": [
"--index",
"1.5",
"--angular-momentum",
"0",
"--surface-prescription",
"harmonics",
"--surface-degree",
"4",
"--max-newton",
"8",
"--restart",
"160",
"--max-linear",
"200",
"--linear-tol",
"0.03",
"--rtol",
"1e-8",
"--atol",
"0",
"--strategy",
"baseline",
"--density-order-final",
"--no-central"
]
},
{
"id": "n3_J0",
"arguments": [
"--index",
"3",
"--angular-momentum",
"0",
"--surface-prescription",
"harmonics",
"--surface-degree",
"4",
"--max-newton",
"8",
"--restart",
"160",
"--max-linear",
"200",
"--linear-tol",
"0.03",
"--rtol",
"1e-8",
"--atol",
"0",
"--strategy",
"baseline",
"--density-order-final"
]
},
{
"id": "n1_J0p01",
"arguments": [
"--index",
"1",
"--angular-momentum",
"0.01",
"--surface-prescription",
"harmonics",
"--surface-degree",
"4",
"--max-newton",
"8",
"--restart",
"160",
"--max-linear",
"200",
"--linear-tol",
"0.03",
"--rtol",
"1e-8",
"--atol",
"0",
"--strategy",
"baseline",
"--density-order-final",
"--no-central"
]
}
]

56
experiments/density_order_frozen_cases_2026-09-19.json Normal file
View File

@@ -0,0 +1,56 @@
[
{
"id": "n1_equilibrium",
"arguments": [
"--index",
"1",
"--angular-momentum",
"0",
"--surface-prescription",
"harmonics",
"--surface-degree",
"4",
"--no-central",
"--max-newton",
"0",
"--density-order-projection",
"--density-snapshot",
"../../diagnostics/spherical_harmonic_surface_2026-09-19/campaign/n1_J0_L4"
]
},
{
"id": "n1p5_equilibrium",
"arguments": [
"--index",
"1.5",
"--angular-momentum",
"0",
"--surface-prescription",
"harmonics",
"--surface-degree",
"4",
"--no-central",
"--max-newton",
"0",
"--density-order-projection",
"--density-snapshot",
"../../diagnostics/spherical_harmonic_surface_2026-09-19/campaign/n1p5_J0_L4"
]
},
{
"id": "n3_seed",
"arguments": [
"--index",
"3",
"--angular-momentum",
"0",
"--surface-prescription",
"harmonics",
"--surface-degree",
"4",
"--max-newton",
"0",
"--density-order-projection"
]
}
]

384
experiments/density_order_projection.hpp Normal file
View File

@@ -0,0 +1,384 @@
#pragma once
// Frozen-field order controls, not solver recipes. All target spaces borrow the
// same mesh and retain the source enthalpy and deformation without modification.
#include "negative_density_projection.hpp"
#include "density_polynomial_regions.hpp"
#include <iostream>
#include <chrono>
#include <map>
#include <span>
namespace experiment::plateau_diagnostics {
namespace density_order_detail {
namespace projection = negative_density_projection_detail;
using ObserverClock = std::chrono::steady_clock;
inline double Seconds(ObserverClock::time_point start, ObserverClock::time_point finish) {
return std::chrono::duration<double>(finish - start).count();
}
inline void WriteTimings(const std::filesystem::path& path, ObserverClock::time_point start,
ObserverClock::time_point prepared, ObserverClock::time_point errors,
ObserverClock::time_point finished) {
auto out = region::Open(path);
out << "phase,seconds\nsetup_and_projection," << Seconds(start, prepared)
<< "\nindependent_error_quadrature," << Seconds(prepared, errors)
<< "\nregion_observation_and_output," << Seconds(errors, finished)
<< "\nobserver_total," << Seconds(start, finished) << '\n';
}
inline int ValidateState(polytrope_validation::PhysicalState& state, double n, double K) {
const auto& fem = state.finiteElements;
if (!std::isfinite(n) || n < 1 || !std::isfinite(K) || K <= 0 || !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 order observations require the current conforming 3D scalar L2 state and finite n>=1,K>0.");
int ranks = 0;
MPI_Comm_size(fem.mesh->GetComm(), &ranks);
if (ranks != 1) throw std::invalid_argument("Density order observations require one MPI rank.");
int degree = -1;
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->GetOrder() > 4 || fe->GetMapType() != mfem::FiniteElement::VALUE ||
fe->GetGeomType() != mfem::Geometry::CUBE || fe->GetDof() != (fe->GetOrder() + 1) * (fe->GetOrder() + 1) * (fe->GetOrder() + 1))
throw std::invalid_argument("Density order observations require scalar VALUE tensor Q2/Q3/Q4 elements.");
if (degree >= 0 && degree != fe->GetOrder()) throw std::invalid_argument("Variable density order is outside this frozen study.");
degree = fe->GetOrder();
}
if (degree < 0) throw std::invalid_argument("Density order observation found no stellar elements.");
return degree;
}
struct SourceCheckpoint final {
mfem::Vector density, enthalpy, potential, gravity, displacement;
explicit SourceCheckpoint(const polytrope_validation::PhysicalState& state)
: density(state.density), enthalpy(state.enthalpy), potential(state.potential),
gravity(state.gravityGradientReference), displacement(state.displacement) {}
void Verify(const polytrope_validation::PhysicalState& state) const {
if (!region::Identical(density, state.density) || !region::Identical(enthalpy, state.enthalpy) ||
!region::Identical(potential, state.potential) || !region::Identical(gravity, state.gravityGradientReference) ||
!region::Identical(displacement, state.displacement))
throw std::logic_error("Density order diagnostic changed a source field.");
}
};
inline void WriteScope(const std::filesystem::path& path, bool projections, int sourceDegree, double n, double K) {
auto out = region::Open(path);
out << "source_fields_bitwise_identical=1\ntrial_preparations=0\nmpi_ranks=1\nsource_density_degree=" << sourceDegree
<< "\nn=" << n << "\nK=" << K << "\nprojection_controls=" << projections << '\n'
<< "All measures use unchanged source enthalpy and physical geometry.\n"
<< "eos_nodal_interpolation is interpolation of rho_EOS(h_h), NOT original Lane-Emden initialization.\n"
<< "eos_l2_q26 uses one fixed order26 rule for both mass and target source at each density order.\n"
<< "Surface records are one-sided samples, not area fractions.\n"
<< "Density coefficients are in different FE spaces; their array differences across orders are not physical norms.\n"
<< "EOS errors and moments are weighted physical integrals; full and negative measures share root-sliced quadrature.\n"
<< "density_order_error_quadrature.csv independently compares unsliced order26/30 EOS error; it should be used for error-norm sensitivity.\n"
<< "Complete EOS integrals are undefined if negative h is sampled; density mass and negative-region measures remain reported.\n"
<< "Projection/interpolation construction rejects negative sampled enthalpy; it does not repair the source.\n"
<< "Root-sliced measures are numerical estimates, not global positivity certificates.\n"
<< "negative_mass_uncertainty_scale is the integral of a coefficient-scaled reconstruction tolerance, not a rigorous bound.\n"
<< "Resolved-negative intervals have a direct FE witness below their reconstruction tolerance; their full volume is still an estimate.\n"
<< "Uncertain shallow intervals and surface values are retained and separately labeled; no density field is clipped.\n"
<< "Near repeated roots and narrow unsampled tangential pockets remain unresolved limitations.\n";
}
struct DegreeFields final {
int degree;
std::unique_ptr<mfem::L2_FECollection> collection;
std::unique_ptr<mfem::ParFiniteElementSpace> space;
std::unique_ptr<mfem::ParGridFunction> interpolation, projection, constant;
mfem::DenseMatrix shapes;
DegreeFields(mfem::ParMesh& mesh, int order, const mfem::IntegrationRule& rule)
: degree(order), collection(std::make_unique<mfem::L2_FECollection>(order, 3, mfem::BasisType::GaussLegendre, mfem::FiniteElement::VALUE)),
space(std::make_unique<mfem::ParFiniteElementSpace>(&mesh, collection.get())),
interpolation(std::make_unique<mfem::ParGridFunction>(space.get())),
projection(std::make_unique<mfem::ParGridFunction>(space.get())),
constant(std::make_unique<mfem::ParGridFunction>(space.get())) {
*interpolation = 0.; *projection = 0.;
mfem::ConstantCoefficient one(1.);
constant->ProjectCoefficient(one);
const auto* fe = collection->FiniteElementForGeometry(mfem::Geometry::CUBE);
shapes.SetSize(fe->GetDof(), rule.GetNPoints());
mfem::Vector shape(fe->GetDof());
for (int q = 0; q < rule.GetNPoints(); ++q) {
fe->CalcShape(rule.IntPoint(q), shape);
for (int i = 0; i < shape.Size(); ++i) shapes(i, q) = shape(i);
}
}
};
struct ErrorField final {
std::string_view label;
int degree;
const mfem::ParGridFunction* density;
};
struct ErrorMeasure final {
Real volume{0}, mass{0}, targetMass{0}, eosSquared{0}, signedEos{0};
Real minimumRho{std::numeric_limits<Real>::infinity()};
double minimumH{std::numeric_limits<double>::infinity()};
std::size_t negativeH{0};
};
inline void ErrorHeader(std::ostream& out) {
out << "control,degree,element,attribute,quadrature_request,quadrature_achieved,quadrature_points,volume,mass,target_mass,"
"eos_squared_integral,signed_eos_integral,minimum_density,minimum_enthalpy,negative_enthalpy_samples,eos_defined_everywhere_sampled\n";
}
inline void WriteError(std::ostream& out, std::string_view label, int degree, int element, int attribute,
int order, const mfem::IntegrationRule& rule, const ErrorMeasure& measure) {
const Real nan = std::numeric_limits<Real>::quiet_NaN();
if (!(measure.volume > 0) || !std::isfinite(measure.volume) || !std::isfinite(measure.mass) ||
!std::isfinite(measure.minimumRho) || !std::isfinite(measure.minimumH) ||
(measure.negativeH == 0 && (!std::isfinite(measure.targetMass) || !std::isfinite(measure.eosSquared) || !std::isfinite(measure.signedEos))))
throw std::runtime_error("Nonfinite independent density error measure.");
out << label << ',' << degree << ',' << element << ',' << attribute << ',' << order << ',' << rule.GetOrder() << ','
<< rule.GetNPoints() << ',' << measure.volume << ',' << measure.mass << ','
<< (measure.negativeH == 0 ? measure.targetMass : nan) << ',' << (measure.negativeH == 0 ? measure.eosSquared : nan) << ','
<< (measure.negativeH == 0 ? measure.signedEos : nan) << ',' << measure.minimumRho << ',' << measure.minimumH << ','
<< measure.negativeH << ',' << (measure.negativeH == 0) << '\n';
}
// A squared EOS error needs more quadrature than either projection moment.
// This pass reuses each mapped sample and reference basis across all supplied
// fields; it does not assemble another matrix or solve another projection.
inline void ObserveErrorRule(polytrope_validation::PhysicalState& state, std::span<const ErrorField> fields,
int order, double n, double K, std::ostream& out) {
const auto& fem = state.finiteElements;
const auto& rule = mfem::IntRules.Get(mfem::Geometry::CUBE, order);
std::map<const mfem::FiniteElement*, std::unique_ptr<mfem::DenseMatrix>> shapeCache;
std::vector<const mfem::DenseMatrix*> shapes(fields.size());
std::vector<mfem::Vector> coefficients(fields.size());
mean_field::mapping::GridFunctionMappingEvaluator mapper(state.domainMapper, state.displacement, *fem.compactificationCoordinate);
mean_field::mapping::MappingPointContext mapped;
for (int element = 0; element < fem.mesh->GetNE(); ++element) {
if (!state.isStellar(element)) continue;
auto* transformation = fem.mesh->GetElementTransformation(element);
for (std::size_t field = 0; field < fields.size(); ++field) {
const auto& density = *fields[field].density;
const auto* fe = density.ParFESpace()->GetFE(element);
if (fe->GetOrder() != fields[field].degree || fe->GetMapType() != mfem::FiniteElement::VALUE || fe->GetGeomType() != mfem::Geometry::CUBE)
throw std::logic_error("Independent density error quadrature received an incompatible scalar field.");
auto& cached = shapeCache[fe];
if (!cached) {
cached = std::make_unique<mfem::DenseMatrix>(fe->GetDof(), rule.GetNPoints());
mfem::Vector shape(fe->GetDof());
for (int q = 0; q < rule.GetNPoints(); ++q) {
fe->CalcShape(rule.IntPoint(q), shape);
for (int i = 0; i < shape.Size(); ++i) (*cached)(i, q) = shape(i);
}
}
shapes[field] = cached.get();
mfem::Array<int> dofs;
auto* dofTransform = density.ParFESpace()->GetElementDofs(element, dofs);
density.GetSubVector(dofs, coefficients[field]);
if (dofTransform) dofTransform->InvTransformPrimal(coefficients[field]);
}
std::vector<ErrorMeasure> measures(fields.size());
for (int q = 0; q < rule.GetNPoints(); ++q) {
const auto& point = rule.IntPoint(q);
if (mapper.EvaluatePoint(*transformation, point, mapped) != mean_field::mapping::MappingStatus::valid || mapped.compactified)
throw std::domain_error("Invalid independent density error mapping.");
const double h = state.enthalpy.GetValue(element, point);
const Real target = EosTarget(h, n, K);
transformation->SetIntPoint(&point);
const double reference = transformation->Jacobian().Det();
const Real weight = static_cast<Real>(point.weight) * reference * mapped.mapping_determinant;
if (!(reference > 0) || !std::isfinite(reference) || !(mapped.mapping_determinant > 0) ||
!std::isfinite(mapped.mapping_determinant) || !(weight > 0) || !std::isfinite(weight))
throw std::domain_error("Invalid independent density error weight.");
for (std::size_t field = 0; field < fields.size(); ++field) {
auto& measure = measures[field];
Real rho = 0;
const auto* shape = shapes[field]->GetColumn(q);
for (int i = 0; i < coefficients[field].Size(); ++i) rho += static_cast<Real>(shape[i]) * coefficients[field](i);
if (!std::isfinite(rho)) throw std::domain_error("Nonfinite independently observed density.");
measure.volume += weight; measure.mass += weight * rho; measure.targetMass += weight * target;
measure.eosSquared += weight * (rho - target) * (rho - target); measure.signedEos += weight * (rho - target);
measure.minimumRho = std::min(measure.minimumRho, rho); measure.minimumH = std::min(measure.minimumH, h);
measure.negativeH += h < 0;
}
}
for (std::size_t field = 0; field < fields.size(); ++field)
WriteError(out, fields[field].label, fields[field].degree, element, fem.mesh->GetAttribute(element), order, rule, measures[field]);
}
out.flush();
}
} // namespace density_order_detail
inline void RecordDensityOrderProjection(polytrope_validation::PhysicalState& state, double n, double K,
const std::filesystem::path& output) {
namespace d = density_order_detail;
const auto observerStart = d::ObserverClock::now();
const int sourceDegree = d::ValidateState(state, n, K);
std::filesystem::create_directories(output);
d::RootSelfChecks(output / "density_order_root_selfchecks.csv");
d::SourceCheckpoint checkpoint(state);
const auto& fem = state.finiteElements;
const auto& rule = mfem::IntRules.Get(mfem::Geometry::CUBE, 26);
std::vector<std::unique_ptr<d::DegreeFields>> controls;
for (int degree : {2, 3, 4}) controls.push_back(std::make_unique<d::DegreeFields>(*fem.mesh, degree, rule));
auto elements = d::region::Open(output / "density_order_projection_elements.csv");
elements << "degree,element,attribute,quadrature_request,quadrature_achieved,quadrature_points,dofs,basis_family,volume,target_mass,"
"interpolation_mass,projection_mass,interpolation_eos_squared,projection_eos_squared,rhs_norm,normal_defect_norm,"
"normal_relative_defect,constant_coefficient_error,constant_sample_error,interpolation_coefficient_minimum,"
"projection_coefficient_minimum,projection_coefficient_maximum\n";
auto nodes = d::region::Open(output / "density_order_nodes.csv");
nodes << "degree,element,local_dof,global_dof,xi,eta,zeta,h,target_rho,interpolation_coefficient,projection_coefficient\n";
auto errorQuadrature = d::region::Open(output / "density_order_error_quadrature.csv");
d::ErrorHeader(errorQuadrature);
mean_field::mapping::GridFunctionMappingEvaluator mapper(state.domainMapper, state.displacement, *fem.compactificationCoordinate);
mean_field::mapping::MappingPointContext mapped;
mfem::Vector weights(rule.GetNPoints()), targets(rule.GetNPoints());
for (int element = 0; element < fem.mesh->GetNE(); ++element) {
if (!state.isStellar(element)) continue;
auto* transformation = fem.mesh->GetElementTransformation(element);
// Geometry and EOS evaluation dominate observer cost; share these
// samples across target orders without sharing their independent M,b.
long double volume = 0, targetMass = 0;
double minimumH = std::numeric_limits<double>::infinity();
for (int q = 0; q < rule.GetNPoints(); ++q) {
const auto& point = rule.IntPoint(q);
if (mapper.EvaluatePoint(*transformation, point, mapped) != mean_field::mapping::MappingStatus::valid || mapped.compactified)
throw std::domain_error("Invalid density order projection mapping.");
const double h = state.enthalpy.GetValue(element, point);
if (!std::isfinite(h) || h < 0) throw std::domain_error("Density order projection encountered inadmissible sampled enthalpy.");
targets(q) = d::EosTarget(h, n, K);
minimumH = std::min(minimumH, h);
transformation->SetIntPoint(&point);
const double reference = transformation->Jacobian().Det();
weights(q) = point.weight * reference * mapped.mapping_determinant;
if (!(reference > 0) || !std::isfinite(reference) || !(mapped.mapping_determinant > 0) ||
!std::isfinite(mapped.mapping_determinant) || !(weights(q) > 0) || !std::isfinite(weights(q)))
throw std::domain_error("Invalid density order projection weight.");
volume += weights(q); targetMass += static_cast<long double>(weights(q)) * targets(q);
}
for (auto& control : controls) {
const auto* fe = control->space->GetFE(element);
const int size = fe->GetDof();
mfem::Array<int> dofs;
if (control->space->GetElementDofs(element, dofs) != nullptr || dofs.Size() != size)
throw std::invalid_argument("Density order projection requires untransformed scalar DG element DOFs.");
for (int i = 0; i < size; ++i) if (dofs[i] < 0) throw std::logic_error("Unexpected signed DG density DOF.");
mfem::Vector interpolated(size), nodalH(size), rhs(size), constantRhs(size), constantExpected;
const auto& nodalRule = fe->GetNodes();
if (nodalRule.GetNPoints() != size) throw std::logic_error("Density target collection is not nodal.");
for (int i = 0; i < size; ++i) {
nodalH(i) = state.enthalpy.GetValue(element, nodalRule.IntPoint(i));
if (!std::isfinite(nodalH(i)) || nodalH(i) < 0)
throw std::domain_error("Density order nodal interpolation encountered inadmissible enthalpy.");
interpolated(i) = d::EosTarget(nodalH(i), n, K);
}
mfem::DenseMatrix mass(size);
mass = 0.; rhs = 0.; constantRhs = 0.;
for (int q = 0; q < rule.GetNPoints(); ++q) {
const auto* shape = control->shapes.GetColumn(q);
for (int j = 0; j < size; ++j) {
const double weighted = weights(q) * shape[j];
rhs(j) += weighted * targets(q); constantRhs(j) += weighted;
auto* column = mass.GetColumn(j);
for (int i = 0; i <= j; ++i) column[i] += weighted * shape[i];
}
}
for (int j = 0; j < size; ++j) for (int i = 0; i < j; ++i) mass(j, i) = mass(i, j);
const auto lower = d::projection::Factor(mass);
const auto projected = d::projection::Solve(lower, rhs), constant = d::projection::Solve(lower, constantRhs);
const auto defect = d::projection::Defect(mass, projected, rhs);
d::projection::VerifySolve(mass, projected, rhs, defect);
d::projection::VerifySolve(mass, constant, constantRhs, d::projection::Defect(mass, constant, constantRhs));
control->constant->GetSubVector(dofs, constantExpected);
double constantError = 0, constantSampleError = 0;
for (int i = 0; i < size; ++i) constantError = std::max(constantError, std::abs(constant(i) - constantExpected(i)));
long double interpolationMass = 0, projectionMass = 0, interpolationError = 0, projectionError = 0;
d::ErrorMeasure interpolationMeasure, projectionMeasure;
interpolationMeasure.volume = projectionMeasure.volume = volume;
interpolationMeasure.targetMass = projectionMeasure.targetMass = targetMass;
interpolationMeasure.minimumH = projectionMeasure.minimumH = minimumH;
for (int q = 0; q < rule.GetNPoints(); ++q) {
long double interpolatedValue = 0, projectedValue = 0, constantValue = 0;
for (int i = 0; i < size; ++i) {
const double shape = control->shapes(i, q);
interpolatedValue += static_cast<long double>(shape) * interpolated(i);
projectedValue += static_cast<long double>(shape) * projected(i);
constantValue += static_cast<long double>(shape) * constant(i);
}
interpolationMass += weights(q) * interpolatedValue; projectionMass += weights(q) * projectedValue;
interpolationError += weights(q) * (interpolatedValue - targets(q)) * (interpolatedValue - targets(q));
projectionError += weights(q) * (projectedValue - targets(q)) * (projectedValue - targets(q));
interpolationMeasure.signedEos += weights(q) * (interpolatedValue - targets(q));
projectionMeasure.signedEos += weights(q) * (projectedValue - targets(q));
interpolationMeasure.minimumRho = std::min(interpolationMeasure.minimumRho, interpolatedValue);
projectionMeasure.minimumRho = std::min(projectionMeasure.minimumRho, projectedValue);
constantSampleError = std::max(constantSampleError, static_cast<double>(std::abs(constantValue - 1)));
}
if (!std::isfinite(constantError) || !std::isfinite(constantSampleError) || constantError > 1e-10 || constantSampleError > 1e-10)
throw std::runtime_error("Density order projection failed independent constant reconstruction.");
control->interpolation->SetSubVector(dofs, interpolated); control->projection->SetSubVector(dofs, projected);
interpolationMeasure.mass = interpolationMass; interpolationMeasure.eosSquared = interpolationError;
projectionMeasure.mass = projectionMass; projectionMeasure.eosSquared = projectionError;
d::WriteError(errorQuadrature, "eos_nodal_interpolation", control->degree, element, fem.mesh->GetAttribute(element), 26, rule, interpolationMeasure);
d::WriteError(errorQuadrature, "eos_l2_q26", control->degree, element, fem.mesh->GetAttribute(element), 26, rule, projectionMeasure);
const double rhsNorm = d::projection::Norm(rhs), defectNorm = d::projection::Norm(defect);
elements << control->degree << ',' << element << ',' << fem.mesh->GetAttribute(element) << ",26," << rule.GetOrder() << ','
<< rule.GetNPoints() << ',' << size << ',' << std::quoted(control->collection->Name()) << ',' << volume << ',' << targetMass << ','
<< interpolationMass << ',' << projectionMass << ',' << interpolationError << ',' << projectionError << ','
<< rhsNorm << ',' << defectNorm << ',' << d::projection::Relative(defectNorm, rhsNorm) << ',' << constantError << ','
<< constantSampleError << ',' << interpolated.Min() << ',' << projected.Min() << ',' << projected.Max() << '\n';
for (int i = 0; i < size; ++i) {
const auto& point = nodalRule.IntPoint(i);
nodes << control->degree << ',' << element << ',' << i << ',' << dofs[i] << ',' << point.x << ',' << point.y << ',' << point.z
<< ',' << nodalH(i) << ',' << interpolated(i) << ',' << interpolated(i) << ',' << projected(i) << '\n';
}
}
}
elements.flush(); nodes.flush(); errorQuadrature.flush();
std::vector<d::ErrorField> errorFields;
for (const auto& control : controls) {
errorFields.push_back({"eos_nodal_interpolation", control->degree, control->interpolation.get()});
errorFields.push_back({"eos_l2_q26", control->degree, control->projection.get()});
}
const auto projectionsFinished = d::ObserverClock::now();
d::ObserveErrorRule(state, errorFields, 30, n, K, errorQuadrature);
const auto errorObservationFinished = d::ObserverClock::now();
d::Outputs observations(output);
for (const auto& control : controls) {
std::cout << "Frozen density order Q" << control->degree << " interpolation/projection observations" << std::endl;
const auto degree = std::to_string(control->degree);
d::projection::WriteVector(output / ("density_order_Q" + degree + "_eos_nodal_interpolation_coefficients.txt"), *control->interpolation);
d::projection::WriteVector(output / ("density_order_Q" + degree + "_eos_l2_q26_coefficients.txt"), *control->projection);
d::Observe(state, *control->interpolation, control->degree, "eos_nodal_interpolation", n, K, observations);
d::Observe(state, *control->projection, control->degree, "eos_l2_q26", n, K, observations);
}
checkpoint.Verify(state);
d::WriteScope(output / "density_order_scope.txt", true, sourceDegree, n, K);
d::WriteTimings(output / "density_order_timings.csv", observerStart, projectionsFinished,
errorObservationFinished, d::ObserverClock::now());
}
inline void RecordDensityOrderState(polytrope_validation::PhysicalState& state, double n, double K,
const std::filesystem::path& output) {
namespace d = density_order_detail;
const auto observerStart = d::ObserverClock::now();
const int degree = d::ValidateState(state, n, K);
std::filesystem::create_directories(output);
d::RootSelfChecks(output / "density_order_root_selfchecks.csv");
d::SourceCheckpoint checkpoint(state);
auto errorQuadrature = d::region::Open(output / "density_order_error_quadrature.csv");
d::ErrorHeader(errorQuadrature);
const std::array<d::ErrorField, 1> errorFields{{{"stored_density", degree, &state.density}}};
const auto setupFinished = d::ObserverClock::now();
for (int order : {26, 30}) d::ObserveErrorRule(state, errorFields, order, n, K, errorQuadrature);
const auto errorObservationFinished = d::ObserverClock::now();
d::Outputs observations(output);
d::Observe(state, state.density, degree, "stored_density", n, K, observations);
d::projection::WriteVector(output / "density_order_stored_density_coefficients.txt", state.density);
checkpoint.Verify(state);
d::WriteScope(output / "density_order_scope.txt", false, degree, n, K);
d::WriteTimings(output / "density_order_timings.csv", observerStart, setupFinished,
errorObservationFinished, d::ObserverClock::now());
}
} // namespace experiment::plateau_diagnostics

383
experiments/density_polynomial_regions.hpp Normal file
View File

@@ -0,0 +1,383 @@
#pragma once
// Diagnostic-only extension of the Q2 slice observer. A monotonic partition
// from derivative roots handles degrees three and four without sampling a sign
// indicator or relying on a particular closed-form cubic/quartic formula.
#include "negative_density_regions.hpp"
#include <string_view>
namespace experiment::plateau_diagnostics::density_order_detail {
namespace region = negative_density_regions_detail;
using Real = region::Real;
struct Polynomial final {
std::array<Real, 5> coefficients{};
int degree{0};
[[nodiscard]] Real Value(Real x) const {
Real value = coefficients[degree], correction = 0;
for (int j = degree - 1; j >= 0; --j) {
const Real product = value * x;
const Real productError = std::fma(value, x, -product);
const Real sum = product + coefficients[j];
const Real sumError = std::abs(product) >= std::abs(coefficients[j])
? (product - sum) + coefficients[j] : (coefficients[j] - sum) + product;
correction = std::fma(correction, x, productError + sumError);
value = sum;
}
// Compensated Horner retains the sign of shallow pockets without
// assuming long double has a wider mantissa on Apple arm64.
return value + correction;
}
[[nodiscard]] Real Scale() const {
Real result = 0;
for (int j = 0; j <= degree; ++j) result += std::abs(coefficients[j]);
return result;
}
[[nodiscard]] Polynomial Derivative() const {
Polynomial result;
result.degree = std::max(0, degree - 1);
for (int j = 1; j <= degree; ++j) result.coefficients[j - 1] = j * coefficients[j];
return result;
}
[[nodiscard]] Real Integral(Real left, Real right) const {
constexpr int binomial[5][5]{{1}, {1, 1}, {1, 2, 1}, {1, 3, 3, 1}, {1, 4, 6, 4, 1}};
const Real width = right - left, middle = left + width / 2;
Real mean = Value(middle);
for (int k = 2; k <= degree; k += 2) {
Real translated = 0;
for (int j = degree; j >= k; --j)
translated = std::fma(translated, middle, binomial[j][k] * coefficients[j]);
mean += translated * std::pow(width / 2, k) / (k + 1);
}
return width * mean;
}
};
inline Polynomial Interpolate(const std::array<double, 5>& values, int degree) {
if (degree < 2 || degree > 4) throw std::invalid_argument("Density slice degree must be two, three, or four.");
std::array<Real, 5> nodes{}, divided{}, basis{};
for (int i = 0; i <= degree; ++i) {
nodes[i] = static_cast<double>(i) / degree;
divided[i] = values[i];
if (!std::isfinite(values[i])) throw std::domain_error("Nonfinite density interpolation value.");
}
for (int width = 1; width <= degree; ++width)
for (int i = degree; i >= width; --i)
divided[i] = (divided[i] - divided[i - 1]) / (nodes[i] - nodes[i - width]);
Polynomial result;
result.degree = degree;
basis[0] = 1;
for (int k = 0; k <= degree; ++k) {
for (int j = 0; j <= k; ++j) result.coefficients[j] += divided[k] * basis[j];
if (k == degree) break;
for (int j = k + 1; j > 0; --j) basis[j] = basis[j - 1] - nodes[k] * basis[j];
basis[0] *= -nodes[k];
}
return result;
}
inline std::vector<Real> Roots(Polynomial polynomial) {
while (polynomial.degree > 0 && polynomial.coefficients[polynomial.degree] == 0) --polynomial.degree;
const Real scale = polynomial.Scale();
if (!std::isfinite(scale)) throw std::domain_error("Nonfinite density polynomial coefficients.");
if (polynomial.degree == 0 || scale == 0) return {};
int exponent = 0;
(void)std::frexp(scale, &exponent);
// A power-of-two rescaling preserves the represented polynomial. Dividing
// by an arbitrary norm can create a pocket at an exact repeated root.
for (int j = 0; j <= polynomial.degree; ++j)
polynomial.coefficients[j] = std::scalbn(polynomial.coefficients[j], -exponent);
if (polynomial.degree == 1) {
const Real root = -polynomial.coefficients[0] / polynomial.coefficients[1];
return std::isfinite(root) && root >= 0 && root <= 1 ? std::vector<Real>{root} : std::vector<Real>{};
}
std::vector<Real> partition{0, 1};
for (Real point : Roots(polynomial.Derivative())) if (point > 0 && point < 1) partition.push_back(point);
std::sort(partition.begin(), partition.end());
partition.erase(std::unique(partition.begin(), partition.end()), partition.end());
std::vector<Real> result;
for (Real point : partition) if (polynomial.Value(point) == 0) result.push_back(point);
for (std::size_t i = 1; i < partition.size(); ++i) {
Real left = partition[i - 1], right = partition[i];
Real fleft = polynomial.Value(left), fright = polynomial.Value(right);
if (fleft == 0 || fright == 0 || (fleft < 0) == (fright < 0)) continue;
// No value tolerance declares a critical point to be zero. Doing so
// could erase both roots of a shallow pocket next to that critical point.
for (int iteration = 0; iteration < 128; ++iteration) {
const Real middle = left + (right - left) / 2;
if (middle == left || middle == right) break;
const Real value = polynomial.Value(middle);
if (value == 0) { left = right = middle; break; }
if ((value < 0) == (fleft < 0)) { left = middle; fleft = value; }
else { right = middle; fright = value; }
}
result.push_back(std::abs(polynomial.Value(left)) <= std::abs(polynomial.Value(right)) ? left : right);
}
std::sort(result.begin(), result.end());
result.erase(std::unique(result.begin(), result.end()), result.end());
return result;
}
inline std::vector<region::Interval> Partition(const Polynomial& polynomial) {
std::vector<Real> points{0, 1};
for (Real root : Roots(polynomial)) if (root > 0 && root < 1) points.push_back(root);
std::sort(points.begin(), points.end());
points.erase(std::unique(points.begin(), points.end()), points.end());
std::vector<region::Interval> result;
for (std::size_t j = 1; j < points.size(); ++j) {
const Real left = points[j - 1], right = points[j];
const bool negative = polynomial.Value(left + (right - left) / 2) < 0;
for (Real fraction : {Real{0.25}, Real{0.75}}) {
const Real value = polynomial.Value(left + fraction * (right - left));
if ((negative && value > 0) || (!negative && value < 0))
throw std::runtime_error("Density polynomial root partition has inconsistent signs.");
}
result.push_back({left, right, negative});
}
return result;
}
inline void RootSelfChecks(const std::filesystem::path& path) {
auto out = region::Open(path);
out << "case,degree,negative_length,expected_negative_length,negative_mass,expected_negative_mass,passed\n";
const auto check = [&](const char* name, Polynomial polynomial, Real expectedLength, Real expectedMass) {
Real length = 0, mass = 0;
for (const auto& interval : Partition(polynomial)) if (interval.negative) {
length += interval.right - interval.left;
mass -= polynomial.Integral(interval.left, interval.right);
}
const bool passed = std::abs(length - expectedLength) < 2e-11L * std::max(Real{1e-8}, expectedLength) &&
std::abs(mass - expectedMass) < 2e-10L * std::max(Real{1e-8}, expectedMass);
out << name << ',' << polynomial.degree << ',' << length << ',' << expectedLength << ',' << mass << ',' << expectedMass << ',' << passed << '\n';
if (!passed) throw std::runtime_error(std::string("General density root self-check failed: ") + name);
};
check("positive_constant", {{{1}}, 0}, 0, 0);
check("negative_constant", {{{-2}}, 0}, 1, 2);
check("linear_crossing", {{{-0.25L, 1}}, 1}, 0.25L, 0.03125L);
check("quadratic_interior_pocket", {{{0.05L, -0.6L, 1.5L}}, 2}, std::sqrt(Real{6}) / 15, std::sqrt(Real{6}) / 2250);
check("cubic_interior_pocket", {{{0, 3.0L / 16, -1, 1}}, 3}, 0.5L, 1.0L / 96);
check("quartic_two_pockets", {{{105.0L / 4096, -11.0L / 32, 43.0L / 32, -2, 1}}, 4}, 0.5L, 19.0L / 15360);
check("quartic_double_roots_nonnegative", {{{9.0L / 256, -3.0L / 8, 11.0L / 8, -2, 1}}, 4}, 0, 0);
// Independently check the power-basis reconstruction, including its highest
// degree, rather than testing only roots of hand-entered coefficient arrays.
for (int degree : {2, 3, 4}) {
std::array<double, 5> values{};
for (int i = 0; i <= degree; ++i) values[i] = std::pow(static_cast<double>(i) / degree, degree);
const auto p = Interpolate(values, degree);
for (Real t : {Real{0.13}, Real{0.37}, Real{0.81}})
if (std::abs(p.Value(t) - std::pow(t, degree)) > 2e-14L)
throw std::runtime_error("Density polynomial interpolation self-check failed.");
out << "reconstruct_t_power_" << degree << ',' << degree << ",0,0,0,0,1\n";
}
}
struct Outputs final {
std::ofstream elements, summary, intervals, surface;
explicit Outputs(const std::filesystem::path& output)
: elements(region::Open(output / "density_order_elements.csv")),
summary(region::Open(output / "density_order_summary.csv")),
intervals(region::Open(output / "density_order_intervals.csv")),
surface(region::Open(output / "density_order_surface.csv")) {
constexpr auto extra = "target_mass,eos_squared_integral,signed_eos_integral,negative_mass_uncertainty_scale,"
"first_moment_x,first_moment_y,first_moment_z,inertia_z,minimum_enthalpy,negative_enthalpy_samples,"
"eos_defined_everywhere_sampled,resolved_negative_intervals,uncertain_negative_intervals,"
"resolved_witness_negative_volume,resolved_witness_negative_mass,uncertain_negative_volume,uncertain_negative_mass,";
elements << "control,degree,tangential_order,interval_order,element,attribute,axis,orientation,surface_face," << extra;
region::Header(elements);
summary << "control,degree,tangential_order,interval_order," << extra;
region::Header(summary);
intervals << "control,degree,element,axis,orientation,u,v,t_left,t_right,c0,c1,c2,c3,c4,negative_volume,negative_mass,"
"minimum_density,witness_h,reconstruction_tolerance,resolved_negative_witness\n";
surface << "control,degree,face,element,u,v,rho,h,rho_eos,error,eos_defined,reconstruction_tolerance,resolved_negative\n";
}
};
struct Extra final {
Real targetMass{0}, eosSquared{0}, signedEos{0}, uncertainty{0};
std::array<Real, 3> firstMoment{};
Real inertiaZ{0};
double minimumH{std::numeric_limits<double>::infinity()};
std::size_t negativeH{0};
std::size_t resolvedNegativeIntervals{0}, uncertainNegativeIntervals{0};
Real resolvedVolume{0}, resolvedMass{0}, uncertainVolume{0}, uncertainMass{0};
void Add(const Extra& other) {
targetMass += other.targetMass; eosSquared += other.eosSquared; signedEos += other.signedEos; uncertainty += other.uncertainty;
inertiaZ += other.inertiaZ; minimumH = std::min(minimumH, other.minimumH); negativeH += other.negativeH;
resolvedNegativeIntervals += other.resolvedNegativeIntervals; uncertainNegativeIntervals += other.uncertainNegativeIntervals;
resolvedVolume += other.resolvedVolume; resolvedMass += other.resolvedMass;
uncertainVolume += other.uncertainVolume; uncertainMass += other.uncertainMass;
for (int d = 0; d < 3; ++d) firstMoment[d] += other.firstMoment[d];
}
void Write(std::ostream& out) const {
const Real nan = std::numeric_limits<Real>::quiet_NaN();
out << (negativeH == 0 ? targetMass : nan) << ',' << (negativeH == 0 ? eosSquared : nan) << ','
<< (negativeH == 0 ? signedEos : nan) << ',' << uncertainty << ',';
for (Real value : firstMoment) out << value << ',';
out << inertiaZ << ',' << minimumH << ',' << negativeH << ',' << (negativeH == 0) << ','
<< resolvedNegativeIntervals << ',' << uncertainNegativeIntervals << ',' << resolvedVolume << ',' << resolvedMass << ','
<< uncertainVolume << ',' << uncertainMass << ',';
}
};
inline double EosTarget(double h, double n, double K) {
if (!std::isfinite(h)) throw std::domain_error("Nonfinite enthalpy in density order observer.");
const double rho = h > 0 ? std::pow(h / ((n + 1) * K), n) : 0;
if (!std::isfinite(rho)) throw std::domain_error("Nonfinite EOS target in density order observer.");
return rho;
}
inline void Observe(polytrope_validation::PhysicalState& state, const mfem::ParGridFunction& density,
int degree, std::string_view label, double n, double K, Outputs& output) {
const auto& fem = state.finiteElements;
const auto axes = region::ChooseAxes(state);
mean_field::mapping::GridFunctionMappingEvaluator mapper(state.domainMapper, state.displacement, *fem.compactificationCoordinate);
mean_field::mapping::MappingPointContext mapped;
const auto& intervalRule = mfem::IntRules.Get(mfem::Geometry::SEGMENT, 14);
for (int order : {14, 22, 30}) {
const auto& tangentRule = mfem::IntRules.Get(mfem::Geometry::SEGMENT, order);
region::Measure global;
Extra globalExtra;
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);
region::Measure local;
Extra localExtra;
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;
std::array<double, 5> values{};
for (int i = 0; i <= degree; ++i) values[i] = density.GetValue(element, region::Point(axis, u, v, static_cast<double>(i) / degree));
const auto polynomial = Interpolate(values, degree);
const Real tolerance = 2048 * std::numeric_limits<double>::epsilon() * polynomial.Scale();
const auto verify = [&](Real t, double rho) {
const Real error = std::abs(polynomial.Value(t) - rho);
local.maximumReconstructionError = std::max(local.maximumReconstructionError, static_cast<double>(error));
local.maximumReconstructionTolerance = std::max(local.maximumReconstructionTolerance, static_cast<double>(tolerance));
if (!std::isfinite(rho) || error > tolerance) throw std::runtime_error("Density order slice reconstruction failed.");
};
for (Real t : {Real{0.07}, Real{0.19}, Real{0.41}, Real{0.63}, Real{0.87}, Real{0.97}})
verify(t, density.GetValue(element, region::Point(axis, u, v, t)));
++local.lines;
const auto critical = Roots(polynomial.Derivative());
const auto evaluate = [&](Real t, region::Measure& piece, bool negative) {
const auto point = region::Point(axis, u, v, t);
if (mapper.EvaluatePoint(*transformation, point, mapped) != mean_field::mapping::MappingStatus::valid || mapped.compactified)
throw std::domain_error("Invalid density order mapping.");
const double rho = density.GetValue(element, point), h = state.enthalpy.GetValue(element, point);
verify(t, rho);
if (!std::isfinite(h)) throw std::domain_error("Nonfinite density order enthalpy.");
if (rho < piece.minimumRho) {
piece.minimumRho = rho; piece.witnessH = h; piece.witnessElement = element;
piece.witnessT = static_cast<double>(t); piece.witnessU = u; piece.witnessV = v;
for (int d = 0; d < 3; ++d) piece.witnessPosition[d] = mapped.physical_position(d);
}
if (negative) {
const double radius = mapped.physical_position.Norml2();
piece.negativeRadiusMin = std::min(piece.negativeRadiusMin, radius);
piece.negativeRadiusMax = std::max(piece.negativeRadiusMax, radius);
}
transformation->SetIntPoint(&point);
const double reference = transformation->Jacobian().Det();
const Real determinant = static_cast<Real>(reference) * mapped.mapping_determinant;
if (!(reference > 0) || !std::isfinite(reference) || !(mapped.mapping_determinant > 0) ||
!std::isfinite(mapped.mapping_determinant) || !std::isfinite(determinant))
throw std::domain_error("Invalid density order physical measure.");
return std::pair{determinant, h};
};
for (const auto& part : Partition(polynomial)) {
region::Measure piece;
(void)evaluate(part.left, piece, part.negative);
(void)evaluate(part.right, piece, part.negative);
for (Real point : critical) if (point > part.left && point < part.right) (void)evaluate(point, piece, part.negative);
for (int iq = 0; iq < intervalRule.GetNPoints(); ++iq) {
const auto& qp = intervalRule.IntPoint(iq);
const Real t = part.left + (part.right - part.left) * qp.x;
const auto [determinant, h] = evaluate(t, piece, part.negative);
const Real weight = static_cast<Real>(up.weight) * vp.weight * qp.weight * (part.right - part.left) * determinant;
const Real rho = polynomial.Value(t), target = EosTarget(h, n, K);
if (!(weight > 0) || !std::isfinite(weight) || !std::isfinite(rho) ||
(part.negative && rho > 0) || (!part.negative && rho < 0))
throw std::runtime_error("Density order interval quadrature has an inconsistent sign.");
piece.volume += weight; piece.mass += weight * rho; ++piece.points;
localExtra.targetMass += weight * target;
localExtra.eosSquared += weight * (rho - target) * (rho - target);
localExtra.signedEos += weight * (rho - target);
localExtra.uncertainty += weight * tolerance;
localExtra.minimumH = std::min(localExtra.minimumH, h);
localExtra.negativeH += h < 0;
for (int d = 0; d < 3; ++d) localExtra.firstMoment[d] += weight * rho * mapped.physical_position(d);
localExtra.inertiaZ += weight * rho * (static_cast<Real>(mapped.physical_position(0)) * mapped.physical_position(0) +
static_cast<Real>(mapped.physical_position(1)) * mapped.physical_position(1));
}
if (part.negative) {
piece.negativeIntervals = 1; piece.negativeVolume = piece.volume; piece.negativeMass = -piece.mass;
piece.negativeTMin = part.left; piece.negativeTMax = part.right;
const bool resolved = piece.minimumRho < -tolerance;
if (resolved) {
++localExtra.resolvedNegativeIntervals;
localExtra.resolvedVolume += piece.negativeVolume; localExtra.resolvedMass += piece.negativeMass;
} else {
++localExtra.uncertainNegativeIntervals;
localExtra.uncertainVolume += piece.negativeVolume; localExtra.uncertainMass += piece.negativeMass;
}
if (order == 30) {
output.intervals << label << ',' << degree << ',' << element << ',' << axis.coordinate << ',' << axis.orientation
<< ',' << u << ',' << v << ',' << part.left << ',' << part.right;
for (Real coefficient : polynomial.coefficients) output.intervals << ',' << coefficient;
output.intervals << ',' << piece.negativeVolume << ',' << piece.negativeMass << ',' << piece.minimumRho << ','
<< piece.witnessH << ',' << tolerance << ',' << resolved << '\n';
}
}
local.Add(piece);
}
}
if (!(local.volume > 0) || !std::isfinite(local.mass) || !std::isfinite(local.negativeMass))
throw std::runtime_error("Nonfinite density order region integral.");
output.elements << label << ',' << degree << ',' << order << ",14," << element << ',' << fem.mesh->GetAttribute(element) << ','
<< axis.coordinate << ',' << axis.orientation << ',' << axis.face << ',';
localExtra.Write(output.elements); region::Write(output.elements, local);
global.Add(local); globalExtra.Add(localExtra);
}
if (!(global.volume > 0)) throw std::runtime_error("Density order observer found no stellar volume.");
output.summary << label << ',' << degree << ',' << order << ",14,";
globalExtra.Write(output.summary); region::Write(output.summary, global);
output.summary.flush(); output.elements.flush(); output.intervals.flush();
}
// These are one-sided field traces on every stellar/nonstellar face, not
// surface-area samples. They expose extrapolation defects missed by nodes.
for (int face = 0; face < fem.mesh->GetNumFaces(); ++face) {
int left = -1, right = -1;
fem.mesh->GetFaceElements(face, &left, &right);
const bool leftStar = left >= 0 && state.isStellar(left), rightStar = right >= 0 && state.isStellar(right);
if (leftStar == rightStar) continue;
const int element = leftStar ? left : right;
auto* transformation = fem.mesh->GetFaceElementTransformations(face);
mfem::Array<int> dofs;
auto* dofTransform = density.ParFESpace()->GetElementDofs(element, dofs);
mfem::Vector coefficients, shape(dofs.Size());
density.GetSubVector(dofs, coefficients);
if (dofTransform) dofTransform->InvTransformPrimal(coefficients);
const auto* fe = density.ParFESpace()->GetFE(element);
for (double u : {0., .25, .5, .75, 1.}) for (double v : {0., .25, .5, .75, 1.}) {
mfem::IntegrationPoint facePoint, point;
facePoint.Set2(u, v);
if (leftStar) transformation->Loc1.Transform(facePoint, point); else transformation->Loc2.Transform(facePoint, point);
const double rho = density.GetValue(element, point), h = state.enthalpy.GetValue(element, point);
if (!std::isfinite(h)) throw std::domain_error("Nonfinite density order surface enthalpy.");
const double target = h >= 0 ? EosTarget(h, n, K) : std::numeric_limits<double>::quiet_NaN();
if (!std::isfinite(rho)) throw std::domain_error("Nonfinite density order surface trace.");
fe->CalcShape(point, shape);
Real magnitude = 0;
for (int i = 0; i < shape.Size(); ++i) magnitude += std::abs(static_cast<Real>(shape(i)) * coefficients(i));
const Real tolerance = 2048 * std::numeric_limits<double>::epsilon() * magnitude;
output.surface << label << ',' << degree << ',' << face << ',' << element << ',' << u << ',' << v << ','
<< rho << ',' << h << ',' << target << ',' << rho - target << ',' << (h >= 0) << ',' << tolerance << ',' << (rho < -tolerance) << '\n';
}
}
output.surface.flush();
}
} // namespace experiment::plateau_diagnostics::density_order_detail

31
experiments/plateau_diagnostic.cpp
View File

@@ -35,6 +35,7 @@ import mean_field;
#include "spherical_harmonic_surface_diagnostics.hpp"
#include "negative_density_diagnostics.hpp"
#include "negative_density_regions.hpp"
#include "density_order_projection.hpp"
using namespace mean_field;
namespace {
@@ -49,6 +50,7 @@ struct Options {
bool reducedSurfaceProbe = false, reducedSurfaceAlgebraCheck = false;
bool densityDiagnostic = false;
bool densityRegions = false;
bool densityOrderProjection = false, densityOrderFinal = false;
int surfaceDegree = 4, surfaceMaximumOrder = -1;
std::string surfacePrescription = "nodal";
bool coupledTrialProbe = false, coupledSeedControl = false;
@@ -135,6 +137,12 @@ Options Parse(int argc, char **argv) {
o.densityDiagnostic = true;
o.densityRegions = true;
}
else if (a == "--density-order-projection") {
o.densityDiagnostic = true;
o.densityOrderProjection = true;
}
else if (a == "--density-order-final")
o.densityOrderFinal = true;
else if (a == "--density-snapshot")
o.densitySnapshot = v();
else if (a == "--surface-prescription")
@@ -209,6 +217,8 @@ Options Parse(int argc, char **argv) {
throw std::invalid_argument("Surface prescriptions are nodal or harmonics");
if (!o.densitySnapshot.empty() && !o.densityDiagnostic)
throw std::invalid_argument("--density-snapshot requires --density-diagnostic");
if ((o.densityOrderProjection && o.densityRegions) || (o.densityDiagnostic && o.densityOrderFinal))
throw std::invalid_argument("Choose one frozen density observation or the coupled final observation");
if (o.densityDiagnostic && (o.maxNewton != 0 || o.strategy != "baseline" || o.probes ||
o.pressureProbe || o.shapeCurvature || o.surfaceTrace || o.strategyProbe || o.positivityProbe ||
o.positivityCertificate || o.preconditionerProbe || o.nonlinearStepProbe || o.reducedSurfaceProbe ||
@@ -687,6 +697,13 @@ int RunContext(Context &context, const Options &o, const Scales &scales) {
<< "\ncontext_gravity_schur_relaxation=" << (contextUsesGravitySchur ? o.schurRelaxation : 0.0)
<< "\nschur_correction_sweeps=" << o.schurCorrectionSweeps << "\nschur_relaxation=" << o.schurRelaxation
<< "\ncompiler=" << __VERSION__ << "\ncompiled=" << __DATE__ << ' ' << __TIME__ << '\n';
metadata << "density_order=" << field::Density::scalarOrder
<< "\nenthalpy_order=" << field::Enthalpy::scalarOrder
<< "\ngravity_potential_order=" << field::Gravity::potentialOrder
<< "\ngravity_flux_order=" << field::Gravity::fluxOrder
<< "\ndisplacement_order=" << field::Displacement::Vector::familyOrder
<< "\ndensity_order_projection=" << o.densityOrderProjection
<< "\ndensity_order_final=" << o.densityOrderFinal << '\n';
if (o.densityDiagnostic) {
metadata << "density_diagnostic=1\ndensity_region_integral=" << o.densityRegions << "\ndensity_source="
<< (o.densitySnapshot.empty() ? "lane_emden_seed" : o.densitySnapshot.string()) << '\n';
@@ -727,7 +744,9 @@ int RunContext(Context &context, const Options &o, const Scales &scales) {
adapter.scatter(closureResidual, closureFull);
SaveVector(o.output / "production_closure_residual_full.txt", closureFull);
SaveVector(o.output / "observed_density_coefficients.txt", physical.density);
if (o.densityRegions)
if (o.densityOrderProjection)
experiment::plateau_diagnostics::RecordDensityOrderProjection(physical, o.n, scales.K, o.output);
else if (o.densityRegions)
experiment::plateau_diagnostics::RecordDensityNegativeRegions(physical, o.n, scales.K, o.output);
else
experiment::plateau_diagnostics::RecordNegativeDensity(physical, o.n, scales.K, o.output);
@@ -1215,6 +1234,16 @@ int RunContext(Context &context, const Options &o, const Scales &scales) {
solver::detail::StellarEquilibriumContextDiagnostics::WithState(context, [&](auto &s, const fem::FEM &f) {
SaveVector(o.output / "accepted_state.txt", s.AcceptedPhysicalState());
physicalMetrics(s, f, "final");
if (o.densityOrderFinal) {
experiment::polytrope_validation::PhysicalState physical(s, f);
experiment::plateau_diagnostics::RecordDensityOrderState(physical, o.n, scales.K, o.output);
auto rules = File(o.output / "final_closure_rules.csv");
rules << "element,quadrature_order,quadrature_points\n";
s.Problem().GetPhysicalOperator().GetBarotropicClosureOperator().VisitMappedGeometryRules(
[&](int element, const mfem::IntegrationRule &rule) {
rules << element << ',' << rule.GetOrder() << ',' << rule.GetNPoints() << '\n';
});
}
const auto &schur = s.physicalInverse->GetGroupedPreconditioner().GetFactorization().GetSchurComplement();
auto sh = File(o.output / "border_schur.csv");
for (int i = 0; i < schur.Height(); ++i) {