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

Add physical-plane FEM model slice diagnostics and plots

Browse Source
This commit is contained in:
Emily Boudreaux
2026-09-19 13:46:11 -04:00
parent 7793fd9111
commit 9bac57a1df
7 changed files with 888 additions and 1 deletions
Download Patch File Download Diff File Expand all files Collapse all files

34
experiments/plateau_diagnostic.cpp
View File

@@ -38,6 +38,7 @@ import mean_field;
#include "density_order_projection.hpp"
#include "continuation_snapshot.hpp"
#include "harmonic_linear_trials.hpp"
#include "model_slice_diagnostics.hpp"
using namespace mean_field;
namespace {
@@ -53,6 +54,8 @@ struct Options {
bool densityDiagnostic = false;
bool densityRegions = false;
bool densityOrderProjection = false, densityOrderFinal = false;
bool modelSlices = false, modelSlicesOnly = false;
experiment::model_slices::Options slices;
int surfaceDegree = 4, surfaceMaximumOrder = -1;
std::string surfacePrescription = "nodal";
bool coupledTrialProbe = false, coupledSeedControl = false;
@@ -154,6 +157,17 @@ Options Parse(int argc, char **argv) {
o.allowUnconvergedRestart = true;
else if (a == "--benchmark-trials")
o.benchmarkTrials = true;
else if (a == "--model-slices")
o.modelSlices = true;
else if (a == "--model-slices-only") {
o.modelSlices = o.modelSlicesOnly = o.densityDiagnostic = true;
}
else if (a == "--slice-plane")
o.slices.plane = v();
else if (a == "--slice-offset")
o.slices.offset = std::stod(v());
else if (a == "--slice-subdivisions")
o.slices.subdivisions = std::stoi(v());
else if (a == "--density-snapshot")
o.densitySnapshot = v();
else if (a == "--surface-prescription")
@@ -231,6 +245,14 @@ Options Parse(int argc, char **argv) {
throw std::invalid_argument("Snapshot restart requires native harmonic coordinates and no legacy replay");
if (o.benchmarkTrials && !o.preconditionerProbe)
throw std::invalid_argument("--benchmark-trials requires --preconditioner-probe");
o.slices.Validate();
if (o.modelSlicesOnly && (o.densityOrderProjection || o.densityRegions || o.densityOrderFinal))
throw std::invalid_argument("--model-slices-only cannot request another density observation");
if (o.modelSlices && (o.probes || o.pressureProbe || o.shapeCurvature || o.surfaceTrace || o.strategyProbe ||
o.positivityProbe || o.positivityCertificate || o.preconditionerProbe || o.nonlinearStepProbe ||
o.reducedSurfaceProbe || o.reducedSurfaceAlgebraCheck || o.weakShapeProbe || o.coupledTrialProbe ||
o.dampedCoupledProbe || !o.probeState.empty() || !o.probeDirection.empty() || !o.replayCandidate.empty()))
throw std::invalid_argument("Model slices require live initial/final diagnostics or density-snapshot replay");
if (o.surfacePrescription != "nodal" && o.surfacePrescription != "harmonics")
throw std::invalid_argument("Surface prescriptions are nodal or harmonics");
if (!o.densitySnapshot.empty() && !o.densityDiagnostic)
@@ -755,6 +777,9 @@ int RunContext(Context &context, const Options &o, const Scales &scales) {
<< "\ndisplacement_order=" << field::Displacement::Vector::familyOrder
<< "\ndensity_order_projection=" << o.densityOrderProjection
<< "\ndensity_order_final=" << o.densityOrderFinal << '\n';
metadata << "model_slices=" << o.modelSlices << "\nmodel_slices_only=" << o.modelSlicesOnly
<< "\nslice_plane=" << o.slices.plane << "\nslice_offset=" << o.slices.offset
<< "\nslice_subdivisions=" << o.slices.subdivisions << '\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';
@@ -795,7 +820,12 @@ 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.densityOrderProjection)
if (o.modelSlices)
experiment::model_slices::Record(physical, o.slices, o.output, "frozen");
if (o.modelSlicesOnly) {
// Replay keeps its ordinary physical metrics and provenance,
// but a plot request need not integrate all negative regions.
} else 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);
@@ -997,6 +1027,8 @@ int RunContext(Context &context, const Options &o, const Scales &scales) {
std::cout << label << " physical diagnostics" << std::endl;
experiment::polytrope_validation::PhysicalState physical(s, f);
auto out = File(o.output / (label + "_physical.csv"));
if (o.modelSlices)
experiment::model_slices::Record(physical, o.slices, o.output, label);
out << "quadrature,metric,value\n";
if (o.surfacePrescription == "harmonics") {
deformation::SphericalHarmonicSurfaceOptions surfaceOptions{.maximumDegree = o.surfaceDegree};