Large Sparse Systems Of Equations Research Articles

A conservative MHD scheme on unstructured Lagrangian grids for Z-pinch hydrodynamic simulations

A new algorithm to model resistive magnetohydrodynamics (MHD) in Z-pinches has been developed. Two-dimensional axisymmetric geometry with azimuthal magnetic field Bθ is considered. Discretization is carried out using unstructured meshes made up of arbitrarily connected polygons. The algorithm is fully conservative for mass, momentum, and energy. Matter energy and magnetic energy are managed separately. The diffusion of magnetic field is solved using a derivative of the Symmetric-Semi-Implicit scheme, Livne et al. (1985) [23], where unconditional stability is obtained without needing to solve large sparse systems of equations. This MHD package has been integrated into the radiation-hydrodynamics code MULTI-2D, Ramis et al. (2009) [20], that includes hydrodynamics, laser energy deposition, heat conduction, and radiation transport. This setup allows to simulate Z-pinch configurations relevant for Inertial Confinement Fusion.

Petascale solvers for anisotropic PDEs in atmospheric modelling on GPU clusters

Multi-GPU parallelisation of scalable CG and multigrid solvers for anisotropic PDEs.Efficient matrix-free CUDA implementation minimises global memory access on GPUs.Excellent weak scaling on up to 16,384 nVidia K20X cards (44 mio. cores, Titan, OLCF).PDEs with 5x1011 unknowns can be solved in 1?s to an accuracy of 10 - 5 .Achieved performance of 0.78PFLOPs and memory bandwidth utilisation of more than 40%. Memory bound applications such as solvers for large sparse systems of equations remain a challenge for GPUs. Fast solvers should be based on numerically efficient algorithms and implemented such that global memory access is minimised. To solve systems with trillions ( O ( 10 12 ) ) unknowns the code has to make efficient use of several million individual processor cores on large GPU clusters.We describe the multi-GPU implementation of two algorithmically optimal iterative solvers for anisotropic PDEs which are encountered in (semi-) implicit time stepping procedures in atmospheric modelling. In this application the condition number is large but independent of the grid resolution and both methods are asymptotically optimal, albeit with different absolute performance. In particular, an important constant in the discretisation is the CFL number; only the multigrid solver is robust to changes in this constant. We parallelise the solvers and adapt them to the specific features of GPU architectures, paying particular attention to efficient global memory access. We achieve a performance of up to 0.78 PFLOPs when solving an equation with 0.55 ? 1012unknowns on 16384 GPUs; this corresponds to about 3% of the theoretical peak performance of the machine and we use more than 40% of the peak memory bandwidth with a Conjugate Gradient (CG) solver. Although the other solver, a geometric multigrid algorithm, has a slightly worse performance in terms of FLOPs per second, overall it is faster as it needs less iterations to converge; the multigrid algorithm can solve a linear PDE with half a trillion unknowns in about one second.

Differentiating the Method of Conjugate Gradients

The method of conjugate gradients (CG) is widely used for the iterative solution of large sparse systems of equations $Ax=b$, where $A\in\Re^{n\times n}$ is symmetric positive definite. Let $x_k$ denote the $k$th iterate of CG. This is a nonlinear differentiable function of $b$. In this paper we obtain expressions for $J_k$, the Jacobian matrix of $x_k$ with respect to $b$. We use these expressions to obtain bounds on $\|J_k\|_2$, the spectral norm condition number of $x_k$, and discuss algorithms to compute or estimate $J_kv$ and $J_k^Tv$ for a given vector $v$.

A Preconditioner for the Finite Element Approximation to the Arbitrary Lagrangian–Eulerian Navier–Stokes Equations

This paper presents an efficient numerical solver for the finite element approximation of the incompressible Navier–Stokes equations within a moving three-dimensional domain. The moving domain is modeled using the arbitrary Lagrangian–Eulerian (ALE) formulation. Applying a finite element approximation leads to the solution of a large sparse system of equations. In this work we look at the application of the $Fp$ preconditioner within GMRES for efficiently solving such systems. Numerical results are presented for tetrahedral and hexahedral elements, and both structured and unstructured meshes. In all cases GMRES convergence rates are seen to be independent of mesh size. Finally, we show an application of this problem for modeling fluid flow within the heart.

Calculation of the eigenfunctions of the two-group neutron diffusion equation and application to modal decomposition of BWR instabilities

