Direct Solver Research Articles

A novel global RBF direct collocation method for solving partial differential equations with variable coefficients

Conventional radial basis function (RBF) approximation techniques face two fundamental challenges: determining the optimal shape parameter and solving the ill-conditioned linear system, which result from the use of typical kernel functions such as the Multiquadric (MQ) RBF. The objective of this study is to address these challenges by employing a novel RBF kernel, termed the coupled MQ RBF (CMQ-RBF). The new kernel combines the MQ with an m-order conical spline, where m denotes the optimal value determined by the so-called back check method (BCM) introduced in this study. Leveraging the CMQ-RBF kernel, we present a novel global RBF collocation method that generates a well-conditioned linear system, enabling direct solver utilization and delivering highly accurate numerical results without the need for laborious shape parameter selection. Consequently, the previously challenging issues can be effectively circumvented by switching to the "better kernel" offered by the CMQ-RBF. The proposed methodology is characterized by its mathematical simplicity and ease of numerical implementation. Several benchmark examples are investigated to verify its performance.

Three-Dimensional Forward Modeling of Transient Electromagnetic Method Considering Induced Polarization Effect Based on Spectral Element Method

The transient electromagnetic method (TEM) is widely used in the exploration of mineral, petroleum, and geothermal resources due to its sensitivity to low-resistivity bodies, limited site constraints, and strong resistance to interference. In practical applications, the TEM often uses a long wire source instead of an idealized horizontal electric dipole (HED) source as the excitation source. This is due to the complex external conditions and the relatively large distance between the receiving zone and the transmitter source. Compared to the HED, the long wire source can provide a larger excitation current, generating stronger signals to meet the requirements of a higher signal-to-noise ratio or deeper exploration. It also produces longer-duration signals, thereby providing better resolution. Additionally, for the interpretation of TEM data, three-dimensional forward modeling plays a crucial role. However, the mature traditional TEM forward method is based on a simple, sometimes inappropriate model, as it is well established that the induced polarization (IP) effect is widely present in the deep earth, especially in oil and gas reservoirs. The presence of the IP effect results in negative responses in field data that do not conform to the traditional theoretical decay law of TEM, which can significantly impact data processing and inversion results. To address this issue, a TEM forward modeling method considering the IP effect based on the spectral element method (SEM) has been developed in this study. Firstly, starting from the Helmholtz equation satisfied by the time domain electric field, we introduce the Debye model with polarization information into the forward modeling by utilizing the differential form of Ohm’s law. As a result, we derive the boundary value problem for the time domain electric field that considers the induced polarization effect. Using Gauss–Lobatto–Legendre (GLL) polynomials as the basis functions, the SEM is employed to discretize the governing equations at each time step and obtain spectral element discretization equations. Then, temporal discretization equations are derived using the second-order backward Euler formula, and the linear system of equations is solved using the Pardiso direct solver. Finally, the electromagnetic responses at any time channel are obtained via SEM interpolation and numerical integration, thereby achieving three-dimensional TEM forward modeling considering the IP effect. The results indicate that this method can effectively reflect the spatial distribution of polarizable subsurface media. It provides valuable references for studying the polarization parameters of subsurface media and performing a three-dimensional inversion of TEM data considering the induced polarization effect.

A Fast Algorithm for Aperiodic Linear Stencil Computation using Fast Fourier Transforms

