Inexact Solvers Research Articles

Double saddle‐point preconditioning for Krylov methods in the inexact sequential homotopy method

AbstractWe derive an extension of the sequential homotopy method that allows for the application of inexact solvers for the linear (double) saddle‐point systems arising in the local semismooth Newton method for the homotopy subproblems. For the class of problems that exhibit (after suitable partitioning of the variables) a zero in the off‐diagonal blocks of the Hessian of the Lagrangian, we propose and analyze an efficient, parallelizable, symmetric positive definite preconditioner based on a double Schur complement approach. For discretized optimal control problems with PDE constraints, this structure is often present with the canonical partitioning of the variables in states and controls. We conclude with numerical results for a badly conditioned and highly nonlinear benchmark optimization problem with elliptic partial differential equations and control bounds. The resulting method allows for the parallel solution of large 3D problems.

On the implementation in Abaqus of the global–local iterative coupling and acceleration techniques

This paper presents results and convergence study of the Global–Local Iterative Coupling through the implementation in the commercial software Abaqus making use of the co-simulation engine. A hierarchical modeling and simulation approach is often required to alleviate modeling burdens. Particular focus has been devoted here on convergence acceleration and performance optimization. Two applications in statics with nonlinear material behavior and geometrically nonlinear formulation are considered here: first a holed curved plate under traction with elastic–plastic material, then a pre-stressed bolted joint connecting two plates between each other and subjected to traction load. Three different convergence acceleration techniques are compared in terms of convergence performance and accuracy. An inexact solver strategy is proposed to improve computing time performance. The results show promising results for the coupling technology and constitute a step forward in the availability of non-intrusive multi-scale modeling capabilities for complex structures and assemblies.

A posteriori error estimates for the Richards equation

The Richards equation is commonly used to model the flow of water and air through soil, and it serves as a gateway equation for multiphase flows through porous media. It is a nonlinear advection–reaction–diffusion equation that exhibits both parabolic–hyperbolic and parabolic–elliptic kind of degeneracies. In this study, we provide reliable, fully computable, and locally space–time efficient a posteriori error bounds for numerical approximations of the fully degenerate Richards equation. For showing global reliability, a nonlocal-in-time error estimate is derived individually for the time-integrated H 1 ( H − 1 ) H^1(H^{-1}) , L 2 ( L 2 ) L^2(L^2) , and the L 2 ( H 1 ) L^2(H^1) errors. A maximum principle and a degeneracy estimator are employed for the last one. Global and local space–time efficiency error bounds are then obtained in a standard H 1 ( H − 1 ) ∩ L 2 ( H 1 ) H^1(H^{-1})\cap L^2(H^1) norm. The reliability and efficiency norms employed coincide when there is no nonlinearity. Moreover, error contributors such as space discretization, time discretization, quadrature, linearization, and data oscillation are identified and separated. The estimates are also valid in a setting where iterative linearization with inexact solvers is considered. Numerical tests are conducted for nondegenerate and degenerate cases having exact solutions, as well as for a realistic case and a benchmark case. It is shown that the estimators correctly identify the errors up to a factor of the order of unity.

Distributed Computation of Stochastic GNE With Partial Information: An Augmented Best-Response Approach

In this paper, we focus on the stochastic generalized Nash equilibrium problem (SGNEP) which is an important and widely-used model in many different fields. In this model, subject to certain global resource constraints, a set of self-interested players aim to optimize their local objectives that depend on their own decisions and the decisions of others and are influenced by some random factors. We propose a distributed stochastic generalized Nash equilibrium seeking algorithm in a partial-decision information setting based on the Douglas-Rachford operator splitting scheme, which relaxes assumptions in the existing literature. The proposed algorithm updates players' local decisions through augmented best-response schemes and subsequent projections onto the local feasible sets, which comprise most of the computational workload. The projected stochastic subgradient method is applied to provide approximate solutions to the augmented best-response subproblems for each player. The Robbins-Siegmund theorem is leveraged to establish the main convergence results to a true Nash equilibrium using the proposed inexact solver. Finally, we illustrate the validity of the proposed algorithm via two numerical examples, i.e., a stochastic Nash-Cournot distribution game and a multi-product assembly problem with the two-stage model.

