Fast Iterative Solvers Research Articles

An improved implicit direct-forcing immersed boundary method (DF-IBM) around arbitrarily moving rigid structures

An improved implicit direct-forcing immersed boundary method (DF-IBM) is presented for simulating incompressible flows around complex rigid structures undergoing arbitrary motion. The current approach harnesses the pressure implicit with splitting of operators algorithm to handle the fluid–solid system's dual constraints in a segregated manner. As a result, the divergence-free condition is preserved throughout the Eulerian domain, and the no-slip velocity boundary condition is exactly enforced on the immersed boundary. A new pressure Poisson equation (PPE) is derived, incorporating the boundary force where the no-slip condition is already met, enabling the use of fast iterative PPE solvers without modifications. The improvement involves integrating Lagrangian weight methods having better reciprocity over the IBM-related linear operators with the implicit formulation. An additional force initialization scheme is introduced to further boost the algorithm's performance. The method's accuracy, efficiency, and capability are verified through various stationary and moving immersed boundary cases. The results are validated against experimental and numerical data from the literature. The proposed improvements seamlessly integrate into existing incompressible fluid solvers with minimal adjustments to the original system equations, highlighting their ease of implementation.

A Fast Multi-Frequency Iterative Acoustic Boundary Element Method Solver Based on a Preconditioned Accelerated Recycling Strategy

Common multi-frequency solution strategies for acoustic systems are either associated with high computational costs or resort to model reduction strategies, which potentially induce considerable approximation errors. In that context, this work proposes a fast multi-frequency solver for acoustic Boundary Element Method (BEM) analyses that combines conventional preconditioners with an accelerated recycling strategy, achieving in that way a fast and accurate iterative solution. Specifically, by approximating both the system and the related preconditioner by individual affine expressions, an upfront Galerkin projection on a global deflation basis is facilitated, therefore enabling the deployment of the accelerated recycling scheme. The latter acts by accelerating the construction of the related deflation projectors, therefore leading to a significantly faster recycling procedure. Finally, the efficiency of the global deflation basis is guaranteed through an Automatic Krylov subspaces Recycling for Deflation (AKR-D) [D. Panagiotopoulos, W. Desmet and E. Deckers, An accelerated subspaces recycling strategy for the deflation of parametric linear systems based on model order reduction, Comput. Methods Appl. Mech. Eng. 403 (2023) 115765] algorithm and thus, the cumulative number of iterations required by an iterative solver within a frequency sweep is drastically decreased. The proposed multi-frequency solver is tested on two industrially relevant examples of an interior and an exterior scattering problem, benchmarked against the traditional solution strategies.

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.

An adaptive solution strategy for Richards' equation

Flow in variably saturated porous media is typically modeled by the Richards equation, a nonlinear elliptic-parabolic equation which is notoriously challenging to solve numerically. In this paper, we propose a robust and fast iterative solver for Richards' equation. The solver relies on an adaptive switching algorithm, based on rigorously derived a posteriori indicators, between two linearization methods: L-scheme and Newton. Although a combined L-scheme/Newton strategy was introduced previously in [1], here, for the first time we propose a reliable and robust criteria for switching between these schemes. The performance of the solver, which can be in principle applied to any spatial discretization and linearization methods, is illustrated through several numerical examples.

Fast iterative solver for 3D finite-element modeling of controlled-source electromagnetic responses with vector and scalar potentials

The controlled-source electromagnetic (CSEM) forward modeling with vector and scalar potentials could provide insights into the inductive and galvanic effects. We have developed a finite-element (FE) forward modeling algorithm for the CSEM with vector and scalar potentials. The first-order vector and nodal shape functions are used for the vector and scalar potentials, respectively. To mitigate the null space and thus the nonuniqueness of the potentials caused by the curl operator, we introduce a new preconditioner constructed from an incomplete Cholesky decomposition with zero fill-ins of the Laplacian approximation to the resulted matrix after the FE discretization. We implemented a preconditioned quasiminimal residual method to iteratively solve the resulting linear system with the new preconditioner. We first verified the accuracy of the algorithm through comparison against the analytic solutions obtained for a three-layered earth model. The new iterative modeling algorithm converges much faster compared with the modeling algorithm with the preconditioner constructed from the original stiffness matrix. The efficiency and accuracy are further illustrated by comparison with the modeling scheme based on coupled potentials with explicit enforcement of the Coulomb gauge condition for the vector potential by modeling three 3D models. In addition, the accuracy of our algorithm is demonstrated via a comparison of our numerical solution for a 3D model with the numerical solution obtained from the integral equation method. We also present the inductive and galvanic components of the horizontal secondary electric field for numerical solutions obtained with direct and iterative solvers. The numerical test demonstrated that, aside from the efficiency improved with the new iterative solver, the unstable and nonuniqueness problem for the potentials is eliminated by the new preconditioner with the Laplacian approximation involved.