Stencil computations are widely used to simulate the change of state of physical systems across a multidimensional grid over multiple timesteps. The state-of-the-art techniques in this area fall into three groups: cache-aware tiled looping algorithms, cache-oblivious divide-and-conquer trapezoidal algorithms, and Krylov subspace methods. In this article, we present two efficient parallel algorithms for performing linear stencil computations. Current direct solvers in this domain are computationally inefficient, and Krylov methods require manual labor and mathematical training. We solve these problems for linear stencils by using discrete Fourier transforms preconditioning on a Krylov method to achieve a direct solver that is both fast and general. Indeed, while all currently available algorithms for solving general linear stencils perform Θ( NT ) work, where N is the size of the spatial grid and T is the number of timesteps, our algorithms perform o ( NT ) work. To the best of our knowledge, we give the first algorithms that use fast Fourier transforms to compute final grid data by evolving the initial data for many timesteps at once. Our algorithms handle both periodic and aperiodic boundary conditions and achieve polynomially better performance bounds (i.e., computational complexity and parallel runtime) than all other existing solutions. Initial experimental results show that implementations of our algorithms that evolve grids of roughly 10 7 cells for around 10 5 timesteps run orders of magnitude faster than state-of-the-art implementations for periodic stencil problems, and 1.3× to 8.5× faster for aperiodic stencil problems. Code Repository: https://github.com/TEAlab/FFTStencils

A matrix solver approach for fracture flow simulation by Analytic Element Method

The manuscript presents a study on 2D groundwater flow through porous and fractured media employing Analytic Element Method. The approximate solution based on a series expansion is presented for the problem with multiple fractures. Influence of each fracture can be expressed in a series satisfying Laplace’s equation perfectly. The unknown coefficients in series expansion are determined by using the discharge potentials of all the elements linked to them. Dimension, location and conductivities were randomly selected for the inhomogeneities. Furthermore, a matrix method as a direct solver for AEM, obtained by implementing the intern boundary conditions as well as convergence analysis of iterative method of Gauss–Seidel type are also discussed.

A study of concurrent multi-frontal solvers for modern massively parallel architectures

Leveraging Trace Theory, we investigate the efficient parallelization of direct solvers for large linear equation systems. Our focus lies on a multi-frontal algorithm, and we present a methodology for achieving near-optimal scheduling on modern massively parallel machines. By employing trace theory with Diekert Graphs and Foata Normal Form, we rigorously validate the correctness of our proposed solution. To establish a strong link between the mesh and elimination tree of the multi-frontal solver, we conduct extensive testing on matrices derived from the Finite Element Method (FEM). Furthermore, we assess the performance of computations on both GPU and CPU platforms, employing practical implementation strategies.

Revising the Boundary Element Method for Thermoviscous Acoustics: An Iterative Approach via Schur Complement

The Helmholtz equation is a reliable model for acoustics in inviscid fluids. Real fluids, however, experience viscous and thermal dissipation that impact the sound propagation dynamics. The viscothermal losses primarily arise in the boundary region between the fluid and solid, the acoustic boundary layers. To preserve model accuracy for structures housing acoustic cavities of comparable size to the boundary layer thickness, meticulous consideration of these losses is essential. Recent research efforts aim to integrate viscothermal effects into acoustic boundary element methods (BEM). While the reduced discretization of BEM is advantageous over finite element methods, it results in fully populated system matrices whose conditioning deteriorates when extended with additional degrees of freedom to account for viscothermal dissipation. Solving such a linear system of equations becomes prohibitively expensive for large-scale applications, as only direct solvers can be used. This work proposes a revised formulation for the viscothermal BEM employing the Schur complement and a change of basis for the boundary coupling. We demonstrate that static condensation significantly improves the conditioning of the coupled problem. When paired with an iterative solution scheme, the approach lowers the algorithmic complexity and thus reduces the computational costs in terms of runtime and storage requirements. The results demonstrate the favorable performance of the new method, indicating its usability for applications of practical relevance in thermoviscous acoustics.

Parallel scalability of sparse solvers for sparsified fractional diffusion problems discretized with FEM

Anomalous diffusion is a non-local process that describes a wide spectrum of natural processes and phenomena with many applications in science and technology. It is described mathematically by the fractional Laplace operator. In this work we examine the integral definition of the fractional Laplacian modeled with the Riesz potential and discretized with the finite element method. The thus obtained system of linear algebraic equations is dense due to the non-local nature of the fractional Laplacian and is computationally complex to solve. With LU factorization, for example, solving the problem has an O(N 3) computational complexity, where N is the number of unknowns. However, it can be observed that a large amount of the off-diagonal coefficients have very small absolute values compared to the diagonal coefficients. Those small off-diagonal coefficients can be lumped (set to zero and added to the diagonal coefficient) without significant loss of accuracy. In this work we employ a direct sparse parallel solver to the resultant sparse matrix. We analyze the parallel performance and speed-up, as well as the accuracy, varying the fractional power and the lumping threshold.

