Krylov Methods Research Articles

Mathematical modelling of nonlinear oscillations of a biodynamical system in the form of a complex cantilever

To ensure the stability of a young seedling of a tree and proper rooting, tree stems are usually staked with one, two or three stakes. In this paper, we present a mathematical model and its approximations of a nonlinear biodynamical system in the form of a complex discrete structure on a cantilever. The nonlinear biodynamical system corresponds to a staked tree with one stake. The geometric nonlinearity of the system is introduced by a spring with cubic nonlinear properties that oscillates in the vertical plane. The influence coefficients of deflection of the cantilever were used for describing single-frequency forced oscillations of a complex biodynamical system under an external single-frequency force. This external single-frequency force has a circular frequency that is close to frequencies from a set of eigen circular frequencies of a corresponding linearized system. The system oscillates in two orthogonal planes - the horizontal and the vertical plane. Forced oscillations, in the vertical direction, of this discrete complex system, are described by subsystems of nonlinear differential equations that are solved in the first approximation by using two methods: the method of variation of constants in combination with the method of averaging and the other extended asymptotic method of nonlinear mechanics Krylov- Bogoliubov- Mitropolsky. The results of the qualitative analysis of the solution of the system of nonlinear differential equations by amplitudes and phases of four nonlinear modes are presented and discussed. In both cases, the frequencies of external force are within the range of resonant frequencies, but constant in time in the case when the generalized method of variation of constants with the method of averaging was used. In the case when the Krylov-Bogoliubov- Mitropolsky asymptotic method of nonlinear mechanics of approximation was used, the frequencies of external force have changeable values within time- frequency changes with different speeds.This system is simpler than a system of governing differential equations along amplitudes and phases of generalized coordinates, and permits a qualitative analysis of nonlinear phenomena in systems with nonlinear dynamics. The method described in this paper is suitable for studying the stability and instability of amplitudes and of phases of nonlinear modes in the first approximation.The asymptotic method of nonlinear mechanics Krylov-Bogoliubov- Mitropolsky allows the analysis of stationary discrete modes of nonlinear small oscillations in the resonant range as well as of non-stationary nonlinear modes. This method can be applied to studying nonlinear dynamics of more complex systems based on cantilevers with complex structure.The most valuable and important elements of the work are emphasized.

High-resolution cerebral blood flow simulation with a domain decomposition method and verified by the TCD measurement

Background: An efficient and accurate blood flow simulation can be useful for understanding many vascular diseases. Accurately resolving the blood flow velocity based on patient-specific geometries and model parameters is still a major challenge because of complex geomerty and turbulence issues. In addition, obtaining results in a short amount of computing time is important so that the simulation can be used in the clinical environment. In this work, we present a parallel scalable method for the patient-specific blood flow simulation with focuses on its parallel performance study and clinical verification.Methods: We adopt a fully implicit unstructured finite element method for a patient-specific simulation of blood flow in a full precerebral artery. The 3D artery is constructed from MRI images, and a parallel Newton–Krylov method preconditioned with a two-level domain decomposition method is adopted to solve the large nonlinear system discretized from the time-dependent 3D Navier-Stokes equations in the artery with an integral outlet boundary condition. The simulated results are verified using the clinical data measured by transcranial Doppler ultrasound, and the parallel performance of the algorithm is studied on a supercomputer.Results: The simulated velocity matches the clinical measured data well. Other simulated blood flow parameters, such as pressure and wall shear stress, are within reasonable ranges. The results show that the parallel algorithm scales up to 2160 processors with a 49% parallel efficiency for solving a problem with over 20 million unstructured elements on a supercomputer. For a standard cerebral blood flow simulation case with approximately 4 million finite elements, the calculation of one cardiac cycle can be finished within one hour with 1000 processors.Conclusion: The proposed method is able to perform high-resolution 3D blood flow simulations in a patient-specific full precerebral artery within an acceptable time, and the simulated results are comparable with the clinical measured data, which demonstrates its high potential for clinical applications.

