Engineering Papers⌕ Search

SEARCH · Engineering Papers

Results for “Sparse direct solver”

Search indexed NASA NTRS and DOE OSTI research on propulsion, heat transfer, battery materials and energy systems. Follow report and document links to the original sources.

Quote a phrase for an exact phrase match. Source license links do not imply unrestricted reuse.

At least 55 records · Page 3
All sourcesNASA NTRSDOE OSTI
NASA NTRS1993

Implementing abstract multigrid or multilevel methods

Multigrid methods can be formulated as an algorithm for an abstract problem that is independent of the partial differential equation, domain, and discretization method. In such an abstract setting, problems not arising from partial differential equations can be treated. A general theory exists for linear problems. The general theory was motivated by a series of abstract solvers (Madpack). The latest version was motivated by the theory. Madpack now allows for a wide variety of iterative and direct solvers, preconditioners, and interpolation and projection schemes, including user callback ones. It allows for sparse, dense, and stencil matrices. Mildly nonlinear problems can be handled. Also, there is a fast, multigrid Poisson solver (two and three dimensions). The type of solvers and design decisions (including language, data structures, external library support, and callbacks) are discussed. Based on the author's experiences with two versions of Madpack, a better approach is proposed. This is based on a mixed language formulation (C and FORTRAN + preprocessor). Reasons for not using FORTRAN, C, or C++ (individually) are given. Implementing the proposed strategy is not difficult.

Douglas, Craig C.↗
DOE OSTI2025

Predicting Flow in Fracture Networks With Quantum Algorithms

Uncertainty quantification plays a crucial role in the modeling of subsurface flow. For instance, uncertainties in the properties of geologic fracture networks significantly impact flow, requiring numerous simulations to accurately estimate quantities of interest. However, each simulation is computationally expensive because it requires solving a large linear system to capture features that involve both small and large fractures. An example is in percolation, where the interaction of many small fractures (which cumulatively can have a large surface area) with the rock matrix must be modeled precisely. Quantum computing is an emerging tool with the potential to address this issue. Quantum algorithms offer a significant speedup in solving linear systems, achieving efficiencies that are challenging to match with classical approaches. These classical approaches include direct solvers, such as LU decomposition, and iterative methods, notably preconditioned conjugate gradient, commonly used in subsurface modeling to solve large sparse systems. However, applying quantum algorithms to geologic fracture flow requires careful attention to algorithmic and problem-specific constraints to fully realize this quantum advantage. In this work we describe a quantum algorithm for generalized Monte Carlo applications with a quadratic speedup over the classical approaches which can be combined with the quantum speedup, currently under investigation, for solving quantum linear systems for subsurface flow. We show that for quantum algorithms the computational cost of estimating a quantity of interest for a statistical ensemble of networks is roughly the same as that of a single realization, essentially implying that one can get uncertainty quantification for free.

58 GEOSCIENCES↗
NASA NTRS2014

Modeling Sound Propagation Through Non-Axisymmetric Jets

A method for computing the far-field adjoint Green's function of the generalized acoustic analogy equations under a locally parallel mean flow approximation is presented. The method is based on expanding the mean-flow-dependent coefficients in the governing equation and the scalar Green's function in truncated Fourier series in the azimuthal direction and a finite difference approximation in the radial direction in circular cylindrical coordinates. The combined spectral/finite difference method yields a highly banded system of algebraic equations that can be efficiently solved using a standard sparse system solver. The method is applied to test cases, with mean flow specified by analytical functions, corresponding to two noise reduction concepts of current interest: the offset jet and the fluid shield. Sample results for the Green's function are given for these two test cases and recommendations made as to the use of the method as part of a RANS-based jet noise prediction code.

Leib, Stewart J.↗
DOE OSTI2022

HyKKT: a hybrid direct-iterative method for solving KKT linear systems

Here, we propose a solution strategy for the large indefinite linear systems arising in interior methods for nonlinear optimization. The method is suitable for implementation on hardware accelerators such as graphical processing units (GPUs). The current gold standard for sparse indefinite systems is the LBLT factorization where L is a lower triangular matrix and B is 1×1 or 2×2 block diagonal. However, this requires pivoting, which substantially increases communication cost and degrades performance on GPUs. Our approach solves a large indefinite system by solving multiple smaller positive definite systems, using an iterative solver on the Schur complement and an inner direct solve (via Cholesky factorization) within each iteration. Cholesky is stable without pivoting, thereby reducing communication and allowing reuse of the symbolic factorization. We demonstrate the practicality of our approach on large optimal power flow problems and show that it can efficiently utilize GPUs and outperform LBL T factorization of the full system.

