Sparse Solver Research Articles

Trefftz finite elements for electromagnetics

It is shown that the method of minimum autonomous blocks (MAB) of Nikol'skii and Nikol'skaia can be reformulated as the Trefftz finite-element method. Solutions of Maxwell's equations in form of plane waves are used to represent fields inside a finite element. Their projections on a set of basis functions on the surface of the element are used to obtain a descriptor of the element in form of an admittance matrix. It is shown that a point-matching projection procedure gives the frequency-domain transmission-line-matrix formulation and Galerkin-type projection leads to the MAB formulation. Admittance matrix representation of the descriptors of the elements makes it possible to use a finite-element-type global matrix assembling procedure and a sparse matrix solver.

Finite‐element modeling of multibody contact and its application to active faults

AbstractEarthquakes have been recognized as resulting from a stick–slip frictional instability along the faults between deformable rocks. An arbitrarily‐shaped contact element strategy, named the node‐to‐point contact element strategy, is proposed, applied with the static‐explicit characters to handle the friction contact between deformable bodies with stick and finite frictional slip and extended here to simulate the active faults in the crust with a more general nonlinear friction law. An efficient contact search algorithm for contact problems among multiple small and finite deformation bodies is also introduced. Moreover, the efficiency of the parallel sparse solver for the nonlinear friction contact problem is investigated. Finally, a model for the plate movement in the north‐east zone of Japan under gravitation is taken as an example to be analyzed with different friction behaviors. Copyright © 2002 John Wiley & Sons, Ltd.

Algorithms and Application of Sparse Matrix Assembly and Equation Solvers for Aeroacoustics

An algorithm for symmetric sparse equation solutions on an unstructured grid is described. Efficient, sequential sparse algorithms for degree-of-freedom reordering, supernodes, symbolic/numerical factorization, and forward backward solution phases are reviewed. Three sparse algorithms for the generation and assembly of symmetric systems of matrix equations are presented. The accuracy and numerical performance of the sequential version of the sparse algorithms are evaluated over the frequency range of interest in a three-dimensional aeroacoustics application. Results show that the solver solutions are accurate using a discretization of 12 points per wavelength. Results also show that the first assembly algorithm is impractical for high-frequency noise calculations. The second and third assembly algorithms have nearly equal performance at low values of source frequencies, but at higher values of source frequencies the third algorithm saves CPU time and RAM. The CPU time and the RAM required by the second and third assembly algorithms are two orders of magnitude smaller than that required by the sparse equation solver. A sequential version of these sparse algorithms can, therefore, be conveniently incorporated into a substructuring for domain decomposition formulation to achieve parallel computation, where different substructures are handles by different parallel processors.

Three-Dimensional Rectangular Duct Code with Application to Impedance Eduction

A zero flow, fully three-dimensional, variable impedance, rectangular duct aeroacoustics code that spans the frequency spectrum of interest in duct liner research is developed. The governing equations and boundary conditions in the duct are solved numerically using the finite element methodology. The methodology makes use of a state-of-the-art, sparse equation solver to obtain the capability to study high-frequency sound waves that may require millions of grid points for resolution. Noise suppression levels predicted from the code are in excellent agreement with those obtained from mode theory. The single-processor performance of the solver, relative to that of the more commonly used band solver, increases with frequency. At a frequency of 17 kHz, the band solver is 4.25 times slower and consumes 2.5 times more memory than the fully sparse equation solver. The duct aeroacoustics code is combined with an optimization algorithm and used successfully to educe the impedance spectrum of a ceramic liner. The primary problem with using the methodology to perform optimization studies at frequencies above 14 kHz is excessive central processor unit time. The results support the recommendation that research be directed toward exploitation of the multiprocessor capability of the solver to further reduce central processor unit time.

Improved Symbolic and Numerical Factorization Algorithms for Unsymmetric Sparse Matrices

We present algorithms for the symbolic and numerical factorization phases in the direct solution of sparse unsymmetric systems of linear equations. We have modified a classical symbolic factorization algorithm for unsymmetric matrices to inexpensively compute minimal elimination structures. We give an efficient algorithm to compute a near-minimal data-dependency graph for unsymmetric multifrontal factorization that is valid irrespective of the amount of dynamic pivoting performed during factorization. Finally, we describe an unsymmetric-pattern multifrontal algorithm for Gaussian elimination with partial pivoting that uses the task- and data-dependency graphs computed during the symbolic phase. These algorithms have been implemented in WSMP---an industrial strength sparse solver package---and have enabled WSMP to significantly outperform other similar solvers. We present experimental results to demonstrate the merits of the new algorithms.