A Generalized Source Integral Equation for Enhanced Compression in Three Dimensions

An integral equation formulation, using generalized directional sources, for three-dimensional scattering by impenetrable and essentially convex bodies, is presented. This formulation increases the rank-deficiency of moment matrix blocks representing interactions between non-overlapping source and observer regions. Thereby, it enables enhanced low-rank approximation based matrix compression and the development of corresponding fast direct solvers. The directional sources are constructed by augmenting the conventional basis functions with spherical absorbing three-dimensional shields, on which auxiliary source distributions are defined. The bottleneck arising from the need to integrate over these distributions upon computing the modified Green’s function is removed by using efficient non-uniform sampling and tabulation of the modified Green’s function or its components, in a region-dependent manner. The formulation is studied and its favorable compressibility is demonstrated, for two fundamental types of compression strategies. The non-uniform sampling approach is employed also to facilitate the rank-revealing analysis of large off-diagonal blocks.

3-D induction log modelling with integral equation method and domain decomposition pre-conditioning

SUMMARY The deployment of electromagnetic (EM) induction tools while drilling is one of the standard routines for assisting the geosteering decision-making process. The conductivity distribution obtained through the inversion of the EM induction log can provide important information about the geological structure around the borehole. To image the 3-D geological structure in the subsurface, 3-D inversion of the EM induction log is required. Because the inversion process is mainly dependent on forward modelling, the use of a fast and accurate forward modelling tool is essential. In this paper, we present an improved version of the integral equation (IE) based modelling technique for general anisotropic media with domain decomposition pre-conditioning. The discretized IE after domain decomposition equals a fixed-point equation that is solved iteratively with either the block Gauss–Seidel or Jacobi pre-conditioning. Within each iteration, the inverse of the block matrix is computed using a Krylov subspace method instead of a direct solver. An additional reduction in computational time is obtained by using an adaptive relative residual stopping criterion in the iterative solver. Using this domain decomposition scheme, numerical experiments show computation time reductions by factors of 1.97–2.84 compared to solving the full-domain IE with a GMRES solver and a contraction IE pre-conditioner. Additionally, the reduction of memory requirement for covering a large area of the induction tool sensitivity enables acceleration with limited GPU memory. Hence, we conclude that the domain decomposition method is improving the efficiency of the IE method by reducing the computation time and memory requirement.

Three-dimensional global instability analysis for high-speed boundary layer flow

Herein, a numerical approach for the three-dimensional global instability analysis of high-speed boundary layer is developed. The approach is based on the Krylov–Schur method for large-scale eigenproblems. A shift-and-invert spectral transformation is employed to accelerate the convergence for eigenvalues of interest. Iterative linear solver is adopted in the transformation to reduce memory usage by avoiding the fill-in problem of direct solver and the explicit construction of matrices. A Robin condition based on BiGlobal instability analysis is imposed at the inflow boundary, and a non-reflecting condition derived from inviscid flow is used for the far-field boundary. The approach is validated in a Blasius boundary layer. According to a comparison between the performances of several algorithms in this case, the conjugate gradient squared method with no preconditioning is chosen as the linear solver. The results provided by three-dimensional global instability analysis consist with those from linear stability theory and BiGlobal instability analysis. A further demonstration of our approach is provided for a three-dimensional roughness-disturbed boundary layer. The fluctuations extracted from direct numerical simulation results illustrate that the symmetric component is the dominant unstable mode in roughness wake, as it grows twice as fast as the antisymmetric component. Two eigenvalues are obtained from the three-dimensional global instability analysis, and these represent the symmetric and antisymmetric unstable modes. The wavelengths and growth rates of both modes delivered by instability analysis agree well with the direct numerical simulation results, and the profiles of disturbance amplitude match the features of the varicose and sinuous modes.

