Low-rank update of preconditioners for the nonlinear Richard's equation

Preconditioners for the Conjugate Gradient method are studied to solve the Newton system with symmetric positive definite (SPD) Jacobian. Following the theoretical work in Bergamaschi et al. (2011) [4] we start from a given approximation of the inverse of the initial Jacobian, and we construct a sequence of preconditioners by means of a low rank update, for the linearized systems arising in the Picard Newton solution of the nonlinear discretized Richards equation. Numerical results onto a very large and realistic test case show that the proposed approach is more efficient, in terms of iteration number and CPU time, as compared to computing the preconditioner of choice at every nonlinear iteration.The support of the CARIPARO Foundation (Grant NPDE: Non-linear Partial Differential Equations: models, analysis, and control - theoretic problems), and of the Spanish DGI grant MTM2010-18674 is acknowledged.Bergamaschi, L.; Bru García, R.; Martínez Calomardo, Á.; Mas Marí, J.; Putti, M. (2013). Low-rank update of preconditioners for the nonlinear Richard's equation. Mathematical and Computer Modelling. 57(7):1933-1941. https://doi.org/10.1016/j.mcm.2012.01.013S1933194157

Low-rank updates of balanced incomplete factorization preconditioners

[EN] Let Ax = b be a large and sparse system of linear equations where A is a nonsingular matrix. An approximate solution is frequently obtained by applying preconditioned terations. Consider the matrix B = A + PQT where P,Q ∈ Rn×k are full rank matrices. In this work, we study the problem of updating a previously computed preconditioner for A in order to solve the updated linear system Bx = b by preconditioned iterations. In particular, we propose a method for updating a Balanced Incomplete Factorization preconditioner. The strategy is based on the computation of an approximate Inverse Sherman-Morrison decomposition for an equivalent augmented linear system. Approximation properties of the preconditioned matrix and an analysis of the computational cost of the algorithm are studied. Moreover the results of the numerical experiments with different types of problems show that the proposed method contributes to accelerate the convergence.This work was supported by the Spanish Ministerio de Economia y Competitividad under grant MTM2014-58159-P.Cerdán Soriano, JM.; Marín Mateos-Aparicio, J.; Mas Marí, J. (2017). Low-rank updates of balanced incomplete factorization preconditioners. Numerical Algorithms. 74(2):337-370. https://doi.org/10.1007/s11075-016-0151-6S337370742Bellavia, S., Bertaccini, D., Morini, B.: Nonsymmetric preconditioner updates in Newton-Krylov methods for nonlinear systems. SIAM J. Sci. Comput. 33 (5), 2595–2619 (2011)Benzi, M., Bertaccini, D.: Approximate inverse preconditioning for shifted linear systems. BIT 43(2), 231–244 (2003)Bergamaschi, L., Bru, R., Martínez, A.: Low-rank update of preconditioners for the inexact Newton method with SPD Jacobian. Math. Comput. Model. 54, 1863–1873 (2011)Bergamaschi, L., Bru, R., Martínez, A., Mas, J., Putti, M.: Low-rank update of preconditioners for the nonlinear Richards Equation. Math. Comput. Model. 57, 1933–1941 (2013)Bergamaschi, L., Gondzio, J., Venturin, M., Zilli, G.: Inexact constraint preconditioners for linear systems arising in interior point methods. Comput. Optim. Appl. 36(2-3), 137–147 (2007)Beroiz, M., Hagstrom, T., Lau, S.R., Price, R.H.: Multidomain, sparse, spectral-tau method for helically symmetric flow. Comput. Fluids 102(0), 250–265 (2014)Bertaccini, D.: Efficient preconditioning for sequences of parametric complex symmetric linear systems. Electron. Trans. Numer. Anal. 18, 49–64 (2004)Bollhöfer, M.: A robust and efficient ILU that incorporates the growth of the inverse triangular factors. SIAM J. Sci. Comput. 25(1), 86–103 (2003)Bollhöfer, M., Saad, Y.: On the relations between ILUs and factored approximate inverses. SIAM. J. Matrix Anal. Appl. 24(1), 219–237 (2002)Bru, R., Cerdán, J., Marín, J., Mas, J.: Preconditioning sparse nonsymmetric linear systems with the Sherman-Morrison formula. SIAM J. Sci. Comput. 25(2), 701–715 (2003)Bru, R., Marín, J., Mas, J., Tůma, M.: Balanced incomplete factorization. SIAM J. Sci. Comput. 30(5), 2302–2318 (2008)Bru, R., Marín, J., Mas, J., Tůma, M.: Improved balanced incomplete factorization. SIAM J. Matrix Anal. Appl. 31(5), 2431–2452 (2010)Cerdán, J., Faraj, T., Malla, N., Marín, J., Mas, J.: Block approximate inverse preconditioners for sparse nonsymmetric linear systems. Electron. Trans. Numer. Anal. 37, 23–40 (2010)Cerdán, J., Marín, J., Mas, J., Tůma, M.: Block balanced incomplete factorization. Technical Report No. TR-IMM2015/04, Polytechnic University of Valencia, Spain (2015)Davis, T.A.: University of Florida Sparse Matrix Collection. available online at http://www.cise.ufl.edu/~davis/sparse/ , NA Digest, vol. 94, issue 42, October 1994.Tebbens, J.D., Tůma, M.: Efficient preconditioning of sequences of nonsymmetric linear systems. SIAM J. Sci Comput. 29(5), 1918–1941 (2007)Tebbens, J.D., Tůma, M.: Preconditioner updates for solving sequences of linear systems in matrix-free environment. Numer Linear Algebra Appl. 17, 997–1019 (2010)Embree, M., Sifuentes, J.A., Soodhalter, K.M., Szyld, D.B., Xue, F.: Short-term recurrence Krylov subspace methods for nearly hermitian matrices. SIAM.J. Matrix Anal. Appl. 33-2, 480–500 (2012)Engquist, B., Ying, L.: Sweeping preconditioner for the Helmholtz equation: hierarchical matrix representation. Commun. Pure Appl. Math. 64, 697–735 (2011)Gatto, P., Christiansen, R.E., Hesthaven, J.S.: A preconditioner based on a low-rank approximation with applications to topology optimization. Technical Report EPFL-ARTICLE-207108, École polytechnique fédérale de Lausanne, EPFL, CH-1015 Lausanne, 2015.Grasedyck, L., Hackbusch, W.: Construction and arithmetics of H-matrices. Computing 70(4), 295–334 (2003)Grasedyck, L., Kressner, D., Tobler, C.: A literature survey of low-rank tensor approximation techniques. GAMM-Mitteilungen 36(1), 53–78 (2013)Greengard, L., Rokhlin, V.: A new version of the fast multipole method for the Laplace equation in three dimensions. Acta Numer. 6(1), 229–269 (1997)Hager, W.W.: Updating the inverse of matrix. SIAM Rev. 31(2), 221–239 (1989)Halko, N., Martinsson, P.G., Tropp, J.A.: Finding structure with randomness: probabilistic algorithms for constructing approximate matrix decompositions. SIAM Rev. 53(2), 217–288 (2011)Kelley, C.T.: Solving Nonlinear Equations with Newton’s Method. Fundamentals of algorithms. Society for Industrial and Applied Mathematics (SIAM), Philadelphia (2003)Saad, Y.: ILUT: a dual threshold incomplete L U factorization. Numer. Linear Algebra Appl. 1(4), 387–402 (1994)Saad, Y., Schulz, M.H.: GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J. Sci. Stat. Comput. 7, 856–869 (1986)Simoncini, V., Szyld, D.B.: The effect of non-optimal bases on the convergence of Krylov subspace methods. Numer Math. 100(4), 711–733 (2005)van der Vorst, H.A.: Bi-CGSTAB: a fast and smoothly converging variant of bi-CG for the solution of non-symmetric linear systems. SIAM J. Sci. Stat. Comput. 12, 631–644 (1992