A 3D incompressible Navier–Stokes velocity–vorticity weak form finite element algorithm

AbstractThe velocity–vorticity formulation is selected to develop a time‐accurate CFD finite element algorithm for the incompressible Navier–Stokes equations in three dimensions.The finite element implementation uses equal order trilinear finite elements on a non‐staggered hexahedral mesh. A second order vorticity kinematic boundary condition is derived for the no slip wall boundary condition which also enforces the incompressibility constraint. A biconjugate gradient stabilized (BiCGSTAB) sparse iterative solver is utilized to solve the fully coupled system of equations as a Newton algorithm. The solver yields an efficient parallel solution algorithm on distributed‐memory machines, such as the IBM SP2. Three dimensional laminar flow solutions for a square channel, a lid‐driven cavity, and a thermal cavity are established and compared with available benchmark solutions. Copyright © 2002 John Wiley & Sons, Ltd.

Analysis and comparison of two general sparse solvers for distributed memory computers

This paper provides a comprehensive study and comparison of two state-of-the-art direct solvers for large sparse sets of linear equations on large-scale distributed-memory computers. One is a multifrontal solver called MUMPS, the other is a supernodal solver called superLU. We describe the main algorithmic features of the two solvers and compare their performance characteristics with respect to uniprocessor speed, interprocessor communication, and memory requirements. For both solvers, preorderings for numerical stability and sparsity play an important role in achieving high parallel efficiency. We analyse the results with various ordering algorithms. Our performance analysis is based on data obtained from runs on a 512-processor Cray T3E using a set of matrices from real applications. We also use regular 3D grid problems to study the scalability of the two solvers.

Wavelet matrix transform approach for electromagnetic scattering by a dielectric cylinder

The wavelet matrix transform approach, in combination with the conventional method of moments, is used to solve the problem of electromagnetic scattering by a dielectric cylinder. First, the problem is formulated in terms of the integral equation for the field of a harmonic source in the presence of the cylinder. Then the equation is discretized to obtain a dense system of linear equations by dividing the cylinder into square cells which are small enough so that the field intensity is nearly uniform in each cell. Finally, the wavelet matrix transform is applied to produce a sparse system of linear equations which is treated effectively by a sparse linear system solver. © 2001 Scripta Technica, Electr Eng Jpn, 137(3): 1–9, 2001

A posteriori finite element output bounds in three space dimensions using the FETI method

A domain decomposition finite element technique for efficiently generating lower and upper bounds to outputs which are linear functionals of the solution to the convection–diffusion equation is presented. The novelty is to utilize a domain decomposition method with Lagrange multipliers and an iterative interface solution scheme – the finite element tearing and interconnecting (FETI) procedure – to extend the bound method to three space dimensions. Numerical results show that the bounds are sharp, their calculation consumes only a fraction of the memory required by a sparse solver, and CPU usage is optimal, i.e. proportional to the total number of d.o.f.

PARDISO: a high-performance serial and parallel sparse linear solver in semiconductor device simulation

The package PARDISO is a high-performance, robust and easy to use software for solving large sparse symmetric or structurally symmetric linear systems of equations on shared memory multiprocessors. PARDISO uses a combination of left- and right-looking Level-3 BLAS supernode techniques to exploit pipelining parallelism. It delivers up to 960 Mflop/s on COMPAQ Alpha ES40 (667 MHz) for irregular problems and sparse matrix factorization has been clocked up at a speedup of 7 on an 8-node SGI Origin 2000. The paper gives an overview of the algorithm, performance results and the integration of the solver into complex industrial simulation tools. Finally, an example is discussed inherently (due to the design goal) producing linear systems close to singularity.

Wavelet matrix transform approach for electromagnetic scattering by a dielectric cylinder

AbstractThe wavelet matrix transform approach, in combination with the conventional method of moments, is used to solve the problem of electromagnetic scattering by a dielectric cylinder. First, the problem is formulated in terms of the integral equation for the field of a harmonic source in the presence of the cylinder. Then the equation is discretized to obtain a dense system of linear equations by dividing the cylinder into square cells which are small enough so that the field intensity is nearly uniform in each cell. Finally, the wavelet matrix transform is applied to produce a sparse system of linear equations which is treated effectively by a sparse linear system solver. © 2001 Scripta Technica, Electr Eng Jpn, 137(3): 1–9, 2001

An MEBDF package for the numerical solution of large sparse systems of stiff initial value problems

