NOX: LOCA Turning Point Continuation Tutorial

Overview

In the LOCA Continuation Tutorial, two turning point bifurcations were discoved in the Chan problem as the parameter $\alpha$ was varied from $0$ to $5$ . Here we provide a brief tutorial on locating one of these turning points and tracking it in the second parameter $\beta$ . There are several methods in LOCA for computing turning points. The example here uses the minimally augmented formulation using the LAPACK direct solver for solving the augmented linear equations. See LOCA::TurningPoint::MinimallyAugmented::ExtendedGroup for a description of this algorithm. Also see LOCA::TurningPoint::MooreSpence::ExtendedGroup for the Moore-Spence turning point formulation. The code fragements discussed below can be found in ChanTPContinuation.C in the Chan subdirectory of the LOCA LAPACK examples directory.

ChanTPContinuation.C

Much of the setup for the turning point continuation problem is the same as for the equilibrium continuation discussed in LOCA Continuation Tutorial. Therefore we will only discuss the differences between the two setups here.

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

int main()
{

  try {
    int n = 100;
    double alpha = 4.0;
    double beta = 0.0;
    double scale = 1.0;
    int maxNewtonIters = 10;

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

By examining the plot in LOCA Continuation Tutorial, a turning point bifurcation occurs near $\alpha=4$ for $\beta = 0$ . We use these values as an initial set of parameter values near the bifurcation. We also set up an output file to store the continuation parameter, solution vector, bifurcation parameter, right null vector, measure of singularity of the Jacobian (sigma) and left null vector at each continuation step. The format is the same as in ChanContinuation.dat, consisting of a series of rows each containing $n+1$ numbers. Three rows are written for each continuation step, the first containing the continuation parameter and the $n$ components of the solution vector, the second containing the bifurcation parameter and the $n$ components of the right null vector, and the third containing sigma and the $n$ components of the left null vector.

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

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

    // Create initial guess for the null vector of jacobian
    Teuchos::RCP<NOX::Abstract::Vector> nullVec =
      Teuchos::rcp(new NOX::LAPACK::Vector(n));
    nullVec->init(1.0);               // initial value 1.0

    // Create initial values for a and b for minimally augmented method
    Teuchos::RCP<NOX::Abstract::Vector> a_vec =
      Teuchos::rcp(new NOX::LAPACK::Vector(n));
    a_vec->init(1.0);

    Teuchos::RCP<NOX::Abstract::Vector> b_vec =
      Teuchos::rcp(new NOX::LAPACK::Vector(n));
    b_vec->init(1.0);

The only additional set up required for turning point tracking in this problem is to compute an initial guess for the null vector of the Jacobian (for the Moore-Spence formulation) or the initial values for the $a$ and $b$ vectors in the minimally augmented formulation. Here we use a vector of all one's in all three cases.

    // Create parameter list
    Teuchos::RCP<Teuchos::ParameterList> paramList =
      Teuchos::rcp(new Teuchos::ParameterList);

    // Create LOCA sublist
    Teuchos::ParameterList& locaParamsList = paramList->sublist("LOCA");

    // Create the stepper sublist and set the stepper parameters
    Teuchos::ParameterList& stepperList = locaParamsList.sublist("Stepper");
    stepperList.set("Continuation Method", "Arc Length");  // Default
    stepperList.set("Continuation Parameter", "beta");     // Must set
    stepperList.set("Initial Value", beta);                // Must set
    stepperList.set("Max Value", 1.0);                     // Must set
    stepperList.set("Min Value", 0.0);                     // Must set
    stepperList.set("Max Steps", 20);                      // Should set
    stepperList.set("Max Nonlinear Iterations", maxNewtonIters); // Should set

    // Create bifurcation sublist
    Teuchos::ParameterList& bifurcationList =
      locaParamsList.sublist("Bifurcation");
    bifurcationList.set("Type", "Turning Point");          // For turning point
    bifurcationList.set("Bifurcation Parameter", "alpha"); // Must set

    // For Moore-Spence formulation w/bordering
    //bifurcationList.set("Formulation", "Moore-Spence");          // Default
    //bifurcationList.set("Solver Method", "Salinger Bordering");  // Default
    //bifurcationList.set("Solver Method", "Phipps Bordering");
    //bifurcationList.set("Bordered Solver Method",
    //                    "LAPACK Direct Solve");   // For Phipps Bordering
    //bifurcationList.set("Length Normalization Vector", nullVec); // Must set
    //bifurcationList.set("Initial Null Vector", nullVec);         // Must set

    // For minimally augmented formulation
    bifurcationList.set("Formulation", "Minimally Augmented");
    bifurcationList.set("Initial A Vector", a_vec);                // Must set
    bifurcationList.set("Initial B Vector", b_vec);                // Must set

    // For minimally augmented method, should set these for good performance
    // Direct solve of bordered equations
    bifurcationList.set("Bordered Solver Method",  "LAPACK Direct Solve");
    // Combine arc-length and turning point bordered rows & columns
    stepperList.set("Bordered Solver Method", "Nested");
    Teuchos::ParameterList& nestedList =
      stepperList.sublist("Nested Bordered Solver");
    // Direct solve of combined bordered system
    nestedList.set("Bordered Solver Method", "LAPACK Direct Solve");

We now set $\beta$ to be the continuation parameter and $\alpha$ to be the bifurcation parameter. We will vary $\beta$ from $0$ to $1$ , computing a value of $\alpha$ for each corresponding value of $\beta$ . The initial value of $\alpha$ is set internally by accessing the component "alpha" in the parameter vector p set below. In the bifurcation sublist, we indicate that we would like to do turning point tracking using the minimally augmented formulation, and pass ref-count pointers to the initial guess for the a and b vectors. Note that these must be casted to NOX::Abstract::Vector ref-count pointers.

Both arc-length continuation and the minimally augmented turning point method add one additional equation to be solved resulting in two nested bordered systems, each adding an additional row and column to the Jacobian matrix. We tell LOCA to combine these rows and columns into one bordered system with two augmented rows and columns (by setting the "Bordered Solver Method" of the "Stepper" sublist to "Nested") and instruct LOCA to use the LAPACK-specific linear solver for solving this 2-bordered system. Note that we must set the "LAPACK Direct Solve" choice twice , once in the "Bifurcation" sublist and once in the "Nested Bordered Solver" sublist. The first specifies the solver for the first and last continuation steps which use natural continuation instead of arc-length continuation. The last specifies the solver for the rest of the continuation steps.

     // Create predictor sublist
    Teuchos::ParameterList& predictorList =
      locaParamsList.sublist("Predictor");
    predictorList.set("Method", "Secant");                         // Default

    // Should use for Moore-Spence w/Salinger Bordering & Secant predictor
    //Teuchos::ParameterList& firstStepPredictor
    //  = predictorList.sublist("First Step Predictor");
    //firstStepPredictor.set("Method", "Random");
    //firstStepPredictor.set("Epsilon", 1.0e-3);

    // Should use for Moore-Spence w/Salinger Bordering & Secant predictor
    //Teuchos::ParameterList& lastStepPredictor
    //  = predictorList.sublist("Last Step Predictor");
    //lastStepPredictor.set("Method", "Random");
    //lastStepPredictor.set("Epsilon", 1.0e-3);

We again use a secant predictor to compute an initial guess at each continuation step. Because we are using the secant predictor, a different predictor must be used for the first step. The default is the "Constant" predictor which is fine for the minimally augmented formulation. However for the Moore-Spence formulation with Salinger bordering, we should use the "Random" predictor because the algorithm is singular. If the problem Jacobian does not depend on the continuation parameter, we can obtain highly ill-conditioned linear solves when not using the random predictor. A random predictor can be chosen for the last step for the same reason.

    // Create step size sublist
    Teuchos::ParameterList& stepSizeList = locaParamsList.sublist("Step Size");
    stepSizeList.set("Method", "Adaptive");                      // Default
    stepSizeList.set("Initial Step Size", 0.1);                  // Should set
    stepSizeList.set("Min Step Size", 1.0e-3);                   // Should set
    stepSizeList.set("Max Step Size", 1.0);                      // Should set

    // Create the "Solver" parameters sublist to be used with NOX Solvers
    Teuchos::ParameterList& nlParams = paramList->sublist("NOX");
    Teuchos::ParameterList& nlPrintParams = nlParams.sublist("Printing");
    nlPrintParams.set("Output Information",
              NOX::Utils::OuterIteration +
              NOX::Utils::OuterIterationStatusTest +
              NOX::Utils::InnerIteration +
              NOX::Utils::Details +
              NOX::Utils::Warning +
              NOX::Utils::StepperIteration +
              NOX::Utils::StepperDetails +
              NOX::Utils::StepperParameters);            // Should set

    // Create LAPACK Factory
    Teuchos::RCP<LOCA::LAPACK::Factory> lapackFactory =
      Teuchos::rcp(new LOCA::LAPACK::Factory);

    // Create global data object
    Teuchos::RCP<LOCA::GlobalData> globalData =
      LOCA::createGlobalData(paramList, lapackFactory);

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

    // Create a group which uses that problem interface. The group will
    // be initialized to contain the default initial guess for the
    // specified problem.
    Teuchos::RCP<LOCA::MultiContinuation::AbstractGroup> grp =
      Teuchos::rcp(new LOCA::LAPACK::Group(globalData, chan));

    grp->setParams(p);

    // Set up the status tests
    Teuchos::RCP<NOX::StatusTest::NormF> statusTestA =
      Teuchos::rcp(new NOX::StatusTest::NormF(1.0e-10,
                          NOX::StatusTest::NormF::Scaled));
    Teuchos::RCP<NOX::StatusTest::MaxIters> statusTestB =
      Teuchos::rcp(new NOX::StatusTest::MaxIters(maxNewtonIters));
    Teuchos::RCP<NOX::StatusTest::Combo> combo =
      Teuchos::rcp(new NOX::StatusTest::Combo(NOX::StatusTest::Combo::OR,
                          statusTestA, statusTestB));

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

    // Solve the nonlinear system
    LOCA::Abstract::Iterator::IteratorStatus status = stepper.run();

    if (status != LOCA::Abstract::Iterator::Finished) {
      if (globalData->locaUtils->isPrintType(NOX::Utils::Error))
    globalData->locaUtils->out()
      << "Stepper failed to converge!" << std::endl;
    }

    // Output the parameter list
    if (globalData->locaUtils->isPrintType(NOX::Utils::StepperParameters)) {
      globalData->locaUtils->out()
    << std::endl << "Final Parameters" << std::endl
    << "****************" << std::endl;
      stepper.getList()->print(globalData->locaUtils->out());
      globalData->locaUtils->out() << std::endl;
    }

    outFile.close();

    LOCA::destroyGlobalData(globalData);
  }

  catch (std::exception& e) {
    cout << e.what() << std::endl;
  }
  catch (const char *s) {
    cout << s << std::endl;
  }
  catch (...) {
    cout << "Caught unknown exception!" << std::endl;
  }

  return 0;
}

