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

experiment: capture frozen live Newton trials without changing solver policy

Browse Source
This commit is contained in:
Emily Boudreaux
2026-09-17 08:06:37 -04:00
parent 95d9fbafe8
commit fd08f7b0e9
3 changed files with 101 additions and 4 deletions
Download Patch File Download Diff File Expand all files Collapse all files

61
experiments/plateau_diagnostic.cpp
View File

@@ -10,6 +10,7 @@
#include <mfem.hpp>
#include <numbers>
#include <span>
#include <sstream>
#include <stdexcept>
#include <string>
#include <vector>
@@ -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<int> 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 <class Context> 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<const mfem::real_t> values) {
mfem::Vector result(static_cast<int>(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<const mfem::real_t> residual) {
for (const auto &b : rows) {
long double ps = 0;
@@ -321,6 +349,14 @@ template <class Context> 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 <class Context> 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<size_t>(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<double>::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);