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 "NOX.H"
00043 #include "Teuchos_ParameterList.hpp"
00044
00045 #include "NOX_Petsc_Group.H"
00046 #include "NOX_Petsc_Interface.H"
00047 #include "NOX_Petsc_Vector.H"
00048 #include "NOX_Petsc_SharedJacobian.H"
00049 #include "Teuchos_ParameterList.hpp"
00050
00051
00052 #include "petscksp.h"
00053
00054 using namespace NOX;
00055 using namespace NOX::Petsc;
00056
00057 Group::Group(Interface& i, Vec& x, Mat& J) :
00058 xVector(x, "Solution"),
00059 RHSVector(x, "RHS", ShapeCopy),
00060 gradVector(x, "Grad", ShapeCopy),
00061 NewtonVector(x, "Newton", ShapeCopy),
00062 sharedJacobianPtr(new SharedJacobian(J)),
00063 sharedJacobian(*sharedJacobianPtr),
00064 jacType("User Supplied"),
00065 userInterface(i)
00066 {
00067 resetIsValid();
00068 }
00069
00070 Group::Group(const Group& source, CopyType type) :
00071 xVector(source.xVector, type),
00072 RHSVector(source.RHSVector, type),
00073 gradVector(source.gradVector, type),
00074 NewtonVector(source.NewtonVector, type),
00075 sharedJacobianPtr(NULL),
00076 sharedJacobian(source.sharedJacobian),
00077 userInterface(source.userInterface),
00078 jacType(source.jacType)
00079 {
00080 switch (type) {
00081
00082 case DeepCopy:
00083
00084 isValidRHS = source.isValidRHS;
00085 isValidGrad = source.isValidGrad;
00086 isValidNewton = source.isValidNewton;
00087 isValidJacobian = source.isValidJacobian;
00088 isValidPreconditioner = source.isValidPreconditioner;
00089 normRHS = source.normRHS;
00090
00091
00092 if (isValidJacobian)
00093 sharedJacobian.getJacobian(this);
00094
00095 break;
00096
00097 case ShapeCopy:
00098 resetIsValid();
00099 break;
00100
00101 default:
00102 cerr << "ERROR: Invalid ConstructorType for group copy constructor." << endl;
00103 throw "NOX Error";
00104 }
00105
00106 }
00107
00108 Group::~Group()
00109 {
00110 delete sharedJacobianPtr;
00111 }
00112
00113 void Group::resetIsValid()
00114 {
00115 isValidRHS = false;
00116 isValidJacobian = false;
00117 isValidGrad = false;
00118 isValidNewton = false;
00119 isValidPreconditioner = false;
00120 }
00121
00122 Teuchos::RCP<NOX::Abstract::Group>
00123 Group::clone(CopyType type) const
00124 {
00125 Teuchos::RCP<NOX::Abstract::Group> newgrp =
00126 Teuchos::rcp(new NOX::Petsc::Group(*this, type));
00127 return newgrp;
00128 }
00129
00130 Abstract::Group&
00131 Group::operator=(const Abstract::Group& source)
00132 {
00133 return operator=(dynamic_cast<const Group&> (source));
00134 }
00135
00136 Abstract::Group&
00137 Group::operator=(const Group& source)
00138 {
00139
00140 xVector = source.xVector;
00141
00142
00143 sharedJacobian = source.sharedJacobian;
00144
00145
00146 isValidRHS = source.isValidRHS;
00147 isValidGrad = source.isValidGrad;
00148 isValidNewton = source.isValidNewton;
00149 isValidJacobian = source.isValidJacobian;
00150 isValidPreconditioner = source.isValidPreconditioner;
00151
00152
00153 if (isValidRHS) {
00154 RHSVector = source.RHSVector;
00155 normRHS = source.normRHS;
00156 }
00157
00158 if (isValidGrad)
00159 gradVector = source.gradVector;
00160
00161 if (isValidNewton)
00162 NewtonVector = source.NewtonVector;
00163
00164
00165 if (isValidJacobian)
00166 sharedJacobian.getJacobian(this);
00167
00168 return *this;
00169 }
00170
00171 void
00172 Group::setX(const Abstract::Vector& y)
00173 {
00174 setX(dynamic_cast<const Vector&> (y));
00175 return;
00176 }
00177
00178 void
00179 Group::setX(const Vector& y)
00180 {
00181 resetIsValid();
00182 xVector = y;
00183 return;
00184 }
00185
00186 void
00187 Group::computeX(const Abstract::Group& grp, const Abstract::Vector& d, double step)
00188 {
00189
00190 const Group& petscgrp = dynamic_cast<const Group&> (grp);
00191 const Vector& petscd = dynamic_cast<const Vector&> (d);
00192 computeX(petscgrp, petscd, step);
00193 return;
00194 }
00195
00196 void
00197 Group::computeX(const Group& grp, const Vector& d, double step)
00198 {
00199 resetIsValid();
00200 xVector.update(1.0, grp.xVector, step, d);
00201 return;
00202 }
00203
00204 Abstract::Group::ReturnType
00205 Group::computeF()
00206 {
00207 if (isF())
00208 return Abstract::Group::Ok;
00209
00210 bool status = false;
00211
00212 status = userInterface.computeF(xVector.getPetscVector(), RHSVector.getPetscVector());
00213
00214 if(status == false) {
00215 cout << "ERROR: Petsc::Group::computeF() - fill failed!!!"
00216 << endl;
00217 throw "NOX Error: RHS Fill Failed";
00218 }
00219
00220 normRHS = RHSVector.norm();
00221
00222 isValidRHS = true;
00223
00224 return Abstract::Group::Ok;
00225 }
00226
00227 Abstract::Group::ReturnType
00228 Group::computeJacobian()
00229 {
00230
00231 if (isJacobian())
00232 return Abstract::Group::Ok;
00233
00234
00235 Mat& Jacobian = sharedJacobian.getJacobian(this);
00236
00237
00238
00239 bool status = false;
00240
00241 if(jacType == "User Supplied") {
00242
00243 status = userInterface.computeJacobian(xVector.getPetscVector(), Jacobian);
00244
00245 if (status == false) {
00246 cout << "ERROR: Petsc::Group::computeJacobian() - fill failed!!!"
00247 << endl;
00248 throw "NOX Error: Jacobian Fill Failed";
00249 }
00250 }
00251 else if(jacType == "Finite Difference") {
00252 cout << "Finite Difference evaluation not yet supported !!\n\n";
00253 throw "NOX Error";
00254 }
00255
00256
00257 isValidJacobian = true;
00258
00259 return Abstract::Group::Ok;
00260 }
00261
00262 Abstract::Group::ReturnType
00263 Group::computeGradient()
00264 {
00265 if (isGradient())
00266 return Abstract::Group::Ok;
00267
00268 if (!isF()) {
00269 cerr << "ERROR: NOX::Petsc::Group::computeGradient() - RHS is out of date wrt X!" << endl;
00270 throw "NOX Error";
00271 }
00272
00273 if (!isJacobian()) {
00274 cerr << "ERROR: NOX::Petsc::Group::computeGradient() - Jacobian is out of date wrt X!" << endl;
00275 throw "NOX Error";
00276 }
00277
00278
00279 const Mat& Jacobian = sharedJacobian.getJacobian();
00280
00281
00282
00283 int ierr = MatMultTranspose(Jacobian, RHSVector.getPetscVector(), gradVector.getPetscVector());
00284
00285
00286 isValidGrad = true;
00287
00288
00289 return Abstract::Group::Ok;
00290 }
00291
00292 Abstract::Group::ReturnType
00293 Group::computeNewton(Teuchos::ParameterList& p)
00294 {
00295 if (isNewton())
00296 return Abstract::Group::Ok;
00297
00298 if (!isF()) {
00299 cerr << "ERROR: NOX::Petsc::Group::computeNewton() - invalid RHS" << endl;
00300 throw "NOX Error";
00301 }
00302
00303 if (!isJacobian()) {
00304 cerr << "ERROR: NOX::Petsc::Group::computeNewton() - invalid Jacobian" << endl;
00305 throw "NOX Error";
00306 }
00307
00308
00309 Mat& Jacobian = sharedJacobian.getJacobian(this);
00310
00311
00312
00313 KSP ksp;
00314 KSPConvergedReason reason;
00315
00316 int ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);
00317
00318
00319
00320
00321
00322 ierr = KSPSetOperators(ksp,Jacobian,Jacobian, DIFFERENT_NONZERO_PATTERN);
00323
00324
00325
00326
00327
00328 ierr = KSPSetFromOptions(ksp);
00329
00330
00331
00332
00333
00334
00335
00336
00337 ierr = KSPSetUp(ksp);
00338
00339 int its = 0;
00340
00341 ierr = KSPSolve( ksp, RHSVector.getPetscVector(), NewtonVector.getPetscVector() );
00342
00343
00344 ierr = KSPGetConvergedReason(ksp,&reason);
00345
00346 if( reason == KSP_DIVERGED_INDEFINITE_PC )
00347 {
00348 cout << "\nDivergence because of indefinite preconditioner;\n";
00349 cout << "Run the executable again but with -pc_ilu_shift option.\n";
00350 }
00351 else if( reason < 0 )
00352 {
00353 cout << "\nOther kind of divergence: this should not happen.\n";
00354 }
00355 else
00356 {
00357 ierr = KSPGetIterationNumber( ksp, &its );
00358 cout << "\nConvergence in " << its << " iterations.\n";
00359 }
00360
00361
00362 NewtonVector.scale(-1.0);
00363
00364
00365 isValidNewton = true;
00366
00367 return Abstract::Group::Ok;
00368 }
00369
00370 Abstract::Group::ReturnType
00371 Group::applyJacobian(const Abstract::Vector& input, Abstract::Vector& result) const
00372 {
00373 const Vector& petscinput = dynamic_cast<const Vector&> (input);
00374 Vector& petscresult = dynamic_cast<Vector&> (result);
00375 return applyJacobian(petscinput, petscresult);
00376 }
00377
00378 Abstract::Group::ReturnType
00379 Group::applyJacobian(const Vector& input, Vector& result) const
00380 {
00381
00382 if (!isJacobian())
00383 return Abstract::Group::BadDependency;
00384
00385
00386 const Mat& Jacobian = sharedJacobian.getJacobian();
00387
00388
00389 MatMult(Jacobian, input.getPetscVector(), result.getPetscVector());
00390
00391 return Abstract::Group::Ok;
00392 }
00393
00394
00395 Abstract::Group::ReturnType
00396 Group::applyRightPreconditioning(Teuchos::ParameterList& params,
00397 const Abstract::Vector& input,
00398 Abstract::Vector& result) const
00399 {
00400 const Vector& petscinput = dynamic_cast<const Vector&> (input);
00401 Vector& petscresult = dynamic_cast<Vector&> (result);
00402 return applyRightPreconditioning(petscinput, petscresult);
00403 }
00404
00405 Abstract::Group::ReturnType
00406 Group::applyRightPreconditioning(const Vector& input, Vector& result) const
00407 {
00408 if (!isJacobian())
00409 return Abstract::Group::BadDependency;
00410
00411
00412 const Mat& Jacobian = sharedJacobian.getJacobian();
00413
00414
00415 Vec& r = result.getPetscVector();
00416
00417
00418 PC pc;
00419
00420 int ierr = PCCreate(PETSC_COMM_WORLD,&pc);
00421
00422
00423
00424
00425
00426
00427 ierr = PCSetFromOptions(pc);
00428
00429
00430
00431
00432
00433 ierr = PCSetOperators(pc,Jacobian,Jacobian, DIFFERENT_NONZERO_PATTERN);
00434 ierr = PCSetUp(pc);
00435
00436
00437 ierr = PCApply(pc,input.getPetscVector(),r);
00438
00439
00440 ierr = PCDestroy(pc);
00441
00442 return Abstract::Group::Ok;
00443 }
00444
00445
00446 Abstract::Group::ReturnType
00447 Group::applyJacobianTranspose(const Abstract::Vector& input, Abstract::Vector& result) const
00448 {
00449 const Vector& petscinput = dynamic_cast<const Vector&> (input);
00450 Vector& petscresult = dynamic_cast<Vector&> (result);
00451 return applyJacobianTranspose(petscinput, petscresult);
00452 }
00453
00454 Abstract::Group::ReturnType
00455 Group::applyJacobianTranspose(const Vector& input, Vector& result) const
00456 {
00457
00458 if (!isJacobian())
00459 return Abstract::Group::BadDependency;
00460
00461
00462 const Mat& Jacobian = sharedJacobian.getJacobian();
00463
00464
00465 int ierr = MatMultTranspose(Jacobian, input.getPetscVector(), result.getPetscVector());
00466
00467 return Abstract::Group::Ok;
00468 }
00469
00470
00471 bool
00472 Group::isF() const
00473 {
00474 return isValidRHS;
00475 }
00476
00477 bool
00478 Group::isJacobian() const
00479 {
00480 return ((sharedJacobian.isOwner(this)) && (isValidJacobian));
00481 }
00482
00483 bool
00484 Group::isGradient() const
00485 {
00486 return isValidGrad;
00487 }
00488
00489 bool
00490 Group::isNewton() const
00491 {
00492 return isValidNewton;
00493 }
00494
00495 bool
00496 Group::isPreconditioner() const
00497 {
00498 return isValidPreconditioner;
00499 }
00500
00501 const Abstract::Vector&
00502 Group::getX() const
00503 {
00504 return xVector;
00505 }
00506
00507 const Abstract::Vector&
00508 Group::getF() const
00509 {
00510 if (!isF()) {
00511 cerr << "ERROR: NOX::Petsc::Group::getF() - invalid RHS" << endl;
00512 throw "NOX Error";
00513 }
00514
00515 return RHSVector;
00516 }
00517
00518 double
00519 Group::getNormF() const
00520 {
00521 if (!isF()) {
00522 cerr << "ERROR: NOX::Petsc::Group::getNormF() - invalid RHS" << endl;
00523 throw "NOX Error";
00524 }
00525
00526 return normRHS;
00527 }
00528
00529 const Abstract::Vector&
00530 Group::getGradient() const
00531 {
00532 if (!isGradient()) {
00533 cerr << "ERROR: NOX::Petsc::Group::getGradient() - invalid gradient" << endl;
00534 throw "NOX Error";
00535 }
00536
00537 return gradVector;
00538 }
00539
00540 const Abstract::Vector&
00541 Group::getNewton() const
00542 {
00543 if (!isNewton()) {
00544 cerr << "ERROR: NOX::Petsc::Group::getNewton() - invalid Newton vector" << endl;
00545 throw "NOX Error";
00546 }
00547
00548 return NewtonVector;
00549 }
00550
00551