658 lines
42 KiB
C++
658 lines
42 KiB
C++
#pragma once
|
|
|
|
#include <algorithm>
|
|
#include <chrono>
|
|
#include <cmath>
|
|
#include <concepts>
|
|
#include <cstring>
|
|
#include <cstdint>
|
|
#include <filesystem>
|
|
#include <fstream>
|
|
#include <iomanip>
|
|
#include <iostream>
|
|
#include <limits>
|
|
#include <memory>
|
|
|
|
#include "plateau_coupled_coarse.hpp"
|
|
#include "plateau_surface_border_schur.hpp"
|
|
#include <stdexcept>
|
|
#include <string>
|
|
#include <type_traits>
|
|
#include <utility>
|
|
|
|
// Include after import mean_field. Numerical recipes remain typed; the runtime
|
|
// selection below belongs only to this single-rank experimental driver.
|
|
namespace experiment::plateau_diagnostics {
|
|
enum class PreconditionerBenchmarkVariant {
|
|
independent,
|
|
material_then_gravity,
|
|
gravity_then_material,
|
|
stellar_ldu,
|
|
material_ldu,
|
|
gravity_schur,
|
|
stellar_ldu_gravity_schur,
|
|
coupled_shape_coarse,
|
|
gravity_schur_coarse,
|
|
surface_border_schur,
|
|
gravity_schur_surface_border_schur
|
|
};
|
|
|
|
inline PreconditionerBenchmarkVariant ParsePreconditionerBenchmarkVariant(const std::string &name) {
|
|
if (name == "independent") return PreconditionerBenchmarkVariant::independent;
|
|
if (name == "material-then-gravity") return PreconditionerBenchmarkVariant::material_then_gravity;
|
|
if (name == "gravity-then-material") return PreconditionerBenchmarkVariant::gravity_then_material;
|
|
if (name == "stellar-ldu") return PreconditionerBenchmarkVariant::stellar_ldu;
|
|
if (name == "material-ldu") return PreconditionerBenchmarkVariant::material_ldu;
|
|
if (name == "gravity-schur") return PreconditionerBenchmarkVariant::gravity_schur;
|
|
if (name == "stellar-ldu-gravity-schur") return PreconditionerBenchmarkVariant::stellar_ldu_gravity_schur;
|
|
if (name == "coupled-shape-coarse") return PreconditionerBenchmarkVariant::coupled_shape_coarse;
|
|
if (name == "gravity-schur-coarse") return PreconditionerBenchmarkVariant::gravity_schur_coarse;
|
|
if (name == "surface-border-schur") return PreconditionerBenchmarkVariant::surface_border_schur;
|
|
if (name == "gravity-schur-surface-border-schur")
|
|
return PreconditionerBenchmarkVariant::gravity_schur_surface_border_schur;
|
|
throw std::invalid_argument("Unknown preconditioner benchmark variant: " + name);
|
|
}
|
|
|
|
struct PreconditionerBenchmarkOptions final {
|
|
std::string variant{"independent"};
|
|
int restart{160};
|
|
int maximumIterations{200};
|
|
double relativeTolerance{0.03};
|
|
int schurCorrectionSweeps{1};
|
|
double schurRelaxation{0.5};
|
|
};
|
|
|
|
namespace preconditioner_benchmark_detail {
|
|
namespace pc = mean_field::preconditioning;
|
|
using Clock = std::chrono::steady_clock;
|
|
|
|
inline double Seconds(Clock::time_point start) {
|
|
return std::chrono::duration<double>(Clock::now() - start).count();
|
|
}
|
|
|
|
inline std::ofstream Open(const std::filesystem::path &path) {
|
|
std::ofstream stream(path);
|
|
stream.exceptions(std::ios::badbit | std::ios::failbit);
|
|
stream << std::setprecision(17);
|
|
return stream;
|
|
}
|
|
|
|
inline void Save(const std::filesystem::path &path, const mfem::Vector &values) {
|
|
auto stream = Open(path);
|
|
stream << values.Size() << '\n';
|
|
values.Print(stream, 1);
|
|
}
|
|
|
|
inline bool Identical(const mfem::Vector &left, const mfem::Vector &right) {
|
|
return left.Size() == right.Size() && (left.Size() == 0 ||
|
|
std::memcmp(left.GetData(), right.GetData(),
|
|
static_cast<std::size_t>(left.Size()) * sizeof(mfem::real_t)) == 0);
|
|
}
|
|
|
|
inline bool Finite(const mfem::Vector &values) {
|
|
for (int i = 0; i < values.Size(); ++i)
|
|
if (!std::isfinite(values(i))) return false;
|
|
return true;
|
|
}
|
|
|
|
inline double BlockNorm(const mfem::Vector &values, int offset, int size) {
|
|
long double sum = 0.0L;
|
|
for (int i = offset; i < offset + size; ++i)
|
|
sum += static_cast<long double>(values(i)) * values(i);
|
|
return std::sqrt(sum);
|
|
}
|
|
|
|
inline double Relative(double numerator, double denominator) {
|
|
return denominator > 0.0 ? numerator / denominator :
|
|
(numerator == 0.0 ? 0.0 : std::numeric_limits<double>::infinity());
|
|
}
|
|
|
|
template <pc::StellarStructureFactorizationPolicy StructurePolicy, bool MaterialLDU = false,
|
|
bool RefinedGravity = false, class Problem>
|
|
auto Prepare(const Problem &problem, const PreconditionerBenchmarkOptions &options) {
|
|
const auto defaults = pc::stellarStructureBlock(problem);
|
|
const auto material = [&]() {
|
|
if constexpr (MaterialLDU) {
|
|
const auto &original = defaults.materialSurfaceComponent();
|
|
return pc::materialSurfaceBlock(problem, original.materialBackend(), original.surfaceBackend(),
|
|
pc::ApproximateMaterialSurfaceLDU{}, original.surfaceSurrogate(), original.diagonalOptions());
|
|
} else {
|
|
return defaults.materialSurfaceComponent();
|
|
}
|
|
}();
|
|
auto gravity = [&]() {
|
|
if constexpr (RefinedGravity) {
|
|
const auto &original = defaults.gravityComponent();
|
|
return pc::GravityFieldBlock(original.massInverseBackend(), original.potentialSchurBackend(),
|
|
pc::GravityConsistentSchurLDU{.correctionSweeps = options.schurCorrectionSweeps,
|
|
.relaxation = options.schurRelaxation});
|
|
} else {
|
|
return defaults.gravityComponent();
|
|
}
|
|
}();
|
|
auto structure = pc::stellarStructureBlock(problem, material, std::move(gravity), StructurePolicy{});
|
|
// Each candidate needs its own completed border columns and scalar Schur;
|
|
// reusing the default border would confound the factorization comparison.
|
|
return pc::prepare(problem, pc::specificationBorderBlock(
|
|
problem, std::move(structure), pc::backend::DenseDirect{}));
|
|
}
|
|
|
|
// Instrumenting J changes its C++ address, whereas the normalization wrapper
|
|
// deliberately permits only its original Jacobian identity. Preserve that
|
|
// binding contract while allowing MFEM to apply the counted forwarding J.
|
|
class InstrumentedBinding final : public mfem::Solver {
|
|
public:
|
|
InstrumentedBinding(mfem::Solver &inverse, const mfem::Operator &original,
|
|
const mfem::Operator &instrumented)
|
|
: mfem::Solver(inverse.Height(), inverse.Width()), m_inverse(inverse),
|
|
m_original(original), m_instrumented(instrumented) {}
|
|
|
|
void SetOperator(const mfem::Operator &operation) override {
|
|
if (&operation != &m_instrumented)
|
|
throw std::invalid_argument("Benchmark preconditioner received a different counted Jacobian.");
|
|
m_inverse.SetOperator(m_original);
|
|
}
|
|
|
|
void Mult(const mfem::Vector &input, mfem::Vector &output) const override {
|
|
m_inverse.Mult(input, output);
|
|
}
|
|
|
|
private:
|
|
mfem::Solver &m_inverse;
|
|
const mfem::Operator &m_original;
|
|
const mfem::Operator &m_instrumented;
|
|
};
|
|
} // namespace preconditioner_benchmark_detail
|
|
|
|
template <class State>
|
|
void BenchmarkPreconditioners(State &state, const std::filesystem::path &output,
|
|
const PreconditionerBenchmarkOptions &options) {
|
|
namespace detail = preconditioner_benchmark_detail;
|
|
namespace pc = mean_field::preconditioning;
|
|
namespace solver = mean_field::solver;
|
|
const auto selected = ParsePreconditionerBenchmarkVariant(options.variant);
|
|
int ranks = 0;
|
|
if (MPI_Comm_size(state.Problem().GetCommunicator(), &ranks) != MPI_SUCCESS || ranks != 1)
|
|
throw std::invalid_argument("Frozen preconditioner benchmarks require one MPI rank.");
|
|
if (options.restart <= 0 || options.maximumIterations <= 0 ||
|
|
!std::isfinite(options.relativeTolerance) || options.relativeTolerance <= 0.0 ||
|
|
options.schurCorrectionSweeps < 0 || !std::isfinite(options.schurRelaxation) ||
|
|
options.schurRelaxation <= 0.0 || !state.normalizedOperator->IsPrepared())
|
|
throw std::invalid_argument("Frozen preconditioner benchmark controls or preparation are invalid.");
|
|
|
|
const mfem::Vector acceptedPhysical(state.AcceptedPhysicalState()),
|
|
acceptedNormalized(state.acceptedNormalizedState), acceptedResidual(state.acceptedNormalizedResidual),
|
|
correction(state.normalizedCorrection), frozenState(state.trialNormalizedState),
|
|
residual(state.trialNormalizedResidual), candidatePhysical(state.candidatePhysicalState);
|
|
const auto generation = state.Problem().GetPreparationGeneration();
|
|
const auto &normalization = state.normalizedOperator->GetNormalization();
|
|
const mfem::Vector stateFactors(normalization.StateFactors()), residualFactors(normalization.ResidualFactors());
|
|
mfem::Vector rhs(residual);
|
|
rhs *= -1.0;
|
|
const double rhsNorm = rhs.Norml2();
|
|
const double threshold = options.relativeTolerance * rhsNorm;
|
|
if (!detail::Finite(rhs) || !std::isfinite(rhsNorm) || rhsNorm <= 0.0 || !std::isfinite(threshold))
|
|
throw std::invalid_argument("A benchmark needs a finite, nonzero frozen Newton right-hand side.");
|
|
mfem::Vector physicalBase, physicalRhs;
|
|
state.normalizedOperator->DenormalizeState(frozenState, physicalBase);
|
|
normalization.DenormalizeResidual(rhs, physicalRhs);
|
|
detail::Save(output / "benchmark_physical_base.txt", physicalBase);
|
|
detail::Save(output / "benchmark_normalized_base.txt", frozenState);
|
|
detail::Save(output / "benchmark_normalized_rhs.txt", rhs);
|
|
detail::Save(output / "benchmark_state_factors.txt", stateFactors);
|
|
detail::Save(output / "benchmark_residual_factors.txt", residualFactors);
|
|
auto layout = detail::Open(output / "layout.csv");
|
|
layout << "kind,block,offset,size,scale\n";
|
|
for (const auto &block : state.Problem().GetManifest().residualBlocks())
|
|
layout << "residual," << block.stableId << ',' << block.offset << ',' << block.size << ',' << block.scale << '\n';
|
|
for (const auto &block : state.Problem().GetManifest().valueBlocks())
|
|
layout << "value," << block.stableId << ',' << block.offset << ',' << block.size << ',' << block.scale << '\n';
|
|
|
|
mfem::Vector direction(rhs.Size()), initialAction(rhs.Size());
|
|
for (int i = 0; i < direction.Size(); ++i)
|
|
direction(i) = std::sin(0.019 * (i + 1)) + 0.25 * std::cos(0.071 * (i + 1));
|
|
direction /= direction.Norml2();
|
|
state.normalizedOperator->Mult(direction, initialAction);
|
|
if (!detail::Finite(initialAction))
|
|
throw std::runtime_error("The frozen operator invariant probe is nonfinite.");
|
|
detail::Save(output / "benchmark_invariant_direction.txt", direction);
|
|
detail::Save(output / "benchmark_invariant_action.txt", initialAction);
|
|
|
|
auto invariants = detail::Open(output / "preconditioner_invariants.csv");
|
|
invariants << "role,metric,value\n";
|
|
const auto verify = [&](const std::string &role) {
|
|
bool valid = true;
|
|
const auto check = [&](const char *name, bool result) {
|
|
invariants << role << ',' << name << ',' << result << '\n';
|
|
valid = valid && result;
|
|
};
|
|
check("accepted_physical_bitwise_identical", detail::Identical(acceptedPhysical, state.AcceptedPhysicalState()));
|
|
check("accepted_normalized_bitwise_identical", detail::Identical(acceptedNormalized, state.acceptedNormalizedState));
|
|
check("accepted_residual_bitwise_identical", detail::Identical(acceptedResidual, state.acceptedNormalizedResidual));
|
|
check("correction_bitwise_identical", detail::Identical(correction, state.normalizedCorrection));
|
|
check("trial_state_bitwise_identical", detail::Identical(frozenState, state.trialNormalizedState));
|
|
check("trial_residual_bitwise_identical", detail::Identical(residual, state.trialNormalizedResidual));
|
|
check("candidate_physical_bitwise_identical", detail::Identical(candidatePhysical, state.candidatePhysicalState));
|
|
check("state_factors_bitwise_identical", detail::Identical(stateFactors, normalization.StateFactors()));
|
|
check("residual_factors_bitwise_identical", detail::Identical(residualFactors, normalization.ResidualFactors()));
|
|
check("preparation_generation_unchanged", generation == state.Problem().GetPreparationGeneration());
|
|
check("normalized_operator_current", state.normalizedOperator->IsPrepared());
|
|
mfem::Vector action(rhs.Size());
|
|
state.normalizedOperator->Mult(direction, action);
|
|
check("jacobian_action_bitwise_identical", detail::Identical(initialAction, action));
|
|
action -= initialAction;
|
|
invariants << role << ",jacobian_action_difference_norm," << action.Norml2() << '\n';
|
|
invariants.flush();
|
|
if (!valid)
|
|
throw std::runtime_error("Frozen preconditioner benchmark changed a context buffer, normalization, or operator.");
|
|
};
|
|
|
|
auto reference = detail::Open(output / "preconditioner_reference.csv");
|
|
reference << "metric,value\n";
|
|
auto start = detail::Clock::now();
|
|
const auto referenceRefresh = state.physicalInverse->Refresh();
|
|
reference << "context_inverse_refresh_seconds," << detail::Seconds(start)
|
|
<< "\ncontext_inverse_refresh_did_work," << referenceRefresh.DidAnyWork() << '\n';
|
|
mfem::Vector referenceCorrection(rhs.Size());
|
|
start = detail::Clock::now();
|
|
state.normalizedInverse->Mult(rhs, referenceCorrection);
|
|
reference << "context_inverse_application_seconds," << detail::Seconds(start)
|
|
<< "\ncontext_inverse_correction_norm," << referenceCorrection.Norml2() << '\n';
|
|
reference.flush();
|
|
if (!detail::Finite(referenceCorrection) || !std::isfinite(referenceCorrection.Norml2()))
|
|
throw std::runtime_error("The context reference inverse action or its norm is nonfinite.");
|
|
detail::Save(output / "context_reference_normalized_correction.txt", referenceCorrection);
|
|
verify("context_reference");
|
|
|
|
auto summary = detail::Open(output / "preconditioner_summary.csv");
|
|
summary << "role,variant,schur_correction_sweeps,schur_relaxation,restart,maximum_iterations,relative_tolerance,rhs_norm,absolute_target,"
|
|
"reported_initial_norm,reported_final_norm,reported_iterations,iterations,solver_converged,"
|
|
"true_residual_norm,true_relative,true_gate_pass,normalized_correction_norm,physical_correction_norm,"
|
|
"construction_seconds,unchanged_state_refresh_seconds,refresh_did_work,normalization_binding_seconds,"
|
|
"krylov_setup_seconds,krylov_seconds,true_check_seconds,setup_plus_solve_seconds,"
|
|
"setup_solve_and_true_check_seconds,measured_candidate_seconds,"
|
|
"jacobian_applications,jacobian_seconds,jacobian_maximum_seconds,"
|
|
"preconditioner_applications,preconditioner_seconds,preconditioner_maximum_seconds,"
|
|
"binding_setups,binding_refreshes,binding_setup_seconds,binding_refresh_seconds,"
|
|
"true_check_jacobian_applications,krylov_other_seconds\n";
|
|
auto coarseReport = detail::Open(output / "preconditioner_coarse.csv");
|
|
coarseReport << "role,variant,stage,metric,value\n";
|
|
auto schurReport = detail::Open(output / "preconditioner_surface_border_schur.csv");
|
|
schurReport << "role,variant,stage,metric,value\n";
|
|
auto schurCoordinates = detail::Open(output / "preconditioner_surface_border_coordinates.csv");
|
|
schurCoordinates << "role,variant,coarse_index,residual_index,correction_index\n";
|
|
auto schurMatrices = detail::Open(output / "preconditioner_surface_border_matrices.csv");
|
|
schurMatrices << "role,variant,row,column,uneliminated,schur\n";
|
|
auto work = detail::Open(output / "preconditioner_work.csv");
|
|
work << "role,variant,stage,metric,count\n";
|
|
auto actions = detail::Open(output / "preconditioner_actions.csv");
|
|
actions << "role,variant,correction_norm,jacobian_action_norm,defect_norm,relative_defect,"
|
|
"application_seconds,jacobian_seconds,context_action_relative_difference,context_equivalence_required\n";
|
|
auto blocks = detail::Open(output / "preconditioner_blocks.csv");
|
|
blocks << "role,variant,block,normalized_rhs_norm,normalized_defect_norm,physical_rhs_norm,"
|
|
"physical_defect_norm,block_relative,global_relative_contribution,squared_defect_fraction\n";
|
|
auto correctionBlocks = detail::Open(output / "preconditioner_correction_blocks.csv");
|
|
correctionBlocks << "role,variant,block,normalized_norm,physical_norm\n";
|
|
auto history = detail::Open(output / "preconditioner_history.csv");
|
|
history << "role,variant,iteration,reported_norm,final\n";
|
|
|
|
const auto run = [&]<pc::StellarStructureFactorizationPolicy Policy, bool MaterialLDU = false,
|
|
bool RefinedGravity = false, bool UseCoarse = false, bool UseSurfaceBorderSchur = false>(
|
|
const std::string &role, const std::string &variant) {
|
|
const auto candidateStart = detail::Clock::now();
|
|
auto physical = detail::Prepare<Policy, MaterialLDU, RefinedGravity>(state.Problem(), options);
|
|
static_assert(pc::PreparedStellarInverseFor<decltype(physical), decltype(state.Problem())>);
|
|
double constructionSeconds = detail::Seconds(candidateStart);
|
|
const auto recordWork = [&](const char *stage) {
|
|
const auto &structure = physical.GetGroupedPreconditioner().GetStructurePreconditioner();
|
|
const auto &stellar = structure.GetFactorization().GetStatistics();
|
|
const auto &gravity = structure.GetGravityPreconditioner().GetFactorization().GetStatistics();
|
|
const auto record = [&](const char *metric, std::uint64_t count) {
|
|
work << role << ',' << variant << ',' << stage << ',' << metric << ',' << count << '\n';
|
|
};
|
|
record("stellar_applications", stellar.applications);
|
|
record("material_inverse_applications", stellar.materialSurfaceInverseApplications);
|
|
record("gravity_inverse_applications", stellar.gravityInverseApplications);
|
|
record("material_to_gravity_applications", stellar.materialToGravityApplications);
|
|
record("gravity_to_material_applications", stellar.gravityToMaterialApplications);
|
|
record("gravity_mass_inverse_applications", gravity.massInverseApplications);
|
|
record("gravity_potential_schur_applications", gravity.potentialSchurApplications);
|
|
record("gravity_divergence_applications", gravity.divergenceApplications);
|
|
record("gravity_transpose_divergence_applications", gravity.transposeDivergenceApplications);
|
|
record("gravity_schur_refinement_sweeps", gravity.schurRefinementSweeps);
|
|
work.flush();
|
|
};
|
|
recordWork("constructed");
|
|
auto sectionStart = detail::Clock::now();
|
|
const auto refresh = physical.Refresh();
|
|
const double refreshSeconds = detail::Seconds(sectionStart);
|
|
sectionStart = detail::Clock::now();
|
|
auto normalized = state.normalizedOperator->MakeScaledPreconditioner(physical);
|
|
const double bindingSeconds = detail::Seconds(sectionStart);
|
|
using NormalizedOperation = std::remove_reference_t<decltype(*state.normalizedOperator)>;
|
|
using Coarse = PreparedCoupledShapeCoarse<NormalizedOperation, decltype(normalized)>;
|
|
using UnborderedStructure = NormalizedUnborderedStructureInverse<NormalizedOperation, decltype(physical)>;
|
|
using SurfaceBorderSchur = PreparedSurfaceBorderSchur<NormalizedOperation, UnborderedStructure>;
|
|
std::unique_ptr<Coarse> coarse;
|
|
std::unique_ptr<UnborderedStructure> unborderedStructure;
|
|
std::unique_ptr<SurfaceBorderSchur> surfaceBorderSchur;
|
|
static_assert(!(UseCoarse && UseSurfaceBorderSchur));
|
|
mfem::Solver *applicationInverse = &normalized;
|
|
const auto recordSurfaceBorderSchur = [&](const char *stage) {
|
|
if (!surfaceBorderSchur) return;
|
|
const auto &stats = surfaceBorderSchur->GetStatistics();
|
|
const auto metric = [&](const char *name, auto value) {
|
|
schurReport << role << ',' << variant << ',' << stage << ',' << name << ',' << value << '\n';
|
|
};
|
|
metric("status", static_cast<int>(stats.status));
|
|
metric("coarse_size", stats.coarseSize);
|
|
metric("reciprocal_condition", stats.reciprocalCondition);
|
|
metric("minimum_scaled_pivot", stats.minimumScaledPivot);
|
|
metric("setup_jacobian_applications", stats.setupOperatorApplications);
|
|
metric("setup_inverse_applications", stats.setupInverseApplications);
|
|
metric("setup_seconds", stats.lastSetupSeconds);
|
|
metric("setup_jacobian_seconds", stats.setupOperatorSeconds);
|
|
metric("setup_inverse_seconds", stats.setupInverseSeconds);
|
|
metric("persistent_vector_bytes", stats.persistentVectorBytes);
|
|
metric("unbordered_adapter_vector_bytes", sizeof(mfem::real_t) *
|
|
4 * static_cast<std::size_t>(state.normalizedOperator->Height()));
|
|
metric("applications", stats.applications);
|
|
metric("internal_jacobian_applications", stats.operatorApplications);
|
|
metric("internal_inverse_applications", stats.inverseApplications);
|
|
metric("internal_jacobian_seconds", stats.operatorSeconds);
|
|
metric("internal_inverse_seconds", stats.inverseSeconds);
|
|
metric("application_seconds", stats.applicationSeconds);
|
|
metric("failed_applications", stats.failedApplications);
|
|
schurReport.flush();
|
|
};
|
|
const auto recordCoarse = [&](const char *stage) {
|
|
if (!coarse) return;
|
|
const auto &stats = coarse->GetStatistics();
|
|
const auto metric = [&](const char *name, auto value) {
|
|
coarseReport << role << ',' << variant << ',' << stage << ',' << name << ',' << value << '\n';
|
|
};
|
|
metric("input_modes", stats.inputModes);
|
|
metric("independent_modes", stats.independentModes);
|
|
metric("active_modes", stats.activeModes);
|
|
metric("fallback", static_cast<int>(stats.fallback));
|
|
metric("minimum_overlap_lower_bound", stats.minimumOverlapLowerBound);
|
|
metric("reciprocal_overlap_condition", stats.reciprocalOverlapCondition);
|
|
metric("setup_jacobian_applications", stats.setupOperatorApplications);
|
|
metric("setup_inverse_applications", stats.setupInverseApplications);
|
|
metric("setup_seconds", stats.lastSetupSeconds);
|
|
metric("setup_jacobian_seconds", stats.lastSetupOperatorSeconds);
|
|
metric("setup_inverse_seconds", stats.lastSetupInverseSeconds);
|
|
metric("persistent_vector_bytes", stats.persistentVectorBytes);
|
|
metric("applications", stats.applications);
|
|
metric("inverse_applications", stats.inverseApplications);
|
|
metric("failed_applications", stats.failedApplications);
|
|
coarseReport.flush();
|
|
};
|
|
if constexpr (UseCoarse) {
|
|
sectionStart = detail::Clock::now();
|
|
const auto modes = MakeCoupledShapeModes(state.Problem(), physical);
|
|
const double basisSeconds = detail::Seconds(sectionStart);
|
|
recordWork("after_basis");
|
|
// Direct basis calls bypass the structure-factorization counters.
|
|
// The material inverse also counts its usual internal couplings.
|
|
for (const char *metric : {"basis_explicit_surface_to_material_applications",
|
|
"basis_material_inverse_applications",
|
|
"basis_material_to_gravity_applications",
|
|
"basis_gravity_inverse_applications"})
|
|
coarseReport << role << ',' << variant << ",constructed," << metric << ',' << modes.size() << '\n';
|
|
const auto coarseStart = detail::Clock::now();
|
|
coarse = std::make_unique<Coarse>(*state.normalizedOperator, normalized, modes,
|
|
state.Problem().GetCommunicator());
|
|
const double completeCoarseSeconds = basisSeconds + detail::Seconds(coarseStart);
|
|
recordWork("after_coarse_setup");
|
|
coarseReport << role << ',' << variant << ",constructed,physical_inverse_construction_seconds,"
|
|
<< constructionSeconds << '\n';
|
|
constructionSeconds += completeCoarseSeconds;
|
|
applicationInverse = coarse.get();
|
|
coarseReport << role << ',' << variant << ",constructed,basis_construction_seconds," << basisSeconds << '\n'
|
|
<< role << ',' << variant << ",constructed,complete_coarse_construction_seconds,"
|
|
<< completeCoarseSeconds << '\n';
|
|
recordCoarse("constructed");
|
|
// These extra checks validate the represented FEM coarse identity;
|
|
// they are diagnostics, separate from construction and solve timing.
|
|
double maximumRelativeError = 0.0;
|
|
for (const auto &source : coarse->GetSourceBasis()) {
|
|
mfem::Vector value(rhs.Size()), image(rhs.Size()), recovered(rhs.Size());
|
|
normalized.Mult(source, value);
|
|
state.normalizedOperator->Mult(value, image);
|
|
coarse->Mult(image, recovered);
|
|
recovered -= value;
|
|
const double valueNorm = value.Norml2(), imageNorm = image.Norml2(), errorNorm = recovered.Norml2();
|
|
if (!detail::Finite(value) || !detail::Finite(image) || !detail::Finite(recovered) ||
|
|
!std::isfinite(valueNorm) || valueNorm <= 0.0 || !std::isfinite(imageNorm) ||
|
|
!std::isfinite(errorNorm))
|
|
throw std::runtime_error("The actual FEM coarse identity produced a nonfinite or null reference.");
|
|
const double relativeError = errorNorm / valueNorm;
|
|
// std::max can discard NaN, so validate each observation first.
|
|
if (!std::isfinite(relativeError))
|
|
throw std::runtime_error("The actual FEM coarse relative identity error is nonfinite.");
|
|
maximumRelativeError = std::max(maximumRelativeError, relativeError);
|
|
}
|
|
coarseReport << role << ',' << variant << ",validation,exact_coarse_modes_checked,"
|
|
<< coarse->GetSourceBasis().size() << '\n'
|
|
<< role << ',' << variant << ",validation,maximum_exact_coarse_relative_error,"
|
|
<< maximumRelativeError << '\n';
|
|
coarseReport.flush();
|
|
if (!std::isfinite(maximumRelativeError) || maximumRelativeError > 1.e-8)
|
|
throw std::runtime_error("The actual FEM coarse inverse failed its exact-subspace check.");
|
|
}
|
|
|
|
if constexpr (UseSurfaceBorderSchur) {
|
|
sectionStart = detail::Clock::now();
|
|
unborderedStructure = std::make_unique<UnborderedStructure>(*state.normalizedOperator, physical);
|
|
surfaceBorderSchur = std::make_unique<SurfaceBorderSchur>(
|
|
*state.normalizedOperator, *unborderedStructure, MakeSurfaceBorderSchurCoordinates(physical),
|
|
state.Problem().GetCommunicator());
|
|
const double extraConstructionSeconds = detail::Seconds(sectionStart);
|
|
constructionSeconds += extraConstructionSeconds;
|
|
const auto &coordinates = surfaceBorderSchur->GetCoordinates();
|
|
const auto &schur = surfaceBorderSchur->GetSchurMatrix();
|
|
const auto &uneliminated = surfaceBorderSchur->GetUneliminatedCoarseMatrix();
|
|
for (int row = 0; row < schur.Height(); ++row) {
|
|
schurCoordinates << role << ',' << variant << ',' << row << ','
|
|
<< coordinates.residualIndices[row] << ','
|
|
<< coordinates.correctionIndices[row] << '\n';
|
|
for (int column = 0; column < schur.Width(); ++column)
|
|
schurMatrices << role << ',' << variant << ',' << row << ',' << column << ','
|
|
<< uneliminated(row, column) << ',' << schur(row, column) << '\n';
|
|
}
|
|
schurCoordinates.flush();
|
|
schurMatrices.flush();
|
|
recordSurfaceBorderSchur("constructed");
|
|
recordWork("after_surface_border_setup");
|
|
if (!surfaceBorderSchur->IsCurrent())
|
|
throw std::runtime_error("The measured surface/scalar Schur matrix is singular or ill-conditioned; no fallback was applied.");
|
|
applicationInverse = surfaceBorderSchur.get();
|
|
}
|
|
|
|
mfem::Vector inverseAction(rhs.Size()), inverseJacobianAction(rhs.Size());
|
|
sectionStart = detail::Clock::now();
|
|
applicationInverse->Mult(rhs, inverseAction);
|
|
const double applicationSeconds = detail::Seconds(sectionStart);
|
|
if (!detail::Finite(inverseAction) || !std::isfinite(inverseAction.Norml2()))
|
|
throw std::runtime_error("The diagnostic inverse action or its norm is nonfinite.");
|
|
detail::Save(output / (role + "_inverse_rhs_action.txt"), inverseAction);
|
|
sectionStart = detail::Clock::now();
|
|
state.normalizedOperator->Mult(inverseAction, inverseJacobianAction);
|
|
const double actionSeconds = detail::Seconds(sectionStart);
|
|
const double actionNorm = inverseJacobianAction.Norml2();
|
|
if (!detail::Finite(inverseJacobianAction) || !std::isfinite(actionNorm))
|
|
throw std::runtime_error("The diagnostic inverse Jacobian image or its norm is nonfinite.");
|
|
inverseJacobianAction -= rhs;
|
|
if (surfaceBorderSchur) {
|
|
double squared = 0.0;
|
|
for (int row : surfaceBorderSchur->GetCoordinates().residualIndices)
|
|
squared += inverseJacobianAction(row) * inverseJacobianAction(row);
|
|
const double coarseRelative = std::sqrt(squared) / rhsNorm;
|
|
schurReport << role << ',' << variant
|
|
<< ",validation,coarse_equation_defect_relative_to_full_rhs," << coarseRelative << '\n';
|
|
schurReport.flush();
|
|
if (!std::isfinite(coarseRelative) || coarseRelative > 1.e-8)
|
|
throw std::runtime_error("The Schur inverse failed the authoritative-J coarse-equation check.");
|
|
}
|
|
mfem::Vector referenceDifference(inverseAction);
|
|
referenceDifference -= referenceCorrection;
|
|
const double referenceRelative = detail::Relative(referenceDifference.Norml2(),
|
|
std::max(inverseAction.Norml2(), referenceCorrection.Norml2()));
|
|
const bool requireEquivalent = std::same_as<Policy, pc::IndependentStellarSubsystems> && !MaterialLDU && !UseCoarse &&
|
|
!UseSurfaceBorderSchur &&
|
|
(!RefinedGravity || options.schurCorrectionSweeps == 0);
|
|
actions << role << ',' << variant << ',' << inverseAction.Norml2() << ',' << actionNorm << ','
|
|
<< inverseJacobianAction.Norml2() << ',' << inverseJacobianAction.Norml2() / rhsNorm << ','
|
|
<< applicationSeconds << ',' << actionSeconds << ',' << referenceRelative << ','
|
|
<< requireEquivalent << '\n';
|
|
actions.flush();
|
|
if (requireEquivalent && (!std::isfinite(referenceRelative) ||
|
|
referenceRelative > 256.0 * std::numeric_limits<double>::epsilon()))
|
|
throw std::runtime_error("Fresh independent inverse differs from the context's default inverse.");
|
|
|
|
sectionStart = detail::Clock::now();
|
|
solver::InstrumentedOperator countedJacobian(*state.normalizedOperator);
|
|
detail::InstrumentedBinding binding(*applicationInverse, *state.normalizedOperator, countedJacobian);
|
|
solver::InstrumentedPreconditioner countedInverse(binding);
|
|
solver::ResidualHistoryMonitor monitor;
|
|
mfem::FGMRESSolver krylov(state.Problem().GetCommunicator());
|
|
krylov.SetPreconditioner(countedInverse);
|
|
krylov.SetOperator(countedJacobian);
|
|
krylov.SetMonitor(monitor);
|
|
krylov.SetKDim(options.restart);
|
|
krylov.SetRelTol(0.0);
|
|
krylov.SetAbsTol(threshold);
|
|
krylov.SetMaxIter(options.maximumIterations);
|
|
krylov.SetPrintLevel(-1);
|
|
krylov.iterative_mode = false;
|
|
mfem::Vector solution(rhs.Size()), defect(rhs.Size());
|
|
solution = 0.0;
|
|
const double krylovSetupSeconds = detail::Seconds(sectionStart);
|
|
recordWork("before_krylov");
|
|
recordCoarse("before_krylov");
|
|
recordSurfaceBorderSchur("before_krylov");
|
|
countedJacobian.ResetStatistics();
|
|
countedInverse.ResetStatistics();
|
|
sectionStart = detail::Clock::now();
|
|
krylov.Mult(rhs, solution);
|
|
const double krylovSeconds = detail::Seconds(sectionStart);
|
|
recordWork("after_krylov");
|
|
recordCoarse("after_krylov");
|
|
recordSurfaceBorderSchur("after_krylov");
|
|
const auto jacobianStats = countedJacobian.GetStatistics();
|
|
const auto inverseStats = countedInverse.GetStatistics();
|
|
const auto lifecycle = countedInverse.GetLifecycleStatistics();
|
|
|
|
// The recurrence's convergence flag is reported, but only this fresh
|
|
// original-operator application determines the linear accuracy gate.
|
|
sectionStart = detail::Clock::now();
|
|
state.normalizedOperator->Mult(solution, defect);
|
|
defect -= rhs;
|
|
const double defectNorm = defect.Norml2();
|
|
const bool finite = detail::Finite(solution) && detail::Finite(defect) && std::isfinite(defectNorm);
|
|
const bool gate = finite && defectNorm <= threshold;
|
|
const double trueCheckSeconds = detail::Seconds(sectionStart);
|
|
const double candidateSeconds = detail::Seconds(candidateStart);
|
|
mfem::Vector physicalSolution, physicalDefect;
|
|
normalization.DenormalizeState(solution, physicalSolution);
|
|
normalization.DenormalizeResidual(defect, physicalDefect);
|
|
const double setupAndSolve = constructionSeconds + bindingSeconds + krylovSetupSeconds + krylovSeconds;
|
|
summary << role << ',' << variant << ',' << (RefinedGravity ? options.schurCorrectionSweeps : 0) << ','
|
|
<< (RefinedGravity ? options.schurRelaxation : 0.0) << ',' << options.restart << ',' << options.maximumIterations << ','
|
|
<< options.relativeTolerance << ',' << rhsNorm << ',' << threshold << ','
|
|
<< krylov.GetInitialNorm() << ',' << krylov.GetFinalNorm() << ',' << krylov.GetNumIterations() << ','
|
|
<< inverseStats.applications << ',' << krylov.GetConverged() << ',' << defectNorm << ','
|
|
<< defectNorm / rhsNorm << ',' << gate << ',' << solution.Norml2() << ',' << physicalSolution.Norml2() << ','
|
|
<< constructionSeconds << ',' << refreshSeconds << ',' << refresh.DidAnyWork() << ',' << bindingSeconds << ','
|
|
<< krylovSetupSeconds << ',' << krylovSeconds << ',' << trueCheckSeconds << ',' << setupAndSolve << ','
|
|
<< setupAndSolve + trueCheckSeconds << ',' << candidateSeconds << ','
|
|
<< jacobianStats.applications << ',' << jacobianStats.totalSeconds << ',' << jacobianStats.maximumSeconds << ','
|
|
<< inverseStats.applications << ',' << inverseStats.totalSeconds << ',' << inverseStats.maximumSeconds << ','
|
|
<< lifecycle.setups << ',' << lifecycle.refreshes << ',' << lifecycle.setupSeconds << ','
|
|
<< lifecycle.refreshSeconds << ",1," <<
|
|
std::max(0.0, krylovSeconds - jacobianStats.totalSeconds - inverseStats.totalSeconds) << '\n';
|
|
summary.flush();
|
|
for (const auto &block : state.Problem().GetManifest().residualBlocks()) {
|
|
const double b = detail::BlockNorm(rhs, block.offset, block.size);
|
|
const double r = detail::BlockNorm(defect, block.offset, block.size);
|
|
blocks << role << ',' << variant << ',' << block.stableId << ',' << b << ',' << r << ','
|
|
<< detail::BlockNorm(physicalRhs, block.offset, block.size) << ','
|
|
<< detail::BlockNorm(physicalDefect, block.offset, block.size) << ','
|
|
<< detail::Relative(r, b) << ',' << r / rhsNorm << ','
|
|
<< (defectNorm > 0.0 ? std::pow(r / defectNorm, 2) : 0.0) << '\n';
|
|
}
|
|
for (const auto &block : state.Problem().GetManifest().valueBlocks())
|
|
correctionBlocks << role << ',' << variant << ',' << block.stableId << ','
|
|
<< detail::BlockNorm(solution, block.offset, block.size) << ','
|
|
<< detail::BlockNorm(physicalSolution, block.offset, block.size) << '\n';
|
|
for (const auto &sample : monitor.GetHistory())
|
|
history << role << ',' << variant << ',' << sample.iteration << ',' << sample.reportedNorm << ','
|
|
<< sample.final << '\n';
|
|
blocks.flush();
|
|
correctionBlocks.flush();
|
|
history.flush();
|
|
detail::Save(output / (role + "_normalized_correction.txt"), solution);
|
|
detail::Save(output / (role + "_physical_correction.txt"), physicalSolution);
|
|
detail::Save(output / (role + "_normalized_linear_defect.txt"), defect);
|
|
detail::Save(output / (role + "_physical_linear_defect.txt"), physicalDefect);
|
|
verify(role);
|
|
std::cout << "Preconditioner " << role << '/' << variant << " iterations=" << inverseStats.applications
|
|
<< " true_relative=" << defectNorm / rhsNorm << " true_gate=" << gate
|
|
<< " krylov_seconds=" << krylovSeconds
|
|
<< " setup_solve_and_true_check_seconds=" << setupAndSolve + trueCheckSeconds << std::endl;
|
|
};
|
|
|
|
run.template operator()<pc::IndependentStellarSubsystems>("baseline", "independent");
|
|
switch (selected) {
|
|
case PreconditionerBenchmarkVariant::independent:
|
|
run.template operator()<pc::IndependentStellarSubsystems>("candidate", options.variant); break;
|
|
case PreconditionerBenchmarkVariant::material_then_gravity:
|
|
run.template operator()<pc::MaterialThenGravityTriangular>("candidate", options.variant); break;
|
|
case PreconditionerBenchmarkVariant::gravity_then_material:
|
|
run.template operator()<pc::GravityThenMaterialTriangular>("candidate", options.variant); break;
|
|
case PreconditionerBenchmarkVariant::stellar_ldu:
|
|
run.template operator()<pc::ApproximateStellarBlockLDU>("candidate", options.variant); break;
|
|
case PreconditionerBenchmarkVariant::material_ldu:
|
|
run.template operator()<pc::IndependentStellarSubsystems, true>("candidate", options.variant); break;
|
|
case PreconditionerBenchmarkVariant::gravity_schur:
|
|
run.template operator()<pc::IndependentStellarSubsystems, false, true>("candidate", options.variant); break;
|
|
case PreconditionerBenchmarkVariant::stellar_ldu_gravity_schur:
|
|
run.template operator()<pc::ApproximateStellarBlockLDU, false, true>("candidate", options.variant); break;
|
|
case PreconditionerBenchmarkVariant::coupled_shape_coarse:
|
|
run.template operator()<pc::IndependentStellarSubsystems, false, false, true>("candidate", options.variant); break;
|
|
case PreconditionerBenchmarkVariant::gravity_schur_coarse:
|
|
run.template operator()<pc::IndependentStellarSubsystems, false, true, true>("candidate", options.variant); break;
|
|
case PreconditionerBenchmarkVariant::surface_border_schur:
|
|
run.template operator()<pc::IndependentStellarSubsystems, false, false, false, true>("candidate", options.variant); break;
|
|
case PreconditionerBenchmarkVariant::gravity_schur_surface_border_schur:
|
|
run.template operator()<pc::IndependentStellarSubsystems, false, true, false, true>("candidate", options.variant); break;
|
|
}
|
|
auto scope = detail::Open(output / "preconditioner_scope.txt");
|
|
scope << "Both solves use fresh physical inverses, separately completed scalar borders, the same physical J,\n"
|
|
<< "the same normalization, the same saved RHS=-F, and zero initial corrections.\n"
|
|
<< "The fresh independent inverse must match the existing context inverse on that RHS.\n"
|
|
<< "construction_seconds includes physical inverse and any selected coarse/basis construction.\n"
|
|
<< "Coarse exact-subspace validation is diagnostic work outside setup/solve totals.\n"
|
|
<< "An inactive coarse fallback has exact_coarse_modes_checked=0; its zero maximum error is vacuous.\n"
|
|
<< "Stellar-level inverse/cross counters exclude direct basis calls, recorded separately; gravity leaf counters include them.\n"
|
|
<< "unchanged_state_refresh_seconds measures an immediate Refresh; it is not a changed-state refresh benchmark.\n"
|
|
<< "Each candidate receives one diagnostic inverse/Jv application before the timed cold Krylov solve.\n"
|
|
<< "Application timers/counts in the summary include only Krylov work; the independent true check is separate.\n"
|
|
<< "Surface/scalar Schur setup is included in construction; its internal J applications are separately counted.\n"
|
|
<< "For that inverse, outer Krylov J counts exclude internal Schur J work; inverse seconds include it.\n"
|
|
<< "Schur status codes are 0=unprepared, 1=ready, 2=singular, 3=ill-conditioned; rejected setup has no fallback.\n"
|
|
<< "The Schur partition uses every native surface parameter and generated scalar block, including a central border if present.\n"
|
|
<< "Its field inverse restricts the unbordered structure inverse, not an exact field elimination.\n"
|
|
<< "Existing scalar-border construction is retained in setup cost but its inverse is not applied by the new wrapper.\n"
|
|
<< "The unbordered adapter is frozen at construction; rebuild it and the Schur wrapper for a new Jacobian.\n"
|
|
<< "setup_plus_solve_seconds includes construction, normalization binding, Krylov setup, and Krylov solve.\n"
|
|
<< "setup_solve_and_true_check_seconds adds the separately measured original-J true-residual gate.\n"
|
|
<< "measured_candidate_seconds also includes diagnostic refresh, coarse exactness checks, and inverse/Jv comparison.\n"
|
|
<< "Common FEM/context construction and frozen-state preparation occur before this helper.\n"
|
|
<< "The helper performs no trial preparation or nonlinear update; WithState restores its accepted seed afterward.\n"
|
|
<< "Block norms are coefficient norms in the stated normalized/physical coordinates, not spatial error norms.\n"
|
|
<< "A linear gate pass does not establish nonlinear acceptance, EOS admissibility, or virial qualification.\n";
|
|
}
|
|
} // namespace experiment::plateau_diagnostics
|