The rest of the driver setup is very similar to ChanContinuation.C We slightly tighten the convergence tolerance to demonstrate the well conditioning of the minimally augmented method. In this case we must use the LAPACK factory since we are using the LAPACK-specific solver method "LAPACK Direct Solve".

After running the example and reading the data file ChanTPContinuation.dat, we can plot the continuation parameter $\beta$ versus the bifurcation parameter $\alpha$ to get a locus of turning point bifurcations in the $(\beta,\alpha)$ parameter space:

There are two branches of the bifurcation curve which come together to form a cusp. Starting at $\beta=0$ on one branch, traversing the cusp, and moving to $\beta=0$ on the second branch connects the two turning points shown in the LOCA Continuation Tutorial.

As in the continuation tutorial, the MATLAB code used to generate this plot is shown below.

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

% read dimension of discretization
n = fscanf(fid, '%d', 1);

beta = [];     % array of continuation parameter values at each step
alpha = [];    % array of bifurcation parameter values at each step
sigma = [];    % array of singular estimates at each step
x = [];        % array of solution components at each step
y = [];        % array of right null vector components at each step
z = [];        % array of left null vector components at each step

while ~feof(fid)

  % read beta
  beta = [beta fscanf(fid, '%g', 1)];

  % read x
  x = [x fscanf(fid, '%g', n)];

  % read alpha
  alpha = [alpha fscanf(fid, '%g', 1)];

  % read y
  y = [y fscanf(fid, '%g', n)];

  % read sigma
  sigma = [sigma fscanf(fid, '%g', 1)];

  % read z
  z = [z fscanf(fid, '%g', n)];

end

% close output file
fclose(fid);

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

figure;
plot(beta,alpha,'bo-');
xlabel('\beta');
ylabel('\alpha    ','Rotation',0);
title('Locus of Turning Points');