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_Interface_xyzt.H"
00040
00041 #ifdef HAVE_NOX_EPETRAEXT
00042
00043 #include "EpetraExt_MatrixMatrix.h"
00044 #include "Epetra_VbrMatrix.h"
00045
00046 LOCA::Epetra::Interface::xyzt::
00047 xyzt( const Teuchos::RCP<LOCA::Epetra::Interface::TimeDependent> &interface_,
00048 const Epetra_MultiVector &splitMultiVec_,
00049 const Teuchos::RCP<Epetra_RowMatrix> &splitJac_,
00050 const Teuchos::RCP<EpetraExt::MultiComm> &globalComm_,
00051 const Epetra_Vector &initialCondVec_,
00052 double dt_,
00053 Teuchos::ParameterList *precPrintParams_,
00054 Teuchos::ParameterList *precLSParams_) :
00055 interface(interface_),
00056 splitJac(splitJac_),
00057 globalComm(globalComm_),
00058 splitVec(*(splitMultiVec_(0))),
00059 splitRes(*(splitMultiVec_(0))),
00060 splitVecOld(*(splitMultiVec_(0))),
00061 initialCondVec(initialCondVec_),
00062 jacobian(0),
00063 solution(0),
00064 solutionOverlap(0),
00065 overlapImporter(0),
00066 timeStepsOnTimeDomain(splitMultiVec_.NumVectors()),
00067 numTimeDomains(globalComm_->NumSubDomains()),
00068 timeDomain(globalComm_->SubDomainRank()),
00069 conStep(0),
00070 rowStencil(0),
00071 rowIndex(0),
00072 precPrintParams(precPrintParams_),
00073 precLSParams(precLSParams_),
00074 isCrsMatrix(true),
00075 floquetFillFlag(false),
00076 savedSplitMassForFloquet(0),
00077 dt(dt_)
00078 {
00079 if (precLSParams)
00080 isPeriodic = precLSParams_->get("Periodic",false);
00081 else
00082 isPeriodic = false;
00083
00084 if (globalComm->MyPID()==0) {
00085
00086 cout << "--------------XYZT Partition Info---------------"
00087 << "\n\tNumProcs = " << globalComm->NumProc()
00088 << "\n\tSpatial Decomposition = " << splitMultiVec_.Comm().NumProc()
00089 << "\n\tNumber of Time Domains = " << numTimeDomains
00090 << "\n\tTime Steps on Domain 0 = " << timeStepsOnTimeDomain
00091 << "\n\tNumber of Time Steps = " << globalComm->NumTimeSteps();
00092 if (isPeriodic) cout << "\n\t-->Solving for a Periodic Orbit!" ;
00093 cout << "\n-----------------------------------------------" << endl;
00094 }
00095
00096
00097
00098
00099
00100 rowStencil = new std::vector< std::vector<int> >(timeStepsOnTimeDomain);
00101 rowIndex = new std::vector<int>;
00102 for (int i=0; i < timeStepsOnTimeDomain; i++) {
00103 if (timeDomain!=0 || i!=0)
00104 (*rowStencil)[i].push_back(-1);
00105 else if (isPeriodic)
00106 (*rowStencil)[i].push_back(globalComm->NumTimeSteps()-1);
00107 (*rowStencil)[i].push_back(0);
00108 (*rowIndex).push_back(i + globalComm->FirstTimeStepOnDomain());
00109 }
00110
00111 jacobian = new EpetraExt::BlockCrsMatrix(*splitJac, *rowStencil,
00112 *rowIndex, *globalComm);
00113
00114
00115
00116 solution = new EpetraExt::BlockVector(splitJac->RowMatrixRowMap(),
00117 jacobian->RowMap());
00118 solutionOverlap = new EpetraExt::BlockVector(splitJac->RowMatrixRowMap(),
00119 jacobian->ColMap());
00120
00121 overlapImporter = new Epetra_Import(solutionOverlap->Map(), solution->Map());
00122
00123
00124
00125 for (int i=0; i < timeStepsOnTimeDomain; i++)
00126 solution->LoadBlockValues(*(splitMultiVec_(i)), (*rowIndex)[i]);
00127
00128
00129 if (precLSParams != 0) {
00130
00131
00132 splitJacCrs = dynamic_cast<Epetra_CrsMatrix *>(splitJac.get());
00133 if (splitJacCrs == NULL) {
00134 isCrsMatrix = false;
00135 cout << "CAST OF splitJacCrs failed!, constructing CRS matrix " << endl;
00136
00137 std::vector< std::vector<int> > row(1); row[0].push_back(0);
00138 std::vector<int> col; col.push_back(0);
00139 splitJacCrs = (Epetra_CrsMatrix *)
00140 new EpetraExt::BlockCrsMatrix(*splitJac, row, col,
00141 splitJac->Comm());
00142 }
00143
00144 preconditioner =
00145 new LOCA::Epetra::xyztPrec(*jacobian, *splitJacCrs, *solution,
00146 *solutionOverlap, *overlapImporter,
00147 *precPrintParams, *precLSParams, globalComm);
00148 }
00149 }
00150
00151 LOCA::Epetra::Interface::xyzt::
00152 ~xyzt()
00153 {
00154 delete solution;
00155 delete solutionOverlap;
00156 delete overlapImporter;
00157 delete jacobian;
00158 delete rowStencil;
00159 delete rowIndex;
00160 }
00161
00162 bool LOCA::Epetra::Interface::xyzt::
00163 computeF(const Epetra_Vector& x, Epetra_Vector& F, const FillType fillFlag)
00164 {
00165 bool stat = true;
00166
00167
00168
00169 solution->Epetra_Vector::operator=(x);
00170 solutionOverlap->Import(*solution, *overlapImporter, Insert);
00171
00172 EpetraExt::BlockVector residual(*solution);
00173
00174 for (int i=0; i < timeStepsOnTimeDomain; i++) {
00175
00176 solution->ExtractBlockValues(splitVec, (*rowIndex)[i]);
00177
00178
00179 if (i==0 && timeDomain==0 && !isPeriodic) {
00180 splitVecOld = initialCondVec;
00181 }
00182 else {
00183 int blockRowOld = (*rowIndex)[i] + (*rowStencil)[i][0];
00184 solutionOverlap->ExtractBlockValues(splitVecOld, blockRowOld);
00185 }
00186
00187
00188 Epetra_Vector& vecDot = splitRes;
00189 vecDot.Update(1.0/dt, splitVec, -1.0/dt, splitVecOld, 0.0);
00190 interface->setXdot(vecDot, dt*(i+globalComm->FirstTimeStepOnDomain()));
00191
00192 splitRes.PutScalar(0.0);
00193 stat = stat && interface->computeF(splitVec, splitRes, fillFlag);
00194
00195 residual.LoadBlockValues(splitRes, (*rowIndex)[i]);
00196 }
00197
00198
00199
00200
00201 F = residual;
00202
00203 return stat;
00204 }
00205
00206 bool LOCA::Epetra::Interface::xyzt::
00207 computeJacobian(const Epetra_Vector& x, Epetra_Operator& Jac)
00208 {
00209 bool stat = true;
00210 jacobian->PutScalar(0.0);
00211
00212
00213
00214 solution->Epetra_Vector::operator=(x);
00215 solutionOverlap->Import(*solution, *overlapImporter, Insert);
00216
00217 for (int i=0; i < timeStepsOnTimeDomain; i++) {
00218
00219 solution->ExtractBlockValues(splitVec, (*rowIndex)[i]);
00220
00221
00222 if (i==0 && timeDomain==0 && !isPeriodic) {
00223 splitVecOld = initialCondVec;
00224 }
00225 else {
00226 int blockRowOld = (*rowIndex)[i] + (*rowStencil)[i][0];
00227 solutionOverlap->ExtractBlockValues(splitVecOld, blockRowOld);
00228 }
00229
00230
00231 Epetra_Vector& vecDot = splitRes;
00232 vecDot.Update(1.0/dt, splitVec, -1.0/dt, splitVecOld, 0.0);
00233 interface->setXdot(vecDot, dt*(i+globalComm->FirstTimeStepOnDomain()));
00234
00235 stat = stat && interface->computeShiftedMatrix(1.0, 1.0/dt, splitVec, *splitJac );
00236
00237
00238 if (i==0 && timeDomain==0 && !isPeriodic) {
00239 jacobian->LoadBlock(*splitJac, i, 0);
00240 }
00241 else {
00242 jacobian->LoadBlock(*splitJac, i, 1);
00243 stat = stat && interface->computeShiftedMatrix(0.0, -1.0/dt, splitVecOld, *splitJac );
00244
00245 if (i==0 && timeDomain==0 && isPeriodic && floquetFillFlag)
00246 (*savedSplitMassForFloquet) = *(splitJac.get());
00247 else jacobian->LoadBlock(*splitJac, i, 0);
00248 }
00249 }
00250
00251 return stat;
00252 }
00253
00254 void LOCA::Epetra::Interface::xyzt::
00255 setParameters(const LOCA::ParameterVector& params)
00256 {
00257 interface->setParameters(params);
00258 }
00259
00260 void LOCA::Epetra::Interface::xyzt::
00261 printSolution(const Epetra_Vector& x, double conParam)
00262 {
00263 solution->Epetra_Vector::operator=(x);
00264
00265
00266 for (int j=0; j<timeDomain; j++) globalComm->Barrier();
00267 for (int i=0; i < timeStepsOnTimeDomain; i++) {
00268 solution->ExtractBlockValues(splitVec, (*rowIndex)[i]);
00269
00270
00271 interface->dataForPrintSolution(conStep, (*rowIndex)[i],
00272 globalComm->NumTimeSteps());
00273 interface->printSolution(splitVec, conParam);
00274 }
00275 for (int j=timeDomain; j<numTimeDomains-1; j++) globalComm->Barrier();
00276
00277 conStep++;
00278 }
00279
00280 void LOCA::Epetra::Interface::xyzt::
00281 setFloquetFillFlag(bool fff)
00282 {
00283 if (fff) {
00284 floquetFillFlag = true;
00285
00286 if (isCrsMatrix)
00287 savedSplitMassForFloquet = new Epetra_CrsMatrix(dynamic_cast<Epetra_CrsMatrix&>(*(splitJac.get())));
00288 else
00289 savedSplitMassForFloquet = new Epetra_VbrMatrix(dynamic_cast<Epetra_VbrMatrix&>(*(splitJac.get())));
00290 }
00291 else {
00292 floquetFillFlag = false;
00293
00294 delete savedSplitMassForFloquet;
00295 }
00296 }
00297
00298 void LOCA::Epetra::Interface::xyzt::
00299 beginFloquetOperatorApplication(Epetra_Vector& v)
00300 {
00301 if (!isPeriodic)
00302 cout << "\n\n\t Must be periodic for Floquet theory\n" << endl;
00303
00304
00305
00306
00307 solution->Epetra_Vector::operator=(v);
00308 solutionOverlap->Import(*solution, *overlapImporter, Insert);
00309
00310
00311 solution->PutScalar(0.0);
00312
00313
00314 if (timeDomain == 0) {
00315 int blockRowOld = (*rowIndex)[0] + (*rowStencil)[0][0];
00316 solutionOverlap->ExtractBlockValues(splitVecOld, blockRowOld);
00317 savedSplitMassForFloquet->Multiply(false, splitVecOld, splitVec);
00318 splitVec.Scale(-1.0);
00319 solution->LoadBlockValues(splitVec, 0);
00320 }
00321
00322 v = *solution;
00323 }
00324
00325 void LOCA::Epetra::Interface::xyzt::
00326 finishFloquetOperatorApplication(Epetra_Vector& v)
00327 {
00328
00329 solution->Epetra_Vector::operator=(v);
00330
00331 if ( timeDomain == numTimeDomains-1 )
00332 solution->ExtractBlockValues(splitVec, globalComm->NumTimeSteps()-1);
00333
00334 solution->PutScalar(0.0);
00335
00336 if ( timeDomain == numTimeDomains-1 )
00337 solution->LoadBlockValues(splitVec, globalComm->NumTimeSteps()-1);
00338
00339 v = *solution;
00340 }
00341
00342 EpetraExt::BlockVector& LOCA::Epetra::Interface::xyzt::
00343 getSolution()
00344 {
00345 return *solution;
00346 }
00347
00348 EpetraExt::BlockCrsMatrix& LOCA::Epetra::Interface::xyzt::
00349 getJacobian()
00350 {
00351 return *jacobian;
00352 }
00353
00354 LOCA::Epetra::xyztPrec& LOCA::Epetra::Interface::xyzt::
00355 getPreconditioner()
00356 {
00357 return *preconditioner;
00358 }
00359
00360 void LOCA::Epetra::Interface::xyzt::
00361 throwError(const string& functionName, const string& errorMsg) const
00362 {
00363 cout << "LOCA::Epetra::Interface::xyzt::" << functionName
00364 << " - " << errorMsg << endl;
00365 throw "LOCA Error";
00366 }
00367
00368 #endif // HAVE_NOX_EPETRAEXT