Block subspace projection preconditioned conjugate gradient method in modal structural analysis
Introduction
Determination of natural frequencies and natural oscillation modes is an important problem that arises in the calculation of buildings and structures subjected to the wind and seismic actions, and other types of dynamic loads. Currently, the finite element method is the main analysis tool; its application reduces the problem of determining the natural frequencies and natural oscillation modes to an algebraic partial generalized eigenvalue problem where is the symmetric positive-definite sparse stiffness matrix, is the symmetric semidefinite sparse or diagonal mass matrix, is the eigenvector and is the eigenvalue, , N is the dimension of the problem, n is the number of required eigenpairs, where N n,
In the case of large-scale design models, the lower part of the spectrum usually includes a significant number of local oscillation modes which do not contribute much to the seismic response of the structure [1], [2], [3]. The number of required eigenpairs for seismic analysis is not known in advance, and the eigenvalue solver extracts eigenpairs until a sufficient amount of modal masses is provided for each of the seismic input directions [1], [2]. Therefore, in the case of seismic analysis, the number of required eigenpairs depends on a large extent on the specifics of the problem and can reach 200–5000 or more. The contribution of ith eigenmode to seismic response is determined by so-called mass participation factor: where , – seismic input direction, – global coordinate system, – vector of the seismic input direction has 1.0 at the degrees of freedom (DOF), corresponding to the seismic input direction, and 0 at the all remaining DOFs, is a mass matrix and is an eigenvector. The modal mass for i mode and input seismic direction dir is defined as: where is a total mass in direction dir. If we take into account all eigenmodes of the discrete model, then But usually, we want to take into account a relatively small number of eigenmodes , and
If we take the less than eigenmodes, we make an error, and internal forces are less than their “exact” values, corresponding to 100% sum of modal masses. That is, we underestimate the magnitudes of internal efforts.
For this reason, seismic design codes require that the analysis takes into account such a number of eigenmodes that provide a sum of modal masses of at least 90% for each of the horizontal seismic directions and 70%–80% for the vertical one. For large-scale design models, it often happens that in the lower part of the spectrum there are a large number of local natural vibration modes that have a small percentage of modal masses (4). Such modes do not make a significant contribution to the seismic response of the structure, which leads to the extraction of a large number of eigenpairs until a sufficient amount of modal masses is achieved. Moreover, often local vibration modes have close frequencies. For instance, a multi-story building has a large number of identical floors, the local vibration modes for which are very close. The following example (Fig. 1) illustrates this.
Fig. 2 shows an increase in the sum of modal masses for the seismic input directions , respectively, with an increase in the number of eigenpairs taken into account in the analysis. Horizontal dot lines correspond to the provisions of most seismic norms requiring at least 90% of the sums of modal masses for the horizontal directions of the seismic input and at least 70% of the sum of modal masses for the vertical direction .
Fig. 3 represents a natural vibration mode to the first frequency, having a local character and giving a small contribution to the movement of the structure in the horizontal direction. The local nature of this eigenmode is because only the floor slab panels oscillate. And only starting from the 523 eigenmode, vibration modes appear that give a significant contribution to the movement of the structure in the horizontal direction (Fig. 4).
Such calculations are usually performed in small and medium design offices using multi-core desktops. Although the majority of the algorithms proposed in this article can be transferred to machines with distributed memory without fundamental modifications, the main development and testing were carried out on multi-core computers with shared memory, since we focus on the capabilities of the average user. For the development of such approaches for multi-core computers of the SMP (Symmetric Multiprocessing) architecture to be effective, it is necessary to take into account their specifics: the limited amount and limited bandwidth of RAM.
Modern software usually use either the block Lanczos method [3], [4], [5] or the block subspace iteration method [2], [6] for solving the problem (1).
Besides, for large problems, the Arnoldi method [7] is used to extract the eigenpairs lying in the lower part of the spectrum. The second-order Arnoldi method [8] is applied to solving the quadratic eigenvalue problem , where is a dissipation matrix. The problem of steady-state damped vibrations of a structure subjected to the action of harmonic loading is reduced to this type. We also note the work [9] in which the linearization-based method converts the rational eigenvalue problem into a linear eigenvalue problem. It is shown that the proposed approach can be used not only to solve the problem of steady-state oscillations of a damped system but also to analyze the natural oscillations of unconstrained systems. In all approaches noted the inverse iteration method [1] is the leading procedure: where is a block of iterated vectors, k is an iteration number. Each step of the procedure (7) approximates the iterated vectors to the eigenvectors of the problem (1) corresponding to the lower part of the spectrum. The system of linear algebraic equations (7) is usually solved using one of the direct methods for sparse matrices, and the matrix is preliminarily factorized according to Cholesky where is a sparse lower triangular matrix. It follows from the substitution of (8) in (7) that the inverse iteration procedure comes down to performing the forward and back substitutions. If the dimension of the problem N is several million equations, then the size of the lower triangular matrix – complete Cholesky factor of the stiffness matrix – is in the range of 10–40 GB for the considered class of problems. In the vast majority of cases, this problem is solved on the widely available multi-core computers with a limited amount of RAM. Therefore, the matrix is written to disk block by block. Thus, the above-mentioned amount of data has to be double read from the disk at each step of the inverse iteration method, which leads to a significant degradation in the performance of the computational process running at the disk access speed, and not at the fast processor speed.
One of the alternative approaches to solving the problem (1), which does not require factorization of a large-scale sparse matrix , is the refined shifted block inverse-free Krylov subspace method [10], [11], [12].
Another alternative approach to solving the problem (1) avoiding a complete factoring of large sparse matrix is the conjugate gradient method, which is based on finding the stationary points of the Rayleigh quotient when using the line search procedure where is the approximation of the ith eigenvector at the iteration step k, is the conjugate direction vector and is the optimization parameter for Rayleigh quotient minimization. In the classical formulation [13], [14] the conjugate gradient method starts with the lowest natural mode () and iterates until the norm of the residual vector , divided by the norm of the vector , becomes less than a certain small number tol: where
Then the second eigenpair (i 2) is determined, and the vectors and are orthogonalized to the convergent eigenvector at each step of the iterations. And so on until the required number of eigenpairs is extracted.
The conjugate gradient method does not require the factorization of the large matrix . Therefore, most problems of structural mechanics can be solved on widely used multi-core shared memory computers without using a disk. Design models of buildings and structures consist of many thin-walled structural elements — plates, rods and shells, for which a significant difference in tensile–compressive stiffness from flexural and shear ones is typical. These and other features of design models of buildings and structures [15], [16] lead to a deterioration in the conditioning number of the problem and a significant slowdown or even locking in the convergence of the above version of conjugate gradient method.
The preconditioning, the shift technique with spectral transformations and simultaneous iteration of several vectors in the block are applied to provide the stability of the computational process and improve the performance of the calculations.
The preconditioning consists of the transition from the original problem (1) to the following one: where is the preconditioning operator which changes the surface given by the Rayleigh quotient (9) to reduce the elongation of the isolines in the vicinity of the stationary points.
Different types of preconditioning have been proposed in [7], [13], [14], [17], [18], [19], [20], [21], [22], [23], [24] and others. When solving the eigenvalue problem (1), the development of efficient preconditioning is necessary, but not sufficient, since when solving the practical structural mechanics problems the convergence is still often locked.
Another technique stabilizing the convergence is to use a shift with spectral transformations. We take preconditioning as where is an incomplete factor of stiffness matrix , or is determined directly as an incomplete factorization , , – is a shift. The value of should be as close as possible to the eigenvalue whose eigenpair is sought for [7], [20]. Moreover, [18] and [20] show the relationship between the PCG method with a shift in preconditioning and the shift in an inverse iteration method.
To better separate close eigenvalues and accelerate the convergence of the conjugate gradient method, spectral transformations are used in the form
The parameters and used here adjust the spectral transformations for a given frequency interval so that the first Ritz pairs, where is a given number, are more isolated than without using spectral transformations. This distinguishes proposed approach from the previously used ones, for example, [4], [5], [25]. The selection of the parameters and will be considered in Section 2.2.4.
The residual vectors of the problem (15) are determined as where is the approximation of the eigenvector at the current step, and is the approximation of the respective eigenvalue. It follows from (16) that: where . Relations (18), (19) connect the residual vector and vector of the problem with spectral transformations (15) with the corresponding vectors of the original preconditioned problem (13).
Another important factor providing the computational stability and high performance of the PCG method is the simultaneous iteration of a block of vectors. One of such approaches, the Local Block PCG (LOBPCG) method, is considered in [26], [27].
Approximations of eigenvectors and conjugate direction vectors are taken as: where is the dimension of the block of iterated vectors, is an iteration number and vectors are defined from solution of the equation set , where is taken from (12).
In the above papers – the dimension of block is taken as the number of required eigenpairs. When this approach is used on multi-core shared memory computers, the number of computed eigenpairs is significantly limited due to the limited amount of RAM. Moreover, when solving structural mechanics problems, we have encountered low computational stability of the above approach: as soon as any Ritz pair begins to converge, the corresponding vectors and tend to zero vectors, and the approximation (20) becomes singular.
The block subspace projection PCG (BSPPCG) method using the expressions (20), is proposed in the previous papers of the authors [28], [29], but unlike the LOBPCG method, the columns of the projection matrix [see (23)] are orthogonalized before they become linearly dependent. Thus, the linear independence of the basis defining the subspace , is provided. Moreover, the dimension of the block m is assumed to be constant and independent of the number of sought eigenpairs n. As soon as the converged eigenpairs are detected in the block, they are removed from the block and placed in the result storage, and new starting vectors are generated in their place in the block, which are then orthogonalized against all the vectors in the block and against all the converged eigenvectors. Since the proposed approach has been developed for multi-core computers with a limited amount of RAM, it is extremely important to retain the relatively small dimension of the block , and usually is much less than .
This paper is a generalization of the approach, described in [28], [29], and focuses on the use of the shift technique with spectral transformations, adjusted to rarefaction a given number of frequencies of current frequency interval, as well as on the error accumulation arising during the orthogonalization of vectors in a block against converged Ritz vectors. In order to effectively use the processor cache and hide the low memory bandwidth, which is typical for most widely used shared memory computers, we tried to coalesce groups of vectors into dense blocks and, wherever possible, pass from BLAS level 1 and 2 procedures to the matrix multiplication procedure dgemm of level 3 BLAS.
Section snippets
Subspace projection approach
We take an approximation similar to (20): where the vectors are determined according to (19). Applying the projection method to the expression in parentheses of Eq. (15) we obtain where:
Convergence in the case of orthogonalization to converged modes
When orthogonalization to converged modes is used, we often encounter the effect depicted in Fig. 6, which shows the process of convergence of the first four Ritz pairs, and the error is determined according to (11). The first and second eigenpairs had converged at the 7th step (, where ), and then they were removed from the block and replaced by new starting vectors. This moment is shown in Fig. 6 as a jump of the error at the 8th iteration step for the curves 1 and 2.
Preconditioning
The incomplete Cholesky factorization “by value” approach based on the sparse matrices technique is used in this paper. The method is based on the fact that all off-diagonal elements satisfying the following inequality are removed from the lower triangular matrix during the incomplete factorization of the stiffness matrix , where is the rejection parameter, . The closer is the parameter to zero, the closer is the lower triangular matrix of incomplete factorization to
Numerical results
The design models are taken from the collection of SCAD Soft, IT company, developer of one of the most popular finite element modeling software, which is used for structural analysis and design, and is certified according to building regulations of the CIS countries [34]. The research was performed on computer A with a 16-core processor AMD Opteron 6276, 2.3/3.2 GHz, 64 GB DDR3 RAM, OS Windows Server 2008 R2 Enterprise SP1, 64 bit, and on computer B with 4-core processor Intel®Core™ i7 – 2760QM
Conclusion
The proposed BSPPCG method with incomplete block factorization preconditioning and shift technique with spectral transformations demonstrates high computational stability and robustness. As example 1 shows, the key point here is the orthogonalization of basis vectors in the iterable block when the loss of linear independence between them occurs. The next important point in this implementation of the method is the constant block size of iterable vectors, which makes it possible to determine
CRediT authorship contribution statement
Sergiy Fialko: Conceptualization, Methodology, Software, Validation, Formal analysis, Writing - original draft, Writing - review & editing, Visualization. Viktor Karpilovskyi: Software, Resources, Supervision, Project administration, Funding acquisition.
Acknowledgment
The authors are deeply grateful to the SCAD Soft IT company, Ukraine for the financial support of this study and the large collection of real-life problems offered by them.
References (37)
- et al.
An eigensolution strategy for large systems
Comput. Struct.
(1983) - et al.
A refined shifted block inverse-free Krylov subspace method for symmetric generalized eigenvalue problems
Comput. Math. Appl.
(2013) - et al.
An improved iterative optimization technique for the leftmost eigenpairs of large symmetric matrices
J. Comput. Phys.
(1988) Iterative methods for solving large-scale problems of structural mechanics using multi-core computers
Arch. Civ. Mech. Eng.
(2014)Preconditioning eigenvalues and some comparison of solvers
J. Comput. Appl. Math.
(2000)Parallel direct solver for solving systems of linear equations resulting from finite element method on multi-core desktops and workstations
Comput. Math. Appl.
(2015)PARFES: A method for solving finite element linear equations on multi-core computers
Adv. Eng. Softw.
(2010)- et al.
Dynamics of Structures
(2003) Three Dimensional Static and Dynamic Analysis of Structures
(2000)- S.Yu. Fialko, E.Z. Kriksunov, V.S. Karpilovskyy, A block Lanczos method with spectral transformations for natural...
The spectral transformation Lanczos method for the numerical solution of large sparse generalize symmetric eigenvalue problem
Math. Comp.
A shifted block Lanczos algorithm for solving sparse symmetric generalized eigenproblems
SIAM J. Matrix Anal. Appl.
Numerical methods for large eigenvalue problems
Solving rational eigenvalue problems via linearization
SIAM J. Matrix Anal. Appl.
SOAR: a second-order Arnoldi method for the solution of the quadratic eigenvalue problem
SIAM J. Matrix Anal. Appl.
A refined Arnoldi type method for large scale eigenvalue problems
Japan J. Ind. Appl. Math.
A refined variant of the inverse-free Krylov subspace method for symmetric generalized eigenvalue problems
Japan J. Ind. Appl. Math.
Solving Large - Scale Problems in Mechanics
Cited by (4)
Application of Acoustic-solid Coupling Theory in New Energy Vehicle Noise Control
2023, WSEAS Transactions on Applied and Theoretical MechanicsBlock Subspace Iteration Method for Structural Analysis on Multicore Computers
2022, Proceedings of the 17th Conference on Computer Science and Intelligence Systems, FedCSIS 2022