Fast Solver Research Articles

A comparison of the stability and performance of depth-integrated ice-dynamics solvers

Abstract. In the last decade, the number of ice-sheet models has increased substantially, in line with the growth of the glaciological community. These models use solvers based on different approximations of ice dynamics. In particular, several depth-integrated dynamics solvers have emerged as fast solvers capable of resolving the relevant physics of ice sheets at the continental scale. However, the numerical stability of these schemes has not been studied systematically to evaluate their effectiveness in practice. Here we focus on three such solvers, the so-called Hybrid, L1L2-SIA and DIVA solvers, as well as the well-known SIA and SSA solvers as boundary cases. We investigate the numerical stability of these solvers as a function of grid resolution and the state of the ice sheet for an explicit time discretization scheme of the mass conservation step. Under simplified conditions with constant viscosity, the maximum stable time step of the Hybrid solver, like the SIA solver, has a quadratic dependence on grid resolution. In contrast, the DIVA solver has a maximum time step that is independent of resolution as the grid becomes increasingly refined, like the SSA solver. A simple 1D implementation of the L1L2-SIA solver indicates that it should behave similarly, but in practice, the complexity of its implementation appears to restrict its stability. In realistic simulations of the Greenland Ice Sheet with a nonlinear rheology, the DIVA and SSA solvers maintain superior numerical stability, while the SIA, Hybrid and L1L2-SIA solvers show markedly poorer performance. At a grid resolution of Δx=4 km, the DIVA solver runs approximately 20 times faster than the Hybrid and L1L2-SIA solvers as a result of a larger stable time step. Our analysis shows that as resolution increases, the ice-dynamics solver can act as a bottleneck to model performance. The DIVA solver emerges as a clear outlier in terms of both model performance and its representation of the ice-flow physics itself.

Preconditioned Legendre Spectral Galerkin Methods for the Non-separable Elliptic Equation

The Legendre spectral Galerkin method of self-adjoint second order elliptic equations usually results in a linear system with a dense and ill-conditioned coefficient matrix. In this paper, the linear system is solved by a preconditioned conjugate gradient (PCG) method where the preconditioner M is constructed by approximating the variable coefficients with a (T+1)-term Legendre series in each direction to a desired accuracy. A feature of the proposed PCG method is that the iteration step increases slightly with the size of the resulting matrix when reaching a certain approximation accuracy. The efficiency of the method lies in that the system with the preconditioner M is approximately solved by a one-step method based on the incomplete LU factorization technique with no fill-in, denoted by ILU(0). The ILU(0) factorization of \(M\in {\mathbb {R}}^{(N-1)^d\times (N-1)^d}\) can be computed using \({\mathcal {O}}(T^{2d} N^d)\) operations, and the number of nonzeros in the factorization factors is of \({\mathcal {O}}(T^{d} N^d)\), \(d=1,2,3\). A conclusion of the algorithm is to fast solve the resulting system from the Legendre Galerkin spectral method for Poisson equations with Dirichlet boundary conditions, which has a complexity of \({\mathcal {O}}(N^d)\). To further speed up the PCG method, an algorithm is developed for fast matrix-vector multiplications by the resulting matrix of Legendre-Galerkin spectral discretization, without the need to explicitly form it. The complexity of the fast matrix-vector multiplications is of \({\mathcal {O}}(N^d (\log _2 N)^{2})\). In view that T is independent of N in one dimension and is set to be of order \({\mathcal {O}}(\log _2 N)\) in two and three dimensions, the PCG method has a \({\mathcal {O}}(N^d (\log _2 N)^{2d})\) quasi-optimal complexity for a d dimensional domain with \((N-1)^d\) unknows, \(d=1,2,3\). In addition, a fast direct solver for the three-dimensional Poisson equation is developed, which is of \({\mathcal {O}}(N^{3} (\log _2 N)^{2})\) and improves the existing results on the computational complexity. Numerical examples are given to demonstrate the efficiency of proposed preconditioners and the algorithm for fast matrix-vector multiplications.

An unconditionally stable fast high order method for thermal phase change models

