# Presentation: 301557

## What’s New in NOX...Plus New Techniques for Solving Large-Scale Steady-State and Transient Stability Analysis Problems with LOCA

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

- Eric Phipps
- Andy Salinger, Roger Pawlowski
- Sandia National Laboratories
- 3rd Annual Trilinos User’s Group Meeting
- November 2, 2005

**What’s New in NOX...Plus New Techniques for Solving Large-Scale Steady-State and Transient Stability Analysis Problems with LOCA**

## Outline

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

## NOX Nonlinear Solver 

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

**NOX Nonlinear Solver**** **

| Group | Vector |
| --- | --- |
| computeF() | innerProduct() |
| computeJacobian() | scale() |
| computeNewton() | norm() |
| applyJacobianInverse() | update() |

## LOCA:  Library of 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**

**Application code provides:**
**Nonlinear steady-state residual and Jacobian fill:**
**Newton-like linear solves:**

- Second parameter,

- 1

- 1

- 3

## Pseudo Arc-length Continuation Solves for Solution and Parameter Simultaneously

**Pseudo Arc-length Continuation Solves for Solution and Parameter Simultaneously**

## Codimension 1 Bifurcations

**Codimension 1 Bifurcations**

**Turning Point**

**Pitchfork**

**Hopf**

**Combustion**
**Buckling of an Arch**
**Buckling of a Beam**
**Pattern formation**
**Cell differentiation (morphogenesis)**
**Vortex Shedding**
**Predator-Prey models**
**Flutter**
**El Ni****ñ****o**

## Super Vectors and Super Groups

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

**Super Vectors and Super Groups**

## NOX and LOCA Software Changes

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

## Memory Management

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

## Epetra Support

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

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

**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 Functions**

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

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

## 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?

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

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

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

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

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

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

**Users implement a simple interface to represent constraints**
**Continuation is a special case of this**

**Constraint Enforcement**

## Bordered Systems of Equations

- Only requires solves of J but
- Requires m+l linear solves
- Has difficulty when J is nearly singular

- Solving bordered systems of equations is a ubiquitous computation:

## Solving Bordered Systems via QR

- 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

- 1H.F Walker, SIAM J. Sci. Comput., 1999
- 2R. Schreiber, SIAM J. Stat. Comput., 1989

- Extension of Householder pseudo-arclength technique by Homer Walker1

- QR Factorization

- Compact WY Representation2

- where

- Rearranged Bordered System

## Turning Point Identification

**... ****but 4 solves of J per Newton iteration are used to drive J singular!**

## Modified Turning Point Bordering Algorithm

**Then **

**Modified Turning Point Bordering Algorithm**

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

## Snap-through Buckling of a Shallow CapTurning Point Identification

**Snap-through Buckling of a Shallow Cap****Turning Point Identification**

- RILU fill factor:  6
- RILU overlap: 6
- Krylov space:  2000
- F = Turning point
- 	residual

## Snap-through Buckling of a Symmetric Cap Turning Point Continuation

- Pseudo-arclength continuation of turning point equations
- Modified bordering method

## 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:**

## Trilinos Implementation of “4D” Capability: creates a “Super” NOX/LOCA interface 

- solve Jacobian, K-1f

**Trilinos Implementation of “4D” Capability: creates a “Super” NOX/LOCA interface **

- NOX/LOCA

- Trilinos’ Linear Solvers:
- Aztec/Ifpack/ML/
- Amesos/Belos/Anasazi

- NOX/LOCA Application Interface

- Application Code

- get Jacobian, J(x,p)

## Space and Time can be partitioned independently, Ex: 4 Time Steps on 4 Procs

**Space and Time can be partitioned independently, ****Ex: 4 Time Steps on 4 Procs**

- Spatial Partition

- Space-Time Partition

- Proc 0:

- Proc 1:

- Proc 3:

- Proc 2:

##  Space-Time solutions can be constrained so that trajectories meet design criteria

- Space-Time: (16K x 32time steps + 2 constraints) unknowns solved once for trajectory which meets design constraint

- Transient: 16K Unknowns solved for every time step, and repeated for every parameter value

- What parameter value leads to a maximum Temperature along the trajectory of 0.7?

- Then, continuation can find curve in parameter space of parameter pairs that meet design criteria

- Rector explosion model in MPSalsa

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

## Points of Contact

- NOX:  Trilinos nonlinear solver package
- Roger Pawlowski (_[rppawlo@sandia.gov](mailto:rppawlo@sandia.gov)_)
- _[www.software.](http://www.software.sandia.gov/nox)__[sandia](http://www.software.sandia.gov/nox)__[.gov/nox](http://www.software.sandia.gov/nox)_
- LOCA:  Trilinos continuation and bifurcation package
- Sub-package of NOX
- Eric Phipps (_[etphipp@sandia.gov](mailto:etphipp@sandia.gov)_)
- Andy Salinger (_[agsalin@sandia.gov](mailto:agsalin@sandia.gov)_)
- _[www.software.sandia.gov/nox](http://www.software.sandia.gov/nox)_
- 4D:  Space-time formulation of transient stability analysis/design problems
- Andy Salinger (_[agsalin@sandia.gov](mailto:agsalin@sandia.gov)_)
- Dan Dunlavy (_dmdunla@sandia.gov_)