附，英文原文Model Order Reduction for Large Scale Engineering
Models Developed in ANSYS
Evgenii B. Rudnyi and Jan G. Korvink IMTEK, Institute of Microsystem Technology
Freiburg University Georges-K¨ohler-Allee, 103 D-79110, Freiburg, Germany {rudnyi,korvink}@imtek.de
http://www.imtek.uni-freiburg.de/simulation/
Abstract. We present the software mor4ansysthat allows engineers to employ modern model reduction
techniques to finite element models developed in ANSYS. We focus on how one extracts the required information from ANSYS and performs model reduction in a C++ implementation that is not dependent on a particular sparse solver. We discuss the computational cost with examples related to structural mechanics and thermal finite element models.
1.Introduction
The model order reduction of linear large-scale dynamic systems is already quite an established area [1]. In many papers (see references in [2]), advantages of model reduction have been demonstrated for a variety of scientific and engineering applications. In the present work, we focus on how engineers can combine this technique with existing commercial finite element software in order to
– Speed up a transient or harmonic analysis,
– Generate automatically compact models for system-level simulation, – Incorporate finite element packages during the design phase.
Model reduction is conventionally applied to a large-scale dynamic system of the first order as follows
Ex˙=Ax+Bu (1.1)
y=Cx
where Aand E are system matrices, Bis the input matrix, Cis the output matrix. The aim of model reduction is to generate a low-dimensional approximation to (1.1) in a similar form
Erz=Arz.+Bru (1.2)
Y=Crz.
that describes well the dependence of the output vector y on the input vector u and so that, at the same time, the dimension of the reduced state vector z is much less than the dimension of the original state vector x
After discretization in space of the partial differential equations describing a user model, a finite element package generally produces a system of ordinary differential equations. At this stage, it is possible to directly apply modern model reduction methods [1]. However, the extraction of the system matrices from a commercial package happens not to be straightforward and here we share our experience on how it can be done with ANSYS [3].
We have chosen the Matrix Market format [4] to represent the reduced model (1.2). We suppose that its simulation will be done in another package, such as Matlab or Mathematica. Functions to work with the reduced model in Mathematica are available at the IMTEK Mathematica Supplement at http://www.imtek.unifreiburg.de/simulation/mathematica/IMSweb/.
The system matrices are high-dimensional and sparse. As a result, the implementation of a model reduction algorithm
usually depends on a particular sparse solver and a storage scheme for sparse matrices. We discuss a C++ interface that allows us to isolate the model reduction and sparse solvers completely for negligible overhead.
Finally, we analyse the computation cost and give the performance results for a few ANSYS models. The comparison of the accuracy of reduced models in respect to the original ANSYS models is given elsewhere [5].
2 mor4ansys
The developed software [6] comprises two almost independent modules (see Fig. 1). The first reads a binary ANSYS file and assembles a dynamic system in the form of Eq (1.1) for first order systems or
Mx¨+Ex˙+Kx=Bu. (2.3)
y=Cx.
for second order systems, where M, E.and K.are the three system matrices. The second module applies the model reduction algorithm to Eq (1.1) or (2.3), that is, it finds a low-dimensional basis V.so that the approximation
x=Vz.+? (2.4)
allows us to reproduce the transient behaviour of the original state vector within the error margin ? . After that, the original equations are projected to the subspace found, for example for Eq (1.2) we have
Er=
VTEV, Ar=
VTAV, Br=
VTB, Cr=CV .
We support three methods to treat second-order systems. When the damping matrix is modeled as Rayleigh damping E=αM+βK, the method from Ref [7] allows us to preserve the coefficients α.and β.as parameters in the reduced model. In the general case, one can choose between the transformation to a first-order system, and second order Arnoldi algorithm (SOAR) [8].
The software can also read, as well as write, the matrices for the original system in the Matrix Market format [4]. A number of model reduction benchmarks has been obtained from ANSYS by means of mor4ansys[9].
Fig. 1. mor4ansysblock-scheme
2.1 Interfacing with ANSYS
The development of the first module happen to be rather difficult because most users of a commercial finite element package do not need the capability to extract the dynamics system in the form of Eq (1.1) or (2.3) and, as a result, this is not a trivial operation.
ANSYS is a huge package and its behavior is not completely consistent. For example, the information described below is not applicable for the fluid dynamics module FLOTRAN.
Our software reads the binary EMATfile with element matrices in order to assemble global system matrices. The file format is documented and ANSYS supplies a library of Fortran subroutines to work with it [10]. An example of how one can use them can be found in the mor4ansyscode [6]. ANSYS has a special command, called a partial solve PSOLVE,
with which one can evaluate element matrices for a given state vector without going through the real solution stage. This allows us to generate an EMAT file efficiently for a given model. However, it was necessary to overcome the following problems:
– The EMAT file does not contain the information about either Dirichlet boundary conditions or equation constraints. They should be extracted separately.
– The EMATfile has a contribution to the load vector from element matrices only. If nodal forces or accelerations are used to apply the load, this information should also be extracted individually. – It is necessary to assemble the global matrices from the element matrices.
During the solution phase, ANSYS can write a binary FULLfile with the assembled system matrices. When we started the development with ANSYS 5.7, this file did not contain the load vector (input matrix). Since then there have been many changes. Since ANSYS 6.0 the FULLfile maintains all the original matrices, the load vector, the Dirichlet and equation constraints in the file. ANSYS 8.0 allows us to make the assembly only and write the FULL file without a real solution phase (equivalent to a partial solution with EMAT). One can now also dump the information from the FULLfile in the Harwell Boeing matrix format. Hence, since ANSYS 8.0, it is possible to use the FULL file efficiently. However, depending on the analysis type the FULLfile may contain not the original stiffness matrix, but rather, a linear combination of system matrices instead.
In the current version of mor4ansys,the EMATfile is employed as the main source to build Eq (1.1) or (2.3). Additional information on the Dirichlet and equation constraints and nodal forces is written in the form of text files by means of ANSYS macros we have developed. The FULLfile can be used to extract the load vector when otherwise this is difficult, for example, as in the case when the acceleration load is used.
ANSYS cannot write several load vectors into the FULL and EMAT files. When multiple-input is to be preserved in Eq (1.1) or (2.3), a user should for each input: – Delete the previously applied load, – Apply a new load, – Generate matrices.
In order to ease this process, the second strategy is also allowed when a user does not have to delete the previous load. In this case, each new load vector contains all the previous vectors and mor4ansyscorrects them at the end of the first phase.
2.2 Running the Model Reduction Algorithm
The Krylov subspaces allow us to obtain a low-dimensional subspace basis for (2.4) with excellent approximating properties by means of a very efficient computation [11,8]. The current version of mor4ansysimplements the block Arnoldi algorithm [11] in order to support multiple inputs, the block size being equal to the number of inputs. Each step of an iterative Krylov subspace algorithm requires us to compute a matrix- vector product, for example, for the first-order system
A?1Eh (2.5)
where h.is some vector. The system matrices are high-dimensional and sparse and one does not compute
A?1explicitly.
The only feasible solution is to solve a linear system of equations for each step as follows
Ag=Eh(2.6)
This constitutes the main computational cost up to the order of the reduced system 30. Later on, the additional cost associated with the orthogonolization process can be also added.
There are many sparse solvers as well as many storage schemes for sparse matrices. Our goal was to implement a model reduction algorithm in a way that does not depend on a particular solver. In addition, we wanted to change solvers at run-time, that is, to allow for run-time polymorphism. As a result, we have chosen the virtual function
Fig. 2. Use of model reduction during design and system-level simulation
mechanism, as its overhead is negligible in our case when the operations by themselves are computationally intensive. Our approach is similar to that in the PETs [12] and Trinilos [13] libraries. The abstract interface is written in terms of relatively low-level functions, as the goal was to cover many different scenarios. The vectors are represented by continuous memory, as they are dense in the case of the Krylov subspaces.
At present, the direct solvers from the TAUCS [14] and UMFPACK [15,16] libraries are supported. The ATLAS library [17] has been used to generate the optimized BLAS. We have found that for many ANSYS models up to 500 000 degrees of freedom the modern direct solvers are quite competitive as the matrix factor fits within 4 Gb of RAM. This allows us to reuse the factorization and achieve good performance.
3 Computational Cost of Model Reduction
We have experimentally observed that for many ANSYS models a reduced model of order 30 is enough to accurately represent the original high-dimensional system [5]. Hence, for simplicity we limit the analysis of the computational cost to this case.
The simulation time of the reduced system comprising 30 equations is very small and we can neglect it. Therefore, for the case when several simulations with different input functions are necessary (the system-level simulation case), the advantage of model reduction is out of the question.
Yet, during the design phase, a reduced model should be generated each time when a user changes the geometry or material properties of the original model. In this case, a reduced model might be used just once. Nevertheless, the model reduction time can be smaller than the simulation time of the original system even in this case. These two different situations are shown in Fig. 2. Below we consider the second case.
Table 1. Computational times on Sun Ultra-80 with 4 Gb of RAM in seconds dimension nnz Stationary Solution in ANSYS 7.0 0.87 2.1 16 304 130 180 590 Stationary solution in ANSYS 8.0 0.63 2.2 15 230 95 150 490 Factoring in TAUCS 0.31 1.3 12 190 91 120 410 generation of the first 30 vectors 0.59 2.7 14 120 120 160 420 4 267 11445 20360 79171 152943 180957 375801 20861 93781 265113 2215638 5887290 7004750 15039875 Let us assume that a direct solver is applicable and the dimension of 30 for the reduced system is sufficient. Then the model reduction time is equal to the time of factoring A in Eq (2.5) and the time required for 30 back substitution steps in Eq (2.6). Table 1 presents computational times for seven ANSYS models where the system matrices are symmetric and positive definite. The first four rows correspond to thermal simulations [18] and the last three to structural mechanics of a bond wire [7].
Each case is specified by its dimension and the number of non zero elements in the stiffness matrix. The time of a stationary solution in ANSYS is given as a reference point. Note that the real simulation time in ANSYS required for the stationary solution is larger than in Table 1 as it includes reading/writing files as well as some other operations. After that is listed the time to factor a matrix by means of a multifrontal solver from the TAUCS library [14] and the time to generate the first 30 vectors. The latter is dominated by the solution of Eq (2.6) by means of back substitution. As the difference to generate the first and thirtieth vectors was less than 10-20%, we can say that the orthogonalization cost was relatively small.
Note that the TAUCS multifrontal solver is even faster than the ANSYS solver. The total time to generate a reduced model is about twice more than that for the stationary solution. At the same time, the reduced model can accurately reproduce any transient and harmonic simulation of the original models within a reasonable frequency range.
The simulation time of a harmonic analysis is the product of solution time for a complex linear system by the number of frequencies needed. The matrix factor cannot be re-used as the linear system to solve depends on frequency. The solution time for a complex linear system is about twice more expensive. Hence model reduction allows us to save simulation time by a factor close to the number of frequencies at which the harmonic response is required. For example, if it is necessary to estimate the transfer function at ten frequencies, then the model reduction plus the simulation of the reduced system is roughly ten times faster than the simulation of the original system.
For the transient simulation, the situation is more difficult to analyse as this depends on the integration strategy. In principle, it is possible to say that the model reduction time above is equivalent to 30 equally spaced timesteps as in this case the same strategy with the re-use of the matrix factor can be applied. However, in our experience, in order to achieve accurate integration results for the examples in Table 1, one either needs at least 600 equally-spaced timesteps or one needs to use adaptive integration schemes where the factor re-use is not possible. In both cases, model reduction plus simulation of the reduced system was more than ten times faster. This shows that model reduction can also be viewed as a fast solver and can be employed even during the optimization phase.
4 Conclusions
We have shown that in the case of the linear dynamics systems (1.1) and (2.3) modern model reduction techniques can speed up finite element transient and harmonic simulation significantly. For nonlinear systems, there are promising theoretical results in the case of polynomial type nonlinearity [19]. Yet, in the nonlinear case in addition to many theoretical problems, it happens that extracting a nonlinear system (1.1) or (2.3) from a commercial finite element tool is a challenge by itself.
Acknowledgment
ANSYS models of the microthruster and the bonded wire have been made by T. Bechtold and J. Lienemann respectively. We would also like to acknowledge an anonymous reviewer for the very helpful comments and suggestions to improve the paper. Partial funding by the DFG project MST-Compact (KO-1883/6), the Italian research council CNR together with the Italian province of Trento PAT, the European Union (grant EU IST-1999-29047, Micropyros) and an operating grant of the University of Freiburg is gratefully acknowledged.
References
1. A. C. Antoulas, D. C. Sorensen. Approximation of Large-Scale Dynamical Systems: An overview. Applied Mathematics & Computer Science, 11(5):1093–1121, 2001.
2. E. B. Rudnyi, J. G. Korvink. Automatic Model Reduction for Transient Simulation of MEMS- based Devices. Sensors Update, 11:3–33, 2002. 3. ANSYS, ANSYS Inc. http://www.ansys.com/
4. R. F. Boisvert, R. Pozo, K. A. Remington. The Matrix Market Exchange Formats: Initial Design. NIST Interim Report 5935, 1996. http://math.nist.gov/MatrixMarket/