97 MATHEMATICS AND COMPUTING↗
NASA NTRS1991

Iterative solution techniques in boundary element analysis

Iterative techniques for the solution of the algebraic equations associated with the direct boundary element analysis (BEA) method are discussed. Continuum structural response analysis problems are considered, employing single- and multizone boundary element models with and without zone condensation. The impact on convergence rate and computer resource requirements associated with the sparse and blocked matrices, resulting in multizone BEA, is studied. Both conjugate gradient and generalized minimum residual preconditioned iterative solvers are applied for these problems and the performance of these algorithms is reported. Included is a quantification of the impact of the preconditioning utilized to render the boundary element matrices solvable by the respective iterative methods in a time competitive with direct methods.

Kane, J. H.↗
NASA NTRS1990

High-performance equation solvers and their impact on finite element analysis

The role of equation solvers in modern structural analysis software is described. Direct and iterative equation solvers which exploit vectorization on modern high-performance computer systems are described and compared. The direct solvers are two Cholesky factorization methods. The first method utilizes a novel variable-band data storage format to achieve very high computation rates and the second method uses a sparse data storage format designed to reduce the number of operations. The iterative solvers are preconditioned conjugate gradient methods. Two different preconditioners are included; the first uses a diagonal matrix storage scheme to achieve high computation rates and the second requires a sparse data storage scheme and converges to the solution in fewer iterations that the first. The impact of using all of the equation solvers in a common structural analysis software system is demonstrated by solving several representative structural analysis problems.

Poole, Eugene L.↗
NASA NTRS1992

High-performance equation solvers and their impact on finite element analysis

The role of equation solvers in modern structural analysis software is described. Direct and iterative equation solvers which exploit vectorization on modern high-performance computer systems are described and compared. The direct solvers are two Cholesky factorization methods. The first method utilizes a novel variable-band data storage format to achieve very high computation rates and the second method uses a sparse data storage format designed to reduce the number od operations. The iterative solvers are preconditioned conjugate gradient methods. Two different preconditioners are included; the first uses a diagonal matrix storage scheme to achieve high computation rates and the second requires a sparse data storage scheme and converges to the solution in fewer iterations that the first. The impact of using all of the equation solvers in a common structural analysis software system is demonstrated by solving several representative structural analysis problems.

Poole, Eugene L.↗
NASA NTRS1994

Time integration algorithms for the two-dimensional Euler equations on unstructured meshes

Explicit and implicit time integration algorithms for the two-dimensional Euler equations on unstructured grids are presented. Both cell-centered and cell-vertex finite volume upwind schemes utilizing Roe's approximate Riemann solver are developed. For the cell-vertex scheme, a four-stage Runge-Kutta time integration, a fourstage Runge-Kutta time integration with implicit residual averaging, a point Jacobi method, a symmetric point Gauss-Seidel method and two methods utilizing preconditioned sparse matrix solvers are presented. For the cell-centered scheme, a Runge-Kutta scheme, an implicit tridiagonal relaxation scheme modeled after line Gauss-Seidel, a fully implicit lower-upper (LU) decomposition, and a hybrid scheme utilizing both Runge-Kutta and LU methods are presented. A reverse Cuthill-McKee renumbering scheme is employed for the direct solver to decrease CPU time by reducing the fill of the Jacobian matrix. A comparison of the various time integration schemes is made for both first-order and higher order accurate solutions using several mesh sizes, higher order accuracy is achieved by using multidimensional monotone linear reconstruction procedures. The results obtained for a transonic flow over a circular arc suggest that the preconditioned sparse matrix solvers perform better than the other methods as the number of elements in the mesh increases.

Slack, David C.↗
DOE OSTI2022

Two-Stage Gauss-Seidel Preconditioners and Smoothers for Krylov Solvers on a GPU Cluster: Preprint