An explicitly-sparse representation for oscillatory kernels with wave atom-like functions

An explicitly-sparse representation for oscillatory kernels is presented in this work by developing a wave atom based method. Multilevel wave atom-like functions are constructed as the transform of the original nodal basis. Then the system matrix in a new non-standard form is derived with respect to the wave atom basis. The wave atom representation is explicitly-sparse in the sense that the system matrix is sparse and computed explicitly. Its sparsity is further enhanced via a-posteriori compression. Finally its log-linear computational complexity with controllable accuracy is demonstrated with numerical results. This explicitly-sparse representation is expected to lay ground to future work related to fast direct solvers and effective preconditioners for high frequency problems. The algorithm may also be viewed as the generalization of wavelet based methods to high frequency cases, and used as a new wideband fast algorithm for wave problems.

Algebraic Inverse Fast Multipole Method: A fast direct solver that is better than HODLR based fast direct solver

This article presents a fast direct solver, termed Algebraic Inverse Fast Multipole Method (from now on abbreviated as AIFMM), for linear systems arising out of N-body problems. AIFMM relies on the following three main ideas: (i) Certain sub-blocks in the matrix corresponding to N-body problems can be efficiently represented as low-rank matrices; (ii) The low-rank sub-blocks in the above matrix are leveraged to construct an extended sparse linear system; (iii) While solving the extended sparse linear system, certain fill-ins that arise in the elimination phase are represented as low-rank matrices and are “redirected” through other variables maintaining zero fill-in sparsity. The main highlights of this article are the following: (i) Our method is completely algebraic (as opposed to the existing Inverse Fast Multipole Method [1–3], from now on abbreviated as IFMM). We rely on our new Nested Cross Approximation [4] (from now on abbreviated as NNCA) to represent the matrix arising out of N-body problems. (ii) A significant contribution is that the algorithm presented in this article is more efficient than the existing IFMMs. In the existing IFMMs, the fill-ins are compressed and redirected as and when they are created. Whereas in this article, we update the fill-ins first without affecting the computational complexity. We then compress and redirect them only once. (iii) Another noteworthy contribution of this article is that we provide a comparison of AIFMM with Hierarchical Off-Diagonal Low-Rank (from now on abbreviated as HODLR) based fast direct solver and NNCA powered GMRES based fast iterative solver. (iv) Additionally, AIFMM is also demonstrated as a preconditioner.

All dielectric highly efficient achromatic meta-lens using inverse design optimization

This work presents a high-efficiency achromatic meta-lens based on inverse design with topology optimization methodology. The meta-lens design with high numerical aperture values (NA = 0.7, NA = 0.8, and NA = 0.9) optimized along wavelength range starts from 550 to 800 nm, then the direct solver along the full extended wavelength band from 400 to 800 nm that applied to the final optimized structures under the three conditions of the high numerical apertures have high focusing efficiency for the all conditions. The optimization problem is based on Kreisselmeier–Steinhauser (k-s) objective function, leading to approximately stable response over the broadband bandwidths of the three designs.

A fast direct solver for surface-based whole-head modeling of transcranial magnetic stimulation

When modeling transcranial magnetic stimulation (TMS) in the brain, a fast and accurate electric field solver can support interactive neuronavigation tasks as well as comprehensive biophysical modeling. We formulate, test, and disseminate a direct (i.e., non-iterative) TMS solver that can accurately determine global TMS fields for any coil type everywhere in a high-resolution MRI-based surface model with ~ 200,000 or more arbitrarily selected observation points within approximately 5 s, with the solution time itself of 3 s. The solver is based on the boundary element fast multipole method (BEM-FMM), which incorporates the latest mathematical advancement in the theory of fast multipole methods—an FMM-based LU decomposition. This decomposition is specific to the head model and needs to be computed only once per subject. Moreover, the solver offers unlimited spatial numerical resolution. Despite the fast execution times, the present direct solution is numerically accurate for the default model resolution. In contrast, the widely used brain modeling software SimNIBS employs a first-order finite element method that necessitates additional mesh refinement, resulting in increased computational cost. However, excellent agreement between the two methods is observed for various practical test cases following mesh refinement, including a biophysical modeling task. The method can be readily applied to a wide range of TMS analyses involving multiple coil positions and orientations, including image-guided neuronavigation. It can even accommodate continuous variations in coil geometry, such as flexible H-type TMS coils. The FMM-LU direct solver is freely available to academic users.

