A new parallel sparse direct solver: Presentation and numerical experiments in large‐scale structural mechanics parallel computing

The main purpose of this work is to present a new parallel direct solver: Dissection solver. It is based on LU factorization of the sparse matrix of the linear system and allows to detect automatically and handle properly the zero‐energy modes, which are important when dealing with DDM. A performance evaluation and comparisons with other direct solvers (MUMPS, DSCPACK) are also given for both sequential and parallel computations. Results of numerical experiments with a two‐level parallelization of large‐scale structural analysis problems are also presented: FETI is used for the global problem parallelization and Dissection for the local multithreading. In this framework, the largest problem we have solved is of an elastic solid composed of 400 subdomains running on 400 computation nodes (3200 cores) and containing about 165 millions dof. The computation of one single iteration consumes less than 20 min of CPU time. Several comparisons to MUMPS are given for the numerical computation of large‐scale linear systems on a massively parallel cluster: performances and weaknesses of this new solver are highlighted. Copyright © 2011 John Wiley & Sons, Ltd.     

SOLUTION OF LARGE LINEAR SYSTEMS ON PIPELINED SIMD MACHINES

We developed a direct out-of-core solver for dense non-symmetric linear systems of ‘arbitrary’ size N×N. The algorithm fully employs the Basic Linear Algebra Subprograms (BLAS), and can therefore easily be adapted to different computer architectures by using the corresponding optimized routines. We used blocked versions of left-looking and right-looking variants of LU decomposition to perform most of the operations in Level 3 BLAS, to reduce the number of I/O operations and to minimize the CPU time usage. The storage requirements of the algorithm are only 2N×NB data elements where NB≪N. Depending on the sustained floating point performance and the sustained I/O rate of the given hardware, we derived formulas that allow for choosing optimal values of NB to balance between CPU time and I/O time. We tested the algorithm by means of linear systems derived from 3D-BEM for strongly and weakly singular integral equations and from interpolation problems for scattered data on closed surfaces in ℝ³. It took only about 2⋅5 CPU minutes on a 5 GFLOPS vector computer SNI S600/20 to solve a linear system of size 10000, which corresponds to a performance of 4⋅3 GFLOPS; a value of NB=650 gives a reasonable I/O time and the necessary main storage size is about 13 Mwords. In addition, we compared the algorithm with (1) an out-of-core version of GMRES and (2) a wavelet transform followed by in-core GMRES after thresholding. At least for boundary integral equations of classical boundary value problems of potential theory, the out-of-core version of GMRES is superior to the direct out-of-core solver and the wavelet transform since the algorithm converged after at most 5 iteration steps. It took about 17 s to solve a system with 8192 unknowns compared with 146 s for direct out-of-core and 402 s for wavelet transform followed by in-core GMRES. © 1997 by John Wiley & Sons, Ltd.     

CMRH method as iterative solver for boundary element acoustic systems

Disceretization of boundary integral equations leads to complex and fully populated linear systems. One inherent drawback of the boundary element method (BEM) is that, the dense linear system has to be constructed and solved for each frequency. For large-scale problems, BEM can be more efficient by improving the construction and solution phases of the linear system. For these problems, the application of common direct solver is inefficient. In this paper, the corresponding linear systems are solved more efficiently than common direct solvers by using the iterative technique called CMRH (Changing Minimal Residual method based on Hessenberg process). In this method, the generation of the basis vectors of the Krylov subspace is based on the Hessenberg process instead of Arnoldi's one that the most known GMRES (Generalized Minimal RESidual) solver uses. Compared to GMRES, the storage requirements are considerably reduced in CMRH.     

Non-intrusive parallelization of multibody system dynamic simulations

This paper evaluates two non-intrusive parallelization techniques for multibody system dynamics: parallel sparse linear equation solvers and OpenMP. Both techniques can be applied to existing simulation software with minimal changes in the code structure; this is a major advantage over Message Passing Interface, the standard parallelization method in multibody dynamics. Both techniques have been applied to parallelize a starting sequential implementation of a global index-3 augmented Lagrangian formulation combined with the trapezoidal rule as numerical integrator, in order to solve the forward dynamics of a variable-loop four-bar mechanism. Numerical experiments have been performed to measure the efficiency as a function of problem size and matrix filling. Results show that the best parallel solver (Pardiso) performs better than the best sequential solver (CHOLMOD) for multibody problems of large and medium sizes leading to matrix fillings above 10. OpenMP also proved to be advantageous even for problems of small sizes. Both techniques delivered speedups above 70% of the maximum theoretical values for a wide range of multibody problems.     

