Computing Hopf Bifurcations in Large-Scale Problems

1. Computing Hopf Bifurcations in Large-Scale Problems Eric Phipps and Andy Salinger Applied Computational Methods Department Sandia National Laboratories 2007 SIAM Conference on Computational Science and Engineering February 23, 2007 Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company,for the United States Department of Energy’s National Nuclear Security Administration under contract DE-AC04-94AL85000.

2. Why Do We Need Stability Analysis Algorithms for Large-Scale Applications? • Nonlinear systems exhibit instabilities, e.g.: • Buckling • Ignition • Onset of Oscillations • Phase Transitions LOCA: Library of Continuation Algorithms We need algorithms, software, and experience to impact ASC- and SciDAC-sized applications (millions of unknowns) These phenomena must be understood in order to perform computational design and optimization. • Established stability/bifurcation analysis libraries exist: • AUTO (Doedel, et al) • CONTENT (Kuznetsov, et al) • MATCONT (Govaerts, et al) • etc… Stability/bifurcation analysis provides qualitative information about time evolution of nonlinear systems by computing families of steady-state solutions.

3. LOCA: Library of Continuation Algorithms • Application code provides: • Nonlinear steady-state residual and Jacobian fill: • Newton-like linear solves: • LOCA provides: • Parameter Continuation: Tracks a family of steady state solutions with parameter • Linear Stability Analysis: Calculates leading eigenvalues via Anasazi (Thornquist, Lehoucq) • Bifurcation Tracking: Locates neutral stability point (x,p) and tracks as a function of a second parameter External force 1 External force 3 1 Second parameter

4. Hopf Bifurcations ODE/DAE Eigenvalue passes through imaginary axis with nonzero imaginary part Stable equilibrium loses stability to a periodic motion Real-equivalent “Moore-Spence” formulation Examples • Vortex Shedding • Predator-Prey models • Flutter • El Niño

5. LOCA Designed for Easy Linking to Existing Newton-based Applications • LOCA targets existing codes that are: • Steady-State, Nonlinear • Newton’s Method • Large-Scale, Parallel • Algorithmic choices for LOCA: • Must work with iterative (approximate) linear solvers on distributed memory machines • Non-Invasive Implementation (e.g. matrix blind) • Should avoid or limit: • Requiring more derivatives • Changing sparsity pattern of matrix • Increasing memory requirements

6. Bordering Algorithms Meet these Requirements Pseudo-Arclength Continuation Bordering Algorithm Full Newton Algorithm

7. Bordering Algorithms Meet These Requirements (“Moore-Spence” Formulation) Hopf Bifurcation Full Newton Algorithm … but complex matrix is singular at the Hopf point! Bordering Algorithm Need to form “complex” matrix

8. LOCA has been rewritten as part of the Trilinos framework • Collection of parallel solver packages, with a common build process, made interoperable where appropriate (Heroux et al) • Better algorithms can be implemented with tighter coupling to linear algebra Trilinos 7.0 Publicly Available • Sandia’s parallel iterative solvers are made available to LOCA users 2004 • LOCA can build upon the existing interface between application codes and the NOX nonlinear solver HPC Software Challenge 2004

9. Bordered Systems of Equations Solving bordered systems of equations is a ubiquitous computation: Pseudo-Arclength Continuation Bordering Algorithm Constraint Following Bifurcation Identification • Only requires solves of J but • Requires m+l linear solves • Has difficulty when J is nearly singular

10. Solving Bordered Systems via QR Extension of Householder pseudo-arclength technique by Homer Walker1 Rearranged Bordered System Compact WY Representation2 QR Factorization where then Write P is nxn, nonsingular, rank m update to J 1H.F Walker, SIAM J. Sci. Comput., 1999 2R. Schreiber, SIAM J. Stat. Comput., 1989

11. Minimally Augmented Hopf Formulation Let At a Hopf point, is singular. Given and , let then is singular if and only if Standard formulation: Note for Newton’s method: 3 linear solves per Newton iteration (5 for Moore-Spence)! For symmetric problems reduces to 2 solves.

12. Minimally Augmented Hopf Formulation for Large-Scale Problems With iterative solvers, so define where Also

13. Real Equivalent Minimally Augmented Hopf Formulation Real-equivalent formulation to avoid complex arithmetic:

14. Vortex Shedding(5K 2D quad elements, 16K unknowns, 16 procs) • Charon (Hennigan, Pawlowski et al., SNL): • Incompressible Navier-Stokes • Heat & mass transfer, reactions • First order, unstructured finite element, pressure stabilization • Analytic, sparse Jacobian (Sacado) • Fully coupled Newton method (NOX) • GMRES (AztecOO) with RILU(k) preconditioner (Ifpack) • Distributed memory parallelism

15. Hopf Bifurcation Near Re = 48.9 • RILU fill factor: 4 • RILU overlap: 2 • Krylov space: 400 • F = Hopf residual (different for Moore-Spence and Minimally Augmented)

16. Pseudo-Arclength Equations Newton Solve Moore-Spence w/Bordering: 7 total linear solves Minimally Augmented: 3 total linear solves Pseudo-Arclength Continuation of Hopfs

17. Vortex SheddingPseudo-arclength Hopf Continuation

18. Summary • Bordered matrix methods provide powerful bifurcation capabilities • Turning points, pitchforks, and Hopfs • No singular matrix solves • Improves robustness, scalability and accuracy • Requires linear algebra-specific implementation • QR approach provides a convenient way to solve bordered systems • Nonsingular • Only involves one linear solve • Only requires simple vector operations • Doesn’t change dimension of the linear system • Future work • Understanding why bordered matrix methods work • Preconditioners for

19. Points of Contact • LOCA: Trilinos continuation and bifurcation package • Sub-package of NOX • Andy Salinger (agsalin@sandia.gov), Eric Phipps (etphipp@sandia.gov) • www.software.sandia.gov/nox • NOX: Trilinos nonlinear solver package • Roger Pawlowski (rppawlo@sandia.gov), Tammy Kolda (tgkolda@sandia.gov) • www.software.sandia.gov/nox • Trilinos: Collection of large-scale linear/nonlinear solvers • Epetra, Ifpack, AztecOO, ML, NOX, LOCA, … • Mike Heroux (maherou@sandia.gov) • www.software.sandia.gov/Trilinos • Trilinos Release 7.0 currently available, includes min. aug. turning point • 8.0 available this fall will include min. aug. pitchfork and Hopf • Charon: Massively Parallel Device Simulator/Reacting Fluid Flows • Rob Hoekstra (rjhoeks@sandia.gov), Roger Pawlowski (rppawlo@sandia.gov)