Fast Iterative Integral Equation Solver for Acoustic Scattering by Inhomogeneous Objects Using the Butterfly Approximation.

The acoustic scattering by highly inhomogeneous objects is analyzed by a method-of-moment solver for the volume integral equation. To enable the treatment of acoustically large scatterers of various topologies, the iterative numerical solution of the resulting system is accelerated via a kernel independent algebraic compression scheme: blocks of the hierarchically partitioned moment stiffness matrix are expressed in butterfly form that, for volume problems, scales favorably compared to the popular low-rank approximation. A detailed description of the algorithm, as implemented in this work, is provided. Validations of the numerical formulation, parameter tuning, and performance study of the fast method for acoustically large objects are presented, in various settings and for a range of examples, representative of biomedical and oceanographic applications.

Fast nonlinear iterative solver for an implicit, energy-conserving, asymptotic-preserving charged-particle orbit integrator

• We developed a new reformulation for energy-conserving AP particle orbit integrator. • New reformulation leads to a simple and efficient Picard iterative solver. • New reformulation shows better performance than an optimized JFNK solver.

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 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.

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.

Solution of a Linear Nondegenerate Matrix Equation Based on the Zero Divisor

New formulas were obtained to solve the linear non-degenerate matrix equations based on zero divisors of numerical matrices. Two theorems were formulated, and a proof to one of them is provided. It is noted that the proof of the second theorem is similar to the proof of the first one. The proved theorem substantiates new formula in solving the equation equivalent in the sense of the solution uniqueness to the known formulas. Its fundamental difference lies in the following: any explicit matrix inversion or determinant calculation is missing; solution is "based" not on the left, but on the right side of the matrix equation; zero divisor method is used (it was never used in classical formulas for solving a matrix equation); zero divisor calculation is reduced to simple operations of permutating the vector elements on the right-hand side of the matrix equation. Examples are provided of applying the proposed method for solving a nondegenerate matrix equation to the numerical matrix equations. High accuracy of the proposed formulas for solving the matrix equations is demonstrated in comparison with standard solvers used in the MATLAB environment. Similar problems arise in the synthesis of fast and ultrafast iterative solvers of linear matrix equations, as well as in nonparametric identification of abnormal (emergency) modes in complex technical systems, for example, in the power systems

A Helmholtz equation solver using unsupervised learning: Application to transcranial ultrasound

Transcranial ultrasound therapy is increasingly used for the non-invasive treatment of brain disorders. However, conventional numerical wave solvers are currently too computationally expensive to be used online during treatments to predict the acoustic field passing through the skull (e.g., to account for subject-specific dose and targeting variations). As a step towards real-time predictions, in the current work, a fast iterative solver for the heterogeneous Helmholtz equation in 2D is developed using a fully-learned optimizer. The lightweight network architecture is based on a modified UNet that includes a learned hidden state. The network is trained using a physics-based loss function and a set of idealized sound speed distributions with fully unsupervised training (no knowledge of the true solution is required). The learned optimizer shows excellent performance on the test set, and is capable of generalization well outside the training examples, including to much larger computational domains, and more complex source and sound speed distributions, for example, those derived from x-ray computed tomography images of the skull.

A multiscale model of terrain dynamics for real-time earthmoving simulation

