DENDRO: /net/www/csela/dendro/releases/Dendro/examples/src/drivers/newElasMesh.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 double gaussian(double mean, double std_deviation);
00075 
00076 int main(int argc, char ** argv ) {     
00077   int size, rank;
00078   bool incCorner = 1;  
00079   unsigned int local_num_pts = 5000;
00080   double gSize[3];
00081   unsigned int ptsLen;
00082   unsigned int maxNumPts = 1;
00083   unsigned int dim = 3;
00084   unsigned int maxDepth = 30;
00085   bool compressLut = false;
00086   double localTime, globalTime;
00087   double startTime, endTime;
00088   DendroIntL localSz, totalSz;
00089   std::vector<ot::TreeNode> linOct, balOct;
00090   std::vector<double> pts;
00091   double mgLoadFac = 2.0;
00092   char partitionFileNameBase[100];
00093   bool dumpPartitions = false;
00094 
00095   PetscInitialize(&argc,&argv,"optionsElasticity",help);
00096   ot::RegisterEvents();
00097 
00098   ot::DAMG_Initialize(MPI_COMM_WORLD);
00099 
00100 #ifdef PETSC_USE_LOG
00101   PetscLogEventRegister(&Jac1DiagEvent,"ODAmatDiag",PETSC_VIEWER_COOKIE);
00102   PetscLogEventRegister(&Jac1MultEvent,"ODAmatMult",PETSC_VIEWER_COOKIE);
00103   PetscLogEventRegister(&Jac1FinestDiagEvent,"ODAmatDiagFinest",PETSC_VIEWER_COOKIE);
00104   PetscLogEventRegister(&Jac1FinestMultEvent,"ODAmatMultFinest",PETSC_VIEWER_COOKIE);
00105 
00106   PetscLogEventRegister(&Jac2DiagEvent,"OMGmatDiag-2",PETSC_VIEWER_COOKIE);
00107   PetscLogEventRegister(&Jac2MultEvent,"OMGmatMult-2",PETSC_VIEWER_COOKIE);
00108   PetscLogEventRegister(&Jac2FinestDiagEvent,"OMGmatDiagFinest-2",PETSC_VIEWER_COOKIE);
00109   PetscLogEventRegister(&Jac2FinestMultEvent,"OMGmatMultFinest-2",PETSC_VIEWER_COOKIE);
00110 
00111   PetscLogEventRegister(&Jac3DiagEvent,"OMGmatDiag-3",PETSC_VIEWER_COOKIE);
00112   PetscLogEventRegister(&Jac3MultEvent,"OMGmatMult-3",PETSC_VIEWER_COOKIE);
00113   PetscLogEventRegister(&Jac3FinestDiagEvent,"OMGmatDiagFinest-3",PETSC_VIEWER_COOKIE);
00114   PetscLogEventRegister(&Jac3FinestMultEvent,"OMGmatMultFinest-3",PETSC_VIEWER_COOKIE);
00115   PetscLogEventRegister(&vecMassDiagEvent,"vMassDiag",PETSC_VIEWER_COOKIE);
00116   PetscLogEventRegister(&vecMassMultEvent,"vMassMult",PETSC_VIEWER_COOKIE);
00117   PetscLogEventRegister(&vecMassFinestDiagEvent,"vMassDiagFinest",PETSC_VIEWER_COOKIE);
00118   PetscLogEventRegister(&vecMassFinestMultEvent,"vMassMultFinest",PETSC_VIEWER_COOKIE);
00119 
00120   PetscLogEventRegister(&elasticityDiagEvent,"elDiag",PETSC_VIEWER_COOKIE);
00121   PetscLogEventRegister(&elasticityMultEvent,"elMult",PETSC_VIEWER_COOKIE);
00122   PetscLogEventRegister(&elasticityFinestDiagEvent,"elDiagFinest",PETSC_VIEWER_COOKIE);
00123   PetscLogEventRegister(&elasticityFinestMultEvent,"elMultFinest",PETSC_VIEWER_COOKIE);
00124 
00125   int stages[3];
00126   PetscLogStageRegister(&stages[0],"P2O.");
00127   PetscLogStageRegister(&stages[1],"Bal");  
00128   PetscLogStageRegister(&stages[2],"Solve");  
00129 #endif
00130 
00131   MPI_Comm_size(MPI_COMM_WORLD,&size);
00132   MPI_Comm_rank(MPI_COMM_WORLD,&rank);
00133 
00134   //  Usage: <exe> local_num_pts[5000] dim[3] maxDepth[30] maxNumPtsPerOctant[1] 
00135   //  incCorner[1] compressLut[0] mgLoadFac[2.0] partitionFileNameBase  
00136 
00137   if(argc > 1) {
00138     local_num_pts = atoi(argv[1]);
00139   }  
00140   if(argc > 2) {
00141     dim = atoi(argv[2]);
00142   }
00143   if(argc > 3) {
00144     maxDepth = atoi(argv[3]);
00145   }
00146   if(argc > 4) {
00147     maxNumPts = atoi(argv[4]);
00148   }
00149   if(argc > 5) {
00150     incCorner = (bool)(atoi(argv[5]));
00151   }  
00152   if(argc > 6) { 
00153     compressLut = (bool)(atoi(argv[6]));
00154   }
00155   if(argc > 7) { 
00156     mgLoadFac = atof(argv[7]);
00157   }
00158   if(argc > 8) {
00159     dumpPartitions = true;
00160     strcpy(partitionFileNameBase, argv[8]);
00161   }
00162 
00163   //Generate Points
00164   if(!rank){
00165     std::cout << "Creating points on the fly... "<<std::endl;
00166   }
00167 
00168   MPI_Barrier(MPI_COMM_WORLD);
00169   startTime = MPI_Wtime();
00170 
00171   srand48(0x12345678 + 76543*rank);
00172 
00173   pts.resize(3*local_num_pts);
00174   for(int i = 0; i < (3*local_num_pts); i++) {
00175     pts[i]= gaussian(0.5, 0.16);
00176   }
00177   MPI_Barrier(MPI_COMM_WORLD);
00178   endTime = MPI_Wtime();
00179   localTime = endTime - startTime;
00180 
00181   if(!rank){
00182     std::cout << " finished creating points  "<<std::endl; // Point size
00183     std::cout <<"point generation time: "<<localTime << std::endl;
00184   }
00185 
00186   ptsLen = pts.size();
00187   std::vector<ot::TreeNode> tmpNodes;
00188   for(int i = 0;i < ptsLen; i+=3) {
00189     if( (pts[i] > 0.0) &&
00190         (pts[i+1] > 0.0)  
00191         && (pts[i+2] > 0.0) &&
00192         ( ((unsigned int)(pts[i]*((double)(1u << maxDepth)))) < (1u << maxDepth))  &&
00193         ( ((unsigned int)(pts[i+1]*((double)(1u << maxDepth)))) < (1u << maxDepth))  &&
00194         ( ((unsigned int)(pts[i+2]*((double)(1u << maxDepth)))) < (1u << maxDepth)) ) {
00195       tmpNodes.push_back( ot::TreeNode((unsigned int)(pts[i]*(double)(1u << maxDepth)),
00196             (unsigned int)(pts[i+1]*(double)(1u << maxDepth)),
00197             (unsigned int)(pts[i+2]*(double)(1u << maxDepth)),
00198             maxDepth,dim,maxDepth) );
00199     }
00200   }
00201   pts.clear();
00202 
00203   par::removeDuplicates<ot::TreeNode>(tmpNodes,false,MPI_COMM_WORLD);   
00204   linOct = tmpNodes;
00205   tmpNodes.clear();
00206   par::partitionW<ot::TreeNode>(linOct, NULL,MPI_COMM_WORLD);
00207   // reduce and only print the total ...
00208   localSz = linOct.size();
00209   par::Mpi_Reduce<DendroIntL>(&localSz, &totalSz, 1, MPI_SUM, 0, MPI_COMM_WORLD);
00210   MPI_Barrier(MPI_COMM_WORLD);  
00211   if(rank==0) {
00212     std::cout<<"# pts= " << totalSz<<std::endl;
00213   }
00214   MPI_Barrier(MPI_COMM_WORLD);  
00215 
00216   pts.resize(3*(linOct.size()));
00217   ptsLen = (3*(linOct.size()));
00218   for(int i = 0; i < linOct.size(); i++) {
00219     pts[3*i] = (((double)(linOct[i].getX())) + 0.5)/((double)(1u << maxDepth));
00220     pts[(3*i)+1] = (((double)(linOct[i].getY())) +0.5)/((double)(1u << maxDepth));
00221     pts[(3*i)+2] = (((double)(linOct[i].getZ())) +0.5)/((double)(1u << maxDepth));
00222   }//end for i
00223   linOct.clear();
00224   gSize[0] = 1.0;
00225   gSize[1] = 1.0;
00226   gSize[2] = 1.0;
00227 
00228   //Points2Octree....
00229   MPI_Barrier(MPI_COMM_WORLD);  
00230 #ifdef PETSC_USE_LOG
00231   PetscLogStagePush(stages[0]);
00232 #endif
00233   startTime = MPI_Wtime();
00234   ot::points2Octree(pts, gSize, linOct, dim, maxDepth, maxNumPts, MPI_COMM_WORLD);
00235   endTime = MPI_Wtime();
00236 #ifdef PETSC_USE_LOG
00237   PetscLogStagePop();
00238 #endif
00239   localTime = endTime - startTime;
00240   par::Mpi_Reduce<double>(&localTime, &globalTime, 1, MPI_MAX, 0, MPI_COMM_WORLD);
00241   if(!rank){
00242     std::cout <<"P2n Time: "<<globalTime << std::endl;
00243   }
00244   // reduce and only print the total ...
00245   localSz = linOct.size();
00246   par::Mpi_Reduce<DendroIntL>(&localSz, &totalSz, 1, MPI_SUM, 0, MPI_COMM_WORLD);
00247   if(rank==0) {
00248     std::cout<<"# of Unbalanced Octants: " << totalSz<<std::endl;
00249   }
00250   pts.clear();
00251 
00252   //Balancing...
00253   MPI_Barrier(MPI_COMM_WORLD);  
00254 #ifdef PETSC_USE_LOG
00255   PetscLogStagePush(stages[1]);
00256 #endif
00257   startTime = MPI_Wtime();
00258   ot::balanceOctree (linOct, balOct, dim, maxDepth, incCorner, MPI_COMM_WORLD, NULL, NULL);
00259   endTime = MPI_Wtime();
00260 #ifdef PETSC_USE_LOG
00261   PetscLogStagePop();
00262 #endif
00263   linOct.clear();
00264   // compute total inp size and output size
00265   localSz = balOct.size();
00266   localTime = endTime - startTime;
00267   par::Mpi_Reduce<DendroIntL>(&localSz, &totalSz, 1, MPI_SUM, 0, MPI_COMM_WORLD);
00268   par::Mpi_Reduce<double>(&localTime, &globalTime, 1, MPI_MAX, 0, MPI_COMM_WORLD);
00269 
00270   if(!rank) {
00271     std::cout << "# of Balanced Octants: "<< totalSz << std::endl;       
00272     std::cout << "bal Time: "<<globalTime << std::endl;
00273   }
00274 
00275   //Solve ...
00276   MPI_Barrier(MPI_COMM_WORLD);  
00277 #ifdef PETSC_USE_LOG
00278   PetscLogStagePush(stages[2]);
00279 #endif
00280 
00281   ot::DAMG       *damg;    
00282   int       nlevels = 1; //number of multigrid levels
00283   PetscInt       numRefinements = 0;
00284   unsigned int   dof = 3; // degrees of freedom per node  
00285 
00286   iC(PetscOptionsGetInt(0,"-numRefinements",&numRefinements,0));
00287   for(int i = 0; i < numRefinements; i++) {
00288     std::vector<ot::TreeNode> tmpOct = balOct;
00289     balOct.clear();
00290     ot::refineOctree(tmpOct, balOct); 
00291   }
00292 
00293   PetscInt nlevelsPetscInt = nlevels;
00294   PetscOptionsGetInt(0, "-nlevels", &nlevelsPetscInt, 0);
00295   nlevels = nlevelsPetscInt;
00296 
00297   // Note: The user context for all levels will be set separately later.
00298   MPI_Barrier(MPI_COMM_WORLD);  
00299 
00300   if(!rank) {
00301     std::cout<<"nlevels initial: "<<nlevels<<std::endl;
00302   }
00303 
00304   ot::DAMGCreateAndSetDA(PETSC_COMM_WORLD, nlevels, NULL, &damg, 
00305       balOct, dof, mgLoadFac, compressLut, incCorner);
00306 
00307   balOct.clear();
00308 
00309   if(!rank) {
00310     std::cout<<"nlevels final: "<<nlevels<<std::endl;
00311   }
00312 
00313   MPI_Barrier(MPI_COMM_WORLD);  
00314 
00315   if(!rank) {
00316     std::cout << "Created DA for all levels."<< std::endl;
00317   }
00318 
00319   MPI_Barrier(MPI_COMM_WORLD);
00320   ot::PrintDAMG(damg);
00321   MPI_Barrier(MPI_COMM_WORLD);
00322 
00323   //Dump Partitions
00324   if(dumpPartitions) {
00325     for(int i = 0; i < nlevels; i++) {
00326       char tmpFileName[200];
00327       sprintf(tmpFileName,"%s_Lev%d.vtk", partitionFileNameBase, i);
00328       ot::writePartitionVTK(damg[i]->da, tmpFileName);
00329       if(damg[i]->da_aux) {
00330         sprintf(tmpFileName,"%s_Lev%d_Aux.vtk", partitionFileNameBase, i);
00331         ot::writePartitionVTK(damg[i]->da_aux, tmpFileName);
00332       }
00333     }
00334   }
00335 
00336   iC(DAMGDestroy(damg));
00337 
00338   if (!rank) {
00339     std::cout << GRN << "Destroyed DAMG" << NRM << std::endl;
00340   }
00341 
00342   endTime = MPI_Wtime();
00343 #ifdef PETSC_USE_LOG
00344   PetscLogStagePop();
00345 #endif
00346 
00347   if (!rank) {
00348     std::cout << GRN << "Finalizing PETSC" << NRM << std::endl;
00349   }
00350   ot::DAMG_Finalize();
00351   PetscFinalize();
00352 }//end main
00353 
00354 double gaussian(double mean, double std_deviation) {
00355   static double t1 = 0, t2=0;
00356   double x1, x2, x3, r;
00357 
00358   using namespace std;
00359 
00360   // reuse previous calculations
00361   if(t1) {
00362     const double tmp = t1;
00363     t1 = 0;
00364     return mean + std_deviation * tmp;
00365   }
00366   if(t2) {
00367     const double tmp = t2;
00368     t2 = 0;
00369     return mean + std_deviation * tmp;
00370   }
00371 
00372   // pick randomly a point inside the unit disk
00373   do {
00374     x1 = 2 * drand48() - 1;
00375     x2 = 2 * drand48() - 1;
00376     x3 = 2 * drand48() - 1;
00377     r = x1 * x1 + x2 * x2 + x3*x3;
00378   } while(r >= 1);
00379 
00380   // Box-Muller transform
00381   r = sqrt(-2.0 * log(r) / r);
00382 
00383   // save for next call
00384   t1 = (r * x2);
00385   t2 = (r * x3);
00386 
00387   return mean + (std_deviation * r * x1);
00388 }//end gaussian
00389 
00390