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

Distinguish strategy failures and expose the actual trial tangent

Browse Source
This commit is contained in:
Emily Boudreaux
2026-09-17 10:09:22 -04:00
parent 5f3eb0b18c
commit 2e31917269
6 changed files with 152 additions and 13 deletions
Download Patch File Download Diff File Expand all files Collapse all files

30
experiments/plateau_diagnostic.cpp
View File

@@ -374,7 +374,7 @@ void FrozenProbes(State &s, const fem::FEM &f, const Options &o, const mfem::Vec
template <class Context> int RunContext(Context &context, const Options &o, const Scales &scales) {
std::vector<operators::RootBlockDescriptor> rows, values;
mfem::Vector residualFactors;
mfem::Vector residualFactors, stateFactors;
std::vector<int> centralRows;
int centralOffset = -1;
auto metadata = File(o.output / "metadata.txt");
@@ -441,6 +441,7 @@ template <class Context> int RunContext(Context &context, const Options &o, cons
rows.assign(r.begin(), r.end());
values.assign(v.begin(), v.end());
residualFactors = s.normalizedOperator->GetNormalization().ResidualFactors();
stateFactors = s.normalizedOperator->GetNormalization().StateFactors();
for (const auto &b : values)
if (b.stableId == "fixed_central_density.border")
centralOffset = b.offset;
@@ -462,6 +463,8 @@ 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";
auto stateConsistency = File(o.output / "trial_state_consistency.csv");
stateConsistency << "iteration,trial,normalized_physical_roundtrip_error,relative_error\n";
mfem::Vector captureBase;
auto captureDirectory = [&](int iteration) {
return o.output / ("capture_" + std::to_string(iteration));
@@ -503,6 +506,20 @@ template <class Context> int RunContext(Context &context, const Options &o, cons
std::cout << "Newton " << e.iteration << " F=" << e.residualNorm << std::endl;
},
[&](const solver::nonlinear::AfterLineSearchTrial &e) {
if (e.candidatePhysicalState.size() != e.candidateNormalizedState.size() ||
e.candidateNormalizedState.size() != static_cast<size_t>(stateFactors.Size()))
throw std::runtime_error("Trial observer exposed incompatible state representations");
long double stateError = 0.0L, stateNorm = 0.0L;
for (int i = 0; i < stateFactors.Size(); ++i) {
const long double difference = e.candidatePhysicalState[i] * stateFactors(i) - e.candidateNormalizedState[i];
stateError += difference*difference;
stateNorm += static_cast<long double>(e.candidateNormalizedState[i])*e.candidateNormalizedState[i];
}
const double relativeStateError = std::sqrt(stateError)/std::max(1.0L,std::sqrt(stateNorm));
stateConsistency << e.iteration << ',' << e.trial << ',' << std::sqrt(stateError) << ',' << relativeStateError << '\n';
stateConsistency.flush();
if (std::isfinite(relativeStateError) && relativeStateError > 64.0*std::numeric_limits<double>::epsilon())
throw std::runtime_error("Trial observer physical state disagrees with its normalized state");
triallog << e.iteration << ',' << e.trial << ',' << e.stepLength << ',' << int(e.disposition) << ','
<< e.rejectionSource << ',' << (e.metric ? e.metric->residualNorm : NAN) << ','
<< e.minimumJacobianDeterminant.value_or(NAN) << '\n';
@@ -520,7 +537,8 @@ template <class Context> int RunContext(Context &context, const Options &o, cons
direction /= e.stepLength;
reconstructed.Add(e.stepLength, direction);
reconstructed -= candidate;
SaveVector(directory / "normalized_direction.txt", direction);
SaveVector(directory / "normalized_effective_trial_direction.txt", direction);
SaveVector(directory / "normalized_direction.txt", copySpan(e.normalizedDirection));
auto info = File(directory / "capture.txt");
info << "iteration=" << e.iteration << "\nfirst_trial_alpha=" << e.stepLength
<< "\ndisposition=" << int(e.disposition) << "\nsource=" << e.rejectionSource
@@ -529,7 +547,9 @@ template <class Context> int RunContext(Context &context, const Options &o, cons
<< "\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"
<< "\ndirection_source=exact_observer_tangent\n"
<< "effective_trial_direction_source=candidate_difference_divided_by_alpha\n"
<< "roundoff_scale_applies_to_effective_trial_direction_only=1\n"
<< "roundoff_scale_is_estimate_not_certificate=1\n";
}
},
@@ -573,7 +593,7 @@ template <class Context> int RunContext(Context &context, const Options &o, cons
const nl::GravitySolveOptions gravity{
.control = {.relativeTolerance = o.gravityTolerance, .absoluteTolerance = 1e-14,
.maximumIterations = o.gravityIterations},
.backend = {.restartLength = o.gravityRestart}};
.backend = {.restartLength = o.gravityRestart, .printLevel = -1}};
auto run = [&](auto numerics) {
auto solve = solver::make(context, nl::Newton(options, nl::NormalizedResidualMetric{}, numerics), observer);
return solve.evaluate();
@@ -597,6 +617,8 @@ template <class Context> int RunContext(Context &context, const Options &o, cons
metric("gravity_solves", strategyDiagnostics.gravitySolves);
metric("gravity_iterations", strategyDiagnostics.gravityIterations);
metric("gravity_failures", strategyDiagnostics.gravityFailures);
metric("gravity_prepared_failures", strategyDiagnostics.gravityPreparedFailures);
metric("linear_report_space", static_cast<int>(strategyDiagnostics.lastLinearReportSpace));
metric("gravity_solve_seconds", strategyDiagnostics.gravitySolveSeconds);
metric("last_gravity_true_residual", strategyDiagnostics.lastGravityTrueResidual);
metric("last_gravity_correction_norm", strategyDiagnostics.lastGravityCorrectionNorm);