Comparative study of inner–outer Krylov solvers for linear systems in structured and high-order unstructured CFD problems

Advanced Krylov subspace methods are investigated for the solution of large sparse linear systems arising from stiff adjoint-based aerodynamic shape optimization problems. A special attention is paid to the flexible inner–outer GMRES strategy combined with most relevant preconditioning and deflation techniques. The choice of this specific class of Krylov solvers for challenging problems is based on its outstanding convergence properties. Typically in our implementation the efficiency of the preconditioner is enhanced with a domain decomposition method with overlapping. However, maintaining the performance of the preconditioner may be challenging since scalability and efficiency of a preconditioning technique are properties often antagonistic to each other. In this paper we demonstrate how flexible inner–outer Krylov methods are able to overcome this critical issue. A numerical study is performed considering either a Finite Volume (FV), or a high-order Discontinuous Galerkin (DG) discretization which affect the arithmetic intensity and memory-bandwidth of the algebraic operations. We consider test cases of transonic turbulent flows with RANS modeling over the two-dimensional supercritical OAT15A airfoil and the three-dimensional ONERA M6 wing. Benefits in terms of robustness and convergence compared to standard GMRES solvers are obtained. Strong scalability analysis shows satisfactory results. Based on these representative problems a discussion of the recommended numerical practices is proposed.

Numerical Simulation of the Kelvin Wake Patterns

The ship wave is of great interest for wave drag and coastal erosion. This paper proposes a mechanism of ship wave transformation to explore the effects of ship speed and ship size on the waveform. Firstly, based on the theory of potential flow, the boundary integral equations for the Kelvin ship waves are obtained by deploying the different Kelvin sources or Rankine sources. Then, these integral equations are numerically discretized to a set of nonlinear equations. Finally, the Jacobian−free Newton–Krylov method with a preconditioner is adopted to solve the nonlinear equations. Though imitating plenty of different Kelvin wave patterns, the mechanism of ship wave transformation is proposed to conveniently generate the polymorphic Kelvin wake patterns. The above numerical simulation scheme is verified by comparing simulation results with real ship waves. After that, the wake angle is discussed with the effects of Froude number, source strength and source type by following the mechanism of ship wave transformation. The results show that the wake angle tends to decrease with ship speed but increase with ship size. In addition, for high ship speeds, the effect on the wake angle can be more dramatic.

EI-NK: A Robust Exponential Integrator Method With Singularity Removal and Newton–Raphson Iterations for Transient Nonlinear Circuit Simulation

In this article, we propose an exponential-integrator-Newton–Krylov (EI-NK) method for transient simulation of large-scale nonlinear circuits. This method aims to address two long-standing problems that affect the performance and robustness of EI in handling nonlinear circuits. First is the numerical instability caused by the singularity in the differential-algebraic equation system. We provide an in-depth analysis of the problem and propose a systematic, algebraic, and sparsity preserving regularization technique to eliminate the unstable modes in the system to be solved. Next, we develop a scheme to enable the iterative Newton–Raphson solution in the EI framework for enhanced nonlinearity handling capability. With the two techniques, we wish to elevate EI’s robustness and performance and make it a competitive alternative to the existing SPICE-type simulators in practical applications.

Preconditioners for fractional diffusion equations based on the spectral symbol

AbstractIt is well known that the discretization of fractional diffusion equations with fractional derivatives , using the so‐called weighted and shifted Grünwald formula, leads to linear systems whose coefficient matrices show a Toeplitz‐like structure. More precisely, in the case of variable coefficients, the related matrix sequences belong to the so‐called generalized locally Toeplitz class. Conversely, when the given FDE has constant coefficients, using a suitable discretization, we encounter a Toeplitz structure associated to a nonnegative function , called the spectral symbol, having a unique zero at zero of real positive order between one and two. For the fast solution of such systems by preconditioned Krylov methods, several preconditioning techniques have been proposed in both the one‐ and two‐dimensional cases. In this article we propose a new preconditioner denoted by which belongs to the ‐algebra and it is based on the spectral symbol . Comparing with some of the previously proposed preconditioners, we show that although the low band structure preserving preconditioners are more effective in the one‐dimensional case, the new preconditioner performs better in the more challenging multi‐dimensional setting.

