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 #include "LOCA_Epetra_xyztPrec.H"
00040 #include "EpetraExt_RowMatrixOut.h"
00041 #include "EpetraExt_VectorOut.h"
00042 #include "EpetraExt_MatrixMatrix.h"
00043 #include "Epetra_Map.h"
00044
00045
00046
00048 LOCA::Epetra::xyztPrec::
00049 xyztPrec(EpetraExt::BlockCrsMatrix &jacobian_,
00050 Epetra_CrsMatrix &splitJac_,
00051 EpetraExt::BlockVector &solution_,
00052 EpetraExt::BlockVector &solutionOverlap_,
00053 Epetra_Import &overlapImporter_,
00054 Teuchos::ParameterList &precPrintParams_,
00055 Teuchos::ParameterList &precLSParams_,
00056 const Teuchos::RCP<EpetraExt::MultiComm> globalComm_) :
00057 jacobian(jacobian_),
00058 splitJac(splitJac_),
00059 solution(solution_),
00060 solutionOverlap(solutionOverlap_),
00061 overlapImporter(overlapImporter_),
00062 printParams(precPrintParams_),
00063 lsParams(precLSParams_),
00064 globalComm(globalComm_),
00065 linSys(std::vector<NOX::Epetra::LinearSystemAztecOO*>(globalComm_->NumTimeStepsOnDomain())),
00066 jacobianBlock(std::vector<Teuchos::RCP<Epetra_CrsMatrix> >(1 + globalComm_->NumTimeStepsOnDomain())),
00067 massBlock(std::vector<Teuchos::RCP<Epetra_CrsMatrix> >(1 + globalComm_->NumTimeStepsOnDomain())),
00068 diagBlockSubdiag(std::vector<Teuchos::RCP<Epetra_Vector> >(globalComm_->NumTimeStepsOnDomain())),
00069 residual(0),
00070 splitVec(0),
00071 splitRes(0),
00072 splitVecOld(0),
00073 isPeriodic(precLSParams_.get("Periodic",false))
00074 {
00075
00076 string prec = lsParams.get("XYZTPreconditioner","None");
00077
00078 Teuchos::RCP<NOX::Epetra::Interface::Required> iReq =
00079 Teuchos::rcp(&((NOX::Epetra::Interface::Required&)*this),false);
00080 Teuchos::RCP<NOX::Epetra::Interface::Jacobian> iJac =
00081 Teuchos::rcp(&((NOX::Epetra::Interface::Jacobian&)*this),false);
00082
00083 if (prec == "Global") {
00084 label = "LOCA::Epetra::xyztPrec::Global";
00085
00086
00087
00088
00089 linSys[0] = new NOX::Epetra::LinearSystemAztecOO(printParams, lsParams,
00090 iReq, iJac, Teuchos::rcp(&jacobian,false), solution);
00091
00092 linSys[0]->setJacobianOperatorForSolve(Teuchos::rcp(&jacobian,false));
00093 }
00094 else if ( prec == "Sequential" || prec == "Parallel" ||
00095 prec == "BlockDiagonal" || prec == "Parareal" || prec == "BDSDT") {
00096 if (prec == "Sequential") {
00097 label = "LOCA::Epetra::xyztPrec::Sequential";
00098
00099
00100 }
00101 else if (prec == "Parallel") {
00102 label = "LOCA::Epetra::xyztPrec::Parallel";
00103
00104
00105 }
00106 else if (prec == "BlockDiagonal") {
00107 label = "LOCA::Epetra::xyztPrec::BlockDiagonal";
00108
00109
00110 }
00111 else if (prec == "BDSDT") {
00112 label = "LOCA::Epetra::xyztPrec::BDSDT";
00113
00114
00115 }
00116 else if (prec == "Parareal") {
00117 label = "LOCA::Epetra::xyztPrec::Parareal";
00118
00119
00120 }
00121
00122
00123 splitVecOld = new Epetra_Vector(splitJac.RowMap());
00124 splitVec = new Epetra_Vector(splitJac.RowMap());
00125 splitRes = new Epetra_Vector(splitJac.RowMap());
00126 splitRes_NEV = new NOX::Epetra::Vector(Teuchos::rcp(splitRes, false),
00127 NOX::Epetra::Vector::CreateView);
00128 splitVec_NEV = new NOX::Epetra::Vector(Teuchos::rcp(splitVec, false),
00129 NOX::Epetra::Vector::CreateView);
00130 residual = new EpetraExt::BlockVector(solution);
00131
00132
00133 int imax = globalComm->NumTimeStepsOnDomain();
00134 if (prec == "Parareal") imax++;
00135
00136 for (int i=0; i < imax; i++) {
00137 jacobianBlock[i] = Teuchos::rcp(new Epetra_CrsMatrix(splitJac));
00138 jacobianBlock[i]->PutScalar(0.0);
00139 if (prec != "BlockDiagonal") {
00140 massBlock[i] = Teuchos::rcp(new Epetra_CrsMatrix(splitJac));
00141 massBlock[i]->PutScalar(0.0);
00142 }
00143 linSys[i] = new NOX::Epetra::LinearSystemAztecOO(printParams, lsParams,
00144 iReq, iJac, jacobianBlock[i], solution);
00145 if (prec == "BDSDT") {
00146 diagBlockSubdiag[i] = Teuchos::rcp(new Epetra_Vector(*splitVec));
00147 }
00148 }
00149
00150 }
00151 else if (prec == "None") {
00152 label = "LOCA::Epetra::xyztPrec::None";
00153
00154
00155 }
00156 else {
00157 string errorMessage = "Option for \"Preconditioner\" is invalid!";
00158 throwError("LOCA::Epetra::xyztPrec::xyztPrec", errorMessage);
00159 }
00160
00161
00162 }
00163
00164 LOCA::Epetra::xyztPrec::
00165 ~xyztPrec()
00166 {
00167 if (splitVec) delete splitVec;
00168 if (splitRes) delete splitRes;
00169 if (splitVecOld) delete splitVecOld;
00170 if (splitRes_NEV) delete splitRes_NEV;
00171 if (splitVec_NEV) delete splitVec_NEV;
00172 if (residual) delete residual;
00173
00174 string prec = lsParams.get("XYZTPreconditioner","None");
00175
00176 int imax=0;
00177 if (prec == "Global")
00178 imax = 1;
00179 else if (prec == "Sequential" || prec == "Parallel" || prec == "BlockDiagonal")
00180 imax = globalComm->NumTimeStepsOnDomain();
00181 else if (prec == "Parareal")
00182 imax = 1 + globalComm->NumTimeStepsOnDomain();
00183
00184 for (int i=0; i < imax; i++) {
00185 linSys[i]->destroyPreconditioner();
00186 delete linSys[i];
00187 }
00188 }
00189
00190
00191 int LOCA::Epetra::xyztPrec::
00192 SetUseTranspose(bool UseTranspose)
00193 {
00194
00195 return false;
00196 }
00197
00198 int LOCA::Epetra::xyztPrec::
00199 Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00200 {
00201
00202 cout << "ERROR: LOCA::Epetra::xyztPrec::Apply() - "
00203 << "method is NOT implemented!! " << endl;
00204 throw "LOCA Error";
00205 return false;
00206 }
00207
00208 int LOCA::Epetra::xyztPrec::
00209 ApplyInverse(const Epetra_MultiVector& input,
00210 Epetra_MultiVector& result) const
00211 {
00212
00213 string prec = lsParams.get("XYZTPreconditioner","None");
00214
00215 if (prec == "None") {
00216 return 0;
00217 }
00218
00219 if (prec == "Global") {
00220
00221 const NOX::Epetra::Vector input_NEV(Teuchos::rcp( const_cast<Epetra_Vector*>(input(0)), false),
00222 NOX::Epetra::Vector::CreateView);
00223 NOX::Epetra::Vector result_NEV(Teuchos::rcp(result(0), false), NOX::Epetra::Vector::CreateView);
00224
00225
00226 bool stat = linSys[0]->applyRightPreconditioning(false, lsParams, input_NEV, result_NEV);
00227
00228
00229 }
00230 else if (prec == "Sequential") {
00231
00232
00233 solution.PutScalar(0.0);
00234 residual->Epetra_Vector::operator=(*input(0));
00235
00236 for (int isd=0; isd < globalComm->NumSubDomains(); isd++) {
00237
00238
00239
00240 solutionOverlap.Import(solution, overlapImporter, Insert);
00241
00242
00243 if ( isd == globalComm->SubDomainRank() ) {
00244
00245 for (int i=0; i < globalComm->NumTimeStepsOnDomain(); i++) {
00246
00247 bool isFirstGlobalTimeStep = (globalComm->FirstTimeStepOnDomain() + i == 0);
00248
00249
00250 linSys[i]->setJacobianOperatorForSolve(jacobianBlock[i]);
00251
00252
00253 residual->ExtractBlockValues(*splitRes, jacobian.RowIndex(i));
00254
00255 if (!isFirstGlobalTimeStep) {
00256
00257
00258
00259 if (i==0)
00260 solutionOverlap.ExtractBlockValues(*splitVecOld, (jacobian.RowIndex(i) - 1));
00261 else
00262 solution.ExtractBlockValues(*splitVecOld, (jacobian.RowIndex(i) - 1));
00263 massBlock[i]->Multiply(false, *splitVecOld, *splitVec);
00264 splitRes->Update(-1.0, *splitVec, 1.0);
00265 }
00266
00267
00268 splitVec->PutScalar(0.0);
00269 bool stat = linSys[i]->applyJacobianInverse(lsParams, *splitRes_NEV, *splitVec_NEV);
00270 solution.LoadBlockValues(*splitVec, jacobian.RowIndex(i));
00271 }
00272 }
00273 }
00274
00275 result = solution;
00276 }
00277 else if (prec == "Parallel") {
00278
00279
00280 solution.PutScalar(0.0);
00281 residual->Epetra_Vector::operator=(*input(0));
00282
00283 for (int i=0; i < globalComm->NumTimeStepsOnDomain(); i++) {
00284
00285
00286 linSys[i]->setJacobianOperatorForSolve(jacobianBlock[i]);
00287
00288
00289 residual->ExtractBlockValues(*splitRes, jacobian.RowIndex(i));
00290
00291 if (i > 0) {
00292
00293
00294 solution.ExtractBlockValues(*splitVecOld, (jacobian.RowIndex(i) - 1));
00295 massBlock[i]->Multiply(false, *splitVecOld, *splitVec);
00296 splitRes->Update(-1.0, *splitVec, 1.0);
00297 }
00298
00299
00300 splitVec->PutScalar(0.0);
00301 bool stat = linSys[i]->applyJacobianInverse(lsParams, *splitRes_NEV, *splitVec_NEV);
00302 solution.LoadBlockValues(*splitVec, jacobian.RowIndex(i));
00303 }
00304
00305 result = solution;
00306 }
00307 else if (prec == "BlockDiagonal") {
00308
00309
00310 solution.PutScalar(0.0);
00311 residual->Epetra_Vector::operator=(*input(0));
00312
00313 for (int i=0; i < globalComm->NumTimeStepsOnDomain(); i++) {
00314
00315
00316 linSys[i]->setJacobianOperatorForSolve(jacobianBlock[i]);
00317
00318
00319 residual->ExtractBlockValues(*splitRes, jacobian.RowIndex(i));
00320
00321
00322 splitVec->PutScalar(0.0);
00323 bool stat = linSys[i]->applyJacobianInverse(lsParams, *splitRes_NEV, *splitVec_NEV);
00324 solution.LoadBlockValues(*splitVec, jacobian.RowIndex(i));
00325 }
00326
00327 result = solution;
00328 }
00329 else if (prec == "BDSDT") {
00330
00331
00332
00333
00334 solution.PutScalar(0.0);
00335 residual->Epetra_Vector::operator=(*input(0));
00336
00337 for (int isd=0; isd < globalComm->NumSubDomains(); isd++) {
00338
00339
00340
00341 solutionOverlap.Import(solution, overlapImporter, Insert);
00342
00343
00344 if ( isd == globalComm->SubDomainRank() ) {
00345
00346 for (int i=0; i < globalComm->NumTimeStepsOnDomain(); i++) {
00347
00348 bool isFirstGlobalTimeStep = (globalComm->FirstTimeStepOnDomain() + i == 0);
00349
00350
00351 residual->ExtractBlockValues(*splitRes, jacobian.RowIndex(i));
00352
00353 if (!isFirstGlobalTimeStep) {
00354
00355
00356
00357 if (i==0)
00358 solutionOverlap.ExtractBlockValues(*splitVecOld, (jacobian.RowIndex(i) - 1));
00359 else
00360 solution.ExtractBlockValues(*splitVecOld, (jacobian.RowIndex(i) - 1));
00361 splitRes->Multiply(-1.0, *diagBlockSubdiag[i], *splitVecOld, 1.0);
00362 }
00363
00364
00365 solution.LoadBlockValues(*splitRes, jacobian.RowIndex(i));
00366 }
00367 }
00368 }
00369
00370
00371
00372
00373
00374 for (int i=0; i < globalComm->NumTimeStepsOnDomain(); i++) {
00375
00376
00377 linSys[i]->setJacobianOperatorForSolve(jacobianBlock[i]);
00378
00379
00380
00381 solution.ExtractBlockValues(*splitRes, jacobian.RowIndex(i));
00382
00383
00384 splitVec->PutScalar(0.0);
00385 bool stat = linSys[i]->applyJacobianInverse(lsParams, *splitRes_NEV, *splitVec_NEV);
00386 solution.LoadBlockValues(*splitVec, jacobian.RowIndex(i));
00387 }
00388
00389 result = solution;
00390 }
00391 else if (prec == "Parareal") {
00392
00393
00394 solution.PutScalar(0.0);
00395 residual->Epetra_Vector::operator=(*input(0));
00396
00397 int N = globalComm->NumTimeStepsOnDomain();
00398
00399
00400 for (int isd=0; isd < globalComm->NumSubDomains(); isd++) {
00401
00402
00403
00404 if (isd>0) solutionOverlap.Import(solution, overlapImporter, Insert);
00405
00406
00407 if ( isd == globalComm->SubDomainRank() &&
00408 ((isd != globalComm->NumSubDomains()-1) && !isPeriodic)) {
00409
00410
00411 linSys[N]->setJacobianOperatorForSolve(jacobianBlock[N]);
00412
00413
00414 residual->ExtractBlockValues(*splitRes, jacobian.RowIndex(N-1));
00415
00416 if (isd > 0) {
00417
00418
00419
00420 solutionOverlap.ExtractBlockValues(*splitVecOld, (jacobian.RowIndex(0) - 1));
00421 massBlock[N]->Multiply(false, *splitVecOld, *splitVec);
00422 splitRes->Update(-1.0, *splitVec, 1.0);
00423 }
00424
00425
00426 splitVec->PutScalar(0.0);
00427 bool stat = linSys[N]->applyJacobianInverse(lsParams, *splitRes_NEV, *splitVec_NEV);
00428 solution.LoadBlockValues(*splitVec, jacobian.RowIndex(N-1));
00429 }
00430 }
00431
00432 solutionOverlap.Import(solution, overlapImporter, Insert);
00433
00434
00435
00436 for (int i=0; i < globalComm->NumTimeStepsOnDomain(); i++) {
00437
00438
00439 linSys[i]->setJacobianOperatorForSolve(jacobianBlock[i]);
00440
00441
00442 residual->ExtractBlockValues(*splitRes, jacobian.RowIndex(i));
00443
00444 bool isFirstGlobalTimeStep = (globalComm->FirstTimeStepOnDomain() + i == 0);
00445 if (!isFirstGlobalTimeStep || isPeriodic) {
00446
00447
00448 if (i==0) {
00449
00450 int offsetToPrevSolution = jacobian.Stencil(i)[0];
00451 cout << "In Parareal, OFFSET " << offsetToPrevSolution << " for time domain " << globalComm->SubDomainRank() << endl;
00452 solutionOverlap.ExtractBlockValues(*splitVecOld, jacobian.RowIndex(i) + offsetToPrevSolution);
00453 }
00454 else
00455 solution.ExtractBlockValues(*splitVecOld, (jacobian.RowIndex(i) - 1));
00456 massBlock[i]->Multiply(false, *splitVecOld, *splitVec);
00457 splitRes->Update(-1.0, *splitVec, 1.0);
00458 }
00459
00460
00461 splitVec->PutScalar(0.0);
00462 bool stat = linSys[i]->applyJacobianInverse(lsParams, *splitRes_NEV, *splitVec_NEV);
00463 solution.LoadBlockValues(*splitVec, jacobian.RowIndex(i));
00464 }
00465
00466 result = solution;
00467 }
00468
00469 return 0;
00470 }
00471
00472 double LOCA::Epetra::xyztPrec::
00473 NormInf() const
00474 {
00475
00476 cout << "ERROR: LOCA::Epetra::xyztPrec::NormInf() - "
00477 << "method is NOT implemented!! " << endl;
00478 throw "LOCA Error";
00479 return 0.0;
00480 }
00481
00482 bool LOCA::Epetra::xyztPrec::
00483 UseTranspose() const
00484 {
00485
00486 cout << "ERROR: LOCA::Epetra::xyztPrec::UseTranspose() - "
00487 << "method is NOT implemented!! " << endl;
00488 throw "LOCA Error";
00489 return false;
00490 }
00491
00492 bool LOCA::Epetra::xyztPrec::
00493 HasNormInf() const
00494 {
00495
00496 return false;
00497 }
00498
00499 const char* LOCA::Epetra::xyztPrec::
00500 Label() const
00501 {
00502 return label.c_str();
00503 }
00504
00505 const Epetra_Comm& LOCA::Epetra::xyztPrec::
00506 Comm() const
00507 {
00508 return jacobian.Comm();
00509 }
00510
00511 const Epetra_Map& LOCA::Epetra::xyztPrec::
00512 OperatorDomainMap () const
00513 {
00514 return jacobian.OperatorDomainMap();
00515 }
00516
00517 const Epetra_Map& LOCA::Epetra::xyztPrec::
00518 OperatorRangeMap () const
00519 {
00520 return jacobian.OperatorRangeMap();
00521 }
00522
00523 bool LOCA::Epetra::xyztPrec::
00524 computeF(const Epetra_Vector&, Epetra_Vector&,
00525 NOX::Epetra::Interface::Required::FillType)
00526 {
00527
00528 return false;
00529 }
00530
00531 bool LOCA::Epetra::xyztPrec::
00532 computeJacobian(const Epetra_Vector&, Epetra_Operator&)
00533 {
00534
00535 return false;
00536 }
00537
00538 bool LOCA::Epetra::xyztPrec::
00539 computePreconditioner(const Epetra_Vector& x,
00540 Epetra_Operator& Prec,
00541 Teuchos::ParameterList* p)
00542 {
00543
00544 string prec = lsParams.get("XYZTPreconditioner","None");
00545 if (prec == "Global") {
00546
00547 linSys[0]->destroyPreconditioner();
00548 linSys[0]->createPreconditioner(x, lsParams, false);
00549 return true;
00550 }
00551 else if (prec == "Sequential" || prec == "Parallel" ||
00552 prec == "BlockDiagonal" || prec == "Parareal" || prec == "BDSDT") {
00553 for (int i=0; i< globalComm->NumTimeStepsOnDomain(); i++ ) {
00554 if (globalComm->FirstTimeStepOnDomain() + i == 0 && !isPeriodic)
00555 jacobian.ExtractBlock(*jacobianBlock[i], 0, 0);
00556 else {
00557 jacobian.ExtractBlock(*jacobianBlock[i], i, 1);
00558 if (prec != "BlockDiagonal") jacobian.ExtractBlock(*massBlock[i], i, 0);
00559 }
00560
00561 linSys[i]->setJacobianOperatorForSolve(jacobianBlock[i]);
00562 linSys[i]->destroyPreconditioner();
00563 linSys[i]->createPreconditioner(x, lsParams, false);
00564 }
00565 if (prec == "Parareal") {
00566 int N = globalComm->NumTimeStepsOnDomain();
00567 jacobian.ExtractBlock(*jacobianBlock[N], N-1, 1);
00568 jacobian.ExtractBlock(*massBlock[N], N-1, 0);
00569
00570
00571 (*massBlock[N]).Scale(1.0 / (double)N);
00572 EpetraExt::MatrixMatrix::Add(*massBlock[N], false, (double)(N-1), *jacobianBlock[N], 1.0);
00573
00574 linSys[N]->setJacobianOperatorForSolve(jacobianBlock[N]);
00575 linSys[N]->destroyPreconditioner();
00576 linSys[N]->createPreconditioner(x, lsParams, false);
00577 }
00578 if (prec == "BDSDT") {
00579 for (int i=0; i< globalComm->NumTimeStepsOnDomain(); i++ ) {
00580 int j=i-1;
00581 if (globalComm->FirstTimeStepOnDomain() + i != 0) {
00582 if (j<0) j=0;
00583 (*jacobianBlock[j]).ExtractDiagonalCopy(*splitVec);
00584 (*massBlock[i]).ExtractDiagonalCopy(*splitRes);
00585 int ierr = (*diagBlockSubdiag[i]).ReciprocalMultiply(1.0, *splitVec, *splitRes, 0.0);
00586 }
00587 }
00588 }
00589 return true;
00590 }
00591 else {
00592
00593 return true;
00594 }
00595 }
00596
00597 void LOCA::Epetra::xyztPrec::
00598 throwError(const string& functionName, const string& errorMsg) const
00599 {
00600 cout << "LOCA::Epetra::xyztPrec" << functionName
00601 << " - " << errorMsg << endl;
00602 throw "LOCA Error";
00603 }
00604