LOCA Continuation Tutorial


Here we provide a brief tutorial for using LOCA to perform a simple continuation using the Chan problem (see ChanProblemInterface). The code fragements discussed below can be found in ChanContinuation.C in the Chan subdirectory of the LOCA LAPACK examples directory.

The ChanProblemInterface implements a 1-D finite-difference discretization of the Chan problem:

\[ \frac{d^2 T}{d x^2} + \alpha s \left(1 + \frac{x + 0.5 x^2}{1 + 0.01 x^2}\right) = 0 \]

subject to the boundary conditions $T(0) = T(1) = \beta$. The parameters are $\alpha$, $\beta$, $s$, and $n$, the size of the discretization. The scaling factor $s$ is used to test continuation at different scales of $\alpha$. Below we will track the solution $T(x,\alpha,\beta)$ as $\alpha$ is varied from $0$ to $5$.

ChanContinuation.C Line by Line

#include "LOCA.H"
#include "LOCA_LAPACK.H"
#include "ChanProblemInterface.H"

To use LOCA in your code, you must always include the header LOCA.H. Since this is a LAPACK problem, we also include LOCA_LAPACK.H. Finally, we include the header for the Chan problem, ChanProblemInterface.H.

int main()
  int n = 100;
  double alpha = 0.0;
  double beta = 0.0;
  double scale = 1.0;
  int maxNewtonIters = 20;

  alpha = alpha / scale;

  try {

Next we set up the basic problem parameters. maxNewtonIters is the maximum number of nonlinear solver iterations we wish to take at each continuation step. NOX and LOCA do throw exceptions (but only when serious errors occur), so all NOX and LOCA calculations should be placed in a try block.

    // Create output file to save solutions
    ofstream outFile("ChanContinuation.dat");
    outFile.setf(ios::scientific, ios::floatfield);

    // Save size of discretizations
    outFile << n << endl;

Here we set up a file for saving the solutions computed at each continuation step. The printSolution method of the ChanProblemInterface is set up in such a way that if a file is provided, the current continuation parameter and solution vector are appended to the file, in addition to printing a portion of the solution vector to the screen. The format of this file is a series of rows, with each row containing $n+1$ numbers, the first is the continuation parameter with the remaining $n$ consisting of each component of the solution vector.

    // Set up the problem interface
    ChanProblemInterface chan(n, alpha, beta, scale, outFile);
    LOCA::ParameterVector p;

Next we instantiate the ChanProblemInterface and create a parameter vector to store the values of the problem parameters. Note that it is not necessary to put every problem parameter into the parameter vector, only those that serve as possible continuation or bifurcation parameters need to be supplied.

    // Create a group which uses that problem interface. The group will
    // be initialized to contain the default initial guess for the
    // specified problem.
    LOCA::LAPACK::Group grp(chan);

Next we instantiate the LAPACK group with the Chan problem and then set the parameter vector in the group. From this point on, the LOCA::Stepper, via the LOCA::LAPACK::Group, will take care of setting parameters in the problem interface.

   // Create parameter list
    NOX::Parameter::List paramList;

    // Create LOCA sublist
    NOX::Parameter::List& locaParamsList = paramList.sublist("LOCA");

    // Create the stepper sublist and set the stepper parameters
    NOX::Parameter::List& stepperList = locaParamsList.sublist("Stepper");
    //stepperList.setParameter("Continuation Method", "Natural");
    stepperList.setParameter("Continuation Method", "Arc Length");
    stepperList.setParameter("Continuation Parameter", "alpha");
    stepperList.setParameter("Initial Value", alpha);
    stepperList.setParameter("Max Value", 5.0/scale);
    stepperList.setParameter("Min Value", 0.0/scale);
    stepperList.setParameter("Max Steps", 50);
    stepperList.setParameter("Max Nonlinear Iterations", maxNewtonIters);
    stepperList.setParameter("Enable Arc Length Scaling", true);
    stepperList.setParameter("Goal Arc Length Parameter Contribution", 0.5);
    stepperList.setParameter("Max Arc Length Parameter Contribution", 0.7);
    stepperList.setParameter("Initial Scale Factor", 1.0);
    stepperList.setParameter("Min Scale Factor", 1.0e-8);
    stepperList.setParameter("Enable Tangent Factor Step Size Scaling",true);
    stepperList.setParameter("Min Tangent Factor", -1.0);
    stepperList.setParameter("Tangent Factor Exponent",1.0);
    stepperList.setParameter("Compute Eigenvalues",false);

    // Create bifurcation sublist
    NOX::Parameter::List& bifurcationList = 
    bifurcationList.setParameter("Method", "None");

    // Create predictor sublist
    NOX::Parameter::List& predictorList = locaParamsList.sublist("Predictor");
    //predictorList.setParameter("Method", "Constant");
    predictorList.setParameter("Method", "Tangent");
    //predictorList.setParameter("Method", "Secant");

    // Create step size sublist
    NOX::Parameter::List& stepSizeList = locaParamsList.sublist("Step Size");
    //stepSizeList.setParameter("Method", "Constant");
    stepSizeList.setParameter("Method", "Adaptive");
    stepSizeList.setParameter("Initial Step Size", 0.1/scale);
    stepSizeList.setParameter("Min Step Size", 1.0e-3/scale);
    stepSizeList.setParameter("Max Step Size", 10.0/scale);
    stepSizeList.setParameter("Aggressiveness", 0.5);

    // Set the LOCA Utilities
    NOX::Parameter::List& locaUtilsList = locaParamsList.sublist("Utilities");
    locaUtilsList.setParameter("Output Information", 
                   LOCA::Utils::Warning +
                   LOCA::Utils::StepperIteration +
                   LOCA::Utils::StepperDetails +
                   LOCA::Utils::Solver +
                   LOCA::Utils::Parameters +

Next we set up the LOCA parameters. We are setting up the problem to perform arc-length continuation in the parameter "alpha" from 0 to 5 with a maximum of 50 continuation steps and maxNewtonIters nonlinear iterations per step. The "Enable Arc Length Scaling", "Goal/Max Arc Length Parameter Contribution" and "Scale Factor" parameters deal with arc-length parameter scaling to equilibrate relative contributions between solution and parameter components in the arc-length equations. The "Tangent Factor" parameters deal with reducing the step size when the solution changes by a large amount from step to step (see LOCA Parameter Reference Page for more details). Since we are doing an equilibrium continuation, we set the bifurcation method to "None". We use a tangent predictor and adaptive step size control with an initial step size of 0.1, maximum of 10.0 and minimum of 0.001. The "Aggressiveness" parameter determines how quickly the step size is increased every successful continuation step. Finally, we specify how much information LOCA should print out during the continuation run.

    // Create the "Solver" parameters sublist to be used with NOX Solvers
    NOX::Parameter::List& nlParams = paramList.sublist("NOX");
    nlParams.setParameter("Nonlinear Solver", "Line Search Based");

    NOX::Parameter::List& nlPrintParams = nlParams.sublist("Printing");
    nlPrintParams.setParameter("Output Information", 
              NOX::Utils::Details +
              NOX::Utils::OuterIteration + 
              NOX::Utils::InnerIteration + 

Next we set up the NOX parameters. We use a simple full-step Newton method for the nonlinear solve at each continuation step.

    // Set up the status tests
    NOX::StatusTest::NormF normF(1.0e-8);
    NOX::StatusTest::MaxIters maxIters(maxNewtonIters);
    NOX::StatusTest::Combo comboOR(NOX::StatusTest::Combo::OR, 

Next we create appropriate status tests for the problem. For convergence at each continuation step, we require the extended (solution and parameter components) residual norm be smaller than 1.0e-8 and the number of nonlinear iterations be smaller than maxNewtonIters.

    // Create the stepper  
    LOCA::Stepper stepper(grp, comboOR, paramList);

    // Perform continuation run
    LOCA::Abstract::Iterator::IteratorStatus status = stepper.run();

    if (status != LOCA::Abstract::Iterator::Finished)
      cout << "Stepper failed to converge!" << endl;

Finally we instantiate the stepper, run the continuation, and check the returned status.

    // Get the final solution from the stepper
    const LOCA::LAPACK::Group& finalGroup = 
      dynamic_cast<const LOCA::LAPACK::Group&>(stepper.getSolutionGroup());
    const NOX::LAPACK::Vector& finalSolution = 
      dynamic_cast<const NOX::LAPACK::Vector&>(finalGroup.getX());

    // Output the parameter list
    if (LOCA::Utils::doPrint(LOCA::Utils::Parameters)) {
      cout << endl << "Final Parameters" << endl
       << "****************" << endl;
      cout << endl;


  catch (string& s) {
    cout << s << endl;
  catch (char *s) {
    cout << s << endl;
  catch (...) {
    cout << "Caught unknown exception!" << endl;

  return 0;

Lastly we copy the final solution out of the stepper, print out the final parameter list, close the output file, and catch any thrown exceptions.

After running the example and plotting the maximum of the temperature versus the continuation parameter $\alpha$ at each step, we obtain the following continuation curve with two turning points:


Turning point bifurcations occur near $\alpha=4$ and $\alpha=3$. For a tutorial on locating these turning point bifurcations and tracking them in the second parameter $\beta$, see LOCA Turning Point Continuation Tutorial.

The above plot was generated via MATLAB using the output file ChanContinuation.dat specified above. For those interested, the MATLAB commands used to generate this plot are shown below.

% open output file
fid = fopen('ChanContinuation.dat');

% read dimension of discretization
n = fscanf(fid, '%d', 1);
alpha = []; % array of continuation parameter values at each step
x = [];     % array of solution components at each step
% read values from output file
while ~feof(fid)
  % read alpha
  alpha = [alpha fscanf(fid, '%g', 1)];
  % read x
  x = [x fscanf(fid, '%g', n)];

% close output file

% compute maximum of each temperature profile
maxT = max(x);

ylabel('T_{max}   ','Rotation',0);
title('Arc-length Continuation:  \beta = 0');

Generated on Thu Sep 18 12:43:29 2008 for NOX by doxygen