A Comparison of Linear Solvers for Resolving Flow in Three‐Dimensional Discrete Fracture Networks

AbstractWe compare various methods for resolving steady flow within three‐dimensional discrete fracture networks, including direct methods, Krylov subspace methods with and without preconditioning, and multi‐grid methods. We compared the performance of the methods based on compute times and scaling of the solution as a function of the number of grid nodes and log‐variance of the hydraulic aperture. The methods are applied to three test cases: (a) variable density of networks with a truncated power‐law distribution of fracture lengths, (b) a fixed network composed of monodisperse fracture sizes but varied permeability/aperture heterogeneity, (c) and a network based on field site in Nevada, US. We chose these cases to allow us to study the impact of the mesh size and flow properties, as well as to demonstrate our conclusions on a large‐scale, realistic problem (more than 40 million mesh nodes). A direct solution using Cholesky factorization outperformed other methods for every example but was closely followed in performance by some algebraic multigrid (AMG) preconditioned Krylov subspace methods. Among the Krylov methods, conjugate gradients (CG) with an AMG preconditioner performs the best. Generally, Cholesky factorization is recommended, but CG with an AMG preconditioner may be suitable for very large problems beyond 40 million nodes where the entire linear system cannot reside in memory.

Interference of internal waves due to two point vortices: linear analytical solution and nonlinear interaction.

In this work, we consider steady two-dimensional interfacial waves in a two-layer stratified fluid, which is induced by a vortex pair located in the lower layer of the fluids. A mathematical model based on the boundary integral equation method and the potential-flow theory is established. The linear analytical solution for the linearized model is given in the form of Cauchy integral and then asymptotic behaviour for large x is presented. The fully nonlinear model is solved by the Jacobian-free Newton–Krylov (JFNK) method numerically. Nonlinear characteristics of wave profiles are identified compared with the linear results under different vortex strengths and the distance between the vortex pair. The amplitude of steady downstream waves is found to vary periodically with respect to the distance of the vortex pair, which can be regarded as the interference between waves produced by each vortex. For equal-strength counter- and co-rotating pairs, the downstream wave heights of linear solutions can be eliminated for some special values of the distance between point vortices, namely, the destructive interference occurs. Meanwhile, the wave only exists between the vortex pair like trapped waves. So does the nonlinear counterpart for counter-rotating pairs, but it could not be diminished with any distance.

Large-Eddy Simulation of Transonic Buffet Using Matrix-Free Discontinuous Galerkin Method

We present an implicit large-eddy simulation of transonic buffet over the OAT15A supercritical airfoil at Mach number 0.73, angle of attack 3.5 deg, and Reynolds number . The simulation is performed using a matrix-free discontinuous Galerkin (DG) method and a diagonally implicit Runge–Kutta scheme on graphics processor units. We propose a Jacobian-free Newton–Krylov method to solve nonlinear systems arising from the discretization of the Navier–Stokes equations. The method successfully predicts the buffet onset, the buffet frequency, and turbulence statistics owing to the high-order DG discretization and an efficient mesh refinement for the laminar and turbulent boundary layers. A number of physical phenomena present in the experiment are captured in our simulation, including periodical low-frequency oscillations of shock wave in the streamwise direction, strong shear layer detached from the shock wave due to shock-wave/boundary-layer interaction and small-scale structures broken down by the shear-layer instability in the transition region, and shock-induced flow separation. The pressure coefficient, the root mean square of the fluctuating pressure, and the streamwise range of the shock wave oscillation agree well with the experimental data. The results suggest that the proposed method can accurately predict the onset of turbulence and buffet phenomena at high Reynolds numbers without a subgrid scale model or a wall model.

