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