Continuation and Bifurcation Methods Using LOCA

1. Continuation and Bifurcation Methods Using LOCA Eric Phipps Andy Salinger, Roger Pawlowski 9233 – Computational Sciences Trilinos User Group Meeting October 16, 2003 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.: • 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.

3. 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

4. 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 • Bifurcation Tracking: Locates neutral stability point (x,p) and tracks as a function of a second parameter

5. Pseudo Arc-length Continuation Solves for Solution and Parameter Simultaneously

6. 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

7. 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

8. Pseudo Arc-length Continuation Bordering Algorithm Full Newton Algorithm Bordering Algorithms Meet these Requirements

9. 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

10. 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

11. 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.

12. 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

13. Arc-length Group applyJacobianInverse() LOCA::Continuation::ArclengthGroup::applyJacobianInverse(constNOX::Abstract::Vector& input,NOX::Abstract::Vector& result)const{ constLOCA::Continuation::ExtendedVector& con_input = dynamic_cast<constLOCA::Continuation::ExtendedVector&>(input); LOCA::Continuation::ExtendedVector& con_result = dynamic_cast<LOCA::Continuation::ExtendedVector&>(result); constNOX::Abstract::Vector& input_x = con_input.getXVec(); double input_p = con_input.getParam(); NOX::Abstract::Vector& result_x = con_result.getXVec(); double& result_p = con_result.getParam(); NOX::Abstract::Vector* b = input_x.clone(NOX::ShapeCopy); underlyingGroupPtr->applyJacobianInverse(input_x, result_x); underlyingGroupPtr->applyJacobianInverse(*dfdpVecPtr, *b); result_p = (predictorVecPtr->getXVec().dot(result_x) – input_p) / (predictorVecPtr->getXVec().dot(*b) – predictorVecPtr->getParam()); result_x.update(-result_p, *b, 1.0); delete b; }

14. 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

15. LOCA::Bifurcation::PitchforkBord::ExtendedGroup Turning Point, Pitchfork Groups NOX::Abstract::Group NOX::Abstract::Group LOCA::Continuation::AbstractGroup LOCA::Continuation::AbstractGroup LOCA::Bifurcation::TPBord::ExtendedGroup LOCA::Bifurcation::TPBord::AbstractGroup • computeDJnDp() • computeDJnDxa() • applySingularJacobianInverse() Concrete group

16. Hopf Groups NOX::Abstract::Group NOX::Abstract::Group LOCA::Continuation::AbstractGroup LOCA::Continuation::AbstractGroup LOCA::Bifurcation::TPBord::AbstractGroup LOCA::Bifurcation::HopfBord::ExtendedGroup LOCA::TimeDependent::AbstractGroup • computeMassMatrix() • applyMassMatrix() LOCA::Bifurcation::HopfBord::AbstractGroup • applyComplex() • applyComplexInverse() • computeDCeDp() • computeDCeDxa() Concrete group

17. 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:

18. 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

19. Salinas Example: 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

20. 3D Rayleigh-Benard Problem in 5x5x1 box MPSalsa 208K unknowns 16 Processors

21. Improve robustness Better step size control Improved bifurcation tracking algorithms Debugging Incorporate Multi-vector support Multi-parameter continuation (Henderson, IBM) Constraint enforcement Automatic differentiation Impact more ASCI codes Structural mechanics Where We’re Going From Here