From fd08f7b0e935883a77e2c85c88d26efd31f8bb04 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Thu, 17 Sep 2026 08:06:37 -0400 Subject: [PATCH] experiment: capture frozen live Newton trials without changing solver policy --- .../plateau_capture_cases_2026-09-17.json | 36 +++++++++++ experiments/plateau_diagnostic.cpp | 61 +++++++++++++++++++ experiments/run_plateau_campaign.py | 8 +-- 3 files changed, 101 insertions(+), 4 deletions(-) create mode 100644 experiments/plateau_capture_cases_2026-09-17.json diff --git a/experiments/plateau_capture_cases_2026-09-17.json b/experiments/plateau_capture_cases_2026-09-17.json new file mode 100644 index 0000000..0f84fb8 --- /dev/null +++ b/experiments/plateau_capture_cases_2026-09-17.json @@ -0,0 +1,36 @@ +[ + { + "id": "n1_J0p10_capture", + "arguments": [ + "--index", + "1", + "--angular-momentum", + "0.1", + "--no-central", + "--restart", + "160", + "--max-newton", + "12", + "--capture-iterations", + "2,3,4,7,11" + ], + "notes": "Repeat prior live trajectory; capture pressure onset and shrinking-step states from observer spans without changing solver policy." + }, + { + "id": "n1p5_J0_capture", + "arguments": [ + "--index", + "1.5", + "--angular-momentum", + "0", + "--no-central", + "--restart", + "160", + "--max-newton", + "3", + "--capture-iterations", + "1,2" + ], + "notes": "Repeat first three guarded Newton attempts; capture first severely damped cold direction and subsequent guarded direction." + } +] diff --git a/experiments/plateau_diagnostic.cpp b/experiments/plateau_diagnostic.cpp index 9e3428a..fef959c 100644 --- a/experiments/plateau_diagnostic.cpp +++ b/experiments/plateau_diagnostic.cpp @@ -10,6 +10,7 @@ #include #include #include +#include #include #include #include @@ -25,6 +26,7 @@ struct Options { double n = 1., angularMomentum = 0., atol = 0., rtol = 1e-8, linearTolerance = .03; int maxNewton = 12, maxLinear = 200, restart = 40, probeRestart = 40; bool probes = false, noCentral = false; + std::vector captureIterations; std::string mesh = "sandbox.smesh"; std::filesystem::path output, probeState; }; @@ -65,6 +67,19 @@ Options Parse(int argc, char **argv) { o.probes = true; else if (a == "--no-central") o.noCentral = true; + else if (a == "--capture-iterations") { + std::istringstream entries(v()); + std::string entry; + while (std::getline(entries, entry, ',')) { + size_t used = 0; + const int iteration = std::stoi(entry, &used); + if (used != entry.size() || iteration < 0) + throw std::invalid_argument("Invalid capture iteration"); + o.captureIterations.push_back(iteration); + } + if (o.captureIterations.empty()) + throw std::invalid_argument("Empty capture iteration list"); + } else throw std::invalid_argument("Unknown option: " + a); } @@ -304,6 +319,19 @@ template int RunContext(Context &context, const Options &o, cons blocklog << "iteration,stage,block,normalized_norm,physical_norm,max_abs_normalized\n"; auto triallog = File(o.output / "trials.csv"); triallog << "iteration,trial,alpha,disposition,source,F,min_jacobian\n"; + mfem::Vector captureBase; + auto captureDirectory = [&](int iteration) { + return o.output / ("capture_" + std::to_string(iteration)); + }; + auto shouldCapture = [&](int iteration) { + return std::find(o.captureIterations.begin(), o.captureIterations.end(), iteration) != o.captureIterations.end(); + }; + auto copySpan = [](std::span values) { + mfem::Vector result(static_cast(values.size())); + for (int i = 0; i < result.Size(); ++i) + result(i) = values[i]; + return result; + }; auto record = [&](int it, const char *stage, std::span residual) { for (const auto &b : rows) { long double ps = 0; @@ -321,6 +349,14 @@ template int RunContext(Context &context, const Options &o, cons auto observer = solver::nonlinear::makeObserver( [&](const solver::nonlinear::BeforeIteration &e) { record(e.iteration, "before", e.normalizedResidual); + if (shouldCapture(e.iteration)) { + const auto directory = captureDirectory(e.iteration); + std::filesystem::create_directory(directory); + captureBase = copySpan(e.normalizedState); + SaveVector(directory / "normalized_base.txt", captureBase); + SaveVector(directory / "physical_base.txt", copySpan(e.physicalState)); + SaveVector(directory / "normalized_residual.txt", copySpan(e.normalizedResidual)); + } std::cout << "Newton " << e.iteration << " F=" << e.residualNorm << std::endl; }, [&](const solver::nonlinear::AfterLineSearchTrial &e) { @@ -328,6 +364,31 @@ template int RunContext(Context &context, const Options &o, cons << e.rejectionSource << ',' << (e.metric ? e.metric->residualNorm : NAN) << ',' << e.minimumJacobianDeterminant.value_or(NAN) << '\n'; triallog.flush(); + if (e.trial == 0 && shouldCapture(e.iteration)) { + if (!(e.stepLength > 0.) || !std::isfinite(e.stepLength) || + e.candidateNormalizedState.size() != static_cast(captureBase.Size())) + throw std::runtime_error("Invalid candidate for diagnostic capture"); + const auto directory = captureDirectory(e.iteration); + const auto candidate = copySpan(e.candidateNormalizedState); + SaveVector(directory / "normalized_candidate.txt", candidate); + SaveVector(directory / "physical_candidate.txt", copySpan(e.candidatePhysicalState)); + mfem::Vector direction(candidate), reconstructed(captureBase); + direction -= captureBase; + direction /= e.stepLength; + reconstructed.Add(e.stepLength, direction); + reconstructed -= candidate; + SaveVector(directory / "normalized_direction.txt", direction); + auto info = File(directory / "capture.txt"); + info << "iteration=" << e.iteration << "\nfirst_trial_alpha=" << e.stepLength + << "\ndisposition=" << int(e.disposition) << "\nsource=" << e.rejectionSource + << "\nbase_norm=" << captureBase.Norml2() << "\ndirection_norm=" << direction.Norml2() + << "\ncandidate_reconstruction_error=" << reconstructed.Norml2() + << "\ndirection_roundoff_scale_estimate=" + << 8. * std::numeric_limits::epsilon() * + (captureBase.Norml2() + candidate.Norml2()) / e.stepLength + << "\ndirection_source=first_trial_difference_divided_by_alpha\n" + << "roundoff_scale_is_estimate_not_certificate=1\n"; + } }, [&](const solver::nonlinear::AfterIteration &e) { record(e.iteration, "after", e.normalizedResidual); diff --git a/experiments/run_plateau_campaign.py b/experiments/run_plateau_campaign.py index 78fcf2f..0b583ba 100644 --- a/experiments/run_plateau_campaign.py +++ b/experiments/run_plateau_campaign.py @@ -124,7 +124,7 @@ def main(): command = [str(executable), *case['arguments'], '--mesh', str(mesh), '--output', str(folder)] state_files = {} for i, word in enumerate(case['arguments']): - if word == '--probe-state': + if word in ('--probe-state', '--probe-direction', '--replay-candidate'): path = Path(case['arguments'][i + 1]) if not path.is_absolute(): path = source / path @@ -144,9 +144,9 @@ def main(): record['exit_code'] = None record.update(seconds=time.monotonic() - start, finished_utc=utc()) record['output_hashes'] = { - name: digest(folder / name) - for name in ('accepted_state.txt', 'layout.csv', 'metadata.txt', 'final_physical.csv') - if (folder / name).is_file() + str(path.relative_to(folder)): digest(path) + for path in sorted(folder.rglob('*')) + if path.is_file() and path.suffix in ('.txt', '.csv', '.json') } write_json(output / (name + '_run.json'), record) summaries.append({**record, **summarize(folder)})