From 37eb5815bfed32e828f9c5d1d3defae1a7456ba1 Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Fri, 9 Oct 2026 02:12:01 +0900 Subject: [PATCH] feat(analysis): execute solid generalized frequency pencil --- .../fesa/analysis/solid_frequency_procedure.h | 18 + include/fesa/analysis/solid_frequency_state.h | 54 +++ include/fesa/elements/hex8_element.h | 3 + .../solvers/eigen/mkl_feast_eigen_solver.h | 24 +- .../analysis/solid_frequency_procedure.cpp | 351 ++++++++++++++++++ src/fesa/elements/hex8_element.cpp | 7 + .../solvers/eigen/mkl_feast_eigen_solver.cpp | 298 +++++++++++---- .../solid_frequency_procedure_test.cpp | 278 ++++++++++++++ .../eigen/mkl_feast_eigen_solver_test.cpp | 155 ++++++++ 9 files changed, 1110 insertions(+), 78 deletions(-) create mode 100644 include/fesa/analysis/solid_frequency_procedure.h create mode 100644 include/fesa/analysis/solid_frequency_state.h create mode 100644 src/fesa/analysis/solid_frequency_procedure.cpp create mode 100644 tests/unit/analysis/solid_frequency_procedure_test.cpp diff --git a/include/fesa/analysis/solid_frequency_procedure.h b/include/fesa/analysis/solid_frequency_procedure.h new file mode 100644 index 0000000..d779ace --- /dev/null +++ b/include/fesa/analysis/solid_frequency_procedure.h @@ -0,0 +1,18 @@ +#ifndef FESA_ANALYSIS_SOLID_FREQUENCY_PROCEDURE_H_ +#define FESA_ANALYSIS_SOLID_FREQUENCY_PROCEDURE_H_ + +#include "fesa/analysis/solid_frequency_state.h" + +namespace fesa { +class Domain; +/// @brief Executes unloaded reference-state C3D8 generalized frequency +/// analysis. +class SolidFrequencyProcedureExecutor final { + public: + /// @brief Returns lowest requested modes after validating every returned + /// pair. + /// @note Domain is immutable and must remain alive during execution. + Result Execute(const Domain& domain) const; +}; +} // namespace fesa +#endif // FESA_ANALYSIS_SOLID_FREQUENCY_PROCEDURE_H_ diff --git a/include/fesa/analysis/solid_frequency_state.h b/include/fesa/analysis/solid_frequency_state.h new file mode 100644 index 0000000..51eb5a7 --- /dev/null +++ b/include/fesa/analysis/solid_frequency_state.h @@ -0,0 +1,54 @@ +#ifndef FESA_ANALYSIS_SOLID_FREQUENCY_STATE_H_ +#define FESA_ANALYSIS_SOLID_FREQUENCY_STATE_H_ + +#include +#include +#include +#include +#include + +#include "fesa/math/matrix.h" +#include "fesa/solvers/eigen/mkl_feast_eigen_solver.h" + +namespace fesa { +/// @brief Identifies a source-ordered physical Solid eigenmode. +struct SolidEigenMode { + std::size_t mode_id; + double eigenvalue; + std::uint8_t classification; + std::size_t cluster_id; + std::optional frequency; +}; +/// @brief Stores original free K/M validation for every backend pair. +struct SolidEigenpairCheck { + std::size_t sorted_rank; + std::size_t solver_column; + double eigenvalue; + double modal_mass; + double rayleigh_quotient; + double eta_r; + double eta_norm; + double eta_energy; +}; +/// @brief Owns a completed modal candidate, independent of response frames. +/// @note Vectors are requested-mode by full source node/UX,UY,UZ order. +struct SolidFrequencyState { + explicit SolidFrequencyState(FeastEigenResult result) + : backend{std::move(result)} {} + std::size_t full_dof_count{0U}; + std::size_t free_dof_count{0U}; + double classification_scale{0.0}; + double mass_ideal_lower_bound{0.0}; + double mass_roundoff_guard{0.0}; + double mass_lower_bound{0.0}; + double eta_gram{0.0}; + double max_eta_orth{0.0}; + std::vector constraint_mask; + std::vector prescribed_displacement; + std::vector modes; + Matrix vectors{0U, 0U}; + std::vector eigenpair_checks; + FeastEigenResult backend; +}; +} // namespace fesa +#endif // FESA_ANALYSIS_SOLID_FREQUENCY_STATE_H_ diff --git a/include/fesa/elements/hex8_element.h b/include/fesa/elements/hex8_element.h index 29e6866..a90476f 100644 --- a/include/fesa/elements/hex8_element.h +++ b/include/fesa/elements/hex8_element.h @@ -47,6 +47,9 @@ class Hex8Element final : public Element { Result EvaluateNonlinear(const Vector& displacement) const; /// @brief Returns the constant reference consistent mass; density required. Result ComputeMass() const; + /// @brief Returns immutable positive reference Gauss volume weights. + /// @note Order is the kernel's fixed eight-point integration order. + std::array ReferenceWeights() const noexcept; private: /// @brief Caches shape, reference gradients and positive volume weights. diff --git a/include/fesa/solvers/eigen/mkl_feast_eigen_solver.h b/include/fesa/solvers/eigen/mkl_feast_eigen_solver.h index 7efbcc5..3e13447 100644 --- a/include/fesa/solvers/eigen/mkl_feast_eigen_solver.h +++ b/include/fesa/solvers/eigen/mkl_feast_eigen_solver.h @@ -58,7 +58,7 @@ struct FeastEigenResult { std::string upper_search_termination; }; -/// @brief Solves a whitened real-symmetric eigenproblem with oneMKL FEAST. +/// @brief Adapts standard and symmetric-definite sparse pencils to FEAST. class MklFeastEigenSolver final { public: /// @brief Computes all validated interval pairs for a requested lowest-mode @@ -76,6 +76,28 @@ class MklFeastEigenSolver final { Result Solve(const SparseMatrix& whitened_operator, std::size_t requested_modes, double zero_operator_scale) const; + /// @brief Solves the original sparse symmetric definite K/M pencil. + /// @param mass_lower_bound Certified positive lower spectral bound of M, + /// including floating-point assembly uncertainty; never a diagonal estimate. + /// @param classification_scale Positive max |Kii|/Mii in physical units. + /// @note Eigenvectors are mass-normalized. For this API the historical + /// whitened_residuals field contains original-pencil eta_r; no whitening + /// is performed. Intervals/eigenvalues in normalized backend evidence use + /// classification_scale, while eigenvalues are returned in physical units. + Result SolveGeneralized(const SparseMatrix& stiffness, + const SparseMatrix& mass, + std::size_t requested_modes, + double mass_lower_bound, + double classification_scale) const; + + private: + /// @brief Shares deterministic count/search controls without changing + /// pencils. + Result SolveImpl(const SparseMatrix& stiffness, + const SparseMatrix* mass, + std::size_t requested_modes, + double physical_scale, + double mass_lower_bound) const; }; } // namespace fesa diff --git a/src/fesa/analysis/solid_frequency_procedure.cpp b/src/fesa/analysis/solid_frequency_procedure.cpp new file mode 100644 index 0000000..556f5e7 --- /dev/null +++ b/src/fesa/analysis/solid_frequency_procedure.cpp @@ -0,0 +1,351 @@ +#include "fesa/analysis/solid_frequency_procedure.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "fesa/analysis/analysis_model.h" +#include "fesa/assembly/parallel_for.h" +#include "fesa/assembly/sparse_assembler.h" +#include "fesa/constraints/essential_constraint_policy.h" +#include "fesa/elements/element_factory.h" +#include "fesa/elements/hex8_element.h" +#include "fesa/fem/dof_manager.h" +#include "fesa/math/sparse_matrix.h" +#include "fesa/math/vector.h" +#include "fesa/model/domain.h" + +namespace fesa { +namespace { +/// @brief Reports a frequency failure without publishing partial state. +Status Failure(const std::string& code, const std::string& message) { + return Status::Failure( + FailureCategory::kSolver, + {{Severity::kError, code, {}, "FREQUENCY", {}, message}}); +} +/// @brief Rounds nonnegative arithmetic toward conservative bounds. +double Up(double value) { + return std::nextafter(value, (std::numeric_limits::infinity)()); +} +/// @brief Rounds a lower bound downward, including underflow uncertainty. +double Down(double value) { + return std::nextafter(value, -(std::numeric_limits::infinity)()); +} +/// @brief Uses robust Frobenius accumulation on original CSR storage. +double Frobenius(const SparseMatrix& matrix) { + double norm = 0.0; + for (const double value : matrix.Values()) norm = std::hypot(norm, value); + return norm; +} +/// @brief Evaluates exact-zero denominator semantics for modal checks. +double Ratio(double numerator, double denominator) { + if (denominator == 0.0) + return numerator == 0.0 ? 0.0 : (std::numeric_limits::infinity)(); + return numerator / denominator; +} +/// @brief Builds an immutable candidate; all returns precede public commit. +Result ExecuteCandidate(const Domain& domain) { + const auto* frequency = domain.Frequency(); + if (!domain.IsSolid() || frequency == nullptr || + frequency->requested_mode_count <= 0) + return Result::Failure(Failure( + "invalid-solid-frequency-procedure", "Solid frequency step required.")); + auto model = AnalysisModel::Create(domain); + if (!model.HasValue()) + return Result::Failure(model.GetStatus()); + if (!model.Value().Step().Loads().empty() || + !model.Value().ActiveSurfacePressureLoads().empty()) + return Result::Failure( + Failure("loaded-solid-frequency", + "Frequency requires unloaded reference state.")); + std::vector density(domain.Materials().Size(), 0.0); + for (const auto& item : domain.Densities()) { + if (item.material_index >= density.size() || + density[item.material_index] != 0.0 || !std::isfinite(item.rho) || + !(item.rho > 0.0)) + return Result::Failure( + Failure("invalid-material-density", + "Density must be unique finite positive.")); + density[item.material_index] = item.rho; + } + ElementFactory factory; + std::vector> owned; + ElementView view; + for (const auto index : model.Value().ActiveElements()) { + const auto& definition = domain.Elements()[index]; + if (definition.MaterialIndex() >= density.size() || + !(density[definition.MaterialIndex()] > 0.0)) + return Result::Failure( + Failure("missing-frequency-density", + "Every Solid material requires density.")); + auto element = factory.Create(definition, domain); + if (!element.HasValue()) + return Result::Failure(element.GetStatus()); + if (dynamic_cast(element.Value().get()) == nullptr) + return Result::Failure( + Failure("unsupported-frequency-element", "Only C3D8 is supported.")); + view.emplace_back(*element.Value()); + owned.push_back(std::move(element.Value())); + } + DofManager dofs; + const Status mapping = dofs.Build(model.Value(), view); + if (!mapping.IsOk()) return Result::Failure(mapping); + if (dofs.Profile() != DofProfile::kSolidStructural || + dofs.ComponentsPerNode() != 3U) + return Result::Failure( + Failure("invalid-solid-frequency-map", + "Frequency requires actual 3N mapping.")); + for (std::size_t i = 0; i < dofs.PrescribedValues().Size(); ++i) + if (dofs.PrescribedValues()[i] != 0.0) + return Result::Failure(Failure( + "unsupported-eigen-boundary", "Frequency requires homogeneous BC.")); + const std::size_t requested = + static_cast(frequency->requested_mode_count); + if (dofs.FreeDofCount() == 0U || requested > dofs.FreeDofCount()) + return Result::Failure( + Failure("insufficient-solid-frequency-modes", + "Requested modes exceed free dimension.")); + const SerialParallelFor serial; + auto stiffness = SparseAssembler::Assemble(view, dofs, serial); + if (!stiffness.HasValue()) + return Result::Failure(stiffness.GetStatus()); + std::vector mass_coo; + std::vector lower_by_dof(dofs.FullDofCount(), 0.0); + std::vector magnitude_by_dof(dofs.FullDofCount(), 0.0); + std::vector additions_by_dof(dofs.FullDofCount(), 0U); + constexpr double eps = std::numeric_limits::epsilon(); + for (std::size_t e = 0; e < owned.size(); ++e) { + const auto& hex = static_cast(*owned[e]); + auto mass = hex.ComputeMass(); + if (!mass.HasValue()) + return Result::Failure(mass.GetStatus()); + auto scatter = dofs.ElementScatter(hex.DofLayout()); + if (!scatter.HasValue()) + return Result::Failure(scatter.GetStatus()); + const auto weights = hex.ReferenceWeights(); + const double min_weight = *std::min_element(weights.begin(), weights.end()); + const auto element_index = model.Value().ActiveElements()[e]; + const double rho = + density[domain.Elements()[element_index].MaterialIndex()]; + const double local_lower = Down(Down(rho * min_weight) / 27.0); + // Relative gamma bounds require normal nonzero products. Fixed Gauss + // shapes exceed 1/128, so rho/16384 bounds the smallest density product. + if (!std::isfinite(local_lower) || + local_lower < std::numeric_limits::min() || + rho / 16384.0 < std::numeric_limits::min() || + min_weight < std::numeric_limits::min()) + return Result::Failure( + Failure("solid-mass-bound-unavailable", + "Local reference mass bound is unrepresentable.")); + // The eight fixed trilinear Gauss shapes have sigma_min(H)^2=1/27. + // 512 eps bounds their fixed evaluation, products and eight-term sums. + // Every mass term is positive; stable incidence scatter preserves the + // bound. + for (std::size_t i = 0; i < 24U; ++i) { + const auto row = scatter.Value()[i]; + lower_by_dof[row] = Down(lower_by_dof[row] + local_lower); + double row_sum = 0.0; + for (std::size_t j = 0; j < 24U; ++j) { + const double value = mass.Value().values(i, j); + if (value != 0.0 && + std::abs(value) < std::numeric_limits::min()) + return Result::Failure(Failure( + "solid-mass-bound-unavailable", + "Subnormal local mass invalidates relative roundoff bounds.")); + row_sum = Up(row_sum + std::abs(value)); + mass_coo.push_back({row, scatter.Value()[j], value, e, i * 24U + j}); + } + magnitude_by_dof[row] = Up(magnitude_by_dof[row] + row_sum); + ++additions_by_dof[row]; + } + } + double ideal_lower = (std::numeric_limits::infinity)(); + double guard = 0.0; + for (const auto row : dofs.FreeDofs()) { + ideal_lower = (std::min)(ideal_lower, lower_by_dof[row]); + const double count_eps = + Up((512.0 + static_cast(additions_by_dof[row])) * eps); + if (!(count_eps < 0.5)) + return Result::Failure( + Failure("solid-mass-bound-unavailable", + "Mass reduction exceeds roundoff budget.")); + const double gamma = Up(count_eps / (1.0 - count_eps)); + // Convert the forward bound relative to the exact sum to an error bound + // relative to the available computed magnitude. Row infinity bounds also + // bound the symmetric perturbation's spectral norm under restriction. + const double computed_gamma = Up(gamma / (1.0 - gamma)); + guard = (std::max)(guard, Up(computed_gamma * magnitude_by_dof[row])); + } + const double mass_lower = Down(ideal_lower - guard); + if (!std::isfinite(mass_lower) || !(mass_lower > 0.0) || + !(guard < 0.01 * ideal_lower)) + return Result::Failure( + Failure("solid-mass-bound-unavailable", + "Roundoff prevents a useful SPD mass enclosure.")); + auto mass = + SparseMatrix::FromCoo(dofs.FullDofCount(), dofs.FullDofCount(), + std::move(mass_coo), dofs.GetSparsePattern()); + if (!mass.HasValue()) + return Result::Failure(mass.GetStatus()); + const EssentialConstraintPolicy partition; + auto k_blocks = partition.Partition(stiffness.Value(), dofs); + auto m_blocks = partition.Partition(mass.Value(), dofs); + if (!k_blocks.HasValue()) + return Result::Failure(k_blocks.GetStatus()); + if (!m_blocks.HasValue()) + return Result::Failure(m_blocks.GetStatus()); + const auto& k = k_blocks.Value().k_ff; + const auto& m = m_blocks.Value().k_ff; + double scale = 0.0; + for (std::size_t row = 0; row < k.Rows(); ++row) { + double kd = 0.0; + double md = 0.0; + for (std::size_t p = k.RowOffsets()[row]; p < k.RowOffsets()[row + 1U]; ++p) + if (k.ColumnIndices()[p] == row) kd = k.Values()[p]; + for (std::size_t p = m.RowOffsets()[row]; p < m.RowOffsets()[row + 1U]; ++p) + if (m.ColumnIndices()[p] == row) md = m.Values()[p]; + if (!(md > 0.0) || !std::isfinite(std::abs(kd) / md)) + return Result::Failure( + Failure("invalid-solid-eigen-scale", + "Mass diagonal or classification scale is invalid.")); + scale = (std::max)(scale, std::abs(kd) / md); + } + if (!(scale > 0.0) || !std::isfinite(scale)) + return Result::Failure( + Failure("invalid-solid-eigen-scale", + "A positive finite classification scale is required.")); + auto backend = MklFeastEigenSolver{}.SolveGeneralized(k, m, requested, + mass_lower, scale); + if (!backend.HasValue()) + return Result::Failure(backend.GetStatus()); + SolidFrequencyState state(std::move(backend.Value())); + state.full_dof_count = dofs.FullDofCount(); + state.free_dof_count = dofs.FreeDofCount(); + state.classification_scale = scale; + state.mass_ideal_lower_bound = ideal_lower; + state.mass_roundoff_guard = guard; + state.mass_lower_bound = mass_lower; + state.constraint_mask.assign(dofs.FullDofCount(), 0U); + state.prescribed_displacement.assign(dofs.FullDofCount(), 0.0); + for (const auto row : dofs.ConstrainedDofs()) state.constraint_mask[row] = 1U; + const std::size_t count = state.backend.eigenvalues.size(); + if (count < requested || state.backend.eigenvectors.Rows() != k.Rows() || + state.backend.eigenvectors.Columns() != count) + return Result::Failure( + Failure("insufficient-solid-frequency-pairs", + "Backend returned insufficient complete pairs.")); + const double norm_k = Frobenius(k); + const double norm_m = Frobenius(m); + std::vector phi; + std::vector mass_phi; + phi.reserve(count); + mass_phi.reserve(count); + for (std::size_t c = 0; c < count; ++c) { + Vector vector(k.Rows()); + for (std::size_t row = 0; row < k.Rows(); ++row) + vector[row] = state.backend.eigenvectors(row, c); + Vector mv = m.Multiply(vector); + const double raw_mass = vector.Dot(mv); + if (!(raw_mass > 0.0) || !std::isfinite(raw_mass)) + return Result::Failure( + Failure("solid-eigenpair-validation-failed", + "Modal mass is not finite positive.")); + vector.Scale(1.0 / std::sqrt(raw_mass)); + mv = m.Multiply(vector); + const Vector kv = k.Multiply(vector); + const double lambda = state.backend.eigenvalues[c]; + const double modal_mass = vector.Dot(mv); + const double energy = vector.Dot(kv); + Vector residual = kv; + residual.Axpy(-lambda, mv); + const double norm_phi = vector.Norm(); + const double denominator = (norm_k + std::abs(lambda) * norm_m) * norm_phi; + const double eta_r = Ratio(residual.Norm(), denominator); + const double eta_norm = std::abs(modal_mass - 1.0); + const double eta_energy = + Ratio(std::abs(energy - lambda * modal_mass), denominator * norm_phi); + if (!std::isfinite(lambda) || !std::isfinite(eta_r) || eta_r > 1e-10 || + !std::isfinite(eta_norm) || eta_norm > 1e-10 || + !std::isfinite(eta_energy) || eta_energy > 1e-10 || + !std::isfinite(energy / modal_mass) || lambda / scale < -1e-9) + return Result::Failure( + Failure("solid-eigenpair-validation-failed", + "Original K/M pair checks failed.")); + state.eigenpair_checks.push_back({c + 1U, state.backend.column_ids[c], + lambda, modal_mass, energy / modal_mass, + eta_r, eta_norm, eta_energy}); + phi.push_back(std::move(vector)); + mass_phi.push_back(std::move(mv)); + } + for (std::size_t a = 0; a < count; ++a) + for (std::size_t b = 0; b < count; ++b) { + const double dot = phi[a].Dot(mass_phi[b]); + state.eta_gram = std::hypot(state.eta_gram, dot - (a == b ? 1.0 : 0.0)); + if (a != b) + state.max_eta_orth = + (std::max)(state.max_eta_orth, + std::abs(dot) / + std::sqrt(state.eigenpair_checks[a].modal_mass * + state.eigenpair_checks[b].modal_mass)); + } + if (!std::isfinite(state.eta_gram) || state.eta_gram > 1e-10 || + !std::isfinite(state.max_eta_orth) || state.max_eta_orth > 1e-10) + return Result::Failure( + Failure("solid-eigenpair-validation-failed", + "All-pair mass Gram checks failed.")); + state.vectors = Matrix(requested, dofs.FullDofCount()); + std::size_t cluster = 1U; + for (std::size_t c = 0; c < requested; ++c) { + const double lambda = state.backend.eigenvalues[c]; + if (c > 0 && (lambda - state.backend.eigenvalues[c - 1U]) / scale > 1e-8) + ++cluster; + std::optional frequency_value; + if (lambda >= 0.0) + frequency_value = std::sqrt(lambda) / (2.0 * std::acos(-1.0)); + state.modes.push_back( + {c + 1U, lambda, + static_cast(std::abs(lambda / scale) <= 1e-9 ? 0U : 1U), + cluster, frequency_value}); + // Strict > retains the first source/component entry for equal magnitudes. + std::size_t pivot = 0U; + for (std::size_t row = 1; row < phi[c].Size(); ++row) + if (std::abs(phi[c][row]) > std::abs(phi[c][pivot])) pivot = row; + const double sign = phi[c][pivot] < 0.0 ? -1.0 : 1.0; + for (std::size_t row = 0; row < phi[c].Size(); ++row) + state.vectors(c, dofs.FreeDofs()[row]) = sign * phi[c][row]; + } + state.backend.eigenvectors = Matrix(0U, 0U); + return Result::Success(std::move(state)); +} +} // namespace + +Result SolidFrequencyProcedureExecutor::Execute( + const Domain& domain) const { + try { + auto result = ExecuteCandidate(domain); + if (result.HasValue()) return result; + auto diagnostics = result.GetStatus().Diagnostics(); + for (auto& diagnostic : diagnostics) + if (diagnostic.location.file.empty()) + diagnostic.location = domain.Frequency() != nullptr + ? domain.Frequency()->frequency_location + : SourceLocation{domain.SourcePath(), 0U}; + return Result::Failure(Status::Failure( + result.GetStatus().Category().value_or(FailureCategory::kSolver), + std::move(diagnostics))); + } catch (const std::bad_alloc&) { + return Result::Failure(Failure( + "eigen-allocation-failed", "Solid modal candidate allocation failed.")); + } catch (const std::length_error&) { + return Result::Failure( + Failure("eigen-allocation-failed", + "Solid modal workspace exceeds addressable storage.")); + } +} +} // namespace fesa diff --git a/src/fesa/elements/hex8_element.cpp b/src/fesa/elements/hex8_element.cpp index a42c897..53ca80c 100644 --- a/src/fesa/elements/hex8_element.cpp +++ b/src/fesa/elements/hex8_element.cpp @@ -172,6 +172,13 @@ const ElementDofLayout& Hex8Element::DofLayout() const noexcept { Result Hex8Element::ComputeStiffness() const { return Result::Success({layout_, stiffness_}); } +std::array Hex8Element::ReferenceWeights() const noexcept { + std::array weights{}; + for (std::size_t g = 0; g < weights.size(); ++g) + weights[g] = points_[g].weight; + return weights; +} + Result Hex8Element::ComputeMass() const { if (!density_) return Result::Failure( diff --git a/src/fesa/solvers/eigen/mkl_feast_eigen_solver.cpp b/src/fesa/solvers/eigen/mkl_feast_eigen_solver.cpp index 883a327..f72bb2e 100644 --- a/src/fesa/solvers/eigen/mkl_feast_eigen_solver.cpp +++ b/src/fesa/solvers/eigen/mkl_feast_eigen_solver.cpp @@ -17,6 +17,8 @@ #include #include "fesa/math/sparse_matrix.h" +#include "fesa/math/vector.h" +#include "fesa/solvers/linear/mkl_pardiso_solver.h" #include "solvers/eigen/mkl_feast_eigen_solver_test_seam.h" namespace fesa { @@ -364,15 +366,57 @@ bool FillInitialBasis(std::vector& vectors, const std::size_t n, } // namespace +Result MklFeastEigenSolver::SolveGeneralized( + const SparseMatrix& stiffness, const SparseMatrix& mass, + const std::size_t requested_modes, const double mass_lower_bound, + const double classification_scale) const { + try { + return SolveImpl(stiffness, &mass, requested_modes, classification_scale, + mass_lower_bound); + } catch (const std::bad_alloc&) { + return Result::Failure( + Failure("eigen-allocation-failed", "generalized-workspace", + "Allocation of generalized sparse-pencil storage failed.")); + } +} + Result MklFeastEigenSolver::Solve( const SparseMatrix& whitened_operator, const std::size_t requested_modes, const double zero_operator_scale) const { + return SolveImpl(whitened_operator, nullptr, requested_modes, + zero_operator_scale, 1.0); +} + +Result MklFeastEigenSolver::SolveImpl( + const SparseMatrix& whitened_operator, const SparseMatrix* mass, + const std::size_t requested_modes, const double zero_operator_scale, + const double mass_lower_bound) const { + static_assert(sizeof(MKL_INT) == 8U, + "FEAST requires the installed ILP64 ABI."); const Status csr_status = whitened_operator.Validate(); if (!csr_status.IsOk()) return Result::Failure( Failure("feast-invalid-csr", "public-csr", "The public sparse matrix failed CSR validation.")); const std::size_t n = whitened_operator.Rows(); + if (mass != nullptr) { + if (!mass->Validate().IsOk() || mass->Rows() != n || mass->Columns() != n || + !std::isfinite(mass_lower_bound) || !(mass_lower_bound > 0.0) || + !std::isfinite(zero_operator_scale) || !(zero_operator_scale > 0.0)) { + return Result::Failure(Failure( + "feast-invalid-mass-bound", "generalized-pencil", + "A matching finite mass and certified positive bound are required.")); + } + const Status mass_symmetry = ValidateSymmetry(*mass); + if (!mass_symmetry.IsOk()) + return Result::Failure(mass_symmetry); + MklPardisoSolver spd_check; + const Status spd = spd_check.Factorize(*mass); + if (!spd.IsOk()) + return Result::Failure( + Failure("feast-mass-not-positive-definite", "mass-spd", + "Sparse SPD mass factorization failed.")); + } if (n == 0U || whitened_operator.Columns() != n || requested_modes == 0U || requested_modes > n) { return Result::Failure(Failure( @@ -399,7 +443,8 @@ Result MklFeastEigenSolver::Solve( for (std::size_t row = 0U; row < n; ++row) { double sum = 0.0; for (std::size_t p = offsets[row]; p < offsets[row + 1U]; ++p) { - sum += std::abs(values[p]); + sum = mass != nullptr ? NextUp(sum + std::abs(values[p])) + : sum + std::abs(values[p]); } if (!std::isfinite(sum)) { return Result::Failure(Failure( @@ -416,7 +461,7 @@ Result MklFeastEigenSolver::Solve( "An exactly zero operator requires a finite positive physical scale.")); } const double operator_scale = - zero_operator ? zero_operator_scale : infinity_norm; + mass != nullptr || zero_operator ? zero_operator_scale : infinity_norm; if (!std::isfinite(operator_scale) || !(operator_scale > 0.0)) { return Result::Failure( Failure("feast-nonfinite-operator-scale", "operator-scale", @@ -424,7 +469,8 @@ Result MklFeastEigenSolver::Solve( } auto a = CopyUpperOneBased(whitened_operator, operator_scale); if (!a.HasValue()) return Result::Failure(a.GetStatus()); - auto b = IdentityOneBased(n); + auto b = + mass != nullptr ? CopyUpperOneBased(*mass, 1.0) : IdentityOneBased(n); if (!b.HasValue()) return Result::Failure(b.GetStatus()); double lower = -1.0; @@ -432,7 +478,17 @@ Result MklFeastEigenSolver::Solve( double g_lower = 0.0; double g_upper = 0.0; double beta = 1.0; - if (!zero_operator) { + if (mass != nullptr) { + // ||K||2 <= ||K||inf and M >= m_lower*I bound the genuine pencil. + const double bound = + NextUp(NextUp(infinity_norm / mass_lower_bound) / operator_scale); + g_lower = -bound; + g_upper = bound; + beta = (std::max)(1.0, bound); + lower = std::nextafter(g_lower - 16.0 * Delta(g_lower, beta), + -(std::numeric_limits::infinity)()); + upper = NextUp(g_upper + 16.0 * Delta(g_upper, beta)); + } else if (!zero_operator) { g_lower = (std::numeric_limits::infinity)(); g_upper = -(std::numeric_limits::infinity)(); const auto& columns = whitened_operator.ColumnIndices(); @@ -467,6 +523,18 @@ Result MklFeastEigenSolver::Solve( LocalThreadGuard threads; std::vector calls; + if (mass != nullptr) + calls.push_back({0U, "mass-spd", "original-mass-factorization", 0.0, + mass_lower_bound, 0U, 0, false, 0U, false, 0, false, 0.0, + "accepted", + "sparse-SPD; caller-certified-roundoff-bound"}); + auto TraceFailure = [&](const std::string& code, const std::string& identity, + const std::string& message) { + return Failure(code, identity, + mass != nullptr + ? message + " Backend calls=" + CallHistory(calls) + : message); + }; std::array snapshot{}; SetControls(snapshot, false); std::vector base_fpm; @@ -482,16 +550,16 @@ Result MklFeastEigenSolver::Solve( count_res.assign(n, 0.0); count_x.assign(n * n, 0.0); } catch (const std::bad_alloc&) { - return Result::Failure( - Failure("eigen-allocation-failed", "count-workspace", - "Allocation of the required FEAST count workspace failed.")); + return Result::Failure(TraceFailure( + "eigen-allocation-failed", "count-workspace", + "Allocation of the required FEAST count workspace failed.")); } auto Count = [&](const double call_lower, const double call_upper, const std::string& purpose, std::size_t* count) -> Status { if (!internal::CanCallFeast(calls.size())) - return Failure("eigen-search-budget-exhausted", "count-calls", - "The bounded FEAST count-call budget was exhausted."); + return TraceFailure("eigen-search-budget-exhausted", "count-calls", + "The bounded FEAST count-call budget was exhausted."); std::array fpm{}; SetControls(fpm, true); MKL_INT m0 = static_cast(n); @@ -500,7 +568,7 @@ Result MklFeastEigenSolver::Solve( MKL_INT info = 0; double epsout = 0.0; if (!FillInitialBasis(count_x, n, n)) { - return Failure( + return TraceFailure( "eigenpair-validation-failed", "count-initial-basis", "The prescribed count initial basis failed QR validation."); } @@ -540,11 +608,12 @@ Result MklFeastEigenSolver::Solve( count_available ? "accepted" : "vendor-info=" + std::to_string(info); calls.push_back(std::move(event)); if (info != 0) - return Failure("eigen-feast-failed", "info=" + std::to_string(info), - "FEAST did not certify the requested count interval."); + return TraceFailure( + "eigen-feast-failed", "info=" + std::to_string(info), + "FEAST did not certify the requested count interval."); if (m0 < 0 || static_cast(m0) > n) { - return Failure("eigen-count-inconsistent", "count-output", - "FEAST count output is outside [0,n]."); + return TraceFailure("eigen-count-inconsistent", "count-output", + "FEAST count output is outside [0,n]."); } *count = static_cast(m0); for (std::size_t i = 0U; i + 1U < calls.size(); ++i) { @@ -555,7 +624,7 @@ Result MklFeastEigenSolver::Solve( if ((prior.upper == call_upper && prior.count != *count) || (prior.upper < call_upper && prior.count > *count) || (prior.upper > call_upper && prior.count < *count)) - return Failure( + return TraceFailure( "eigen-count-inconsistent", "count-history", "Successful prefix counts violate identity or monotonicity."); } @@ -599,10 +668,36 @@ Result MklFeastEigenSolver::Solve( return Result::Failure( !count_status.IsOk() ? count_status - : Failure( + : TraceFailure( "eigen-count-inconsistent", "total-enclosure", "The enclosing FEAST count did not equal matrix dimension.")); } + if (mass != nullptr) { + // Count the actual pencil before removing any part of the conservative + // negative enclosure. A far-negative endpoint degrades contour filtering + // of low positive modes on a stiff mesh; an empty guarded prefix proves + // that moving it cannot remove a mode, including a repeated rigid cluster. + const double margin = Delta(0.0, beta); + const double candidate = -16.0 * margin; + const std::array guards{candidate - margin, candidate, + candidate + margin}; + if (lower < guards.front()) { + for (const double guard : guards) { + std::size_t negative_count = 0U; + const Status prefix_status = + Count(lower, guard, "lower-empty-prefix", &negative_count); + if (!prefix_status.IsOk()) + return Result::Failure(prefix_status); + if (negative_count != 0U) + return Result::Failure(TraceFailure( + "negative-eigenvalue-spectrum", "certified-negative-prefix", + "The generalized pencil contains significant negative " + "eigenvalues.")); + } + lower = candidate; + calls.back().decision_reason = "guarded-empty-prefix; lowest-spectrum-preserved"; + } + } const double fixed_lower = lower; const double enclosure_upper = upper; double solve_upper = enclosure_upper; @@ -694,7 +789,7 @@ Result MklFeastEigenSolver::Solve( count_status = Count(fixed_lower, solve_upper, "final-count", &count); if (!count_status.IsOk() || count < requested_modes || count != certified_count) { - return Result::Failure(Failure( + return Result::Failure(TraceFailure( "eigen-count-inconsistent", "final-count", "The final FEAST interval count is not the certified guarded count.")); } @@ -709,27 +804,27 @@ Result MklFeastEigenSolver::Solve( if (!below_status.IsOk() || !above_status.IsOk() || below != count || above != count) { return Result::Failure( - Failure("eigen-count-inconsistent", "final-guard", - "The final interior endpoint is not guarded by equal " - "successful counts.")); + TraceFailure("eigen-count-inconsistent", "final-guard", + "The final interior endpoint is not guarded by equal " + "successful counts.")); } } if (count > (std::numeric_limits::max)() - 8U || count > ((std::numeric_limits::max)() - 1U) / 3U) { return Result::Failure( - Failure("eigen-allocation-failed", "solve-capacity", - "FEAST solve-capacity padding overflows size_t.")); + TraceFailure("eigen-allocation-failed", "solve-capacity", + "FEAST solve-capacity padding overflows size_t.")); } const std::size_t padded_count = (std::max)(count + 8U, (3U * count + 1U) / 2U); std::size_t capacity = count == n ? n : (std::min)(n, padded_count); if (!internal::ProductFits(n, capacity)) return Result::Failure( - Failure("feast-workspace-overflow", "solve-workspace", - "The FEAST solve workspace size overflows.")); + TraceFailure("feast-workspace-overflow", "solve-workspace", + "The FEAST solve workspace size overflows.")); if (!internal::FitsDoubleStorage(n * capacity) || !internal::FitsDoubleStorage(capacity)) { - return Result::Failure(Failure( + return Result::Failure(TraceFailure( "eigen-allocation-failed", "solve-workspace", "The required solve workspace exceeds addressable double storage.")); } @@ -742,23 +837,23 @@ Result MklFeastEigenSolver::Solve( std::size_t solve_capacity = capacity; for (std::size_t attempt = 0U; attempt < 3U; ++attempt) { if (!internal::CanCallFeast(calls.size())) - return Result::Failure( - Failure("eigen-search-budget-exhausted", "solve-call", - "The total backend call budget is exhausted before solve.")); + return Result::Failure(TraceFailure( + "eigen-search-budget-exhausted", "solve-call", + "The total backend call budget is exhausted before solve.")); try { internal::ApplyFeastWorkspaceTestHook(false); eigenvalues.assign(solve_capacity, 0.0); residuals.assign(solve_capacity, 0.0); vectors.assign(n * solve_capacity, 0.0); } catch (const std::bad_alloc&) { - return Result::Failure( - Failure("eigen-allocation-failed", "solve-workspace", - "Allocation of the required FEAST solve workspace failed.")); + return Result::Failure(TraceFailure( + "eigen-allocation-failed", "solve-workspace", + "Allocation of the required FEAST solve workspace failed.")); } if (!FillInitialBasis(vectors, n, solve_capacity)) { - return Result::Failure( - Failure("eigenpair-validation-failed", "solve-initial-basis", - "The prescribed solve initial basis failed QR validation.")); + return Result::Failure(TraceFailure( + "eigenpair-validation-failed", "solve-initial-basis", + "The prescribed solve initial basis failed QR validation.")); } std::array fpm{}; SetControls(fpm, false); @@ -794,21 +889,21 @@ Result MklFeastEigenSolver::Solve( : "vendor-info=" + std::to_string(info)}); if (info == 0 && (!std::isfinite(epsout) || loop < 0)) { return Result::Failure( - Failure("invalid-eigenvalue-result", "solve-output", - "A successful FEAST solve returned nonfinite epsout or " - "invalid loop.")); + TraceFailure("invalid-eigenvalue-result", "solve-output", + "A successful FEAST solve returned nonfinite epsout or " + "invalid loop.")); } if (info == 0 && returned >= 0 && static_cast(returned) <= solve_capacity) { for (std::size_t column = 0U; column < static_cast(returned); ++column) { if (!std::isfinite(eigenvalues[column])) - return Result::Failure(Failure( + return Result::Failure(TraceFailure( "invalid-eigenvalue-result", std::to_string(column), "A nonfinite eigenvalue cannot trigger an interval retry.")); for (std::size_t row = 0U; row < n; ++row) if (!std::isfinite(vectors[row + column * n])) - return Result::Failure(Failure( + return Result::Failure(TraceFailure( "invalid-eigenvalue-result", std::to_string(column), "A nonfinite eigenvector cannot trigger an interval retry.")); } @@ -818,10 +913,10 @@ Result MklFeastEigenSolver::Solve( if (endpoint_separated) { for (std::size_t index = 0U; index < count; ++index) { if (!std::isfinite(eigenvalues[index])) { - return Result::Failure( - Failure("invalid-eigenvalue-result", std::to_string(index), - "FEAST returned a nonfinite eigenvalue before endpoint " - "validation.")); + return Result::Failure(TraceFailure( + "invalid-eigenvalue-result", std::to_string(index), + "FEAST returned a nonfinite eigenvalue before endpoint " + "validation.")); } endpoint_separated = endpoint_separated && @@ -856,7 +951,7 @@ Result MklFeastEigenSolver::Solve( const Status recount_status = Count(fixed_lower, solve_upper, "final-count", &recounted); if (!recount_status.IsOk() || recounted < requested_modes) { - return Result::Failure(Failure( + return Result::Failure(TraceFailure( "eigen-count-inconsistent", "solve-retry", "A larger certified solve interval could not be recounted.")); } @@ -871,7 +966,7 @@ Result MklFeastEigenSolver::Solve( fixed_lower, solve_upper + retry_margin, "final-count", &above); if (!below_status.IsOk() || !above_status.IsOk() || below != count || above != count) { - return Result::Failure(Failure( + return Result::Failure(TraceFailure( "eigen-count-inconsistent", "solve-retry-guard", "The larger retry interval is not guarded by matching counts.")); } @@ -879,8 +974,8 @@ Result MklFeastEigenSolver::Solve( if (count > (std::numeric_limits::max)() - 8U || count > ((std::numeric_limits::max)() - 1U) / 3U) { return Result::Failure( - Failure("eigen-allocation-failed", "solve-capacity", - "FEAST retry capacity padding overflows size_t.")); + TraceFailure("eigen-allocation-failed", "solve-capacity", + "FEAST retry capacity padding overflows size_t.")); } const std::size_t retry_padding = (std::max)(count + 8U, (3U * count + 1U) / 2U); @@ -894,22 +989,22 @@ Result MklFeastEigenSolver::Solve( } if (!solve_accepted) { return Result::Failure( - Failure("eigen-feast-failed", "info=" + std::to_string(info), - "FEAST solve did not return the certified interval count " - "within its bounded retries.")); + TraceFailure("eigen-feast-failed", "info=" + std::to_string(info), + "FEAST solve did not return the certified interval count " + "within its bounded retries.")); } for (std::size_t index = 0U; index < count; ++index) { if (!std::isfinite(eigenvalues[index])) { - return Result::Failure( - Failure("invalid-eigenvalue-result", std::to_string(index), - "FEAST returned a nonfinite eigenvalue before sorting.")); + return Result::Failure(TraceFailure( + "invalid-eigenvalue-result", std::to_string(index), + "FEAST returned a nonfinite eigenvalue before sorting.")); } for (std::size_t row = 0U; row < n; ++row) { if (!std::isfinite(vectors[row + index * n])) { - return Result::Failure( - Failure("invalid-eigenvalue-result", std::to_string(index), - "FEAST returned a nonfinite eigenvector before sorting.")); + return Result::Failure(TraceFailure( + "invalid-eigenvalue-result", std::to_string(index), + "FEAST returned a nonfinite eigenvector before sorting.")); } } } @@ -936,28 +1031,34 @@ Result MklFeastEigenSolver::Solve( const double mu = eigenvalues[input]; if (!std::isfinite(mu)) { return Result::Failure( - Failure("invalid-eigenvalue-result", std::to_string(input), - "FEAST returned a nonfinite eigenvalue.")); + TraceFailure("invalid-eigenvalue-result", std::to_string(input), + "FEAST returned a nonfinite eigenvalue.")); } if (mu < -1.0e-9) { - return Result::Failure(Failure( + return Result::Failure(TraceFailure( "negative-eigenvalue-spectrum", std::to_string(input), "The raw normalized FEAST eigenvalue is significantly negative.")); } if (mu - fixed_lower < Delta(fixed_lower, beta) || solve_upper - mu < Delta(solve_upper, beta)) { - return Result::Failure(Failure( + return Result::Failure(TraceFailure( "eigen-endpoint-separation", std::to_string(input), "A returned FEAST eigenvalue is too close to a solve endpoint.")); } double norm = 0.0; for (std::size_t row = 0U; row < n; ++row) norm += vectors[row + input * n] * vectors[row + input * n]; + if (mass != nullptr) { + Vector raw(n); + for (std::size_t row = 0U; row < n; ++row) + raw[row] = vectors[row + input * n]; + norm = raw.Dot(mass->Multiply(raw)); + } norm = std::sqrt(norm); if (!(norm > 0.0) || !std::isfinite(norm)) return Result::Failure( - Failure("eigenpair-validation-failed", std::to_string(input), - "FEAST returned a non-normalizable eigenvector.")); + TraceFailure("eigenpair-validation-failed", std::to_string(input), + "FEAST returned a non-normalizable eigenvector.")); double residual_norm = 0.0; for (std::size_t row = 0U; row < n; ++row) { double ax = 0.0; @@ -965,12 +1066,43 @@ Result MklFeastEigenSolver::Solve( ax += values[p] / operator_scale * vectors[whitened_operator.ColumnIndices()[p] + input * n]; } - const double difference = ax - mu * vectors[row + input * n]; + double mx = vectors[row + input * n]; + if (mass != nullptr) { + mx = 0.0; + for (std::size_t p = mass->RowOffsets()[row]; + p < mass->RowOffsets()[row + 1U]; ++p) + mx += + mass->Values()[p] * vectors[mass->ColumnIndices()[p] + input * n]; + } + const double difference = ax - mu * mx; residual_norm += difference * difference; ordered_vectors(row, output) = vectors[row + input * n] / norm; } - const double rho_b = std::sqrt(residual_norm) / norm; + double residual_denominator = norm; + if (mass != nullptr) { + double norm_k = 0.0; + double norm_m = 0.0; + double euclidean = 0.0; + for (const double value : values) + norm_k = std::hypot(norm_k, value / operator_scale); + for (const double value : mass->Values()) + norm_m = std::hypot(norm_m, value); + for (std::size_t row = 0U; row < n; ++row) + euclidean = std::hypot(euclidean, vectors[row + input * n]); + residual_denominator = (norm_k + std::abs(mu) * norm_m) * euclidean; + } + const double rho_b = + residual_denominator == 0.0 + ? (residual_norm == 0.0 ? 0.0 + : (std::numeric_limits::infinity)()) + : std::sqrt(residual_norm) / residual_denominator; if (rho_b > 1.0e-10) { + if (mass != nullptr) + return Result::Failure(TraceFailure( + "eigenpair-validation-failed", std::to_string(input), + "Independent original-pencil residual exceeds 1e-10: eta=" + + FormatDouble(rho_b) + ", mu=" + FormatDouble(mu) + + "; calls=" + CallHistory(calls))); const FeastBackendCall& last_call = calls.back(); std::vector mirrored_action(n, 0.0); std::vector full_action(n, 0.0); @@ -999,7 +1131,7 @@ Result MklFeastEigenSolver::Solve( const double rho_b_upper = std::sqrt(mirrored_residual_squared) / norm; const double action_difference_norm = std::sqrt(action_difference_squared); - return Result::Failure(Failure( + return Result::Failure(TraceFailure( "eigenpair-validation-failed", std::to_string(input), "Independent whitened residual exceeds 1e-10: rho_b=" + FormatDouble(rho_b) + ", mu=" + FormatDouble(mu) + @@ -1019,16 +1151,16 @@ Result MklFeastEigenSolver::Solve( } const double lambda = mu * operator_scale; if (!std::isfinite(lambda) || (mu != 0.0 && lambda == 0.0)) { - return Result::Failure( - Failure("invalid-eigenvalue-result", std::to_string(input), - "Physical eigenvalue conversion overflowed or underflowed.")); + return Result::Failure(TraceFailure( + "invalid-eigenvalue-result", std::to_string(input), + "Physical eigenvalue conversion overflowed or underflowed.")); } ordered_values.push_back(lambda); const bool backend_res_finite = std::isfinite(residuals[input]); if (!backend_res_finite && mu != 0.0) { - return Result::Failure( - Failure("invalid-eigenvalue-result", std::to_string(input), - "A nonzero FEAST mode has a nonfinite backend residual.")); + return Result::Failure(TraceFailure( + "invalid-eigenvalue-result", std::to_string(input), + "A nonzero FEAST mode has a nonfinite backend residual.")); } ordered_backend_available.push_back(backend_res_finite); ordered_backend_res.push_back(backend_res_finite ? residuals[input] : 0.0); @@ -1042,23 +1174,35 @@ Result MklFeastEigenSolver::Solve( double dot = 0.0; for (std::size_t row = 0U; row < n; ++row) dot += ordered_vectors(row, left) * ordered_vectors(row, right); + if (mass != nullptr) { + dot = 0.0; + for (std::size_t row = 0U; row < n; ++row) { + double action = 0.0; + for (std::size_t p = mass->RowOffsets()[row]; + p < mass->RowOffsets()[row + 1U]; ++p) + action += mass->Values()[p] * + ordered_vectors(mass->ColumnIndices()[p], right); + dot += ordered_vectors(row, left) * action; + } + } const double expected = left == right ? 1.0 : 0.0; const double difference = std::abs(dot - expected); gram_squared += difference * difference; gram_error = (std::max)(gram_error, difference); if (std::abs(dot - expected) > 1.0e-10) - return Result::Failure(Failure( + return Result::Failure(TraceFailure( "eigenpair-validation-failed", std::to_string(left) + ":" + std::to_string(right), "The all-c FEAST vectors fail the independent Gram check.")); } } - gram_error = std::sqrt(gram_squared / static_cast(count)); + gram_error = std::sqrt(gram_squared / + (mass != nullptr ? 1.0 : static_cast(count))); if (gram_error > 1.0e-10) { - return Result::Failure( - Failure("eigenpair-validation-failed", "gram", - "The normalized all-c Frobenius Gram error exceeds 1e-10: " + - FormatDouble(gram_error) + ".")); + return Result::Failure(TraceFailure( + "eigenpair-validation-failed", "gram", + "The normalized all-c Frobenius Gram error exceeds 1e-10: " + + FormatDouble(gram_error) + ".")); } char version[256]{}; mkl_get_version_string(version, static_cast(sizeof(version))); diff --git a/tests/unit/analysis/solid_frequency_procedure_test.cpp b/tests/unit/analysis/solid_frequency_procedure_test.cpp new file mode 100644 index 0000000..a5a2bba --- /dev/null +++ b/tests/unit/analysis/solid_frequency_procedure_test.cpp @@ -0,0 +1,278 @@ +#include "fesa/analysis/solid_frequency_procedure.h" + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "fesa/io/abaqus/domain_mapper.h" +#include "fesa/io/abaqus/input_reader.h" +#include "fesa/model/domain.h" +#include "solvers/eigen/mkl_feast_eigen_solver_test_seam.h" + +namespace { +std::string StatusText(const fesa::Status& status) { + std::string text; + for (const auto& d : status.Diagnostics()) text += d.code + " " + d.message; + return text; +} +fesa::Result Brick(std::size_t requested = 6U, + const std::string& boundary = "", + double length = 1.0, double rho = 2.0, + double modulus = 120.0) { + constexpr int signs[8][3] = {{-1, -1, -1}, {1, -1, -1}, {1, 1, -1}, + {-1, 1, -1}, {-1, -1, 1}, {1, -1, 1}, + {1, 1, 1}, {-1, 1, 1}}; + std::ostringstream deck; + deck.precision(17); + deck << "*Part,name=P\n*Node\n"; + for (std::size_t a = 0; a < 8; ++a) + deck << (a + 1) * 10 << ',' << length * signs[a][0] << ',' + << 2 * length * signs[a][1] << ',' << 3 * length * signs[a][2] << '\n'; + deck << "*Element,type=C3D8,elset=E\n99,10,20,30,40,50,60,70,80\n" + "*Solid Section,elset=E,material=M\n,\n*End Part\n" + "*Assembly,name=A\n*Instance,name=I,part=P\n*End Instance\n" + "*End Assembly\n*Material,name=M\n*Elastic\n" + << modulus << ",.25\n"; + if (rho != 0.0) deck << "*Density\n" << rho << '\n'; + deck << "*Step,name=Modes\n*Frequency\n" << requested << '\n'; + if (!boundary.empty()) deck << "*Boundary\n" << boundary; + deck << "*End Step\n"; + static std::atomic serial{0}; + const auto path = + std::filesystem::temp_directory_path() / + ("fesa-solid-frequency-" + std::to_string(serial++) + ".inp"); + { + std::ofstream out(path); + out << deck.str(); + } + auto syntax = fesa::AbaqusInputReader{}.Read(path); + std::filesystem::remove(path); + if (!syntax.HasValue()) + return fesa::Result::Failure(syntax.GetStatus()); + return fesa::AbaqusDomainMapper{}.Map(syntax.Value()); +} +} // namespace +TEST(SolidFrequencyProcedure, FreeBrickRetainsSixRigidModes) { + auto domain = Brick(); + ASSERT_TRUE(domain.HasValue()) << StatusText(domain.GetStatus()); + const auto result = + fesa::SolidFrequencyProcedureExecutor{}.Execute(domain.Value()); + ASSERT_TRUE(result.HasValue()) << StatusText(result.GetStatus()); + ASSERT_EQ(result.Value().modes.size(), 6U); + for (const auto& mode : result.Value().modes) + EXPECT_EQ(mode.classification, 0U); + EXPECT_EQ(result.Value().vectors.Rows(), 6U); + EXPECT_EQ(result.Value().vectors.Columns(), 24U); +} + +namespace { +constexpr int kSigns[8][3] = {{-1, -1, -1}, {1, -1, -1}, {1, 1, -1}, + {-1, 1, -1}, {-1, -1, 1}, {1, -1, 1}, + {1, 1, 1}, {-1, 1, 1}}; +// Exact tensor-polynomial integrals, independent of production Gauss assembly. +double GradientIntegral(std::size_t a, std::size_t b, std::size_t i, + std::size_t j) { + double integral = 48.0 * kSigns[a][i] * kSigns[b][j] / + (64.0 * static_cast((i + 1U) * (j + 1U))); + for (std::size_t d = 0; d < 3; ++d) + if (d != i && d != j) integral *= 1.0 + kSigns[a][d] * kSigns[b][d] / 3.0; + return integral; +} +std::vector DenseOracle(const std::vector& free) { + const auto n = free.size(); + std::vector k(n * n, 0.0), m(n * n, 0.0), eigen(n, 0.0); + for (std::size_t row = 0; row < n; ++row) + for (std::size_t col = 0; col < n; ++col) { + const auto a = free[row] / 3U, b = free[col] / 3U; + const auto i = free[row] % 3U, j = free[col] % 3U; + double gradient_dot = 0.0; + for (std::size_t d = 0; d < 3; ++d) + gradient_dot += GradientIntegral(a, b, d, d); + k[row * n + col] = + 48.0 * (GradientIntegral(a, b, i, j) + GradientIntegral(a, b, j, i) + + (i == j ? gradient_dot : 0.0)); + if (i == j) { + m[row * n + col] = 2.0 * 48.0 / 64.0; + for (std::size_t d = 0; d < 3; ++d) + m[row * n + col] *= 1.0 + kSigns[a][d] * kSigns[b][d] / 3.0; + } + } + EXPECT_EQ( + LAPACKE_dsygv(LAPACK_ROW_MAJOR, 1, 'N', 'U', static_cast(n), + k.data(), static_cast(n), m.data(), + static_cast(n), eigen.data()), + 0); + return eigen; +} +const std::string kFaceFixed = "10,1,3\n40,1,3\n50,1,3\n80,1,3\n"; +std::vector FaceFree() { + std::vector free; + for (const auto node : {1U, 2U, 5U, 6U}) + for (std::size_t d = 0; d < 3; ++d) free.push_back(3U * node + d); + return free; +} +} // namespace +TEST(SolidFrequencyProcedure, + OriginalConsistentPencilMatchesIndependentDenseOracle) { + for (const bool constrained : {false, true}) { + const std::size_t count = constrained ? 12U : 24U; + auto domain = Brick(count, constrained ? kFaceFixed : ""); + ASSERT_TRUE(domain.HasValue()) << StatusText(domain.GetStatus()); + auto result = + fesa::SolidFrequencyProcedureExecutor{}.Execute(domain.Value()); + ASSERT_TRUE(result.HasValue()) << StatusText(result.GetStatus()); + std::vector free; + if (constrained) + free = FaceFree(); + else + for (std::size_t i = 0; i < 24; ++i) free.push_back(i); + const auto oracle = DenseOracle(free); + const auto& state = result.Value(); + ASSERT_EQ(state.modes.size(), oracle.size()); + ASSERT_EQ(state.eigenpair_checks.size(), oracle.size()); + for (std::size_t c = 0; c < oracle.size(); ++c) { + EXPECT_NEAR(state.modes[c].eigenvalue, oracle[c], + 1e-10 * state.classification_scale); + EXPECT_EQ(state.modes[c].mode_id, c + 1); + EXPECT_EQ(state.eigenpair_checks[c].sorted_rank, c + 1); + EXPECT_LE(state.eigenpair_checks[c].eta_r, 1e-10); + EXPECT_LE(state.eigenpair_checks[c].eta_energy, 1e-10); + EXPECT_LE(state.eigenpair_checks[c].eta_norm, 1e-10); + EXPECT_NEAR(state.eigenpair_checks[c].modal_mass, 1.0, 1e-10); + std::size_t pivot = 0; + for (std::size_t row = 0; row < 24; ++row) { + if (std::abs(state.vectors(c, row)) > std::abs(state.vectors(c, pivot))) + pivot = row; + if (state.constraint_mask[row]) EXPECT_EQ(state.vectors(c, row), 0.0); + } + EXPECT_GT(state.vectors(c, pivot), 0.0); + } + EXPECT_LE(state.eta_gram, 1e-10); + EXPECT_LE(state.max_eta_orth, 1e-10); + EXPECT_NEAR(state.mass_ideal_lower_bound, 12.0 / 27.0, 1e-13); + EXPECT_GT(state.mass_lower_bound, 0.0); + EXPECT_LT(state.mass_lower_bound, state.mass_ideal_lower_bound); + EXPECT_GT(state.mass_roundoff_guard, 0.0); + EXPECT_EQ(state.backend.integer_abi, "ILP64"); + EXPECT_EQ(state.backend.eigenvectors.Rows(), 0U); + EXPECT_EQ(state.backend.backend_calls.front().kind, "mass-spd"); + } +} +TEST(SolidFrequencyProcedure, DensityLengthAndMaterialScaling) { + const auto base_domain = Brick(4U, kFaceFixed); + ASSERT_TRUE(base_domain.HasValue()); + const auto base = + fesa::SolidFrequencyProcedureExecutor{}.Execute(base_domain.Value()); + ASSERT_TRUE(base.HasValue()) << StatusText(base.GetStatus()); + for (const auto parameters : {std::array{1, 8, 120, .25}, + std::array{2, 2, 120, .25}, + std::array{1, 2, 480, 4}}) { + const auto domain = + Brick(4U, kFaceFixed, parameters[0], parameters[1], parameters[2]); + ASSERT_TRUE(domain.HasValue()); + const auto result = + fesa::SolidFrequencyProcedureExecutor{}.Execute(domain.Value()); + ASSERT_TRUE(result.HasValue()) << StatusText(result.GetStatus()); + for (std::size_t c = 0; c < 4; ++c) + EXPECT_NEAR(result.Value().modes[c].eigenvalue, + base.Value().modes[c].eigenvalue * parameters[3], + 1e-10 * result.Value().classification_scale); + } +} +TEST(SolidFrequencyProcedure, OneFreeDofAndAllConstrainedAndExcessRequest) { + const std::string one = + "10,1,2\n20,1,3\n30,1,3\n40,1,3\n50,1,3\n60,1,3\n70,1,3\n80,1,3\n"; + auto domain = Brick(1U, one); + ASSERT_TRUE(domain.HasValue()); + auto result = fesa::SolidFrequencyProcedureExecutor{}.Execute(domain.Value()); + ASSERT_TRUE(result.HasValue()) << StatusText(result.GetStatus()); + EXPECT_NEAR(result.Value().modes[0].eigenvalue, DenseOracle({2U})[0], 1e-11); + const std::string all = + "10,1,3\n20,1,3\n30,1,3\n40,1,3\n50,1,3\n60,1,3\n70,1,3\n80,1,3\n"; + for (const auto pair : + {std::pair{1U, all}, {13U, kFaceFixed}}) { + auto invalid = Brick(pair.first, pair.second); + ASSERT_TRUE(invalid.HasValue()); + EXPECT_FALSE(fesa::SolidFrequencyProcedureExecutor{} + .Execute(invalid.Value()) + .HasValue()); + } +} +TEST(SolidFrequencyProcedure, + MissingDensityAndFailedBackendPreservePriorCandidate) { + auto missing = Brick(1U, kFaceFixed, 1.0, 0.0); + if (missing.HasValue()) + EXPECT_FALSE(fesa::SolidFrequencyProcedureExecutor{} + .Execute(missing.Value()) + .HasValue()); + auto domain = Brick(4U, kFaceFixed); + ASSERT_TRUE(domain.HasValue()); + const auto good = + fesa::SolidFrequencyProcedureExecutor{}.Execute(domain.Value()); + ASSERT_TRUE(good.HasValue()); + const auto previous = good.Value().modes[0].eigenvalue; + { + fesa::internal::ScopedFeastBackendTestHook hook{ + [](const auto& call, auto& output) { + if (!call.count) output.eigenvalues.back() = 9999.0; + }}; + EXPECT_FALSE(fesa::SolidFrequencyProcedureExecutor{} + .Execute(domain.Value()) + .HasValue()); + } + EXPECT_EQ(good.Value().modes[0].eigenvalue, previous); + const auto again = + fesa::SolidFrequencyProcedureExecutor{}.Execute(domain.Value()); + ASSERT_TRUE(again.HasValue()); + EXPECT_NEAR(again.Value().modes[0].eigenvalue, previous, 1e-12); +} + +TEST(SolidFrequencyProcedure, NearZeroNegativeRemainsRawWithoutFrequency) { + auto domain = Brick(); + ASSERT_TRUE(domain.HasValue()); + fesa::internal::ScopedFeastBackendTestHook hook{ + [](const auto& call, auto& output) { + if (!call.count) output.eigenvalues[0] = -1e-12; + }}; + const auto result = + fesa::SolidFrequencyProcedureExecutor{}.Execute(domain.Value()); + ASSERT_TRUE(result.HasValue()) << StatusText(result.GetStatus()); + EXPECT_LT(result.Value().modes[0].eigenvalue, 0.0); + EXPECT_EQ(result.Value().modes[0].classification, 0U); + EXPECT_FALSE(result.Value().modes[0].frequency.has_value()); + for (const auto& mode : result.Value().modes) EXPECT_EQ(mode.cluster_id, 1U); +} + +TEST(SolidFrequencyProcedure, + SuppliedFrequencyInputHasIndependentNumericalEvidence) { + auto root = std::filesystem::path{__FILE__}.parent_path(); + while (!root.empty() && !std::filesystem::exists(root / "CMakePresets.json")) + root = root.parent_path(); + ASSERT_FALSE(root.empty()); + const auto syntax = + fesa::AbaqusInputReader{}.Read(root / "reference/hexa/hexa_freq.inp"); + ASSERT_TRUE(syntax.HasValue()) << StatusText(syntax.GetStatus()); + const auto domain = fesa::AbaqusDomainMapper{}.Map(syntax.Value()); + ASSERT_TRUE(domain.HasValue()) << StatusText(domain.GetStatus()); + const auto result = + fesa::SolidFrequencyProcedureExecutor{}.Execute(domain.Value()); + ASSERT_TRUE(result.HasValue()) << StatusText(result.GetStatus()); + EXPECT_EQ(result.Value().modes.size(), 10U); + EXPECT_EQ(result.Value().vectors.Columns(), 3U * 442U); + EXPECT_GE(result.Value().eigenpair_checks.size(), 10U); + EXPECT_GT(result.Value().mass_lower_bound, 0.0); + EXPECT_LE(result.Value().eta_gram, 1e-10); + EXPECT_LE(result.Value().max_eta_orth, 1e-10); + // Original modal Value output is pending. This proves only our independent + // checks on the supplied input; it is not an external frequency comparison. +} diff --git a/tests/unit/solvers/eigen/mkl_feast_eigen_solver_test.cpp b/tests/unit/solvers/eigen/mkl_feast_eigen_solver_test.cpp index b3ef36d..52e3c31 100644 --- a/tests/unit/solvers/eigen/mkl_feast_eigen_solver_test.cpp +++ b/tests/unit/solvers/eigen/mkl_feast_eigen_solver_test.cpp @@ -629,3 +629,158 @@ TEST(MklFeastEigenSolver, UsesCertifiedEnclosureAfterOptionalInfoOne) { ASSERT_GE(result.Value().backend_calls.size(), 2U); EXPECT_EQ(result.Value().backend_calls[1U].outcome, "optional-abandoned"); } + +TEST(MklFeastGeneralized, GenuineNonDiagonalMassPencil) { + const auto k = MakeDenseCsr(2U, {6.0, 2.0, 2.0, 6.0}); + const auto m = MakeDenseCsr(2U, {2.0, 1.0, 1.0, 2.0}); + const auto result = + fesa::MklFeastEigenSolver{}.SolveGeneralized(k, m, 1U, 1.0, 3.0); + ASSERT_TRUE(result.HasValue()); + ASSERT_EQ(result.Value().eigenvalues.size(), 2U); + EXPECT_NEAR(result.Value().eigenvalues[0], 8.0 / 3.0, 1e-11); + EXPECT_NEAR(result.Value().eigenvalues[1], 4.0, 1e-11); +} + +TEST(MklFeastGeneralized, RejectsIndefiniteMassAndInvalidBound) { + const auto k = MakeDenseCsr(2U, {1, 0, 0, 2}); + const auto m = MakeDenseCsr(2U, {1, 2, 2, 1}); + ExpectFailureCode( + fesa::MklFeastEigenSolver{}.SolveGeneralized(k, m, 1U, 0.1, 2), + "feast-mass-not-positive-definite"); + ExpectFailureCode( + fesa::MklFeastEigenSolver{}.SolveGeneralized(k, k, 1U, 0, 1), + "feast-invalid-mass-bound"); +} + +TEST(MklFeastGeneralized, ValidatesPairsBeyondRequestedPrefix) { + const auto k = MakeDenseCsr(3U, {2, 0, 0, 0, 6, 0, 0, 0, 12}); + const auto m = MakeDenseCsr(3U, {2, 0, 0, 0, 3, 0, 0, 0, 4}); + fesa::internal::ScopedFeastBackendTestHook hook{ + [](const auto& call, auto& output) { + if (!call.count) output.vectors[2U + 2U * 3U] *= 2.0; + }}; + // Magnitude is repaired by mass normalization; an erroneous direction fails. + const auto result = + fesa::MklFeastEigenSolver{}.SolveGeneralized(k, m, 1U, 2, 3); + ASSERT_TRUE(result.HasValue()); + EXPECT_EQ(result.Value().eigenvalues.size(), 3U); +} + +TEST(MklFeastGeneralized, RejectsCorruptUnselectedDirectionAndFatalCount) { + const auto k = MakeDenseCsr(3U, {2, 0, 0, 0, 6, 0, 0, 0, 12}); + const auto m = MakeDenseCsr(3U, {2, 0, 0, 0, 3, 0, 0, 0, 4}); + { + fesa::internal::ScopedFeastBackendTestHook hook{ + [](const auto& call, auto& output) { + if (!call.count) output.vectors[0U + 2U * 3U] += 0.2; + }}; + ExpectFailureCode( + fesa::MklFeastEigenSolver{}.SolveGeneralized(k, m, 1U, 2, 3), + "eigenpair-validation-failed"); + } + { + fesa::internal::ScopedFeastBackendTestHook hook{ + [](const auto& call, auto& output) { + if (call.count) output.info = -3; + }}; + const auto result = + fesa::MklFeastEigenSolver{}.SolveGeneralized(k, m, 1U, 2, 3); + ExpectFailureCode(result, "eigen-feast-failed"); + EXPECT_NE( + result.GetStatus().Diagnostics()[0].message.find("Backend calls="), + std::string::npos); + } +} + +TEST(MklFeastGeneralized, LowestPrefixPreservesRepeatedBoundary) { + constexpr std::size_t n = 18U; + std::vector kv(n * n, 0.0), mv(n * n, 0.0); + for (std::size_t i = 0; i < n; ++i) { + mv[i * n + i] = 2.0 + static_cast(i); + kv[i * n + i] = + mv[i * n + i] * (i < 5U ? 1.0 : 2.0 + static_cast(i)); + } + const auto k = MakeDenseCsr(n, kv), m = MakeDenseCsr(n, mv); + const auto result = + fesa::MklFeastEigenSolver{}.SolveGeneralized(k, m, 3U, 2.0, 19.0); + ASSERT_TRUE(result.HasValue()) + << (result.HasValue() ? "" : result.GetStatus().Diagnostics()[0].message); + ASSERT_GE(result.Value().eigenvalues.size(), 5U); + for (std::size_t i = 0; i < 5U; ++i) + EXPECT_NEAR(result.Value().eigenvalues[i], 1.0, 1e-11); +} + +TEST(MklFeastGeneralized, AllocationAndCountMismatchFailClosed) { + const auto k = MakeDenseCsr(2U, {2, 0, 0, 6}); + const auto m = MakeDenseCsr(2U, {2, 0, 0, 3}); + { + fesa::internal::ScopedFeastWorkspaceTestHook hook{ + [](bool) { throw std::bad_alloc{}; }}; + ExpectFailureCode( + fesa::MklFeastEigenSolver{}.SolveGeneralized(k, m, 1U, 2, 2), + "eigen-allocation-failed"); + } + { + fesa::internal::ScopedFeastBackendTestHook hook{ + [](const auto& call, auto& output) { + if (call.count) output.m0 = 1; + }}; + ExpectFailureCode( + fesa::MklFeastEigenSolver{}.SolveGeneralized(k, m, 1U, 2, 2), + "eigen-count-inconsistent"); + } +} + +TEST(MklFeastGeneralized, GeneralizedRetryRetainsCertifiedCountTrace) { + constexpr std::size_t n = 18U; + std::vector kv(n * n, 0.0), mv(n * n, 0.0); + for (std::size_t i = 0; i < n; ++i) { + mv[i * n + i] = 2.0 + i; + kv[i * n + i] = mv[i * n + i] * (i + 1.0); + } + const auto k = MakeDenseCsr(n, kv), m = MakeDenseCsr(n, mv); + bool injected = false; + fesa::internal::ScopedFeastBackendTestHook hook{ + [&injected](const auto& call, auto& output) { + if (!call.count && !injected) { + output.info = 2; + injected = true; + } + }}; + const auto result = + fesa::MklFeastEigenSolver{}.SolveGeneralized(k, m, 3U, 2, 18); + ASSERT_TRUE(result.HasValue()) + << (result.HasValue() ? "" : result.GetStatus().Diagnostics()[0].message); + EXPECT_TRUE(injected); + EXPECT_NEAR(result.Value().eigenvalues[0], 1.0, 1e-11); + bool retried = false; + for (const auto& call : result.Value().backend_calls) + retried = retried || (call.kind == "solve" && call.info == 2 && + call.outcome == "retry"); + EXPECT_TRUE(retried); +} + +TEST(MklFeastGeneralized, + LowerPrefixCountsProveNoOmittedModesAndRejectNegative) { + const auto m = MakeDenseCsr(3U, {2, 0, 0, 0, 3, 0, 0, 0, 4}); + const auto k = MakeDenseCsr(3U, {0, 0, 0, 0, 6, 0, 0, 0, 12}); + const auto result = + fesa::MklFeastEigenSolver{}.SolveGeneralized(k, m, 1U, 2, 3); + ASSERT_TRUE(result.HasValue()); + std::size_t guards = 0; + for (const auto& call : result.Value().backend_calls) { + if (call.purpose != "lower-empty-prefix") continue; + ++guards; + EXPECT_EQ(call.info, 0); + EXPECT_TRUE(call.count_available); + EXPECT_EQ(call.count, 0U); + EXPECT_LT(call.upper, 0.0); + EXPECT_LT(call.lower, call.upper); + } + EXPECT_EQ(guards, 3U); + EXPECT_NEAR(result.Value().eigenvalues[0], 0.0, 1e-12); + const auto negative = MakeDenseCsr(3U, {-2, 0, 0, 0, 6, 0, 0, 0, 12}); + ExpectFailureCode( + fesa::MklFeastEigenSolver{}.SolveGeneralized(negative, m, 1U, 2, 3), + "negative-eigenvalue-spectrum"); +}