Composing Scalable Nonlinear Algebraic Solvers

Most efficient linear solvers use composable algorithmic components, with the most common model being the combination of a Krylov accelerator and one or more preconditioners. A similar set of concepts may be used for nonlinear algebraic systems, where nonlinear composition of different nonlinear solvers may significantly improve the time to solution. We describe the basic concepts of nonlinear composition and preconditioning and present a number of solvers applicable to nonlinear partial differential equations. We have developed a software framework in order to easily explore the possible combinations of solvers. We show that the performance gains from using composed solvers can be substantial compared with gains from standard Newton-Krylov methods.Comment: 29 pages, 14 figures, 13 table

Metric based up-scaling

We consider divergence form elliptic operators in dimension $n\geq 2$ with $L^\infty$ coefficients. Although solutions of these operators are only H\"{o}lder continuous, we show that they are differentiable ($C^{1,\alpha}$) with respect to harmonic coordinates. It follows that numerical homogenization can be extended to situations where the medium has no ergodicity at small scales and is characterized by a continuum of scales by transferring a new metric in addition to traditional averaged (homogenized) quantities from subgrid scales into computational scales and error bounds can be given. This numerical homogenization method can also be used as a compression tool for differential operators.Comment: Final version. Accepted for publication in Communications on Pure and Applied Mathematics. Presented at CIMMS (March 2005), Socams 2005 (April), Oberwolfach, MPI Leipzig (May 2005), CIRM (July 2005). Higher resolution figures are available at http://www.acm.caltech.edu/~owhadi

Homogenization of Parabolic Equations with a Continuum of Space and Time Scales

This paper addresses the issue of the homogenization of linear divergence form parabolic operators in situations where no ergodicity and no scale separation in time or space are available. Namely, we consider divergence form linear parabolic operators in $\Omega \subset \mathbb{R}^n$ with $L^\infty(\Omega \times (0,T))$-coefficients. It appears that the inverse operator maps the unit ball of $L^2(\Omega\times (0,T))$ into a space of functions which at small (time and space) scales are close in $H^1$ norm to a functional space of dimension $n$. It follows that once one has solved these equations at least $n$ times it is possible to homogenize them both in space and in time, reducing the number of operation counts necessary to obtain further solutions. In practice we show under a Cordes-type condition that the first order time derivatives and second order space derivatives of the solution of these operators with respect to caloric coordinates are in $L^2$ (instead of $H^{-1}$ with Euclidean coordinates). If the medium is time-independent, then it is sufficient to solve $n$ times the associated elliptic equation in order to homogenize the parabolic equation

Preconditioners for Krylov subspace methods: An overview

When simulating a mechanism from science or engineering, or an industrial process, one is frequently required to construct a mathematical model, and then resolve this model numerically. If accurate numerical solutions are necessary or desirable, this can involve solving large-scale systems of equations. One major class of solution methods is that of preconditioned iterative methods, involving preconditioners which are computationally cheap to apply while also capturing information contained in the linear system. In this article, we give a short survey of the field of preconditioning. We introduce a range of preconditioners for partial differential equations, followed by optimization problems, before discussing preconditioners constructed with less standard objectives in mind

An analysis for some methods and algorithms of quantum chemistry