Numerical finite element formulation of the Schapery non‐linear viscoelastic material model

This study presents a numerical integration method for the non‐linear viscoelastic behaviour of isotropic materials and structures. The Schapery's three‐dimensional (3D) non‐linear viscoelastic material model is integrated within a displacement‐based finite element (FE) environment. The deviatoric and volumetric responses are decoupled and the strain vector is decomposed into instantaneous and hereditary parts. The hereditary strains are updated at the end of each time increment using a recursive formulation. The constitutive equations are expressed in an incremental form for each time step, assuming a constant incremental strain rate. A new iterative procedure with predictor–corrector type steps is combined with the recursive integration method. A general polynomial form for the parameters of the non‐linear Schapery model is proposed. The consistent algorithmic tangent stiffness matrix is realized and used to enhance convergence and help achieve a correct convergent state. Verifications of the proposed numerical formulation are performed and compared with a previous work using experimental data for a glassy amorphous polymer PMMA. Copyright © 2003 John Wiley & Sons, Ltd.     

High-performance multilevel iterative aggregation solver for large finite-element structural analysis problems

The present study is aimed to overcome difficulties faced with industrial applications of multilevel iterative methods to arbitrary finite element (FE) structural analysis problems. The coarse grid concept, used in multigrid methods, is substituted with an aggregation coarse model based on the mechanical principle. On the base of this approach together with previously developed multilevel preconditioner, an efficient iterative equation solver FEAGS was developed for using in standard comprehensive finite element software systems. Numerical examples for analyses of three-dimensional (3-D) frame structures demonstrate the efficiency of the method. Comparison with incomplete Cholesky conjugate gradient (ICCG) method and direct methods is presented.     

A Generalised Conjugate Residual method for the solution of non‐symmetric systems of equations with multiple right‐hand sides

This paper presents an iterative algorithm for solving non‐symmetric systems of equations with multiple right‐hand sides. The algorithm is an extension of the Generalised Conjugate Residual method (GCR) and combines the advantages of a direct solver with those of an iterative solver: it does not have to restart from scratch for every right‐hand side, it tends to require less memory than a direct solver, and it can be implemented efficiently on a parallel computer. We will show that the extended GCR algorithm can be competitive with a direct solver when running on a single processor. We will also show that the algorithm performs well on a Cray T3E parallel computer. Copyright © 1999 John Wiley & Sons, Ltd.     

The bipenalty method for arbitrary multipoint constraints

In finite element (FE) analysis, traditional penalty methods impose constraints by adding virtual stiffness to the FE system. In dynamics, this can decrease the critical time step of the system when conditionally stable time integration schemes are used by introducing spurious modes with high eigenfrequencies. Recent studies have shown that using mass penalties alongside traditional stiffness penalties can mitigate this effect for systems with a one single‐point constraint. In the present work, we extend this finding to include systems with an arbitrary set of multipoint constraints. By analysing the generalised eigenvalue problem, we show that the values of spurious eigenfrequencies may be controlled by the choice of stiffness and mass penalty parameters. The method is demonstrated using numerical examples, including a one‐dimensional contact–impact formulation and a two‐dimensional crack propagation analysis. The results show that constraint imposition using the bipenalty method can be employed such that the critical time step of an analysis is unaffected, whereas also displaying superiority over the mass penalty method in terms of accuracy and versatility. Copyright © 2012 John Wiley & Sons, Ltd.     

Numerical implementation of a neural network based material model in finite element analysis

Neural network (NN) based constitutive models can capture non‐linear material behaviour. These models are versatile and have the capacity to continuously learn as additional material response data becomes available. NN constitutive models are increasingly used within the finite element (FE) method for the solution of boundary value problems. NN constitutive models, unlike commonly used plasticity models, do not require special integration procedures for implementation in FE analysis. NN constitutive model formulation does not use a material stiffness matrix concept in contrast to the elasto‐plastic matrix central to conventional plasticity based models. This paper addresses numerical implementation issues related to the use of NN constitutive models in FE analysis. A consistent material stiffness matrix is derived for the NN constitutive model that leads to efficient convergence of the FE Newton iterations. The proposed stiffness matrix is general and valid regardless of the material behaviour represented by the NN constitutive model. Two examples demonstrate the performance of the proposed NN constitutive model implementation. Copyright © 2004 John Wiley & Sons, Ltd.     

基于动刚度实验的发动机悬置系统计算模型及模态验证