An efficient algorithm for the numerical integration of large sparse systems of stiff initial value ordinary differential equations (ODEs) and differential-algebraic equations (DAEs) is described. The algorithm is constructed by embedding a standard sparse linear algebraic equation solver into a suitably modified MEBDF code. An important practical application of this algorithm is in the numerical solution of time dependent partial differential equations (PDEs), particularly in two or more space dimensions, using the method of lines (MOL). A code based on this algorithm is illustrated by application to several problems of practical interest and its performance is compared to that of the standard code LSODES.

Advanced FI2TD algorithms for transient eddy current problems

Transient eddy current formulations based on the Finite Integration Technique (FIT) for the magneto‐quasistatic regime are extended to include motional induction effects of moving conductors with simple geometries by different approaches. A new regularization of the formulation using discrete grad‐div augmentation of the curlcurl formulation is presented and tested. To improve the implicit time integration process, several schemes for an error controlled variable time step selection are presented and for the repetitive solution of the arising large sparse systems of equations a sparse direct solver is compared to iterative methods such as a preconditioned conjugate gradient method and a new algebraic multigrid solver, which is aware of the curlcurl nullspace.

Solving Large Nonconvex Water Resources Management Models Using Generalized Benders Decomposition

Nonconvex nonlinear programming (NLP) problems arise frequently in water resources management, e.g., reservoir operations, groundwater remediation, and integrated water quantity and quality management. Such problems are usually large and sparse. Existing software for global optimization cannot cope with problems of this size, while current local sparse NLP solvers, e.g., MINOS (Murtagh and Saunders 1987), or CONOPT (Drud 1994) cannot guarantee a global solution. In this paper, we apply the Generalized Benders Decomposition (GBD) algorithm to two large nonconvex water resources models involving reservoir operations and water allocation in a river basin, using an approximation to the GBD cuts proposed by Floudas et al. (1989) and Floudas (1995). To ensure feasibility of the GBD subproblem, we relax its constraints by introducing elastic slack variables, penalizing these slacks in the objective function. This approach leads to solutions with excellent objective values in run times much less than the GAMS NLP solvers MINOS5 and CONOPT2, if the complicating variables are carefully selected. Using these solutions as initial points for MINOS5 or CONOPT2 often leads to further improvements.

Chronopotentiometry at a microband electrode: simulation study using a Rosenbrock time integration scheme for differential–algebraic equations and a direct sparse solver

Two-dimensional digital simulation of chronopotentiometry at a microband electrode cannot be performed by conventional alternating direction implicit finite-difference simulation methods, as a result of non-local boundary conditions at the electrode. It is also potentially troublesome for iterative Krylov algorithms for solving linear algebraic equations that result from implicit temporal integration schemes. These difficulties can be avoided by representing the spatially discretised initial boundary value problem in the form of a set of differential–algebraic equations, to which a suitable integration scheme, such as the Rosenbrock ROWDA3 scheme, can be applied. The resulting linear algebraic equations can be solved unfailingly by a general direct sparse algorithm such as Y12M. Using this technique it has been found that the chronopotentiometric potential–time curves and transition time characteristics of a single charge transfer reaction at a microband electrode resemble those for a microhemicylinder electrode of the same surface area, although transition times are somewhat shorter. The simulation reveals also that the assumption of a uniform flux of the concentration at the electrode, encountered in the literature in the context of the two-dimensional theories of chronopotentiometry at microelectrodes, is not adequate, owing to non-negligible edge effects observed.

The use of new SQP methods for the optimization of utility systems

An equation-oriented (EO) mathematical model for the optimization of utility systems has been developed. The model includes an accurate physical description of water/steam and air process streams. These together with EO models for the required unit operations were used to optimize two utility systems: a small model containing two steam turbines and a larger model of a combined heat and power plant. The applicability of the model was established using a novel SQP (Sequential Quadratic Programming) method which employs a ‘tolerance tube’ approach (Zoppke-Donaldson , A Tolerance-Tube Approach to Sequential Quadratic Programming with Applications, PhD thesis, 1995) to decide whether to accept a step from the QP subproblem. Comparison was then made with a further new penalty-free optimization method named filterSQP (Nonlinear Programming without a Penalty Function, 1997), which was used to solve the model containing two steam turbines. The tuning of both the model and the use of filterSQP were further explored with numerical experiments on the effect of scaling; comparison of dense and sparse versions of filterSQP; and different initial sizes for the trust region. Also, the correctness is demonstrated of the ‘shadow price’ given by the SQP code for the steam generation temperature. The results show that the sparse filterSQP solver works very efficiently and provides a basis for modelling and optimization of larger, industrial size problems which are reported in a companion paper.

A hybrid BEM/WTM approach for analysis of the EM scattering from large open-ended cavities

