LOCA Turning Point Continuation Tutorial

Overview

In 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$. This example uses the turning point bordering method to formulate and solve the extended set of nonlinear equations specifying a turning point bifurcation. Please see the documentation for LOCA::Bifurcation::TPBord::ExtendedGroup for a description of the turning point formulation and how it is solved. 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 = 20;

    // 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, and 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. Two rows are written for each continuation step, the first containing the continuation parameter and the $n$ components of the solution vector, and the second containing the bifurcation parameter and the $n$ components of the null vector.

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

    // Set up the problem interface
    ChanProblemInterface chan(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.
     LOCA::LAPACK::Group grp(chan);
    grp.setParams(p);

    // Create initial guess for the null vector of jacobian
    NOX::LAPACK::Vector nullVec(n);  // length n
    nullVec.init(1.0);             // initial value 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. Here we use a vector of all one's. We will also use this vector as the length scaling vector $\phi$ in the turning point tracking algorithm.

    // 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", "beta");
    stepperList.setParameter("Initial Value", beta);
    stepperList.setParameter("Max Value", 1.0);
    stepperList.setParameter("Min Value", 0.0);
    stepperList.setParameter("Max Steps", 20);
    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",false);
    stepperList.setParameter("Min Tangent Factor", -1.0);
    stepperList.setParameter("Tangent Factor Exponent",1.0);

    // Create bifurcation sublist
    NOX::Parameter::List& bifurcationList = 
      locaParamsList.sublist("Bifurcation");
    bifurcationList.setParameter("Method", "Turning Point");
    bifurcationList.setParameter("Bifurcation Parameter", "alpha");
    bifurcationList.setParameter("Length Normalization Vector", 
             dynamic_cast<NOX::Abstract::Vector*>(&nullVec));
    bifurcationList.setParameter("Initial Null Vector",
             dynamic_cast<NOX::Abstract::Vector*>(&nullVec));

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 above. In the bifurcation sublist, we indicate that we would like to do turning point tracking using the bordering method (LOCA::Bifurcation::TPBord::ExtendedGroup), and pass pointers to the initial guess for the null vector and length scaling vector. Note that these must be casted to NOX::Abstract::Vector pointers.

    // Create predictor sublist
    NOX::Parameter::List& predictorList = locaParamsList.sublist("Predictor");
    //predictorList.setParameter("Method", "Constant");
    predictorList.setParameter("Method", "Secant");
    //predictorList.setParameter("Method", "Random");
    //predictorList.setParameter("Epsilon", 1.0e-3);

    NOX::Parameter::List& firstStepPredictor 
      = predictorList.sublist("First Step Predictor");
    firstStepPredictor.setParameter("Method", "Random");
    firstStepPredictor.setParameter("Epsilon", 1.0e-3);

    NOX::Parameter::List& lastStepPredictor 
      = predictorList.sublist("Last Step Predictor");
    lastStepPredictor.setParameter("Method", "Random");
    lastStepPredictor.setParameter("Epsilon", 1.0e-3);

We now use a secant predictor to compute an initial guess at each continuation step. In general, this is a better choice for bifurcation tracking than the tangent predictor because of the bordering techniques used. The tangent predictor requires a solve of the underlying Jacobian matrix which is singular at the turning point. Because we are using the secant predictor, we must supply a predictor to use for the first step. Here we choose the random predictor. This ensures that the system is perturbed off of the singularity before computing the first step. We choose a random predictor for the last step for the same reason.

    // Create step size sublist
    NOX::Parameter::List& stepSizeList = locaParamsList.sublist("Step Size");
    stepSizeList.setParameter("Method", "Adaptive");
    stepSizeList.setParameter("Initial Step Size", 0.1);
    stepSizeList.setParameter("Min Step Size", 1.0e-3);
    stepSizeList.setParameter("Max Step Size", 1.0);
    stepSizeList.setParameter("Aggressiveness", 0.5);
    stepSizeList.setParameter("Failed Step Reduction Factor", 0.5);
    stepSizeList.setParameter("Successful Step Increase Factor", 1.26); // for constant

    // 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::SolverDetails);

    // 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::OuterIteration + 
                   //NOX::Utils::OuterIterationStatusTest + 
                   //NOX::Utils::InnerIteration +
                   //NOX::Utils::Parameters +
                   //NOX::Utils::Details + 
                   NOX::Utils::Warning);

    // Create the "Line Search" sublist for the "Line Search Based" solver
    NOX::Parameter::List& searchParams = nlParams.sublist("Line Search");
    searchParams.setParameter("Method", "Full Step");

    // Create the newton and  linear solver parameters sublist
  NOX::Parameter::List& directionParameters = nlParams.sublist("Direction");
  NOX::Parameter::List& newtonParameters = directionParameters.sublist("Newton");
  NOX::Parameter::List& linearSolverParameters = newtonParameters.sublist("Linear Solver");

    // Set up the status tests
    NOX::StatusTest::NormF statusTestA(1.0e-5, NOX::StatusTest::NormF::Scaled);
    NOX::StatusTest::MaxIters statusTestB(maxNewtonIters);
    NOX::StatusTest::Combo combo(NOX::StatusTest::Combo::OR, statusTestA, statusTestB);

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

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

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

    // Output the parameter list
    if (NOX::Utils::doPrint(NOX::Utils::Parameters)) {
      cout << endl << "Final Parameters" << endl
       << "****************" << endl;
      stepper.getParameterList().print(cout);
      cout << endl;
    }

    outFile.close();
  }

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

  return 0;
}

The rest of the driver setup is very similar to ChanContinuation.C The only difference being we use a slightly looser nonlinear solver tolerance and suppress some of NOX's screen output.

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:

chan_tp.png

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 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
x = [];        % array of solution components at each step
y = [];        % array of 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)];
  
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');

Generated on Thu Sep 18 12:40:41 2008 for NOX by doxygen 1.3.9.1