diff --git a/experiments/density_order_coupled_cases_2026-09-19.json b/experiments/density_order_coupled_cases_2026-09-19.json new file mode 100644 index 0000000..93fbd44 --- /dev/null +++ b/experiments/density_order_coupled_cases_2026-09-19.json @@ -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" + ] + } +] diff --git a/experiments/density_order_frozen_cases_2026-09-19.json b/experiments/density_order_frozen_cases_2026-09-19.json new file mode 100644 index 0000000..a6c97f0 --- /dev/null +++ b/experiments/density_order_frozen_cases_2026-09-19.json @@ -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" + ] + } +] diff --git a/experiments/density_order_projection.hpp b/experiments/density_order_projection.hpp new file mode 100644 index 0000000..c54e325 --- /dev/null +++ b/experiments/density_order_projection.hpp @@ -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 +#include +#include +#include + +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(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(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 collection; + std::unique_ptr space; + std::unique_ptr interpolation, projection, constant; + mfem::DenseMatrix shapes; + DegreeFields(mfem::ParMesh& mesh, int order, const mfem::IntegrationRule& rule) + : degree(order), collection(std::make_unique(order, 3, mfem::BasisType::GaussLegendre, mfem::FiniteElement::VALUE)), + space(std::make_unique(&mesh, collection.get())), + interpolation(std::make_unique(space.get())), + projection(std::make_unique(space.get())), + constant(std::make_unique(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::infinity()}; + double minimumH{std::numeric_limits::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::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 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> shapeCache; + std::vector shapes(fields.size()); + std::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(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 dofs; + auto* dofTransform = density.ParFESpace()->GetElementDofs(element, dofs); + density.GetSubVector(dofs, coefficients[field]); + if (dofTransform) dofTransform->InvTransformPrimal(coefficients[field]); + } + std::vector 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(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(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> controls; + for (int degree : {2, 3, 4}) controls.push_back(std::make_unique(*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::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(weights(q)) * targets(q); + } + for (auto& control : controls) { + const auto* fe = control->space->GetFE(element); + const int size = fe->GetDof(); + mfem::Array 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(shape) * interpolated(i); + projectedValue += static_cast(shape) * projected(i); + constantValue += static_cast(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(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 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 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 diff --git a/experiments/density_polynomial_regions.hpp b/experiments/density_polynomial_regions.hpp new file mode 100644 index 0000000..f27ad28 --- /dev/null +++ b/experiments/density_polynomial_regions.hpp @@ -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 + +namespace experiment::plateau_diagnostics::density_order_detail { +namespace region = negative_density_regions_detail; +using Real = region::Real; + +struct Polynomial final { + std::array 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& values, int degree) { + if (degree < 2 || degree > 4) throw std::invalid_argument("Density slice degree must be two, three, or four."); + std::array nodes{}, divided{}, basis{}; + for (int i = 0; i <= degree; ++i) { + nodes[i] = static_cast(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 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{root} : std::vector{}; + } + std::vector 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 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 Partition(const Polynomial& polynomial) { + std::vector 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 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 values{}; + for (int i = 0; i <= degree; ++i) values[i] = std::pow(static_cast(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 firstMoment{}; + Real inertiaZ{0}; + double minimumH{std::numeric_limits::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::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(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 values{}; + for (int i = 0; i <= degree; ++i) values[i] = density.GetValue(element, region::Point(axis, u, v, static_cast(i) / degree)); + const auto polynomial = Interpolate(values, degree); + const Real tolerance = 2048 * std::numeric_limits::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(error)); + local.maximumReconstructionTolerance = std::max(local.maximumReconstructionTolerance, static_cast(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(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(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(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(mapped.physical_position(0)) * mapped.physical_position(0) + + static_cast(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 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::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(shape(i)) * coefficients(i)); + const Real tolerance = 2048 * std::numeric_limits::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 diff --git a/experiments/plateau_diagnostic.cpp b/experiments/plateau_diagnostic.cpp index 0a1a256..8d86bcd 100644 --- a/experiments/plateau_diagnostic.cpp +++ b/experiments/plateau_diagnostic.cpp @@ -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) {