DENDRO: /net/www/csela/dendro/releases/Dendro/examples/src/drivers/elasticitySolver.C Source File

00001 
00007 #include "mpi.h"
00008 #include "petsc.h"
00009 #include "sys.h"
00010 #include <vector>
00011 #include "TreeNode.h"
00012 #include "parUtils.h"
00013 #include "omg.h"
00014 #include "handleStencils.h"
00015 #include "elasticityJac.h"
00016 #include "omgJac.h"
00017 #include "colors.h"
00018 #include <cstdlib>
00019 #include <cstring>
00020 #include "externVars.h"
00021 #include "dendro.h"
00022 
00023 static char help[] = "Solves static Navier-Lame equations\
00024                       using FEM, MG, DA and Matrix-Free methods.\n\n";
00025 
00026 //Don't want time to be synchronized. Need to check load imbalance.
00027 #ifdef MPI_WTIME_IS_GLOBAL
00028 #undef MPI_WTIME_IS_GLOBAL
00029 #endif
00030 
00031 #ifndef iC
00032 #define iC(fun) {CHKERRQ(fun);}
00033 #endif
00034 
00035 #ifdef PETSC_USE_LOG
00036 int elasticityDiagEvent;
00037 int elasticityMultEvent;
00038 int elasticityFinestDiagEvent;
00039 int elasticityFinestMultEvent;
00040 
00041 int vecMassDiagEvent;
00042 int vecMassMultEvent;
00043 int vecMassFinestDiagEvent;
00044 int vecMassFinestMultEvent;
00045 #endif
00046 
00047 #ifdef PETSC_USE_LOG
00048 //user-defined variables
00049 int Jac1DiagEvent;
00050 int Jac1MultEvent;
00051 int Jac1FinestDiagEvent;
00052 int Jac1FinestMultEvent;
00053 
00054 int Jac2DiagEvent;
00055 int Jac2MultEvent;
00056 int Jac2FinestDiagEvent;
00057 int Jac2FinestMultEvent;
00058 
00059 int Jac3DiagEvent;
00060 int Jac3MultEvent;
00061 int Jac3FinestDiagEvent;
00062 int Jac3FinestMultEvent;
00063 #endif
00064 
00065 double***** LaplacianType1Stencil; 
00066 double**** LaplacianType2Stencil; 
00067 double***** MassType1Stencil; 
00068 double**** MassType2Stencil; 
00069 double****** ShapeFnStencil;
00070 
00071 double**** GradDivType2Stencil; 
00072 
00073 int main(int argc, char ** argv ) {     
00074   int size, rank;
00075   bool incCorner = 1;  
00076   unsigned int numPts;
00077   unsigned int solveU = 0;
00078   unsigned int writeB = 0;  
00079   char Kstr[20];
00080   char pFile[50],bFile[50],uFile[50];
00081   double gSize[3];
00082   unsigned int ptsLen;
00083   unsigned int maxNumPts= 1;
00084   unsigned int dim=3;
00085   unsigned int maxDepth=30;
00086   bool compressLut=true;
00087   double localTime, totalTime;
00088   double startTime, endTime;
00089   DendroIntL locSz, totalSz;
00090   std::vector<ot::TreeNode> linOct, balOct;
00091   std::vector<double> pts;
00092   double mgLoadFac = 2.0;
00093 
00094   PetscInitialize(&argc,&argv,"optionsElasticity",help);
00095   ot::RegisterEvents();
00096 
00097   ot::DAMG_Initialize(MPI_COMM_WORLD);
00098 
00099 #ifdef PETSC_USE_LOG
00100   PetscLogEventRegister(&Jac1DiagEvent,"ODAmatDiag",PETSC_VIEWER_COOKIE);
00101   PetscLogEventRegister(&Jac1MultEvent,"ODAmatMult",PETSC_VIEWER_COOKIE);
00102   PetscLogEventRegister(&Jac1FinestDiagEvent,"ODAmatDiagFinest",PETSC_VIEWER_COOKIE);
00103   PetscLogEventRegister(&Jac1FinestMultEvent,"ODAmatMultFinest",PETSC_VIEWER_COOKIE);
00104 
00105   PetscLogEventRegister(&Jac2DiagEvent,"OMGmatDiag-2",PETSC_VIEWER_COOKIE);
00106   PetscLogEventRegister(&Jac2MultEvent,"OMGmatMult-2",PETSC_VIEWER_COOKIE);
00107   PetscLogEventRegister(&Jac2FinestDiagEvent,"OMGmatDiagFinest-2",PETSC_VIEWER_COOKIE);
00108   PetscLogEventRegister(&Jac2FinestMultEvent,"OMGmatMultFinest-2",PETSC_VIEWER_COOKIE);
00109 
00110   PetscLogEventRegister(&Jac3DiagEvent,"OMGmatDiag-3",PETSC_VIEWER_COOKIE);
00111   PetscLogEventRegister(&Jac3MultEvent,"OMGmatMult-3",PETSC_VIEWER_COOKIE);
00112   PetscLogEventRegister(&Jac3FinestDiagEvent,"OMGmatDiagFinest-3",PETSC_VIEWER_COOKIE);
00113   PetscLogEventRegister(&Jac3FinestMultEvent,"OMGmatMultFinest-3",PETSC_VIEWER_COOKIE);
00114   PetscLogEventRegister(&vecMassDiagEvent,"vMassDiag",PETSC_VIEWER_COOKIE);
00115   PetscLogEventRegister(&vecMassMultEvent,"vMassMult",PETSC_VIEWER_COOKIE);
00116   PetscLogEventRegister(&vecMassFinestDiagEvent,"vMassDiagFinest",PETSC_VIEWER_COOKIE);
00117   PetscLogEventRegister(&vecMassFinestMultEvent,"vMassMultFinest",PETSC_VIEWER_COOKIE);
00118 
00119   PetscLogEventRegister(&elasticityDiagEvent,"elDiag",PETSC_VIEWER_COOKIE);
00120   PetscLogEventRegister(&elasticityMultEvent,"elMult",PETSC_VIEWER_COOKIE);
00121   PetscLogEventRegister(&elasticityFinestDiagEvent,"elDiagFinest",PETSC_VIEWER_COOKIE);
00122   PetscLogEventRegister(&elasticityFinestMultEvent,"elMultFinest",PETSC_VIEWER_COOKIE);
00123 
00124   int stages[3];
00125   PetscLogStageRegister(&stages[0],"P2O.");
00126   PetscLogStageRegister(&stages[1],"Bal");  
00127   PetscLogStageRegister(&stages[2],"Solve");  
00128 #endif
00129 
00130   MPI_Comm_size(MPI_COMM_WORLD,&size);
00131   MPI_Comm_rank(MPI_COMM_WORLD,&rank);
00132   if(argc < 2) {
00133     std::cerr << "Usage: " << argv[0] << "inpfile  maxDepth[30] solveU[0]\
00134       writeB[0] dim[3] maxNumPtsPerOctant[1] incCorner[1] compressLut[1] mgLoadFac[2.0] " << std::endl;
00135     return -1;
00136   }
00137   if(argc > 2) {
00138     maxDepth = atoi(argv[2]);
00139   }
00140   if(argc > 3) {
00141     solveU = atoi(argv[3]);
00142   }
00143   if(argc > 4) {
00144     writeB = atoi(argv[4]);
00145   }
00146   if(argc > 5) {
00147     dim = atoi(argv[5]);
00148   }
00149   if(argc > 6) {
00150     maxNumPts = atoi(argv[6]);
00151   }
00152   if(argc > 7) { incCorner = (bool)(atoi(argv[7]));}  
00153   if(argc > 8) { compressLut = (bool)(atoi(argv[8]));}
00154   if(argc > 9) { mgLoadFac = atof(argv[9]); }
00155 
00156   strcpy(bFile,argv[1]);
00157   ot::int2str(rank,Kstr);
00158   strcat(bFile,Kstr);
00159   strcat(bFile,"_\0");
00160   ot::int2str(size,Kstr);
00161   strcat(bFile,Kstr);
00162   strcpy(pFile,bFile);
00163   strcpy(uFile,bFile);
00164   strcat(bFile,"_Bal.ot\0");
00165   strcat(pFile,".pts\0");
00166   strcat(uFile,".sol\0");
00167 
00168   //Points2Octree....
00169   MPI_Barrier(MPI_COMM_WORLD);  
00170   if(!rank){
00171     std::cout << " reading  "<<pFile<<std::endl; // Point size
00172   }
00173   ot::readPtsFromFile(pFile, pts);
00174   if(!rank){
00175     std::cout << " finished reading  "<<pFile<<std::endl; // Point size
00176   }
00177   MPI_Barrier(MPI_COMM_WORLD);  
00178   ptsLen = pts.size();
00179   std::vector<ot::TreeNode> tmpNodes;
00180   for(int i=0;i<ptsLen;i+=3) {
00181     if( (pts[i] > 0.0) &&
00182         (pts[i+1] > 0.0)  
00183         && (pts[i+2] > 0.0) &&
00184         ( ((unsigned int)(pts[i]*((double)(1u << maxDepth)))) < (1u << maxDepth))  &&
00185         ( ((unsigned int)(pts[i+1]*((double)(1u << maxDepth)))) < (1u << maxDepth))  &&
00186         ( ((unsigned int)(pts[i+2]*((double)(1u << maxDepth)))) < (1u << maxDepth)) ) {
00187       tmpNodes.push_back( ot::TreeNode((unsigned int)(pts[i]*(double)(1u << maxDepth)),
00188             (unsigned int)(pts[i+1]*(double)(1u << maxDepth)),
00189             (unsigned int)(pts[i+2]*(double)(1u << maxDepth)),
00190             maxDepth,dim,maxDepth) );
00191     }
00192   }
00193   pts.clear();
00194   par::removeDuplicates<ot::TreeNode>(tmpNodes,false,MPI_COMM_WORLD);   
00195   linOct = tmpNodes;
00196   tmpNodes.clear();
00197   par::partitionW<ot::TreeNode>(linOct, NULL,MPI_COMM_WORLD);
00198   // reduce and only print the total ...
00199   locSz = linOct.size();
00200   par::Mpi_Reduce<DendroIntL>(&locSz, &totalSz, 1, MPI_SUM, 0, MPI_COMM_WORLD);
00201   MPI_Barrier(MPI_COMM_WORLD);  
00202   if(rank==0) {
00203     std::cout<<"# pts= " << totalSz<<std::endl;
00204   }
00205   MPI_Barrier(MPI_COMM_WORLD);  
00206 
00207   pts.resize(3*(linOct.size()));
00208   ptsLen = (3*(linOct.size()));
00209   for(int i=0;i<linOct.size();i++) {
00210     pts[3*i] = (((double)(linOct[i].getX())) + 0.5)/((double)(1u << maxDepth));
00211     pts[(3*i)+1] = (((double)(linOct[i].getY())) +0.5)/((double)(1u << maxDepth));
00212     pts[(3*i)+2] = (((double)(linOct[i].getZ())) +0.5)/((double)(1u << maxDepth));
00213   }//end for i
00214   linOct.clear();
00215   gSize[0] = 1.;
00216   gSize[1] = 1.;
00217   gSize[2] = 1.;
00218 
00219   MPI_Barrier(MPI_COMM_WORLD);  
00220 #ifdef PETSC_USE_LOG
00221   PetscLogStagePush(stages[0]);
00222 #endif
00223   startTime = MPI_Wtime();
00224   ot::points2Octree(pts, gSize, linOct, dim, maxDepth, maxNumPts, MPI_COMM_WORLD);
00225   endTime = MPI_Wtime();
00226 #ifdef PETSC_USE_LOG
00227   PetscLogStagePop();
00228 #endif
00229   localTime = endTime - startTime;
00230   par::Mpi_Reduce<double>(&localTime, &totalTime, 1, MPI_MAX, 0, MPI_COMM_WORLD);
00231   if(!rank){
00232     std::cout <<"P2n Time: "<<totalTime << std::endl;
00233   }
00234   // reduce and only print the total ...
00235   locSz = linOct.size();
00236   par::Mpi_Reduce<DendroIntL>(&locSz, &totalSz, 1, MPI_SUM, 0, MPI_COMM_WORLD);
00237   if(rank==0) {
00238     std::cout<<"# of Unbalanced Octants: " << totalSz<<std::endl;
00239   }
00240   pts.clear();
00241 
00242   //Balancing...
00243   MPI_Barrier(MPI_COMM_WORLD);  
00244 #ifdef PETSC_USE_LOG
00245   PetscLogStagePush(stages[1]);
00246 #endif
00247   startTime = MPI_Wtime();
00248   ot::balanceOctree (linOct, balOct, dim, maxDepth, incCorner, MPI_COMM_WORLD, NULL, NULL);
00249   endTime = MPI_Wtime();
00250 #ifdef PETSC_USE_LOG
00251   PetscLogStagePop();
00252 #endif
00253   linOct.clear();
00254   if(writeB) { 
00255     ot::writeNodesToFile(bFile,balOct);
00256   }
00257   // compute total inp size and output size
00258   locSz = balOct.size();
00259   localTime = endTime - startTime;
00260   par::Mpi_Reduce<DendroIntL>(&locSz, &totalSz, 1, MPI_SUM, 0, MPI_COMM_WORLD);
00261   par::Mpi_Reduce<double>(&localTime, &totalTime, 1, MPI_MAX, 0, MPI_COMM_WORLD);
00262 
00263   if(!rank) {
00264     std::cout << "# of Balanced Octants: "<< totalSz << std::endl;       
00265     std::cout << "bal Time: "<<totalTime << std::endl;
00266   }
00267 
00268   //Solve ...
00269   MPI_Barrier(MPI_COMM_WORLD);  
00270 #ifdef PETSC_USE_LOG
00271   PetscLogStagePush(stages[2]);
00272 #endif
00273 
00274   ot::DAMG       *damg;    
00275   int       nlevels = 1; //number of multigrid levels
00276   PetscInt       numRefinements = 0;
00277   unsigned int   dof = 3; // degrees of freedom per node  
00278 
00279   iC(PetscOptionsGetInt(0,"-numRefinements",&numRefinements,0));
00280   for(int i = 0; i < numRefinements; i++) {
00281     std::vector<ot::TreeNode> tmpOct = balOct;
00282     balOct.clear();
00283     ot::refineOctree(tmpOct, balOct); 
00284   }
00285 
00286   PetscInt nlevelsPetscInt = nlevels;
00287   PetscOptionsGetInt(0, "-nlevels", &nlevelsPetscInt, 0);
00288   nlevels = nlevelsPetscInt;
00289 
00290   // Note: The user context for all levels will be set separately later.
00291   MPI_Barrier(MPI_COMM_WORLD);  
00292 
00293   if(!rank) {
00294     std::cout<<"nlevels initial: "<<nlevels<<std::endl;
00295   }
00296 
00297   ot::DAMGCreateAndSetDA(PETSC_COMM_WORLD, nlevels, NULL, &damg, 
00298       balOct, dof, mgLoadFac, compressLut, incCorner);
00299 
00300   if(!rank) {
00301     std::cout<<"nlevels final: "<<nlevels<<std::endl;
00302   }
00303 
00304   MPI_Barrier(MPI_COMM_WORLD);  
00305 
00306   if(!rank) {
00307     std::cout << "Created DA for all levels."<< std::endl;
00308   }
00309 
00310   MPI_Barrier(MPI_COMM_WORLD);
00311   ot::PrintDAMG(damg);
00312   MPI_Barrier(MPI_COMM_WORLD);
00313 
00314   MPI_Barrier(MPI_COMM_WORLD);
00315   if(!rank) {
00316     std::cout << "Creating stencils..."<< std::endl;
00317   }
00318   MPI_Barrier(MPI_COMM_WORLD);
00319 
00320   createLmatType2(LaplacianType2Stencil);
00321   createMmatType2(MassType2Stencil);
00322   createGDmatType2(GradDivType2Stencil);
00323 
00324   MPI_Barrier(MPI_COMM_WORLD);
00325   if(!rank) {
00326     std::cout << "Created all stencils."<< std::endl;
00327   }
00328   MPI_Barrier(MPI_COMM_WORLD);
00329 
00330   ot::DAMGCreateSuppressedDOFs(damg);
00331 
00332   SetElasticityContexts(damg);
00333 
00334   MPI_Barrier(MPI_COMM_WORLD);
00335   if(!rank) {
00336     std::cout << "Set Elasticity Contexts all levels."<< std::endl;
00337   }
00338 
00339   //Functions for using KSP_Shell (will be used @ the coarsest grid if not all
00340 //processors are active on the coarsest grid)
00341 
00342   ot::getPrivateMatricesForKSP_Shell = getPrivateMatricesForKSP_Shell_Elas;
00343   
00344   //Set function pointers so that PC_BlockDiag could be used.
00345 
00346   ot::getDofAndNodeSizeForPC_BlockDiag = getDofAndNodeSizeForElasticityMat;
00347 
00348   ot::computeInvBlockDiagEntriesForPC_BlockDiag = computeInvBlockDiagEntriesForElasticityMat;
00349 
00350   PetscInt rhsType = 2;
00351   iC(PetscOptionsGetInt(0,"-rhsType",&rhsType,0));
00352 
00353   if(rhsType == 4) {
00354     ot::DAMGSetKSP(damg, CreateElasticityMat,
00355         ComputeElasticityMat, ComputeRHS4);
00356   }else {
00357     ot::DAMGSetKSP(damg, CreateElasticityMat,
00358         ComputeElasticityMat, ComputeElasticityRHS);
00359   }
00360 
00361   PetscReal norm2;
00362   PetscReal normInf;
00363 
00364   MPI_Barrier(MPI_COMM_WORLD);
00365   if(!rank) {
00366     std::cout << "Solving with 0-intial guess"<< std::endl;
00367   }
00368   MPI_Barrier(MPI_COMM_WORLD);
00369 
00370   startTime = MPI_Wtime();
00371   iC(ot::DAMGSolve(damg));
00372   endTime = MPI_Wtime();
00373   localTime = endTime - startTime;
00374   par::Mpi_Reduce<double>(&localTime, &totalTime, 1, MPI_MAX, 0, MPI_COMM_WORLD);
00375 
00376   MPI_Barrier(MPI_COMM_WORLD);
00377 
00378   if (!rank) {
00379     std::cout << GRN << "Done Solve" << NRM << std::endl;
00380     std::cout << "Solve Time: "<<totalTime << std::endl;
00381   }
00382 
00383   MPI_Barrier(MPI_COMM_WORLD);
00384 
00385   destroyLmatType2(LaplacianType2Stencil);
00386   destroyMmatType2(MassType2Stencil);
00387   destroyGDmatType2(GradDivType2Stencil);
00388 
00389   if (!rank) {
00390     std::cout << GRN << "Destroyed Stencils" << NRM << std::endl;
00391   }
00392 
00393   DestroyElasticityContexts(damg);
00394 
00395   if (!rank) {
00396     std::cout << GRN << "Destroyed Elasticity Contexts." << NRM << std::endl;
00397   }
00398 
00399   MPI_Barrier(MPI_COMM_WORLD);
00400 
00401   iC(DAMGDestroy(damg));
00402 
00403   if (!rank) {
00404     std::cout << GRN << "Destroyed DAMG" << NRM << std::endl;
00405   }
00406 
00407   endTime = MPI_Wtime();
00408 #ifdef PETSC_USE_LOG
00409   PetscLogStagePop();
00410 #endif
00411   balOct.clear();
00412 
00413   if (!rank) {
00414     std::cout << GRN << "Finalizing PETSC" << NRM << std::endl;
00415   }
00416   ot::DAMG_Finalize();
00417   PetscFinalize();
00418 }//end function
00419