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