472 lines
30 KiB
C++
472 lines
30 KiB
C++
#pragma once
|
|
|
|
// Included after import mean_field. The experiment owns observations only;
|
|
// the prepared production transaction owns correction and fallback decisions.
|
|
#include "plateau_nonlinear_step_diagnostics.hpp"
|
|
|
|
#include <optional>
|
|
#include <vector>
|
|
|
|
namespace experiment::plateau_diagnostics {
|
|
namespace coupled_trial_detail {
|
|
// Expose only the production workspace operations needed by this frozen probe.
|
|
// Keeping the affine gravity gate here would let the comparison silently drift
|
|
// away from the trial path used by the nonlinear solver.
|
|
template <class State>
|
|
struct GravityWorkspace final : mean_field::solver::detail::StellarGravityWorkspace<State> {
|
|
using Base = mean_field::solver::detail::StellarGravityWorkspace<State>;
|
|
using Base::Base;
|
|
using Base::ApplyCorrection;
|
|
using Base::BeginGravityAttempt;
|
|
using Base::Refresh;
|
|
using Base::desiredResidual;
|
|
using Base::fullCorrection;
|
|
using Base::gravityCorrection;
|
|
using Base::lastAffineDisagreement;
|
|
using Base::lastPreparedDefect;
|
|
using Base::lastPreparedTargetThreshold;
|
|
using Base::lastSolve;
|
|
using Base::subsystem;
|
|
using Base::trialSolveSucceeded;
|
|
};
|
|
}
|
|
|
|
struct CoupledTrialProbeOptions final {
|
|
double outerLinearTolerance{.03};
|
|
mean_field::solver::nonlinear::GravitySolveOptions gravity{};
|
|
mean_field::solver::CoupledTrialCorrectionOptions coupled{};
|
|
bool gravityPreconditionPath{true};
|
|
bool seedControl{false};
|
|
double sufficientDecrease{1.e-4};
|
|
};
|
|
|
|
template <class State>
|
|
void ProbeCoupledTrial(State &state, const mfem::Vector &normalizedBase,
|
|
const mfem::Vector &suppliedDirection, const std::filesystem::path &output,
|
|
CoupledTrialProbeOptions options = {}) {
|
|
using namespace mean_field;
|
|
namespace d = nonlinear_step_detail;
|
|
using GravityWorkspace = coupled_trial_detail::GravityWorkspace<State>;
|
|
using GravityCorrection = solver::PreparedGravityCurvatureCorrection<typename GravityWorkspace::Subsystem>;
|
|
using Transaction = solver::detail::PreparedCoupledTrialTransaction<State>;
|
|
const auto communicator = state.Problem().GetCommunicator();
|
|
int ranks = 0;
|
|
MPI_Comm_size(communicator, &ranks);
|
|
const int size = state.Problem().StateSize();
|
|
if (ranks != 1 || size != state.Problem().EquationSize() || normalizedBase.Size() != size ||
|
|
suppliedDirection.Size() != size || !d::Finite(normalizedBase) || !d::Finite(suppliedDirection) ||
|
|
!std::isfinite(options.outerLinearTolerance) || options.outerLinearTolerance <= 0. ||
|
|
options.outerLinearTolerance >= 1. || !std::isfinite(options.sufficientDecrease) ||
|
|
options.sufficientDecrease <= 0. || options.sufficientDecrease >= 1.)
|
|
throw std::invalid_argument("Coupled-trial probes require one rank, compatible finite vectors and valid controls.");
|
|
options.gravity.Validate();
|
|
options.coupled.Validate();
|
|
const auto start = d::Clock::now();
|
|
const double nan = std::numeric_limits<double>::quiet_NaN();
|
|
const mfem::Vector base(normalizedBase), direction(suppliedDirection);
|
|
const auto &normalization = state.normalizedOperator->GetNormalization();
|
|
const mfem::Vector factors(normalization.StateFactors()), residualFactors(normalization.ResidualFactors());
|
|
const mfem::Vector savedAccepted(state.acceptedNormalizedState),
|
|
savedAcceptedResidual(state.acceptedNormalizedResidual), savedPhysical(state.AcceptedPhysicalState()),
|
|
savedCorrection(state.normalizedCorrection), savedPhysicalCorrection(state.physicalCorrection),
|
|
savedVolumeDirection(state.volumeDisplacementDirection);
|
|
const double savedMinimumJacobian = state.acceptedMinimumJacobianDeterminant;
|
|
std::filesystem::create_directories(output);
|
|
const auto file = [&](const std::string &name) { return output / ("coupled_trial_" + name); };
|
|
auto summary = d::Open(file("summary.csv"));
|
|
auto preparations = d::Open(file("preparations.csv"));
|
|
auto trials = d::Open(file("trials.csv"));
|
|
auto blocks = d::Open(file("blocks.csv"));
|
|
auto manifest = d::Open(file("layout.csv"));
|
|
summary << "metric,value\n";
|
|
preparations << "sequence,kind,sample,alpha,prepared,rejection_reason,rejection_stage,seconds,residual_norm\n";
|
|
trials << "sample,alpha,gravity_path,raw_prepared,raw_rejection_reason,raw_rejection_stage,raw_norm,raw_merit,"
|
|
"gravity_attempted,gravity_reported_gate,gravity_prepared_gate,gravity_target_defect,"
|
|
"gravity_affine_disagreement,gravity_target_threshold,before_prepared,before_usable,before_norm,"
|
|
"before_merit,before_accepted,soc_called,soc_status,soc_attempted,aux_status,aux_reported_gate,"
|
|
"aux_iterations,aux_rhs_norm,aux_reported_true_norm,aux_verified_true_norm,aux_verified_gate,"
|
|
"aux_threshold,aux_verification_applied,correction_norm,correction_cap,correction_ratio,"
|
|
"corrected_metric_observed,corrected_norm,corrected_merit,corrected_accepted,selected_prepared,"
|
|
"selected_norm,selected_merit,selected_accepted,selected_is_original,fallback_residual_difference,"
|
|
"raw_seconds,gravity_seconds,soc_seconds,soc_refresh_seconds,soc_solve_seconds,"
|
|
"soc_verification_seconds,soc_preparation_seconds,soc_preparations,soc_restorations,"
|
|
"aux_operator_applications,aux_inverse_applications,verified_action_defect_norm,"
|
|
"capture_seconds,total_sample_seconds,aux_core_status,rhs_saved,correction_saved,verified_vectors_saved\n";
|
|
blocks << "sample,alpha,stage,block,offset,size,residual_norm\n";
|
|
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 : state.Problem().GetManifest().residualBlocks())
|
|
manifest << "residual," << row.stableId << ',' << row.offset << ',' << row.size << '\n';
|
|
manifest.flush();
|
|
const auto metric = [&](const char *name, auto value) { summary << name << ',' << value << '\n'; };
|
|
metric("schema_version", 1);
|
|
metric("independent_frozen_alphas", 1);
|
|
metric("accepted_commits", 0);
|
|
metric("outer_linear_tolerance", options.outerLinearTolerance);
|
|
metric("sufficient_decrease", options.sufficientDecrease);
|
|
metric("gravity_path", options.gravityPreconditionPath);
|
|
metric("seed_control", options.seedControl);
|
|
metric("coupled_maximum_correction_ratio", options.coupled.maximumCorrectionRatio);
|
|
metric("coupled_relative_tolerance", options.coupled.control.relativeTolerance);
|
|
metric("coupled_absolute_tolerance", options.coupled.control.absoluteTolerance);
|
|
metric("coupled_maximum_iterations", options.coupled.control.maximumIterations);
|
|
metric("coupled_restart", options.coupled.backend.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", options.gravity.backend.restartLength);
|
|
metric("polynomial_direction_repairs", 0);
|
|
metric("sampled_enthalpy_policy", 1);
|
|
metric("coupled_attempt_budget_per_independent_alpha", 1);
|
|
metric("soc_status_not_attempted", int(solver::CoupledTrialAttemptStatus::not_attempted));
|
|
metric("soc_status_original_accepted", int(solver::CoupledTrialAttemptStatus::original_accepted));
|
|
metric("soc_status_budget_exhausted", int(solver::CoupledTrialAttemptStatus::budget_exhausted));
|
|
metric("soc_status_kept", int(solver::CoupledTrialAttemptStatus::kept));
|
|
metric("soc_status_linear_failure", int(solver::CoupledTrialAttemptStatus::linear_failure));
|
|
metric("soc_status_norm_limit", int(solver::CoupledTrialAttemptStatus::norm_limit));
|
|
metric("soc_status_non_finite", int(solver::CoupledTrialAttemptStatus::non_finite));
|
|
metric("soc_status_inadmissible", int(solver::CoupledTrialAttemptStatus::inadmissible));
|
|
metric("soc_status_not_improved", int(solver::CoupledTrialAttemptStatus::not_improved));
|
|
std::vector<double> alphas = options.seedControl ? std::vector<double>{1.}
|
|
: std::vector<double>{1. / 16., 1. / 64., 1. / 256., 1. / 1024.};
|
|
for (std::size_t sample = 0; sample < alphas.size(); ++sample)
|
|
summary << "alpha_" << sample << ',' << alphas[sample] << '\n';
|
|
mfem::Vector residual(size), baseAction(size);
|
|
bool captured = false, haveBaseAction = false, finished = false;
|
|
int preparationsCount = 0, refreshes = 0, directActions = 0;
|
|
double preparationSeconds = 0., refreshSeconds = 0., actionSeconds = 0., restorationError = 0.;
|
|
const auto prepare = [&](const mfem::Vector &candidate, const char *kind, int sample, double alpha) {
|
|
state.trialNormalizedState = candidate;
|
|
state.normalizedOperator->DenormalizeState(candidate, state.candidatePhysicalState);
|
|
const auto started = d::Clock::now();
|
|
auto result = state.PrepareTrial();
|
|
const double seconds = d::Seconds(started);
|
|
preparationSeconds += seconds;
|
|
preparations << ++preparationsCount << ',' << kind << ',' << sample << ',' << alpha << ',' << bool(result)
|
|
<< ',' << (result ? -1 : int(result.error().reason)) << ','
|
|
<< (result ? -1 : int(result.error().stage)) << ',' << seconds << ','
|
|
<< (result ? d::Norm(state.trialNormalizedResidual) : nan) << '\n';
|
|
preparations.flush();
|
|
if (result && !d::Finite(state.trialNormalizedResidual))
|
|
throw std::runtime_error("Coupled-trial preparation produced a nonfinite residual.");
|
|
return result;
|
|
};
|
|
const auto refresh = [&]() {
|
|
const auto started = d::Clock::now();
|
|
state.physicalInverse->Refresh();
|
|
refreshSeconds += d::Seconds(started);
|
|
++refreshes;
|
|
};
|
|
const auto restore = [&](const char *kind) {
|
|
if (!prepare(base, kind, -1, 0.)) throw std::runtime_error("The coupled-trial base cannot be restored.");
|
|
if (!captured) { residual = state.trialNormalizedResidual; captured = true; }
|
|
restorationError = std::max(restorationError, d::Difference(residual, state.trialNormalizedResidual));
|
|
refresh();
|
|
};
|
|
const auto action = [&](const mfem::Vector &input, mfem::Vector &result) {
|
|
const auto started = d::Clock::now();
|
|
state.normalizedOperator->Mult(input, result);
|
|
actionSeconds += d::Seconds(started);
|
|
++directActions;
|
|
if (!d::Finite(result)) throw std::runtime_error("The coupled-trial full-J action is nonfinite.");
|
|
};
|
|
const auto evaluate = [&](const mfem::Vector &value) {
|
|
return solver::nonlinear::getMetric(solver::nonlinear::NormalizedResidualMetric{}, value, communicator);
|
|
};
|
|
const auto blockRows = [&](std::size_t sample, double alpha, const char *stage, const mfem::Vector &value) {
|
|
blocks << sample << ',' << alpha << ',' << stage << ",all,0," << size << ',' << d::Norm(value) << '\n';
|
|
for (const auto &row : state.Problem().GetManifest().residualBlocks())
|
|
blocks << sample << ',' << alpha << ',' << stage << ',' << row.stableId << ',' << row.offset << ','
|
|
<< row.size << ',' << d::Norm(value, row.offset, row.size) << '\n';
|
|
blocks.flush();
|
|
};
|
|
const auto finish = [&]() {
|
|
// Output is append-only evidence. A failed cleanup must not retry a
|
|
// previously written filename or replace the original failure record.
|
|
finished = true;
|
|
const bool directionsNeverChanged = d::Identical(savedCorrection, state.normalizedCorrection) &&
|
|
d::Identical(savedPhysicalCorrection, state.physicalCorrection) &&
|
|
d::Identical(savedVolumeDirection, state.volumeDisplacementDirection);
|
|
state.normalizedCorrection = savedCorrection;
|
|
state.physicalCorrection = savedPhysicalCorrection;
|
|
state.volumeDisplacementDirection = savedVolumeDirection;
|
|
restore("final_base_restore");
|
|
d::Save(file("restored_residual.txt"), state.trialNormalizedResidual);
|
|
mfem::Vector restoredAction(size);
|
|
action(direction, restoredAction);
|
|
d::Save(file("restored_action.txt"), restoredAction);
|
|
const bool acceptedIdentical = d::Identical(savedAccepted, state.acceptedNormalizedState) &&
|
|
d::Identical(savedAcceptedResidual, state.acceptedNormalizedResidual) &&
|
|
d::Identical(savedPhysical, state.AcceptedPhysicalState()) &&
|
|
savedMinimumJacobian == state.acceptedMinimumJacobianDeterminant;
|
|
const bool factorsIdentical = d::Identical(factors, normalization.StateFactors()) &&
|
|
d::Identical(residualFactors, normalization.ResidualFactors());
|
|
const bool directionsIdentical = d::Identical(savedCorrection, state.normalizedCorrection) &&
|
|
d::Identical(savedPhysicalCorrection, state.physicalCorrection) &&
|
|
d::Identical(savedVolumeDirection, state.volumeDisplacementDirection);
|
|
metric("accepted_buffers_bitwise_identical", acceptedIdentical);
|
|
metric("normalization_factors_bitwise_identical", factorsIdentical);
|
|
metric("direction_buffers_bitwise_identical", directionsIdentical);
|
|
metric("direction_buffers_unchanged_before_restoration", directionsNeverChanged);
|
|
metric("maximum_base_restoration_error", restorationError);
|
|
metric("restored_base_action_error", haveBaseAction ? d::Difference(baseAction, restoredAction) : nan);
|
|
metric("fresh_operator_and_inverse_current", state.normalizedOperator->IsPrepared() &&
|
|
state.physicalInverse->IsCurrent() && state.normalizedInverse->IsCurrent());
|
|
metric("explicit_probe_preparations", preparationsCount);
|
|
metric("explicit_probe_preparation_seconds", preparationSeconds);
|
|
metric("explicit_probe_inverse_refreshes", refreshes);
|
|
metric("explicit_probe_inverse_refresh_seconds", refreshSeconds);
|
|
metric("explicit_probe_operator_applications", directActions);
|
|
metric("explicit_probe_operator_seconds", actionSeconds);
|
|
metric("total_seconds", d::Seconds(start));
|
|
summary.flush();
|
|
if (!acceptedIdentical || !factorsIdentical || !directionsIdentical || !directionsNeverChanged)
|
|
throw std::runtime_error("The coupled-trial probe changed accepted buffers, directions or normalization.");
|
|
};
|
|
|
|
try {
|
|
restore("initial_base_restore");
|
|
mfem::Vector physicalBase(size);
|
|
state.normalizedOperator->DenormalizeState(base, physicalBase);
|
|
action(direction, baseAction);
|
|
haveBaseAction = true;
|
|
mfem::Vector defect(baseAction);
|
|
defect += residual;
|
|
const auto baseMetric = evaluate(residual);
|
|
const double outerThreshold = options.outerLinearTolerance * baseMetric.residualNorm;
|
|
const bool outerGate = d::Norm(defect) <= outerThreshold;
|
|
d::Save(file("base_normalized.txt"), base);
|
|
d::Save(file("base_physical.txt"), physicalBase);
|
|
d::Save(file("base_residual.txt"), residual);
|
|
d::Save(file("state_factors.txt"), factors);
|
|
d::Save(file("residual_factors.txt"), residualFactors);
|
|
d::Save(file("direction.txt"), direction);
|
|
d::Save(file("base_action.txt"), baseAction);
|
|
d::Save(file("base_defect.txt"), defect);
|
|
metric("base_residual_norm", baseMetric.residualNorm);
|
|
metric("base_merit", baseMetric.merit);
|
|
metric("direction_norm", d::Norm(direction));
|
|
metric("base_merit_slope", d::Dot(residual, baseAction, 0, size));
|
|
metric("outer_verified_defect_norm", d::Norm(defect));
|
|
metric("outer_verified_relative_defect", d::Ratio(d::Norm(defect), baseMetric.residualNorm));
|
|
metric("outer_verified_threshold", outerThreshold);
|
|
metric("outer_verified_gate", outerGate);
|
|
metric("outer_gate_failure_skips_trials", !outerGate);
|
|
summary.flush();
|
|
if (!outerGate) { metric("probe_exception", 0); finish(); return; }
|
|
|
|
Transaction transaction(state, options.coupled, 1);
|
|
std::optional<GravityWorkspace> gravity;
|
|
std::optional<GravityCorrection> gravityCorrection;
|
|
if (options.gravityPreconditionPath) {
|
|
gravity.emplace(state, options.gravity);
|
|
gravityCorrection.emplace(gravity->subsystem);
|
|
metric("compact_gravity_operator", gravity->subsystem.UsesCompactOperator());
|
|
metric("compact_gravity_inverse", gravity->subsystem.UsesCompactPreconditioner());
|
|
}
|
|
for (std::size_t sample = 0; sample < alphas.size(); ++sample) {
|
|
const auto sampleStart = d::Clock::now();
|
|
const double alpha = alphas[sample];
|
|
const std::string suffix = "_a" + std::to_string(sample) + ".txt";
|
|
restore("capture_base");
|
|
const auto captureStart = d::Clock::now();
|
|
// These are independent frozen trials. Resetting here does not
|
|
// represent several attempts inside one production line search.
|
|
transaction.CaptureBase(residual, direction);
|
|
transaction.BeginTrial();
|
|
if (gravity) {
|
|
gravity->subsystem.Refresh();
|
|
gravityCorrection->CaptureBase(residual, direction);
|
|
gravity->BeginGravityAttempt();
|
|
}
|
|
const double captureSeconds = d::Seconds(captureStart);
|
|
mfem::Vector prediction(residual), rawState(base);
|
|
prediction.Add(alpha, baseAction);
|
|
rawState.Add(alpha, direction);
|
|
d::Save(file("prediction" + suffix), prediction);
|
|
d::Save(file("raw_state" + suffix), rawState);
|
|
const auto rawStart = d::Clock::now();
|
|
auto preparation = prepare(rawState, "raw", static_cast<int>(sample), alpha);
|
|
const double rawSeconds = d::Seconds(rawStart);
|
|
const bool rawPrepared = bool(preparation);
|
|
const int rawReason = rawPrepared ? -1 : int(preparation.error().reason);
|
|
const int rawStage = rawPrepared ? -1 : int(preparation.error().stage);
|
|
solver::nonlinear::MetricEvaluation rawMetric{nan, nan};
|
|
if (rawPrepared) {
|
|
rawMetric = evaluate(state.trialNormalizedResidual);
|
|
d::Save(file("raw_residual" + suffix), state.trialNormalizedResidual);
|
|
blockRows(sample, alpha, "raw", state.trialNormalizedResidual);
|
|
}
|
|
bool gravityAttempted = false, gravityReportedGate = !gravity, gravityPreparedGate = !gravity;
|
|
double gravitySeconds = 0.;
|
|
if (rawPrepared && gravity) {
|
|
const auto gravityStart = d::Clock::now();
|
|
gravity->Refresh(state);
|
|
gravityCorrection->BuildLinearPrediction(alpha, gravity->desiredResidual);
|
|
gravityAttempted = true;
|
|
gravity->lastSolve = gravityCorrection->SolveCorrection(
|
|
alpha, state.trialNormalizedResidual, gravity->gravityCorrection, options.gravity.control);
|
|
gravity->trialSolveSucceeded = gravity->lastSolve.Converged();
|
|
gravityReportedGate = gravity->trialSolveSucceeded;
|
|
if (gravityReportedGate) {
|
|
// This is the same production affine-gravity assessment,
|
|
// including its contribution-aware arithmetic allowance.
|
|
preparation = gravity->ApplyCorrection(state);
|
|
gravityPreparedGate = bool(preparation) && gravity->TrialSolveSucceeded();
|
|
d::Save(file("gravity_correction" + suffix), gravity->fullCorrection);
|
|
}
|
|
gravitySeconds = d::Seconds(gravityStart);
|
|
}
|
|
const bool beforePrepared = bool(preparation);
|
|
const bool usable = beforePrepared && gravityReportedGate && gravityPreparedGate;
|
|
mfem::Vector beforeState(state.trialNormalizedState), beforeResidual;
|
|
solver::nonlinear::MetricEvaluation beforeMetric{nan, nan};
|
|
if (beforePrepared) {
|
|
beforeResidual = state.trialNormalizedResidual;
|
|
beforeMetric = evaluate(beforeResidual);
|
|
d::Save(file("before_state" + suffix), beforeState);
|
|
d::Save(file("before_residual" + suffix), beforeResidual);
|
|
blockRows(sample, alpha, "before", beforeResidual);
|
|
}
|
|
const auto accepts = [&](const solver::nonlinear::MetricEvaluation &candidate) {
|
|
return std::isfinite(candidate.merit) && candidate.merit >= 0. &&
|
|
candidate.merit <= (1. - options.sufficientDecrease * alpha) * baseMetric.merit;
|
|
};
|
|
const auto better = [](const auto &candidate, const auto &original) {
|
|
return std::isfinite(candidate.merit) && candidate.merit >= 0. &&
|
|
std::isfinite(candidate.residualNorm) && candidate.residualNorm >= 0. &&
|
|
candidate.merit < original.merit;
|
|
};
|
|
const bool beforeAccepted = usable && accepts(beforeMetric);
|
|
auto selectedMetric = beforeMetric;
|
|
bool correctedObserved = false;
|
|
mfem::Vector observedCorrectedState, observedCorrectedResidual;
|
|
solver::nonlinear::MetricEvaluation observedCorrectedMetric{nan, nan};
|
|
const auto observeMetric = [&](const mfem::Vector &value) {
|
|
const auto result = evaluate(value);
|
|
// This callback observes the actual prepared correction before
|
|
// a possible merit-based rollback. It neither reconstructs nor
|
|
// reapplies the remedy and does not change its selection rule.
|
|
if (!d::Identical(state.trialNormalizedState, beforeState)) {
|
|
observedCorrectedState = state.trialNormalizedState;
|
|
observedCorrectedResidual = value;
|
|
observedCorrectedMetric = result;
|
|
correctedObserved = true;
|
|
}
|
|
return result;
|
|
};
|
|
const auto previousDiagnostics = transaction.GetDiagnostics();
|
|
double socSeconds = 0.;
|
|
if (usable) {
|
|
const auto socStart = d::Clock::now();
|
|
transaction.TryImprove(state, alpha, preparation, selectedMetric, observeMetric,
|
|
accepts, better, [](const mfem::Vector &) { return true; });
|
|
socSeconds = d::Seconds(socStart);
|
|
}
|
|
const auto &diagnostics = transaction.GetDiagnostics();
|
|
const auto &report = diagnostics.currentSolve;
|
|
const bool attempted = report.has_value();
|
|
const bool rhsSaved = attempted && d::Finite(transaction.GetCore().GetRightHandSide());
|
|
const bool correctionSaved = attempted && d::Finite(transaction.GetCorrection());
|
|
const bool verifiedVectorsSaved = attempted && report->verificationApplied && rhsSaved &&
|
|
d::Finite(transaction.GetCore().GetVerifiedAction()) && d::Finite(transaction.GetCore().GetVerifiedDefect());
|
|
double savedDefectNorm = nan;
|
|
if (attempted) {
|
|
if (rhsSaved) d::Save(file("soc_rhs" + suffix), transaction.GetCore().GetRightHandSide());
|
|
if (correctionSaved) d::Save(file("soc_correction" + suffix), transaction.GetCorrection());
|
|
if (verifiedVectorsSaved) {
|
|
d::Save(file("soc_action" + suffix), transaction.GetCore().GetVerifiedAction());
|
|
mfem::Vector savedDefect(transaction.GetCore().GetVerifiedAction());
|
|
savedDefect -= transaction.GetCore().GetRightHandSide();
|
|
savedDefectNorm = d::Norm(savedDefect);
|
|
d::Save(file("soc_defect" + suffix), savedDefect);
|
|
}
|
|
}
|
|
if (correctedObserved) {
|
|
d::Save(file("corrected_state" + suffix), observedCorrectedState);
|
|
d::Save(file("corrected_residual" + suffix), observedCorrectedResidual);
|
|
mfem::Vector physical(size);
|
|
state.normalizedOperator->DenormalizeState(observedCorrectedState, physical);
|
|
d::Save(file("corrected_physical" + suffix), physical);
|
|
blockRows(sample, alpha, "corrected", observedCorrectedResidual);
|
|
}
|
|
const bool selectedPrepared = bool(preparation);
|
|
const bool selectedIsOriginal = beforePrepared && d::Identical(beforeState, state.trialNormalizedState);
|
|
double fallbackDifference = nan;
|
|
if (selectedPrepared) {
|
|
selectedMetric = evaluate(state.trialNormalizedResidual);
|
|
d::Save(file("selected_state" + suffix), state.trialNormalizedState);
|
|
d::Save(file("selected_residual" + suffix), state.trialNormalizedResidual);
|
|
d::Save(file("selected_physical" + suffix), state.candidatePhysicalState);
|
|
blockRows(sample, alpha, "selected", state.trialNormalizedResidual);
|
|
if (selectedIsOriginal) fallbackDifference = d::Difference(beforeResidual, state.trialNormalizedResidual);
|
|
}
|
|
const bool selectedAccepted = usable && selectedPrepared && accepts(selectedMetric);
|
|
trials << sample << ',' << alpha << ',' << options.gravityPreconditionPath << ',' << rawPrepared << ','
|
|
<< rawReason << ',' << rawStage << ',' << rawMetric.residualNorm << ',' << rawMetric.merit << ','
|
|
<< gravityAttempted << ',' << gravityReportedGate << ',' << gravityPreparedGate << ','
|
|
<< (gravityAttempted ? gravity->lastPreparedDefect : nan) << ','
|
|
<< (gravityAttempted ? gravity->lastAffineDisagreement : nan) << ','
|
|
<< (gravityAttempted ? gravity->lastPreparedTargetThreshold : nan) << ','
|
|
<< beforePrepared << ',' << usable << ',' << beforeMetric.residualNorm << ',' << beforeMetric.merit
|
|
<< ',' << beforeAccepted << ',' << usable << ',' << int(diagnostics.status) << ',' << attempted
|
|
<< ',' << (attempted ? int(report->linearSolve.status) : -1) << ','
|
|
<< (attempted && report->linearSolve.Converged()) << ','
|
|
<< (attempted ? report->linearSolve.iterations : 0) << ','
|
|
<< (attempted ? report->linearSolve.rightHandSideNorm : nan) << ','
|
|
<< (attempted ? report->linearSolve.trueResidualNorm : nan) << ','
|
|
<< (attempted ? report->verifiedTrueResidualNorm : nan) << ','
|
|
<< (attempted && report->verifiedGate) << ','
|
|
<< (attempted ? report->convergenceThreshold : nan) << ','
|
|
<< (attempted && report->verificationApplied) << ','
|
|
<< (attempted ? report->correctionNorm : nan) << ','
|
|
<< (attempted ? report->maximumCorrectionNorm : nan) << ','
|
|
<< (attempted ? d::Ratio(report->correctionNorm, alpha * d::Norm(direction)) : nan) << ','
|
|
<< correctedObserved << ',' << observedCorrectedMetric.residualNorm << ','
|
|
<< observedCorrectedMetric.merit << ',' << (correctedObserved && accepts(observedCorrectedMetric)) << ','
|
|
<< selectedPrepared << ',' << selectedMetric.residualNorm << ',' << selectedMetric.merit << ','
|
|
<< selectedAccepted << ',' << selectedIsOriginal << ',' << fallbackDifference << ','
|
|
<< rawSeconds << ',' << gravitySeconds << ',' << socSeconds << ','
|
|
<< diagnostics.refreshSeconds - previousDiagnostics.refreshSeconds << ','
|
|
<< diagnostics.solveSeconds - previousDiagnostics.solveSeconds << ','
|
|
<< diagnostics.verificationSeconds - previousDiagnostics.verificationSeconds << ','
|
|
<< diagnostics.preparationSeconds - previousDiagnostics.preparationSeconds << ','
|
|
<< diagnostics.preparations - previousDiagnostics.preparations << ','
|
|
<< diagnostics.restorations - previousDiagnostics.restorations << ','
|
|
<< (attempted ? report->linearSolve.operatorApplications : 0) << ','
|
|
<< (attempted ? report->linearSolve.inversePreconditionerApplications : 0) << ','
|
|
<< savedDefectNorm << ',' << captureSeconds << ',' << d::Seconds(sampleStart) << ','
|
|
<< (attempted ? int(report->status) : -1) << ',' << rhsSaved << ',' << correctionSaved << ','
|
|
<< verifiedVectorsSaved << '\n';
|
|
trials.flush();
|
|
summary.flush();
|
|
}
|
|
const auto &totals = transaction.GetDiagnostics();
|
|
metric("soc_attempts", totals.attempts);
|
|
metric("soc_kept", totals.kept);
|
|
metric("soc_fallbacks", totals.fallbacks);
|
|
metric("soc_iterations", totals.iterations);
|
|
metric("soc_solve_seconds", totals.solveSeconds);
|
|
metric("soc_verification_seconds", totals.verificationSeconds);
|
|
metric("soc_preparation_seconds", totals.preparationSeconds);
|
|
metric("soc_refresh_seconds", totals.refreshSeconds);
|
|
metric("soc_preparations", totals.preparations);
|
|
metric("soc_restorations", totals.restorations);
|
|
if (gravity) {
|
|
metric("gravity_solves", gravity->subsystem.GetSolveCount());
|
|
metric("gravity_iterations", gravity->subsystem.GetIterationCount());
|
|
metric("gravity_failures", gravity->subsystem.GetFailureCount());
|
|
metric("gravity_solve_seconds", gravity->subsystem.GetSolveSeconds());
|
|
}
|
|
} catch (...) {
|
|
const auto failure = std::current_exception();
|
|
metric("probe_exception", 1);
|
|
if (!finished) finish();
|
|
std::rethrow_exception(failure);
|
|
}
|
|
metric("probe_exception", 0);
|
|
finish();
|
|
}
|
|
} // namespace experiment::plateau_diagnostics
|