Parallel Newton--Krylov BDDC and FETI-DP Deluxe Solvers for Implicit Time discretizations of the Cardiac Bidomain Equations

Two novel parallel Newton-Krylov Balancing Domain Decomposition by Constraints (BDDC) and Dual-Primal Finite Element Tearing and Interconnecting (FETI-DP) solvers are here constructed, analyzed and tested numerically for implicit time discretizations of the three-dimensional Bidomain system of equations. This model represents the most advanced mathematical description of the cardiac bioelectrical activity and it consists of a degenerate system of two non-linear reaction-diffusion partial differential equations (PDEs), coupled with a stiff system of ordinary differential equations (ODEs). A finite element discretization in space and a segregated implicit discretization in time, based on decoupling the PDEs from the ODEs, yields at each time step the solution of a non-linear algebraic system. The Jacobian linear system at each Newton iteration is solved by a Krylov method, accelerated by BDDC or FETI-DP preconditioners, both augmented with the recently introduced {\em deluxe} scaling of the dual variables. A polylogarithmic convergence rate bound is proven for the resulting parallel Bidomain solvers. Extensive numerical experiments on linux clusters up to two thousands processors confirm the theoretical estimates, showing that the proposed parallel solvers are scalable and quasi-optimal.

Low-CP-Rank Tensor Completion via Practical Regularization