Multiple local minima of PDE-constrained optimisation problems via deflation

Nonconvex optimisation problems constrained by partial differential equations (PDEs) may permit distinct local minima. In this paper we present a numerical technique, called deflation, for computing multiple local solutions of such optimisation problems. The basic approach is to apply a nonlinear transformation to the Karush-Kuhn-Tucker optimality conditions that eliminates previously found solutions from consideration. Starting from some initial guess, Newton's method is used to find a stationary point of the Lagrangian; this solution is then deflated away, and Newton's method is initialised from the same initial guess to find other solutions. In this paper, we investigate how the Schur complement preconditioners widely used in PDE-constrained optimisation perform after deflation. We prove an upper bound on the number of new distinct eigenvalues of a matrix after an arbitrary additive perturbation; from this it follows that for diagonalisable operators the number of Krylov iterations required for exact convergence of the Newton step at most doubles compared to the undeflated problem. While deflation is not guaranteed to converge to all minima, these results indicate the approach scales to arbitrary-dimensional problems if a scalable Schur complement pre-conditioner is available. The technique is demonstrated on a discretised nonconvex PDE-constrained optimisation problem with approximately ten million degrees of freedom

Scalable parallel simulation of variably saturated flow

In this thesis we develop highly accurate simulation tools for variably saturated flow through porous media able to take advantage of the latest supercomputing resources. Hence, we aim for parallel scalability to very large compute resources of over 105 CPU cores. Our starting point is the parallel subsurface flow simulator ParFlow. This library is of widespread use in the hydrology community and known to have excellent parallel scalability up to 16k processes. We first investigate the numerical tools this library implements in order to perform the simulations it was designed for. ParFlow solves the governing equation for subsurface flow with a cell centered finite difference (FD) method. The code targets high performance computing (HPC) systems by means of distributed memory parallelism. We propose to reorganize ParFlow's mesh subsystem by using fast partitioning algorithms provided by the parallel adaptive mesh refinement (AMR) library p4est. We realize this in a minimally invasive manner by modifying selected parts of the code to reinterpret the existing mesh data structures. Furthermore, we evaluate the scaling performance of the modified version of ParFlow, demonstrating excellent weak and strong scaling up to 458k cores of the Juqueen supercomputer at the Jülich Supercomputing Centre. The above mentioned results were obtained for uniform meshes and hence without explicitly exploiting the AMR capabilities of the p4est library. A natural extension of our work is to activate such functionality and make ParFlow a true AMR application. Enabling ParFlow to use AMR is challenging for several reasons: It may be based on assumptions on the parallel partition that cannot be maintained with AMR, it may use mesh-related metadata that is replicated on all CPUs, and it may assume uniform meshes in the construction of mathematical operators. Additionally, the use of locally refined meshes will certainly change the spectral properties of these operators. In this work, we develop an algorithmic approach to activate the usage of locally refined grids in ParFlow. AMR allows meshes where elements of different size neighbor each other. In this case, ParFlow may incur erroneous results when it attempts to communicate data between inter-element boundaries. We propose and discuss two solutions to this issue operating at two different levels: The first manipulates the indices of the degrees of freedom, While the second operates directly on the degrees of freedom. Both approaches aim to introduce minimal changes to the original ParFlow code. In an AMR framework, the FD method taken by ParFlow will require modifications to correctly deal with different size elements. Mixed finite elements (MFE) are on the other hand better suited for the usage of AMR. It is known that the cell centered FD method used in ParFlow might be reinterpreted as a MFE discretization using Raviart-Thomas elements of lower order. We conclude this thesis presenting a block preconditioner for saddle point problems arising from a MFE on locally refined meshes. We evaluate its robustness with respect to various classes of coefficients for uniform and locally refined meshes

