Solvers and preconditioners (`pyfe3d.solver`)#

Helpers for the eigenvalue problems that the element matrices of this package are usually assembled for, namely linear buckling

\[([K_{C_0}] + \lambda [K_G])\{u\} = \{0\}\]

and natural frequency

\[([K_{C_0}] - \omega_n^2 [M])\{u\} = \{0\}\]

The reason this module exists is that neither problem is safe to hand to scipy.sparse.linalg.eigsh() without preparation, and the two traps are easy to fall into and silent when they bite:

The matrices are badly scaled. A shell element mixes translations and rotations, so the diagonal of \([K_{C_0}]\) spans many orders of magnitude, and the drilling stiffness makes that worse: with the fictitious penalty of K6ROT the drilling entries are some six orders of magnitude below the membrane ones, and a condition number of \(10^{14}\) is easy to reach. diagonal_preconditioner() equilibrates the diagonal, which does not change the eigenvalues at all and usually decides whether the iterative solver converges.
The shift of the Cayley mode has to be chosen, not guessed. With eigsh(A=KG, M=K, sigma=sigma, mode='cayley', which='SM') the returned eigenvalues are the ones with the smallest \(|(\mu + \sigma)/(\mu - \sigma)|\), which are the most negative \(\mu\), and therefore the lowest positive load multipliers \(\lambda = -1/\mu\), only while \(\sigma\) is larger than \(|\mu|\) of those critical eigenvalues. Otherwise the eigenvalues closest to \(-\sigma\) come back instead, with no warning. estimate_cayley_sigma() estimates a safe shift from the matrices themselves.

On top of that, check_eigenpairs() verifies the result: it checks the residual of every eigenpair and, when \([K_{C_0}]\) is positive definite, that no load multiplier was missed below the lowest one found, by testing whether \([K_{C_0}] + s [K_G]\) is still positive definite just below it. That second check is the only cheap way to know that the lowest mode really is the lowest, and it is what distinguishes a converged answer from a plausible-looking one.

linear_buckling() and natural_frequency() put the three together and are what most users should call.

The mathematics of the shift estimate and of the verification follow the implementation of structsolve.linear_buckling.lb by the same author,

saullocastro/structsolve

which should be preferred when its additional solvers, such as the static condensation of the degrees-of-freedom where \([K_G]\) vanishes, are wanted. The versions here are self-contained so that this package keeps depending only on NumPy and SciPy.

Note

ARPACK, used by scipy.sparse.linalg.eigsh(), has been reported to return wrong eigenpairs at random, without raising, when linked against Intel MKL 2024.2.0 to 2025.0.0, which includes some Anaconda builds of SciPy. The dsteqr routine of those versions returns wrong eigenvectors for matrices larger than 32 by 32. MKL 2025.0.1 or newer is not affected. The verification of check_eigenpairs() catches such a failure, which is another reason to keep it on.

pyfe3d.solver.check_eigenpairs(K, KG, eigvals, eigvecs, rtol=0.001, min_rel_gap=0.0001)#

Verify the eigenpairs of \(([K_{C_0}] + \lambda [K_G])\{u\} = \{0\}\)

Two checks are made:

the relative residual \(||K u + \lambda K_G u||/(||K u|| + |\lambda| ||K_G u||)\) of every eigenpair with a finite \(\lambda\) must not exceed rtol;
no load multiplier is missing below the lowest positive \(\lambda_1\) returned: when \([K_{C_0}]\) is positive definite, \([K_{C_0}] + s [K_G]\) with \(s = \lambda_1 (1 - min\_rel\_gap)\) must also be positive definite.

The second check is the one that catches a solver that converged to an interior part of the spectrum, which is silent otherwise.

Parameters:

K, KGsparse_matrix: Constitutive and geometric stiffness matrices.
eigvalsarray_like: Load multipliers \(\lambda\).
eigvecsarray_like: Eigenvectors, one per column.
rtolfloat, optional: Largest acceptable relative residual.
min_rel_gapfloat, optional: Relative margin below \(\lambda_1\) at which the completeness of the spectrum is tested.

Returns:

errorstr or None: Description of the check that failed, or None when all passed.

pyfe3d.solver.diagonal_preconditioner(K, floor=1e-30)#

Jacobi preconditioner of a stiffness matrix

Returns the diagonal matrix \([D]\) with

\[D_{ii} = \frac{1}{\sqrt{max({K_{C_0}}_{ii}, floor)}}\]

such that \([D][K_{C_0}][D]\) has a unit diagonal. For a generalized eigenvalue problem both matrices are transformed, \([D][A][D]\) and \([D][B][D]\), which leaves the eigenvalues unchanged and maps the eigenvectors to \([D]^{-1}\{u\}\), so the eigenvectors of the original problem are recovered with u = D @ u_scaled.

Parameters:

