NOX Development
Loading...
Searching...
No Matches
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);
Implementation of NOX::Abstract::Vector for STL std::vector<double> (using LAPACK for some computatio...
Definition NOX_LAPACK_Vector.H:78

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",
// 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);
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,
Teuchos::RCP<NOX::StatusTest::MaxIters> statusTestB =
Teuchos::rcp(new NOX::StatusTest::MaxIters(maxNewtonIters));
Teuchos::RCP<NOX::StatusTest::Combo> combo =
statusTestA, statusTestB));
// Create the stepper
LOCA::Stepper stepper(globalData, grp, combo, paramList);
// Solve the nonlinear system
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;
}
IteratorStatus
Enumerated type for status of the iterator.
Definition LOCA_Abstract_Iterator.H:101
@ Finished
The iterator is finished.
Definition LOCA_Abstract_Iterator.H:103
Implementation of the LOCA::Abstract::Factory for LAPACK groups.
Definition LOCA_LAPACK_Factory.H:68
Extension of the NOX::LAPACK::Group to LOCA.
Definition LOCA_LAPACK_Group.H:106
LOCA's container for holding a set of parameters that are used by the LOCA continuation routines.
Definition LOCA_Parameter_Vector.H:73
int addParameter(std::string label, double value=0.0)
Adds a parameter to the list. Returns the index value assigned to the parameter.
Definition LOCA_Parameter_Vector.C:79
Implementation of LOCA::Abstract::Iterator for computing points along a continuation curve.
Definition LOCA_Stepper.H:110
Arbitrary combination of status tests.
Definition NOX_StatusTest_Combo.H:84
@ OR
Logically "OR" together the results of the tests in this combination.
Definition NOX_StatusTest_Combo.H:96
Failure test based on the maximum number of nonlinear solver iterations.
Definition NOX_StatusTest_MaxIters.H:76
Various convergence tests based on the norm of the residual.
Definition NOX_StatusTest_NormF.H:119
@ Scaled
Scale the norm by the length of the vector.
Definition NOX_StatusTest_NormF.H:128
@ InnerIteration
2^2
Definition NOX_Utils.H:152
@ OuterIteration
2^1
Definition NOX_Utils.H:150
@ Warning
2^0
Definition NOX_Utils.H:148
@ StepperDetails
2^9 – For LOCA
Definition NOX_Utils.H:166
@ StepperIteration
2^8 – For LOCA
Definition NOX_Utils.H:164
@ Error
Errors are always printed.
Definition NOX_Utils.H:146
@ StepperParameters
2^10 – For LOCA
Definition NOX_Utils.H:168
@ OuterIterationStatusTest
2^5
Definition NOX_Utils.H:158
@ Details
2^4
Definition NOX_Utils.H:156

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