Dimension reduction is analytical methods for reconstructing high-order tensors that the intrinsic rank of these tensor data is relatively much smaller than the dimension of the ambient measurement space. Typically, this is the case for most real world datasets in signals, images and machine learning. The CANDECOMP/PARAFAC (CP, aka Canonical Polyadic) tensor completion is a widely used approach to find a low-rank approximation for a given tensor. In the tensor model (Sanogo and Navasca in 2018 52nd Asilomar conference on signals, systems, and computers, pp 845–849, https://doi.org/10.1109/ACSSC.2018.8645405, 2018), a sparse regularization minimization problem via ell _1 norm was formulated with an appropriate choice of the regularization parameter. The choice of the regularization parameter is important in the approximation accuracy. Due to the emergence of the massive data, one is faced with an onerous computational burden for computing the regularization parameter via classical approaches (Gazzola and Sabaté Landman in GAMM-Mitteilungen 43:e202000017, 2020) such as the weighted generalized cross validation (WGCV) (Chung et al. in Electr Trans Numer Anal 28:2008, 2008), the unbiased predictive risk estimator (Stein in Ann Stat 9:1135–1151, 1981; Vogel in Computational methods for inverse problems, 2002), and the discrepancy principle (Morozov in Doklady Akademii Nauk, Russian Academy of Sciences, pp 510–512, 1966). In order to improve the efficiency of choosing the regularization parameter and leverage the accuracy of the CP tensor, we propose a new algorithm for tensor completion by embedding the flexible hybrid method (Gazzola in Flexible krylov methods for lp regularization) into the framework of the CP tensor. The main benefits of this method include incorporating the regularization automatically and efficiently as well as improving accuracy in the reconstruction and algorithmic robustness. Numerical examples from image reconstruction and model order reduction demonstrate the efficacy of the proposed algorithm.

QForte: An Efficient State-Vector Emulator and Quantum Algorithms Library for Molecular Electronic Structure.

We introduce a novel open-source software package QForte, a comprehensive development tool for new quantum simulation algorithms. QForte incorporates functionality for handling molecular Hamiltonians, Fermionic encoding, ansatz construction, time evolution, and state-vector emulation, requiring only a classical electronic structure package as a dependency. QForte also contains black-box implementations of a wide variety of quantum algorithms, including variational and projective quantum eigensolvers, adaptive eigensolvers, quantum imaginary time evolution, and quantum Krylov methods. We highlight two features of QForte: (i) how the Python class structure of QForte enables the facile implementation of new algorithms, and (ii) how existing algorithms can be executed in just a few lines of code.

Steady-State Analysis of Electrical Networks in Pandapower Software: Computational Performances of Newton–Raphson, Newton–Raphson with Iwamoto Multiplier, and Gauss–Seidel Methods

At the core of every system for the efficient control of the network steady-state operation is the AC-power-flow problem solver. For local distribution networks to continue to operate effectively, it is necessary to use the most powerful and numerically stable AC-power-flow problem solvers within the software that controls the power flows in these networks. This communication presents the results of analyses of the computational performance and stability of three methods for solving the AC-power-flow problem. Specifically, this communication compares the robustness and speed of execution of the Gauss–Seidel (G–S), Newton–Raphson (N–R), and Newton–Raphson method with Iwamoto multipliers (N–R–I), which were tested in open-source pandapower software using a meshed electrical network model of various topologies. The test results show that the pandapower implementations of the N–R method and the N–R–I method are significantly more robust and faster than the G–S method, regardless of the network topology. In addition, a generalized Python interface between the pandapower and the SciPy package was implemented and tested, and results show that the hybrid Powell, Levenberg–Marquardt, and Krylov methods, a quasilinearization algorithm, and the continuous Newton method can sometimes achieve better results than the classical N–R method.

Jacobian-free Newton Krylov two-node coarse mesh finite difference based on nodal expansion method

A Jacobian-Free Newton Krylov Two-Nodal Coarse Mesh Finite Difference algorithm based on Nodal Expansion Method (NEM_TNCMFD_JFNK) is successfully developed and proposed to solve the three-dimensional (3D) and multi-group reactor physics models. In the NEM_TNCMFD_JFNK method, the efficient JFNK method with the Modified Incomplete LU (MILU) preconditioner is integrated and applied into the discrete systems of the NEM-based two-node CMFD method by constructing the residual functions of only the nodal average fluxes and the eigenvalue. All the nonlinear corrective nodal coupling coefficients are updated on the basis of two-nodal NEM formulation including the discontinuity factor in every few newton steps. All the expansion coefficients and interface currents of the two-node NEM need not be chosen as the solution variables to evaluate the residual functions of the NEM_TNCMFD_JFNK method, therefore, the NEM_TNCMFD_JFNK method can greatly reduce the number of solution variables and the computational cost compared with the JFNK based on the conventional NEM. Finally the NEM_TNCMFD_JFNK code is developed and then analyzed by simulating the representative PWR MOX/UO2 core benchmark, the popular NEACRP 3D core benchmark and the complicated full-core pin-by-pin homogenous core model. Numerical solutions show that the proposed NEM_TNCMFD_JFNK method with the MILU preconditioner has the good numerical accuracy and can obtain higher computational efficiency than the NEM-based two-node CMFD algorithm with the power method in the outer iteration and the Krylov method using the MILU preconditioner in the inner iteration, which indicates the NEM_TNCMFD_JFNK method can serve as a potential and efficient numerical tool for reactor neutron diffusion analysis module in the JFNK-based multiphysics coupling application.

Fast Solution of Fully Implicit Runge--Kutta and Discontinuous Galerkin in Time for Numerical PDEs, Part I: the Linear Setting

Fully implicit Runge-Kutta (IRK) methods have many desirable properties as time integration schemes in terms of accuracy and stability, but are rarely used in practice with numerical PDEs due to the difficulty of solving the stage equations. This paper introduces a theoretical and algorithmic preconditioning framework for solving the systems of equations that arise from IRK methods applied to linear numerical PDEs (without algebraic constraints). This framework also naturally applies to discontinuous Galerkin discretizations in time. Under quite general assumptions on the spatial discretization that yield stable time integration, the preconditioned operator is proven to have conditioning ~O(1), with only weak dependence on number of stages/polynomial order; for example, the preconditioned operator for 10th-order Gauss IRK has condition number less than two, independent of the spatial discretization. The new method can be used with arbitrary existing preconditioners for backward Euler-type time stepping schemes, and is amenable to the use of three-term recursion Krylov methods when the underlying spatial discretization is symmetric. The new method is demonstrated to be effective on various high-order finite-difference and finite-element discretizations of linear parabolic and hyperbolic problems, demonstrating fast, scalable solution of up to 10th order accuracy. The new method consistently outperforms existing block preconditioning approaches, and in several cases, the new method can achieve 4th-order accuracy using Gauss integration with roughly half the number of preconditioner applications and wallclock time as required using standard diagonally implicit RK methods.

The Effect of Trim on Stability and Seakeeping of Tanker, Container, and Bulk Carrier

Stability and seakeeping are important factor that must be owned by a ship because it related to the safety of the ship, meanwhile the trim also affecting the ship operation. The study has a purpose to find out the impact of the trim to the stability and seaworthiness of the ship. The analysis is carried out on 3 types of ships, namely tankers, containers, and bulk carriers. A.N. Krylov method is using to calculate the stability, whereas the seakeeping is analysed by strip theory. Both of the results for stability and seakeeping are calculated by Maxsurf Software. Stability and seakeeping analysis was carried out when the ship is in even keel condition, trim by bow, and trim by stern in maximum and extreme condition. The results shows that the trim by stern of the ship can increased the ship stability rate by 0.5~5.4%. While the effect of trim on ship motions varies due to the differences characteristics of each ship. The results shows due to the trim impact by values of RAO’s, which some of ship motions have increased meanwhile for some ship motion have decreased of ship motion condition.

Fast 3-D Controlled-Source Electromagnetic Modeling Combining UPML and Rational Krylov Method

Controlled-source electromagnetic (CSEM) surveying is a critical tool for sensing and locating underground anomalies and structures. In this letter, based on uniaxial perfectly matched layer (UPML) and rational Krylov (RK) method, we propose a fast algorithm for 3-D Multifrequency CSEM modeling. We use the frequency-independent UPML to truncate the boundaries and adopt an RK method to rapidly solve the 3-D multifrequency CSEM problems. The accuracy and efficiency of our algorithm are verified by two examples, i.e., a two-layer model and a 3-D model. Numerical experiments indicate that our algorithm is computationally efficient, obtaining nearly 20-fold speedup on a laptop compared with the conventional 3-D CSEM using finite element method (3DCSEM) modeling.

3-D Large-Scale TEM Modeling Using Restarting Polynomial Krylov Method

The transient electromagnetic (TEM) method is widely used in near-surface geophysical prospecting. The high-precision forward and inversion of large-scale complex models is a challenging problem. It is difficult to solve the large-scale problems using the direct solvers due to the large memory requirements. In this article, we present a restarting polynomial Krylov method for modeling 3-D large-scale TEM responses. The mimetic finite volume method is carried out for spatial discretization. The step-off TEM response then can be expressed as a matrix exponential function. The restarting polynomial Krylov method is used to solve the matrix exponential function. For a given restart subspace dimension, the residual is used to obtain the forward response at any time that meets the given accuracy. This method does not need to solve large-scale linear equations. Furthermore, the memory usage is mainly determined by the number of spatial discrete grids and the restart subspace order. The numerical experiments of large-scale model demonstrated that the method is accurate and uses limited memory.

Jacobian-Free Newton–Krylov Method for the Simulation of Non-Thermal Plasma Discharges with High-Order Time Integration and Physics-Based Preconditioning

Jacobian-Free Newton–Krylov Method for the Simulation of Non-Thermal Plasma Discharges with High-Order Time Integration and Physics-Based Preconditioning

A Compact Heart Iteration for Large Eigenvalues Problems

In this paper, we present a compact version of the Heart iteration. One that requires less matrix-vector products per iteration and attains faster convergence. The Heart iteration is a new type of Restarted Krylov methods for calculating peripheral eigenvalues of symmetric matrices. The new framework avoids the Lanczos tridiagonalization process and the use of implicit restarts. This simplifies the restarting mechanism and allows the introduction of several modifications. Convergence is assured by a monotonicity property that pushes the computed Ritz values toward their limits. Numerical experiments illustrate the usefulness of the proposed approach.