Thermal phase change problems arise in a large number of applications. In this paper, we consider a phase field model instead of the classical Stefan model to describe phenomena, which may appear in some complex phase change problems such as dendritic crystal growth, phase transformations in metallic alloys, etc. Our aim is to propose efficient and accurate schemes for the model, which is the coupling of a heat transfer equation and a phase field equation. The schemes are constructed based on an auxiliary variable approach for the phase field equation and semi-implicit treatment for the heat transfer equation. The main novelty of the paper consists in: (i) construction of the efficient schemes, which only requires solving several second-order elliptic problems with constant coefficients; (ii) proof of the unconditional stability of the schemes; (iii) fast high order solver for the resulting equations at each time step. A series of numerical examples are presented to verify the theoretical claims and to illustrate the efficiency of our method. As far as we know, it seems this is the first attempt made for the thermal phase change model of this type.

A New Directional Algebraic Fast Multipole Method Based Iterative Solver for the Lippmann-Schwinger Equation Accelerated with HODLR Preconditioner

We present a fast iterative solver for scattering problems in 2D, where a penetrable object with compact support is considered. By representing the scattered field as a volume potential in terms of the Green's function, we arrive at the Lippmann-Schwinger equation in integral form, which is then discretized using an appropriate quadrature technique. The discretized linear system is then solved using an iterative solver accelerated by Directional Algebraic Fast Multipole Method (DAFMM). The DAFMM presented here relies on the directional admissibility condition of the 2D Helmholtz kernel. And the construction of low-rank factorizations of the appropriate low-rank matrix sub-blocks is based on our new Nested Cross Approximation (NCA)~\cite{ arXiv:2203.14832 [math.NA]}. The advantage of our new NCA is that the search space of so-called far-field pivots is smaller than that of the existing NCAs. Another significant contribution of this work is the use of HODLR based direct solver as a preconditioner to further accelerate the iterative solver. In one of our numerical experiments, the iterative solver does not converge without a preconditioner. We show that the HODLR preconditioner is capable of solving problems that the iterative solver can not. Another noteworthy contribution of this article is that we perform a comparative study of the HODLR based fast direct solver, DAFMM based fast iterative solver, and HODLR preconditioned DAFMM based fast iterative solver for the discretized Lippmann-Schwinger problem. To the best of our knowledge, this work is one of the first to provide a systematic study and comparison of these different solvers for various problem sizes and contrast functions. In the spirit of reproducible computational science, the implementation of the algorithms developed in this article is made available at \url{https://github.com/vaishna77/Lippmann_Schwinger_Solver}.

Fast Solver for Implicit Continuous Set Model Predictive Control of Electric Drives

This paper proposes a fast and accurate solver for implicit Continuous Set Model Predictive Control for the current control loop of synchronous motor drives with input constraints, allowing for reaching the maximum voltage feasible set. The related control problem requires an iterative solver to find the optimal solution. The real-time certification of the algorithm is of paramount importance to move the technology toward industrial-scale applications. A relevant feature of the proposed solver is that the total number of operations can be computed in the worst-case scenario. Thus, the maximum computational time is known a priori. The solver is deeply illustrated, showing its feasibility for real-time applications in the microseconds range by means of experimental tests. The proposed method outperforms general-purpose algorithms in terms of computation time, while keeping the same accuracy.

Low-dose Cone-Beam Computed Tomography Reconstruction through a fast Three-Dimensional Compressed Sensing Method Based on the Three-Dimensional Pseudo-polar Fourier Transform.

Background:Reconstruction of high quality two dimensional images from fan beam computed tomography (CT) with a limited number of projections is already feasible through Fourier based iterative reconstruction method. However, this article is focused on a more complicated reconstruction of three dimensional (3D) images in a sparse view cone beam computed tomography (CBCT) by utilizing Compressive Sensing (CS) based on 3D pseudo polar Fourier transform (PPFT).Method:In comparison with the prevalent Cartesian grid, PPFT re gridding is potent to remove rebinning and interpolation errors. Furthermore, using PPFT based radon transform as the measurement matrix, reduced the computational complexity.Results:In order to show the computational efficiency of the proposed method, we compare it with an algebraic reconstruction technique and a CS type algorithm. We observed convergence in <20 iterations in our algorithm while others would need at least 50 iterations for reconstructing a qualified phantom image. Furthermore, using a fast composite splitting algorithm solver in each iteration makes it a fast CBCT reconstruction algorithm. The algorithm will minimize a linear combination of three terms corresponding to a least square data fitting, Hessian (HS) Penalty and l1 norm wavelet regularization. We named it PP-based compressed sensing-HS-W. In the reconstruction range of 120 projections around the 360° rotation, the image quality is visually similar to reconstructed images by Feldkamp-Davis-Kress algorithm using 720 projections. This represents a high dose reduction.Conclusion:The main achievements of this work are to reduce the radiation dose without degrading the image quality. Its ability in removing the staircase effect, preserving edges and regions with smooth intensity transition, and producing high-resolution, low-noise reconstruction results in low-dose level are also shown.