A primal-dual active-set method for non-negativity constrained total variation deblurring problems

Master'sMASTER OF SCIENC

The workshop on iterative methods for large scale nonlinear problems

The aim of the workshop was to bring together researchers working on large scale applications with numerical specialists of various kinds. Applications that were addressed included reactive flows (combustion and other chemically reacting flows, tokamak modeling), porous media flows, cardiac modeling, chemical vapor deposition, image restoration, macromolecular modeling, and population dynamics. Numerical areas included Newton iterative (truncated Newton) methods, Krylov subspace methods, domain decomposition and other preconditioning methods, large scale optimization and optimal control, and parallel implementations and software. This report offers a brief summary of workshop activities and information about the participants. Interested readers are encouraged to look into an online proceedings available at http://www.usi.utah.edu/logan.proceedings. In this, the material offered here is augmented with hypertext abstracts that include links to locations such as speakers` home pages, PostScript copies of talks and papers, cross-references to related talks, and other information about topics addresses at the workshop

Spatio-temporal integral equation methods with applications

Electromagnetic interactions are vital in many applications including physics, chemistry, material sciences and so on. Thus, a central problem in physical modeling is the electromagnetic analysis of materials. Here, we consider the numerical solution of the Maxwell equation for the evolution of the electromagnetic field given the charges, and the Newton or Schr\\"odinger equation for the evolution of particles. By combining integral equation techniques with new spectral deferred correction algorithms in time and hierarchical methods in space, we develop fast solvers for the calculation of electromagnetism with relaxations of the model in different scenarios. The dissertation consists of two parts, aiming to resolve the challenges in the temporal and spatial direction, respectively. In the first part, we study a new class of time stepping methods for time-dependent differential equations. The core algorithm uses the pseudo-spectral collocation formulation to discretize the Picard type integral equation reformulation, producing a highly accurate and stable representation, which is then solved via the deferred correction technique. By exploiting the mathematical properties of the formulation and the convergence procedure, we develop some new preconditioning techniques from different perspectives that are accurate, robust, and can be much more efficient than existing methods. As is typical of spectral methods, the solution to the discretization is spectral accurate and the time step-size is optimal, though the cost of solving the system can be high. Thus, the solver is particularly suited to problems where very accurate solutions are sought or large time-step is required, e.g., chaotic systems or long-time simulation. In the second part, we study the hierarchical methods with emphasis on the spatial integral equations. In the first application, we implement a parallel version of the adaptive recursive solver for two-point boundary value problem by Cilk multithreaded runtime system based on the integral equation formulation. In the second application, we apply the hierarchical method to two-layered media Helmholtz equations in the acoustic and electromagnetic scattering problems. With the method of images and integral representations, the spatially heterogeneous translation operators are derived with rigorous error analysis, and the information is then compressed and spread in a fashion similar to fast multipole methods. The preliminary results suggest that our approach can be faster than existing algorithms with several orders of magnitude. We demonstrate our solver on a number of examples and discuss various useful extensions. Preliminary results are favorable and show the viability of our techniques for integral equations. Such integral equation methods could well have a broad impact on many areas of computational science and engineering. We describe further applications in biology, chemistry, and physics, and outline some directions for future work.Doctor of Philosoph