Two-point Landweber-type method with convex penalty terms for nonsmooth nonlinear inverse problems

Abstract In this work, we propose a two-point Landweber-type method with general convex penalty terms for solving nonsmooth nonlinear inverse problems. The design of our method can cope with nonsmooth nonlinear inverse problems or the nonlinear inverse problems whose data are contaminated by various types of noise. The method consists of the two-point acceleration strategy and inner solvers. Inner solvers are used to solve the minimization problems with respect to the penalty terms in each steps. If the minimization problems can be solved explicitly, the inner solvers will be chosen to be the exact solvers. Otherwise, we will use inexact solvers as inner solvers. Convergence results are given without utilizing the Gâteaux differentiability of the forward operator or the reflexivity of the image space. Numerical simulations are given to test the performance of the proposed method.

A recycling preconditioning method with auxiliary tip subspace for elastic crack propagation simulation using XFEM

In this paper, we propose a recycling preconditioning method with auxiliary tip subspace for solving a sequence of highly ill-conditioned linear systems of equations of different sizes arising from elastic crack propagation problems discretized by the extended finite element method. To construct a Schwarz type preconditioner, the finite element mesh is decomposed into crack tip subdomains, which contain all the degrees of freedom (DOFs) of the branch enrichment functions, and regular subdomains, which contain the standard DOFs and the DOFs of the Heaviside and the Junction enrichment functions. As cracks propagate these subdomains are modified accordingly, and the subdomain matrices are constructed as the restriction of the global matrix to the subdomains. In the overlapping Schwarz preconditioners, the crack tip subproblems are solved exactly and the regular subproblems are solved by some inexact solvers, such as ILU. We consider problems with and without crack intersections and develop a simple scheme to update, instead of re-computing, the subdomain problems as cracks propagate, in which only crack tip subdomains are updated around the new crack tips and all the regular subdomains remain unchanged. Therefore, no extra search is required, and the sizes of crack tip subproblems do not increase as cracks propagate, which greatly saves the computational cost. Moreover, starting from the second system, the Krylov subspace method uses a nontrivial initial guess constructed using the solution of the previous system with a modification around the new crack tips. The strategy accelerates the convergence remarkably. Numerical experiments demonstrate the efficiency of the proposed algorithms applied to problems with several types of cracks.

A domain decomposition method for Isogeometric multi-patch problems with inexact local solvers

In Isogeometric Analysis, the computational domain is often described as multi-patch, where each patch is given by a tensor product spline/NURBS parametrization. In this work we propose a FETI-like solver where local inexact solvers exploit the tensor product structure at the patch level. To this purpose, we extend to the isogeometric framework the so-called All-Floating variant of FETI, that allows us to use the Fast Diagonalization method at the patch level. We construct then a preconditioner for the whole system and prove its quasi-robustness with respect to the local mesh-size h and patch-size H: precisely the condition number of the preconditioned system is bounded by the square of the logarithm of H∕h. Our numerical tests confirm the theory and also show a favourable dependence of the computational cost of the method from the spline degree p.

Progressive Iterative Approximation with Preconditioners

The progressive iterative approximation (PIA) plays an important role in curve and surface fitting. By using the diagonally compensated reduction of the collocation matrix, we propose the preconditioned progressive iterative approximation (PPIA) to improve the convergence rate of PIA. For most of the normalized totally positive bases, we show that the presented PPIA can accelerate the convergence rate significantly in comparison with the weighted progressive iteration approximation (WPIA) and the progressive iterative approximation with different weights (DWPIA). Furthermore, we propose an inexact variant of the PPIA (IPPIA) to reduce the computational complexity of the PPIA. We introduce the inexact solver of the preconditioning system by employing some state-of-the-art iterative methods. Numerical results show that both the PPIA and the IPPIA converge faster than the WPIA and DWPIA, while the elapsed CPU times of the PPIA and IPPIA are less than those of the WPIA and DWPIA.

