Emily Boudreaux
75cc638739
This commit uses global pre allocated work space to dramatically reduce memory usage and allocation time
446 lines
28 KiB
C++
446 lines
28 KiB
C++
#include <algorithm>
|
|
#include <array>
|
|
#include <chrono>
|
|
#include <cmath>
|
|
#include <cstdlib>
|
|
#include <filesystem>
|
|
#include <fstream>
|
|
#include <iomanip>
|
|
#include <iostream>
|
|
#include <limits>
|
|
#include <map>
|
|
#include <numbers>
|
|
#include <stdexcept>
|
|
#include <string>
|
|
#include <vector>
|
|
#include <mfem.hpp>
|
|
|
|
import mean_field;
|
|
|
|
#include "polytrope_analytic_self_checks.hpp"
|
|
#include "polytrope_validation_measurements.hpp"
|
|
#include "polytrope_radial_profiles.hpp"
|
|
|
|
namespace {
|
|
using namespace mean_field;
|
|
using namespace experiment::polytrope_validation;
|
|
using Clock = std::chrono::steady_clock;
|
|
|
|
struct Options final {
|
|
std::string mesh{"sandbox.smesh"};
|
|
std::filesystem::path output{"polytrope_validation_results"};
|
|
std::filesystem::path replay;
|
|
std::string mode{"solve"};
|
|
double absoluteTolerance{1.0e-8};
|
|
double relativeTolerance{1.0e-8};
|
|
double linearTolerance{0.03};
|
|
int maximumNewtonIterations{8};
|
|
int maximumLinearIterations{80};
|
|
int quadratureOrder{14};
|
|
int checkQuadratureOrder{18};
|
|
bool profiles{true};
|
|
RadialProfileOptions radial;
|
|
};
|
|
|
|
Options Parse(const int argc, char **argv) {
|
|
Options options;
|
|
for (int index = 1; index < argc; ++index) {
|
|
const std::string argument = argv[index];
|
|
auto value = [&]() -> std::string {
|
|
if (++index >= argc) throw std::invalid_argument("Missing value after " + argument);
|
|
return argv[index];
|
|
};
|
|
if (argument == "--mesh") options.mesh = value();
|
|
else if (argument == "--output") options.output = value();
|
|
else if (argument == "--self-check") options.mode = "self-check";
|
|
else if (argument == "--analytic-mesh") options.mode = "analytic-mesh";
|
|
else if (argument == "--solve") options.mode = "solve";
|
|
else if (argument == "--replay") {
|
|
options.mode = "replay";
|
|
options.replay = value();
|
|
options.mesh = (options.replay / "input.smesh").string();
|
|
}
|
|
else if (argument == "--absolute-tolerance") options.absoluteTolerance = std::stod(value());
|
|
else if (argument == "--relative-tolerance") options.relativeTolerance = std::stod(value());
|
|
else if (argument == "--linear-tolerance") options.linearTolerance = std::stod(value());
|
|
else if (argument == "--max-newton") options.maximumNewtonIterations = std::stoi(value());
|
|
else if (argument == "--max-linear-iterations") options.maximumLinearIterations = std::stoi(value());
|
|
else if (argument == "--quadrature-order") options.quadratureOrder = std::stoi(value());
|
|
else if (argument == "--check-quadrature-order") options.checkQuadratureOrder = std::stoi(value());
|
|
else if (argument == "--skip-profiles") options.profiles = false;
|
|
else if (argument == "--mu-points") options.radial.muPointCount = std::stoi(value());
|
|
else if (argument == "--phi-points") options.radial.azimuthPointCount = std::stoi(value());
|
|
else if (argument == "--exterior-shells") options.radial.exteriorShellCount = std::stoi(value());
|
|
else if (argument == "--help") {
|
|
std::cout << "polytrope_validation_experiment [--self-check | --analytic-mesh | --solve]\n"
|
|
" [--mesh sandbox.smesh] [--output NEW_DIRECTORY] [--skip-profiles]\n"
|
|
" [--absolute-tolerance 1e-8] [--relative-tolerance 1e-8]\n"
|
|
" [--linear-tolerance 0.03] [--max-newton 8] [--max-linear-iterations 80]\n"
|
|
" [--quadrature-order 14] [--check-quadrature-order 18]\n"
|
|
" [--replay SOLVER_OUTPUT_DIRECTORY] (saved nonrotating fields; no Newton context)\n"
|
|
" [--mu-points 6] [--phi-points 12] [--exterior-shells 8]\n"
|
|
"Single MPI rank. Default: nonrotating n=1 production solve.\n"
|
|
"Exit codes: 0 all requested checks pass; 1 nonlinear failure; 2 execution error;\n"
|
|
"3 verification failure. Existing output directories are never overwritten.\n";
|
|
options.mode = "help";
|
|
return options;
|
|
} else throw std::invalid_argument("Unknown option: " + argument);
|
|
}
|
|
if (!std::isfinite(options.absoluteTolerance) || options.absoluteTolerance < 0.0 ||
|
|
!std::isfinite(options.relativeTolerance) || options.relativeTolerance < 0.0 ||
|
|
!(options.absoluteTolerance > 0.0 || options.relativeTolerance > 0.0) ||
|
|
!(options.linearTolerance > 0.0 && options.linearTolerance < 1.0) ||
|
|
options.maximumNewtonIterations < 1 || options.maximumLinearIterations < 1 ||
|
|
options.quadratureOrder < 2 || options.checkQuadratureOrder <= options.quadratureOrder ||
|
|
options.radial.muPointCount < 2 || options.radial.azimuthPointCount < 4 || options.radial.exteriorShellCount < 1) {
|
|
throw std::invalid_argument("Invalid tolerance, iteration limit, or quadrature orders.");
|
|
}
|
|
if (options.mode == "replay") options.mesh = (options.replay / "input.smesh").string();
|
|
return options;
|
|
}
|
|
|
|
std::ofstream File(const std::filesystem::path &path) {
|
|
std::ofstream stream(path);
|
|
stream.exceptions(std::ios::failbit | std::ios::badbit);
|
|
stream << std::setprecision(17);
|
|
return stream;
|
|
}
|
|
|
|
bool SelfChecks(const Options &options) {
|
|
const auto report = RunAnalyticSelfChecks();
|
|
auto stream = File(options.output / "analytic_self_checks.csv");
|
|
stream << "check,observed,expected,scale,absolute_error,scaled_error,tolerance,passed\n";
|
|
for (const auto &check : report.checks) {
|
|
stream << check.name << ',' << check.observed << ',' << check.expected << ',' << check.scale << ','
|
|
<< check.AbsoluteError() << ',' << check.ScaledError() << ',' << check.tolerance << ',' << check.passed << '\n';
|
|
if (!check.passed) std::cerr << "Analytic check failed: " << check.name << " scaled error=" << check.ScaledError() << '\n';
|
|
}
|
|
std::cout << "Independent analytic checks: " << report.checks.size() << ", passed=" << report.Passed() << std::endl;
|
|
return report.Passed();
|
|
}
|
|
|
|
void WriteMetrics(const std::filesystem::path &path, const Measurements &measurements) {
|
|
auto stream = File(path);
|
|
stream << "metric,value\n";
|
|
for (const auto &[name, value] : measurements) stream << name << ',' << value << '\n';
|
|
}
|
|
|
|
std::map<std::string, std::string> ReadMetadata(const std::filesystem::path &directory) {
|
|
std::ifstream stream(directory / "metadata.txt");
|
|
if (!stream) throw std::runtime_error("Cannot read replay metadata.");
|
|
std::map<std::string, std::string> result;
|
|
for (std::string line; std::getline(stream, line);) {
|
|
const auto separator = line.find('=');
|
|
if (separator != std::string::npos) result[line.substr(0, separator)] = line.substr(separator+1);
|
|
}
|
|
return result;
|
|
}
|
|
|
|
Measurements ReadSavedSolverMetrics(const std::filesystem::path &directory) {
|
|
std::ifstream stream(directory / "physical_metrics.csv");
|
|
if (!stream) throw std::runtime_error("Replay needs completed physical_metrics.csv to preserve solver/border diagnostics.");
|
|
Measurements result;
|
|
std::string line;
|
|
std::getline(stream, line);
|
|
while (std::getline(stream, line)) {
|
|
const auto separator = line.find(',');
|
|
if (separator == std::string::npos) throw std::runtime_error("Malformed saved physical metrics.");
|
|
const auto name = line.substr(0, separator);
|
|
if (name == "normalized_bordered_residual" || name == "normalized_unbordered_residual" ||
|
|
name == "normalized_central_border_action" || name == "central_border" ||
|
|
name == "bernoulli_constant" || name == "angular_velocity_norm") {
|
|
result[name] = std::stod(line.substr(separator+1));
|
|
}
|
|
}
|
|
if (result.size() != 6 || result.at("angular_velocity_norm") != 0.0) {
|
|
throw std::runtime_error("Replay currently requires complete saved diagnostics and exactly zero rotation.");
|
|
}
|
|
return result;
|
|
}
|
|
|
|
template <typename SampleState>
|
|
bool Measure(SampleState &state, const N1Reference &reference, const Options &options,
|
|
Measurements additional = {}) {
|
|
std::cout << "Physical volume integration, order=" << options.quadratureOrder << std::endl;
|
|
const auto base = MeasureVolumes(state, reference, options.quadratureOrder);
|
|
WriteMetrics(options.output / "volume_metrics_base.csv", base);
|
|
std::cout << "Independent higher-order integration, order=" << options.checkQuadratureOrder << std::endl;
|
|
auto metrics = MeasureVolumes(state, reference, options.checkQuadratureOrder);
|
|
auto quadrature = File(options.output / "quadrature_comparison.csv");
|
|
quadrature << "metric,base,check,absolute_difference,relative_difference\n";
|
|
for (const auto &[name, high] : metrics) {
|
|
const double low = base.at(name);
|
|
quadrature << name << ',' << low << ',' << high << ',' << std::abs(high-low) << ','
|
|
<< std::abs(high-low) / std::max(std::abs(high), 1.0e-300) << '\n';
|
|
}
|
|
metrics["quadrature_virial_absolute_change"] = std::abs(metrics.at("virial_signed") - base.at("virial_signed"));
|
|
metrics["quadrature_binding_relative_change"] = std::abs(metrics.at("binding_energy") - base.at("binding_energy")) / std::abs(reference.BindingEnergy());
|
|
std::cout << "Stellar surface and all-element corner sampling" << std::endl;
|
|
metrics.merge(MeasureSurfaceAndCorners(state, reference, options.checkQuadratureOrder));
|
|
metrics.merge(additional);
|
|
if (options.profiles) {
|
|
std::cout << "Physical radial projection (stellar interior and finite exterior)" << std::endl;
|
|
const auto profiles = WriteRadialProfiles(state, reference, options.output, options.radial);
|
|
metrics["profile_requested_points"] = profiles.requestedPoints;
|
|
metrics["profile_missing_points"] = profiles.requestedPoints - profiles.locatedPoints;
|
|
metrics["profile_material_points"] = profiles.materialPoints;
|
|
metrics["profile_locator_element_attempts"] = profiles.locatorElementAttempts;
|
|
metrics["profile_maximum_location_error"] = profiles.maximumLocationError;
|
|
metrics["profile_angular_moment_error"] = profiles.angularMomentError;
|
|
metrics["profile_maximum_scaled_error"] = profiles.maximumScaledFieldError;
|
|
metrics["profile_maximum_angular_rms_scaled"] = profiles.maximumAngularRmsScaled;
|
|
}
|
|
metrics["negative_density_maximum_scaled"] = std::max(0.0, -metrics.at("minimum_density")) / reference.CentralDensity();
|
|
metrics["negative_enthalpy_maximum_scaled"] = std::max(0.0, -metrics.at("minimum_enthalpy")) / reference.CentralEnthalpy();
|
|
WriteMetrics(options.output / "physical_metrics.csv", metrics);
|
|
|
|
// Initial screening budgets, declared before running the solver. A pass
|
|
// is not a mesh-convergence certificate. The complete errors are saved.
|
|
std::map<std::string, double> budgets{
|
|
{"mass_relative_error", 1.0e-4}, {"volume_radius_relative_error", 1.0e-4},
|
|
{"surface_radius_relative_rms_error", 1.0e-4},
|
|
{"density_relative_l2_error", 1.0e-4}, {"enthalpy_relative_l2_error", 1.0e-4},
|
|
{"potential_relative_l2_error", 1.0e-4}, {"gravity_gradient_relative_l2_error", 1.0e-4},
|
|
{"binding_relative_error", 1.0e-4}, {"pressure_integral_relative_error", 1.0e-4},
|
|
{"moment_of_inertia_relative_error", 1.0e-4}, {"virial_error", 1.0e-6},
|
|
{"force_virial_error", 1.0e-6}, {"gravity_energy_consistency", 1.0e-6},
|
|
{"eos_enthalpy_scaled_rms", 1.0e-6}, {"bernoulli_scaled_rms_variation", 1.0e-4},
|
|
{"quadrature_virial_absolute_change", 1.0e-8}, {"quadrature_binding_relative_change", 1.0e-8},
|
|
{"invalid_stellar_corner_samples", 0.0}, {"kinetic_energy", 1.0e-14},
|
|
{"negative_density_maximum_scaled", 1.0e-8}, {"negative_enthalpy_maximum_scaled", 1.0e-8}
|
|
};
|
|
if (options.profiles) budgets["profile_missing_points"] = 0.0;
|
|
if (options.profiles && options.mode == "analytic-mesh") budgets["profile_maximum_scaled_error"] = 1.0e-8;
|
|
if (options.profiles && options.mode == "analytic-mesh") budgets["profile_maximum_angular_rms_scaled"] = 1.0e-10;
|
|
if (metrics.contains("normalized_unbordered_residual")) budgets["normalized_unbordered_residual"] = 1.0e-8;
|
|
auto checks = File(options.output / "verification_checks.csv");
|
|
checks << "metric,observed,maximum_allowed,passed\n";
|
|
bool passed = true;
|
|
for (const auto &[name, budget] : budgets) {
|
|
const double value = metrics.at(name);
|
|
const bool okay = std::isfinite(value) && std::abs(value) <= budget;
|
|
checks << name << ',' << value << ',' << budget << ',' << okay << '\n';
|
|
passed = passed && okay;
|
|
if (!okay) std::cout << "Screen failed: " << name << '=' << value << " budget=" << budget << '\n';
|
|
}
|
|
for (const std::string name : {"mass", "binding_energy", "pressure_integral", "virial_ratio", "virial_error",
|
|
"force_virial_error", "density_relative_l2_error", "enthalpy_relative_l2_error",
|
|
"potential_relative_l2_error", "surface_radius_relative_rms_error"}) {
|
|
std::cout << name << '=' << metrics.at(name) << '\n';
|
|
}
|
|
std::cout << "Physical screening passed=" << passed << " (not a resolution-convergence claim)" << std::endl;
|
|
return passed;
|
|
}
|
|
|
|
double BlockNorm(const mfem::Vector &vector, const int offset, const int size) {
|
|
long double sum = 0.0L;
|
|
for (int index = offset; index < offset + size; ++index) sum += static_cast<long double>(vector(index)) * vector(index);
|
|
return std::sqrt(sum);
|
|
}
|
|
|
|
int Run(const Options &options) {
|
|
if (options.mode == "help") return 0;
|
|
if (!std::filesystem::create_directory(options.output)) {
|
|
throw std::invalid_argument("Output directory already exists; choose a new --output directory.");
|
|
}
|
|
auto metadata = File(options.output / "metadata.txt");
|
|
const N1Reference reference{utils::G, utils::MASS, utils::RADIUS};
|
|
reference.Validate();
|
|
metadata << "mode=" << options.mode << "\nmesh=" << std::filesystem::absolute(options.mesh).string()
|
|
<< "\ncompiled=" << __DATE__ << ' ' << __TIME__ << "\ncompiler=" << __VERSION__
|
|
<< "\nmpi_ranks=1\nmodel=nonrotating_n1_fixed_mass_fixed_central_density_zero_surface_pressure"
|
|
<< "\nG=" << reference.gravitationalConstant << "\nM=" << reference.mass << "\nR=" << reference.radius
|
|
<< "\nK=" << reference.PolytropicConstant() << "\nrho_c=" << reference.CentralDensity()
|
|
<< "\nabsolute_tolerance=" << options.absoluteTolerance << "\nrelative_tolerance=" << options.relativeTolerance
|
|
<< "\nlinear_tolerance=" << options.linearTolerance << "\nmax_newton=" << options.maximumNewtonIterations
|
|
<< "\nmax_linear_iterations=" << options.maximumLinearIterations
|
|
<< "\npolynomial_increment=" << MEAN_FIELD_UNIFORM_POLYNOMIAL_ORDER_INCREMENT
|
|
<< "\nnormalization=production_frozen_physical_Riesz_diagonal\nprofiles=" << options.profiles << '\n';
|
|
metadata << "mu_points=" << options.radial.muPointCount << "\nphi_points=" << options.radial.azimuthPointCount
|
|
<< "\nexterior_shells=" << options.radial.exteriorShellCount << '\n';
|
|
metadata.flush();
|
|
if (!SelfChecks(options)) return 3;
|
|
if (options.mode == "self-check") return 0;
|
|
|
|
const auto meshSnapshot = options.output / "input.smesh";
|
|
std::filesystem::copy_file(options.mesh, meshSnapshot);
|
|
metadata << "mesh_snapshot=" << std::filesystem::absolute(meshSnapshot).string() << '\n';
|
|
utils::Args arguments;
|
|
arguments.mesh_file = meshSnapshot.string();
|
|
arguments.p.rtol = arguments.p.atol = 1.0e-12;
|
|
auto finiteElements = fem::setup_fem(arguments.mesh_file, arguments, 0);
|
|
if (!finiteElements.okay()) throw std::runtime_error("Could not construct finite elements.");
|
|
metadata << "mesh_bytes=" << std::filesystem::file_size(meshSnapshot)
|
|
<< "\nelements=" << finiteElements.mesh->GetNE()
|
|
<< "\ndensity_order=" << finiteElements.densityFes->GetMaxElementOrder()
|
|
<< "\nenthalpy_order=" << finiteElements.enthalpyFes->GetMaxElementOrder()
|
|
<< "\npotential_order=" << finiteElements.gravityPotentialFes->GetMaxElementOrder()
|
|
<< "\ngravity_flux_order=" << finiteElements.gravityFluxFes->GetMaxElementOrder()
|
|
<< "\ndisplacement_order=" << finiteElements.displacementFes->GetMaxElementOrder() << '\n';
|
|
metadata.flush();
|
|
if (options.mode == "analytic-mesh") {
|
|
AnalyticMeshState analytic(finiteElements, reference);
|
|
const bool passed = Measure(analytic, reference, options);
|
|
metadata << "physical_screen_passed=" << passed << '\n';
|
|
return passed ? 0 : 3;
|
|
}
|
|
if (options.mode == "replay") {
|
|
const auto source = ReadMetadata(options.replay);
|
|
if (source.at("mode") != "solve" || source.at("model") != "nonrotating_n1_fixed_mass_fixed_central_density_zero_surface_pressure" ||
|
|
std::stod(source.at("G")) != reference.gravitationalConstant || std::stod(source.at("M")) != reference.mass ||
|
|
std::stod(source.at("R")) != reference.radius) {
|
|
throw std::runtime_error("Replay source does not match this nonrotating n=1 experiment.");
|
|
}
|
|
auto savedMetrics = ReadSavedSolverMetrics(options.replay);
|
|
PhysicalState physical(finiteElements, options.replay);
|
|
const bool converged = source.at("solver_converged") == "1";
|
|
metadata << "replay_source=" << std::filesystem::absolute(options.replay).string()
|
|
<< "\nsolver_converged=" << converged
|
|
<< "\nsolver_diagnostics=copied_from_source_not_recomputed\n";
|
|
if (!converged) metadata << "solver_failure=" << source.at("solver_failure") << '\n';
|
|
metadata.flush();
|
|
const bool passed = Measure(physical, reference, options, std::move(savedMetrics));
|
|
metadata << "physical_screen_passed=" << passed << '\n';
|
|
return !converged ? 1 : (passed ? 0 : 3);
|
|
}
|
|
|
|
auto stellarModel = model::StellarModel(
|
|
eos::Polytrope({.n = 1.0, .K = reference.PolytropicConstant()}),
|
|
surface::Isobaric({.Psurf = dimensions::PressureValue{0.0}}),
|
|
integral::FixedTotalMass({.Mtotal = dimensions::MassValue{reference.mass}}),
|
|
integral::FixedAngularMomentum({.Jtotal = dimensions::AngularMomentumValue{0.0},
|
|
.axis = {0.0, 0.0, 1.0}, .center = {0.0, 0.0, 0.0}}),
|
|
constraint::FixedCentralDensity({.RhoC = dimensions::DensityValue{reference.CentralDensity()}})
|
|
);
|
|
auto discretization = equilibrium::makeStellarDiscretization(std::move(finiteElements),
|
|
normalization::PhysicalRieszDiagonal{dimensions::LengthValue{reference.radius}, reference.gravitationalConstant});
|
|
std::cout << "Constructing production context; nonrotating n=1, nonlinear atol=" << options.absoluteTolerance << std::endl;
|
|
const auto start = Clock::now();
|
|
auto context = solver::makeContext(std::move(stellarModel), std::move(discretization),
|
|
preconditioning::makePreconditioner(), solver::linear::FGMRES({.restartLength = 40, .printLevel = -1}));
|
|
metadata << "context_seconds=" << std::chrono::duration<double>(Clock::now()-start).count() << '\n';
|
|
metadata.flush();
|
|
// Exercise accepted-field reconstruction before the expensive solve,
|
|
// and retain a seed baseline to detect physical degradation by Newton.
|
|
std::cout << "Measuring production seed before Newton" << std::endl;
|
|
solver::detail::StellarEquilibriumContextDiagnostics::WithState(context, [&](auto &state, const fem::FEM &fem) {
|
|
PhysicalState physical(state, fem);
|
|
auto seedMetrics = MeasureVolumes(physical, reference, options.quadratureOrder);
|
|
seedMetrics["normalized_bordered_residual"] = physical.normalizedBorderedResidualNorm;
|
|
seedMetrics["normalized_unbordered_residual"] = physical.normalizedUnborderedResidualNorm;
|
|
seedMetrics["central_border"] = physical.centralBorder;
|
|
WriteMetrics(options.output / "seed_physical_metrics.csv", seedMetrics);
|
|
std::cout << "Seed |F|=" << physical.normalizedBorderedResidualNorm
|
|
<< ", density relative L2=" << seedMetrics.at("density_relative_l2_error")
|
|
<< ", virial error=" << seedMetrics.at("virial_error") << std::endl;
|
|
});
|
|
auto trajectory = File(options.output / "newton_history.csv");
|
|
trajectory << "iteration,residual,step,trials,iteration_seconds,linear_iterations,linear_relative_residual,linear_seconds\n";
|
|
auto observer = solver::nonlinear::makeObserver(
|
|
[](const solver::nonlinear::BeforeIteration &event) {
|
|
std::cout << "Newton " << event.iteration << ": |F|=" << event.residualNorm << std::endl;
|
|
},
|
|
[&](const solver::nonlinear::AfterIteration &event) {
|
|
trajectory << event.iteration << ',' << event.residualNorm << ',' << event.acceptedStepLength << ','
|
|
<< event.lineSearchTrials << ',' << event.iterationSeconds << ','
|
|
<< (event.linearSolve ? event.linearSolve->iterations : 0) << ','
|
|
<< (event.linearSolve ? event.linearSolve->relativeTrueResidualNorm : 0.0) << ','
|
|
<< (event.linearSolve ? event.linearSolve->solveSeconds : 0.0) << '\n';
|
|
trajectory.flush();
|
|
std::cout << " step=" << event.acceptedStepLength << ", |F|=" << event.residualNorm
|
|
<< ", seconds=" << event.iterationSeconds;
|
|
if (event.linearSolve) std::cout << ", linear_iterations=" << event.linearSolve->iterations
|
|
<< ", true_relative=" << event.linearSolve->relativeTrueResidualNorm;
|
|
std::cout << std::endl;
|
|
}
|
|
);
|
|
auto nonlinear = solver::nonlinear::Newton(solver::nonlinear::NewtonOptions{
|
|
.relativeTolerance = options.relativeTolerance, .absoluteTolerance = options.absoluteTolerance,
|
|
.maximumIterations = options.maximumNewtonIterations,
|
|
.linearSolve = {.relativeTolerance = options.linearTolerance, .absoluteTolerance = 0.0,
|
|
.maximumIterations = options.maximumLinearIterations}, .backtracking = {}
|
|
});
|
|
const auto report = [&]() {
|
|
auto equilibriumSolver = solver::make(context, nonlinear, observer);
|
|
return equilibriumSolver.evaluate();
|
|
}(); // Release the solver before the guarded diagnostic callback.
|
|
metadata << "solver_converged=" << report.converged()
|
|
<< "\naccepted_steps=" << report.completedNonlinearIterations() << '\n';
|
|
if (!report.converged()) metadata << "solver_failure=" << report.failure().message << '\n';
|
|
metadata.flush();
|
|
std::cout << "Solver converged=" << report.converged() << "; measuring last accepted state" << std::endl;
|
|
bool passed = false;
|
|
solver::detail::StellarEquilibriumContextDiagnostics::WithState(context, [&](auto &state, const fem::FEM &fem) {
|
|
// Preserve expensive accepted coefficients before postprocessing.
|
|
// This is experiment data, not a versioned production checkpoint.
|
|
auto accepted = File(options.output / "accepted_state.txt");
|
|
accepted << state.AcceptedPhysicalState().Size() << '\n';
|
|
state.AcceptedPhysicalState().Print(accepted, 1);
|
|
accepted.close();
|
|
auto layout = File(options.output / "state_layout.csv");
|
|
layout << "block,offset,size\n";
|
|
for (const auto &block : state.Problem().GetManifest().valueBlocks()) {
|
|
layout << block.stableId << ',' << block.offset << ',' << block.size << '\n';
|
|
}
|
|
PhysicalState physical(state, fem);
|
|
auto saveField = [&](const char *name, const mfem::ParGridFunction &field) {
|
|
auto stream = File(options.output / (std::string(name) + ".gf"));
|
|
field.Save(stream);
|
|
};
|
|
saveField("density", physical.density);
|
|
saveField("enthalpy", physical.enthalpy);
|
|
saveField("potential", physical.potential);
|
|
saveField("gravity_gradient_reference", physical.gravityGradientReference);
|
|
saveField("displacement", physical.displacement);
|
|
auto residuals = File(options.output / "residual_blocks.csv");
|
|
residuals << "block,size,physical_bordered_l2,physical_unbordered_l2,normalized_bordered_l2,normalized_unbordered_l2\n";
|
|
for (const auto &block : state.Problem().GetManifest().residualBlocks()) {
|
|
residuals << block.stableId << ',' << block.size << ','
|
|
<< BlockNorm(physical.physicalResidualBordered, block.offset, block.size) << ','
|
|
<< BlockNorm(physical.physicalResidualUnbordered, block.offset, block.size) << ','
|
|
<< BlockNorm(physical.normalizedBorderedResidual, block.offset, block.size) << ','
|
|
<< BlockNorm(physical.normalizedUnborderedResidual, block.offset, block.size) << '\n';
|
|
}
|
|
auto reconstruction = File(options.output / "field_reconstruction.csv");
|
|
reconstruction << "field,reduced_size,full_true_size,maximum_round_trip_error\n";
|
|
for (const auto &field : physical.reconstructionReports) {
|
|
reconstruction << field.field << ',' << field.reducedSize << ',' << field.fullTrueSize << ',' << field.maximumRoundTripError << '\n';
|
|
}
|
|
metadata << "state_size=" << state.AcceptedPhysicalState().Size()
|
|
<< "\naccepted_normalized_residual=" << physical.normalizedBorderedResidualNorm << '\n';
|
|
if (physical.centralDensityReport) {
|
|
metadata << "central_constraint_density_inferred_from_h=" << physical.centralDensityReport->achievedDensity
|
|
<< "\ncentral_constraint_h=" << physical.centralDensityReport->achievedEnthalpy << '\n';
|
|
}
|
|
metadata.flush();
|
|
passed = Measure(physical, reference, options, {
|
|
{"normalized_bordered_residual", physical.normalizedBorderedResidualNorm},
|
|
{"normalized_unbordered_residual", physical.normalizedUnborderedResidualNorm},
|
|
{"normalized_central_border_action", physical.normalizedCentralBorderActionNorm},
|
|
{"central_border", physical.centralBorder}, {"bernoulli_constant", physical.bernoulliConstant},
|
|
{"angular_velocity_norm", physical.rotation.angular_velocity().Norml2()}
|
|
});
|
|
});
|
|
metadata << "physical_screen_passed=" << passed << "\ntotal_seconds=" << std::chrono::duration<double>(Clock::now()-start).count() << '\n';
|
|
return !report.converged() ? 1 : (passed ? 0 : 3);
|
|
}
|
|
} // namespace
|
|
|
|
int main(int argc, char **argv) {
|
|
mfem::Mpi::Init(argc, argv);
|
|
int result = 0;
|
|
try {
|
|
int ranks = 0;
|
|
MPI_Comm_size(MPI_COMM_WORLD, &ranks);
|
|
if (ranks != 1) throw std::invalid_argument("Run this verification on exactly one MPI rank.");
|
|
mfem::Device device("cpu");
|
|
std::cout << std::setprecision(12);
|
|
result = Run(Parse(argc, argv));
|
|
} catch (const std::exception &error) {
|
|
std::cerr << "polytrope verification failure: " << error.what() << std::endl;
|
|
result = 2;
|
|
}
|
|
mfem::Mpi::Finalize();
|
|
return result;
|
|
}
|