Experimental and Numerical Investigation of Contact Parameters in a Dovetail Type of Blade Root Joints

This paper focuses on the contact characteristics of the blade root joints subjected to the dry friction damping under periodic excitation. The numerical method and experimental procedure are combined to trace the contact behavior in the nonlinear vibration conditions. In experimental procedure, a novel excitation method alongside the accurate measurements is used to determine the frequencies of the blade under different axial loads. In numerical simulations, local behavior of contact areas is investigated using the reduction method as a reliable and fast solver. Subsequently, by using both experimental measurements and numerical outcomes in a developed code, the global stiffness matrix is calculated. This leads to find the normal and tangential stiffness in the contact areas of a dovetail blade root joints. The results indicate that the proposed method can provide an accurate quantitative assessment for investigation the dynamic response of the joints with focusing the contact areas.

Fast Solver for J2-Perturbed Lambert Problem Using Deep Neural Network

This paper presents a novel and fast solver for the J2-perturbed Lambert problem. The solver consists of an intelligent initial guess generator combined with a differential correction procedure. The intelligent initial guess generator is a deep neural network that is trained to correct the initial velocity vector coming from the solution of the unperturbed Lambert problem. The differential correction module takes the initial guess and uses a forward shooting procedure to further update the initial velocity and exactly meet the terminal conditions. Eight sample forms are analyzed and compared to find the optimum form to train the neural network on the J2-perturbed Lambert problem. The accuracy and performance of this novel approach will be demonstrated on a representative test case: the solution of a multirevolution J2-perturbed Lambert problem in the Jupiter system. The performance of the proposed approach will be compared against a classical standard shooting method and a homotopy-based perturbed Lambert algorithm. It will be shown that, for a comparable level of accuracy, the proposed method is significantly faster than the other two.

Fast iterative solver for the optimal control of time‐dependent PDEs with Crank–Nicolson discretization in time

AbstractIn this article, we derive a new, fast, and robust preconditioned iterative solution strategy for the all‐at‐once solution of optimal control problems with time‐dependent PDEs as constraints, including the heat equation and the non‐steady convection–diffusion equation. After applying an optimize‐then‐discretize approach, one is faced with continuous first‐order optimality conditions consisting of a coupled system of PDEs. As opposed to most work in preconditioning the resulting discretized systems, where a (first‐order accurate) backward Euler method is used for the discretization of the time derivative, we employ a (second‐order accurate) Crank–Nicolson method in time. We apply a carefully tailored invertible transformation for symmetrizing the matrix, and then derive an optimal preconditioner for the saddle‐point system obtained. The key components of this preconditioner are an accurate mass matrix approximation, a good approximation of the Schur complement, and an appropriate multigrid process to apply this latter approximation—these are constructed using our work in transforming the matrix system. We prove the optimality of the approximation of the Schur complement through bounds on the eigenvalues, and test our solver against a widely‐used preconditioner for the linear system arising from a backward Euler discretization. These demonstrate the effectiveness and robustness of our solver with respect to mesh‐sizes, regularization parameter, and diffusion coefficient.

Fast Power Series Solution of Large 3-D Electrodynamic Integral Equation for PEC Scatterers

This paper presents a new fast power series solution method to solve the Hierarchal Method of Moment (MoM) matrix for a large complex, perfectly electric conducting (PEC) 3D structures. The proposed power series solution converges in just two (2) iterations which is faster than the conventional fast solver–based iterative solution. The method is purely algebraic in nature and, as such applicable to existing conventional methods. The method uses regular fast solver Hierarchal Matrix (H-Matrix) and can also be applied to Multilevel Fast Multipole Method Algorithm (MLFMA). In the proposed method, we use the scaling of the symmetric near-field matrix to develop a diagonally dominant overall matrix to enable a power series solution. Left and right block scaling coefficients are required for scaling near-field blocks to diagonal blocks using Schur’s complement method. However, only the right-hand scaling coefficients are computed for symmetric near-field matrix leading to saving of computation time and memory. Due to symmetric property, the left side-block scaling coefficients are just the transpose of the right-scaling blocks. Next, the near-field blocks are replaced by scaled near-field diagonal blocks. Now the scaled near-field blocks in combination with far-field and scaling coefficients are subjected to power series solution terminating after only two terms. As all the operations are performed on the near-field blocks, the complexity of scaling coefficient computation is retained as O(N)O⁢(N). The power series solution only involves the matrix-vector product of the far-field, scaling coefficients blocks, and inverse of scaled near-field blocks. Hence, the solution cost remains O(NlogN)O⁢(N⁢l⁢o⁢g⁢N). Several numerical results are presented to validate the efficiency and robustness of the proposed numerical method.