Gauss-Seidel (GS) relaxation is often employed as a preconditioner for a Krylov solver or as a smoother for Algebraic Multigrid (AMG). However, the requisite sparse triangular solve is difficult to parallelize on many-core architectures such as graphics processing units (GPUs). In the present study, the performance of the sequential GS relaxation based on a triangular solve is compared with two-stage variants, replacing the direct triangular solve with a fixed number of inner Jacobi-Richardson (JR) iterations. When a small number of inner iterations is sufficient to maintain the Krylov convergence rate, the two-stage GS (GS2) often outperforms the sequential algorithm on many-core architectures. The GS2 algorithm is also compared with JR. When they perform the same number of ops for SpMV (e.g. three JR sweeps compared to two GS sweeps with one inner JR sweep), the GS2 iterations, and the Krylov solver preconditioned with GS2, may converge faster than the JR iterations. Moreover, for some problems (e.g. elasticity), it was found that JR may diverge with a damping factor of one, whereas two-stage GS may improve the convergence with more inner iterations. Finally, to study the performance of the two-stage smoother and preconditioner for a practical problem, these were applied to incompressible uid ow simulations on GPUs.

algebraic multigrid↗
NASA NTRS1989

Application of direct solvers to unstructured meshes for the Euler and Navier-Stokes equations using upwind schemes

The application of Newton iteration to inviscid and viscous airfoil calculations on unstructured meshes is examined. A cell-centered finite volume scheme is employed on an unstructured mesh consisting of triangles. Roe's flux difference splitting scheme is used to compute the inviscid fluxes. Higher order accuracy is achieved by an interpolation procedure that makes use of auxiliary gradients. The efficient solution of the sparse linear system of equations which arises upon linearization in time is addressed. Results are presented for inviscid and viscous test cases. The complications which arise due to the introduction of nonlinear limiters are addressed.

Venkatakrishnan, V.↗
DOE OSTI2021

Fast and Scalable Sparse Triangular Solver for Multi-GPU Based HPC Architectures

Designing efficient and scalable sparse linear algebra kernels on modern multi-GPU based HPC systems is a daunting task due to significant irregular memory references and workload imbalance across the GPUs. This is particularly the case for \textit{Sparse Triangular Solver (SpTRSV)} which introduces additional two-dimensional computation dependencies among subsequent computation steps. Dependency information is exchanged and shared among GPUs, thus warrant for efficient memory allocation, data partitioning, and workload distribution as well as fine-grained communication and synchronization support. In this work, we demonstrate that directly adopting unified memory can adversely affect the performance of SpTRSV on multi-GPU architectures, despite linking via fast interconnect like NVLinks and NVSwitches. Alternatively, we employ the latest NVSHMEM technology based on Partitioned Global Address Space programming model to enable efficient fine-grained communication and drastic synchronization overhead reduction. Furthermore, to handle workload imbalance, we propose a malleable task-pool execution model which can further enhance the utilization of GPUs. By applying these techniques, our experiments on the NVIDIA multi-GPU supernode V100-DGX-1 and DGX-2 systems demonstrate that our design can achieve on average 3.53x (up to 9.86x) speedup on a DGX-1 system and 3.66x (up to 9.64x) speedup on a DGX-2 system with 4-GPUs over the Unified-Memory design. The comprehensive sensitivity and scalability studies also show that the proposed zero-copy SpTRSV is able to fully utilize the computing and communication resources of the multi-GPU system.

Xie, Chenhao↗
DOE OSTI2024

Sparse Cholesky factorization for solving nonlinear PDEs via Gaussian processes

In recent years, there has been widespread adoption of machine learning-based approaches to automate the solving of partial differential equations (PDEs). Among these approaches, Gaussian processes (GPs) and kernel methods have garnered considerable interest due to their flexibility, robust theoretical guarantees, and close ties to traditional methods. They can transform the solving of general nonlinear PDEs into solving quadratic optimization problems with nonlinear, PDE-induced constraints. However, the complexity bottleneck lies in computing with dense kernel matrices obtained from pointwise evaluations of the covariance kernel, and its partial derivatives, a result of the PDE constraint and for which fast algorithms are scarce. The primary goal of this paper is to provide a near-linear complexity algorithm for working with such kernel matrices. We present a sparse Cholesky factorization algorithm for these matrices based on the near-sparsity of the Cholesky factor under a novel ordering of pointwise and derivative measurements. The near-sparsity is rigorously justified by directly connecting the factor to GP regression and exponential decay of basis functions in numerical homogenization. We then employ the Vecchia approximation of GPs, which is optimal in the Kullback-Leibler divergence, to compute the approximate factor. This enables us to compute ϵ-approximate inverse Cholesky factors of the kernel matrices with complexity O(N log d (N/ϵ)) in space and O(N log 2d (N/ϵ)) in time. We integrate sparse Cholesky factorizations into optimization algorithms to obtain fast solvers of the nonlinear PDE. We numerically illustrate our algorithm’s near-linear space/time complexity for a broad class of nonlinear PDEs such as the nonlinear elliptic, Burgers, and Monge-Ampère equations. In summary, we provide a fast, scalable, and accurate method for solving general PDEs with GPs and kernel methods.

