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_MinimallyAugmented_ModifiedConstraint.H"
00043 #include "LOCA_TurningPoint_MinimallyAugmented_AbstractGroup.H"
00044 #include "LOCA_BorderedSolver_AbstractStrategy.H"
00045 #include "LOCA_Parameter_SublistParser.H"
00046 #include "LOCA_GlobalData.H"
00047 #include "LOCA_ErrorCheck.H"
00048 #include "LOCA_Factory.H"
00049 #include "NOX_Utils.H"
00050 #include "Teuchos_ParameterList.hpp"
00051 #include "LOCA_BorderedSolver_JacobianOperator.H"
00052
00053 LOCA::TurningPoint::MinimallyAugmented::ModifiedConstraint::
00054 ModifiedConstraint(
00055 const Teuchos::RCP<LOCA::GlobalData>& global_data,
00056 const Teuchos::RCP<LOCA::Parameter::SublistParser>& topParams,
00057 const Teuchos::RCP<Teuchos::ParameterList>& tpParams,
00058 const Teuchos::RCP<LOCA::TurningPoint::MinimallyAugmented::AbstractGroup>& g,
00059 bool is_symmetric,
00060 const NOX::Abstract::Vector& a,
00061 const NOX::Abstract::Vector* b,
00062 int bif_param) :
00063 Constraint(global_data, topParams, tpParams, g, is_symmetric, a, b,
00064 bif_param),
00065 w_vector_update(a.createMultiVector(1, NOX::ShapeCopy)),
00066 v_vector_update(a.createMultiVector(1, NOX::ShapeCopy)),
00067 w_residual(a.createMultiVector(1, NOX::ShapeCopy)),
00068 v_residual(a.createMultiVector(1, NOX::ShapeCopy)),
00069 deltaX(a.createMultiVector(1, NOX::ShapeCopy)),
00070 sigma1(1, 1),
00071 sigma2(1, 1),
00072 deltaP(0),
00073 isFirstSolve(true),
00074 includeNewtonTerms(false)
00075 {
00076
00077 w_vector_update->init(0.0);
00078 v_vector_update->init(0.0);
00079
00080 includeNewtonTerms = tpParams->get("Include Newton Terms", false);
00081 }
00082
00083 LOCA::TurningPoint::MinimallyAugmented::ModifiedConstraint::
00084 ModifiedConstraint(
00085 const LOCA::TurningPoint::MinimallyAugmented::ModifiedConstraint& source,
00086 NOX::CopyType type) :
00087 Constraint(source, type),
00088 w_vector_update(source.w_vector_update->clone(type)),
00089 v_vector_update(source.v_vector_update->clone(type)),
00090 w_residual(source.w_residual->clone(type)),
00091 v_residual(source.v_residual->clone(type)),
00092 deltaX(source.deltaX->clone(type)),
00093 sigma1(source.sigma1),
00094 sigma2(source.sigma2),
00095 deltaP(source.deltaP),
00096 isFirstSolve(source.isFirstSolve),
00097 includeNewtonTerms(source.includeNewtonTerms)
00098 {
00099 }
00100
00101 LOCA::TurningPoint::MinimallyAugmented::ModifiedConstraint::
00102 ~ModifiedConstraint()
00103 {
00104 }
00105
00106 void
00107 LOCA::TurningPoint::MinimallyAugmented::ModifiedConstraint::
00108 copy(const LOCA::MultiContinuation::ConstraintInterface& src)
00109 {
00110 const LOCA::TurningPoint::MinimallyAugmented::ModifiedConstraint& source =
00111 dynamic_cast<const LOCA::TurningPoint::MinimallyAugmented::ModifiedConstraint&>(src);
00112
00113 if (this != &source) {
00114 Constraint::copy(source);
00115
00116 *w_vector_update = *(source.w_vector_update);
00117 *v_vector_update = *(source.v_vector_update);
00118 *w_residual = *(source.w_residual);
00119 *v_residual = *(source.v_residual);
00120 *deltaX = *(source.deltaX);
00121 sigma1.assign(source.sigma1);
00122 sigma2.assign(source.sigma2);
00123 deltaP = source.deltaP;
00124 isFirstSolve = source.isFirstSolve;
00125 includeNewtonTerms = source.includeNewtonTerms;
00126 }
00127 }
00128
00129 Teuchos::RCP<LOCA::MultiContinuation::ConstraintInterface>
00130 LOCA::TurningPoint::MinimallyAugmented::ModifiedConstraint::
00131 clone(NOX::CopyType type) const
00132 {
00133 return Teuchos::rcp(new ModifiedConstraint(*this, type));
00134 }
00135
00136 NOX::Abstract::Group::ReturnType
00137 LOCA::TurningPoint::MinimallyAugmented::ModifiedConstraint::
00138 computeConstraints()
00139 {
00140 if (isValidConstraints)
00141 return NOX::Abstract::Group::Ok;
00142
00143 string callingFunction =
00144 "LOCA::TurningPoint::MinimallyAugmented::ModifiedConstraint::computeConstraints()";
00145 NOX::Abstract::Group::ReturnType status;
00146 NOX::Abstract::Group::ReturnType finalStatus = NOX::Abstract::Group::Ok;
00147
00148
00149 status = grpPtr->computeJacobian();
00150 finalStatus =
00151 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00152 finalStatus,
00153 callingFunction);
00154
00155
00156 Teuchos::RCP<const LOCA::BorderedSolver::JacobianOperator> op =
00157 Teuchos::rcp(new LOCA::BorderedSolver::JacobianOperator(grpPtr));
00158 borderedSolver->setMatrixBlocksMultiVecConstraint(op,
00159 a_vector,
00160 b_vector,
00161 Teuchos::null);
00162
00163
00164 Teuchos::RCP<Teuchos::ParameterList> linear_solver_params =
00165 parsedParams->getSublist("Linear Solver");
00166
00167
00168 if (isFirstSolve) {
00169
00170 cout << "solving for base w,v..." << endl;
00171
00172
00173 NOX::Abstract::MultiVector::DenseMatrix one(1,1);
00174 one(0,0) = dn;
00175
00176
00177 status = borderedSolver->initForSolve();
00178 finalStatus =
00179 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00180 finalStatus,
00181 callingFunction);
00182 status = borderedSolver->applyInverse(*linear_solver_params,
00183 NULL,
00184 &one,
00185 *v_vector,
00186 sigma1);
00187 finalStatus =
00188 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00189 finalStatus,
00190 callingFunction);
00191
00192
00193 if (!isSymmetric) {
00194 status = borderedSolver->initForTransposeSolve();
00195 finalStatus =
00196 globalData->locaErrorCheck->combineAndCheckReturnTypes(
00197 status,
00198 finalStatus,
00199 callingFunction);
00200 status = borderedSolver->applyInverseTranspose(*linear_solver_params,
00201 NULL,
00202 &one,
00203 *w_vector,
00204 sigma2);
00205 finalStatus =
00206 globalData->locaErrorCheck->combineAndCheckReturnTypes(
00207 status,
00208 finalStatus,
00209 callingFunction);
00210
00211 }
00212 else {
00213 *w_vector = *v_vector;
00214 sigma2.assign(sigma1);
00215 }
00216
00217 isFirstSolve = false;
00218 }
00219
00220
00221 else {
00222
00223 cout << "solving for updates..." << endl;
00224
00225
00226 status = grpPtr->applyJacobianMultiVector(*v_vector, *v_residual);
00227 finalStatus =
00228 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00229 finalStatus,
00230 callingFunction);
00231 v_residual->update(Teuchos::NO_TRANS, 1.0, *a_vector, sigma1, 0.0);
00232
00233
00234 NOX::Abstract::MultiVector::DenseMatrix sigma1_residual(1,1);
00235 v_vector->multiply(1.0, *b_vector, sigma1_residual);
00236 sigma1_residual(0,0) -= dn;
00237
00238 if (includeNewtonTerms) {
00239
00240
00241 Teuchos::RCP<NOX::Abstract::MultiVector> Jv_x_dx =
00242 deltaX->clone(NOX::ShapeCopy);
00243 status = grpPtr->computeDJnDxaMulti((*v_vector)[0], *deltaX, *Jv_x_dx);
00244 finalStatus =
00245 globalData->locaErrorCheck->combineAndCheckReturnTypes(
00246 status,
00247 finalStatus,
00248 callingFunction);
00249
00250
00251 Teuchos::RCP<NOX::Abstract::MultiVector> Jv_p1 =
00252 deltaX->clone(2);
00253 vector<int> idx(1); idx[0] = 0;
00254 Teuchos::RCP<NOX::Abstract::MultiVector> Jv_p =
00255 Jv_p1->subView(idx);
00256 status = grpPtr->computeDJnDpMulti(bifParamID, (*v_vector)[0], *Jv_p1,
00257 false);
00258 finalStatus =
00259 globalData->locaErrorCheck->combineAndCheckReturnTypes(
00260 status,
00261 finalStatus,
00262 callingFunction);
00263
00264 v_residual->update(1.0, *Jv_x_dx, deltaP, *Jv_p, 1.0);
00265
00266
00267 status = grpPtr->computeJacobian();
00268 finalStatus =
00269 globalData->locaErrorCheck->combineAndCheckReturnTypes(
00270 status,
00271 finalStatus,
00272 callingFunction);
00273 }
00274
00275
00276 NOX::Abstract::MultiVector::DenseMatrix sigma1_update(1,1);
00277 status = borderedSolver->initForSolve();
00278 finalStatus =
00279 globalData->locaErrorCheck->combineAndCheckReturnTypes(
00280 status,
00281 finalStatus,
00282 callingFunction);
00283 status = borderedSolver->applyInverse(*linear_solver_params,
00284 v_residual.get(),
00285 &sigma1_residual,
00286 *v_vector_update,
00287 sigma1_update);
00288 finalStatus =
00289 globalData->locaErrorCheck->combineAndCheckReturnTypes(
00290 status,
00291 finalStatus,
00292 callingFunction);
00293
00294
00295 v_vector->update(-1.0, *v_vector_update, 1.0);
00296 sigma1(0,0) -= sigma1_update(0,0);
00297
00298 if (!isSymmetric) {
00299
00300 status = grpPtr->applyJacobianTransposeMultiVector(*w_vector,
00301 *w_residual);
00302 finalStatus =
00303 globalData->locaErrorCheck->combineAndCheckReturnTypes(
00304 status,
00305 finalStatus,
00306 callingFunction);
00307 w_residual->update(Teuchos::NO_TRANS, 1.0, *b_vector, sigma2, 0.0);
00308
00309
00310 NOX::Abstract::MultiVector::DenseMatrix sigma2_residual(1,1);
00311 w_vector->multiply(1.0, *a_vector, sigma2_residual);
00312 sigma2_residual(0,0) -= dn;
00313
00314 if (includeNewtonTerms) {
00315
00316
00317 Teuchos::RCP<NOX::Abstract::MultiVector> Jtw_x_dx =
00318 deltaX->clone(NOX::ShapeCopy);
00319 status = grpPtr->computeDwtJnDx((*w_vector)[0], (*deltaX)[0],
00320 (*Jtw_x_dx)[0]);
00321 finalStatus =
00322 globalData->locaErrorCheck->combineAndCheckReturnTypes(
00323 status,
00324 finalStatus,
00325 callingFunction);
00326
00327
00328 Teuchos::RCP<NOX::Abstract::MultiVector> Jtw_p1 =
00329 deltaX->clone(2);
00330 vector<int> idx(1); idx[0] = 0;
00331 Teuchos::RCP<NOX::Abstract::MultiVector> Jtw_p =
00332 Jtw_p1->subView(idx);
00333 status = grpPtr->computeDwtJDp(bifParamID, (*w_vector)[0], *Jtw_p1,
00334 false);
00335 finalStatus =
00336 globalData->locaErrorCheck->combineAndCheckReturnTypes(
00337 status,
00338 finalStatus,
00339 callingFunction);
00340
00341 w_residual->update(1.0, *Jtw_x_dx, deltaP, *Jtw_p, 1.0);
00342
00343
00344 status = grpPtr->computeJacobian();
00345 finalStatus =
00346 globalData->locaErrorCheck->combineAndCheckReturnTypes(
00347 status,
00348 finalStatus,
00349 callingFunction);
00350 }
00351
00352
00353 NOX::Abstract::MultiVector::DenseMatrix sigma2_update(1,1);
00354 status = borderedSolver->initForTransposeSolve();
00355 finalStatus =
00356 globalData->locaErrorCheck->combineAndCheckReturnTypes(
00357 status,
00358 finalStatus,
00359 callingFunction);
00360 status = borderedSolver->applyInverseTranspose(*linear_solver_params,
00361 w_residual.get(),
00362 &sigma2_residual,
00363 *w_vector_update,
00364 sigma2_update);
00365 finalStatus =
00366 globalData->locaErrorCheck->combineAndCheckReturnTypes(
00367 status,
00368 finalStatus,
00369 callingFunction);
00370
00371
00372 w_vector->update(-1.0, *w_vector_update, 1.0);
00373 sigma2(0,0) -= sigma2_update(0,0);
00374
00375 }
00376 else {
00377 *w_vector = *v_vector;
00378 sigma2.assign(sigma1);
00379 }
00380
00381 }
00382
00383
00384 status = grpPtr->applyJacobianMultiVector(*v_vector, *Jv_vector);
00385 finalStatus =
00386 globalData->locaErrorCheck->combineAndCheckReturnTypes(status,
00387 finalStatus,
00388 callingFunction);
00389 Jv_vector->multiply(-1.0, *w_vector, constraints);
00390
00391
00392
00393
00394
00395 sigma_scale = dn;
00396 constraints.scale(1.0/sigma_scale);
00397
00398 if (globalData->locaUtils->isPrintType(NOX::Utils::OuterIteration)) {
00399 globalData->locaUtils->out() <<
00400 "\n\tEstimate for singularity of Jacobian (sigma1) = " <<
00401 globalData->locaUtils->sciformat(sigma1(0,0));
00402 globalData->locaUtils->out() <<
00403 "\n\tEstimate for singularity of Jacobian (sigma2) = " <<
00404 globalData->locaUtils->sciformat(sigma2(0,0));
00405 globalData->locaUtils->out() <<
00406 "\n\tEstimate for singularity of Jacobian (sigma) = " <<
00407 globalData->locaUtils->sciformat(constraints(0,0)) << std::endl <<
00408 "\tScale factor = " <<
00409 globalData->locaUtils->sciformat(sigma_scale) << std::endl;
00410 }
00411
00412 isValidConstraints = true;
00413
00414
00415 if (updateVectorsEveryIteration) {
00416 if (globalData->locaUtils->isPrintType(NOX::Utils::OuterIteration)) {
00417 globalData->locaUtils->out() <<
00418 "\n\tUpdating null vectors for the next nonlinear iteration" <<
00419 std::endl;
00420 }
00421 *a_vector = *w_vector;
00422 *b_vector = *v_vector;
00423
00424 a_vector->scale(std::sqrt(dn) / (*a_vector)[0].norm());
00425 b_vector->scale(std::sqrt(dn) / (*b_vector)[0].norm());
00426 }
00427
00428 return finalStatus;
00429 }
00430
00431 void
00432 LOCA::TurningPoint::MinimallyAugmented::ModifiedConstraint::
00433 preProcessContinuationStep(LOCA::Abstract::Iterator::StepStatus stepStatus)
00434 {
00435 Constraint::preProcessContinuationStep(stepStatus);
00436
00437
00438 w_vector_update->init(0.0);
00439 v_vector_update->init(0.0);
00440
00441 isFirstSolve = true;
00442 }
00443
00444 void
00445 LOCA::TurningPoint::MinimallyAugmented::ModifiedConstraint::
00446 postProcessContinuationStep(LOCA::Abstract::Iterator::StepStatus stepStatus)
00447 {
00448 Constraint::postProcessContinuationStep(stepStatus);
00449
00450 if (stepStatus == LOCA::Abstract::Iterator::Successful) {
00451
00452
00453 w_vector_update->init(0.0);
00454 v_vector_update->init(0.0);
00455
00456 isFirstSolve = true;
00457 }
00458 }
00459
00460 void
00461 LOCA::TurningPoint::MinimallyAugmented::ModifiedConstraint::
00462 setNewtonUpdates(const NOX::Abstract::Vector& dx, double dp, double step)
00463 {
00464 (*deltaX)[0].update(step, dx, 0.0);
00465 deltaP = step*dp;
00466 }
00467
00468