In der theoretischen Berechnung der Eigenschaften von Atomen, Molekülen und Festkörpern spielt die Lösung der elektronischen Schrödingergleichung, einer Operatoreigenwertgleichung für den Hamiltonoperator H des jeweiligen Systems, eine zentrale Rolle. Besondere Bedeutung kommt hierbei dem kleinsten Eigenwert von H zu, der die Grundzustandsenergie des Systems angibt. Um den unterschiedlichen Anforderungen in der Fülle von Anwendungsgebieten der elektronischen Schrödingergleichung gerecht zu werden, wurden in den letzten Jahrzehnten verschiedenste Näherungsverfahren entwickelt, um die Lösung dieses extrem hochdimensionalen Minimierungsproblems zu approximieren. Das Ziel der vorliegenden Arbeit ist es, eine (mathematische) Analysis für Aspekte einiger der verwendeten Methoden der Quantenchemie zu liefern. Zu diesem Zweck gliedert sich die Arbeit in vier Teile: Der erste Teil gibt eine Einführung in den mathematischen, hauptsächlich der Funktionalanalysis zuzuschreibenden Hintergrund, der bei der Behandlung der elektronischen Schrödingergleichung als Operatoreigenwertgleichung notwendig ist, und stellt viele der in den späteren Kapiteln benötigten Handwerkszeuge zur Verfügung. Der zweite Teil beschäftigt sich mit einem Gradientenalgorithmus mit Orthogonalitätsnebenbedingungen, der zur der Lösung der in der Beschreibung größerer Systeme wichtigen Hartree-Fock- und Kohn-Sham-Gleichungen, aber auch zur algorithmischen Behandlung der CI-Methode und außerhalb der Elektronenstrukturberechnung in der Berechnung invarianter Unterräume verwendet wird. Wir formulieren den Algorithmus allgemeiner als Verfahren zur Behandlung von Minimierungsproblemen auf der sogenannten Grassmann-Mannigfaltigkeit [1] und beweisen mit Hilfe dieses Formalismus unter anderem lineare Konvergenz des Algorithmus und die quadratische Konvergenz der zugeh√∂rigen Energien. Im dritten Teil der Arbeit wird die in der Praxis für hochgenaue Rechnungen bedeutsame Coupled-Cluster-Methode, traditionell ein Ansatz zur Approximation der Galerkinlösung der Schrödingergleichung innerhalb einer gegebenen Diskretisierung [2], als Verfahren im unendlichdimensionalen, undiskretisierten Raum formuliert. Zu diesem Zweck wird zunächst die Stetigkeit des Clusteroperators T als Operator vom Sobolevraum H1 in sich bewiesen: hieraus lässt sich dann die (unendlichdimensionale Verallgemeinerung der bekannten) Nullstellengleichung für die Coupled-Cluster-Funktion formulieren. Wir zeigen die lokale starke Monotonie der CC-Funktion, mit deren Hilfe wir Existenz- und Eindeutigkeitsaussagen und einen zielorientierten Fehlerschätzer nach [3] beweisen. Schließlich diskutieren wir die algorithmische Behandlung der oben genannten Nullstellengleichung. Teil 4 beschäftigt sich mit der DIIS-Methode, einem im Rahmen der Quantenchemie standardmäßig verwendeten Verfahren zur Konvergenzbeschleuningung iterativer Algorithmen. Wir identifizieren DIIS mit einer Variante des projezierten Broyden-Verfahrens [4] und zeigen, dass sich das Verfahren, angewandt auf lineare Probleme, als Variante des GMRES-Verfahrens auffassen lässt. Für den allgemeinen Fall beweisen wir schließlich zwei lokale Konvergenzaussagen und diskutieren die Umstände, unter denen DIIS superlineares Konvergenzverhalten zeigen kann. [1] T. A. Arias, A. Edelman, S. T. Smith, SIAM J. Matrix Anal. and Appl. 20, 2, 1999. [2] R. Schneider, Num. Math. 113, 3, 2009. [3] R. Becker, R. Rannacher, Acta Numerica 2000 (A. Iserlet, ed.), Cambridge University Press, 2001. [4] D. M. Gay, R. B. Schnabel, Nonlinear Programming 3, Academic Press, 1978.In the field of ab-initio calculation of the properties of atoms, molecules and solids, the solution of the electronic Schrödinger equation, an operator eigenvalue equation for the Hamiltonian of the system, plays a major role. Of utmost significance is the lowest eigenvalue of H, representing the ground state energy of the system. To meet the requirements of the multitude of possible applications of the elctronic Schrödinger equation, the last decades have seen the development of a variety of different methods designed to approximate the solution of this extremely high-dimensional minimization problem. The aim of the present work is to deliver a (mathematical) analysis for some aspects of some of these methods used in the context of quantum chemistry. The work consists of four parts: The first part gives an introduction to the mathematical background, mainly belonging to the field of functional analysis, that is needed for the rigogous treatment of the electronic Schrödinger equation as an operator eigenvalue equation, and provides many of the technical means needed in the later chapters. The second part is concerned with a gradient algorithm with orthogonality constraints, which is used for the solution Hartree-Fock and Kohn-Sham equations playing an important role in the description of larger systems and which also serves for the algorithmic treatment of the CI method and - outside of the field of electronic structure calculation - for the calculation of invariant subspaces. The algorithm is formulated as an abstract method for the treatment of minimization problems on the so-called Grassmann manifold [1]; with the help of this formalism, linear convergence of the algorithm and quadratic convergence of the corresponding eigenvalues is proven. The third part of the work is concerned with the Coupled Cluster method, being of high practical significance in calculations where high accuracy is demanded. We lift the method, usually formulated as an ansatz for the approximation of the Galerkin solution in a finite dimensional, discretised subspace [2] to the continuous, undiscretised space, resulting in what we will call the continuous Coupled Cluster method. To define the continuous method, we first prove the continuity of the cluster operator T as an operator mapping the Sobolev space H1 to itself; with the help of this result, the infinite dimensional globalization of the canonical) Coupled Cluster equations can be formulated. Afterwards, we prove local strong monotonicity of the CC function, from which we derive existence and (local) uniqueness statements and a goal-oriented a-posteriori error estimator in the fashion of [3]. Finally, we discuss the algorithmic treatment of the CC root equation. The last part of this work features an analysis for the acceleration technique DIIS that is commonly used in quantum chemistry codes. We identify DIIS with a variant of a projected Broyden's method [4] and show that when applied to linear systems, the method can be interpreted as a variant of the well-known GMRES method. For the global nonlinear case, we finally prove two local convergence results and discuss the circumstances under which DIIS can show superlinear convergence. [1] T. A. Arias, A. Edelman, S. T. Smith, SIAM J. Matrix Anal. and Appl. 20, 2, 1999. [2] R. Schneider, Num. Math. 113, 3, 2009. [3] R. Becker, R. Rannacher, Acta Numerica 2000 (A. Iserlet, ed.), Cambridge University Press, 2001. [4] D. M. Gay, R. B. Schnabel, Nonlinear Programming 3, Academic Press, 1978

