|
|
|
|
@@ -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 <algorithm>
|
|
|
|
|
#include <array>
|
|
|
|
|
#include <chrono>
|
|
|
|
|
#include <cmath>
|
|
|
|
|
#include <cstdint>
|
|
|
|
|
#include <cstring>
|
|
|
|
|
#include <exception>
|
|
|
|
|
#include <filesystem>
|
|
|
|
|
#include <fstream>
|
|
|
|
|
#include <iomanip>
|
|
|
|
|
#include <limits>
|
|
|
|
|
#include <stdexcept>
|
|
|
|
|
#include <string>
|
|
|
|
|
#include <string_view>
|
|
|
|
|
#include <type_traits>
|
|
|
|
|
|
|
|
|
|
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<double>(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<long double>(a(i)) * b(i);
|
|
|
|
|
return result;
|
|
|
|
|
}
|
|
|
|
|
inline double Norm(const mfem::Vector &v, int offset, int size) {
|
|
|
|
|
return std::sqrt(static_cast<double>(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<double>::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<double>::quiet_NaN()};
|
|
|
|
|
double targetDefect{std::numeric_limits<double>::quiet_NaN()};
|
|
|
|
|
double affineDisagreement{std::numeric_limits<double>::quiet_NaN()};
|
|
|
|
|
double arithmeticAllowance{std::numeric_limits<double>::quiet_NaN()};
|
|
|
|
|
double targetThreshold{std::numeric_limits<double>::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 <class State>
|
|
|
|
|
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<decltype(state.Problem())>::FormType;
|
|
|
|
|
using namespace utils::blocks;
|
|
|
|
|
constexpr value_block<type_index_v<surface_deformation::parameters::value, typename Form::value_blocks>> shape{};
|
|
|
|
|
constexpr value_block<type_index_v<gravity::gradient::value, typename Form::value_blocks>> gradient{};
|
|
|
|
|
constexpr value_block<type_index_v<gravity::poisson::value, typename Form::value_blocks>> 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<double>::quiet_NaN();
|
|
|
|
|
const double eps = std::numeric_limits<mfem::real_t>::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<const char *, 3> names{"full", "shape", "fields"};
|
|
|
|
|
constexpr std::array<double, 8> 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<mfem::Vector, 3> directions{mfem::Vector(size), mfem::Vector(size), mfem::Vector(size)};
|
|
|
|
|
std::array<mfem::Vector, 3> 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<Form>(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<std::array<d::Observation, 3>, 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<double, 8> 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<double>(d::Dot(originalResidual, actions[k], 0, size));
|
|
|
|
|
const double projection = observation.prepared
|
|
|
|
|
? static_cast<double>(d::Dot(predicted, observation.remainder, 0, size)) : nan;
|
|
|
|
|
const auto identityError = [&](int offset, int count) {
|
|
|
|
|
return observation.prepared ? static_cast<double>(
|
|
|
|
|
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<double>(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<double, 4> norms{nan, nan, nan, nan};
|
|
|
|
|
std::array<long double, 4> 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<double>(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
|