97 MATHEMATICS AND COMPUTING↗
NASA NTRS1992

An incremental strategy for calculating consistent discrete CFD sensitivity derivatives

In this preliminary study involving advanced computational fluid dynamic (CFD) codes, an incremental formulation, also known as the 'delta' or 'correction' form, is presented for solving the very large sparse systems of linear equations which are associated with aerodynamic sensitivity analysis. For typical problems in 2D, a direct solution method can be applied to these linear equations which are associated with aerodynamic sensitivity analysis. For typical problems in 2D, a direct solution method can be applied to these linear equations in either the standard or the incremental form, in which case the two are equivalent. Iterative methods appear to be needed for future 3D applications; however, because direct solver methods require much more computer memory than is currently available. Iterative methods for solving these equations in the standard form result in certain difficulties, such as ill-conditioning of the coefficient matrix, which can be overcome when these equations are cast in the incremental form; these and other benefits are discussed. The methodology is successfully implemented and tested in 2D using an upwind, cell-centered, finite volume formulation applied to the thin-layer Navier-Stokes equations. Results are presented for two laminar sample problems: (1) transonic flow through a double-throat nozzle; and (2) flow over an isolated airfoil.

Korivi, Vamshi Mohan↗
DOE OSTI2021

Evidence for the 3D Radiative Effects of Boundary-Layer Clouds From Observations of Direct and Diffuse Surface Solar Fluxes

Numerical experiments have revealed the importance of horizontal transport of light in the presence of clouds (“3D effects”), with consequences for climate, weather, and solar resource availability predictions. Yet, analysis of 3D effects from observations remain sparse because of the difficulty to isolate the effect of horizontal transport in radiation measurements. In this study, we provide observational evidence for 3D effects based on the direct-diffuse partition of surface solar fluxes. It is compared to outputs from the ecRad radiative transfer scheme run on retrieved cloud profiles. Furthermore, the direct-beam calculation takes careful account of the field-of-view of the pyrheliometer to ensure consistency between observed and modeled direct fluxes. Only the solver that accounts for 3D effects is able to reproduce the observed mean direct-diffuse partition as a function of solar zenith angle and cloud cover, in particular at large solar zenith angles where cloud sides intercept most of the direct beam.

54 ENVIRONMENTAL SCIENCES↗
NASA NTRSDate not supplied

Performance Optimization Methods for a Memory-Bound, Unstructured-Grid CFD Application on Massively Parallel GPU Platforms

Computational performance of the FUN3D unstructured-grid computational fluid dynamics (CFD) application on massively parallel GPU environments is memory-bound and highly dependent upon efficient reads from and atomic updates to the irregular cell-, edge-, and node-based data structures. In this talk, we present recent efforts into optimizing select performance-critical kernels on NVIDIA Tesla V100 and A100 GPUs and AMD CDNA MI100 GPUs. A novel use of L2 cache residency controls and asynchronous loads into on-chip shared memory are explored on the A100 GPU for the sparse iterative solver, which is dominated by mixed-precision, sparse matrix vector multiplication. Demonstrations show that these methods improve global memory bandwidth utilization by 13.5% on the A100 GPU. Several techniques are also presented that use registers and/or shared memory to facilitate array transposition and aggregation which combine to reduce the frequency and increase the cache efficiency of floating-point atomic updates to the irregular data structures. These methods are demonstrated to improve the kernel throughput by nearly 500% on select kernels on the AMD MI100 over atomic updates directly to global memory. Overall, both V100 and A100 GPUs outperformed the MI100 GPU on kernels dominated by double-precision atomic updates; however, the techniques demonstrated here reduced the performance gap and improved the MI100 performance.

