Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages

createCoarser.C

Go to the documentation of this file.
00001 
00002 #include "mpi.h"
00003 #include "petsc.h"
00004 #include "sys.h"
00005 #include "TreeNode.h"
00006 #include "parUtils.h"
00007 #include "omg.h"
00008 #include "oda.h"
00009 #include <cstdlib>
00010 #include "colors.h"
00011 #include "externVars.h"
00012 #include "dendro.h"
00013 
00014 static char help[] = "Solves 3D PBE using FEM, MG, DA and Matrix-Free methods.\n\n";
00015 //Don't want time to be synchronized. Need to check load imbalance.
00016 #ifdef MPI_WTIME_IS_GLOBAL
00017 #undef MPI_WTIME_IS_GLOBAL
00018 #endif
00019 
00020 int main(int argc, char ** argv ) {     
00021   int size, rank;
00022   unsigned int numPts;
00023   char pFile[50],oFile[50];
00024   double gSize[3];
00025   unsigned int ptsLen;
00026   unsigned int maxNumPts= 1;
00027   unsigned int dim=3;
00028   unsigned int maxDepth=30;
00029   DendroIntL locSz, totalSz;
00030   int nlevels;
00031   std::vector<ot::TreeNode> linOct, balOct;
00032   std::vector<double> pts;
00033 
00034   PetscInitialize(&argc,&argv,0,help);
00035   ot::RegisterEvents();
00036 
00037   ot::DAMG_Initialize(MPI_COMM_WORLD);
00038 
00039 #ifdef PETSC_USE_LOG
00040   int stages[3];
00041   PetscLogStageRegister(&stages[0],"P2O.");
00042   PetscLogStageRegister(&stages[1],"Bal");  
00043   PetscLogStageRegister(&stages[2],"Solve");  
00044 #endif
00045 
00046   MPI_Comm_size(MPI_COMM_WORLD,&size);
00047   MPI_Comm_rank(MPI_COMM_WORLD,&rank);
00048   if(argc < 2) {
00049     std::cerr << "Usage: " << argv[0] << "inpfile  maxDepth[30] \
00050       dim[3] maxNumPtsPerOctant[1] nlevels[2] " << std::endl;
00051     return -1;
00052   }
00053   if(argc > 2) {
00054     maxDepth = atoi(argv[2]);
00055   }
00056   if(argc > 3) {
00057     dim = atoi(argv[3]);
00058   }
00059   if(argc > 4) {
00060     maxNumPts = atoi(argv[4]);
00061   }
00062   if(argc > 5) {
00063     nlevels = atoi(argv[5]);
00064   }
00065 
00066   sprintf(pFile,"%s%d_%d.pts",argv[1],rank,size);
00067   sprintf(oFile,"%s_Out%d_%d.ot",argv[1],rank,size);
00068 
00069   //Points2Octree....
00070   MPI_Barrier(MPI_COMM_WORLD);  
00071   if(!rank){
00072     std::cout << " reading  "<<pFile<<std::endl; // Point size
00073   }
00074   ot::readPtsFromFile(pFile, pts);
00075   if(!rank){
00076     std::cout << " finished reading  "<<pFile<<std::endl; // Point size
00077   }
00078   MPI_Barrier(MPI_COMM_WORLD);  
00079   ptsLen = pts.size();
00080   std::vector<ot::TreeNode> tmpNodes;
00081   for(int i=0;i<ptsLen;i+=3) {
00082     if( (pts[i] > 0.0) &&
00083         (pts[i+1] > 0.0)  
00084         && (pts[i+2] > 0.0) &&
00085         ( ((unsigned int)(pts[i]*((double)(1u << maxDepth)))) < (1u << maxDepth))  &&
00086         ( ((unsigned int)(pts[i+1]*((double)(1u << maxDepth)))) < (1u << maxDepth))  &&
00087         ( ((unsigned int)(pts[i+2]*((double)(1u << maxDepth)))) < (1u << maxDepth)) ) {
00088       tmpNodes.push_back( ot::TreeNode((unsigned int)(pts[i]*(double)(1u << maxDepth)),
00089             (unsigned int)(pts[i+1]*(double)(1u << maxDepth)),
00090             (unsigned int)(pts[i+2]*(double)(1u << maxDepth)),
00091             maxDepth,dim,maxDepth) );
00092     }
00093   }
00094   pts.clear();
00095   par::removeDuplicates<ot::TreeNode>(tmpNodes,false,MPI_COMM_WORLD);   
00096   linOct = tmpNodes;
00097   tmpNodes.clear();
00098   par::partitionW<ot::TreeNode>(linOct, NULL,MPI_COMM_WORLD);
00099   // reduce and only print the total ...
00100   locSz = linOct.size();
00101   par::Mpi_Reduce<DendroIntL>(&locSz, &totalSz, 1, MPI_SUM, 0, MPI_COMM_WORLD);
00102   MPI_Barrier(MPI_COMM_WORLD);  
00103   if(rank==0) {
00104     std::cout<<"# pts= " << totalSz<<std::endl;
00105   }
00106   MPI_Barrier(MPI_COMM_WORLD);  
00107 
00108   pts.resize(3*(linOct.size()));
00109   ptsLen = (3*(linOct.size()));
00110   for(int i=0;i<linOct.size();i++) {
00111     pts[3*i] = (((double)(linOct[i].getX())) + 0.5)/((double)(1u << maxDepth));
00112     pts[(3*i)+1] = (((double)(linOct[i].getY())) +0.5)/((double)(1u << maxDepth));
00113     pts[(3*i)+2] = (((double)(linOct[i].getZ())) +0.5)/((double)(1u << maxDepth));
00114   }//end for i
00115   linOct.clear();
00116   gSize[0] = 1.;
00117   gSize[1] = 1.;
00118   gSize[2] = 1.;
00119 
00120   MPI_Barrier(MPI_COMM_WORLD);  
00121 #ifdef PETSC_USE_LOG
00122   PetscLogStagePush(stages[0]);
00123 #endif
00124   ot::points2Octree(pts, gSize, linOct, dim, maxDepth, maxNumPts, MPI_COMM_WORLD);
00125 #ifdef PETSC_USE_LOG
00126   PetscLogStagePop();
00127 #endif
00128   // reduce and only print the total ...
00129   locSz = linOct.size();
00130   par::Mpi_Reduce<DendroIntL>(&locSz, &totalSz, 1, MPI_SUM, 0, MPI_COMM_WORLD);
00131   if(rank==0) {
00132     std::cout<<"# of Unbalanced Octants: " << totalSz<<std::endl;
00133   }
00134   pts.clear();
00135 
00136   //Balancing...
00137   MPI_Barrier(MPI_COMM_WORLD);  
00138 #ifdef PETSC_USE_LOG
00139   PetscLogStagePush(stages[1]);
00140 #endif
00141   ot::balanceOctree (linOct, balOct, dim, maxDepth, true, MPI_COMM_WORLD, NULL, NULL);
00142 #ifdef PETSC_USE_LOG
00143   PetscLogStagePop();
00144 #endif
00145   linOct.clear();
00146   // compute total inp size and output size
00147   locSz = balOct.size();
00148   par::Mpi_Reduce<DendroIntL>(&locSz, &totalSz, 1, MPI_SUM, 0, MPI_COMM_WORLD);
00149 
00150   if(!rank) {
00151     std::cout << "# of Balanced Octants: "<< totalSz << std::endl;       
00152   }
00153 
00154   MPI_Barrier(MPI_COMM_WORLD);  
00155 #ifdef PETSC_USE_LOG
00156   PetscLogStagePush(stages[2]);
00157 #endif
00158 
00159 
00160   ot::DAMG        *damg;    
00161   unsigned int       dof =1;// degrees of freedom per node  
00162 
00163   // Note: The user context for all levels will be set separately later.
00164   MPI_Barrier(MPI_COMM_WORLD);  
00165 
00166   if(!rank) {
00167     std::cout<<"nlevels initial: "<<nlevels<<std::endl;
00168   }
00169 
00170   ot::DAMGCreateAndSetDA(PETSC_COMM_WORLD, nlevels, NULL, &damg,
00171       balOct, dof, 1.5, false, true);
00172 
00173   if(!rank) {
00174     std::cout<<"nlevels final: "<<nlevels<<std::endl;
00175   }
00176 
00177   MPI_Barrier(MPI_COMM_WORLD);  
00178 
00179   if(!rank) {
00180     std::cout << "Created DA for all levels."<< std::endl;
00181   }
00182 
00183   MPI_Barrier(MPI_COMM_WORLD);
00184   ot::PrintDAMG(damg);
00185   MPI_Barrier(MPI_COMM_WORLD);
00186 
00187   ot::DA* dac = damg[0]->da;
00188 
00189   std::vector<ot::TreeNode> outputOct;
00190 
00191   if(dac->iAmActive()) {
00192     double factor = (1.0/(static_cast<double>(1u << maxDepth)));
00193     for(dac->init<ot::DA_FLAGS::WRITABLE>();
00194         dac->curr() < dac->end<ot::DA_FLAGS::WRITABLE>();
00195         dac->next<ot::DA_FLAGS::WRITABLE>() ) {
00196       Point currPt = dac->getCurrentOffset();
00197       unsigned int lev = dac->getLevel(dac->curr());
00198       ot::TreeNode tmpNode(currPt.xint(), currPt.yint(), currPt.zint(), (lev-1), dim, maxDepth);
00199       outputOct.push_back(tmpNode); 
00200     } 
00201   }
00202 
00203   ot::writeNodesToFile(oFile, outputOct);
00204   outputOct.clear();
00205 
00206   iC(DAMGDestroy(damg));
00207 
00208   if (!rank) {
00209     std::cout << GRN << "Destroyed DAMG" << NRM << std::endl;
00210   }
00211 
00212 #ifdef PETSC_USE_LOG
00213   PetscLogStagePop();
00214 #endif
00215   balOct.clear();
00216 
00217   if (!rank) {
00218     std::cout << GRN << "Finalizing ..." << NRM << std::endl;
00219   }
00220 
00221   ot::DAMG_Finalize();
00222   PetscFinalize();
00223 
00224 }//end function
00225 

Generated on Tue Mar 24 16:13:58 2009 for DENDRO by  doxygen 1.3.9.1