Robust fast direct integral equation solver for three-dimensional doubly periodic scattering problems with a large number of layers

Frequency-domain wave scattering problems that arise in acoustics and electromagnetism can be often described by the Helmholtz equation. This work presents a boundary integral equation method for the Helmholtz equation in 3-D multilayered media with many doubly periodic smooth layers. Compared with conventional quasi-periodic Green's function method, the new method is robust at all scattering parameters. A periodizing scheme is used to decompose the solution into near- and far-field contributions. The near-field contribution uses the free-space Green's function in an integral equation on the interface in the unit cell and its immediate eight neighbors; the far-field contribution uses proxy point sources that enclose the unit cell. A specialized high-order quadrature is developed to discretize the underlying surface integral operators to keep the number of unknowns per layer small. We achieve overall linear computational complexity in the number of layers by reducing the linear system into block tridiagonal form and then solving the system directly via block LU decomposition. The new solver is capable of handling a 100-interface structure with 961.3k unknowns to 10−5 accuracy in less than 2 hours on a desktop workstation.

Fast 3D simulation of magnetotelluric data in anisotropic media using a rational Krylov method

This study introduces a computational method aimed at accelerating multifrequency 3D magnetotelluric (MT) forward modeling in 3D conductivity structures with general anisotropy. The approach combines an edge-based finite-element (FE) method with a rational Krylov subspace method. In this method, the MT source term is expressed as a planar current source, and the frequency-dependent electric field response is the product of a transfer function and a constant vector of the current source. Consequently, the rapid approximate electric field for multifrequency calculations can be obtained by constructing the orthogonal basis of the rational Krylov subspace. By applying this technique, the large sparse matrix resulting from FE discretization can be projected onto a much lower-order matrix, significantly reducing its size to tens or hundreds. Incorporating the direct solver PARDISO to construct the orthogonal basis enables a considerably faster solution compared with the previous FE method. Moreover, the algorithm is implemented using Julia, a high-level programming language known for its readability, maintainability, and extensibility. The validity of the algorithm is demonstrated through a 1D anisotropic model, and two additional 3D models are designed to further evaluate the applicability and effectiveness of our algorithm. Numerical experiments illustrate a notable increase in modeling speed compared with the traditional edge-based FE method.

GPU-resident sparse direct linear solvers for alternating current optimal power flow analysis

Integrating renewable resources within the transmission grid at a wide scale poses significant challenges for economic dispatch as it requires analysis with more optimization parameters, constraints, and sources of uncertainty. This motivates the investigation of more efficient computational methods, especially those for solving the underlying linear systems, which typically take more than half of the overall computation time. In this paper, we present our work on sparse linear solvers that take advantage of hardware accelerators, such as graphical processing units (GPUs), and improve the overall performance when used within economic dispatch computations. We treat the problems as sparse, which allows for faster execution but also makes the implementation of numerical methods more challenging. We present the first GPU-native sparse direct solver that can execute on both AMD and NVIDIA GPUs. We demonstrate significant performance improvements when using high-performance linear solvers within alternating current optimal power flow (ACOPF) analysis. Furthermore, we demonstrate the feasibility of getting significant performance improvements by executing the entire computation on GPU-based hardware. Finally, we identify outstanding research issues and opportunities for even better utilization of heterogeneous systems, including those equipped with GPUs.

Virtual battery electric vehicle development via 1D tools

