Eric Phipps Andy Salinger, Roger Pawlowski Sandia National Laboratories

1. What’s New in NOX…Plus New Techniques for Solving Large-Scale Steady-State and Transient Stability Analysis Problems with LOCA Eric Phipps Andy Salinger, Roger Pawlowski Sandia National Laboratories 3rd Annual Trilinos User’s Group Meeting November 2, 2005 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. Quick overview of the NOX and LOCA Software changes in NOX and LOCA Memory management Epetra support Merit functions Strategies and factories Multi-vectors Miscellaneous New algorithms and techniques in LOCA Constraint enforcement Solving bordered systems via QR Modified turning point bordering Introduction to the “4D” project: Turning transient stability analysis problems in to “steady-state” problems Outline

3. NOX Nonlinear Solver 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.

4. 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 1 3 1 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 • Flutter • El Niño Pitchfork Hopf

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

8. NOX Major refactor started just after Trilinos 6.0 release. The aim was to clean-up memory management, scaling issues and to address a variety of user requests. LOCA Refactor started after Trilinos 5.0 release, but in a parallel source code tree. Refactor will be complete for Trilinos 7.0 release and old code will be removed. Aim was to increase flexibility of software design and reduce code maintenance burden. With the release of Trilinos 7.0, all codes will have to make many changes when calling NOX/LOCA. This was timed to give us a large window between releases to update all codes that rely on NOX/LOCA. NOX and LOCA Software Changes

9. All NOX/LOCA objects now use Teuchos::RefCountPtr whenever a persisting relationship between objects is established. This includes StatusTest, Solver, Stepper, concrete groups and vectors, and parameter lists. All clone() methods return Teuchos::RefCountPtr. User supplied objects are now passed into NOX/LOCA as Teuchos::RefCountPtr objects in the parameter list. Teuchos is now a required dependency, it is no longer optional. This requires NOX to be built as part of Trilinos. Independent builds are still supported, but discouraged.(The new test harness lives at the top Trilinos level and users need these scripts to run the tests.) Memory Management

10. Weighted norm support has been improved. The dot() method of the NOX::Abstract::Vector has been replaced by innerProduct() A vector space object may be provided to NOX::Epetra::Vector to override the default inner product. The old Epetra support is gone and EpetraNew is now Epetra. The NOX/LOCA Epetra interface no longer assumes only one Jacobian matrix exists in the application code. Epetra interface computeJacobian() now accepts a Jacobian argument. Epetra Support

11. NOX now has a pure virtual class that defines a merit function. Can be used in the line searches, trust region, and steepest descent direction calculations. Users can supply a merit function through the parameter list. Merit function Merit Functions

12. Previously, any application/linear algebra specific operation LOCA needed was added to the group Creates a “blob” anti-pattern Difficult for a group to support multiple possible implementations Have moved instead to a strategy interface and factory pattern Users can provide their own strategies and a factory to instantiate them Instantiated strategy is controlled through parameter lists Users can even supply constructor arguments through parameter lists (using Teuchos::RefCountPtr) Makes it easy to provide multiple implementations for a single linear algebra interface (e.g., Epetra) Drastically increases flexibility of LOCA However, requires users to implement more classes Strategy Interfaces & Factories

13. All LOCA algorithm implementations have been rewritten to use multi-vectors instead of single vectors: Same implementation for single parameter continuation versus multi-parameter continuation Block elimination (i.e. bordering) algorithms always involve linear solves with multiple right-hand sides Default multi-vector implementation is provided as an array of vectors LOCA must interact with NOX, which doesn’t use multi-vectors What is the best relationship of a vector to a multi-vector? Multi-Vectors

14. Output All output can be redirected to an output stream and error stream of the user’s choice Note: For Epetra support, we can not redirect the Aztec output to a user ostream Parallel issues have been streamlined and hidden by using the Teuchos::oblackholestream The following will be moved from prerelease to release status in Trilinos 7.0 Broyden Direction Tensor Solver/Direction/Curvilinear Line Search Nonlinear-CG Direction/Line Search Inexact Trust Region Solver and related utilities Multi-physics drivers Many tests have been added and/or expanded Miscellaneous NOX Changes