In this paper, numerical methods aiming at determining the eigenfunctions, their adjoint and the corresponding eigenvalues of the two-group neutron diffusion equations representing any heterogeneous system are investigated. First, the classical power iteration method is modified so that the calculation of modes higher than the fundamental mode is possible. Thereafter, the explicitly-restarted Arnoldi method, belonging to the class of Krylov subspace methods, is touched upon. Although the modified power iteration method is a computationally-expensive algorithm, its main advantage is its robustness, i.e. the method always converges to the desired eigenfunctions without any need from the user to set up any parameter in the algorithm. On the other hand, the Arnoldi method, which requires some parameters to be defined by the user, is a very efficient method for calculating eigenfunctions of large sparse system of equations with a minimum computational effort. These methods are thereafter used for off-line analysis of the stability of boiling water reactors in a two-dimensional representation of the core. Since several oscillation modes are usually excited (global and regional oscillations) when unstable conditions are encountered, the characterization of the stability of the reactor using for instance the Decay Ratio as a stability indicator might be difficult if the contribution from each of the modes are not separated from each other. Such a modal decomposition is applied to a stability test performed at the Swedish Ringhals-1 unit in September 2002, after the use of the Arnoldi method for pre-calculating the different eigenmodes of the neutron flux throughout the reactor. The modal decomposition clearly demonstrates the excitation of both the global and regional oscillations. Furthermore, such oscillations are found to be intermittent with a time-varying phase shift between the first and second azimuthal modes.

Generalized Rigid and Generalized Affine Image Registration and Interpolation by Geometric Multigrid

Generalized rigid and generalized affine registration and interpolation obtained by finite displacements and by optical flow are here developed variationally and numerically as well as with respect to a geometric multigrid solution process. For high order optimality systems under natural boundary conditions, it is shown that the convergence criteria of Hackbusch (Iterative Solution of Large Sparse Systems of Equations. Springer, Berlin, 1993) are met. Specifically, the Galerkin formalism is used together with a multi-colored ordering of unknowns to permit vectorization of a symmetric successive over-relaxation on image processing systems. The geometric multigrid procedure is situated as an inner iteration within an outer Newton or lagged diffusivity iteration, which in turn is embedded within a pyramidal scheme that initializes each outer iteration from predictions obtained on coarser levels. Differences between results obtainable by finite displacements and by optical flows are elucidated. Specifically, independence of image order can be shown for optical flow but in general not for finite displacements. Also, while autonomous optical flows are used in practice, it is shown explicitly that finite displacements generate a broader class of registrations. This work is motivated by applications in histological reconstruction and in dynamic medical imaging, and results are shown for such realistic examples.

Geometric multigrid for high-order regularizations of early vision problems

The surface estimation problem is used as a model to demonstrate a framework for solving early vision problems by high-order regularization with natural boundary conditions. Because the application of algebraic multigrid is usually constrained by an M-matrix condition which does not hold for discretizations of high-order problems, a geometric multigrid framework is developed for the efficient solution of the associated optimality systems. It is shown that the convergence criteria of Hackbusch [W. Hackbusch, Iterative Solution of Large Sparse Systems of Equations, Springer, 1993] are met, and in particular the general elliptic regularity required is proved. Further, the Galerkin formalism is used together with a multicolored ordering of unknowns to permit vectorization of a symmetric Gauss–Seidel relaxation in image processing systems. The implementation is analyzed computationally and inaccuracies are corrected by lumping and by proper floating point representations. Direct one-dimensional calculations are used to estimate the effect of regularization order, regularization strength, relaxation, and data support on the multigrid reduction factor. A finite difference formulation is ruled out in favor of a finite element formulation. A representative problem from magnetic resonance coil sensitivity estimation is solved using increasingly higher orders of regularization, and the results are compared in terms of accuracy and multigrid convergence.

Iterative solution of fuzzy linear systems

Linear systems have important applications to many branches of science and engineering. In many applications, at least some of the parameters of the system are represented by fuzzy rather than crisp numbers. So it is immensely important to develop numerical procedures that would appropriately treat general fuzzy linear systems (will be denoted by FLS) and solve them. In this paper firstly a general fuzzy linear system using the embedding approach, has been investigated and then several well-known numerical algorithms for solving system of linear equations such as Richardson, Extrapolated Richardson, Jacobi, JOR, Gaus †–Seidel, EGS, SOR, AOR, ESOR, SSOR, USSOR, EMA and MSOR are extended for solving FLS. The iterative methods are followed by convergence theorems and the presented algorithms are tested by solving some numerical examples. ( †Hackbusch noticed that Gauss spelling is less correct than Gaus [W. Hackbusch, Iterative Solution of Large Sparse Systems of Equations, Springer-Verlag Inc., New York, 1994 [9]].)