A high-precision low-dispersive nearly analytic difference method with its application in frequency-domain seismic waveform inversion

ABSTRACTTo improve the efficiency of seismic waveform inversion, we construct a novel high-order (6th-order) frequency-domain nearly analytic discrete method for forward modelling in inversion processes. Compared with some existing numerical schemes, this method is more powerful in suppressing numerical dispersion and enhancing the accuracy of computational results based on coarse discrete grids. We describe the discretisation of the frequency-domain wave equation, investigate the efficiency of wave-field simulation and perform numerical dispersion analysis. The inversion efficiency can also be improved when the inexact rotated block triangular preconditioned Krylov iteration solvers are incorporated into the corresponding linear system. Moreover, the inversion behaviour of frequency-domain NAD methods is also systematically discussed as one of its important applications. We take three representative media models as examples and the reliable inversion results are obtained reflecting the effectiveness of the proposed methods.

An adaptive [formula omitted]-refinement strategy with inexact solvers and computable guaranteed bound on the error reduction factor

In this work we extend our recently proposed adaptive refinement strategy for hp-finite element approximations of elliptic problems by taking into account an inexact algebraic solver. Namely, on each level of refinement and on each iteration of an (arbitrary) iterative algebraic solver, we compute guaranteed a posteriori error bounds on the algebraic and the total errors in energy norm. The algebraic error is the difference between the inexact discrete solution obtained by an iterative algebraic solver and the (unavailable) exact discrete solution. On the other hand, the total error stands for the difference between the inexact discrete solution and the (unavailable) exact solution of the partial differential equation. For the algebraic error upper bound, we crucially exploit the whole nested hierarchy of hp-finite element spaces created by the adaptive algorithm, whereas the remaining parts of the total error upper and lower bounds are computed using the finest space only. These error bounds allow us to formulate adaptive stopping criteria for the algebraic solver ensuring that the algebraic error does not significantly contribute to the total error. Next, we use the total error bound to mark mesh vertices for refinement via Dörfler’s bulk-chasing criterion. On patches associated with marked vertices only, we solve two separate primal finite element problems with homogeneous Dirichlet (Neumann) boundary conditions, which serve to decide between h-, p-, or hp-refinement. Altogether, we show that these ingredients lead to a computable guaranteed bound on the ratio of the total errors of the inexact approximations between successive refinements (the error reduction factor), when the stopping criteria are satisfied. Finally, in a series of numerical experiments, we investigate the practicality of the proposed adaptive solver, the accuracy of our bound on the reduction factor, and show that exponential convergence rates are achieved even in the presence of an inexact algebraic solver.

Goal-oriented a posteriori error estimation for conforming and nonconforming approximations with inexact solvers

We derive a unified framework for goal-oriented a posteriori estimation covering in particular higher-order conforming, nonconforming, and discontinuous Galerkin finite element methods, as well as the finite volume method. The considered problem is a model linear second-order elliptic equation with inhomogeneous Dirichlet and Neumann boundary conditions and the quantity of interest is given by an arbitrary functional composed of a volumetric weighted mean value (source) term and a surface weighted mean (Dirichlet boundary) flux term. We specifically do not request the primal and dual discrete problems to be resolved exactly, allowing for inexact solves. Our estimates are based on H(div)-conforming flux reconstructions and H1-conforming potential reconstructions and provide a guaranteed upper bound on the goal error. The overall estimator is split into components corresponding to the primal and dual discretization and algebraic errors, which are then used to prescribe efficient stopping criteria for the employed iterative algebraic solvers. Numerical experiments are performed for the finite volume method applied to the Darcy porous media flow problem in two and three space dimensions. They show excellent effectivity indices even in presence of primal and dual algebraic errors and enable to spare a large percentage of unnecessary algebraic iterations.

