Solving LargeScale Continuation and Bifurcation Problems with LOCA

1
Solving Large-Scale Continuation and Bifurcation
Problems with LOCA
Eric Phipps Andy Salinger, Roger Pawlowski 9233
Computational Sciences Trilinos User Group
Meeting November 3, 2004
Sandia is a multiprogram laboratory operated by
Sandia Corporation, a Lockheed Martin
Company,for the United States Department of
Energys 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
• Buckling
• 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
Outline

• What LOCA is/does
• Structural mechanics examples in Salinas
• Snap-through buckling of an arch
• Euler buckling of a beam
• LOCA software and configure options
• Overview of LOCAs design and implementation
• Bordering algorithms
• Super groups/vectors
• Summary of LOCAs current capabilities
• New capabilities since last TUG
• Multi-parameter continuation
• Householder arclength continuation
• Modified turning point bordering algorithm

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

Pitchfork
Hopf
7
Snap-through Buckling of a Symmetric Arch

• Unstressed state of beam is flat
• Negative gravity used to bend beam into an arch
• Ends are hinged
• Beam is loaded in center
• 100 Salinas beam elements
• Continuation parameter is center load

8
Snap-through Buckling of a Symmetric ArchChange
in stability at the turning point
9
Snap-through Buckling of a Symmetric
ArchTracking the turning point in a second
parameter

• Pseudo arc-length continuation on turning point
  equations
• Bending moment as continuation parameter
• Solving for load

10
Euler Buckling of a 3D Beam

• 1x1x50 solid aluminum beam
• 4x4x200 Salinas hex8 elements
• Ends are hinged, right end constrained to move
  along x-axis
• Continuation parameter is horizontal load at
  right end
• Buckles at load 3770, beam theory predicts 3290.

11
LOCA v2.0

• Complete rewrite of LOCA v1.0 (C-library) around
  NOX in C
• Trilinos package inside NOX subdirectory
• Trilinos/packages/nox/src-loca
• Completely dependent on NOX (i.e., you cant
  build LOCA without NOX)
• Leverages NOX interface to application code
• Uses design of NOX to implement continuation and
  bifurcation tracking in a generic way
• All LOCA SQA tools are combined with NOX
• Autoconf/automake
• Documentation
• Bugzilla, Bonsai
• Mail lists

12
Summary of Relevant Configuration Options

• Top-level option to build LOCA in Trilinos
• --enable-loca (Default is on)
• Most configuration options mirror NOX
• --enable-loca-lapack Enable LOCA LAPACK support
  (automatically enabled if NOX LAPACK support is
  enabled)
• --enable-loca-epetra Enable LOCA Epetra support
  (automatically enabled if NOX Epetra support is
  enabled)
• --enable-loca-lapack-examples Build LOCA LAPACK
  examples (automatically enabled if NOX LAPACK
  examples are enabled)
• --enable-loca-epetra-examples Build LOCA Epetra
  examples (automatically enabled if NOX Epetra
  examples are enabled)
• Other options
• --with-loca-anasazi Build LOCA-Anasazi
  interface (for automated eigen-analysis during
  continuation run)
• --with-loca-mf Build LOCA interface to MF
  (multi-parameter continuation)
• Other Trilinos options that must be in place
• --enable-teuchos, --enable-teuchos-complex,
  --enable-anasazi (if LOCA Anasazi support is
  enabled)

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

• LOCA targets existing codes that are
• Steady-State, Nonlinear
• Newtons 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

14
Bordering Algorithms Meet these Requirements
Full Newton Algorithm
15
Bordering Algorithms Meet these Requirements
Full Newton Algorithm
Turning Point Bifurcation
Bordering Algorithm

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

16
Abstraction of Continuation Process
LOCA Stepper

• 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

Predictor modules
NOX continuation/ bifurcation groups
Step size modules
17
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.
18
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
19
Continuation Groups
NOXAbstractGroup
LOCAContinuationExtendedGroup
LOCAContinuationNaturalGroup
LOCAContinuationArclengthGroup
NOXAbstractGroup

• LOCAContinuationAbstractGroup
• setParam()
• getParam()
• operator ()
• computeDfDp()
• computeEigenvalues()
• printSolution()

