BiCGSTAB Solver Research Articles

Flexible GMRES solver for interpolating MLPG analysis of heat conduction

This article evaluates the performance of flexible generalized minimum residual (FGMRES) solver in meshless local Petrov–Galerkin (MLPG) analysis of heat conduction. Interpolating MLPG formulation has been used for the meshless discretization, resulting in an asymmetric algebraic equation system. FGMRES is a variant of GMRES, which allows the preconditioner to vary from step to step. Also, the GMRES solver can be used as a preconditioner. The FGMRES method has been used with Jacobi, SOR(k), and GMRES preconditioners for MLPG analysis of heat conduction in two-dimensional and three-dimensional geometries. The performance of the FGMRES solver is compared with the restarted GMRES solver (and BiCGSTAB solver in 3-D complex geometry). Jacobi and SOR(k) preconditioners are used in the restarted GMRES and BICGSTAB solvers. Results demonstrate that the FGMRES solver performs better than the GMRES and BiCGSTAB solvers, except for two-dimensional complex geometry.

Reliability Investigation of BiCGStab and IDR Solvers for the Advection-Diffusion-Reaction Equation

Reliability Investigation of BiCGStab and IDR Solvers for the Advection-Diffusion-Reaction Equation

Hardware Acceleration of High-Performance Computational Flow Dynamics Using High-Bandwidth Memory-Enabled Field-Programmable Gate Arrays

Scientific computing is at the core of many High-Performance Computing applications, including computational flow dynamics. Because of the utmost importance to simulate increasingly larger computational models, hardware acceleration is receiving increased attention due to its potential to maximize the performance of scientific computing. Field-Programmable Gate Arrays could accelerate scientific computing because of the possibility to fully customize the memory hierarchy important in irregular applications such as iterative linear solvers. In this article, we study the potential of using Field-Programmable Gate Arrays in High-Performance Computing because of the rapid advances in reconfigurable hardware, such as the increase in on-chip memory size, increasing number of logic cells, and the integration of High-Bandwidth Memories on board. To perform this study, we propose a novel Sparse Matrix-Vector multiplication unit and an ILU0 preconditioner tightly integrated with a BiCGStab solver kernel. We integrate the developed preconditioned iterative solver in Flow from the Open Porous Media project, a state-of-the-art open source reservoir simulator. Finally, we perform a thorough evaluation of the FPGA solver kernel in both stand-alone mode and integrated in the reservoir simulator, using the NORNE field, a real-world case reservoir model using a grid with more than 10 5 cells and using three unknowns per cell.

Algebraic multigrid for the finite pointset method

We investigate algebraic multigrid (AMG) methods for the linear systems arising from the discretization of Navier–Stokes equations via the finite pointset method. In the segregated approach, three pressure systems and one velocity system need to be solved. In the coupled approach, one of the pressure systems is coupled with the velocity system, leading to a coupled velocity-pressure saddle point system. The discretization of the differential operators used in FPM leads to non-symmetric matrices that do not have the M-matrix property. Even though the theoretical framework for many AMG methods requires these properties, our AMG methods can be successfully applied to these matrices and show a robust and scalable convergence when compared to a BiCGStab(2) solver.

An Efficient Preconditioner for 3-D Finite Difference Modeling of the Electromagnetic Diffusion Process in the Frequency Domain

Krylov subspace solvers for frequency-domain electromagnetic forward modeling problems converge remarkably more slowly as the period increases. In this article, we present an efficient four-color cellblock Gauss Seidel (GS) preconditioner for finite-difference (FD) electromagnetic modeling in geophysical applications. Rather than updating the FD electromagnetic (EM) equation edge by edge, as in a traditional GS scheme, we renew six edge components attached to one node simultaneously (i.e., in cellblock manner) effectively enforcing a local divergence free condition for currents. To improve implementation efficiency, we reorder the nodes on the FD grid into four colors so that nodes in each color are uncoupled, allowing the use of highly parallel vectorized algorithms. The four-color cellblock GS preconditioner is implemented in the MATLAB code, in conjunction with a BiCGstab solver. It is compared, in terms of iteration number and computing time, with other three commonly used preconditioners [GS, symmetric successive overrelaxation (SSOR) and incomplete lower and upper triangular matrix decomposition (ILU)] on three models—two synthetic and one modified from the version of real data inversion. The comparison indicates that the proposed algorithm is extremely stable and efficient compared with the other three preconditioners tested, over a range of periods (1-1000 s). Especially at long periods, the improvement of our proposed algorithm is substantial. In addition, a parallel implementation of the cellblock GS preconditioner is straightforward due to the independence of nodes in each color.