Localization of the W-1,q norm for local a posteriori efficiency

Abstract This paper gives a direct proof of localization of dual norms of bounded linear functionals on the Sobolev space ${W^{1,p}_0(\varOmega )}$, $1 \leq p \leq \infty $. The basic condition is that the functional in question vanishes over locally supported test functions from ${W^{1,p}_0(\varOmega )}$ which form a partition of unity in $\varOmega $, apart from close to the boundary $\partial \varOmega $. We also study how to weaken this condition. The results allow in particular to establish local efficiency and robustness with respect to the exponent $p$ of a posteriori estimates for nonlinear partial differential equations in divergence form, including the case of inexact solvers. Numerical illustrations support the theory.

Local Fourier analysis for mixed finite-element methods for the Stokes equations

In this paper, we develop a local Fourier analysis of multigrid methods based on block-structured relaxation schemes for stable and stabilized mixed finite-element discretizations of the Stokes equations, to analyze their convergence behavior. Three relaxation schemes are considered: distributive, Braess-Sarazin, and Uzawa relaxation. From this analysis, parameters that minimize the local Fourier analysis smoothing factor are proposed for the stabilized methods with distributive and Braess-Sarazin relaxation. Considering the failure of the local Fourier analysis smoothing factor in predicting the true two-grid convergence factor for the stable discretization, we numerically optimize the two-grid convergence predicted by local Fourier analysis in this case. We also compare the efficiency of the presented algorithms with variants using inexact solvers. Finally, some numerical experiments are presented to validate the two-grid and multigrid convergence factors.

Vertex-Based Preconditioners for the Coarse Problems of BDDC

We present a family of approximate BDDC preconditioners based on inexact solvers for the coarse problem. The basic idea is to replace the direct solver for a standard BDDC coarse problem by a preco...

Adaptive BEM with inexact PCG solver yields almost optimal computational costs

We consider the preconditioned conjugate gradient method (PCG) with optimal preconditioner in the frame of the boundary element method for elliptic first-kind integral equations. Our adaptive algorithm steers the termination of PCG as well as the local mesh-refinement. Besides convergence with optimal algebraic rates, we also prove almost optimal computational complexity. In particular, we provide an additive Schwarz preconditioner which can be computed in linear complexity and which is optimal in the sense that the condition numbers of the preconditioned systems are uniformly bounded. As model problem serves the 2D or 3D Laplace operator and the associated weakly-singular integral equation with energy space widetilde{H}^{-1/2}(Gamma ). The main results also hold for the hyper-singular integral equation with energy space H^{1/2}(Gamma ).

Nonlinear BDDC Methods with approximate solvers

New nonlinear BDDC (Balancing Domain Decomposition by Constraints) domain decomposition methods using inexact solvers for the subdomains and the coarse problem are proposed. In nonlinear domain decomposition methods, the nonlinear problem is decomposed before linearization to improve concurrency and robustness. For linear problems, the new methods are equivalent to known inexact BDDC methods. The new approaches are therefore discussed in the context of other known inexact BDDC methods for linear problems. Relations are pointed out, and the advantages of the approaches chosen here are highlighted. For the new approaches, using an algebraic multigrid method as a building block, parallel scalability is shown for more than half a million (524288) MPI ranks on the JUQUEEN IBM BG/Q supercomputer (JSC Jülich, Germany) and on up to 193600 cores of the Theta Xeon Phi supercomputer (ALCF, Argonne National Laboratory, USA), which is based on the recent Intel Knights Landing (KNL) many-core architecture. One of our nonlinear inexact BDDC domain decomposition methods is also applied to three-dimensional plasticity problems. Comparisons to standard Newton-Krylov-BDDC methods are provided.

Development of an efficient tightly coupled method for multiphysics reactor transient analysis