An hp-adaptive discontinuous Galerkin method for phase field fracture

The phase field method is becoming the de facto choice for the numerical analysis of complex problems that involve multiple initiating, propagating, interacting, branching and merging fractures. However, within the context of finite element modelling, the method requires a fine mesh in regions where fractures will propagate, in order to capture sharp variations in the phase field representing the fractured/damaged regions. This means that the method can become computationally expensive when the fracture propagation paths are not known a priori. This paper presents a 2D -adaptive discontinuous Galerkin finite element method for phase field fracture that includes a posteriori error estimators for both the elasticity and phase field equations, which drive mesh adaptivity for static and propagating fractures. This combination means that it is possible to be reliably and efficiently solve phase field fracture problems with arbitrary initial meshes, irrespective of the initial geometry or loading conditions. This ability is demonstrated on several example problems, which are solved using a light-BFGS (Broyden–Fletcher–Goldfarb–Shanno) quasi-Newton algorithm. The examples highlight the importance of driving mesh adaptivity using both the elasticity and phase field errors for physically meaningful, yet computationally tractable, results. They also reveal the importance of including -refinement, which is typically not included in existing phase field literature. The above features provide a powerful and general tool for modelling fracture propagation with controlled errors and degree-of-freedom optimised meshes

International Conference on Continuous Optimization (ICCOPT) 2019 Conference Book