Fast prediction of transonic aeroelasticity using computational fluid dynamics

The exploitation of computational fluid dynamics for non linear aeroelastic simulations is mainly based on time domain simulations of the Euler and Navier-Stokes equations coupled with structural models. Current industrial practice relies heavily on linear methods which can lead to conservative design and flight envelope restrictions. The significant aeroelastic effects caused by nonlinear aerodynamics include the transonic flutter dip and limit cycle oscillations. An intensive research effort is underway to account for aerodynamic nonlinearity at a practical computational cost.To achieve this a large reduction in the numbers of degrees of freedoms is required and leads to the construction of reduced order models which provide compared with CFD simulations an accurate description of the dynamical system at much lower cost. In this thesis we consider limit cycle oscillations as local bifurcations of equilibria which are associated with degenerate behaviour of a system of linearised aeroelastic equations. This extra information can be used to formulate a method for the augmented solve of the onset point of instability - the flutter point. This method contains all the fidelity of the original aeroelastic equations at much lower cost as the stability calculation has been reduced from multiple unsteady computations to a single steady state one. Once the flutter point has been found, the centre manifold theory is used to reduce the full order system to two degrees of freedom. The thesis describes three methods for finding stability boundaries, the calculation of a reduced order models for damping and for limit cycle oscillations predictions. Results are shown for aerofoils, and the AGARD, Goland, and a supercritical transport wing. It is shown that the methods presented allow results comparable to the full order system predictions to be obtained with CPU time reductions of between one and three orders of magnitude