Overview

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);
outFile.precision(14);
 
// Save size of discretizations
outFile << n << std::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.

// 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 Method", "Natural");
stepperList.set("Continuation Parameter", "alpha");  // Must set
stepperList.set("Initial Value", alpha);             // Must set
stepperList.set("Max Value", 5.0/scale);             // Must set
stepperList.set("Min Value", 0.0/scale);             // Must set
stepperList.set("Max Steps", 50);                    // Should set
stepperList.set("Max Nonlinear Iterations", maxNewtonIters); // Should set
stepperList.set("Compute Eigenvalues",false);        // Default
 
// Create bifurcation sublist
Teuchos::ParameterList& bifurcationList =
  locaParamsList.sublist("Bifurcation");
bifurcationList.set("Type", "None");                 // Default
 
// Create predictor sublist
Teuchos::ParameterList& predictorList =
  locaParamsList.sublist("Predictor");
predictorList.set("Method", "Secant");               // Default
//predictorList.set("Method", "Constant");
//predictorList.set("Method", "Tangent");
 
// Create step size sublist
Teuchos::ParameterList& stepSizeList = locaParamsList.sublist("Step Size");
stepSizeList.set("Method", "Adaptive");             // Default
stepSizeList.set("Initial Step Size", 0.1/scale);   // Should set
stepSizeList.set("Min Step Size", 1.0e-3/scale);    // Should set
stepSizeList.set("Max Step Size", 10.0/scale);      // Should set

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. Since we are doing an equilibrium continuation, we set the bifurcation method to "None". We use a secant predictor and adaptive step size control with an initial step size of 0.1, maximum of 10.0 and minimum of 0.001.

// 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::Details +
          NOX::Utils::OuterIteration +
          NOX::Utils::InnerIteration +
          NOX::Utils::Warning +
          NOX::Utils::StepperIteration +
          NOX::Utils::StepperDetails +
          NOX::Utils::StepperParameters);  // Should set

Next we set up the NOX parameters. We use a simple full-step Newton method for the nonlinear solve at each continuation step (the default), and specify what output we want printed to the screen.

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

We now create the LAPACK factory and the global data object. The factory is optional, and allows the use of LAPACK-specific strategies (e.g., the LAPACK eigensolver DGGEV). However the global data object is manditory. If no factory is provided, the second argument to createGlobalData is empty.

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

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.
Teuchos::RCP<LOCA::MultiContinuation::AbstractGroup> grp =
  Teuchos::rcp(new LOCA::LAPACK::Group(globalData, chan));
 
grp->setParams(p);

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.

// Set up the status tests
Teuchos::RCP<NOX::StatusTest::NormF> normF =
  Teuchos::rcp(new NOX::StatusTest::NormF(1.0e-8));
Teuchos::RCP<NOX::StatusTest::MaxIters> maxIters =
  Teuchos::rcp(new NOX::StatusTest::MaxIters(maxNewtonIters));
Teuchos::RCP<NOX::StatusTest::Generic> comboOR =
  Teuchos::rcp(new NOX::StatusTest::Combo(NOX::StatusTest::Combo::OR,
                      normF,
                      maxIters));

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(globalData, grp, comboOR, paramList);
 
// Perform continuation run
LOCA::Abstract::Iterator::IteratorStatus status = stepper.run();
 
// Check for convergence
if (status != LOCA::Abstract::Iterator::Finished) {
  if (globalData->locaUtils->isPrintType(NOX::Utils::Error))
globalData->locaUtils->out()
  << "Stepper failed to converge!" << std::endl;
}

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

    // Get the final solution from the stepper
    Teuchos::RCP<const LOCA::LAPACK::Group> finalGroup =
      Teuchos::rcp_dynamic_cast<const LOCA::LAPACK::Group>(stepper.getSolutionGroup());
    const NOX::LAPACK::Vector& finalSolution =
      dynamic_cast<const NOX::LAPACK::Vector&>(finalGroup->getX());
 
    // Print final solution
    globalData->locaUtils->out()
                << std::endl << "Final solution is " << std::endl;
    finalGroup->printSolution(finalSolution,
                      finalGroup->getParam("alpha"));
 
    // 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;
}

Lastly we copy the final solution out of the stepper, print out the final parameter list, close the output file, destroy the global data object, 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 the 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)];
 
end
 
% close output file
fclose(fid);
 
% compute maximum of each temperature profile
maxT = max(x);
 
plot(alpha,maxT,'bo-');
xlabel('\alpha');
ylabel('T_{max}   ','Rotation',0);
title('Arc-length Continuation:  \beta = 0');