GPU CPU unstructured CFD memory↗
DOE OSTI2023

The Kokkos OpenMPTarget Backend: Implementation and Lessons Learned

As the supercomputing landscape diversifies, solutions such as Kokkos to write vendor agnostic applications and libraries have risen in popularity. Kokkos provides a programming model designed for performance portability, which allows developers to write a single source implementation that can run efficiently on various architectures. At its heart, Kokkos maps parallel algorithms to architecture and vendor specific backends written in lower level programming models such as CUDA and HIP. Another approach to writing vendor agnostic parallel code is using OpenMP’s directives based approach, which lets developers annotate code to express parallelism. It is implemented at the compiler level and is supported by all major high performance computing vendors, as well as the primary Open Source toolchains GNU and LLVM. Since its inception, Kokkos has used OpenMP to parallelize on CPU architectures. In this paper, we explore leveraging OpenMP for a GPU backend and discuss the challenges we encountered when mapping the Kokkos APIs and semantics to OpenMP target constructs. As an exemplar workload we chose a simple conjugate gradient solver for sparse matrices. We find that performance on NVIDIA and AMD GPUs varies widely based on details of the implementation strategy and the chosen compiler. Furthermore, the performance of the OpenMP implementations decreases with increasing complexity of the investigated algorithms.

Gayatri, Rahulkumar↗
NASA NTRS1998

Parallel Implicit Algorithms for CFD

The main goal of this project was efficient distributed parallel and workstation cluster implementations of Newton-Krylov-Schwarz (NKS) solvers for implicit Computational Fluid Dynamics (CFD.) "Newton" refers to a quadratically convergent nonlinear iteration using gradient information based on the true residual, "Krylov" to an inner linear iteration that accesses the Jacobian matrix only through highly parallelizable sparse matrix-vector products, and "Schwarz" to a domain decomposition form of preconditioning the inner Krylov iterations with primarily neighbor-only exchange of data between the processors. Prior experience has established that Newton-Krylov methods are competitive solvers in the CFD context and that Krylov-Schwarz methods port well to distributed memory computers. The combination of the techniques into Newton-Krylov-Schwarz was implemented on 2D and 3D unstructured Euler codes on the parallel testbeds that used to be at LaRC and on several other parallel computers operated by other agencies or made available by the vendors. Early implementations were made directly in Massively Parallel Integration (MPI) with parallel solvers we adapted from legacy NASA codes and enhanced for full NKS functionality. Later implementations were made in the framework of the PETSC library from Argonne National Laboratory, which now includes pseudo-transient continuation Newton-Krylov-Schwarz solver capability (as a result of demands we made upon PETSC during our early porting experiences). A secondary project pursued with funding from this contract was parallel implicit solvers in acoustics, specifically in the Helmholtz formulation. A 2D acoustic inverse problem has been solved in parallel within the PETSC framework.

Keyes, David E.↗
DOE OSTI2021

Scalable Techniques for Stochastic Power Flow Problems (Final Report)

The proposed research focuses on developing scalable algorithms for two-stage security-constrained OPF problems with AC power flow constraints, a class of problems complicated by (i) scale arising from a scenario representation; and (ii) the presence of nonlinearity, nonconvexity, and possibly second-stage discreteness or complementarity. Unfortunately, most existing solvers cannot contend with both challenges simultaneously; accordingly, the proposed research focuses on developing solution techniques that can both scale with the number of scenarios and contend with nonconvexity and second-stage complementarity. We consider three avenues for addressing such problems: (i) Variable sample-size SQP (VS-SQP) methods that combine sparse Quasi-Newton updates with a scalable variance-reduced stochastic gradient scheme for stochastic QP subproblems, allowing for contending with second-stage complementarity via regularization; (ii) Variable sample-size stochastic Interior-point (VS-sIP) schemes that propose a sampling-based regularized (to allow for contending with complementarity) interior-point schemes in which a Schur-complement technique is employed for decomposing the Newton direction computation step; (iii) Variable sample-size tractable ADMM (VS-tADMM) schemes combine variable sample-sizes with carefully designed techniques for resolving each of the nonconvex updates (by leveraging the QCQP structures). We intend to compare the three schemes using performance profiles in terms of solution quality, scalability, etc. and then select one scheme which will then be developed and further refined in Python for purposes of the GO competition.

42 ENGINEERING↗