Ksparse_matrix: Matrix whose diagonal sets the scaling, normally the constitutive stiffness matrix with the boundary conditions already applied.
floorfloat, optional: Lower bound for the diagonal entries, which guards against the null diagonal of an unconstrained degree-of-freedom.

Returns:

Dsparse_matrix: Diagonal preconditioner.

Examples

>>> D = diagonal_preconditioner(Kuu)
>>> eigvals, eigvecs = eigsh(A=D@KGuu@D, M=D@Kuu@D)
>>> eigvecs = D @ eigvecs

pyfe3d.solver.estimate_cayley_sigma(K, KG, safety=10.0, max_iter=50, rel_tol=0.001)#

Safe shift for the Cayley mode of scipy.sparse.linalg.eigsh()

The shift has to exceed \(|\mu|\) of the critical eigenvalues of \([K_G]\{u\} = \mu [K_{C_0}]\{u\}\), see the module documentation. The largest \(|\mu|\) is estimated with power iterations on \([K_{C_0}]^{-1}[K_G]\), whose ratio of norms in the \([K_{C_0}]\)-norm grows towards it, and the result is multiplied by safety.

The fallback value 1. is returned when \([K_{C_0}]\) is singular or not positive definite, or when the linear solution is inaccurate, since in those cases the power iteration says nothing.

Parameters:

K, KGsparse_matrix: Constitutive and geometric stiffness matrices, with the boundary conditions already applied.
safetyfloat, optional: Multiplier applied to the estimated largest \(|\mu|\).
max_iterint, optional: Maximum number of power iterations.
rel_tolfloat, optional: Relative increment below which the power iteration is stopped.

Returns:

sigmafloat: Shift to be passed to eigsh(..., sigma=sigma, mode='cayley').

pyfe3d.solver.is_positive_definite(A)#

Tells whether the symmetric sparse matrix A is positive definite

A is factorized with a symmetric permutation and without row interchanges. The Gaussian elimination of a positive definite matrix never meets a zero or a negative pivot, and by Sylvester’s law of inertia a negative pivot means a negative eigenvalue.

pyfe3d.solver.linear_buckling(K, KG, num_eigvalues=6, tol=0, sigma=None, precondition=True, check=True, check_rtol=0.001)#

Linear buckling analysis with a verified spectrum

Solves \(([K_{C_0}] + \lambda [K_G])\{u\} = \{0\}\) with scipy.sparse.linalg.eigsh() in Cayley mode, using the shift of estimate_cayley_sigma(), the preconditioner of diagonal_preconditioner() and the verification of check_eigenpairs().

Parameters:

K, KGsparse_matrix: Constitutive and geometric stiffness matrices, with the boundary conditions already applied, i.e. only the free degrees-of-freedom.
num_eigvaluesint, optional: Number of load multipliers to compute.
tolfloat, optional: Tolerance passed to scipy.sparse.linalg.eigsh(), 0 meaning machine precision.
sigmafloat or None, optional: Shift of the Cayley mode. With None it is estimated from the matrices, which is the recommended use.
preconditionbool, optional: Whether to equilibrate the diagonal before solving. This does not change the eigenvalues.
checkbool, optional: Whether to verify the eigenpairs and raise on failure.
check_rtolfloat, optional: Largest acceptable relative residual of the eigenpairs.

Returns:

eigvalsnp.ndarray: Load multipliers \(\lambda\), the positive ones first in increasing order, followed by the negative ones.
eigvecsnp.ndarray: Eigenvectors, one per column, in the same order.

Raises:

RuntimeError: When check=True and the verification fails, which usually means that the mesh, and not the solver, needs attention: a spectrum with many nearly coincident eigenvalues, as in the buckling of a cylindrical shell, is the typical case.

pyfe3d.solver.natural_frequency(K, M, num_eigvalues=6, tol=0, precondition=True)#

Natural frequencies of \(([K_{C_0}] - \omega_n^2 [M])\{u\} = \{0\}\)

Parameters:

K, Msparse_matrix: Constitutive stiffness and mass matrices, with the boundary conditions already applied.
num_eigvaluesint, optional: Number of frequencies to compute.
tolfloat, optional: Tolerance passed to scipy.sparse.linalg.eigsh().
preconditionbool, optional: Whether to equilibrate the diagonal of \([K_{C_0}]\) before solving.

Returns:

omegannp.ndarray: Circular frequencies in rad/s, in increasing order.
eigvecsnp.ndarray: Eigenvectors, one per column.

Notes

A mass matrix without rotary inertia, and one without inertia associated with the drilling rotation, is singular, so the eigenvalues associated with those degrees-of-freedom are infinite. The shift-invert mode used here looks for the eigenvalues closest to zero and is not troubled by them, but a lumped mass matrix will place the massless degrees-of-freedom wherever the factorization puts them, so prefer a consistent mass matrix when the lowest frequencies matter.

Solvers and preconditioners (pyfe3d.solver)#

Solvers and preconditioners (`pyfe3d.solver`)#