Identification of boundary heat fluxes in a falling film experiment using high resolution temperature measurements

In this paper, we consider a three-dimensional inverse heat conduction problem (IHCP) in a falling film experiment. The wavy film is heated electrically by a thin constantan foil and the temperature on the back side of this foil is measured by high resolution infrared (IR) thermography. The transient heat flux at the inaccessible film side of the foil is determined from the IR data and the electrical heating power. The IHCP is formulated as a mathematical optimization problem, which is solved with the conjugate gradient method. In each step of the iterative process two direct transient heat conduction problems must be solved. We apply a one step θ-method and piecewise linear finite elements on a tetrahedral grid for the time and space discretization, respectively. The resulting large sparse system of equations is solved with a preconditioned Krylov subspace method. We give results of simulated experiments, which illustrate the performance and tuning of the solution method, and finally present the estimation results from temperature measurements obtained during falling film experiments.

On the Approximate Cyclic Reduction Preconditioner

We present a preconditioning method for the iterative solution of large sparse systems of equations. The preconditioner is based on ideas both from ILU preconditioning and from multigrid. The resulting preconditioning technique requires the matrix only. A multilevel structure is obtained by using maximal independent sets for graph coarsening. A Schur complement approximation is constructed using a sequence of point-Gaussian elimination steps. The resulting preconditioner has a transparent modular structure similar to the algorithmic structure of a multigrid V-cycle.

Application of a factored newton-relaxation scheme to calculation of discrete aerodynamic sensitivity derivatives*

The objective of this paper is to describe the modifications and extensions required to apply an existing Euler solver to the calculation of discrete aerodynamic sensitivity derivatives and to evaluate the perfomance of a Newton-relaxation scheme as applied to the very large sparse systems of equations resulting from aerodynamics sensitivity analysis. The preponderance of prior work has utilized direct methods to solve these linear equations resulting in techniques requiring large amounts of computer memory making application to complex 2D or 3D configurations difficult, if not impractical. The algorithm applied to invert the aerodynamic sensitivity system of equations was originally developed for the Euler/Navier-Stokes CFD codes and applied to numerous complex two- and three-dimensional configurations. The delta form of the aero-dynamic sensitivity equations was used to enable efficient coupling with CFD codes using high spatial resolution schemes, viscous effects, and multi-block decomposition. Sample computations are included for estimation of sensitivity derivatives on a 2D airfoil at various points in design space. The design space gradients are compared with gradient estimations from numerical derivatives. Coupling calculated design space gradients with nonlinear optimization algorithms will produce a design optimization tool applicable to various aerodynamic devices including airfoils, compressor cascades, and propulsion test facility forebody simulators. Development, application and evaluation of design optimization capability using these gradients will be the subject of future reports.

Iterative Solution of Large Sparse Systems of Equations.

Iterative Solution of Large Sparse Systems of Equations.

A multifrontal algorithm for the solution of large systems of equations using network-based parallel computing

A multifrontal algorithm for solving large sparse systems of equations has been developed and implemented in a distributed computing environment. The Parallel Virtual Machine (PVM) software was used to combine several networked workstations into a single parallel computational resource. The motivation for this work was the need for solving efficiently the large systems of equations arising in Finite Element Analysis of transport and reaction processes underlying the chemical vapor deposition of thin films. Both out-of-core and in-core computations produce significant speedups for a flow and heat transfer problem using up to 7 networked workstations in a parallel mode. The performance of the algorithm was found to deteriorate rapidly as the granularity of the problem decreased. Thus, it is desirable to scale up the problem size with the number of processors to achieve respectable speedups with the proposed multifrontal technique. The intrinsic fault tolerance capabilities of the PVM software were implemented for completing the computations when one or more workstations became unavailable during a run.

Solution of large sparse equation systems by parallel-section method on es computers

Solution of large sparse equation systems by parallel-section method on es computers

Conjugate direction relaxation solutions for 3-D magnetotelluric modeling

