330 likes | 559 Views
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.
E N D
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 Energy’s National Nuclear Security Administration under contract DE-AC04-94AL85000.
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.
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 LOCA’s design and implementation Bordering algorithms Super groups/vectors Summary of LOCA’s current capabilities New capabilities since last TUG Multi-parameter continuation Householder arclength continuation Modified turning point bordering algorithm Outline
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,
Pseudo Arc-length Continuation Solves for Solution and Parameter Simultaneously
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
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
Snap-through Buckling of a Symmetric ArchChange in stability at the turning point
Locus of turning points 1 solution 3 solutions 1 solution 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
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.
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 can’t 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
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)
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
Pseudo Arc-length Continuation Bordering Algorithm Bordering Algorithms Meet these Requirements Full Newton Algorithm
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
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
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.
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
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
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; }
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
Interfacing Application Codes to LOCA • Can overload many additional methods if better techniques are available • block solves • singular matrix solves • estimating derivatives:
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 LOCA’s Current Capabilities(New since last TUG in red)
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 • LOCA’s continuation groups implement MF’sinterface • 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)
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
Bordering Algorithm for Newton Updates Turning Point Equations Improving the Turning Point Bordering Algorithm blow up in the direction of as However, don’t. Idea: restrict to be orthogonal to and adjust bordering algorithm appropriately, e.g., solve where
Modified Turning Point Bordering Algorithm Solve where Then
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
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 Where We’re Going From Here