The Sixth International Conference on Continuous Optimization took place on the campus of the Technical University of Berlin, August 3-8, 2019. The ICCOPT is a flagship conference of the Mathematical Optimization Society (MOS), organized every three years. ICCOPT 2019 was hosted by the Weierstrass Institute for Applied Analysis and Stochastics (WIAS) Berlin. It included a Summer School and a Conference with a series of plenary and semi-plenary talks, organized and contributed sessions, and poster sessions. This book comprises the full conference program. It contains, in particular, the scientific program in survey style as well as with all details, and information on the social program, the venue, special meetings, and more

Computational modeling and simulation of nonlinear electromagnetic and multiphysics problems

In this dissertation, nonlinear electromagnetic and multiphysics problems are modeled and simulated using various three-dimensional full-wave methods in the time domain. The problems under consideration fall into two categories. One is nonlinear electromagnetic problems with the nonlinearity embedded in either the permeability or the conductivity of the material's constitutive properties. The other is multiphysics problems that involve interactions between electromagnetic and other physical phenomena. A numerical solution of nonlinear magnetic problems is formulated using the three-dimensional time-domain finite element method (TDFEM) combined with the inverse Jiles-Atherton vector hysteresis model. A second-order nonlinear partial differential equation (PDE) that governs the nonlinear magnetic problem is constructed through the magnetic vector potential in the time domain, which is solved by applying the Newton-Raphson method. To solve the ordinary differential equation (ODE) representing the magnetic hysteresis accurately and efficiently, several ODE solvers are specifically designed and investigated. To improve the computational efficiency of the Newton-Raphson method, the multi-dimensional secant methods are incorporated in the nonlinear TDFEM solver. A nonuniform time-stepping scheme is also developed using the weighted residual approach to remove the requirement of a uniform time-step size during the simulation. Breakdown phenomena during high-power microwave (HPM) operation are investigated using different physical and mathematical models. During the breakdown process, the bound charges in solid dielectrics and air molecules break free and are pushed to move by the Lorentz force produced by the electromagnetic fields. The motion of free electrons produces plasma currents, which generate secondary electromagnetic fields that couple back to the externally applied fields and interact with the free electrons. When the incident field intensity is high enough, this will lead to an exponential increase of the charged particles known as breakdown. Such a process is first described by a nonlinear conductivity of the solid dielectric as a function of the electric field to model the dielectric breakdown phenomenon. The air breakdown problem encountered with HPM operation is then simulated with the plasma current modeled by a simplified plasma fluid equation. Both the dielectric and air breakdown problems are solved with the TDFEM together with a Newton's method, where the dielectric breakdown is treated as a pure nonlinear electromagnetic problem, while the air breakdown is treated as a multiphysics problem. To describe the plasma behavior more accurately, the plasma density and velocity are modeled by the equations of diffusion and motion, respectively. This results in a multiphysics and multiscale system depicted by the nonlinearly coupled full-wave Maxwell and plasma fluid equations, which are solved by a nodal discontinuous Galerkin time-domain (DGTD) method in three dimensions. The air breakdown during the HPM operation and the resulting plasma formation and shielding are modeled and simulated. Several important numerical issues in the simulation of nonlinear electromagnetic and multiphysics problems have been investigated and discussed. A continuity-preserving and divergence-cleaning scheme for electromagnetic problems involving inhomogeneous materials has been proposed based on the purely and damped hyperbolic Maxwell equations. A divergence-cleaning method is presented to enforce Gauss's laws and normal flux continuity by introducing auxiliary variables and damping terms into the original Maxwell's equations, which result in artificial propagation and dissipation of the numerical errors. Based on the DGTD method, dynamic h- and p-adaptation algorithms are developed for a full-wave analysis of electromagnetic and multiphysics problems. The dynamic h-adaptation algorithm can dynamically refine the mesh to resolve the local variation of the fields during the wave propagation, while the dynamic p-adaptation algorithm can determine and adjust the basis order in real time during the simulation. Both algorithms developed and investigated in this dissertation are highly flexible and efficient, and are powerful simulation tools in the solution of nonlinear electromagnetic and multiphysics problems