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_MultiPoint.H"
00040
00041 #ifdef HAVE_NOX_EPETRAEXT
00042
00043 LOCA::Epetra::Interface::MultiPoint::
00044 MultiPoint(
00045 const Teuchos::RCP<LOCA::Epetra::Interface::Required> &iReq_,
00046 const Teuchos::RCP< NOX::Epetra::Interface::Jacobian> &iJac_,
00047 const Epetra_MultiVector &splitMultiVec_,
00048 const Teuchos::RCP<Epetra_RowMatrix> &splitJac_,
00049 const Teuchos::RCP<EpetraExt::MultiComm> &globalComm_) :
00050 iReq(iReq_),
00051 iJac(iJac_),
00052 splitJac(splitJac_),
00053 globalComm(globalComm_),
00054 splitVec(*(splitMultiVec_(0))),
00055 splitRes(*(splitMultiVec_(0))),
00056 jacobian(0),
00057 solution(0),
00058 solutionOverlap(0),
00059 overlapImporter(0),
00060 timeStepsOnTimeDomain(splitMultiVec_.NumVectors()),
00061 numTimeDomains(globalComm_->NumSubDomains()),
00062 timeDomain(globalComm_->SubDomainRank()),
00063 conStep(0),
00064 rowStencil(0),
00065 rowIndex(0)
00066 {
00067
00068 if (globalComm->MyPID()==0) {
00069
00070 cout << "----------MultiPoint Partition Info------------"
00071 << "\n\tNumProcs = " << globalComm->NumProc()
00072 << "\n\tSpatial Decomposition = " << splitMultiVec_.Comm().NumProc()
00073 << "\n\tNumber of Domains = " << numTimeDomains
00074 << "\n\tSteps on Domain 0 = " << timeStepsOnTimeDomain
00075 << "\n\tTotal Number of Steps = " << globalComm->NumTimeSteps();
00076 cout << "\n-----------------------------------------------" << endl;
00077 }
00078
00079
00080
00081
00082 rowStencil = new std::vector< std::vector<int> >(timeStepsOnTimeDomain);
00083 rowIndex = new std::vector<int>;
00084 for (int i=0; i < timeStepsOnTimeDomain; i++) {
00085 (*rowStencil)[i].push_back(0);
00086 (*rowIndex).push_back(i + globalComm->FirstTimeStepOnDomain());
00087 }
00088
00089 jacobian = new EpetraExt::BlockCrsMatrix(*splitJac, *rowStencil,
00090 *rowIndex, *globalComm);
00091
00092
00093
00094 solution = new EpetraExt::BlockVector(splitJac->RowMatrixRowMap(),
00095 jacobian->RowMap());
00096 solutionOverlap = new EpetraExt::BlockVector(splitJac->RowMatrixRowMap(),
00097 jacobian->ColMap());
00098
00099 overlapImporter = new Epetra_Import(solutionOverlap->Map(), solution->Map());
00100
00101
00102
00103 for (int i=0; i < timeStepsOnTimeDomain; i++)
00104 solution->LoadBlockValues(*(splitMultiVec_(i)), (*rowIndex)[i]);
00105 }
00106
00107 LOCA::Epetra::Interface::MultiPoint::
00108 ~MultiPoint()
00109 {
00110 delete solution;
00111 delete solutionOverlap;
00112 delete overlapImporter;
00113 delete jacobian;
00114 delete rowStencil;
00115 delete rowIndex;
00116 }
00117
00118 bool LOCA::Epetra::Interface::MultiPoint::
00119 computeF(const Epetra_Vector& x, Epetra_Vector& F, const FillType fillFlag)
00120 {
00121 bool stat = true;
00122
00123
00124
00125 solution->Epetra_Vector::operator=(x);
00126 solutionOverlap->Import(*solution, *overlapImporter, Insert);
00127
00128 EpetraExt::BlockVector residual(*solution);
00129
00130 for (int i=0; i < timeStepsOnTimeDomain; i++) {
00131
00132 solution->ExtractBlockValues(splitVec, (*rowIndex)[i]);
00133
00134 splitRes.PutScalar(0.0);
00135
00136
00137 iReq->setMultiPointParameter((*rowIndex)[i]);
00138 stat = iReq->computeF(splitVec, splitRes, fillFlag);
00139
00140 residual.LoadBlockValues(splitRes, (*rowIndex)[i]);
00141 }
00142
00143
00144
00145
00146 F = residual;
00147
00148 return stat;
00149 }
00150
00151 bool LOCA::Epetra::Interface::MultiPoint::
00152 computeJacobian(const Epetra_Vector& x, Epetra_Operator& Jac)
00153 {
00154 bool stat = true;
00155 jacobian->PutScalar(0.0);
00156
00157
00158
00159 solution->Epetra_Vector::operator=(x);
00160
00161 for (int i=0; i < timeStepsOnTimeDomain; i++) {
00162
00163 solution->ExtractBlockValues(splitVec, (*rowIndex)[i]);
00164
00165
00166 iReq->setMultiPointParameter((*rowIndex)[i]);
00167 stat = iJac->computeJacobian(splitVec, *splitJac );
00168
00169 jacobian->LoadBlock(*splitJac, i, 0);
00170 }
00171
00172 return stat;
00173 }
00174
00175 void LOCA::Epetra::Interface::MultiPoint::
00176 setParameters(const LOCA::ParameterVector& params)
00177 {
00178 iReq->setParameters(params);
00179 }
00180
00181 void LOCA::Epetra::Interface::MultiPoint::
00182 printSolution(const Epetra_Vector& x, double conParam)
00183 {
00184 solution->Epetra_Vector::operator=(x);
00185
00186
00187 for (int j=0; j<timeDomain; j++) globalComm->Barrier();
00188 for (int i=0; i < timeStepsOnTimeDomain; i++) {
00189 solution->ExtractBlockValues(splitVec, (*rowIndex)[i]);
00190
00191
00192 iReq->dataForPrintSolution(conStep, (*rowIndex)[i],
00193 globalComm->NumTimeSteps());
00194 iReq->printSolution(splitVec, conParam);
00195 }
00196 for (int j=timeDomain; j<numTimeDomains-1; j++) globalComm->Barrier();
00197
00198 conStep++;
00199 }
00200
00201 EpetraExt::BlockVector& LOCA::Epetra::Interface::MultiPoint::
00202 getSolution()
00203 {
00204 return *solution;
00205 }
00206
00207 EpetraExt::BlockCrsMatrix& LOCA::Epetra::Interface::MultiPoint::
00208 getJacobian()
00209 {
00210 return *jacobian;
00211 }
00212
00213 void LOCA::Epetra::Interface::MultiPoint::
00214 throwError(const string& functionName, const string& errorMsg) const
00215 {
00216 cout << "LOCA::Epetra::Interface::MultiPoint::" << functionName
00217 << " - " << errorMsg << endl;
00218 throw "LOCA Error";
00219 }
00220
00221 #endif // HAVE_NOX_EPETRAEXT