Mandatory
Default implementation available
Optional
Concrete group
20
Arc-length Group applyJacobianInverse()
LOCAContinuationArclengthGroupapplyJacobianI
nverse(const NOXAbstractVector input,

NOXAbstractVector result) const const
LOCAContinuationExtendedVector con_input
dynamic_castltconst LOCAContinuationExtended
Vectorgt(input) LOCAContinuationExtendedVect
or con_result dynamic_castltLOCAContinuation
ExtendedVectorgt(result) const
NOXAbstractVector input_x
con_input.getXVec() double input_p
con_input.getParam() NOXAbstractVector
result_x con_result.getXVec() double
result_p con_result.getParam() NOXAbstract
Vector b input_x.clone(NOXShapeCopy) unde
rlyingGroupPtr-gtapplyJacobianInverse(input_x,
result_x) underlyingGroupPtr-gtapplyJacobianInver
se(dfdpVecPtr, b) result_p
(predictorVecPtr-gtgetXVec().dot(result_x)
input_p) / (predictorVecPtr-gtg
etXVec().dot(b) predictorVecPtr-gtgetParam())
result_x.update(-result_p, b, 1.0) delete
b
21
Turning Point, Pitchfork Groups
NOXAbstractGroup
NOXAbstractGroup
LOCAContinuationAbstractGroup
LOCAContinuationAbstractGroup
LOCABifurcationTPBordExtendedGroup

• LOCABifurcationTPBordAbstractGroup
• computeDJnDp()
• computeDJnDxa()
• applySingularJacobianInverse()

Concrete group
22
Interfacing Application Codes to LOCA

• Can overload many additional methods if better
  techniques are available
• block solves
• singular matrix solves
• estimating derivatives

23
LOCAs Current Capabilities(New since last TUG
in red)

• Single parameter continuation
• Natural
• Pseudo Arc-length
• Householder arc-length
• Multi-parameter continuation
• Bifurcations
• Turning point
• Modified turning point
• Pitchfork
• Hopf
• Predictors
• Constant (i.e., Euler)
• Tangent
• Secant
• Random

• Restart
• Step size control
• Constant
• Adaptive
• Status tests for bifurcations
• Natural artificial homotopy
• Computing eigenvalues with Anasazi
• Jacobian inverse
• Shift-Invert
• Cayley
• Native support for
• LAPACK
• Epetra

24
Multi-Parameter Continuation

• Multi-parameter continuation supplied through
  Multifario (MF) code (Mike Henderson, IBM)
• General purpose code for covering an implicitly
  defined manifold
• Generic link through LOCA
• LOCA stepper wraps MF driver
• LOCAs continuation groups implement MFs
  interface
• Uses new NOX multi-vector support
• MF library in Trilinos3PL
• Two examples in LOCA repository
• Chan (LAPACK),
• Tcubed (Epetra)
• Resulting data files best visualized with OpenDX
  (www.opendx.org)

25
Multi-Parameter Continuation ExampleChan Problem
26
Householder Pseudo Arc-Length Continuation (New)
Newton solve for pseudo arc-length continuation

• Idea of Homer Walker Solve
• Q is given by a Householder transformation

• Disadvantage Non-generic
• Currently only have Epetra implementation
• Belos implementation coming soon

• Advantage Nearly twice as fast
• Eliminates 2nd solve of bordering method
• Applying Q only requires dot product saxpy
• No change in preconditioning

H.F Walker, SIAM J. Sci. Comput., 1999
27
Improving the Turning Point Bordering Algorithm
blow up in the direction of
as However,
dont. Idea restrict to be
orthogonal to and adjust bordering algorithm
appropriately, e.g., solve
where
28
Modified Turning Point Bordering Algorithm
Solve where
Then
29
3D Rayleigh-Benard Problem in 5x5x1 box (208K
unknowns, 16 processors)

• Salsa application code
• Aztec GMRES solver
• Ifpack RILU preconditioner
• RILU fill factor 2
• RILU overlap 2
• Krylov space 500
• F Turning point
• residual

30
Where Were Going From Here

• Improve robustness/user interface
• More step size control algorithms for homotopy
  problems
• Continued work on improved bifurcation tracking
  algorithms
• Tests for automatic bifurcation location
• Automatic branch switching???
• Incorporate
• High order predictors
• Constraint enforcement
• Complete transition to a multivector-based
  implementation
• Finish off the Belos and Epetra group
  implementations
• Implement LOCA-TSF adaptor
• Tests