Picard Iteration is a widely used coupling method for multiphysics simulations. This method allows one to directly leverage existing and well-developed single-physics programs without re-writing large portions of the codes. In Picard Iteration, single-physics codes just iteratively pass solutions to each other as inputs until each code has reached a converged solution. However, multiphysics computation linked by Picard Iteration is susceptible to over-solving, which can make the overall computation much less efficient. Over-solving means that each single-physics code provides an accurate solution in each Picard Iteration, which is not necessary in practice. Solving the single-physics codes in an inexact manner, i.e. with relaxed termination criteria, can help avoid this problem. This work develops a modified Picard Iteration coupling method with adaptive, inexact termination criteria for the underlying single-physics codes. Also, nested within the inexact Picard Iteration, inexact Newton methods were applied in the single-physics codes. The effect on the overall computation efficiency due to the inexact (relaxed) termination criteria at both levels is investigated by applying them to solve reactor transient problems. A reactor dynamics problem with temperature feedback in one-dimensional slab geometry is used to scope the behavior of nested inexact solvers. Then these methods are applied to a larger two-dimensional Boiling Water Reactor (BWR) problem. Computational time savings reach 55% for the two-dimensional problem. Additionally, applying an inexact termination criterion (inexact Newton method) to each single-physics code results in a further time savings of up to 18%.

Rate optimal adaptive FEM with inexact solver for nonlinear operators

We prove convergence with optimal algebraic rates for an adaptive finite element method for nonlinear equations with strongly monotone operator. Unlike prior works, our analysis also includes the iterative and inexact solution of the arising nonlinear systems by means of the Picard iteration. Using nested iteration, we prove, in particular, that the number of of Picard iterations is uniformly bounded in generic cases, and the overall computational cost is (almost) optimal. Numerical experiments confirm the theoretical results.

Iterative minimization algorithm for efficient calculations of transition states

This paper presents an efficient algorithmic implementation of the iterative minimization formulation (IMF) for fast local search of transition state on potential energy surface. The IMF is a second order iterative scheme providing a general and rigorous description for the eigenvector-following (min-mode following) methodology. We offer a unified interpretation in numerics via the IMF for existing eigenvector-following methods, such as the gentlest ascent dynamics, the dimer method and many other variants. We then propose our new algorithm based on the IMF. The main feature of our algorithm is that the translation step is replaced by solving an optimization subproblem associated with an auxiliary objective function which is constructed from the min-mode information. We show that using an efficient scheme for the inexact solver and enforcing an adaptive stopping criterion for this subproblem, the overall computational cost will be effectively reduced and a super-linear rate between the accuracy and the computational cost can be achieved. A series of numerical tests demonstrate the significant improvement in the computational efficiency for the new algorithm.

On the scalability of inexact balancing domain decomposition by constraints with overlapped coarse/fine corrections

In this work, we analyze the scalability of inexact two-level balancing domain decomposition by constraints (BDDC) preconditioners for Krylov subspace iterative solvers, when using a highly scalable asynchronous parallel implementation where fine and coarse correction computations are overlapped in time. This way, the coarse-grid problem can be fully overlapped by fine-grid computations (which are embarrassingly parallel) in a wide range of cases. Further, we consider inexact solvers to reduce the computational cost/complexity and memory consumption of coarse and local problems and boost the scalability of the solver. Out of our numerical experimentation, we conclude that the BDDC preconditioner is quite insensitive to inexact solvers. In particular, one cycle of algebraic multigrid (AMG) is enough to attain algorithmic scalability. Further, the clear reduction of computing time and memory requirements of inexact solvers compared to sparse direct ones makes possible to scale far beyond state-of-the-art BDDC implementations. Excellent weak scalability results have been obtained with the proposed inexact/overlapped implementation of the two-level BDDC preconditioner, up to 93,312 cores and 20 billion unknowns on JUQUEEN. Further, we have also applied the proposed setting to unstructured meshes and partitions for the pressure Poisson solver in the backward-facing step benchmark domain.