diff --git a/experiments/NONLINEAR_STEP_INVESTIGATION_2026-09-18.md b/experiments/NONLINEAR_STEP_INVESTIGATION_2026-09-18.md index 46181ac..0988afb 100644 --- a/experiments/NONLINEAR_STEP_INVESTIGATION_2026-09-18.md +++ b/experiments/NONLINEAR_STEP_INVESTIGATION_2026-09-18.md @@ -106,6 +106,23 @@ introduced by gravity correction; it is not identified solely with the raw physical shape Hessian. Below-floor differences are unresolved, not evidence of exact linearity. +The EOS provides an additional interpretation limit: rho(h) is proportional +to max(h,0)^1.5 here. At an exact zero with a positive enthalpy direction, +the closure remainder can scale as alpha^1.5 while its first derivative is +correct. A superlinear slope different from two does not alone establish a +Jacobian defect; the direct tangent check is whether R/alpha decreases in a +resolved regime. The strong zero-pressure boundary row is already expressed +as h-h(P_surface), with derivative delta h, so it does not inherit a vanishing +dP/dh boundary row. + +The paired directions also permit a conditioning diagnostic without another +operator solve: w = p_tight - p_loose and Jw = Jp_tight - Jp_loose. A large +change in direction with a small action demonstrates a weak response in that +direction, not a proved nullspace or a particular symmetry. Measure the +shape/fields action cancellation as well. Fixed-K, fixed-M n = 1.5 does not +have the familiar exact n = 3 mass-preserving homology; no additional gauge +or central-density constraint is justified by that analogy. + ## Execution and validation Run a late-state probe first, then use its evidence to choose the matched diff --git a/experiments/nonlinear_step_controls_2026-09-18.json b/experiments/nonlinear_step_controls_2026-09-18.json index 4fc0eea..ff08604 100644 --- a/experiments/nonlinear_step_controls_2026-09-18.json +++ b/experiments/nonlinear_step_controls_2026-09-18.json @@ -6,7 +6,6 @@ "--angular-momentum", "0", "--no-central", "--max-newton", "0", - "--strategy", "curvature", "--nonlinear-step-probe", "--probe-state", "/Users/tboudreaux/.codex/.chatgpt-projects/g-p-6a6b3f2194ec8191bebcd59075897f1d/diagnostics/preconditioner_improvement_2026-09-18/newton_v1/n1p5_J0_curvature_gravity_schur/accepted_state.txt", "--restart", "160", @@ -27,7 +26,6 @@ "--angular-momentum", "0", "--no-central", "--max-newton", "0", - "--strategy", "curvature", "--nonlinear-step-probe", "--restart", "160", "--max-linear", "600", diff --git a/experiments/nonlinear_step_late_case_2026-09-18.json b/experiments/nonlinear_step_late_case_2026-09-18.json index a52cb11..cc9b84a 100644 --- a/experiments/nonlinear_step_late_case_2026-09-18.json +++ b/experiments/nonlinear_step_late_case_2026-09-18.json @@ -6,7 +6,6 @@ "--angular-momentum", "0", "--no-central", "--max-newton", "0", - "--strategy", "curvature", "--nonlinear-step-probe", "--probe-state", "/Users/tboudreaux/.codex/.chatgpt-projects/g-p-6a6b3f2194ec8191bebcd59075897f1d/diagnostics/preconditioner_improvement_2026-09-18/newton_v1/n1p5_J0_curvature_gravity_schur/accepted_state.txt", "--restart", "160", diff --git a/experiments/plateau_diagnostic.cpp b/experiments/plateau_diagnostic.cpp index e2d7f43..c0723bf 100644 --- a/experiments/plateau_diagnostic.cpp +++ b/experiments/plateau_diagnostic.cpp @@ -23,6 +23,7 @@ import mean_field; #include "plateau_pressure_obstruction_diagnostics.hpp" #include "plateau_shape_curvature_diagnostics.hpp" #include "plateau_numerical_strategy_diagnostics.hpp" +#include "plateau_nonlinear_step_diagnostics.hpp" #include "plateau_gravity_application_diagnostics.hpp" #include "plateau_positivity_strategy_diagnostics.hpp" #include "plateau_polynomial_certificate.hpp" @@ -35,7 +36,7 @@ struct Options { int maxNewton = 12, maxLinear = 200, restart = 40, probeRestart = 40; bool probes = false, noCentral = false; bool strategyProbe = false, positivityProbe = false, positivityCertificate = false; - bool preconditionerProbe = false; + bool preconditionerProbe = false, nonlinearStepProbe = false; std::string preconditionerVariant = "independent"; int schurCorrectionSweeps = 1; double schurRelaxation = 0.5; @@ -100,6 +101,8 @@ Options Parse(int argc, char **argv) { o.positivityCertificate = true; else if (a == "--preconditioner-probe") o.preconditionerProbe = true; + else if (a == "--nonlinear-step-probe") + o.nonlinearStepProbe = true; else if (a == "--preconditioner-variant") o.preconditionerVariant = v(); else if (a == "--schur-sweeps") @@ -141,24 +144,33 @@ Options Parse(int argc, char **argv) { if (o.positivityCertificate && o.probeState.empty()) throw std::invalid_argument("--positivity-certificate requires --probe-state"); if (o.positivityCertificate && (o.probes || o.pressureProbe || o.shapeCurvature || o.surfaceTrace || - o.strategyProbe || o.positivityProbe || o.preconditionerProbe || !o.probeDirection.empty() || !o.replayCandidate.empty())) + o.strategyProbe || o.positivityProbe || o.preconditionerProbe || o.nonlinearStepProbe || !o.probeDirection.empty() || !o.replayCandidate.empty())) throw std::invalid_argument("--positivity-certificate is a certificate-only mode; omit other probe options"); if (o.preconditionerProbe) { experiment::plateau_diagnostics::ParsePreconditionerBenchmarkVariant(o.preconditionerVariant); if (o.maxNewton != 0 || o.probes || o.pressureProbe || o.shapeCurvature || o.surfaceTrace || - o.strategyProbe || o.positivityProbe || o.positivityCertificate || !o.probeDirection.empty() || + o.strategyProbe || o.positivityProbe || o.positivityCertificate || o.nonlinearStepProbe || !o.probeDirection.empty() || !o.replayCandidate.empty() || !o.captureIterations.empty() || o.strategy != "baseline") throw std::invalid_argument("--preconditioner-probe requires --max-newton 0 and no nonlinear/probe strategies"); if (!std::isfinite(o.linearTolerance) || o.linearTolerance <= 0.0) throw std::invalid_argument("A preconditioner benchmark requires a finite positive --linear-tol"); } + if (o.nonlinearStepProbe) { + if (o.maxNewton != 0 || o.probes || o.pressureProbe || o.shapeCurvature || o.surfaceTrace || + o.strategyProbe || o.positivityProbe || o.positivityCertificate || o.preconditionerProbe || + !o.probeDirection.empty() || !o.replayCandidate.empty() || !o.captureIterations.empty() || + o.strategy != "baseline") + throw std::invalid_argument("--nonlinear-step-probe requires --max-newton 0 and no other probes or strategy; it compares raw and curvature paths itself"); + if (!std::isfinite(o.linearTolerance) || o.linearTolerance <= 0.0) + throw std::invalid_argument("A nonlinear-step probe requires a finite positive --linear-tol"); + } if (!o.preconditionerProbe) { if (o.preconditionerVariant != "independent" && o.preconditionerVariant != "gravity-schur" && o.preconditionerVariant != "stellar-ldu" && o.preconditionerVariant != "stellar-ldu-gravity-schur") throw std::invalid_argument("Newton preconditioner variants are independent, gravity-schur, stellar-ldu, or stellar-ldu-gravity-schur"); if (o.preconditionerVariant != "independent") { if (o.probes || o.pressureProbe || o.shapeCurvature || o.surfaceTrace || o.strategyProbe || - o.positivityProbe || o.positivityCertificate || !o.probeState.empty() || !o.probeDirection.empty() || + o.positivityProbe || o.positivityCertificate || (!o.probeState.empty() && !o.nonlinearStepProbe) || !o.probeDirection.empty() || !o.replayCandidate.empty()) throw std::invalid_argument("A custom Newton preconditioner cannot be combined with frozen or certificate probes"); if (o.strategy != "baseline" && o.strategy != "positivity" && o.strategy != "curvature" && @@ -441,6 +453,7 @@ int RunContext(Context &context, const Options &o, const Scales &scales) { << "\nM=1\nR_reference=1\nxi1=" << scales.xi << "\nmu1=" << scales.mu << "\ncentral_phase=" << !o.noCentral << "\nrestart=" << o.restart << "\nprobe_restart=" << o.probeRestart << "\nfrozen_probe_backend=" << (o.preconditionerProbe ? "direct_mfem_fgmres_preconditioner_benchmark" : + o.nonlinearStepProbe ? "prepared_fgmres_cold_with_independent_true_gate" : o.positivityCertificate ? "polynomial_certificate_no_linear_solve" : o.probeDirection.empty() ? "direct_mfem_fgmres_cold" : "supplied_direction_no_linear_solve") << "\nlinear_tol=" << o.linearTolerance << "\nmax_linear=" << o.maxLinear @@ -448,6 +461,7 @@ int RunContext(Context &context, const Options &o, const Scales &scales) { << "\nstrategy=" << o.strategy << "\ngravity_tolerance=" << o.gravityTolerance << "\ngravity_iterations=" << o.gravityIterations << "\ngravity_restart=" << o.gravityRestart << "\npreconditioner_probe=" << o.preconditionerProbe + << "\nnonlinear_step_probe=" << o.nonlinearStepProbe << "\npreconditioner_variant=" << o.preconditionerVariant << "\ncontext_preconditioner_variant=" << (o.preconditionerProbe ? "independent" : o.preconditionerVariant) << "\ncontext_gravity_schur_sweeps=" @@ -456,6 +470,27 @@ int RunContext(Context &context, const Options &o, const Scales &scales) { << "\nschur_correction_sweeps=" << o.schurCorrectionSweeps << "\nschur_relaxation=" << o.schurRelaxation << "\ncompiler=" << __VERSION__ << "\ncompiled=" << __DATE__ << ' ' << __TIME__ << '\n'; metadata.flush(); + if (o.nonlinearStepProbe) { + metadata << "nonlinear_step_input=" << (o.probeState.empty() ? "lane_emden_seed" : "physical_snapshot") + << "\nnonlinear_step_paths=raw,curvature\nnonlinear_step_positivity_repair=0\n"; + if (!o.probeState.empty()) metadata << "probe_state=" << o.probeState << '\n'; + metadata.flush(); + solver::detail::StellarEquilibriumContextDiagnostics::WithState(context, [&](auto &s, const fem::FEM &) { + mfem::Vector base(s.acceptedNormalizedState); + if (!o.probeState.empty()) { + const auto snapshot = LoadVector(o.probeState, s.AcceptedPhysicalState().Size()); + s.normalizedOperator->NormalizeState(snapshot, base); + } + experiment::plateau_diagnostics::ProbeNonlinearStep(s, base, o.output, { + .linear = {.relativeTolerance = o.linearTolerance, .absoluteTolerance = 0., + .maximumIterations = o.maxLinear}, + .restartLength = o.restart, + .gravity = {.control = {.relativeTolerance = o.gravityTolerance, .absoluteTolerance = 1e-14, + .maximumIterations = o.gravityIterations}, + .backend = {.restartLength = o.gravityRestart, .printLevel = -1}}}); + }); + return 0; + } // 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'; diff --git a/experiments/plateau_nonlinear_step_diagnostics.hpp b/experiments/plateau_nonlinear_step_diagnostics.hpp new file mode 100644 index 0000000..622bbb5 --- /dev/null +++ b/experiments/plateau_nonlinear_step_diagnostics.hpp @@ -0,0 +1,590 @@ +#pragma once + +// Experiment-only: include after import mean_field. The caller holds the state +// lease. A physical snapshot uses the context's frozen normalization throughout. +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace experiment::plateau_diagnostics { +struct NonlinearStepProbeOptions final { + mean_field::solver::LinearSolveControl linear{.relativeTolerance = .03, .absoluteTolerance = 0., + .maximumIterations = 200}; + int restartLength{40}; + mean_field::solver::nonlinear::GravitySolveOptions gravity{}; +}; + +namespace nonlinear_step_detail { +using Clock = std::chrono::steady_clock; +inline double Seconds(Clock::time_point start) { + return std::chrono::duration(Clock::now() - start).count(); +} +inline bool Finite(const mfem::Vector &v) { + for (int i = 0; i < v.Size(); ++i) if (!std::isfinite(v(i))) return false; + return true; +} +inline bool Identical(const mfem::Vector &a, const mfem::Vector &b) { + return a.Size() == b.Size() && + (a.Size() == 0 || std::memcmp(a.GetData(), b.GetData(), sizeof(mfem::real_t) * a.Size()) == 0); +} +inline long double Dot(const mfem::Vector &a, const mfem::Vector &b, int offset, int size) { + long double result = 0.; + for (int i = offset; i < offset + size; ++i) result += static_cast(a(i)) * b(i); + return result; +} +inline double Norm(const mfem::Vector &v, int offset, int size) { + return std::sqrt(static_cast(Dot(v, v, offset, size))); +} +inline double Norm(const mfem::Vector &v) { return Norm(v, 0, v.Size()); } +inline double Difference(const mfem::Vector &a, const mfem::Vector &b) { + mfem::Vector difference(a); + difference -= b; + return Norm(difference); +} +inline double Ratio(double numerator, double denominator) { + return std::isfinite(numerator) && std::isfinite(denominator) && denominator > 0. + ? numerator / denominator : std::numeric_limits::quiet_NaN(); +} +inline std::ofstream Open(const std::filesystem::path &path) { + if (std::filesystem::exists(path)) + throw std::invalid_argument("Nonlinear-step output already exists: " + path.string()); + std::ofstream file(path); + file.exceptions(std::ios::failbit | std::ios::badbit); + file << std::setprecision(17); + return file; +} +inline void Save(const std::filesystem::path &path, const mfem::Vector &v) { + auto file = Open(path); + file << v.Size() << '\n'; + for (int i = 0; i < v.Size(); ++i) file << v(i) << '\n'; +} +struct Observation final { + bool prepared{false}, innerConverged{true}, preparedGate{true}, usable{false}; + mfem::Vector residual, remainder; + double correctionNorm{std::numeric_limits::quiet_NaN()}; + double targetDefect{std::numeric_limits::quiet_NaN()}; + double affineDisagreement{std::numeric_limits::quiet_NaN()}; + double arithmeticAllowance{std::numeric_limits::quiet_NaN()}; + double targetThreshold{std::numeric_limits::quiet_NaN()}; + double seconds{0.}; +}; +} // namespace nonlinear_step_detail + +// The SOC target is G0 + alpha J_G p, not zero gravity residual. Its retained +// rows may have substantial curvature even when that target is met accurately. +// Full and split vectors let independent analysis test this trial map's tangent +// without mistaking scalar norm cancellation for a correct vector derivative. +template +void ProbeNonlinearStep(State &state, const mfem::Vector &normalizedBase, + const std::filesystem::path &output, NonlinearStepProbeOptions options = {}) { + using namespace mean_field; + namespace d = nonlinear_step_detail; + using Form = typename std::remove_cvref_t::FormType; + using namespace utils::blocks; + constexpr value_block> shape{}; + constexpr value_block> gradient{}; + constexpr value_block> potential{}; + options.linear.Validate(); + options.gravity.Validate(); + int ranks = 0; + MPI_Comm_size(state.Problem().GetCommunicator(), &ranks); + const int size = state.Problem().StateSize(); + if (ranks != 1 || size != state.Problem().EquationSize() || normalizedBase.Size() != size || + !d::Finite(normalizedBase) || options.restartLength < 1) + throw std::invalid_argument("Nonlinear-step probes require one rank, compatible finite vectors and valid controls."); + const auto start = d::Clock::now(); + const double nan = std::numeric_limits::quiet_NaN(); + const double eps = std::numeric_limits::epsilon(); + const mfem::Vector base(normalizedBase); + const auto &normalization = state.normalizedOperator->GetNormalization(); + const mfem::Vector stateFactors(normalization.StateFactors()), residualFactors(normalization.ResidualFactors()); + const mfem::Vector savedAcceptedState(state.acceptedNormalizedState), + savedAcceptedResidual(state.acceptedNormalizedResidual), savedAcceptedPhysical(state.AcceptedPhysicalState()); + const mfem::Vector savedCorrection(state.normalizedCorrection), savedPhysicalCorrection(state.physicalCorrection), + savedVolumeDirection(state.volumeDisplacementDirection); + const auto &layout = state.Problem().GetManifest().layout(); + const auto rows = state.Problem().GetManifest().residualBlocks(); + const int shapeOffset = layout.offset(shape), shapeSize = layout.size(shape); + if (shapeSize <= 0) throw std::invalid_argument("The nonlinear-step probe requires a nonempty shape block."); + constexpr std::array names{"full", "shape", "fields"}; + constexpr std::array alphas{1., .25, .0625, .015625, .00390625, .0009765625, + .000244140625, .00006103515625}; + std::filesystem::create_directories(output); + const auto file = [&](const std::string &name) { return output / ("nonlinear_step_" + name); }; + auto summary = d::Open(file("summary.csv")); + auto linear = d::Open(file("linear.csv")); + auto preparations = d::Open(file("preparations.csv")); + auto paths = d::Open(file("paths.csv")); + auto blocks = d::Open(file("blocks.csv")); + auto mixed = d::Open(file("mixed.csv")); + auto inner = d::Open(file("inner.csv")); + auto restorations = d::Open(file("restorations.csv")); + summary << "metric,value\n"; + linear << "status,converged,iterations,rhs_norm,reported_true_norm,verified_true_norm,verified_relative_true," + "threshold,verified_gate,operator_applications,inverse_applications,seconds\n"; + preparations << "sequence,kind,direction,alpha,prepared,rejection_reason,rejection_stage,seconds,residual_norm\n"; + paths << "sample,alpha,direction,path,prepared,inner_converged,prepared_gate,usable,base_norm,action_norm," + "residual_norm,merit,slope,linear_prediction_norm,remainder_norm,remainder_over_alpha," + "remainder_over_alpha_squared,relative_derivative_error,linear_dot_remainder,norm_identity_error," + "expected_step_norm,actual_step_norm,addition_error_norm,increment_resolved,gravity_correction_norm," + "gravity_target_defect,gravity_affine_disagreement,gravity_arithmetic_allowance,gravity_target_threshold," + "seconds,heuristic_roundoff_scale\n"; + blocks << "sample,alpha,direction,path,block,offset,size,base_norm,action_norm,linear_prediction_norm,residual_norm," + "remainder_norm,remainder_over_alpha,remainder_over_alpha_squared,linear_dot_remainder," + "norm_identity_error,heuristic_roundoff_scale\n"; + mixed << "sample,alpha,path,block,offset,size,usable,full_remainder_norm,shape_remainder_norm,fields_remainder_norm," + "mixed_remainder_norm,linear_dot_full,linear_dot_shape,linear_dot_fields,linear_dot_mixed," + "projection_identity_error,norm_identity_error\n"; + inner << "sample,alpha,direction,status,converged,iterations,rhs_norm,true_norm,relative_true,threshold," + "operator_applications,inverse_applications,seconds,compact_operator,compact_inverse\n"; + restorations << "sequence,kind,block,offset,size,base_norm,difference_norm\n"; + auto metric = [&](const char *name, auto value) { summary << name << ',' << value << '\n'; }; + metric("schema_version", 1); + metric("linear_relative_tolerance", options.linear.relativeTolerance); + metric("linear_absolute_tolerance", options.linear.absoluteTolerance); + metric("linear_maximum_iterations", options.linear.maximumIterations); + metric("linear_restart_length", options.restartLength); + metric("gravity_relative_tolerance", options.gravity.control.relativeTolerance); + metric("gravity_absolute_tolerance", options.gravity.control.absoluteTolerance); + metric("gravity_maximum_iterations", options.gravity.control.maximumIterations); + metric("gravity_restart_length", options.gravity.backend.restartLength); + metric("shape_offset", shapeOffset); + metric("shape_size", shapeSize); + metric("machine_epsilon", eps); + metric("polynomial_direction_repairs", 0); + metric("accepted_commits", 0); + metric("linear_gate_failure_skips_trials", 0); + metric("path_seconds_include_shared_raw_preparation", 1); + metric("roundoff_scale_is_subtraction_heuristic", 1); + for (std::size_t sample = 0; sample < alphas.size(); ++sample) + summary << "alpha_" << sample << ',' << alphas[sample] << '\n'; + mfem::Vector originalResidual(size); + bool haveResidual = false; + int preparationCount = 0, refreshCount = 0, directFullActions = 0, captureFullActions = 0; + int assessmentGravityActions = 0, successfulPreparedPaths = 0; + double preparationSeconds = 0., refreshSeconds = 0., directActionSeconds = 0., captureSeconds = 0.; + double assessmentGravitySeconds = 0.; + double restorationError = 0.; + auto eachBlock = [&](auto &&operation) { + operation("all", 0, size); + for (const auto &row : rows) operation(row.stableId, row.offset, row.size); + }; + auto prepare = [&](const mfem::Vector &candidate, const char *kind, const char *direction, double alpha) { + state.trialNormalizedState = candidate; + const auto preparationStart = d::Clock::now(); + const auto result = state.PrepareTrial(); + const double elapsed = d::Seconds(preparationStart); + preparationSeconds += elapsed; + ++preparationCount; + preparations << preparationCount << ',' << kind << ',' << direction << ',' << alpha << ',' << bool(result) + << ',' << (result ? -1 : int(result.error().reason)) << ',' + << (result ? -1 : int(result.error().stage)) << ',' << elapsed << ',' + << (result ? d::Norm(state.trialNormalizedResidual) : nan) << '\n'; + preparations.flush(); + if (result && !d::Finite(state.trialNormalizedResidual)) + throw std::runtime_error("Nonlinear-step preparation produced a nonfinite residual."); + return bool(result); + }; + auto refresh = [&]() { + const auto refreshStart = d::Clock::now(); + state.physicalInverse->Refresh(); + refreshSeconds += d::Seconds(refreshStart); + ++refreshCount; + }; + auto restore = [&](const char *kind) { + if (!prepare(base, kind, "base", 0.)) + throw std::runtime_error("The nonlinear-step probe could not restore its base."); + if (!haveResidual) { originalResidual = state.trialNormalizedResidual; haveResidual = true; } + mfem::Vector difference(state.trialNormalizedResidual); + difference -= originalResidual; + restorationError = std::max(restorationError, d::Norm(difference)); + eachBlock([&](auto name, int offset, int count) { + restorations << preparationCount << ',' << kind << ',' << name << ',' << offset << ',' << count << ',' + << d::Norm(originalResidual, offset, count) << ',' << d::Norm(difference, offset, count) << '\n'; + }); + restorations.flush(); + }; + auto action = [&](const mfem::Vector &direction, mfem::Vector &image) { + const auto actionStart = d::Clock::now(); + state.normalizedOperator->Mult(direction, image); + directActionSeconds += d::Seconds(actionStart); + ++directFullActions; + if (!d::Finite(image)) throw std::runtime_error("Nonlinear-step J action is nonfinite."); + }; + auto finish = [&]() { + restore("final_restore"); + d::Save(file("restored_residual.txt"), state.trialNormalizedResidual); + refresh(); + state.normalizedCorrection = savedCorrection; + state.physicalCorrection = savedPhysicalCorrection; + state.volumeDisplacementDirection = savedVolumeDirection; + const bool acceptedUnchanged = d::Identical(savedAcceptedState, state.acceptedNormalizedState) && + d::Identical(savedAcceptedResidual, state.acceptedNormalizedResidual) && + d::Identical(savedAcceptedPhysical, state.AcceptedPhysicalState()); + const bool factorsUnchanged = d::Identical(stateFactors, normalization.StateFactors()) && + d::Identical(residualFactors, normalization.ResidualFactors()); + metric("maximum_original_restoration_error", restorationError); + metric("accepted_state_change", d::Difference(savedAcceptedState, state.acceptedNormalizedState)); + metric("accepted_residual_change", d::Difference(savedAcceptedResidual, state.acceptedNormalizedResidual)); + metric("accepted_physical_change", d::Difference(savedAcceptedPhysical, state.AcceptedPhysicalState())); + metric("accepted_buffers_bitwise_identical", acceptedUnchanged); + metric("normalization_factors_bitwise_identical", factorsUnchanged); + metric("preparations", preparationCount); + metric("preparation_seconds", preparationSeconds); + metric("inverse_refreshes", refreshCount); + metric("inverse_refresh_seconds", refreshSeconds); + metric("direct_full_operator_applications", directFullActions); + metric("direct_full_operator_seconds", directActionSeconds); + metric("capture_full_operator_applications", captureFullActions); + metric("capture_seconds", captureSeconds); + metric("assessment_gravity_operator_applications", assessmentGravityActions); + metric("assessment_gravity_operator_seconds", assessmentGravitySeconds); + metric("prepared_paths", successfulPreparedPaths); + metric("total_seconds", d::Seconds(start)); + summary.flush(); + if (!acceptedUnchanged || !factorsUnchanged) + throw std::runtime_error("The nonlinear-step probe changed accepted buffers or normalization."); + }; + + try { + restore("initial_restore"); + refresh(); + mfem::Vector physicalBase; + state.normalizedOperator->DenormalizeState(base, physicalBase); + d::Save(file("physical_base.txt"), physicalBase); + d::Save(file("normalized_base.txt"), base); + d::Save(file("residual.txt"), originalResidual); + d::Save(file("state_factors.txt"), stateFactors); + d::Save(file("residual_factors.txt"), residualFactors); + auto manifest = d::Open(file("layout.csv")); + manifest << "kind,block,offset,size\n"; + for (const auto &value : state.Problem().GetManifest().valueBlocks()) + manifest << "value," << value.stableId << ',' << value.offset << ',' << value.size << '\n'; + for (const auto &row : rows) + manifest << "residual," << row.stableId << ',' << row.offset << ',' << row.size << '\n'; + mfem::Vector rhs(originalResidual); + rhs *= -1.; + d::Save(file("rhs.txt"), rhs); + std::array directions{mfem::Vector(size), mfem::Vector(size), mfem::Vector(size)}; + std::array actions{mfem::Vector(size), mfem::Vector(size), mfem::Vector(size)}; + directions[0] = 0.; + auto backend = solver::linear::prepareLinearBackend( + solver::linear::FGMRES({.restartLength = options.restartLength, .printLevel = -1}), + *state.normalizedOperator, *state.normalizedInverse, state.Problem().GetCommunicator()); + const auto solve = backend.Solve(rhs, directions[0], options.linear); + // Reprepare before the independent certificate: this also tests that a + // solve has not changed the operator's promised frozen linearization. + restore("verify_linear_base"); + action(directions[0], actions[0]); + mfem::Vector defect(actions[0]); + defect -= rhs; + const double trueNorm = d::Norm(defect), rhsNorm = d::Norm(rhs); + const double threshold = options.linear.ConvergenceThreshold(rhsNorm); + const bool trueGate = d::Finite(directions[0]) && std::isfinite(trueNorm) && trueNorm <= threshold; + linear << int(solve.status) << ',' << solve.Converged() << ',' << solve.iterations << ',' << rhsNorm << ',' + << solve.trueResidualNorm << ',' << trueNorm << ',' << d::Ratio(trueNorm, rhsNorm) << ',' << threshold + << ',' << trueGate << ',' << solve.operatorApplications << ',' + << solve.inversePreconditionerApplications << ',' << solve.solveSeconds << '\n'; + linear.flush(); + metric("linear_verified_gate", trueGate); + metric("linear_solve_operator_applications", solve.operatorApplications); + metric("linear_solve_inverse_applications", solve.inversePreconditionerApplications); + metric("linear_solve_seconds", solve.solveSeconds); + metric("linear_solve_operator_seconds", solve.operatorSeconds); + metric("linear_solve_inverse_seconds", solve.inversePreconditionerSeconds); + if (!d::Finite(directions[0])) throw std::runtime_error("A frozen direction is nonfinite."); + directions[1] = 0.; + directions[2] = directions[0]; + for (int i = shapeOffset; i < shapeOffset + shapeSize; ++i) { + directions[1](i) = directions[0](i); + directions[2](i) = 0.; + } + for (int k = 1; k < 3; ++k) { restore("split_action_base"); action(directions[k], actions[k]); } + mfem::Vector splitError(actions[0]); + splitError -= actions[1]; + splitError -= actions[2]; + metric("action_split_error_norm", d::Norm(splitError)); + metric("base_residual_norm", d::Norm(originalResidual)); + metric("full_direction_norm", d::Norm(directions[0])); + metric("full_merit_slope", d::Dot(originalResidual, actions[0], 0, size)); + for (int k = 0; k < 3; ++k) { + d::Save(file(std::string("direction_") + names[k] + ".txt"), directions[k]); + d::Save(file(std::string("action_") + names[k] + ".txt"), actions[k]); + } + restore("gravity_setup_base"); + refresh(); + auto partition = solver::MakeGravityBlockPartition
(layout); + solver::PreparedGravitySubsystem gravitySubsystem(*state.normalizedOperator, *state.normalizedInverse, + partition, state.Problem().GetCommunicator(), solver::linear::FGMRES(options.gravity.backend)); + solver::PreparedGravityCurvatureCorrection correction(gravitySubsystem); + metric("compact_gravity_operator", gravitySubsystem.UsesCompactOperator()); + metric("compact_gravity_inverse", gravitySubsystem.UsesCompactPreconditioner()); + mfem::Vector compactCorrection(partition.GravityValueSize()), gravityAction(partition.GravityResidualSize()); + mfem::Vector target(partition.GravityResidualSize()), beforeGravity(partition.GravityResidualSize()); + mfem::Vector source(partition.GravityResidualSize()), compactValues(partition.GravityValueSize()); + mfem::Vector scaleAction(partition.GravityResidualSize()), rootValues(size), fullCorrection(size); + std::uint64_t innerOperatorApplications = 0, innerInverseApplications = 0; + double innerOperatorSeconds = 0., innerInverseSeconds = 0.; + const auto assessmentAction = [&](const mfem::Vector &values, mfem::Vector &image) { + const auto actionStart = d::Clock::now(); + gravitySubsystem.GetOperator().Mult(values, image); + assessmentGravitySeconds += d::Seconds(actionStart); + ++assessmentGravityActions; + }; + + for (std::size_t sample = 0; sample < alphas.size(); ++sample) { + const double alpha = alphas[sample]; + std::array, 2> observations; + for (int k = 0; k < 3; ++k) { + restore("capture_base"); + refresh(); + gravitySubsystem.Refresh(); + const auto captureStart = d::Clock::now(); + correction.CaptureBase(originalResidual, directions[k]); + captureSeconds += d::Seconds(captureStart); + ++captureFullActions; + mfem::Vector candidate(base), predicted(originalResidual), expectedIncrement(directions[k]); + candidate.Add(alpha, directions[k]); + predicted.Add(alpha, actions[k]); + d::Save(file(std::string("prediction_") + names[k] + "_a" + std::to_string(sample) + ".txt"), predicted); + expectedIncrement *= alpha; + mfem::Vector actualIncrement(candidate); + actualIncrement -= base; + auto &raw = observations[0][k]; + auto &curved = observations[1][k]; + const auto rawStart = d::Clock::now(); + raw.prepared = prepare(candidate, "raw", names[k], alpha); + raw.seconds = d::Seconds(rawStart); + raw.usable = raw.prepared; + if (raw.prepared) { + raw.residual = state.trialNormalizedResidual; + raw.remainder = raw.residual; + raw.remainder -= predicted; + } + const auto correctedStart = d::Clock::now(); + curved.innerConverged = curved.preparedGate = false; + if (raw.prepared) { + refresh(); + gravitySubsystem.Refresh(); + correction.BuildLinearPrediction(alpha, target); + const auto report = correction.SolveCorrection(alpha, raw.residual, compactCorrection, + options.gravity.control); + innerOperatorApplications += report.operatorApplications; + innerInverseApplications += report.inversePreconditionerApplications; + innerOperatorSeconds += report.operatorSeconds; + innerInverseSeconds += report.inversePreconditionerSeconds; + inner << sample << ',' << alpha << ',' << names[k] << ',' << int(report.status) << ',' + << report.Converged() << ',' << report.iterations << ',' << report.rightHandSideNorm << ',' + << report.trueResidualNorm << ',' << report.relativeTrueResidualNorm << ',' + << (std::isfinite(report.rightHandSideNorm) && report.rightHandSideNorm >= 0. + ? options.gravity.control.ConvergenceThreshold(report.rightHandSideNorm) : nan) << ',' + << report.operatorApplications << ',' << report.inversePreconditionerApplications << ',' + << report.solveSeconds << ',' << gravitySubsystem.UsesCompactOperator() << ',' + << gravitySubsystem.UsesCompactPreconditioner() << '\n'; + inner.flush(); + curved.innerConverged = report.Converged(); + if (report.Converged()) { + partition.InjectGravityValues(compactCorrection, fullCorrection); + d::Save(file(std::string("correction_curvature_") + names[k] + "_a" + + std::to_string(sample) + ".txt"), fullCorrection); + assessmentAction(compactCorrection, gravityAction); + partition.GatherGravityResiduals(raw.residual, beforeGravity); + source = beforeGravity; + std::array contributions{}; + // Match production ApplyCorrection's arithmetic guard: + // separate the large cancelling g/Phi terms, rather than + // estimating preparation error from a tiny gravity defect. + const auto measure = [&](const mfem::Vector &values, auto block, bool subtractSource) { + rootValues = 0.; + for (int i = layout.offset(block); i < layout.offset(block) + layout.size(block); ++i) + rootValues(i) = values(i); + partition.GatherGravityValues(rootValues, compactValues); + assessmentAction(compactValues, scaleAction); + if (subtractSource) source -= scaleAction; + return d::Norm(scaleAction); + }; + contributions[0] = measure(candidate, gradient, true); + contributions[1] = measure(candidate, potential, true); + contributions[2] = d::Norm(source); + contributions[3] = measure(fullCorrection, gradient, false); + contributions[4] = measure(fullCorrection, potential, false); + contributions[5] = d::Norm(beforeGravity); + contributions[6] = d::Norm(target); + curved.correctionNorm = d::Norm(compactCorrection); + partition.AddGravityValues(compactCorrection, candidate); + // Recompute the full prepared residual after correction; + // algebraic A delta alone cannot certify this trial map. + state.normalizedOperator->DenormalizeState(candidate, state.candidatePhysicalState); + curved.prepared = prepare(candidate, "curvature", names[k], alpha); + if (curved.prepared) { + curved.residual = state.trialNormalizedResidual; + curved.remainder = curved.residual; + curved.remainder -= predicted; + mfem::Vector actualGravity(partition.GravityResidualSize()); + partition.GatherGravityResiduals(curved.residual, actualGravity); + contributions[7] = d::Norm(actualGravity); + mfem::Vector affine(beforeGravity); + affine += gravityAction; + affine -= actualGravity; + actualGravity -= target; + curved.targetDefect = d::Norm(actualGravity); + curved.affineDisagreement = d::Norm(affine); + const auto assessment = solver::detail::AssessGravityPreparedResidual( + contributions, curved.affineDisagreement, curved.targetDefect, report.trueResidualNorm); + curved.preparedGate = assessment.admissible; + curved.arithmeticAllowance = assessment.arithmeticAllowance; + curved.targetThreshold = assessment.targetThreshold; + curved.usable = curved.preparedGate; + } + } + } + // Curvature cost includes its shared raw preparation. Summing + // path seconds double-counts that work; total_seconds does not. + curved.seconds = raw.seconds + d::Seconds(correctedStart); + for (int path = 0; path < 2; ++path) { + const auto &observation = observations[path][k]; + const char *pathName = path == 0 ? "raw" : "curvature"; + mfem::Vector intendedIncrement(expectedIncrement), pathIncrement(actualIncrement); + if (path == 1 && observation.innerConverged) { + intendedIncrement += fullCorrection; + pathIncrement = candidate; + pathIncrement -= base; + } + mfem::Vector pathAdditionError(pathIncrement); + pathAdditionError -= intendedIncrement; + const bool haveIncrement = path == 0 || observation.innerConverged; + const double expectedPathStep = haveIncrement ? d::Norm(intendedIncrement) : nan; + const double actualPathStep = haveIncrement ? d::Norm(pathIncrement) : nan; + const bool pathResolved = haveIncrement && expectedPathStep > 64. * eps * (1. + d::Norm(base)) && + actualPathStep > 0. && d::Norm(pathAdditionError) <= .01 * expectedPathStep; + const double residualNorm = observation.prepared ? d::Norm(observation.residual) : nan; + const double remainderNorm = observation.prepared ? d::Norm(observation.remainder) : nan; + const double slope = static_cast(d::Dot(originalResidual, actions[k], 0, size)); + const double projection = observation.prepared + ? static_cast(d::Dot(predicted, observation.remainder, 0, size)) : nan; + const auto identityError = [&](int offset, int count) { + return observation.prepared ? static_cast( + d::Dot(observation.residual, observation.residual, offset, count) - + d::Dot(predicted, predicted, offset, count) - + 2. * d::Dot(predicted, observation.remainder, offset, count) - + d::Dot(observation.remainder, observation.remainder, offset, count)) : nan; + }; + const auto arithmeticFloor = [&](int offset, int count) { + // Only a subtraction heuristic; production's gravity + // gate separately measures assembled cancellation scales. + return observation.prepared ? std::max(8. * restorationError, 128. * eps * ( + 1. + d::Norm(originalResidual, offset, count) + alpha * d::Norm(actions[k], offset, count) + + d::Norm(observation.residual, offset, count))) : nan; + }; + paths << sample << ',' << alpha << ',' << names[k] << ',' << pathName << ',' + << observation.prepared << ',' << observation.innerConverged << ',' + << observation.preparedGate << ',' << observation.usable << ',' << d::Norm(originalResidual) + << ',' << d::Norm(actions[k]) << ',' << residualNorm << ',' << .5 * residualNorm * residualNorm + << ',' << slope << ',' << d::Norm(predicted) << ',' << remainderNorm << ',' + << remainderNorm / alpha << ',' << remainderNorm / (alpha * alpha) << ',' + << d::Ratio(remainderNorm / alpha, d::Norm(actions[k])) << ',' << projection << ',' + << identityError(0, size) << ',' << expectedPathStep << ',' << actualPathStep << ',' + << (haveIncrement ? d::Norm(pathAdditionError) : nan) << ',' << pathResolved << ',' << observation.correctionNorm << ',' + << observation.targetDefect << ',' << observation.affineDisagreement << ',' + << observation.arithmeticAllowance << ',' << observation.targetThreshold << ',' + << observation.seconds << ',' << arithmeticFloor(0, size) << '\n'; + if (observation.prepared) { + ++successfulPreparedPaths; + const std::string suffix = std::string(pathName) + '_' + names[k] + "_a" + std::to_string(sample) + ".txt"; + d::Save(file("residual_" + suffix), observation.residual); + // The raw state must not accidentally inherit the SOC + // update from the shared trial buffer. + mfem::Vector recordedState(base); + recordedState.Add(alpha, directions[k]); + if (path == 1) recordedState = candidate; + d::Save(file("state_" + suffix), recordedState); + } + eachBlock([&](auto name, int offset, int count) { + const double norm = observation.prepared ? d::Norm(observation.remainder, offset, count) : nan; + blocks << sample << ',' << alpha << ',' << names[k] << ',' << pathName << ',' << name << ',' + << offset << ',' << count << ',' << d::Norm(originalResidual, offset, count) << ',' + << d::Norm(actions[k], offset, count) << ',' << d::Norm(predicted, offset, count) << ',' + << (observation.prepared ? d::Norm(observation.residual, offset, count) : nan) << ',' + << norm << ',' << norm / alpha << ',' << norm / (alpha * alpha) << ',' + << (observation.prepared ? static_cast(d::Dot(predicted, observation.remainder, offset, count)) : nan) + << ',' << identityError(offset, count) << ',' << arithmeticFloor(offset, count) << '\n'; + }); + } + paths.flush(); blocks.flush(); + } + // For corrected trials this is an interaction of the corrected + // trial map, not solely the Hessian of the original affine path. + mfem::Vector fullPrediction(originalResidual); + fullPrediction.Add(alpha, actions[0]); + for (int path = 0; path < 2; ++path) { + const auto &o = observations[path]; + const bool allPrepared = o[0].prepared && o[1].prepared && o[2].prepared; + const bool usable = o[0].usable && o[1].usable && o[2].usable; + mfem::Vector interaction(size); + if (allPrepared) { + interaction = o[0].remainder; + interaction -= o[1].remainder; + interaction -= o[2].remainder; + } + eachBlock([&](auto name, int offset, int count) { + std::array norms{nan, nan, nan, nan}; + std::array projections{nan, nan, nan, nan}; + double identity = nan; + if (allPrepared) { + for (int k = 0; k < 3; ++k) { + norms[k] = d::Norm(o[k].remainder, offset, count); + projections[k] = d::Dot(fullPrediction, o[k].remainder, offset, count); + } + norms[3] = d::Norm(interaction, offset, count); + projections[3] = d::Dot(fullPrediction, interaction, offset, count); + identity = static_cast(d::Dot(o[0].residual, o[0].residual, offset, count) - + d::Dot(fullPrediction, fullPrediction, offset, count) - + 2. * (projections[1] + projections[2] + projections[3]) - + d::Dot(o[0].remainder, o[0].remainder, offset, count)); + } + mixed << sample << ',' << alpha << ',' << (path == 0 ? "raw" : "curvature") << ',' << name << ',' + << offset << ',' << count << ',' << usable; + for (double norm : norms) mixed << ',' << norm; + for (long double projection : projections) mixed << ',' << projection; + mixed << ',' << projections[0] - projections[1] - projections[2] - projections[3] << ',' << identity << '\n'; + }); + } + mixed.flush(); summary.flush(); + } + metric("gravity_solves", gravitySubsystem.GetSolveCount()); + metric("gravity_iterations", gravitySubsystem.GetIterationCount()); + metric("gravity_failures", gravitySubsystem.GetFailureCount()); + metric("gravity_solve_seconds", gravitySubsystem.GetSolveSeconds()); + metric("gravity_solve_operator_applications", innerOperatorApplications); + metric("gravity_solve_inverse_applications", innerInverseApplications); + metric("gravity_solve_operator_seconds", innerOperatorSeconds); + metric("gravity_solve_inverse_seconds", innerInverseSeconds); + double maximumActionChange = 0.; + for (int k = 0; k < 3; ++k) { + restore("final_split_action_base"); + mfem::Vector restoredAction(size); + action(directions[k], restoredAction); + maximumActionChange = std::max(maximumActionChange, d::Difference(actions[k], restoredAction)); + d::Save(file(std::string("restored_action_") + names[k] + ".txt"), restoredAction); + } + metric("maximum_restored_action_difference", maximumActionChange); + } catch (...) { + const auto failure = std::current_exception(); + finish(); + std::rethrow_exception(failure); + } + finish(); +} +} // namespace experiment::plateau_diagnostics