15. Available in Trilinos 6.0: Constraint enforcement Solving bordered systems via QR Nonsingular turning point bordering methods Will be available in Trilinos 7.0 Nonsingular pitchfork and Hopf bordering methods Minimally augment turning point formulations New Algorithms/Techniques in LOCA

16. Users can provide additional algebraic constraints to be solved along with application code’s residual (e.g., boundary conditions, etc…) Users implement a simple interface to represent constraints Continuation is a special case of this Constraint Enforcement

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

18. Solving Bordered Systems via QR Extension of Householder pseudo-arclength technique by Homer Walker1 Rearranged Bordered System Compact WY Representation2 QR Factorization • Advantages • 1 solve per nonlinear iteration • P is nonsingular even when J is singular • Applying P only requires dot products + saxpy’s • Preconditioner for J works well as preconditioner for P • Does not change the dimension of linear system • Disadvantages • Requires solver specific implementation, but have general Epetra implementation where 1H.F Walker, SIAM J. Sci. Comput., 1999 2R. Schreiber, SIAM J. Stat. Comput., 1989

19. Turning Point Bifurcation Full Newton Algorithm Bordering Algorithm Turning Point Identification … but 4 solves of J per Newton iteration are used to drive J singular!

20. Modified Turning Point Bordering Algorithm Solve 5 bordered systems of equations using QR approach Then

21. Snap-through Buckling of a Shallow Cap • Salinas structural mechanics app • 66K tri-shell elements, 200K unkowns, 30 procs • FEI/Epetra interface • AztecOO GMRES solver • Ifpack RILU(k) preconditioner

22. Snap-through Buckling of a Shallow CapTurning Point Identification • RILU fill factor: 6 • RILU overlap: 6 • Krylov space: 2000 • F = Turning point • residual

23. Snap-through Buckling of a Symmetric Cap Turning Point Continuation • Pseudo-arclength continuation of turning point equations • Modified bordering method

24. Instead, solve for all solutions at once: where Sophisticated Steady-State Design/Analysis Tools can Impact Transient Problems with Space-Time “4D’ Formulation Transient Simulation of: First solve: Then solve: Then solve: … and with Newton solve:

25. get Jacobian, K(y) “4D” Super-Interface get Jacobian, Ji(xi,p) get Mass Matrix, Bi Trilinos Implementation of “4D” Capability: creates a “Super” NOX/LOCA interface NOX/LOCA solve Jacobian, J-1f solve Jacobian, K-1f get Jacobian, J(x,p) Trilinos’ Linear Solvers: Aztec/Ifpack/ML/ Amesos/Belos/Anasazi NOX/LOCA Application Interface Application Code

26. Each processor owns 1 time step for the entire spatial domain proc 0: proc 0: proc 0: proc 1: proc 0: proc 1: proc 0: proc 0: proc 2: Each processor owns 4 time steps for ¼ of the spatial domain proc 0: proc 3: proc 1: Each processor owns 2 time steps for ½ of the spatial domain Space and Time can be partitioned independently, Ex: 4 Time Steps on 4 Procs Spatial Partition Space-Time Partition Proc 0: Proc 1: Proc 2: Proc 3:

27. g Space-Time solutions can be constrained so that trajectories meet design criteria Transient: 16K Unknowns solved for every time step, and repeated for every parameter value Then, continuation can find curve in parameter space of parameter pairs that meet design criteria What parameter value leads to a maximum Temperature along the trajectory of 0.7? Rector explosion model in MPSalsa Space-Time: (16K x 32time steps + 2 constraints) unknowns solved once for trajectory which meets design constraint

28. Future Capabilities and Algorithms for Space-Time Solutions of large problems • Current and Future Work: • Preconditioning for efficient solves of block matrix (Dunlavy) • Periodic orbit tracking

29. Points of Contact • NOX: Trilinos nonlinear solver package • Roger Pawlowski (rppawlo@sandia.gov) • www.software.sandia.gov/nox • LOCA: Trilinos continuation and bifurcation package • Sub-package of NOX • Eric Phipps (etphipp@sandia.gov) • Andy Salinger (agsalin@sandia.gov) • www.software.sandia.gov/nox • 4D: Space-time formulation of transient stability analysis/design problems • Andy Salinger (agsalin@sandia.gov) • Dan Dunlavy (dmdunla@sandia.gov)