A matrix-theoretic spectral analysis of incompressible Navier–Stokes staggered DG approximations and a related spectrally based preconditioning approach

The incompressible Navier-Stokes equations are solved in a channel, using a Discontinuous Galerkin method over staggered grids. The resulting linear systems are studied both in terms of the structure and in terms of the spectral features of the related coefficient matrices. In fact, the resulting matrices are of block type, each block showing Toeplitz-like, band, and tensor structure at the same time. Using this rich matrix-theoretic information and the Toeplitz, Generalized Locally Toeplitz technology, a quite complete spectral analysis is presented, with the target of designing and analyzing fast iterative solvers for the associated large linear systems. Quite promising numerical results are presented, commented, and critically discussed for elongated two- and three-dimensional geometries.

A second-order numerical method for Landau-Lifshitz-Gilbert equation with large damping parameters

A second order accurate numerical scheme is proposed and implemented for the Landau-Lifshitz-Gilbert equation, which models magnetization dynamics in ferromagnetic materials, with large damping parameters. The main advantages of this method are associated with the following features: (1) It only solves linear systems of equations with coefficient matrices independent of the magnetization, and fast solvers are available, so that the numerical efficiency has been greatly improved, in comparison with the existing Gauss-Seidel project method. (2) The second-order accuracy in time is achieved, and it is unconditionally stable for large damping parameters. Moreover, both the second-order accuracy and the great efficiency improvement will be verified by several numerical examples in the 1D and 3D simulations. In the presence of large damping parameters, it is observed that this method is unconditionally stable and finds physically reasonable structures while many existing methods have failed. For the domain wall dynamics, the linear dependence of wall velocity with respect to the damping parameter and the external magnetic field will be obtained through the reported simulations.

ALGAMES: a fast augmented Lagrangian solver for constrained dynamic games

Dynamic games are an effective paradigm for dealing with the control of multiple interacting actors. This paper introduces augmented Lagrangian GAME-theoretic solver (ALGAMES), a solver that handles trajectory-optimization problems with multiple actors and general nonlinear state and input constraints. Its novelty resides in satisfying the first-order optimality conditions with a quasi-Newton root-finding algorithm and rigorously enforcing constraints using an augmented Lagrangian method. We evaluate our solver in the context of autonomous driving on scenarios with a strong level of interactions between the vehicles. We assess the robustness of the solver using Monte Carlo simulations. It is able to reliably solve complex problems like ramp merging with three vehicles three times faster than a state-of-the-art DDP-based approach. A model-predictive control (MPC) implementation of the algorithm, running at more than 60 Hz, demonstrates ALGAMES’ ability to mitigate the “frozen robot” problem on complex autonomous driving scenarios like merging onto a crowded highway.

A variational framework for computing Wannier functions using dictionary learning

This article introduces a data-driven variational framework for computing Wannier functions by using basis pursuit to transfer information from a general feature dictionary. General features are learned by applying dictionary learning to a dataset of Wannier functions. Our approach displays a systematically controllable energy-localization trade-off, an objective functional allowing for the use of fast numerical solvers, and is capable of being used in highly-efficient self-consistent algorithms.

Hybrid Method of ACA and Interpolation-based Algorithms for Analysis of Eddy Current Nondestructive Evaluations

ABSTRACT In this article, the hybridization of adaptive cross approximation (ACA) algorithm and interpolation-based separation of the kernel function is proposed to accelerate solving the matrix equations resulted in the boundary element method (BEM) for 3D arbitrary-shaped eddy current nondestructive evaluation problems. The hybrid method combines the advantages of both the ACA algorithm and the interpolation-based methods, and resolves the shortcoming of pure ACA method, when modeling the planar eddy current nondestructive evaluation problems, that it cannot compress the null entries the BEM generated when the testing and basis patches are co-planar. In the proposed method, the submatrices associated with the null entries are compressed by the interpolation-based method, while the others are compressed by the ACA algorithm. Several benchmarks are shown to demonstrate both the robustness and efficiency of the proposed fast and general solver.

Novel accurate modeling of dust loaded wire-duct precipitators using FDM-FMG method on one fine computational domains

