# Presentation: 297195

## Continuation and Bifurcation Methods Using 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
- 9233 – Computational Sciences
**Trilinos User Group Meeting**
- October 16, 2003

**Continuation and Bifurcation Methods ****Using LOCA**

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

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

- We need algorithms, software, and experience to impact  ASCI- and SciDAC-sized applications.

- Stability/bifurcation analysis provides qualitative information about time evolution of nonlinear systems by computing families of steady-state solutions.

- LOCA:  Library of Continuation Algorithms

## History

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

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

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

## LOCA Designed for Easy Linking to Existing Newton-based Applications

- 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

- LOCA targets existing codes that are:
- Steady-State, Nonlinear
- Newton’s Method
- Large-Scale, Parallel

## Bordering Algorithms Meet these Requirements

## Bordering Algorithms Meet these Requirements

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

**Turning Point Bifurcation**

**Full Newton Algorithm**

**Bordering Algorithm**

## Abstraction of Continuation Process

**Given initial guess            , step size**
**Solve nonlinear equations to find 1****st**** 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**
**	****	**

**LOCA Stepper**

**Predictor modules**

**Step size modules**

**NOX + **
**continuation/**
**bifurcation **
**groups**

## NOX Nonlinear Solver (Kolda, Pawlowski, Hooper, Shadid)

**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 ****(Kolda, Pawlowski, Hooper, Shadid)**

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

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

## Arc-length Group applyJacobianInverse()

- LOCA::Continuation::ArclengthGroup::applyJacobianInverse(const NOX::Abstract::Vector& input,                  					    	              NOX::Abstract::Vector& result)  const {
- 	const LOCA::Continuation::ExtendedVector& con_input =
- 		dynamic_cast<const LOCA::Continuation::ExtendedVector&>(input);
- 	LOCA::Continuation::ExtendedVector& con_result =
- 		dynamic_cast<LOCA::Continuation::ExtendedVector&>(result);
- 	const NOX::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; }

**Arc-length Group applyJacobianInverse()**

## Continuation Groups

**Continuation Groups**

- LOCA::Continuation::ExtendedGroup

- _LOCA::Continuation::NaturalGroup_

- _LOCA::Continuation::ArclengthGroup_

- _LOCA::Continuation::AbstractGroup_
-  setParam()
-  getParam()
-  operator = ()
-  computeDfDp()
-  computeEigenvalues()
-  printSolution()

- Mandatory

- Default implementation
- available

- Optional

- NOX::Abstract::Group

- NOX::Abstract::Group

- Concrete group

## Turning Point, Pitchfork Groups

- LOCA::Continuation::AbstractGroup

- _LOCA::Bifurcation::TPBord::ExtendedGroup_

- _LOCA::Bifurcation::TPBord::AbstractGroup_
-  computeDJnDp()
-  computeDJnDxa()
-  applySingularJacobianInverse()

- NOX::Abstract::Group

- NOX::Abstract::Group

- Concrete group

- LOCA::Continuation::AbstractGroup

## Hopf Groups

- LOCA::Continuation::AbstractGroup

- _LOCA::Bifurcation::HopfBord::ExtendedGroup_

- _LOCA::Bifurcation::HopfBord::AbstractGroup_
-  applyComplex()
-  applyComplexInverse()
-  computeDCeDp()
-  computeDCeDxa()

- NOX::Abstract::Group

- NOX::Abstract::Group

- Concrete group

- _LOCA::TimeDependent::AbstractGroup_
-  computeMassMatrix()
-  applyMassMatrix()

- LOCA::Continuation::AbstractGroup

- _LOCA::Bifurcation::TPBord::AbstractGroup_

## Interfacing Application Codes to LOCA v2.0

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

**Interfacing Application Codes to LOCA v2.0**

## LOCA’s Current Capabilities

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

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

## 3D Rayleigh-Benard Problem in 5x5x1 box

**3D Rayleigh-Benard Problem in 5**x**5**x**1 box**

-  MPSalsa
-  208K unknowns
-  16 Processors

## Where We’re Going From Here

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