diff --git a/CMakeLists.txt b/CMakeLists.txt index c57ee40..9bc0c81 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -410,6 +410,14 @@ add_executable(density_order_contract_tests EXCLUDE_FROM_ALL tests/operators/prepared_barotropic_closure.cpp) target_link_libraries(density_order_contract_tests PRIVATE mean_field test_mod Catch2::Catch2 Boost::boost) +# Continuation changes owned runtime state; exercise its transactional boundary +# without rebuilding unrelated geometry and equation-of-state test suites. +add_executable(continuation_contract_tests EXCLUDE_FROM_ALL + tests/test_main.cpp + tests/solver/stellar_numerical_strategies.cpp + tests/solver/continuation_snapshot.cpp) +target_link_libraries(continuation_contract_tests PRIVATE mean_field test_mod Catch2::Catch2 Boost::boost) + # A deliberately opt-in API workbench. It is compiled explicitly during API # verification but temporary user experiments do not break the default build. add_executable(sandbox EXCLUDE_FROM_ALL sandbox.cpp) diff --git a/experiments/continuation_snapshot.hpp b/experiments/continuation_snapshot.hpp new file mode 100644 index 0000000..7136878 --- /dev/null +++ b/experiments/continuation_snapshot.hpp @@ -0,0 +1,259 @@ +#pragma once + +// Include after import mean_field. These helpers transfer an initial guess; +// they do not certify that it solves the destination model. The caller binds +// model, mesh, field orders and surface-basis provenance before installation. +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace experiment::continuation { +using SnapshotMetadata = std::map; + +struct SnapshotBlock final { + std::string name; + int offset; + int size; +}; + +namespace detail { +inline constexpr std::string_view centralBorder = "fixed_central_density.border"; + +inline bool Finite(const mfem::Vector &vector) { + for (int i = 0; i < vector.Size(); ++i) + if (!std::isfinite(vector(i))) return false; + return true; +} + +inline bool Identical(const mfem::Vector &left, const mfem::Vector &right) { + return left.Size() == right.Size() && + (left.Size() == 0 || std::memcmp(left.GetData(), right.GetData(), + sizeof(mfem::real_t) * static_cast(left.Size())) == 0); +} + +inline void RequireSerial(MPI_Comm communicator) { + int ranks = 0; + if (communicator == MPI_COMM_NULL || MPI_Comm_size(communicator, &ranks) != MPI_SUCCESS || ranks != 1) + throw std::invalid_argument("Continuation snapshots currently require one MPI rank."); +} + +inline int NonnegativeInteger(std::string_view text) { + int value = -1; + const auto result = std::from_chars(text.data(), text.data() + text.size(), value); + if (text.empty() || result.ec != std::errc{} || result.ptr != text.data() + text.size() || value < 0) + throw std::invalid_argument("A snapshot layout offset or size is not a nonnegative integer."); + return value; +} + +inline std::vector Fields(std::string line) { + if (!line.empty() && line.back() == '\r') line.pop_back(); + if (line.find('"') != std::string::npos) + throw std::invalid_argument("Snapshot layout requires the producer's unquoted CSV format."); + std::vector result; + std::size_t first = 0; + while (true) { + const auto comma = line.find(',', first); + result.push_back(line.substr(first, comma == std::string::npos ? comma : comma - first)); + if (comma == std::string::npos) return result; + first = comma + 1; + } +} + +// Contiguity establishes an exact coordinate partition independently of file +// order. Equal total sizes alone cannot detect permuted or missing physics. +inline int ValidatePartition(std::span blocks) { + if (blocks.empty()) throw std::invalid_argument("Snapshot layout has no physical value blocks."); + std::vector ordered(blocks.begin(), blocks.end()); + std::set names; + for (const auto &block : ordered) { + if (block.name.empty() || !names.insert(block.name).second || block.offset < 0 || block.size < 0 || + block.offset > std::numeric_limits::max() - block.size) + throw std::invalid_argument("Snapshot layout has invalid or repeated block coordinates."); + if (block.name == centralBorder && block.size != 1) + throw std::invalid_argument("A central-density border must contain exactly one scalar."); + } + std::sort(ordered.begin(), ordered.end(), [](const auto &a, const auto &b) { + return a.offset < b.offset || (a.offset == b.offset && a.size < b.size); + }); + int size = 0; + for (const auto &block : ordered) { + if (block.offset != size) + throw std::invalid_argument("Snapshot layout contains a gap or overlapping coordinates."); + size += block.size; + } + if (size == 0) throw std::invalid_argument("Snapshot state is empty."); + return size; +} + +inline std::vector ReadValueBlocks(const std::filesystem::path &path) { + std::ifstream input(path); + std::string line; + if (!input || !std::getline(input, line)) + throw std::invalid_argument("Cannot read snapshot layout: " + path.string()); + const auto header = Fields(line); + std::map columns; + for (std::size_t i = 0; i < header.size(); ++i) + if (header[i].empty() || !columns.emplace(header[i], i).second) + throw std::invalid_argument("Snapshot layout has an invalid CSV header."); + for (const auto *name : {"kind", "block", "offset", "size"}) + if (!columns.contains(name)) throw std::invalid_argument("Snapshot layout lacks a required column."); + std::vector values, residuals; + while (std::getline(input, line)) { + const auto fields = Fields(line); + if (fields.size() != header.size()) + throw std::invalid_argument("Snapshot layout has an incomplete CSV row."); + const auto &kind = fields[columns.at("kind")]; + SnapshotBlock block{fields[columns.at("block")], NonnegativeInteger(fields[columns.at("offset")]), + NonnegativeInteger(fields[columns.at("size")])}; + if (kind == "value") values.push_back(std::move(block)); + else if (kind == "residual") residuals.push_back(std::move(block)); + else throw std::invalid_argument("Snapshot layout has an unknown coordinate kind."); + } + if (input.bad()) throw std::runtime_error("Snapshot layout read failed."); + ValidatePartition(values); + if (!residuals.empty()) ValidatePartition(residuals); + return values; +} + +inline mfem::Vector ReadVector(const std::filesystem::path &path, int expectedSize) { + std::ifstream input(path); + long long count = -1; + if (!(input >> count) || count != expectedSize) + throw std::invalid_argument("Snapshot vector length disagrees with its named layout: " + path.string()); + mfem::Vector values(expectedSize); + for (int i = 0; i < expectedSize; ++i) + if (!(input >> values(i)) || !std::isfinite(values(i))) + throw std::invalid_argument("Snapshot vector is incomplete or nonfinite: " + path.string()); + std::string trailing; + if (input >> trailing || input.bad()) + throw std::invalid_argument("Snapshot vector has trailing data or a read error: " + path.string()); + return values; +} +} // namespace detail + +inline SnapshotMetadata ParseSnapshotMetadata(const std::filesystem::path &directory) { + const auto path = directory / "metadata.txt"; + std::ifstream input(path); + if (!input) throw std::invalid_argument("Continuation requires snapshot metadata: " + path.string()); + SnapshotMetadata metadata; + std::string line; + while (std::getline(input, line)) { + if (!line.empty() && line.back() == '\r') line.pop_back(); + if (line.empty()) continue; + const auto separator = line.find('='); + if (separator == std::string::npos || separator == 0 || + !metadata.emplace(line.substr(0, separator), line.substr(separator + 1)).second) + throw std::invalid_argument("Snapshot metadata has a malformed or repeated key."); + } + if (input.bad() || metadata.empty()) throw std::invalid_argument("Snapshot metadata is empty or unreadable."); + return metadata; +} + +inline mfem::Vector MapNamedPhysicalBlocks( + const mfem::Vector &source, std::span sourceBlocks, + std::span targetBlocks, bool allowCentralBorderChange = false +) { + if (source.Size() != detail::ValidatePartition(sourceBlocks) || !detail::Finite(source)) + throw std::invalid_argument("Snapshot source does not match its finite physical coordinate layout."); + mfem::Vector target(detail::ValidatePartition(targetBlocks)); + target = 0.0; + std::map sourceByName, targetByName; + for (const auto &block : sourceBlocks) sourceByName.emplace(block.name, &block); + for (const auto &block : targetBlocks) targetByName.emplace(block.name, &block); + for (const auto &block : sourceBlocks) + if (!targetByName.contains(block.name) && !(allowCentralBorderChange && block.name == detail::centralBorder)) + throw std::invalid_argument("Snapshot contains an unknown destination block: " + block.name); + for (const auto &block : targetBlocks) { + const auto found = sourceByName.find(block.name); + if (found == sourceByName.end()) { + if (allowCentralBorderChange && block.name == detail::centralBorder) continue; + throw std::invalid_argument("Snapshot is missing a required destination block: " + block.name); + } + const auto &from = *found->second; + if (from.size != block.size) + throw std::invalid_argument("Snapshot block has a different coordinate count: " + block.name); + for (int i = 0; i < block.size; ++i) target(block.offset + i) = source(from.offset + i); + } + return target; +} + +template +mfem::Vector LoadMappedSnapshot(const Problem &problem, const std::filesystem::path &path, + bool allowCentralBorderChange = false, bool allowUnconvergedRestart = false) { + detail::RequireSerial(problem.GetCommunicator()); + const auto statePath = std::filesystem::is_directory(path) ? path / "accepted_state.txt" : path; + const auto directory = std::filesystem::exists(statePath.parent_path() / "metadata.txt") ? + statePath.parent_path() : statePath.parent_path().parent_path(); + const bool isAcceptedEndpoint = statePath.filename() == "accepted_state.txt" && statePath.parent_path() == directory; + if (!isAcceptedEndpoint && !allowUnconvergedRestart) + throw std::invalid_argument("Captured iterates require the explicit unconverged-restart override."); + const auto metadata = ParseSnapshotMetadata(directory); + const auto converged = metadata.find("converged"); + if (converged != metadata.end() && converged->second != "0" && converged->second != "1") + throw std::invalid_argument("Snapshot metadata has an invalid convergence flag."); + if (!allowUnconvergedRestart && (converged == metadata.end() || converged->second != "1")) + throw std::invalid_argument("Continuation refuses an unverified or unconverged predecessor; an explicit restart override is required."); + const auto sourceBlocks = detail::ReadValueBlocks(directory / "layout.csv"); + const auto source = detail::ReadVector(statePath, detail::ValidatePartition(sourceBlocks)); + std::vector targetBlocks; + for (const auto &block : problem.GetManifest().valueBlocks()) + targetBlocks.push_back({std::string(block.stableId), block.offset, block.size}); + auto mapped = MapNamedPhysicalBlocks(source, sourceBlocks, targetBlocks, allowCentralBorderChange); + if (mapped.Size() != problem.StateSize()) + throw std::logic_error("Destination manifest does not cover its physical state."); + return mapped; +} + +// No mesh or FEM preparation is needed: these cases catch the silent block +// reorder/border mistakes that a length-only snapshot reader would permit. +inline int RunSnapshotCoordinateSelfChecks() { + int checks = 0; + const auto require = [&](bool condition) { + if (!condition) throw std::logic_error("Continuation coordinate self-check failed."); + ++checks; + }; + const auto rejects = [&](auto action) { + bool rejected = false; + try { action(); } catch (const std::invalid_argument &) { rejected = true; } + require(rejected); + }; + const std::vector source{{"density", 0, 2}, {"enthalpy", 2, 1}}; + const std::vector reordered{{"enthalpy", 0, 1}, {"density", 1, 2}}; + const std::vector bordered{{"enthalpy", 0, 1}, {std::string(detail::centralBorder), 1, 1}, {"density", 2, 2}}; + mfem::Vector values(3); + values(0) = 2.; values(1) = 3.; values(2) = 7.; + const auto mapped = MapNamedPhysicalBlocks(values, source, reordered); + require(mapped(0) == 7. && mapped(1) == 2. && mapped(2) == 3.); + const auto added = MapNamedPhysicalBlocks(values, source, bordered, true); + require(added.Size() == 4 && added(0) == 7. && added(1) == 0. && added(2) == 2. && added(3) == 3.); + require(detail::Identical(MapNamedPhysicalBlocks(added, bordered, source, true), values)); + rejects([&] { MapNamedPhysicalBlocks(values, source, bordered); }); + rejects([&] { MapNamedPhysicalBlocks(added, bordered, source); }); + auto bad = source; bad[1].name = "unknown"; + rejects([&] { MapNamedPhysicalBlocks(values, bad, source); }); + bad = source; bad[1].offset = 1; + rejects([&] { MapNamedPhysicalBlocks(values, bad, source); }); + bad = source; bad[1].offset = 3; + rejects([&] { MapNamedPhysicalBlocks(values, bad, source); }); + bad = source; bad[1].name = "density"; + rejects([&] { MapNamedPhysicalBlocks(values, bad, source); }); + bad = source; bad[0].size = 1; bad[1].offset = 1; bad[1].size = 2; + rejects([&] { MapNamedPhysicalBlocks(values, source, bad); }); + values(0) = std::numeric_limits::quiet_NaN(); + rejects([&] { MapNamedPhysicalBlocks(values, source, reordered); }); + return checks; +} +} // namespace experiment::continuation diff --git a/experiments/harmonic_linear_trials.hpp b/experiments/harmonic_linear_trials.hpp new file mode 100644 index 0000000..062869f --- /dev/null +++ b/experiments/harmonic_linear_trials.hpp @@ -0,0 +1,78 @@ +#pragma once + +// Include after the benchmark and physical-snapshot helpers. Each observation +// starts from the same prepared physical state; nothing is accepted or clipped. +namespace experiment::plateau_diagnostics { +template +void RecordHarmonicLinearTrials(State &state, const mean_field::fem::FEM &fem, + const std::filesystem::path &output, double index) { + namespace d = preconditioner_benchmark_detail; + const mfem::Vector base(state.trialNormalizedState), residual(state.trialNormalizedResidual), + savedCorrection(state.normalizedCorrection), accepted(state.AcceptedPhysicalState()); + mfem::Vector physicalBase; + state.normalizedOperator->DenormalizeState(base, physicalBase); + const auto restore = [&] { + state.trialNormalizedState = base; + if (!state.PrepareTrial()) throw std::runtime_error("Linear trial base cannot be restored"); + }; + auto trials = d::Open(output / "linear_direction_trials.csv"); + trials << "role,alpha,prepared,reason,stage,base_norm,trial_norm,relative_norm,linear_defect_relative," + "direction_norm,action_norm,merit_slope,remainder_norm,remainder_over_alpha," + "predicted_norm,signed_remainder_projection,merit_identity_error,decrease_screen\n"; + auto blocks = d::Open(output / "linear_direction_trial_blocks.csv"); + blocks << "role,alpha,block,base_norm,action_norm,trial_norm,remainder_norm\n"; + auto bounds = d::Open(output / "linear_direction_enthalpy_bounds.csv"); + bounds << "role,metric,value\n"; + for (const std::string role : {"baseline", "candidate"}) { + restore(); + const auto correction = coupled_physical_detail::Load(output / (role + "_normalized_correction.txt"), base.Size()); + mfem::Vector action(base.Size()), defect(base.Size()); + state.normalizedOperator->Mult(correction, action); + defect = action; defect += residual; + const double relative = defect.Norml2() / residual.Norml2(); + const double slope = residual * action; + state.normalizedCorrection = correction; + for (const auto &[key, value] : MeasureEnthalpyStep(state, fem, index, &physicalBase)) + bounds << role << ',' << key << ',' << value << '\n'; + d::Save(output / (role + "_trial_jacobian_action.txt"), action); + for (const double alpha : {1., .25, .0625, .015625, .00390625, .0009765625, .000244140625}) { + restore(); + state.trialNormalizedState = base; + state.trialNormalizedState.Add(alpha, correction); + const auto prepared = state.PrepareTrial(); + trials << role << ',' << alpha << ',' << prepared.has_value() << ',' + << (prepared ? -1 : int(prepared.error().reason)) << ',' + << (prepared ? -1 : int(prepared.error().stage)) << ',' << residual.Norml2() << ','; + if (!prepared) { + trials << "nan,nan," << relative << ',' << correction.Norml2() << ',' << action.Norml2() + << ',' << slope << ",nan,nan,nan,nan,nan,0\n"; + continue; + } + const mfem::Vector trial(state.trialNormalizedResidual); + mfem::Vector prediction(residual), remainder(trial); + prediction.Add(alpha, action); + remainder -= prediction; + const double projection = prediction * remainder; + const double identity = trial * trial - prediction * prediction - 2. * projection - remainder * remainder; + trials << trial.Norml2() << ',' << trial.Norml2() / residual.Norml2() << ',' << relative << ',' + << correction.Norml2() << ',' << action.Norml2() << ',' << slope << ',' << remainder.Norml2() + << ',' << remainder.Norml2() / alpha << ',' << prediction.Norml2() << ',' << projection << ',' + << identity << ',' << (trial * trial <= (1. - 1e-4 * alpha) * (residual * residual)) << '\n'; + for (const auto &block : state.Problem().GetManifest().residualBlocks()) + blocks << role << ',' << alpha << ',' << block.stableId << ',' + << d::BlockNorm(residual, block.offset, block.size) << ',' + << d::BlockNorm(action, block.offset, block.size) << ',' + << d::BlockNorm(trial, block.offset, block.size) << ',' + << d::BlockNorm(remainder, block.offset, block.size) << '\n'; + } + trials.flush(); blocks.flush(); bounds.flush(); + } + restore(); + state.normalizedCorrection = savedCorrection; + if (!d::Identical(accepted, state.AcceptedPhysicalState()) || + !d::Identical(residual, state.trialNormalizedResidual)) + throw std::runtime_error("Linear trial observations altered their base or accepted state"); + auto invariants = d::Open(output / "linear_direction_trial_invariants.csv"); + invariants << "metric,value\naccepted_unchanged,1\nbase_residual_restored,1\n"; +} +} // namespace experiment::plateau_diagnostics diff --git a/experiments/plateau_diagnostic.cpp b/experiments/plateau_diagnostic.cpp index 8d86bcd..6cf9b07 100644 --- a/experiments/plateau_diagnostic.cpp +++ b/experiments/plateau_diagnostic.cpp @@ -36,6 +36,8 @@ import mean_field; #include "negative_density_diagnostics.hpp" #include "negative_density_regions.hpp" #include "density_order_projection.hpp" +#include "continuation_snapshot.hpp" +#include "harmonic_linear_trials.hpp" using namespace mean_field; namespace { @@ -70,7 +72,8 @@ struct Options { std::vector captureIterations; std::string mesh = "sandbox.smesh"; std::filesystem::path output, probeState, probeDirection, replayCandidate; - std::filesystem::path densitySnapshot; + std::filesystem::path densitySnapshot, restartSnapshot; + bool continuation = false, allowUnconvergedRestart = false, benchmarkTrials = false; }; Options Parse(int argc, char **argv) { Options o; @@ -143,6 +146,14 @@ Options Parse(int argc, char **argv) { } else if (a == "--density-order-final") o.densityOrderFinal = true; + else if (a == "--restart-snapshot") + o.restartSnapshot = v(); + else if (a == "--continuation") + o.continuation = true; + else if (a == "--allow-unconverged-restart") + o.allowUnconvergedRestart = true; + else if (a == "--benchmark-trials") + o.benchmarkTrials = true; else if (a == "--density-snapshot") o.densitySnapshot = v(); else if (a == "--surface-prescription") @@ -213,6 +224,13 @@ Options Parse(int argc, char **argv) { } if (o.output.empty() || o.n < 1 || o.n >= 5 || o.maxNewton < 0 || o.maxLinear < 1 || o.restart < 1) throw std::invalid_argument("Invalid options"); + if ((o.continuation || o.allowUnconvergedRestart) && o.restartSnapshot.empty()) + throw std::invalid_argument("Continuation/restart permission requires --restart-snapshot"); + if (!o.restartSnapshot.empty() && (o.densityDiagnostic || !o.probeState.empty() || + !o.probeDirection.empty() || !o.replayCandidate.empty() || o.surfacePrescription != "harmonics")) + 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"); if (o.surfacePrescription != "nodal" && o.surfacePrescription != "harmonics") throw std::invalid_argument("Surface prescriptions are nodal or harmonics"); if (!o.densitySnapshot.empty() && !o.densityDiagnostic) @@ -234,10 +252,10 @@ Options Parse(int argc, char **argv) { // size match is insufficient to establish modal snapshot provenance. if (o.reducedSurfaceProbe || o.reducedSurfaceAlgebraCheck || o.weakShapeProbe || o.weakShapeFinal || o.probes || o.pressureProbe || o.shapeCurvature || o.surfaceTrace || o.strategyProbe || - o.positivityProbe || o.positivityCertificate || o.preconditionerProbe || o.nonlinearStepProbe || + o.positivityProbe || o.positivityCertificate || o.coupledTrialProbe || o.dampedCoupledProbe || !o.probeState.empty() || !o.probeDirection.empty() || - !o.replayCandidate.empty() || o.strategy != "baseline") - throw std::invalid_argument("Native harmonic diagnostics currently support live baseline Newton trajectories only"); + !o.replayCandidate.empty() || (o.strategy != "baseline" && o.strategy != "curvature")) + throw std::invalid_argument("Native harmonic diagnostics exclude legacy nodal probes and unsupported strategies"); } else if (o.surfaceMaximumOrder != -1) { throw std::invalid_argument("--surface-maximum-order requires --surface-prescription harmonics"); } @@ -656,6 +674,36 @@ int RunContext(Context &context, const Options &o, const Scales &scales) { mfem::Vector residualFactors, stateFactors; std::vector centralRows; int centralOffset = -1; + if (!o.restartSnapshot.empty()) { + mfem::Vector restartPhysical; + solver::detail::StellarEquilibriumContextDiagnostics::WithState(context, [&](auto &s, const fem::FEM &) { + const auto directory = std::filesystem::is_directory(o.restartSnapshot) ? o.restartSnapshot : + (std::filesystem::exists(o.restartSnapshot.parent_path() / "metadata.txt") ? + o.restartSnapshot.parent_path() : o.restartSnapshot.parent_path().parent_path()); + const auto saved = experiment::continuation::ParseSnapshotMetadata(directory); + const auto number = [&](const char *key, double expected) { + const double actual = std::stod(saved.at(key)); + if (!std::isfinite(actual) || std::abs(actual - expected) > + 1e-13 * std::max(1., std::abs(expected))) + throw std::invalid_argument(std::string("Restart metadata mismatch: ") + key); + }; + if (saved.at("surface_prescription") != o.surfacePrescription) + throw std::invalid_argument("Restart surface coordinate prescription differs"); + number("surface_degree", o.surfaceDegree); number("surface_maximum_order", o.surfaceMaximumOrder); + number("density_order", field::Density::scalarOrder); number("enthalpy_order", field::Enthalpy::scalarOrder); + number("gravity_potential_order", field::Gravity::potentialOrder); + number("gravity_flux_order", field::Gravity::fluxOrder); + number("displacement_order", field::Displacement::Vector::familyOrder); + number("M", 1.); number("R_reference", 1.); + if (!o.continuation) { + number("n", o.n); number("K", scales.K); number("J", o.angularMomentum); + number("central_phase", !o.noCentral); + } + restartPhysical = experiment::continuation::LoadMappedSnapshot(s.Problem(), + o.restartSnapshot, o.continuation, o.allowUnconvergedRestart); + }); + solver::detail::StellarEquilibriumContextDiagnostics::InitializePhysicalState(context, restartPhysical); + } const bool contextUsesGravitySchur = !o.preconditionerProbe && (o.preconditionerVariant == "gravity-schur" || o.preconditionerVariant == "stellar-ldu-gravity-schur"); auto metadata = File(o.output / "metadata.txt"); @@ -697,6 +745,9 @@ 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 << "restart_snapshot=" << o.restartSnapshot.string() << "\ncontinuation=" << o.continuation + << "\nallow_unconverged_restart=" << o.allowUnconvergedRestart + << "\nbenchmark_trials=" << o.benchmarkTrials << '\n'; metadata << "density_order=" << field::Density::scalarOrder << "\nenthalpy_order=" << field::Enthalpy::scalarOrder << "\ngravity_potential_order=" << field::Gravity::potentialOrder @@ -880,9 +931,10 @@ int RunContext(Context &context, const Options &o, const Scales &scales) { } // Custom Newton recipes do not instantiate the nine-variant frozen suite. if constexpr (FullDefaultStack) if (o.preconditionerProbe) { - metadata << "benchmark_input=" << (o.probeState.empty() ? "lane_emden_seed" : "physical_snapshot") << '\n'; + metadata << "benchmark_input=" << (!o.restartSnapshot.empty() ? "checked_restart_snapshot" : + o.probeState.empty() ? "lane_emden_seed" : "physical_snapshot") << '\n'; if (!o.probeState.empty()) metadata << "probe_state=" << o.probeState << '\n'; - solver::detail::StellarEquilibriumContextDiagnostics::WithState(context, [&](auto &s, const fem::FEM &) { + solver::detail::StellarEquilibriumContextDiagnostics::WithState(context, [&](auto &s, const fem::FEM &f) { if (o.probeState.empty()) { s.trialNormalizedState = s.acceptedNormalizedState; } else { @@ -899,6 +951,8 @@ int RunContext(Context &context, const Options &o, const Scales &scales) { {.variant = o.preconditionerVariant, .restart = o.restart, .maximumIterations = o.maxLinear, .relativeTolerance = o.linearTolerance, .schurCorrectionSweeps = o.schurCorrectionSweeps, .schurRelaxation = o.schurRelaxation}); + if (o.benchmarkTrials) + experiment::plateau_diagnostics::RecordHarmonicLinearTrials(s, f, o.output, o.n); }); // A missed linear gate is a measured outcome, not failure to execute // the experiment. The CSV retains both recurrence and true gate results. diff --git a/experiments/run_plateau_campaign.py b/experiments/run_plateau_campaign.py index 6b580e2..a00cbdb 100644 --- a/experiments/run_plateau_campaign.py +++ b/experiments/run_plateau_campaign.py @@ -132,9 +132,29 @@ def main(): path = source / path path = path.resolve() state_files[str(path)] = digest(path) - elif word == '--density-snapshot': + elif word in ('--density-snapshot', '--restart-snapshot'): path = Path(case['arguments'][i + 1]) path = (path if path.is_absolute() else source / path).resolve() + if word == '--restart-snapshot': + snapshot_file = path / 'accepted_state.txt' if path.is_dir() else path + folder = snapshot_file.parent + if not (folder / 'metadata.txt').is_file(): + folder = folder.parent + prior_manifest = folder.parent / 'provenance.json' + prior_run = folder.parent / (folder.name + '_run.json') + prior = json.loads(prior_manifest.read_text()) + saved_run = json.loads(prior_run.read_text()) + if prior['mesh_sha256'] != provenance['mesh_sha256']: + raise RuntimeError('Restart mesh differs from the snapshot campaign.') + for saved_file in (snapshot_file, folder / 'metadata.txt', folder / 'layout.csv'): + value = digest(saved_file) + key = str(saved_file.relative_to(folder)) + if value != saved_run['output_hashes'][key]: + raise RuntimeError('Restart input no longer matches its original run: ' + key) + state_files[str(saved_file)] = value + state_files[str(prior_manifest)] = digest(prior_manifest) + state_files[str(prior_run)] = digest(prior_run) + continue # Modal coordinates also depend on the mesh and named block # layout. Bind replay to its original recorded campaign rather # than trusting a coincidentally matching vector length. diff --git a/libmeanfield/interface/solver/stellar_context.cppm b/libmeanfield/interface/solver/stellar_context.cppm index f2a629a..a49469d 100644 --- a/libmeanfield/interface/solver/stellar_context.cppm +++ b/libmeanfield/interface/solver/stellar_context.cppm @@ -840,6 +840,48 @@ export namespace mean_field::solver { export namespace mean_field::solver::detail { struct StellarEquilibriumContextDiagnostics final { + // Research continuation supplies physical coefficients in a newly + // compiled problem. Keep this explicit mutation separate from WithState: + // a prepared initial guess has no equilibrium certification, and a + // correction recycled from the preceding problem has no valid tangent. + template + static void InitializePhysicalState(Context &context, const mfem::Vector &physicalState) { + if (context.hasActiveSolver() || !context.isReady()) + throw std::logic_error("Initial-state installation requires a ready, inactive context."); + auto &state = *context.m_state; + bool valid = physicalState.Size() == state.Problem().StateSize(); + for (int i = 0; i < physicalState.Size(); ++i) + valid = valid && std::isfinite(physicalState(i)); + if (!AllStellarRanksSatisfy(valid, state.Problem().GetCommunicator())) + throw std::invalid_argument("Initial physical coefficients have invalid size or values."); + context.AcquireSolver(); + try { + state.BeginEvaluation(); + state.normalizedOperator->NormalizeState(physicalState, state.trialNormalizedState); + const auto preparation = state.PrepareTrial(); + if (!preparation) + throw std::invalid_argument("Initial physical coefficients fail operator preparation."); + bool finiteResidual = true; + for (int i = 0; i < state.trialNormalizedResidual.Size(); ++i) + finiteResidual = finiteResidual && std::isfinite(state.trialNormalizedResidual(i)); + if (!AllStellarRanksSatisfy(finiteResidual, state.Problem().GetCommunicator())) + throw std::invalid_argument("Initial physical coefficients produce a nonfinite residual."); + state.PrepareTrialCommit(); + state.CommitPreparedTrial(preparation->minimumJacobianDeterminant); + state.normalizedCorrection = 0.; + state.physicalCorrection = 0.; + state.volumeDisplacementDirection = 0.; + state.linearRightHandSide = 0.; + context.ReleaseSolver(); + } catch (...) { + const auto original = std::current_exception(); + try { state.RestoreAccepted(); } + catch (...) { context.ReleaseSolver(); throw; } + context.ReleaseSolver(); + std::rethrow_exception(original); + } + } + // The callback must not retain references to runtime storage. It may // prepare trial states, but accepted vectors must remain unchanged. // Restore the production preparation and correction on every exit. diff --git a/tests/solver/continuation_snapshot.cpp b/tests/solver/continuation_snapshot.cpp new file mode 100644 index 0000000..6e20373 --- /dev/null +++ b/tests/solver/continuation_snapshot.cpp @@ -0,0 +1,65 @@ +#include +#include +#include +#include +#include +#include +#include +#include +import mean_field; +#include "../../experiments/continuation_snapshot.hpp" + +namespace { +struct Manifest { + struct Block { std::string stableId; int offset, size; }; + std::vector blocks{{"specific_enthalpy", 0, 1}, {"density", 1, 2}}; + const auto &valueBlocks() const { return blocks; } +}; +struct Problem { + Manifest manifest; + MPI_Comm GetCommunicator() const { return MPI_COMM_SELF; } + const auto &GetManifest() const { return manifest; } + int StateSize() const { return 3; } +}; +struct SnapshotFiles { + const std::filesystem::path directory = std::filesystem::temp_directory_path() / + ("meanfield-continuation-" + std::to_string(std::chrono::steady_clock::now().time_since_epoch().count())); + SnapshotFiles() { + std::filesystem::create_directories(directory / "capture_0"); + write("metadata.txt", "converged=1\n"); + write("layout.csv", "kind,block,offset,size,scale\nvalue,density,0,2,1\nvalue,specific_enthalpy,2,1,1\n"); + write("accepted_state.txt", "3\n2\n3\n7\n"); + write("capture_0/physical_base.txt", "3\n2\n3\n7\n"); + } + ~SnapshotFiles() { std::error_code error; std::filesystem::remove_all(directory, error); } + void write(const char *name, const std::string &contents) { std::ofstream(directory / name) << contents; } +}; +} + +TEST_CASE("Continuation file loading checks provenance and distinguishes captures from endpoints", + "[solver][continuation_initial_state]") { + using experiment::continuation::LoadMappedSnapshot; + SnapshotFiles files; + const Problem problem; + auto mapped = LoadMappedSnapshot(problem, files.directory); + REQUIRE(mapped.Size() == 3); + CHECK(mapped(0) == 7.); CHECK(mapped(1) == 2.); CHECK(mapped(2) == 3.); + CHECK_NOTHROW(LoadMappedSnapshot(problem, files.directory / "accepted_state.txt")); + CHECK_THROWS_AS(LoadMappedSnapshot(problem, files.directory / "capture_0/physical_base.txt"), std::invalid_argument); + CHECK_NOTHROW(LoadMappedSnapshot(problem, files.directory / "capture_0/physical_base.txt", false, true)); + files.write("metadata.txt", "converged=0\n"); + CHECK_THROWS_AS(LoadMappedSnapshot(problem, files.directory), std::invalid_argument); + CHECK_NOTHROW(LoadMappedSnapshot(problem, files.directory, false, true)); + files.write("metadata.txt", "n=2\n"); + CHECK_THROWS_AS(LoadMappedSnapshot(problem, files.directory), std::invalid_argument); + files.write("metadata.txt", "converged=1\nconverged=1\n"); + CHECK_THROWS_AS(LoadMappedSnapshot(problem, files.directory), std::invalid_argument); + files.write("metadata.txt", "converged=1\n"); + for (const auto contents : {"2\n2\n3\n", "3\n2\n3\n7\nextra", "3\n2\n3\nnan\n"}) { + files.write("accepted_state.txt", contents); + CHECK_THROWS_AS(LoadMappedSnapshot(problem, files.directory), std::invalid_argument); + } + files.write("accepted_state.txt", "3\n2\n3\n7\n"); + files.write("layout.csv", "kind,block,offset,size,scale\nvalue,density,0,3,1\n"); + CHECK_THROWS_AS(LoadMappedSnapshot(problem, files.directory), std::invalid_argument); +} diff --git a/tests/solver/stellar_numerical_strategies.cpp b/tests/solver/stellar_numerical_strategies.cpp index 680fa5a..3bd2405 100644 --- a/tests/solver/stellar_numerical_strategies.cpp +++ b/tests/solver/stellar_numerical_strategies.cpp @@ -14,6 +14,7 @@ import mean_field; import test_helpers; +#include "../../experiments/continuation_snapshot.hpp" namespace stellar_numerical_strategy_test { using namespace mean_field; @@ -194,6 +195,54 @@ namespace stellar_numerical_strategy_test { } } // namespace stellar_numerical_strategy_test +TEST_CASE("Continuation snapshot transfers named coordinates and explicitly handles the central border", + "[solver][continuation_initial_state]") { + CHECK(experiment::continuation::RunSnapshotCoordinateSelfChecks() == 11); +} + +TEST_CASE("Explicit physical restart prepares a new initial guess and rolls back invalid input", + "[solver][continuation_initial_state]") { + using namespace stellar_numerical_strategy_test; + auto context = makeContext(); + using Diagnostics = solver::detail::StellarEquilibriumContextDiagnostics; + const auto original = snapshot(context); + mfem::Vector changed(original.physical); + int enthalpyOffset = -1; + Diagnostics::WithState(context, [&](auto &state, const auto &) { + for (const auto &block : state.Problem().GetManifest().valueBlocks()) { + if (block.stableId == "density") + for (int i = block.offset; i < block.offset + block.size; ++i) changed(i) *= 1.001; + if (block.stableId == "specific_enthalpy") enthalpyOffset = block.offset; + } + }); + REQUIRE(relativeDifference(original.physical, changed) > 0.); + REQUIRE_NOTHROW(Diagnostics::InitializePhysicalState(context, changed)); + const auto installed = snapshot(context); + CHECK(relativeDifference(changed, installed.physical) < 1e-14); + CHECK(relativeDifference(original.residual, installed.residual) > 1e-12); + Diagnostics::WithState(context, [&](auto &state, const auto &) { + CHECK(state.normalizedCorrection.Norml2() == 0.); + CHECK(state.physicalInverse->IsCurrent()); + CHECK(state.normalizedInverse->IsCurrent()); + }); + mfem::Vector wrongSize(changed.Size() - 1); + CHECK_THROWS_AS(Diagnostics::InitializePhysicalState(context, wrongSize), std::invalid_argument); + mfem::Vector nonfinite(changed); + nonfinite(0) = std::numeric_limits::quiet_NaN(); + CHECK_THROWS_AS(Diagnostics::InitializePhysicalState(context, nonfinite), std::invalid_argument); + REQUIRE(enthalpyOffset >= 0); + mfem::Vector inadmissible(changed); + inadmissible(enthalpyOffset) = -1e10; + CHECK_THROWS_AS(Diagnostics::InitializePhysicalState(context, inadmissible), std::invalid_argument); + checkPreserved(installed, snapshot(context)); + { + auto active = solver::make(context, nonlinear::Newton(oneAttempt())); + CHECK_THROWS_AS(Diagnostics::InitializePhysicalState(context, changed), std::logic_error); + } + REQUIRE_NOTHROW(Diagnostics::InitializePhysicalState(context, original.physical)); + checkPreserved(original, snapshot(context)); +} + TEST_CASE("Typed numerical policies preserve the existing Newton configuration", "[solver][numerical_strategies]") { using namespace stellar_numerical_strategy_test; const nonlinear::Newton original;