A multiscale model for real-time simulation of terrain dynamics is explored. To represent the dynamics on different scales the model combines the description of soil as a continuous solid, as distinct particles and as rigid multibodies. The models are dynamically coupled to each other and to the earthmoving equipment. Agitated soil is represented by a hybrid of contacting particles and continuum solid, with the moving equipment and resting soil as geometric boundaries. Each zone of active soil is aggregated into distinct bodies, with the proper mass, momentum and frictional-cohesive properties, which constrain the equipment’s multibody dynamics. The particle model parameters are pre-calibrated to the bulk mechanical parameters for a wide range of different soils. The result is a computationally efficient model for earthmoving operations that resolve the motion of the soil, using a fast iterative solver, and provide realistic forces and dynamic for the equipment, using a direct solver for high numerical precision. Numerical simulations of excavation and bulldozing operations are performed to test the model and measure the computational performance. Reference data is produced using coupled discrete element and multibody dynamics simulations at relatively high resolution. The digging resistance and soil displacements with the real-time multiscale model agree with the reference model up to 10–25%, and run more than three orders of magnitude faster.

Fast iterative solvers for the two-dimensional spatial fractional Ginzburg–Landau equations

In this work, we propose an alternating direction implicit (ADI) scheme to discrete the two-dimensional spatial fractional Ginzburg–Landau equations. Meanwhile, a matrix splitting iteration method that preserves the Toeplitz structure is proposed for solving the resulting complex linear system, in which the circulant preconditioning technique and fast Fourier transform (FFT) can be adopted to improve computing efficiency. Convergence properties of the corresponding method are derived under some conditions. Numerical experiments demonstrate that our method outperforms some existing iteration methods in iterative steps and computing time.

A Fast Algorithm for the Variable-Order Spatial Fractional Advection-Diffusion Equation

We propose a fast algorithm for the variable-order (VO) space-fractional advection-diffusion equations with nonlinear source terms on a finite domain. Due to the impact of the space-dependent the VO, the resulting coefficient matrices arising from the finite difference discretization of the fractional advection-diffusion equation are dense without Toeplitz-like structure. By the properties of the elements of coefficient matrices, we show that the off-diagonal blocks can be approximated by low-rank matrices. Then we present a fast algorithm based on the polynomial interpolation to approximate the coefficient matrices. The approximation can be constructed in $${\mathcal {O}}(kN)$$ operations and requires $${\mathcal {O}}(kN)$$ storage with N and k being the number of unknowns and the approximants, respectively. Moreover, the matrix-vector multiplication can be implemented in $${\mathcal {O}} (kN\log N)$$ complexity, which leads to a fast iterative solver for the resulting linear systems. The stability and convergence of the new scheme are also studied. Numerical tests are carried out to exemplify the accuracy and efficiency of the proposed method.

A Segregated Approach for Modeling the Electrochemistry in the 3-D Microstructure of Li-Ion Batteries and Its Acceleration Using Block Preconditioners

Battery performance is strongly correlated with electrode microstructure. Electrode materials for lithium-ion batteries have complex microstructure geometries that require millions of degrees of freedom to solve the electrochemical system at the microstructure scale. A fast-iterative solver with an appropriate preconditioner is then required to simulate large representative volume in a reasonable time. In this work, a finite element electrochemical model is developed to resolve the concentration and potential within the electrode active materials and the electrolyte domains at the microstructure scale, with an emphasis on numerical stability and scaling performances. The block Gauss-Seidel (BGS) numerical method is implemented because the system of equations within the electrodes is coupled only through the nonlinear Butler–Volmer equation, which governs the electrochemical reaction at the interface between the domains. The best solution strategy found in this work consists of splitting the system into two blocks—one for the concentration and one for the potential field—and then performing block generalized minimal residual preconditioned with algebraic multigrid, using the FEniCS and the Portable, Extensible Toolkit for Scientific Computation libraries. Significant improvements in terms of time to solution (six times faster) and memory usage (halving) are achieved compared with the MUltifrontal Massively Parallel sparse direct Solver. Additionally, BGS experiences decent strong parallel scaling within the electrode domains. Last, the system of equations is modified to specifically address numerical instability induced by electrolyte depletion, which is particularly valuable for simulating fast-charge scenarios relevant for automotive application.

A Novel h–φ Approach for Solving Eddy–Current Problems in Multiply Connected Regions