Domain decomposition based parallel computing for multi-scale coronary blood flow simulations

In the present study, we demonstrated a parallel P2P1 finite element scheme with a four-step fractional splitting approach to conduct a multiscale coronary flow simulation. The three-dimensional (3D) computational fluid dynamics (CFD) for patient-specific coronary artery flow and zero-dimensional (0D) lumped-parameter network (LPN) modeling for distal coronary beds was fully coupled, and an MPI parallel algorithm based on domain decomposition was applied. A parallel conjugate gradient (CG)-LPN subroutine for the 3D-0D coupled system with a monolithic scheme was derived, and it provides a correct pressure solution that may not be obtained by the conventional CG solver, particularly when the subdomain division intersects the 3D-0D coupling outlet. The overall computing time for parallel CG-LPN does not show a noticeable difference from the conventional CG solver, despite the extra MPI calls for data transfer at the interfaces of subdomains.For the BiCGSTAB solver, the block ILU(0) preconditioner showed favorable performance for a high density mesh compared with the simple Jacobi preconditioner. MPI_COMM is a major bottleneck that saturates the overall parallel performance at high core count, but a computing time of less than 10 min per cardiac cycle on a medium density mesh could be attained for a patient-specific coronary flow simulation using 60 CPU cores run in parallel, which is in an acceptable range for clinical practice. Further tests of accuracy are needed with a large set of patients to enable wide use of the proposed technique in helping to make interventional decisions in routine clinical practice for coronary stenotic lesions.

Stability and convergence of a Crank–Nicolson finite volume method for space fractional diffusion equations

We analyze a Crank–Nicolson finite volume method (CN-FVM) for the time-dependent two-sided conservative space fractional diffusion equation (of order 2−α). We prove that the proposed method is unconditionally stable in a weighted discrete norm and has a convergence rate of order O(τ2+h1+α), where τ and h are the time step size and spatial mesh size, respectively. In addition, we present a matrix-free preconditioned fast BiCGSTAB solver for the discrete linear algebraic system, which has a linear memory requirement and almost linear computational complexity. Numerical experiments show strong potential of the fast CN-FVM.

On preconditioned BiCGSTAB solver for MLPG method applied to heat conduction in 3D complex geometry

Development of meshfree methods is in full swing in the last few decades as these methods overcome the shortcomings of mesh based methods. Development of meshfree methods for large scale problems requires fast iterative solvers. Preconditioner plays an important role in improving convergence rate of an iterative solver. Recently proposed SRJ method can be used as a stand-alone solver and as a preconditioner in CG type methods. Objective of the present work is to evaluate the performance of SRJ method as a preconditioner in BiCGSTAB method in connection with the MLPG method. Three dimensional Poisson equation and heat conduction equation are solved in regular and irregular domain respectively. The SRJ method is tested as a solver and as a preconditioner in BiCGSTAB method. Its performance as a preconditioner is compared with Jacobi and SOR preconditioners.

On preconditioned BiCGSTAB solver for MLPG method applied to heat conduction in complex geometry

ABSTRACTMeshless local Petrov–Galerkin (MLPG) method is a promising meshfree method for continuum problems in complex domains, especially for large deformation, moving boundary and phase change problems. For large-scale problems, iterative methods for solving the discretized equations are more suitable than direct methods. Krylov subspace solvers of conjugate gradient type are the most preferred iterative solvers. The convergence rate of these methods depends on preconditioner used. Recently, proposed schedule relaxation Jacobi (SRJ) method can be used as a stand-alone solver and as a preconditioner. In the present work, the SRJ method is tested as a stand-alone solver and as a preconditioner for BiCGSTAB solver in the MLPG method, and its performance has been compared with successive overrelaxation (k) preconditioner. Two-dimensional linear steady-state heat conduction in complex shape geometry has been used as the model test problem.

Preconditioned Krylov solvers on GPUs

Abstract In this paper, we study the effect of enhancing GPU-accelerated Krylov solvers with preconditioners. We consider the BiCGSTAB, CGS, QMR, and IDR(s) Krylov solvers. For a large set of test matrices, we assess the impact of Jacobi and incomplete factorization preconditioning on the solvers’ numerical stability and time-to-solution performance. We also analyze how the use of a preconditioner impacts the choice of the fastest solver.

