00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042 #include "LOCA_Pitchfork_MooreSpence_ExtendedGroup.H"
00043 #include "LOCA_Pitchfork_MooreSpence_AbstractGroup.H"
00044 #include "LOCA_Pitchfork_MooreSpence_SolverStrategy.H"
00045 #include "LOCA_Parameter_Vector.H"
00046 #include "Teuchos_ParameterList.hpp"
00047 #include "LOCA_GlobalData.H"
00048 #include "LOCA_Factory.H"
00049 #include "LOCA_Parameter_SublistParser.H"
00050 #include "NOX_Utils.H"
00051 #include "LOCA_ErrorCheck.H"
00052
00053 LOCA::Pitchfork::MooreSpence::ExtendedGroup::ExtendedGroup(
00054 const Teuchos::RCP<LOCA::GlobalData>& global_data,
00055 const Teuchos::RCP<LOCA::Parameter::SublistParser>& topParams,
00056 const Teuchos::RCP<Teuchos::ParameterList>& tpParams,
00057 const Teuchos::RCP<LOCA::Pitchfork::MooreSpence::AbstractGroup>& g)
00058 : LOCA::Extended::MultiAbstractGroup(),
00059 LOCA::MultiContinuation::AbstractGroup(),
00060 globalData(global_data),
00061 parsedParams(topParams),
00062 pitchforkParams(tpParams),
00063 grpPtr(g),
00064 xMultiVec(globalData, g->getX(), 1),
00065 fMultiVec(globalData, g->getX(), 2),
00066 newtonMultiVec(globalData, g->getX(), 1),
00067 asymMultiVec(),
00068 lengthMultiVec(),
00069 xVec(),
00070 fVec(),
00071 ffMultiVec(),
00072 dfdpMultiVec(),
00073 newtonVec(),
00074 asymVec(),
00075 lengthVec(),
00076 solverStrategy(),
00077 index_f(1),
00078 index_dfdp(1),
00079 bifParamID(1),
00080 isValidF(false),
00081 isValidJacobian(false),
00082 isValidNewton(false)
00083 {
00084 const char *func = "LOCA::Pitchfork::MooreSpence::ExtendedGroup()";
00085
00086
00087 *(xMultiVec.getColumn(0)->getXVec()) = g->getX();
00088
00089 if (!pitchforkParams->isParameter("Bifurcation Parameter")) {
00090 globalData->locaErrorCheck->throwError(func,
00091 "\"Bifurcation Parameter\" name is not set!");
00092 }
00093 string bifParamName = pitchforkParams->get(
00094 "Bifurcation Parameter",
00095 "None");
00096 const ParameterVector& p = grpPtr->getParams();
00097 bifParamID[0] = p.getIndex(bifParamName);
00098
00099 if (!pitchforkParams->isParameter("Antisymmetric Vector")) {
00100 globalData->locaErrorCheck->throwError(func,
00101 "\"Antisymmetric Vector\" is not set!");
00102 }
00103 Teuchos::RCP<NOX::Abstract::Vector> asymVecPtr =
00104 (*pitchforkParams).INVALID_TEMPLATE_QUALIFIER
00105 get< Teuchos::RCP<NOX::Abstract::Vector> >("Antisymmetric Vector");
00106
00107 if (!pitchforkParams->isParameter("Length Normalization Vector")) {
00108 globalData->locaErrorCheck->throwError(func,
00109 "\"Length Normalization Vector\" is not set!");
00110 }
00111 Teuchos::RCP<NOX::Abstract::Vector> lenVecPtr =
00112 (*pitchforkParams).INVALID_TEMPLATE_QUALIFIER
00113 get< Teuchos::RCP<NOX::Abstract::Vector> >("Length Normalization Vector");
00114
00115 if (!pitchforkParams->isParameter("Initial Null Vector")) {
00116 globalData->locaErrorCheck->throwError(func,
00117 "\"Initial Null Vector\" is not set!");
00118 }
00119 Teuchos::RCP<NOX::Abstract::Vector> nullVecPtr =
00120 (*pitchforkParams).INVALID_TEMPLATE_QUALIFIER
00121 get< Teuchos::RCP<NOX::Abstract::Vector> >("Initial Null Vector");
00122
00123 bool perturbSoln = pitchforkParams->get(
00124 "Perturb Initial Solution",
00125 false);
00126 double perturbSize = pitchforkParams->get(
00127 "Relative Perturbation Size",
00128 1.0e-3);
00129 asymMultiVec =
00130 asymVecPtr->createMultiVector(1, NOX::DeepCopy);
00131 lengthMultiVec =
00132 lenVecPtr->createMultiVector(1, NOX::DeepCopy);
00133 *(xMultiVec.getColumn(0)->getNullVec()) = *nullVecPtr;
00134
00135
00136 solverStrategy =
00137 globalData->locaFactory->createMooreSpencePitchforkSolverStrategy(
00138 parsedParams,
00139 pitchforkParams);
00140
00141
00142 setupViews();
00143
00144 init(perturbSoln, perturbSize);
00145 }
00146
00147 LOCA::Pitchfork::MooreSpence::ExtendedGroup::ExtendedGroup(
00148 const LOCA::Pitchfork::MooreSpence::ExtendedGroup& source,
00149 NOX::CopyType type)
00150 : globalData(source.globalData),
00151 parsedParams(source.parsedParams),
00152 pitchforkParams(source.pitchforkParams),
00153 grpPtr(Teuchos::rcp_dynamic_cast<LOCA::Pitchfork::MooreSpence::AbstractGroup>(source.grpPtr->clone(type))),
00154 xMultiVec(source.xMultiVec, type),
00155 fMultiVec(source.fMultiVec, type),
00156 newtonMultiVec(source.newtonMultiVec, type),
00157 asymMultiVec(source.asymMultiVec->clone(type)),
00158 lengthMultiVec(source.lengthMultiVec->clone(type)),
00159 xVec(),
00160 fVec(),
00161 ffMultiVec(),
00162 dfdpMultiVec(),
00163 newtonVec(),
00164 asymVec(),
00165 lengthVec(),
00166 solverStrategy(source.solverStrategy),
00167 index_f(1),
00168 index_dfdp(1),
00169 bifParamID(source.bifParamID),
00170 isValidF(source.isValidF),
00171 isValidJacobian(source.isValidJacobian),
00172 isValidNewton(source.isValidNewton)
00173 {
00174
00175
00176 solverStrategy =
00177 globalData->locaFactory->createMooreSpencePitchforkSolverStrategy(
00178 parsedParams,
00179 pitchforkParams);
00180
00181
00182 setupViews();
00183
00184 if (type == NOX::ShapeCopy) {
00185 isValidF = false;
00186 isValidJacobian = false;
00187 isValidNewton = false;
00188 }
00189 }
00190
00191 LOCA::Pitchfork::MooreSpence::ExtendedGroup::~ExtendedGroup()
00192 {
00193 }
00194
00195 NOX::Abstract::Group&
00196 LOCA::Pitchfork::MooreSpence::ExtendedGroup::operator=(
00197 const NOX::Abstract::Group& source)
00198 {
00199 copy(source);
00200 return *this;
00201 }
00202
00203 Teuchos::RCP<NOX::Abstract::Group>
00204 LOCA::Pitchfork::MooreSpence::ExtendedGroup::clone(NOX::CopyType type) const
00205 {
00206 return
00207 Teuchos::rcp(new LOCA::Pitchfork::MooreSpence::ExtendedGroup(*this,
00208 type));
00209 }
00210
00211 void
00212 LOCA::Pitchfork::MooreSpence::ExtendedGroup::setX(
00213 const NOX::Abstract::Vector& y)
00214 {
00215 const LOCA::Pitchfork::MooreSpence::ExtendedVector& yy =
00216 dynamic_cast<const LOCA::Pitchfork::MooreSpence::ExtendedVector&>(y);
00217 grpPtr->setX( *yy.getXVec() );
00218 *xVec = y;
00219 setBifParam(xVec->getBifParam());
00220
00221 isValidF = false;
00222 isValidJacobian = false;
00223 isValidNewton = false;
00224 }
00225
00226 void
00227 LOCA::Pitchfork::MooreSpence::ExtendedGroup::computeX(
00228 const NOX::Abstract::Group& g,
00229 const NOX::Abstract::Vector& d,
00230 double step)
00231 {
00232 const LOCA::Pitchfork::MooreSpence::ExtendedGroup& gg =
00233 dynamic_cast<const LOCA::Pitchfork::MooreSpence::ExtendedGroup&>(g);
00234 const LOCA::Pitchfork::MooreSpence::ExtendedVector& dd =
00235 dynamic_cast<const LOCA::Pitchfork::MooreSpence::ExtendedVector&>(d);
00236
00237 grpPtr->computeX(*(gg.grpPtr), *dd.getXVec(), step);
00238 xVec->update(1.0, gg.getX(), step, dd, 0.0);
00239 setBifParam(xVec->getBifParam());
00240
00241 isValidF = false;
00242 isValidJacobian = false;
00243 isValidNewton = false;
00244 }
00245
00246 NOX::Abstract::Group::ReturnType
00247 LOCA::Pitchfork::MooreSpence::ExtendedGroup::computeF()
00248 {
00249 if (isValidF)
00250 return NOX::Abstract::Group::Ok;
00251
00252 string callingFunction =
00253 "LOCA::Pitchfork::MooreSpence::ExtendedGroup::computeF()";
00254 NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
00255 NOX::Abstract::Group::ReturnType status;
00256
00257 double sigma = xVec->getSlack();
00258
00259
00260 if (!grpPtr->isF()) {
00261 status = grpPtr->computeF();
00262 finalStatus =
00263 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00264 finalStatus,
00265 callingFunction);
00266 }
00267 fVec->getXVec()->update(1.0, grpPtr->getF(), sigma, *asymVec, 0.0);
00268
00269
00270 if (!grpPtr->isJacobian()) {
00271 status = grpPtr->computeJacobian();
00272 finalStatus =
00273 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00274 finalStatus,
00275 callingFunction);
00276 }
00277
00278
00279 status = grpPtr->applyJacobian(*(xVec->getNullVec()),
00280 *(fVec->getNullVec()));
00281 finalStatus =
00282 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00283 callingFunction);
00284
00285
00286 fVec->getSlack() = grpPtr->innerProduct(*(xVec->getXVec()), *asymVec);
00287
00288
00289 fVec->getBifParam() = lTransNorm(*(xVec->getNullVec())) - 1.0;
00290
00291 isValidF = true;
00292
00293 return finalStatus;
00294 }
00295
00296 NOX::Abstract::Group::ReturnType
00297 LOCA::Pitchfork::MooreSpence::ExtendedGroup::computeJacobian()
00298 {
00299 if (isValidJacobian)
00300 return NOX::Abstract::Group::Ok;
00301
00302 string callingFunction =
00303 "LOCA::Pitchfork::MooreSpence::ExtendedGroup::computeJacobian()";
00304 NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
00305 NOX::Abstract::Group::ReturnType status;
00306
00307
00308
00309
00310
00311 status = grpPtr->computeDfDpMulti(bifParamID,
00312 *fMultiVec.getXMultiVec(),
00313 false);
00314 finalStatus =
00315 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00316 callingFunction);
00317
00318
00319 double sigma = xVec->getSlack();
00320 fVec->getXVec()->update(sigma, *asymVec, 1.0);
00321
00322
00323 status = grpPtr->computeDJnDpMulti(bifParamID,
00324 *(xVec->getNullVec()),
00325 *fMultiVec.getNullMultiVec(),
00326 isValidF);
00327
00328 finalStatus =
00329 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00330 callingFunction);
00331
00332
00333 status = grpPtr->computeJacobian();
00334 finalStatus =
00335 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00336 callingFunction);
00337
00338 solverStrategy->setBlocks(
00339 grpPtr,
00340 Teuchos::rcp(this, false),
00341 asymMultiVec,
00342 xVec->getNullVec(),
00343 fVec->getNullVec(),
00344 fMultiVec.getColumn(1)->getXVec(),
00345 fMultiVec.getColumn(1)->getNullVec());
00346
00347 isValidJacobian = true;
00348
00349 return finalStatus;
00350 }
00351
00352 NOX::Abstract::Group::ReturnType
00353 LOCA::Pitchfork::MooreSpence::ExtendedGroup::computeGradient()
00354 {
00355 return NOX::Abstract::Group::NotDefined;
00356 }
00357
00358 NOX::Abstract::Group::ReturnType
00359 LOCA::Pitchfork::MooreSpence::ExtendedGroup::computeNewton(
00360 Teuchos::ParameterList& params)
00361 {
00362 if (isValidNewton)
00363 return NOX::Abstract::Group::Ok;
00364
00365 string callingFunction =
00366 "LOCA::Pitchfork::MooreSpence::ExtendedGroup::computeNewton()";
00367 NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
00368 NOX::Abstract::Group::ReturnType status;
00369
00370
00371 if (!isF()) {
00372 status = computeF();
00373 finalStatus =
00374 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00375 finalStatus,
00376 callingFunction);
00377 }
00378
00379
00380 if (!isJacobian()) {
00381 status = computeJacobian();
00382 finalStatus =
00383 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00384 finalStatus,
00385 callingFunction);
00386 }
00387
00388
00389 newtonMultiVec.init(0.0);
00390
00391
00392 status = solverStrategy->solve(params, *ffMultiVec, newtonMultiVec);
00393 finalStatus =
00394 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00395 callingFunction);
00396
00397 newtonMultiVec.scale(-1.0);
00398
00399 isValidNewton = true;
00400
00401 return finalStatus;
00402 }
00403
00404 NOX::Abstract::Group::ReturnType
00405 LOCA::Pitchfork::MooreSpence::ExtendedGroup::applyJacobian(
00406 const NOX::Abstract::Vector& input,
00407 NOX::Abstract::Vector& result) const
00408 {
00409
00410 Teuchos::RCP<NOX::Abstract::MultiVector> mv_input =
00411 input.createMultiVector(1, NOX::DeepCopy);
00412 Teuchos::RCP<NOX::Abstract::MultiVector> mv_result =
00413 result.createMultiVector(1, NOX::DeepCopy);
00414
00415
00416 NOX::Abstract::Group::ReturnType status =
00417 applyJacobianMultiVector(*mv_input, *mv_result);
00418
00419
00420 result = (*mv_result)[0];
00421
00422 return status;
00423 }
00424
00425 NOX::Abstract::Group::ReturnType
00426 LOCA::Pitchfork::MooreSpence::ExtendedGroup::applyJacobianTranspose(
00427 const NOX::Abstract::Vector& input,
00428 NOX::Abstract::Vector& result) const
00429 {
00430
00431 Teuchos::RCP<NOX::Abstract::MultiVector> mv_input =
00432 input.createMultiVector(1, NOX::DeepCopy);
00433 Teuchos::RCP<NOX::Abstract::MultiVector> mv_result =
00434 result.createMultiVector(1, NOX::DeepCopy);
00435
00436
00437 NOX::Abstract::Group::ReturnType status =
00438 applyJacobianTransposeMultiVector(*mv_input, *mv_result);
00439
00440
00441 result = (*mv_result)[0];
00442
00443 return status;
00444 }
00445
00446 NOX::Abstract::Group::ReturnType
00447 LOCA::Pitchfork::MooreSpence::ExtendedGroup::applyJacobianInverse(
00448 Teuchos::ParameterList& params,
00449 const NOX::Abstract::Vector& input,
00450 NOX::Abstract::Vector& result) const
00451 {
00452
00453 Teuchos::RCP<NOX::Abstract::MultiVector> mv_input =
00454 input.createMultiVector(1, NOX::DeepCopy);
00455 Teuchos::RCP<NOX::Abstract::MultiVector> mv_result =
00456 result.createMultiVector(1, NOX::DeepCopy);
00457
00458
00459 NOX::Abstract::Group::ReturnType status =
00460 applyJacobianInverseMultiVector(params, *mv_input, *mv_result);
00461
00462
00463 result = (*mv_result)[0];
00464
00465 return status;
00466 }
00467
00468 NOX::Abstract::Group::ReturnType
00469 LOCA::Pitchfork::MooreSpence::ExtendedGroup::applyJacobianMultiVector(
00470 const NOX::Abstract::MultiVector& input,
00471 NOX::Abstract::MultiVector& result) const
00472 {
00473 string callingFunction =
00474 "LOCA::Pitchfork::MooreSpence::ExtendedGroup::applyJacobianMultiVector()";
00475 NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
00476 NOX::Abstract::Group::ReturnType status;
00477
00478 if (!isJacobian()) {
00479 globalData->locaErrorCheck->throwError(callingFunction,
00480 "Called with invalid Jacobian!");
00481 }
00482
00483
00484 const LOCA::Pitchfork::MooreSpence::ExtendedMultiVector& pf_input =
00485 dynamic_cast<const LOCA::Pitchfork::MooreSpence::ExtendedMultiVector&>(input);
00486 LOCA::Pitchfork::MooreSpence::ExtendedMultiVector& pf_result =
00487 dynamic_cast<LOCA::Pitchfork::MooreSpence::ExtendedMultiVector&>(result);
00488
00489
00490 Teuchos::RCP<const NOX::Abstract::MultiVector> input_x =
00491 pf_input.getXMultiVec();
00492 Teuchos::RCP<const NOX::Abstract::MultiVector> input_null =
00493 pf_input.getNullMultiVec();
00494 Teuchos::RCP<const NOX::Abstract::MultiVector::DenseMatrix> input_slack = pf_input.getSlacks();
00495 Teuchos::RCP<const NOX::Abstract::MultiVector::DenseMatrix> input_param = pf_input.getBifParams();
00496
00497
00498 Teuchos::RCP<NOX::Abstract::MultiVector> result_x =
00499 pf_result.getXMultiVec();
00500 Teuchos::RCP<NOX::Abstract::MultiVector> result_null =
00501 pf_result.getNullMultiVec();
00502 Teuchos::RCP<NOX::Abstract::MultiVector::DenseMatrix> result_slack =
00503 pf_result.getSlacks();
00504 Teuchos::RCP<NOX::Abstract::MultiVector::DenseMatrix> result_param =
00505 pf_result.getBifParams();
00506
00507
00508 Teuchos::RCP<NOX::Abstract::MultiVector> tmp =
00509 input_null->clone(NOX::ShapeCopy);
00510
00511
00512 if (!grpPtr->isJacobian()) {
00513 status = grpPtr->computeJacobian();
00514 finalStatus =
00515 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00516 finalStatus,
00517 callingFunction);
00518 }
00519
00520
00521 status = grpPtr->applyJacobianMultiVector(*input_x, *result_x);
00522 finalStatus =
00523 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00524 finalStatus,
00525 callingFunction);
00526
00527
00528 result_x->update(Teuchos::NO_TRANS, 1.0, *asymMultiVec, *input_slack);
00529
00530
00531
00532 result_x->update(Teuchos::NO_TRANS, 1.0, *(dfdpMultiVec->getXMultiVec()),
00533 *input_param);
00534
00535
00536 status = grpPtr->applyJacobianMultiVector(*input_null, *result_null);
00537 finalStatus =
00538 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00539 callingFunction);
00540
00541
00542 result_null->update(Teuchos::NO_TRANS, 1.0,
00543 *(dfdpMultiVec->getNullMultiVec()),
00544 *input_param);
00545
00546
00547 status = grpPtr->computeDJnDxaMulti(*(xVec->getNullVec()),
00548 *(fVec->getNullVec()),
00549 *input_x, *tmp);
00550 finalStatus =
00551 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00552 callingFunction);
00553
00554
00555 result_null->update(1.0, *tmp, 1.0);
00556
00557
00558 grpPtr->innerProduct(*asymMultiVec, *input_x, *result_slack);
00559
00560
00561 lTransNorm(*input_null, *result_param);
00562
00563 return finalStatus;
00564 }
00565
00566 NOX::Abstract::Group::ReturnType
00567 LOCA::Pitchfork::MooreSpence::ExtendedGroup::applyJacobianTransposeMultiVector(
00568 const NOX::Abstract::MultiVector& input,
00569 NOX::Abstract::MultiVector& result) const
00570 {
00571 globalData->locaErrorCheck->throwError(
00572 "LOCA::Pitchfork::MooreSpence::ExtendedGroup::applyJacobianTransposeMultiVector()",
00573 "Method not implemented!");
00574
00575 return NOX::Abstract::Group::NotDefined;
00576 }
00577
00578 NOX::Abstract::Group::ReturnType
00579 LOCA::Pitchfork::MooreSpence::ExtendedGroup::applyJacobianInverseMultiVector(
00580 Teuchos::ParameterList& params,
00581 const NOX::Abstract::MultiVector& input,
00582 NOX::Abstract::MultiVector& result) const
00583 {
00584
00585 const LOCA::Pitchfork::MooreSpence::ExtendedMultiVector& pf_input =
00586 dynamic_cast<const LOCA::Pitchfork::MooreSpence::ExtendedMultiVector&>(input);
00587 LOCA::Pitchfork::MooreSpence::ExtendedMultiVector& pf_result =
00588 dynamic_cast<LOCA::Pitchfork::MooreSpence::ExtendedMultiVector&>(result);
00589
00590 NOX::Abstract::Group::ReturnType res = solverStrategy->solve(params, pf_input, pf_result);
00591
00592 return res;
00593 }
00594
00595 bool
00596 LOCA::Pitchfork::MooreSpence::ExtendedGroup::isF() const
00597 {
00598 return isValidF;
00599 }
00600
00601 bool
00602 LOCA::Pitchfork::MooreSpence::ExtendedGroup::isJacobian() const
00603 {
00604 return isValidJacobian;
00605 }
00606
00607 bool
00608 LOCA::Pitchfork::MooreSpence::ExtendedGroup::isGradient() const
00609 {
00610 return false;
00611 }
00612
00613 bool
00614 LOCA::Pitchfork::MooreSpence::ExtendedGroup::isNewton() const
00615 {
00616 return isValidNewton;
00617 }
00618
00619 const NOX::Abstract::Vector&
00620 LOCA::Pitchfork::MooreSpence::ExtendedGroup::getX() const
00621 {
00622 return *xVec;
00623 }
00624
00625 const NOX::Abstract::Vector&
00626 LOCA::Pitchfork::MooreSpence::ExtendedGroup::getF() const
00627 {
00628 return *fVec;
00629 }
00630
00631 double
00632 LOCA::Pitchfork::MooreSpence::ExtendedGroup::getNormF() const
00633 {
00634 return fVec->norm();
00635 }
00636
00637 const NOX::Abstract::Vector&
00638 LOCA::Pitchfork::MooreSpence::ExtendedGroup::getGradient() const
00639 {
00640 globalData->locaErrorCheck->throwError(
00641 "LOCA::Pitchfork::MooreSpence::ExtendedGroup::getGradient()",
00642 " - not implemented");
00643 return getNewton();
00644 }
00645
00646 const NOX::Abstract::Vector&
00647 LOCA::Pitchfork::MooreSpence::ExtendedGroup::getNewton() const
00648 {
00649 return *newtonVec;
00650 }
00651
00652 double
00653 LOCA::Pitchfork::MooreSpence::ExtendedGroup::getNormNewtonSolveResidual() const
00654 {
00655 string callingFunction =
00656 "LOCA::Pitchfork::MooreSpence::ExtendedGroup::getNormNewtonSolveResidual()";
00657 NOX::Abstract::Group::ReturnType finalStatus;
00658 LOCA::Pitchfork::MooreSpence::ExtendedVector residual = *fVec;
00659
00660 finalStatus = applyJacobian(*newtonVec, residual);
00661 globalData->locaErrorCheck->checkReturnType(finalStatus, callingFunction);
00662
00663 residual.update(1.0, *fVec, 1.0);
00664 return residual.norm();
00665 }
00666
00667 Teuchos::RCP<const LOCA::MultiContinuation::AbstractGroup>
00668 LOCA::Pitchfork::MooreSpence::ExtendedGroup::getUnderlyingGroup() const
00669 {
00670 return grpPtr;
00671 }
00672
00673 Teuchos::RCP<LOCA::MultiContinuation::AbstractGroup>
00674 LOCA::Pitchfork::MooreSpence::ExtendedGroup::getUnderlyingGroup()
00675 {
00676 return grpPtr;
00677 }
00678
00679 void
00680 LOCA::Pitchfork::MooreSpence::ExtendedGroup::copy(
00681 const NOX::Abstract::Group& src)
00682 {
00683 const LOCA::Pitchfork::MooreSpence::ExtendedGroup& source =
00684 dynamic_cast<const LOCA::Pitchfork::MooreSpence::ExtendedGroup&>(src);
00685
00686
00687 if (this != &source) {
00688
00689
00690 globalData = source.globalData;
00691 parsedParams = source.parsedParams;
00692 pitchforkParams = source.pitchforkParams;
00693 grpPtr->copy(*(source.grpPtr));
00694 xMultiVec = source.xMultiVec;
00695 fMultiVec = source.fMultiVec;
00696 newtonMultiVec = source.newtonMultiVec;
00697 *asymMultiVec = *source.asymMultiVec;
00698 *lengthMultiVec = *source.lengthMultiVec;
00699 index_f = source.index_f;
00700 index_dfdp = source.index_dfdp;
00701 bifParamID = source.bifParamID;
00702 isValidF = source.isValidF;
00703 isValidJacobian = source.isValidJacobian;
00704 isValidNewton = source.isValidNewton;
00705
00706
00707 setupViews();
00708
00709
00710 solverStrategy =
00711 globalData->locaFactory->createMooreSpencePitchforkSolverStrategy(
00712 parsedParams,
00713 pitchforkParams);
00714 }
00715 }
00716
00717 void
00718 LOCA::Pitchfork::MooreSpence::ExtendedGroup::setParamsMulti(
00719 const vector<int>& paramIDs,
00720 const NOX::Abstract::MultiVector::DenseMatrix& vals)
00721 {
00722 grpPtr->setParamsMulti(paramIDs, vals);
00723 for (unsigned int i=0; i<paramIDs.size(); i++) {
00724 if (paramIDs[i] == bifParamID[0])
00725 setBifParam(vals(i,0));
00726 }
00727 }
00728
00729 void
00730 LOCA::Pitchfork::MooreSpence::ExtendedGroup::setParams(
00731 const LOCA::ParameterVector& p)
00732 {
00733 isValidF = false;
00734 isValidJacobian = false;
00735 isValidNewton = false;
00736
00737 grpPtr->setParams(p);
00738 setBifParam(p[bifParamID[0]]);
00739 }
00740
00741 void
00742 LOCA::Pitchfork::MooreSpence::ExtendedGroup::setParam(int paramID,
00743 double val)
00744 {
00745 if (paramID == bifParamID[0])
00746 setBifParam(val);
00747 else
00748 grpPtr->setParam(paramID, val);
00749 }
00750
00751 void
00752 LOCA::Pitchfork::MooreSpence::ExtendedGroup::setParam(string paramID,
00753 double val)
00754 {
00755 const LOCA::ParameterVector& pVec = grpPtr->getParams();
00756 if (pVec.getIndex(paramID) == bifParamID[0])
00757 setBifParam(val);
00758 else
00759 grpPtr->setParam(paramID, val);
00760 }
00761
00762 const LOCA::ParameterVector&
00763 LOCA::Pitchfork::MooreSpence::ExtendedGroup::getParams() const
00764 {
00765 return grpPtr->getParams();
00766 }
00767
00768 double
00769 LOCA::Pitchfork::MooreSpence::ExtendedGroup::getParam(int paramID) const
00770 {
00771 return grpPtr->getParam(paramID);
00772 }
00773
00774 double
00775 LOCA::Pitchfork::MooreSpence::ExtendedGroup::getParam(string paramID) const
00776 {
00777 return grpPtr->getParam(paramID);
00778 }
00779
00780 NOX::Abstract::Group::ReturnType
00781 LOCA::Pitchfork::MooreSpence::ExtendedGroup::computeDfDpMulti(
00782 const vector<int>& paramIDs,
00783 NOX::Abstract::MultiVector& dfdp,
00784 bool isValid_F)
00785 {
00786 string callingFunction =
00787 "LOCA::Pitchfork::MooreSpence::ExtendedGroup::computeDfDpMulti()";
00788 NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
00789 NOX::Abstract::Group::ReturnType status;
00790
00791
00792 LOCA::Pitchfork::MooreSpence::ExtendedMultiVector& pf_dfdp =
00793 dynamic_cast<LOCA::Pitchfork::MooreSpence::ExtendedMultiVector&>(dfdp);
00794
00795
00796
00797
00798
00799 status = grpPtr->computeDfDpMulti(paramIDs, *pf_dfdp.getXMultiVec(),
00800 false);
00801 finalStatus =
00802 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00803 callingFunction);
00804
00805
00806 double sigma = xVec->getSlack();
00807 pf_dfdp.getColumn(0)->getXVec()->update(sigma, *asymVec, 1.0);
00808
00809
00810 status = grpPtr->computeDJnDpMulti(paramIDs, *(xVec->getNullVec()),
00811 *pf_dfdp.getNullMultiVec(), isValid_F);
00812 finalStatus =
00813 globalData->locaErrorCheck->combineAndCheckReturnTypes(status, finalStatus,
00814 callingFunction);
00815
00816
00817 if (!isValid_F) {
00818 pf_dfdp.getScalar(0,0) = grpPtr->innerProduct(*(xVec->getXVec()),
00819 *asymVec);
00820 pf_dfdp.getScalar(1,0) = lTransNorm(*(xVec->getNullVec()));
00821 }
00822 for (int i=0; i<dfdp.numVectors()-1; i++) {
00823 pf_dfdp.getScalar(0,i+1) = 0.0;
00824 pf_dfdp.getScalar(1,i+1) = 0.0;
00825 }
00826
00827 return finalStatus;
00828 }
00829
00830 void
00831 LOCA::Pitchfork::MooreSpence::ExtendedGroup::preProcessContinuationStep(
00832 LOCA::Abstract::Iterator::StepStatus stepStatus)
00833 {
00834 grpPtr->preProcessContinuationStep(stepStatus);
00835 }
00836
00837 void
00838 LOCA::Pitchfork::MooreSpence::ExtendedGroup::postProcessContinuationStep(
00839 LOCA::Abstract::Iterator::StepStatus stepStatus)
00840 {
00841 grpPtr->postProcessContinuationStep(stepStatus);
00842 }
00843
00844 void
00845 LOCA::Pitchfork::MooreSpence::ExtendedGroup::projectToDraw(
00846 const NOX::Abstract::Vector& x,
00847 double *px) const
00848 {
00849 const LOCA::Pitchfork::MooreSpence::ExtendedVector& mx =
00850 dynamic_cast<const LOCA::Pitchfork::MooreSpence::ExtendedVector&>(x);
00851
00852 grpPtr->projectToDraw(*(mx.getXVec()), px);
00853 px[grpPtr->projectToDrawDimension()] = mx.getBifParam();
00854 }
00855
00856 int
00857 LOCA::Pitchfork::MooreSpence::ExtendedGroup::projectToDrawDimension() const
00858 {
00859 return grpPtr->projectToDrawDimension() + 1;
00860 }
00861
00862 void
00863 LOCA::Pitchfork::MooreSpence::ExtendedGroup::printSolution(
00864 const double conParam) const
00865 {
00866 if (globalData->locaUtils->isPrintType(NOX::Utils::StepperDetails)) {
00867 globalData->locaUtils->out() <<
00868 "LOCA::Pitchfork::MooreSpence::ExtendedGroup::printSolution\n";
00869
00870 globalData->locaUtils->out() << "Pitchfork located at: " <<
00871 globalData->locaUtils->sciformat(conParam) << " " <<
00872 globalData->locaUtils->sciformat(getBifParam()) << std::endl;
00873
00874 globalData->locaUtils->out() <<
00875 "\tSlack variable sigma = " <<
00876 globalData->locaUtils->sciformat(xVec->getSlack()) << std::endl;
00877
00878 globalData->locaUtils->out() <<
00879 "\tPrinting Solution Vector for conParam = " <<
00880 globalData->locaUtils->sciformat(conParam) << std::endl;
00881 }
00882 grpPtr->printSolution(conParam);
00883 if (globalData->locaUtils->isPrintType(NOX::Utils::StepperDetails)) {
00884 globalData->locaUtils->out() <<
00885 "\tPrinting Null Vector for bif param = " <<
00886 globalData->locaUtils->sciformat(getBifParam()) << std::endl;
00887 }
00888 grpPtr->printSolution(*(xVec->getNullVec()), xVec->getBifParam());
00889 }
00890
00891 void
00892 LOCA::Pitchfork::MooreSpence::ExtendedGroup::printSolution(
00893 const NOX::Abstract::Vector& x_,
00894 const double conParam) const
00895 {
00896 const LOCA::Pitchfork::MooreSpence::ExtendedVector& pf_x =
00897 dynamic_cast<const LOCA::Pitchfork::MooreSpence::ExtendedVector&>(x_);
00898
00899 if (globalData->locaUtils->isPrintType(NOX::Utils::StepperDetails)) {
00900 globalData->locaUtils->out() <<
00901 "LOCA::Pitchfork::MooreSpence::ExtendedGroup::printSolution\n";
00902
00903 globalData->locaUtils->out() << "Pitchfork located at: " <<
00904 globalData->locaUtils->sciformat(conParam) << " " <<
00905 globalData->locaUtils->sciformat(pf_x.getBifParam()) << std::endl;
00906
00907 globalData->locaUtils->out() <<
00908 "\tSlack variable sigma = " <<
00909 globalData->locaUtils->sciformat(pf_x.getSlack()) << std::endl;
00910
00911 globalData->locaUtils->out() <<
00912 "\tPrinting Solution Vector for conParam = " <<
00913 globalData->locaUtils->sciformat(conParam) << std::endl;
00914 }
00915 grpPtr->printSolution(*pf_x.getXVec(), conParam);
00916 if (globalData->locaUtils->isPrintType(NOX::Utils::StepperDetails)) {
00917 globalData->locaUtils->out() <<
00918 "\tPrinting Null Vector for bif param = " <<
00919 globalData->locaUtils->sciformat(pf_x.getBifParam()) << std::endl;
00920 }
00921 grpPtr->printSolution(*pf_x.getNullVec(), pf_x.getBifParam());
00922 }
00923
00924 double
00925 LOCA::Pitchfork::MooreSpence::ExtendedGroup::getBifParam() const
00926 {
00927 return grpPtr->getParam(bifParamID[0]);
00928 }
00929
00930 double
00931 LOCA::Pitchfork::MooreSpence::ExtendedGroup::lTransNorm(
00932 const NOX::Abstract::Vector& n) const
00933 {
00934 return lengthVec->innerProduct(n) / lengthVec->length();
00935 }
00936
00937 void
00938 LOCA::Pitchfork::MooreSpence::ExtendedGroup::lTransNorm(
00939 const NOX::Abstract::MultiVector& n,
00940 NOX::Abstract::MultiVector::DenseMatrix& result) const
00941 {
00942 n.multiply(1.0 / lengthVec->length(), *lengthMultiVec, result);
00943 }
00944
00945 void
00946 LOCA::Pitchfork::MooreSpence::ExtendedGroup::setBifParam(double param)
00947 {
00948 grpPtr->setParam(bifParamID[0], param);
00949 xVec->getBifParam() = param;
00950
00951 isValidF = false;
00952 isValidJacobian = false;
00953 isValidNewton = false;
00954 }
00955
00956 void
00957 LOCA::Pitchfork::MooreSpence::ExtendedGroup::setupViews()
00958 {
00959 index_f[0] = 0;
00960 index_dfdp[0] = 1;
00961
00962 xVec = xMultiVec.getColumn(0);
00963 fVec = fMultiVec.getColumn(0);
00964 newtonVec = newtonMultiVec.getColumn(0);
00965 asymVec = Teuchos::rcp(&(*asymMultiVec)[0], false);
00966 lengthVec = Teuchos::rcp(&(*lengthMultiVec)[0],false);
00967
00968 ffMultiVec = Teuchos::rcp_dynamic_cast<LOCA::Pitchfork::MooreSpence::ExtendedMultiVector>(fMultiVec.subView(index_f),true);
00969
00970 dfdpMultiVec = Teuchos::rcp_dynamic_cast<LOCA::Pitchfork::MooreSpence::ExtendedMultiVector>(fMultiVec.subView(index_dfdp),true);
00971
00972 }
00973
00974 void
00975 LOCA::Pitchfork::MooreSpence::ExtendedGroup::init(bool perturbSoln,
00976 double perturbSize)
00977 {
00978 xVec->getBifParam() = getBifParam();
00979
00980
00981 double lVecDotNullVec = lTransNorm(*(xVec->getNullVec()));
00982
00983 if (fabs(lVecDotNullVec) < 1.0e-8) {
00984 globalData->locaErrorCheck->throwError(
00985 "LOCA::Pitchfork::MooreSpence::ExtendedGroup::init()",
00986 "null vector cannot be orthogonal to length-scaling vector: ");
00987 }
00988 if (globalData->locaUtils->isPrintType(NOX::Utils::StepperDetails)) {
00989 globalData->locaUtils->out() <<
00990 "\tIn LOCA::Pitchfork::MooreSpence::ExtendedGroup::init(), " <<
00991 "scaling null vector by:" <<
00992 globalData->locaUtils->sciformat(1.0 / lVecDotNullVec) << std::endl;
00993 }
00994 xVec->getNullVec()->scale(1.0/lVecDotNullVec);
00995
00996
00997 double psi_norm = sqrt( grpPtr->innerProduct(*asymVec, *asymVec) );
00998 if (globalData->locaUtils->isPrintType(NOX::Utils::StepperDetails)) {
00999 globalData->locaUtils->out() <<
01000 "\tIn LOCA::Pitchfork::MooreSpence::ExtendedGroup::init(), " <<
01001 "scaling asymmetric vector by:" <<
01002 globalData->locaUtils->sciformat(1.0 / psi_norm) << std::endl;
01003 }
01004 asymVec->scale(1.0 / psi_norm);
01005
01006 if (perturbSoln) {
01007 if (globalData->locaUtils->isPrintType(NOX::Utils::StepperDetails)) {
01008 globalData->locaUtils->out() <<
01009 "\tIn LOCA::Pitchfork::MooreSpence::ExtendedGroup::init(), " <<
01010 "applying random perturbation to initial solution of size: " <<
01011 globalData->locaUtils->sciformat(perturbSize) << endl;
01012 }
01013 Teuchos::RCP<NOX::Abstract::Vector> perturb =
01014 xVec->getXVec()->clone(NOX::ShapeCopy);
01015 perturb->random();
01016 perturb->scale(*(xVec->getXVec()));
01017 xVec->getXVec()->update(perturbSize, *perturb, 1.0);
01018 grpPtr->setX(*(xVec->getXVec()));
01019 }
01020 }
01021