A novel h-φ approach for solving 3-D time-harmonic eddy current problems is presented. It makes it possible to limit the number of degrees of freedom required for the discretization such as the T-Ω method, while overcoming topological issues related to it when multiply connected domains are considered. Global basis functions, needed for representing magnetic field in the insulating region, are obtained by a fast iterative solver. The computation of thick cuts by high-complexity computational topology tools, typically required by the T-Ω method, is thus avoided. The final matrix system turns out to be symmetric and full-rank unlike the more classical A-A method, which requires gauging of magnetic vector potential to ensure uniqueness. Numerical tests show that the proposed method is accurate and the field problem solution is obtained in a reasonable computational time even for 3-D models with millions of mesh elements.

Fast upwind and Eulerian-Lagrangian control volume schemes for time-dependent directional space-fractional advection-dispersion equations

We develop control volume methods for two-dimensional time-dependent advection-dominated directional space-fractional advection-dispersion equations with the directional space-fractional derivative weighted in all the directions by a probability measure in the unit circle, which are used to model the anisotropic superdiffusive transport of solutes in groundwater moving through subsurface heterogeneous porous media.We develop a fast upwind control volume method for the governing equation to eliminate the spurious numerical oscillations that often occur in space-centered numerical discretizations of advection term, which are relatively straightforward to implement. We also develop a Eulerian-Lagrangian control-volume method for the governing equation, which symmetrizes the governing equation by combining the time-derivative term and the advection term into a material derivative term along characteristic curves. Both methods are locally mass-conservative, which are essential in these applications.Due to the nonlocal nature of the directional space-fractional differential operators, corresponding numerical discretizations usually generate full stiffness matrices. Conventional direct solvers tend to require O(N2) memory requirement and have O(N3) computational complexity per time step, where N is the number of spatial unknowns, which is computationally significantly more expensive than the numerical approximations of integer-order advection-diffusion equations. Based on the analysis of the structure of stiffness matrix, we propose a fast Krylov subspace iterative solver to accelerate the numerical approximations of both the upwind and Eulerian-Lagrangian control volume methods, which reduce computational complexity from O(N3) by a direct solver to O(Nlog⁡N) per Krylov subspace iteration per time step and a memory requirement from O(N2) to O(N). Numerical results are presented to show the utility of the methods.

An embedded–hybridized discontinuous Galerkin finite element method for the Stokes equations

We present and analyze a new embedded–hybridized discontinuous Galerkin finite element method for the Stokes problem. The method has the attractive properties of full hybridized methods, namely an H(div)-conforming velocity field, pointwise satisfaction of the continuity equation and a priori error estimates for the velocity that are independent of the pressure. The embedded–hybridized formulation has advantages over a full hybridized formulation in that it has fewer global degrees-of-freedom for a given mesh and the algebraic structure of the resulting linear system is better suited to fast iterative solvers. The analysis results are supported by a range of numerical examples that demonstrate rates of convergence, and which show computational efficiency gains over a full hybridized formulation.

Optimization of Patterned Surfaces for Improved Superhydrophobicity through Cost-Effective Large-Scale Computations.

The pattern design of superhydrophobic surfaces can be significantly aided by computations that predict the Cassie-Baxter (CB) to Wenzel (W) transition, which is responsible for the break-down of superhydrophobic behavior. We present a computational framework for the optimization of patterned surfaces based on the energy barriers of the CB-W transitions which comprises the following elements: (a) design of structured surface patterns, for example, arrays of pillars, with parameterized geometric features such as size, pitch, slope, and roundness. (b) Computation of the wetting states with a modified Young-Laplace equation that facilitates the introduction of solid/liquid interactions for complex surface patterns and has significantly lower computational cost than other commonly used methods, such as the volume-of-fluid, phase-field, and so forth. (c) Incorporation of the modified Young-Laplace in the simplified string method, allowing the calculation of the minimum energy paths of wetting transitions which, apart from the energy barriers, also reveal the transition mechanisms (CB failure modes). (d) Accommodation of large-scale problems with good parallel performance and scalability on multicore-distributed memory systems using fast iterative solvers and the Message Passing Interface communication protocol. We demonstrate the computational framework with a shape optimization study of inverted conical frustum pillars. The optimization objective function is the resistance to the CB-W transition, which is quantified by the energy barrier-a relatively large energy barrier suggests improved superhydrophobicity. We also report the parallel performance, in terms of parallel speedup for problems ranging from three hundred thousands to 12 million degrees of freedom, solved using up to 40 processing cores.