OpenSees » SRC » system_of_eqn » linearSOE » petsc
PetscSolver.cpp

00001 /* ****************************************************************** **
00002 **    OpenSees - Open System for Earthquake Engineering Simulation    **
00003 **          Pacific Earthquake Engineering Research Center            **
00004 **                                                                    **
00005 **                                                                    **
00006 ** (C) Copyright 1999, The Regents of the University of California    **
00007 ** All Rights Reserved.                                               **
00008 **                                                                    **
00009 ** Commercial use of this program without express permission of the   **
00010 ** University of California, Berkeley, is strictly prohibited.  See   **
00011 ** file 'COPYRIGHT'  in main directory for information on usage and   **
00012 ** redistribution,  and for a DISCLAIMER OF ALL WARRANTIES.           **
00013 **                                                                    **
00014 ** Developed by:                                                      **
00015 **   Frank McKenna (fmckenna@ce.berkeley.edu)                         **
00016 **   Gregory L. Fenves (fenves@ce.berkeley.edu)                       **
00017 **   Filip C. Filippou (filippou@ce.berkeley.edu)                     **
00018 **                                                                    **
00019 ** ****************************************************************** */
00020                                                                         
00021 // $Revision: 1.5 $
00022 // $Date: 2006/01/13 00:02:03 $
00023 // $Source: /usr/local/cvs/OpenSees/SRC/system_of_eqn/linearSOE/petsc/PetscSolver.cpp,v $
00024                                                                         
00025 // Written: fmk & om
00026 // Created: 7/98
00027 // Revision: A
00028 //
00029 // Description: This file contains the class implementation for 
00030 // FullGenLinLapackSolver. It solves the FullGenLinSOE object by calling
00031 // Lapack routines.
00032 //
00033 // What: "@(#) FullGenLinLapackSolver.h, revA"
00034 
00035 #include <PetscSolver.h>
00036 #include <PetscSOE.h>
00037 #include <f2c.h>
00038 #include <math.h>
00039 #include <Channel.h>
00040 #include <FEM_ObjectBroker.h>
00041 #include <Timer.h>
00042 #include <ID.h>
00043 #include <Vector.h>
00044 
00045 int PetscSolver::numSolver=0;
00046 
00047 PetscSolver::PetscSolver()
00048   :LinearSOESolver(SOLVER_TAGS_PetscSolver), 
00049    rTol(PETSC_DEFAULT), aTol(PETSC_DEFAULT), dTol(PETSC_DEFAULT), maxIts(PETSC_DEFAULT)
00050 {
00051   numSolver++;
00052 }
00053 
00054 PetscSolver::PetscSolver(KSPType meth, PCType pre)
00055   :LinearSOESolver(SOLVER_TAGS_PetscSolver), method(meth), preconditioner(pre),
00056    rTol(PETSC_DEFAULT), aTol(PETSC_DEFAULT), dTol(PETSC_DEFAULT), maxIts(PETSC_DEFAULT)
00057 {
00058   numSolver++;
00059 }
00060 
00061 PetscSolver::PetscSolver(KSPType meth, PCType pre, double relTol, double absTol, double divTol, int maxIterations)
00062   :LinearSOESolver(SOLVER_TAGS_PetscSolver), method(meth), preconditioner(pre),
00063    rTol(relTol), aTol(absTol), dTol(divTol), maxIts(maxIterations)
00064 {
00065   numSolver++;
00066 }
00067 
00068 PetscSolver::~PetscSolver()
00069 {
00070 
00071   KSPDestroy(ksp);
00072   numSolver--;
00073   /*
00074   if (numSolver == 0)
00075     PetscFinalize();
00076   */
00077 }
00078 
00079 
00080 
00081 int
00082 PetscSolver::solve(void)
00083 {
00084   int size = theSOE->size;
00085   int factored = theSOE->isFactored;
00086   
00087   int numProcesses = theSOE->numProcesses;
00088   int processID = theSOE->processID;
00089   
00090   int ierr;
00091   ierr = MatAssemblyBegin(theSOE->A, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 
00092   ierr = MatAssemblyEnd(theSOE->A, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 
00093 
00094   //
00095   // if parallel, we must zero X & form the total B: each processor has own contributions
00096   //
00097   
00098   static Vector recvVector(1);
00099 
00100   if (numProcesses > 1) {
00101     Vector *vectX = theSOE->vectX;
00102     Vector *vectB = theSOE->vectB;
00103 
00104     // zero X
00105     vectX->Zero();
00106 
00107     //
00108     // form B on each
00109     //
00110 
00111     int numChannels = theSOE->numChannels;
00112     Channel **theChannels = theSOE->theChannels;
00113     
00114     if (processID != 0) {
00115       Channel *theChannel = theChannels[0];
00116       
00117       theChannel->sendVector(0, 0, *vectB);
00118       theChannel->recvVector(0, 0, *vectB);
00119       
00120     } else {
00121       
00122       if (recvVector.Size() != size) 
00123         recvVector.resize(size);
00124       for (int j=0; j<numChannels; j++) {
00125         Channel *theChannel = theChannels[j];
00126         theChannel->recvVector(0, 0, recvVector);
00127         *vectB += recvVector;
00128       }
00129       for (int j=0; j<numChannels; j++) {
00130         Channel *theChannel = theChannels[j];
00131         theChannel->sendVector(0, 0, *vectB);
00132       }
00133     }
00134   }
00135 
00136   //
00137   // solve and mark as having been solved
00138   //
00139 
00140   ierr = KSPSolve(ksp, theSOE->b, theSOE->x); CHKERRQ(ierr); 
00141   theSOE->isFactored = 1;
00142 
00143   //
00144   // if parallel, we must form the total X: each processor has startRow through endRow-1
00145   //
00146 
00147   if (numProcesses > 1) {
00148     Vector *vectX = theSOE->vectX;
00149 
00150     int numChannels = theSOE->numChannels;
00151     Channel **theChannels = theSOE->theChannels;
00152     
00153     if (processID != 0) {
00154       Channel *theChannel = theChannels[0];
00155       
00156       theChannel->sendVector(0, 0, *vectX);
00157       theChannel->recvVector(0, 0, *vectX);
00158       
00159     } else {
00160       
00161       if (recvVector.Size() != size) 
00162         recvVector.resize(size);
00163 
00164       for (int j=0; j<numChannels; j++) {
00165         Channel *theChannel = theChannels[j];
00166         theChannel->recvVector(0, 0, recvVector);
00167         *vectX += recvVector;
00168       }
00169       for (int j=0; j<numChannels; j++) {
00170         Channel *theChannel = theChannels[j];
00171         theChannel->sendVector(0, 0, *vectX);
00172       }
00173     }
00174   }
00175 
00176   return ierr;
00177 }
00178 
00179 
00180 int
00181 PetscSolver::setSize()
00182 {
00183   /* 
00184    * Create linear solver context
00185    */
00186   
00187    KSPCreate(PETSC_COMM_WORLD, &ksp);
00188 
00189    /* 
00190     *  Set operators. NOTE: matrix that defines the linear system
00191     *  also serves as the preconditioning matrix.
00192     */
00193 
00194    KSPSetOperators(ksp, theSOE->A, theSOE->A, DIFFERENT_NONZERO_PATTERN);
00195 
00196    /* 
00197     *  Set solution scheme & tolerances
00198     */
00199    
00200    int ierr;
00201    ierr = KSPSetType(ksp, method); CHKERRQ(ierr); 
00202    ierr = KSPSetTolerances(ksp, rTol, aTol, dTol, maxIts); 
00203 
00204    /* 
00205     *  Set preconditioning scheme
00206     */
00207 
00208    KSPGetPC(ksp, &pc);
00209    ierr = PCSetType(pc,  preconditioner); CHKERRQ(ierr); 
00210 
00211    /* 
00212     *  Finally mark so that uses last solution as initial guess in each solution
00213     *    NOTE: maybe change this as another user supplied option
00214     */
00215 
00216    //   ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); CHKERRQ(ierr); 
00217    return ierr;
00218 }
00219 
00220 
00221 int 
00222 PetscSolver::setLinearSOE(PetscSOE &theSys)
00223 {
00224   theSOE = &theSys;
00225   return 0;
00226 }
00227 
00228 
00229 int
00230 PetscSolver::sendSelf(int cTag, Channel &theChannel)
00231 {
00232   static ID idData(3);
00233   idData(0) = maxIts;
00234   if (strcmp(method, KSPCG) == 0) 
00235     idData(1) = 0;
00236   else if (strcmp(method, KSPBICG) == 0)
00237     idData(1) = 1;
00238   else if (strcmp(method, KSPRICHARDSON) == 0)
00239     idData(1) = 2;
00240   else if (strcmp(method, KSPCHEBYCHEV) == 0)
00241     idData(1) = 3;
00242   else if (strcmp(method, KSPGMRES) ==0) 
00243     idData(1) = 4;
00244   else {
00245     opserr << "PetscSolver::sendSelf() - unknown method set\n";
00246     return -1;
00247   }
00248 
00249   if (strcmp(preconditioner, PCJACOBI) == 0)
00250     idData(2) = 0;
00251   else if (strcmp(preconditioner, PCILU) == 0)
00252     idData(2) = 1;
00253   else if (strcmp(preconditioner, PCICC) == 0)
00254     idData(2) = 2;
00255   else if (strcmp(preconditioner, PCBJACOBI) == 0)
00256     idData(2) = 3;
00257   else if (strcmp(preconditioner, PCNONE) == 0)
00258     idData(2) = 4;
00259   else {
00260     opserr << "PetscSolver::sendSelf() - unknown preconditioner set\n";
00261     return -1;
00262   }
00263 
00264   theChannel.sendID(0, cTag, idData);
00265 
00266   static Vector data(3);
00267   data(0) = rTol;
00268   data(1) = aTol;
00269   data(2) = dTol;
00270 
00271   theChannel.sendVector(0, cTag, data);
00272   return 0;
00273 }
00274 
00275 int
00276 PetscSolver::recvSelf(int cTag, Channel &theChannel, 
00277                       FEM_ObjectBroker &theBroker)
00278 {
00279   static ID idData(3);
00280   theChannel.recvID(0, cTag, idData);
00281   maxIts = idData(0);
00282   if (idData(1) == 0) 
00283     method = KSPCG;
00284   else if (idData(1) == 1)
00285     method = KSPBICG;
00286   else if (idData(1) == 2) 
00287     method = KSPRICHARDSON;
00288   else if (idData(1) == 3)
00289     method = KSPCHEBYCHEV;
00290   else if (idData(1) == 4)
00291     method = KSPGMRES;
00292   else {
00293     opserr << "PetscSolver::recvSelf() - unknown method recvd\n";
00294     return -1;
00295   }
00296 
00297   if (idData(2) == 0)
00298     preconditioner = PCJACOBI;
00299   else if (idData(2) == 1)
00300     preconditioner = PCILU;
00301   else if (idData(2) == 2)
00302     preconditioner = PCICC;
00303   else if (idData(2) == 3)
00304     preconditioner = PCBJACOBI;
00305   else if (idData(2) == 4)
00306     preconditioner = PCNONE;
00307   else {
00308     opserr << "PetscSolver::sendSelf() - unknown preconditioner set\n";
00309     return -1;
00310   }
00311 
00312 
00313   static Vector data(3);
00314   theChannel.recvVector(0, cTag, data);
00315   rTol = data(0);
00316   aTol = data(1);
00317   dTol = data(2);
00318 
00319   return 0;
00320 }
00321 
00322 
00323
Generated on Mon Oct 23 15:05:29 2006 for OpenSees by 1.5.0
opensees-support @ berkeley.edu
Supported by the National Science Foundation