First Application of Lattice QCD to Pezy-SC Processor

Pezy-SC processor is a novel new architecture developed by Pezy Computing K. K. that has achieved large computational power with low electric power consumption. It works as an accelerator device similarly to GPGPUs. A programming environment that resembles OpenCL is provided. Using a hybrid parallel system “Suiren” installed at KEK, we port and tune a simulation code of lattice QCD, which is computational elementary particle physics based on Monte Carlo method. We offload an iterative solver of a linear equation for a fermion matrix, which is in general the most time consuming part of the lattice QCD simulations. On single and multiple Pezy-SC devices, the sustained performance is measured for the matrix multiplications and a BiCGStab solver. We examine how the data layout affects the performance. The results demonstrate that the Pezy-SC processors provide a feasible environment to perform numerical lattice QCD simulations.

Radiation magnetohydrodynamics in global simulations of protoplanetary discs

Our aim is to study the thermal and dynamical evolution of protoplanetary disks in global simulations, including the physics of radiation transfer and magneto-hydrodynamic (MHD) turbulence caused by the magneto-rotational instability. We develop a radiative transfer method based on the flux-limited diffusion approximation that includes frequency dependent irradiation by the central star. This hybrid scheme is implemented in the PLUTO code. The focus of our implementation is on the performance of the radiative transfer method. Using an optimized Jacobi preconditioned BiCGSTAB solver, the radiative module is three times faster than the MHD step for the disk setup we consider. We obtain weak scaling efficiencies of 70% up to 1024 cores. We present the first global 3D radiation MHD simulations of a stratified protoplanetary disk. The disk model parameters are chosen to approximate those of the system AS 209 in the star-forming region Ophiuchus. Starting the simulation from a disk in radiative and hydrostatic equilibrium, the magnetorotational instability quickly causes MHD turbulence and heating in the disk. For the disk parameters we use, turbulent dissipation heats the disk midplane and raises the temperature by about 15% compared to passive disk models. A roughly flat vertical temperature profile establishes in the disk optically thick region close to the midplane. We reproduce the vertical temperature profile with a viscous disk models for which the stress tensor vertical profile is flat in the bulk of the disk and vanishes in the disk corona. The present paper demonstrates for the first time that global radiation MHD simulations of turbulent protoplanetary disks are feasible with current computational facilities. This opens up the windows to a wide range of studies of the dynamics of protoplanetary disks inner parts, for which there are significant observational constraints.

Evaluation and FPGA Implementation of Sparse Linear Solvers for Video Processing Applications

Sparse linear systems are commonly used in video processing applications, such as edge-aware filtering or video retargeting. Due to the 2-D nature of images, the involved problem sizes are large and thus solving such systems is computationally challenging. In this paper, we address sparse linear solvers for real-time video applications. We investigate several solver techniques, discuss hardware trade-offs, and provide field-programmable gate array (FPGA) architectures and implementation results of a Cholesky direct solver and of an iterative BiCGSTAB solver. The FPGA implementations solve 32 k × 32 k matrices at up to 50 f/s and outperform software implementations by at least one order of magnitude.

An efficient block parallel AMR method for two phase interfacial flow simulations

The direct numerical simulation of two phase interfacial flows can be computationally challenging, as the strong resolution needed to follow the deformations of the interface leads to a lot of time spent solving the whole computation domain. Efficient solution of such problems requires an adaptive mesh refinement capability to concentrate computational effort where it is most needed. In this paper a parallel adaptive algorithm to solve incompressible two-phase flows with surface tension is presented: the AMR is handled with the help of the PARAMESH package. The free interface between fluids is tracked via Level Set approach; the jump conditions at the interface for pressure and velocity are imposed by the Ghost Fluid method. A multigrid preconditioned BiCG-stab solver adapted to the AMR data structure has been developed to allow high density ratio computations (up to 1:1000). Special treatment has been done at the refinement jumps to maintain the fine mesh accuracy. Computational results are compared in different test cases with analytical solutions or literature, and show very good agreement with the references. The effectiveness of PARAMESH parallelization has been quite well maintained, as shown in the strong and weak scaling tests. Speed-up capabilities of the AMR are demonstrated.

Higher order Galerkin time discretizations and fast multigrid solvers for the heat equation