An effective hybrid boundary-element method (BEM) and wavelet-transform method (WTM) is proposed to analyze electromagnetic scattering from three-dimensional (3-D) open-ended cavities with arbitrary shapes. This hybrid technique formulates the original cavity problems by a magnetic field integral equation. The BEM is employed to establish the mapping between the original complex integral surface and the unit square. The WTM is used to reduce the density of the moment matrix. Since a surface integral equation has to be solved, the WTM requires a two-dimensional (2-D) wavelet basis to implement the numerical computation. The previous fast iterative algorithm for 2-D wavelets has been extended for efficiently constructing various 2-D wavelet basis functions by a tensorial product from two one-dimensional (1-D) regular multiresolution analyses. Unlike the conventional method of moments, the proposed hybrid technique can always obtain sparse moment matrix equations, which can be efficiently solved by sparse solvers. As the level scales for numerical discretization of cavities increase, larger compression rates can be obtained, which makes it possible for the hybrid BEM/WTM technique to efficiently solve scattering from large open-ended cavities with complex terminations. Numerical results are presented to demonstrate the merits of the proposed method.

Adapting EMAP3D to parallel processing [for EM launcher modelling

EMAP3D (Electromechanical Analysis Program in Three Dimensions) is a single-processor (serial) program that uses finite element (FE) methods to solve coupled electromagnetic, thermal and structural problems for high velocity conductors in transient electromagnetic fields. Its primary application has been the simulation of electromagnetic launchers and pulsed rotating machines. While ETVAP3D has been applied successfully to a wide range of problems, its serial execution time limits the achievement of finer detail in large problems, even on high performance vector machines. The present class of production simulations, involving 100,000 unknowns, can take a week to complete on time-sharing machines at computation centers. To reduce simulation time, a parallel implementation of the EMAPSD program has been undertaken. While vector parallel programming is straightforward and a reasonable approach, access to more than 16 vector processors is limited. The availability, low cost, and scalability of massively parallel processing (MPP) make MPP computing more attractive than vector-parallel processing. The number of MPP processors on PC Beowulf clusters usually ranges from 8 to 100 and can be as high as several thousand on IBM (SP) and Gray (T3E) systems. The authors have decoupled the matrix generation and components of the FE algorithm, thereby allowing us to use any of the new MPP-parallel solvers on the matrix equations. Since the number of zero matrix elements is high for EMAPSD problems, a sparse matrix solver is ideal. Hence, their parallel implementation uses the sparse iterative solvers of PETSc (Portable Extensible Toolkit for Scientific Computing). Here, they report the performance, scalability and use of PETSc preconditioners and solver algorithms as a solution engine for real EMAP3D simulations and test cases.

Design of generic direct sparse linear system solver in C++ for power system analysis

Design of generic direct sparse linear system solver in C++ for power system analysis

MPI-based implementation of a PCG solver using an EBE architecture and preconditioner for implicit, 3-D finite element analysis

This work describes a coarse-grain parallel implementation of a linear preconditioned conjugate gradient solver using an element-by-element architecture and preconditioner for computation. The solver, implemented within a nonlinear, implicit finite element code, uses an MPI-based message-passing approach to provide portable parallel execution on shared, distributed, and distributed-shared memory computers. The flexibility of the element-by-element approach permits a dual-level mesh decomposition; a coarse, domain-level decomposition creates a load-balanced domain for each processor for parallel computation, while a second level decomposition breaks each domain into blocks of similar elements (same constitutive model, order of integration, element type) for fine-grained parallel computation on each processor. The key contribution here is a new parallel implementation of the Hughes–Winget (HW) element-by-element preconditioner suitable for arbitrary, unstructured meshes. The implementation couples an unstructured dependency graph with a new balanced graph-coloring algorithm to schedule parallel computations within and across domains. The code also includes the diagonal preconditioner and a modern parallel (threaded) sparse direct solver for comparison. Three example problems with up to 158,000 elements and 180,000 nodes analyzed on an SGI/Cray Origin 2000 illustrate the parallel performance of the algorithms and preconditioners. Analyses with varying block sizes illustrate that the two-level decomposition improves overall execution speed with the block size tuned for the cache memory architecture of the executing platform. This implementation of the HW preconditioner shows reasonable parallel efficiency – typically 80% on 48 processors. Efficiency for the diagonal preconditioner is also high, with total speedups reaching 86% on 48 CPUs. Calculation of the tangent element stiffnesses shows superlinear speedups for each of the test problems, while the computation of strains/stresses/residual forces shows 80% parallel efficiency on 48 processors.