DENDRO: /net/www/csela/dendro/releases/Dendro/examples/src/drivers/pts2Mesh.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 "odaUtils.h"
00015 #include "handleStencils.h"
00016 #include "elasticityJac.h"
00017 #include "omgJac.h"
00018 #include <cstdlib>
00019 #include <cstring>
00020 #include "colors.h"
00021 #include "externVars.h"
00022 #include "dendro.h"
00023 
00024 static char help[] = "Solves static Navier-Lame equations\
00025                       using FEM, MG, DA and Matrix-Free methods.\n\n";
00026 
00027 //Don't want time to be synchronized. Need to check load imbalance.
00028 #ifdef MPI_WTIME_IS_GLOBAL
00029 #undef MPI_WTIME_IS_GLOBAL
00030 #endif
00031 
00032 #ifndef iC
00033 #define iC(fun) {CHKERRQ(fun);}
00034 #endif
00035 
00036 #ifdef PETSC_USE_LOG
00037 int elasticityDiagEvent;
00038 int elasticityMultEvent;
00039 int elasticityFinestDiagEvent;
00040 int elasticityFinestMultEvent;
00041 
00042 int vecMassDiagEvent;
00043 int vecMassMultEvent;
00044 int vecMassFinestDiagEvent;
00045 int vecMassFinestMultEvent;
00046 #endif
00047 
00048 #ifdef PETSC_USE_LOG
00049 //user-defined variables
00050 int Jac1DiagEvent;
00051 int Jac1MultEvent;
00052 int Jac1FinestDiagEvent;
00053 int Jac1FinestMultEvent;
00054 
00055 int Jac2DiagEvent;
00056 int Jac2MultEvent;
00057 int Jac2FinestDiagEvent;
00058 int Jac2FinestMultEvent;
00059 
00060 int Jac3DiagEvent;
00061 int Jac3MultEvent;
00062 int Jac3FinestDiagEvent;
00063 int Jac3FinestMultEvent;
00064 #endif
00065 
00066 double***** LaplacianType1Stencil; 
00067 double**** LaplacianType2Stencil; 
00068 double***** MassType1Stencil; 
00069 double**** MassType2Stencil; 
00070 double****** ShapeFnStencil;
00071 
00072 double**** GradDivType2Stencil; 
00073 
00074 int main(int argc, char ** argv ) {     
00075   int size, rank;
00076   bool incCorner = 1;  
00077   double gSize[3];
00078   unsigned int ptsLen;
00079   unsigned int maxNumPts = 1;
00080   unsigned int dim = 3;
00081   unsigned int maxDepth = 30;
00082   bool compressLut = false;
00083   double localTime, globalTime;
00084   double startTime, endTime;
00085   DendroIntL localSz, totalSz;
00086   std::vector<ot::TreeNode> linOct, balOct;
00087   std::vector<double> pts;
00088   double mgLoadFac = 2.0;
00089   char partitionFileNameBase[100];
00090   bool dumpPartitions = false;
00091   char ptsFileName[256];
00092 
00093   PetscInitialize(&argc,&argv,"optionsElasticity",help);
00094   ot::RegisterEvents();
00095 
00096   ot::DAMG_Initialize(MPI_COMM_WORLD);
00097 
00098   MPI_Comm_size(MPI_COMM_WORLD,&size);
00099   MPI_Comm_rank(MPI_COMM_WORLD,&rank);
00100 
00101   if(argc < 2) {
00102     std::cout<<"  Usage: <exe> ptsFilePrefix dim[3] maxDepth[30]"<<
00103       " maxNumPtsPerOctant[1] incCorner[1] compressLut[0]"<<
00104       " mgLoadFac[2.0] partitionFileNameBase "<<std::endl; 
00105     ot::DAMG_Finalize();
00106     PetscFinalize();
00107   }
00108 
00109   if(argc > 1) {
00110     sprintf(ptsFileName, "%s%d_%d.pts", argv[1], rank, size);
00111   }  
00112   if(argc > 2) {
00113     dim = atoi(argv[2]);
00114   }
00115   if(argc > 3) {
00116     maxDepth = atoi(argv[3]);
00117   }
00118   if(argc > 4) {
00119     maxNumPts = atoi(argv[4]);
00120   }
00121   if(argc > 5) {
00122     incCorner = (bool)(atoi(argv[5]));
00123   }  
00124   if(argc > 6) { 
00125     compressLut = (bool)(atoi(argv[6]));
00126   }
00127   if(argc > 7) { 
00128     mgLoadFac = atof(argv[7]);
00129   }
00130   if(argc > 8) {
00131     dumpPartitions = true;
00132     strcpy(partitionFileNameBase, argv[8]);
00133   }
00134 
00135 #ifdef PETSC_USE_LOG
00136   PetscLogEventRegister(&Jac1DiagEvent,"ODAmatDiag",PETSC_VIEWER_COOKIE);
00137   PetscLogEventRegister(&Jac1MultEvent,"ODAmatMult",PETSC_VIEWER_COOKIE);
00138   PetscLogEventRegister(&Jac1FinestDiagEvent,"ODAmatDiagFinest",PETSC_VIEWER_COOKIE);
00139   PetscLogEventRegister(&Jac1FinestMultEvent,"ODAmatMultFinest",PETSC_VIEWER_COOKIE);
00140 
00141   PetscLogEventRegister(&Jac2DiagEvent,"OMGmatDiag-2",PETSC_VIEWER_COOKIE);
00142   PetscLogEventRegister(&Jac2MultEvent,"OMGmatMult-2",PETSC_VIEWER_COOKIE);
00143   PetscLogEventRegister(&Jac2FinestDiagEvent,"OMGmatDiagFinest-2",PETSC_VIEWER_COOKIE);
00144   PetscLogEventRegister(&Jac2FinestMultEvent,"OMGmatMultFinest-2",PETSC_VIEWER_COOKIE);
00145 
00146   PetscLogEventRegister(&Jac3DiagEvent,"OMGmatDiag-3",PETSC_VIEWER_COOKIE);
00147   PetscLogEventRegister(&Jac3MultEvent,"OMGmatMult-3",PETSC_VIEWER_COOKIE);
00148   PetscLogEventRegister(&Jac3FinestDiagEvent,"OMGmatDiagFinest-3",PETSC_VIEWER_COOKIE);
00149   PetscLogEventRegister(&Jac3FinestMultEvent,"OMGmatMultFinest-3",PETSC_VIEWER_COOKIE);
00150   PetscLogEventRegister(&vecMassDiagEvent,"vMassDiag",PETSC_VIEWER_COOKIE);
00151   PetscLogEventRegister(&vecMassMultEvent,"vMassMult",PETSC_VIEWER_COOKIE);
00152   PetscLogEventRegister(&vecMassFinestDiagEvent,"vMassDiagFinest",PETSC_VIEWER_COOKIE);
00153   PetscLogEventRegister(&vecMassFinestMultEvent,"vMassMultFinest",PETSC_VIEWER_COOKIE);
00154 
00155   PetscLogEventRegister(&elasticityDiagEvent,"elDiag",PETSC_VIEWER_COOKIE);
00156   PetscLogEventRegister(&elasticityMultEvent,"elMult",PETSC_VIEWER_COOKIE);
00157   PetscLogEventRegister(&elasticityFinestDiagEvent,"elDiagFinest",PETSC_VIEWER_COOKIE);
00158   PetscLogEventRegister(&elasticityFinestMultEvent,"elMultFinest",PETSC_VIEWER_COOKIE);
00159 
00160   int stages[3];
00161   PetscLogStageRegister(&stages[0],"P2O.");
00162   PetscLogStageRegister(&stages[1],"Bal");  
00163   PetscLogStageRegister(&stages[2],"Solve");  
00164 #endif
00165 
00166 
00167   //Read pts from files
00168   MPI_Barrier(MPI_COMM_WORLD);  
00169   if(!rank){
00170     std::cout << " reading  "<<ptsFileName<<std::endl; // Point size
00171   }
00172   ot::readPtsFromFile(ptsFileName, pts);
00173   if(!rank){
00174     std::cout << " finished reading  "<<ptsFileName<<std::endl; // Point size
00175   }
00176   MPI_Barrier(MPI_COMM_WORLD);  
00177 
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 
00195   par::removeDuplicates<ot::TreeNode>(tmpNodes,false,MPI_COMM_WORLD);   
00196   linOct = tmpNodes;
00197   tmpNodes.clear();
00198   par::partitionW<ot::TreeNode>(linOct, NULL,MPI_COMM_WORLD);
00199   // reduce and only print the total ...
00200   localSz = linOct.size();
00201   par::Mpi_Reduce<DendroIntL>(&localSz, &totalSz, 1, MPI_SUM, 0, MPI_COMM_WORLD);
00202   MPI_Barrier(MPI_COMM_WORLD);  
00203   if(rank==0) {
00204     std::cout<<"# pts= " << totalSz<<std::endl;
00205   }
00206   MPI_Barrier(MPI_COMM_WORLD);  
00207 
00208   pts.resize(3*(linOct.size()));
00209   ptsLen = (3*(linOct.size()));
00210   for(int i = 0; i < linOct.size(); i++) {
00211     pts[3*i] = (((double)(linOct[i].getX())) + 0.5)/((double)(1u << maxDepth));
00212     pts[(3*i)+1] = (((double)(linOct[i].getY())) +0.5)/((double)(1u << maxDepth));
00213     pts[(3*i)+2] = (((double)(linOct[i].getZ())) +0.5)/((double)(1u << maxDepth));
00214   }//end for i
00215   linOct.clear();
00216   gSize[0] = 1.0;
00217   gSize[1] = 1.0;
00218   gSize[2] = 1.0;
00219 
00220   //Points2Octree....
00221   MPI_Barrier(MPI_COMM_WORLD);  
00222 #ifdef PETSC_USE_LOG
00223   PetscLogStagePush(stages[0]);
00224 #endif
00225   startTime = MPI_Wtime();
00226   ot::points2Octree(pts, gSize, linOct, dim, maxDepth, maxNumPts, MPI_COMM_WORLD);
00227   endTime = MPI_Wtime();
00228 #ifdef PETSC_USE_LOG
00229   PetscLogStagePop();
00230 #endif
00231   localTime = endTime - startTime;
00232   par::Mpi_Reduce<double>(&localTime, &globalTime, 1, MPI_MAX, 0, MPI_COMM_WORLD);
00233   if(!rank){
00234     std::cout <<"P2n Time: "<<globalTime << std::endl;
00235   }
00236   // reduce and only print the total ...
00237   localSz = linOct.size();
00238   par::Mpi_Reduce<DendroIntL>(&localSz, &totalSz, 1, MPI_SUM, 0, MPI_COMM_WORLD);
00239   if(rank==0) {
00240     std::cout<<"# of Unbalanced Octants: " << totalSz<<std::endl;
00241   }
00242   pts.clear();
00243 
00244   //Balancing...
00245   MPI_Barrier(MPI_COMM_WORLD);  
00246 #ifdef PETSC_USE_LOG
00247   PetscLogStagePush(stages[1]);
00248 #endif
00249   startTime = MPI_Wtime();
00250   ot::balanceOctree (linOct, balOct, dim, maxDepth, incCorner, MPI_COMM_WORLD, NULL, NULL);
00251   endTime = MPI_Wtime();
00252 #ifdef PETSC_USE_LOG
00253   PetscLogStagePop();
00254 #endif
00255   linOct.clear();
00256   // compute total inp size and output size
00257   localSz = balOct.size();
00258   localTime = endTime - startTime;
00259   par::Mpi_Reduce<DendroIntL>(&localSz, &totalSz, 1, MPI_SUM, 0, MPI_COMM_WORLD);
00260   par::Mpi_Reduce<double>(&localTime, &globalTime, 1, MPI_MAX, 0, MPI_COMM_WORLD);
00261 
00262   if(!rank) {
00263     std::cout << "# of Balanced Octants: "<< totalSz << std::endl;       
00264     std::cout << "bal Time: "<<globalTime << std::endl;
00265   }
00266 
00267   //Solve ...
00268   MPI_Barrier(MPI_COMM_WORLD);  
00269 #ifdef PETSC_USE_LOG
00270   PetscLogStagePush(stages[2]);
00271 #endif
00272 
00273   ot::DAMG       *damg;    
00274   int       nlevels = 1; //number of multigrid levels
00275   PetscInt       numRefinements = 0;
00276   unsigned int   dof = 3; // degrees of freedom per node  
00277 
00278   iC(PetscOptionsGetInt(0,"-numRefinements",&numRefinements,0));
00279   for(int i = 0; i < numRefinements; i++) {
00280     std::vector<ot::TreeNode> tmpOct = balOct;
00281     balOct.clear();
00282     ot::refineOctree(tmpOct, balOct); 
00283   }
00284 
00285   PetscInt nlevelsPetscInt = nlevels;
00286   PetscOptionsGetInt(0, "-nlevels", &nlevelsPetscInt, 0);
00287   nlevels = nlevelsPetscInt;
00288 
00289   // Note: The user context for all levels will be set separately later.
00290   MPI_Barrier(MPI_COMM_WORLD);  
00291 
00292   if(!rank) {
00293     std::cout<<"nlevels initial: "<<nlevels<<std::endl;
00294   }
00295 
00296   ot::DAMGCreateAndSetDA(PETSC_COMM_WORLD, nlevels, NULL, &damg, 
00297       balOct, dof, mgLoadFac, compressLut, incCorner);
00298 
00299   balOct.clear();
00300 
00301   if(!rank) {
00302     std::cout<<"nlevels final: "<<nlevels<<std::endl;
00303   }
00304 
00305   MPI_Barrier(MPI_COMM_WORLD);  
00306 
00307   if(!rank) {
00308     std::cout << "Created DA for all levels."<< std::endl;
00309   }
00310 
00311   MPI_Barrier(MPI_COMM_WORLD);
00312   ot::PrintDAMG(damg);
00313   MPI_Barrier(MPI_COMM_WORLD);
00314 
00315   //Dump Partitions
00316   if(dumpPartitions) {
00317     for(int i = 0; i < nlevels; i++) {
00318       char tmpFileName[200];
00319       sprintf(tmpFileName,"%s_Lev%d.vtk", partitionFileNameBase, i);
00320       ot::writePartitionVTK(damg[i]->da, tmpFileName);
00321       if(damg[i]->da_aux) {
00322         sprintf(tmpFileName,"%s_Lev%d_Aux.vtk", partitionFileNameBase, i);
00323         ot::writePartitionVTK(damg[i]->da_aux, tmpFileName);
00324       }
00325     }
00326   }
00327 
00328   iC(DAMGDestroy(damg));
00329 
00330   if (!rank) {
00331     std::cout << GRN << "Destroyed DAMG" << NRM << std::endl;
00332   }
00333 
00334   endTime = MPI_Wtime();
00335 #ifdef PETSC_USE_LOG
00336   PetscLogStagePop();
00337 #endif
00338 
00339   if (!rank) {
00340     std::cout << GRN << "Finalizing PETSC" << NRM << std::endl;
00341   }
00342   ot::DAMG_Finalize();
00343   PetscFinalize();
00344 }//end main
00345