Actual testing and prototyping costs make up a significant portion of engineering budgets. Virtual demonstration mainly relies on fast and accurate models with robust performance prediction capability as a cost-effective solution. In this manuscript, GT-Suite, a one-dimensional simulation tool is preferred for developing an isothermal battery model. The developed battery model is implemented to a Light Commercial Vehicle model to run Worldwide Harmonized Light Vehicles Test Cycle. The critical outputs such as state of charge, energy depletion, heat rejection, and generated power are reported. Histograms such as pack current, charge, and discharge current with respect to its nominal capacity (C-rate) are also created to examine the operation capability of battery. In the results, it is seen that the state of charge diminishes from 90% to 82%, as an expected behavior. It is also found that the Parallel Sparse Direct and Multi-Recursive Iterative Linear Solvers methodology reduces the simulation duration by approximately 100 times, in comparison with Singular Value Decomposition Electrical Inversion Scheme. The runtime of the battery pack modeled with the cellular approach combined with the SVD method is more than 90 h. However, the runtime drops to 0.75 h when the PARDISO technique is applied. The method developed within this study can be used for rapid and accurate development of batteries. The verification can be completed in virtual environment and a vital reduction in engineering/prototype/test costs can be guaranteed. The innovation the developed methodology propose is ensuring a fast and reliable performance assessment via taking the holistic effect of integrated models into consideration. Hence it is possible for original equipment manufacturer to completely or gradually eliminate the actual tests and hence prototype costs, test system, and engineer allocations.

3D inversion of multifrequency controlled-source electromagnetic data and its application to geothermal exploration in the Tianzhen region of the northern Datong Basin, China

Owing to its sensitivity to hydrothermal alteration and geothermal fluids, the distribution of subsurface electrical conductivity provides critical information for characterizing geothermal systems. For geothermal exploration, the controlled-source electromagnetic (CSEM) technique provides a crucial tool to image the subsurface resistivity structures, especially in areas with strong cultural noise. Usually, land-based CSEM surveys are carried out with dozens of operating frequencies to enhance spatial resolution. However, the 3D inversion of multifrequency CSEM data using the commonly used direct solver is challenging with limited computational resources. In this study, we implement a practical inversion strategy for interpreting 3D multifrequency CSEM data. By combining a hybrid direct-iterative solver with the inexact Gauss-Newton optimization, 3D inversion of CSEM data involving multiple frequencies can be performed effectively on typical workstations. First, we test the effectiveness of the developed approach using synthetic multifrequency 3D CSEM data generated for a simplified geothermal model. The inversion strategy is then applied to the multifrequency CSEM field data collected for the geothermal exploration at Tianzhen region in Shanxi Province, China. The comparison between the inversion results using the sparse data set (six frequencies) and the dense data set (12 frequencies) highlights the necessity of using data from more frequencies in the inversion to improve the resolution. The resulting 3D resistivity model clearly delineates the clay alteration layers and shallow thermal reservoirs of the potential high-temperature geothermal system within the study area, while its deep heat source is not revealed owing to the limited investigation depth of the survey.

A multigrid assisted reanalysis method for accelerated heat transfer topology optimization

Heat sinks are widely applied in industry, and heat transfer topology optimization can improve their performance while saving materials. However, the prohibitive computational cost hinders its development. Therefore, the Multigrid Assisted Reanalysis (MGAR) method is proposed to address the issue. The Combined Approximate Reduced Order Model (CAROM) is employed to reduce the dimensionality of governing equations, and computational efficiency is significantly improved. In addition, CAROM should be dynamically updated to ensure the accuracy of solutions, thus constructing the CAROM also take up a major computational cost. For this purpose, the V-cycle multigrid is employed to construct basis vectors of CAROM instead of direct solvers. Since the MGAR involves almost no matrix decomposition, except for solving the reduced linear systems, the MGAR is equipped to handle large-scale problems due to the advantage of high computational efficiency and low memory requirements. The competence of the MGAR is verified by several cases including thermal compliance minimization and maximum temperature minimization.