计算模型的真实性、可靠性是模型建立的关键问题,线性系统动力学模型不能适应参数呈非线性系统的建模。以发动机悬置系统为例,通过实验获得刚度和阻尼参数,用实验的质量矩阵、刚度矩阵建立系统动力学计算模型,通过实验运行模态与计算模态相关分析得到验证模型。结果表明考虑动刚度参数的计算模型比传统静刚度模型更接近真实情形。     

Non-proportional loading in sequentially linear analysis for 3D stress states

Sequentially linear analysis (SLA), a non-incremental-iterative approach towards finite element simulation of quasi-brittle materials, is based on sequentially identifying a critical integration point in the model, to reduce its strength and stiffness, and the associated critical load multiplier (λ_crit), to scale the linear analysis results. In this article, two novel methods are presented to enable SLA simulations for non-proportional loading situations in a three-dimensional fixed smeared crack framework. In the first approach, the cubic function in the load multiplier is analytically solved for real roots using trigonometric solutions or the Cardano method. In the second approach, the load multiplier is expressed as a function of the inclination of a potential damage plane and is deduced using a constrained optimization approach. The first method is preferred over the second for the validation studies due to computational efficiency and accuracy reasons. A three-point bending beam test, with and without prestress, and an RC slab tested in shear, with and without axial loads, are used as benchmarks. The proposed solution method shows good agreement with the experiments in terms of force-displacement curves and damage evolution.     

Application of piece‐wise linear weight functions for 2D 8‐node quadrilateral element in contact problems

The present study is a continuation of our previous work with the aim to reduce problems caused by standard higher order elements in contact problems. The difficulties can be attributed to the inherent property of the Galerkin method which gives uneven distributions of nodal forces resulting in oscillating contact pressures. The proposed remedy is use of piece‐wise linear weight functions. The methods to establish stiffness and/or mass matrix for 8‐node quadrilateral element in 2D are presented, i.e. the condensing and direct procedures. The energy and nodal displacement error norms are also checked to establish the convergence ratio. Interpretation of calculated contact pressures is discussed. Two new 2D 8‐node quadrilateral elements, QUAD8C and QUAD8D, are derived and tested in many examples, which show their good performance in contact problems. Copyright © 2004 John Wiley & Sons, Ltd.     

Out-of-core solver for large,multi-zone boundary element matrices

In this paper an out-of-core solver is developed for a three-dimensional elastostatic boundary element program. This program includes multi-regions which enables a large model to be divided into a number of homogeneous sub-models, each with their own material properties. The system matrix produced by multi-regions is sparse, blocked and unsymmetric in character and so reduces disc space and solution time. The solver presented in this paper utilizes efficiently the structure of this matrix by holding only non-zero parts of this matrix in a sequential, unformatted direct access file, reading in as much of the file as possible into the working space, performing Gauss elimination, then writing part of the matrix back to file. The maximum utilization of working space is particularly important on vector machines as vector activity is maximized whilst record read/write is minimized. The effect of multi-regions on CPU is demonstrated on both the CONVEX and CRAY machines.     

Hierarchical subspace evolution method for super large parallel computing: A linear solver and an eigensolver as examples

Solving linear equations and finding eigenvalues are essential tasks in many simulations for engineering applications, but these tasks often cause performance bottlenecks. In this work, the hierarchical subspace evolution method (HiSEM), a hierarchical iteration framework for solving scientific computing problems with solution locality, is proposed. In HiSEM, the original problem is converted to a corresponding minimization function. The problem is decomposed into a series of subsystems. Subspaces and their weights are established for the subsystems and evolve in each iteration. The subspaces are calculated based on local equations and knowledge of physical problems. A small-scale minimization problem determines the weights of the subspaces. The solution system can be hierarchically established based on the subspaces. As the iterations continue, the degrees of freedom gradually converge to an accurate solution. Two parallel algorithms are derived from HiSEM. One algorithm is designed for symmetric positive definite linear equations, and the other is designed for generalized eigenvalue problems. The linear solver and eigensolver performance is evaluated using a series of benchmarks and a tower model with a complex topology. Algorithms derived from HiSEM can solve a super large-scale problem with high performance and good scalability.     

An incremental 2D constitutive model accounting for linear viscoelasticity and damage development in short fibre composites

A model accounting for linear viscoelasticity and microdamage evolution in short fibre composites is described. An incremental 2D formulation suitable for FE‐simulation is derived and implemented in FE‐solver ABAQUS. The implemented subroutine allows for simulation close to the final failure of the material. The formulation and subroutine is validated with analytical results and experimental data in a tensile test with constant strain rate using sheet moulding compound composites. FE‐simulation of a four‐point bending test is performed using shell elements. The result is compared with linear elastic solution and test data using a plot of maximum surface strain in compression and tension versus applied force. The model accounts for damage evolution due to tensile loading and neglects any damage evolution in compression, where the material has higher strength. Simulation and test results are in very good agreement regarding the slope of the load–strain curve and the slope change. Copyright © 2005 John Wiley & Sons, Ltd.     