Worldwide, electrostatic precipitators (ESP) have been extensively utilized to separate fine particles for diverse large-scale industrial applications. In this regard, this paper presents a novel approach for modeling the dust-loaded ESP on the fine computational domain where the need for a fast solver arises. Unlike the previously published numerical techniques, the finite difference method (FDM) integrated with a full multi-grid method (FMG), labeled FDM-FMG, is developed to resolve Poisson and continuity equations on one fine computational domain. For clean and dust-loaded ESP, the proposed FMG is checked versus successive over-relaxation (SOR) on fine domains where the proposed one is greatly transcendent in terms of convergence characteristics and hence the computational performance (CPU time). For the first time, two major issues are highlighted and solved: the first concerning issue is the chosen ion mobility as an important factor in the simulation results and the second one is choosing an optimal computational grid for dust loaded precipitators that grantees both low truncation and roundoff errors, results in well-matched with experimental measurements nominated in the previous publishing. The novel idea of working on various grid sizes and tracking the optimal ones gives the FDM-FMG an advantage of predicting a precise picture for the electrical situations in industrial ESP over the other numerical techniques. After all, the impact of changing the spacing between the different wires and the height of the ionized wires on the distributions of current, ion, and particle charge densities on the ground are deeply simulated and presented in dust-loaded ESP. The proposed FDM-FMG can be a promising tool for the designers and manufacturers of precipitators, thanks to its superior computational performance.

A new FV scheme and fast cell-centered multigrid solver for 3D anisotropic diffusion equations with discontinuous coefficients

In this paper, an efficient cell-centered extrapolation cascadic multigrid (CEXCMG) method is proposed for solving large linear system of equations resulting from finite volume (FV) discretizations of three dimensional (3D) anisotropic diffusion equations with discontinuous coefficients. For cell-centered FV schemes, the values at vertex need to be approximated often by a linear combination of neighboring cell-centered values. In the literature, the weighted coefficients are obtained by solving local linear system of equations which is costly in 3D. One of the novelties of this paper is a new approach for obtaining vertex values by interpolating the cell-centered ones, which avoids solving local linear system of equations even with arbitrary diffusion tensors. Another main novelty of this paper is a new cascadic multigrid solver based on a prolongation operator, the newly developed explicit gradient transfer method, and a splitting extrapolation operator for solving 3D anisotropic diffusion equations with discontinuous coefficients. Numerical experiments are presented to demonstrate the efficiency and robustness of the CEXCMG method in terms of the mesh size and the contrast in the coefficients of the anisotropic diffusion tensor.

Fast image inpainting strategy based on the space-fractional modified Cahn-Hilliard equations

The solution strategy of the space-fractional modified Cahn-Hilliard equation as a tool for the gray value image inpainting model is studied. The existing strategies solve the convexity splitting scheme of the vector-valued Cahn-Hilliard model by Fourier spectral method. In this paper, we constructed a fast solver for the discretized linear systems possessing the saddle-point structure within block-Toeplitz-Toeplitz-block (BTTB) structure arising from the 2D space-fractional modified Cahn-Hilliard equation. The new solver enjoys computational advantage since circulant approximation and fast Fourier transforms (FFTs) can be used for solving the involved linear subsystems. Theoretical analysis shows the spectrum of the preconditioned matrix clusters around 1, which implies the fast convergence rate of the proposed preconditioner. Numerical examples are given to confirm the effectiveness of our method.

Multiscale modeling and simulation methods for electromagnetic and multiphysics problems

Multiscale modeling and simulation methods for electromagnetic and multiphysics problems

Discretization error estimates for discontinuous Galerkin isogeometric analysis

Isogeometric analysis is a spline-based discretization method to partial differential equations which show the approximation power of a high-order method. The number of degrees of freedom, however, is as small as the number of degrees of freedom of a low-order method. This does not come for free as the original formulation of isogeometric analysis requires a global geometry function. Since this is too restrictive for many kinds of applications, the domain is usually decomposed into patches, where each patch is parameterized with its own geometry function. In simpler cases, the patches can be combined in a conforming way. However, for non-matching discretizations or for varying coefficients, a non-conforming discretization is desired. An symmetric interior penalty discontinuous Galerkin method for isogeometric analysis has been previously introduced. In the present paper, we give error estimates that are explicit in the spline degree. This opens the door towards the construction and the analysis of fast linear solvers, particularly multigrid solvers for non-conforming multipatch isogeometric analysis.