Continuation and Bifurcation Methods Using LOCA Eric Phipps Andy Salinger, Roger Pawlowski Sandia National Laboratories Trilinos Workshop at Copper Mountain March 30, 2004 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.
Why Do We Need Stability Analysis Algorithms for Large-Scale Applications? Nonlinear systems exhibit instabilities, e.g.: • Multiple steady states • Ignition • Symmetry Breaking • Onset of Oscillations • Phase Transitions LOCA: Library of Continuation Algorithms We need algorithms, software, and experience to impact ASCI- and SciDAC-sized applications. These phenomena must be understood in order to perform computational design and optimization. Established stability/bifurcation analysis libraries exist: • AUTO (Doedel) • CONTENT (Kuznetsov) • MATCONT (Govaerts) Stability/bifurcation analysis provides qualitative information about time evolution of nonlinear systems by computing families of steady-state solutions.
LOCA library grew out of continuation code in MPSalsa Andy Salinger, John Shadid, Roger Pawlowski, Louis Romero, Rich Lehoucq, Ed Wilkes, Beth Burroughs, Nawaf Bou-Rabee LOCA 1.0 released April 2002 Written in C with wrapper functions for linking to application code ~200 downloads Complete rewrite in C++ around NOX framework began September 2002, part of Trilinos release September 2003. History
r Tmax 1 3 1 Reaction Rate, r Second parameter, h LOCA: Library or Continuation Algorithms 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
Pseudo Arc-length Continuation Solves for Solution and Parameter Simultaneously
Codimension 1 Bifurcations Turning Point • Combustion • Buckling of an Arch • Buckling of a Beam • Pattern formation • Cell differentiation (morphogenesis) • Vortex Shedding • Predator-Prey models • Puberty Pitchfork Hopf
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
Pseudo Arc-length Continuation Bordering Algorithm Full Newton Algorithm Bordering Algorithms Meet these Requirements
Bordering Algorithms Meet these Requirements Full Newton Algorithm Turning Point Bifurcation … but 4 solves of per Newton Iteration are used to drive singular! Bordering Algorithm
Given initial guess , step size Solve nonlinear equations to find 1st point on curve while !stop Compute predictor Compute predicted point Solve continuation equations for using as initial guess If successful Postprocess (e.g., compute eigenvalues, output data) Increase step size Else Decrease step size Restore previous solution End if If or or stop = true End while Abstraction of Continuation Process LOCA Stepper Predictor modules NOX + continuation/ bifurcation groups Step size modules
NOX Nonlinear Solver (Kolda, Pawlowski, Hooper, Shadid) NOX implements various methods for solving Code to evaluate is encapsulated in a Group. NOX solver methods are generic, and implemented in terms of group/vector abstract interfaces: NOX solvers will work with any group/vector that implements these interfaces.
Super Vectors and Super Groups Idea: Given a vector to store and a group representing the equations , build an extended (“super”) group representing, e.g., pseudo arc-length continuation equations: and a super vector to store the solution component and parameter component . Super groups/vectors are generic: All abstract group/vector methods for super groups/vectors implemented in terms of methods of the underlying groups/vectors. Super groups are NOX groups: Extended nonlinear equations solved by most NOX solvers
Continuation Groups NOX::Abstract::Group LOCA::Continuation::ExtendedGroup LOCA::Continuation::NaturalGroup LOCA::Continuation::ArclengthGroup NOX::Abstract::Group LOCA::Continuation::AbstractGroup • setParam() • getParam() • operator = () • computeDfDp() • computeEigenvalues() • printSolution() Mandatory Default implementation available Optional Concrete group
Interfacing Application Codes to LOCA v2.0 • Interfacing NOX to the application code is 90% of the work! For • continuation • turning point tracking • pitchfork tracking at very minimum must be able to additionally set/retrieve parameter values, save complete state of system by copying group. • For Hopf tracking, must implement a complex solve: • Can overload many additional methods if better techniques are available • block solves • singular matrix solves • estimating derivatives:
Single parameter continuation Natural Pseudo Arc-length Bifurcations Turning point Pitchfork Hopf Predictors Constant Tangent Secant Random Step size control Constant Adaptive Artificial Homotopy Generic interface to Anasazi Native support for LAPACK (all intefaces) Epetra (all except Hopf) LOCA’s Current Capabilities
Continuation Example: Chan Problem ChanProblemInterface.H ChanProblemInterface.C ChanContinuation.C ChanContinuation.txt
Turning Point Continuation Example ChanTPContinuation.C ChanTPContinuation.txt
Structural Mechanics Example: Salinas Bending a 1D Beam Example problem from Salinas test suite Original continuation run with 50 load steps NOX/LOCA interface written by Russell Hooper Variable step size algorithm reduced to 37 load steps
3D Rayleigh-Bénard Problem in 5x5x1 box(208K unknowns, 16 processors) MPSalsa (Shadid et al., SNL): • Incompressible Navier-Stokes • Heat and Mass Transfer, Reactions • Unstrucured Finite Element (Galerkin/Least-Squares)
Improve robustness Better step size control Improved bifurcation tracking algorithms Debugging More features for homotopy Incorporate Multi-vector support Multi-parameter continuation (Henderson, IBM) Constraint enforcement Automatic differentiation Future Development