An implementation of the stiffness derivative method as a discrete analytical sensitivity analysis and its application to mixed mode in LEFM

In this work, an improvement in the stiffness derivative method based on a shape design sensitivity analysis is proposed, so that the error inherent in the finite difference procedure is avoided. For a global estimation of G from a given finite element solution, this approach is shown to be equivalent to the well-known J-integral when the latter is numerically implemented through its equivalent domain integral. However, it is verified that its direct application to 2D mixed mode problems of linear elastic fracture mechanics through the field decomposition technique yields estimates for G_I and G_II which are in general more accurate for the proposed method. The importance of the velocity field is also remarked and some suggestions for its choice are given.     

An algorithm for the matrix‐free solution of quasistatic frictional contact problems

A contact enforcement algorithm has been developed for matrix‐free quasistatic finite element techniques. Matrix‐free (iterative) solution algorithms such as non‐linear conjugate gradients (CG) and dynamic relaxation (DR) are desirable for large solid mechanics applications where direct linear equation solving is prohibitively expensive, but in contrast to more traditional Newton–Raphson and quasi‐Newton iteration strategies, the number of iterations required for convergence is typically of the same order as the number of degrees of freedom of the model. It is therefore crucial that each of these iterations be inexpensive to per‐form, which is of course the essence of a matrix free method. In applying such methods to contact problems we emphasize here two requirements: first, that the treatment of the contact should not make an average equilibrium iteration considerably more expensive; and second, that the contact constraints should be imposed in such a way that they do not introduce spurious energy that acts against the iterative solver. These practical concerns are utilized to develop an iterative technique for accurate constraint enforcement that is suitable for non‐linear conjugate gradient and dynamic relaxation iterative schemes. Copyright © 1999 John Wiley & Sons, Ltd.     

Experimental observations and finite element modelling of damage initiation and evolution in carbon/epoxy non-crimp fabric composites

Damage in carbon/epoxy non-crimp stitched fabric (NCF) reinforced composites, produced by the resin transfer moulding (RTM) process is described. Formation of the stitching loop results in a certain disturbance of the uniform placement of the fibres. These deviations in fibre placement produce resin-rich zones that can influence the mechanical behaviour of the composite part. Tensile tests on quadriaxial (45°/90°-45°/0°)_s laminates are performed accompanied by acoustic emission (AE) registration and X-ray imaging. Early initiation of damage (matrix cracking) in plies with different fibre orientation has been detected. Damage sites correlate with the resin-rich zones created by the stitching. Finite element (FE) analysis is carried out to develop a model that describes damage of the NCF composites. Numerical multi-level FE homogenization is performed to obtain effective elastic orthotropic properties of NCF composite at micro (unit cell of unidirectional tow) and meso (fabric unit cell) levels. A hierarchical sequence of FE models of different scales is created to analyze in detail the 3D stress state of the NCF composite (meso unit cell). A multi-level submodeling approach is applied during FE analysis. Zones of matrix-dominated damage are predicted. A comparison of non-destructive testing results with computational model is performed. Fracture mechanics parameters of matrix crack are computed and cracks growth stability is studied.     

A new monolithic Newton‐multigrid‐based FEM solution scheme for large strain dynamic poroelasticity problems

This paper presents a new efficient monolithic finite element solution scheme to treat the set of PDEs governing a 2D, biphasic, saturated theory of porous media model with intrinsically coupled and incompressible solid and fluid constituents for infinitesimal and large elastic deformation. Our approach, which inherits some of its techniques from CFD, is characterized by the following aspects: (1) it only performs operator evaluation with no additional Gateaux derivatives. In particular, the computations of the time‐consuming material tangent matrix are not involved here; (2) it solves the non‐linear dynamic problem with no restriction on the strength of coupling; (3) it is more efficient than the linear u v p solver discussed in previous works; (4) it requires weaker derivatives, and hence, lower‐order FE can be tested; and (5) the boundary conditions are reduced, solution independent and more convenient to apply than in the old u v p formulation. For the purpose of validation and comparison, prototypical simulations including analytical solutions are carried out, and at the end, an adaptive time stepping procedure is introduced to handle the rapid change in the numbers of nonlinear iterations that may occur. Copyright © 2016 John Wiley & Sons, Ltd.