We discuss numerical properties of continuous Galerkin–Petrov and discontinuous Galerkin time discretizations applied to the heat equation as a prototypical example for scalar parabolic partial differential equations. For the space discretization, we use biquadratic quadrilateral finite elements on general two-dimensional meshes. We discuss implementation aspects of the time discretization as well as efficient methods for solving the resulting block systems. Here, we compare a preconditioned BiCGStab solver as a Krylov space method with an adapted geometrical multigrid solver. Only the convergence of the multigrid method is almost independent of the mesh size and the time step leading to an efficient solution process. By means of numerical experiments we compare the different time discretizations with respect to accuracy and computational costs.

Parallel Newton-Krylov solvers for modeling of a navigation lock filling system

The Galerkin least-squares finite element method for solving the Reynolds-averaged incompressible turbulent 3-D Navier-Stokes equations is employed to simulate a navigation lock filling system in the numerical code Adaptive Hydraulics (ADH). The linear system is solved at each nonlinear iteration within every time-step using biconjugate gradient stabilized (BiCGstab) in combination with block-Jacobi (bjacobi) preconditioners, as it failed to solve the linear system because of dramatic changes in flow velocity and pressure early in the simulation. To overcome this problem, we used the Portable Extensible Toolkit for Scientific Computation (PETSc), a numerical library that provides multiple types of linear solvers. PETSc has been incorporated into the ADH code. The ADH-PETSc interface helps to systematically investigate the best linear solver for an ADH simulation. We found that a variant, known as enhanced BiCGstab(l) in combination with the additive Schwarz method (ASM), made it possible to simulate the John Day lock filling system. The BiCGstab(l) solver improved the rate of convergence because of a more reliable update strategy for the residuals. In addition, the simulation was run with various numbers of processors. The result shows good scaling of solution time as the number of processors increases

A sharp interface finite volume method for elliptic equations on Cartesian grids

We present a second order sharp interface finite volume method for the solution of the three-dimensional elliptic equation ∇·(β(x→)∇u(x→))=f(x→) with variable coefficients on Cartesian grids. In particular, we focus on interface problems with discontinuities in the coefficient, the source term, the solution, and the fluxes across the interface. The method uses standard piecewise trilinear finite elements for normal cells and a double piecewise trilinear ansatz for the solution on cells intersected by the interface resulting always in a compact 27-point stencil. Singularities associated with vanishing partial volumes of intersected grid cells are removed by a two-term asymptotic approach. In contrast to the 2D method presented by two of the authors in [M. Oevermann, R. Klein, A Cartesian grid finite volume method for elliptic equations with variable coefficients and embedded interfaces, Journal of Computational Physics 219 (2006) 749–769] we use a minimization technique to determine the unknown coefficients of the double trilinear ansatz. This simplifies the treatment of the different cut-cell types and avoids additional special operations for degenerated interface topologies. The resulting set of linear equations has been solved with a BiCGSTAB solver preconditioned with an algebraic multigrid. In various testcases – including large β-ratios and non-smooth interfaces – the method achieves second order of accuracy in the L∞ and L2 norm.

Three-dimensional magnetotellurics modeling using edgebased finite-element unstructured meshes

Three-dimensional forward modeling is a challenge for geometrically complex magnetotellurics (MT) problems. We present a new edge-based finite-element algorithm using an unstructured mesh for accurately and efficiently simulating 3D MT responses. The electric field curl-curl equation in the frequency domain was used to deduce the H (curl) variation weak form of the MT forward problem, the Galerkin rule was used to derive a linear finite-element equation on the linear-edge tetrahedroid space, and, finally, a BI-CGSTAB solver was used to estimate the unknown electric fields. A local mesh refinement technique in the neighbor of the measuring MT stations was used to greatly improve the accuracies of the numerical solutions. Four synthetic models validated the powerful performance of our algorithms. We believe that our method will effectively contribute to processing more complex MT studies.

A near-field preconditioner and its performance in conjunction with the BiCGstab(ell) solver

This paper describes a simple physically-motivated near-field preconditioning scheme that is effective in accelerating convergence of surface, volume, and combined surface/volume integral equations for a broad variety of electromagnetic scattering problems. It can be easily implemented numerically in method of moment (MoM) solvers (both conventional and those employing matrix-compression techniques), irrespective of the analytical form of the integral-equation kernel. It has low memory and CPU requirements, both of which scale linearly with the number of unknowns, and is easily amenable to efficient parallelization. We demonstrate the preconditioner's performance (in conjunction with the BiCGstab(ell) iterative solver) on two representative geometries, and observe a significant reduction in the number of iterations required for convergence.