In recent years, there has been a tremendous amount of progress made in three‐dimensional (3-D) magnetotelluric modeling algorithms. Much of this work has been devoted to the integral equation technique (e.g., Hohmann, 1975; Weidelt, 1975; Wannamaker et al., 1984; Wannamaker, 1991). This method has contributed significantly to our understanding of electromagnetic field behavior in 3-D models. However, some of the very earliest work in 3-D modeling concentrated on differential methods (e.g., Jones and Pascoe, 1972; Reddy et al., 1977). It is generally recognized that differential methods are better suited than integral equation methods to model arbitrarily complex geometries, and consequently this area has recently been receiving a great deal of attention (e.g., Madden and Mackie, 1989; Xinghua et al., 1991; Mackie et al., 1993; Smith, 1992, personnal communication). Differential methods lead to large sparse systems of equations to be solved for the unknown field values. It is possible to use relaxation algorithms to quickly obtain approximate solutions to these systems of equations without resorting to standard matrix inversion routines or sparse matrix solvers.

Parallel ICCG on a hierarchical memory multiprocessor — Addressing the triangular solve bottleneck

The incomplete Cholesky conjugate gradient (ICCG) algorithm is a commonly used iterative method for solving large sparse systems of equations. In this paper, we study the parallel solution of sparse triangular systems of equations, the most difficult aspect of implementing the ICCG method on a multiprocessor. We focus on shared-memory multiprocessor architectures with deep memory hierarchies. On such architectures we find that previously proposed parallelization approaches result in little or no speedup. The reason is that these approaches cause significant increases in the amount of memory system traffic as compared to a sequential approach. Increases of as much as a factor of 10 on four processors were observed. In this paper we propose new techniques for limiting these increases, including data remappings to increase spatial locality, new processor synchronization techniques to decrease the use of auxiliary data structures, and data partitioning techniques to reduce the amount of interprocessor communication. With these techniques, memory system traffic is reduced to as little as one-sixth of its previous volume. The resulting speedups are greatly improved as well, although they are still much less than linear. We discuss the factors that limit further speedups. We present both simulation results and results of experiments on an SGI 4D/340 multiprocessor.

Numerical simulation of mechanical systems using methods for differential‐algebraic equations

AbstractA simple approach to the numerical simulation of mechanical systems consisting of rigid and flexible bodies is presented. The mechanical system may consist of rigid bodies, different types of flexible bodies, joints and actuators and may have arbitrary topological structure with kinematical loops. The equations of motion are formulated as a large sparse system of equations in absolute co‐ordinates as well as relative co‐ordinates. These equations are numerically integrated as a system of differential‐algebraic equations using modern numerical methods.

Studies of finite element procedures—the conjugate gradient and GMRES methods in adina and ADINA-F

In this paper we present some experiences with the use of the conjugate gradient and GMRES iterative methods for the solution of large sparse systems of equations as recently implemented in ADINA and ADINA-F for structural and fluid flow problems. The conjugate gradient method preconditioned with the incomplete Cholesky decomposition is used for the symmetric positive definite systems resulting from structural problems. For fluid flow problems, the systems of equations are nonsymmetric and indefinite. For these systems, the biconjugate gradient and GMRES algorithms with preconditioning by incomplete LU factorization are used. The performance of the iterative methods is compared with the direct solution methods. The results from our numerical experiments show that the use of these iterative methods for large sparse systems can lead to significant reductions in storage requirements and computation times, especially for nonlinear structural dynamics problems and three-dimensional problems in general.

Some issues in solving large sparse systems of equations

In Don and Gallo (1987) it was argued that a large sparse system of equations can be solved efficiently if the equations are ordered in such a way that the set of feedback variables is small. Also an algorithm to find such an ordering was presented. In this paper we discuss some complications that may arise in applying these techniques. The issues considered have come up in applied econometric models and relate to nonnormalized equations, implicit equations, and vector equations. It is argued that a good normalization does not always exist (even if a bad one does). We discuss how implicit equations, whether single noninvertible equations or system equilibrium conditions, can be handled. Finally it is shown that vector equations pose no real problem if they are appropriately represented in the structural graph.

Some Issues in Solving Large Sparse Systems of Equations

In Don and Gallo (1987) it was argued that a large sparse; system of equations can be solved efficiently if the equations are ordered in such a way that the set of feedback variables is small. Also an algorithm to find such an ordering was presented.In this paper we discuss some complications that may arise in applying these techniques. The issues considered have come up in applied econometric models and relate to non-normalized equations, implicit equations and vector equations. It will be argued that a good normalization does not always exist (even if a bad one does), and a respecification may be required. We discuss how implicit equations, whether single non-invertible equations or system equilibrium conditions, can be handled. Finally it is shown that vector equations pose no real problem if they are appropriately represented in the structural graph.