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

checkError.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 "omgJac.h"
00010 #include "handleStencils.h"
00011 #include "colors.h"
00012 #include <cstdlib>
00013 #include <cstring>
00014 #include "externVars.h"
00015 #include "dendro.h"
00016 
00017 static char help[] = "Scalar Problem";
00018 
00019 //Don't want time to be synchronized. Need to check load imbalance.
00020 #ifdef MPI_WTIME_IS_GLOBAL
00021 #undef MPI_WTIME_IS_GLOBAL
00022 #endif
00023 
00024 #ifdef PETSC_USE_LOG
00025 //user-defined variables
00026 int Jac1DiagEvent;
00027 int Jac1MultEvent;
00028 int Jac1FinestDiagEvent;
00029 int Jac1FinestMultEvent;
00030 
00031 int Jac2DiagEvent;
00032 int Jac2MultEvent;
00033 int Jac2FinestDiagEvent;
00034 int Jac2FinestMultEvent;
00035 
00036 int Jac3DiagEvent;
00037 int Jac3MultEvent;
00038 int Jac3FinestDiagEvent;
00039 int Jac3FinestMultEvent;
00040 #endif
00041 
00042 double***** LaplacianType1Stencil; 
00043 double**** LaplacianType2Stencil; 
00044 double***** MassType1Stencil; 
00045 double**** MassType2Stencil; 
00046 double****** ShapeFnStencil;
00047 
00048 int main(int argc, char ** argv ) {     
00049   int size, rank;
00050   bool incCorner = 1;  
00051   unsigned int numPts;
00052   unsigned int solveU = 0;
00053   unsigned int writeB = 0;  
00054   char Kstr[20];
00055   char pFile[50],bFile[50],uFile[50];
00056   double gSize[3];
00057   unsigned int ptsLen;
00058   unsigned int maxNumPts= 1;
00059   unsigned int dim=3;
00060   unsigned int maxDepth=30;
00061   bool compressLut=true;
00062   double localTime, totalTime;
00063   double startTime, endTime;
00064   DendroIntL locSz, totalSz;
00065   std::vector<ot::TreeNode> linOct, balOct;
00066   std::vector<double> pts;
00067   double mgLoadFac = 2.0;
00068 
00069   PetscInitialize(&argc,&argv,"options",help);
00070   ot::RegisterEvents();
00071 
00072   ot::DAMG_Initialize(MPI_COMM_WORLD);
00073 
00074 #ifdef PETSC_USE_LOG
00075   PetscLogEventRegister(&Jac1DiagEvent,"ODAmatDiag",PETSC_VIEWER_COOKIE);
00076   PetscLogEventRegister(&Jac1MultEvent,"ODAmatMult",PETSC_VIEWER_COOKIE);
00077   PetscLogEventRegister(&Jac1FinestDiagEvent,"ODAmatDiagFinest",PETSC_VIEWER_COOKIE);
00078   PetscLogEventRegister(&Jac1FinestMultEvent,"ODAmatMultFinest",PETSC_VIEWER_COOKIE);
00079 
00080   PetscLogEventRegister(&Jac2DiagEvent,"OMGmatDiag-2",PETSC_VIEWER_COOKIE);
00081   PetscLogEventRegister(&Jac2MultEvent,"OMGmatMult-2",PETSC_VIEWER_COOKIE);
00082   PetscLogEventRegister(&Jac2FinestDiagEvent,"OMGmatDiagFinest-2",PETSC_VIEWER_COOKIE);
00083   PetscLogEventRegister(&Jac2FinestMultEvent,"OMGmatMultFinest-2",PETSC_VIEWER_COOKIE);
00084 
00085   PetscLogEventRegister(&Jac3DiagEvent,"OMGmatDiag-3",PETSC_VIEWER_COOKIE);
00086   PetscLogEventRegister(&Jac3MultEvent,"OMGmatMult-3",PETSC_VIEWER_COOKIE);
00087   PetscLogEventRegister(&Jac3FinestDiagEvent,"OMGmatDiagFinest-3",PETSC_VIEWER_COOKIE);
00088   PetscLogEventRegister(&Jac3FinestMultEvent,"OMGmatMultFinest-3",PETSC_VIEWER_COOKIE);
00089 
00090   int stages[3];
00091   PetscLogStageRegister(&stages[0],"P2O.");
00092   PetscLogStageRegister(&stages[1],"Bal");  
00093   PetscLogStageRegister(&stages[2],"Solve");  
00094 #endif
00095 
00096   MPI_Comm_size(MPI_COMM_WORLD,&size);
00097   MPI_Comm_rank(MPI_COMM_WORLD,&rank);
00098   if(argc < 3) {
00099     std::cerr << "Usage: " << argv[0] << "inpfile  maxDepth[30] solveU[0]\
00100       writeB[0] dim[3] maxNumPtsPerOctant[1] incCorner[1] compressLut[1] mgLoadFac[2.0] " << std::endl;
00101     return -1;
00102   }
00103   if(argc > 2) {
00104     maxDepth = atoi(argv[2]);
00105   }
00106   if(argc > 3) {
00107     solveU = atoi(argv[3]);
00108   }
00109   if(argc > 4) {
00110     writeB = atoi(argv[4]);
00111   }
00112   if(argc > 5) {
00113     dim = atoi(argv[5]);
00114   }
00115   if(argc > 6) {
00116     maxNumPts = atoi(argv[6]);
00117   }
00118   if(argc > 7) { incCorner = (bool)(atoi(argv[7]));}  
00119   if(argc > 8) { compressLut = (bool)(atoi(argv[8]));}
00120   if(argc > 9) { mgLoadFac = atof(argv[9]); }
00121 
00122   strcpy(bFile,argv[1]);
00123   ot::int2str(rank,Kstr);
00124   strcat(bFile,Kstr);
00125   strcat(bFile,"_\0");
00126   ot::int2str(size,Kstr);
00127   strcat(bFile,Kstr);
00128   strcpy(pFile,bFile);
00129   strcpy(uFile,bFile);
00130   strcat(bFile,"_Bal.ot\0");
00131   strcat(pFile,".pts\0");
00132   strcat(uFile,".sol\0");
00133 
00134   //Points2Octree....
00135   MPI_Barrier(MPI_COMM_WORLD);  
00136   if(!rank){
00137     std::cout << " reading  "<<pFile<<std::endl; // Point size
00138   }
00139   ot::readPtsFromFile(pFile, pts);
00140   if(!rank){
00141     std::cout << " finished reading  "<<pFile<<std::endl; // Point size
00142   }
00143   MPI_Barrier(MPI_COMM_WORLD);  
00144   ptsLen = pts.size();
00145   std::vector<ot::TreeNode> tmpNodes;
00146   for(int i=0;i<ptsLen;i+=3) {
00147     if( (pts[i] > 0.0) &&
00148         (pts[i+1] > 0.0)  
00149         && (pts[i+2] > 0.0) &&
00150         ( ((unsigned int)(pts[i]*((double)(1u << maxDepth)))) < (1u << maxDepth))  &&
00151         ( ((unsigned int)(pts[i+1]*((double)(1u << maxDepth)))) < (1u << maxDepth))  &&
00152         ( ((unsigned int)(pts[i+2]*((double)(1u << maxDepth)))) < (1u << maxDepth)) ) {
00153       tmpNodes.push_back( ot::TreeNode((unsigned int)(pts[i]*(double)(1u << maxDepth)),
00154             (unsigned int)(pts[i+1]*(double)(1u << maxDepth)),
00155             (unsigned int)(pts[i+2]*(double)(1u << maxDepth)),
00156             maxDepth,dim,maxDepth) );
00157     }
00158   }
00159   pts.clear();
00160   par::removeDuplicates<ot::TreeNode>(tmpNodes,false,MPI_COMM_WORLD);   
00161   linOct = tmpNodes;
00162   tmpNodes.clear();
00163   par::partitionW<ot::TreeNode>(linOct, NULL,MPI_COMM_WORLD);
00164   // reduce and only print the total ...
00165   locSz = linOct.size();
00166   DendroIntL totalNumPts;
00167   par::Mpi_Reduce<DendroIntL>(&locSz, &totalNumPts, 1, MPI_SUM, 0, MPI_COMM_WORLD);
00168   MPI_Barrier(MPI_COMM_WORLD);  
00169   if(rank==0) {
00170     std::cout<<"# pts= " << totalNumPts<<std::endl;
00171   }
00172   MPI_Barrier(MPI_COMM_WORLD);  
00173 
00174   PetscTruth setMatPropsUsingPts;
00175   PetscOptionsHasName(0,"-setMatPropsUsingPts",&setMatPropsUsingPts);
00176 
00177   std::vector<ot::TreeNode> matPropNodes;
00178   if(setMatPropsUsingPts) {
00179     matPropNodes = linOct;
00180   }
00181 
00182   pts.resize(3*(linOct.size()));
00183   ptsLen = (3*(linOct.size()));
00184 
00185   for(int i=0;i<linOct.size();i++) {
00186     pts[3*i] = (((double)(linOct[i].getX())) + 0.5)/((double)(1u << maxDepth));
00187     pts[(3*i)+1] = (((double)(linOct[i].getY())) +0.5)/((double)(1u << maxDepth));
00188     pts[(3*i)+2] = (((double)(linOct[i].getZ())) +0.5)/((double)(1u << maxDepth));
00189   }//end for i
00190   linOct.clear();
00191   gSize[0] = 1.;
00192   gSize[1] = 1.;
00193   gSize[2] = 1.;
00194 
00195   PetscTruth usingRegularOctree;
00196   PetscInt regLev;
00197   PetscOptionsHasName(0,"-useRegularOctreeAtLevel",&usingRegularOctree);
00198   PetscOptionsGetInt(0,"-useRegularOctreeAtLevel",&regLev,0);
00199 
00200   if(usingRegularOctree) {
00201     if(!rank) {
00202       std::cout<<"Creating Regular Fine Grid Octree."<<std::endl;
00203     }
00204     balOct.clear();
00205     createRegularOctree(balOct,regLev,3,maxDepth,MPI_COMM_WORLD);
00206   } else {   
00207     MPI_Barrier(MPI_COMM_WORLD);        
00208 #ifdef PETSC_USE_LOG
00209     PetscLogStagePush(stages[0]);
00210 #endif
00211     startTime = MPI_Wtime();
00212     ot::points2Octree(pts, gSize, linOct, dim, maxDepth, maxNumPts, MPI_COMM_WORLD);
00213     endTime = MPI_Wtime();
00214 #ifdef PETSC_USE_LOG
00215     PetscLogStagePop();
00216 #endif
00217     localTime = endTime - startTime;
00218     par::Mpi_Reduce<double>(&localTime, &totalTime, 1, MPI_MAX, 0, MPI_COMM_WORLD);
00219     if(!rank){
00220       std::cout <<"P2n Time: "<<totalTime << std::endl;
00221     }
00222     // reduce and only print the total ...
00223     locSz = linOct.size();
00224     par::Mpi_Reduce<DendroIntL>(&locSz, &totalSz, 1, MPI_SUM, 0, MPI_COMM_WORLD);
00225     if(rank==0) {
00226       std::cout<<"# of Unbalanced Octants: " << totalSz<<std::endl;
00227     }
00228     pts.clear();
00229 
00230     //Balancing...
00231     MPI_Barrier(MPI_COMM_WORLD);        
00232 #ifdef PETSC_USE_LOG
00233     PetscLogStagePush(stages[1]);
00234 #endif
00235     startTime = MPI_Wtime();
00236     ot::balanceOctree (linOct, balOct, dim, maxDepth, incCorner, MPI_COMM_WORLD, NULL, NULL);
00237     endTime = MPI_Wtime();
00238 #ifdef PETSC_USE_LOG
00239     PetscLogStagePop();
00240 #endif
00241     linOct.clear();
00242   }
00243 
00244   PetscInt minCoarsestLevel = 0;
00245   PetscOptionsGetInt(0,"-minCoarsestLevel",&minCoarsestLevel,0);
00246 
00247   assert(minCoarsestLevel <= maxDepth);
00248 
00249   if(!rank) {
00250     std::cout<<"min coarsest level reset to: "<<minCoarsestLevel<<std::endl;
00251   }
00252 
00253   //Since we only refine the coarse octants in an already balanced octree, the
00254   //2:1 balance constraint will be preserved after the following operation.
00255   //The following loop preserves the sorted order.
00256   //Linearity is preserved as well.
00257   std::vector<ot::TreeNode> tmpBalOct;
00258   for(unsigned int i = 0; i < balOct.size(); i++) {
00259     if(balOct[i].getLevel() >= minCoarsestLevel) {
00260       tmpBalOct.push_back(balOct[i]);
00261     } else {     
00262       ot::TreeNode currOct(balOct[i].getX(), balOct[i].getY(),
00263           balOct[i].getZ(), (minCoarsestLevel-1), dim, maxDepth);
00264       while(true) {
00265         currOct.addChildren(tmpBalOct);
00266         unsigned char cNum = currOct.getChildNumber();       
00267         while((cNum == 7) && (balOct[i].isAncestor(currOct))) {
00268           ot::TreeNode tmpNode = currOct.getParent();
00269           currOct = tmpNode;
00270           cNum = currOct.getChildNumber();
00271         }
00272         //At this point, we need the first decendant (at level =
00273         //minCoarsestLevel-1) of the next brother of currOct 
00274         if(balOct[i].isAncestor(currOct)) {
00275           ot::TreeNode parOfCurr = currOct.getParent();
00276           std::vector<ot::TreeNode> tmpChildren;
00277           parOfCurr.addChildren(tmpChildren);
00278           currOct = ot::TreeNode(tmpChildren[cNum+1].getX(), tmpChildren[cNum+1].getY(),
00279               tmpChildren[cNum+1].getZ(), (minCoarsestLevel-1), dim, maxDepth);
00280         } else {
00281           break;
00282         }
00283       }
00284     }
00285   }//end for i
00286   balOct = tmpBalOct;
00287   tmpBalOct.clear();
00288 
00289   if(writeB) { 
00290     ot::writeNodesToFile(bFile,balOct);
00291   }
00292   // compute total inp size and output size
00293   locSz = balOct.size();
00294   localTime = endTime - startTime;
00295   par::Mpi_Reduce<DendroIntL>(&locSz, &totalSz, 1, MPI_SUM, 0, MPI_COMM_WORLD);
00296   par::Mpi_Reduce<double>(&localTime, &totalTime, 1, MPI_MAX, 0, MPI_COMM_WORLD);
00297 
00298   if(!rank) {
00299     std::cout << "# of Balanced Octants: "<< totalSz << std::endl;       
00300     std::cout << "bal Time: "<<totalTime << std::endl;
00301   }
00302 
00303   //Solve ...
00304   MPI_Barrier(MPI_COMM_WORLD);  
00305 #ifdef PETSC_USE_LOG
00306   PetscLogStagePush(stages[2]);
00307 #endif
00308 
00309   ot::DAMG        *damg;    
00310   int       nlevels = 1; //number of multigrid levels
00311   PetscInt       numRefinements = 0;
00312   unsigned int       dof =1;// degrees of freedom per node  
00313 
00314   PetscInt nlevelsPetscInt = nlevels;
00315   PetscOptionsGetInt(0, "-nlevels", &nlevelsPetscInt,0);
00316   nlevels = nlevelsPetscInt;
00317 
00318   PetscOptionsGetInt(0,"-numRefinements",&numRefinements,0);
00319   for(int i = 0; i < numRefinements; i++) {
00320     std::vector<ot::TreeNode> tmpOct = balOct;
00321     balOct.clear();
00322     ot::refineOctree(tmpOct, balOct); 
00323   }
00324 
00325   MPI_Barrier(MPI_COMM_WORLD);  
00326 
00327   if(!rank) {
00328     std::cout<<"nlevels initial: "<<nlevels<<std::endl;
00329   }
00330 
00331   // Note: The user context for all levels will be set separately later.
00332   ot::DAMGCreateAndSetDA(PETSC_COMM_WORLD, nlevels, NULL, &damg,
00333       balOct, dof, mgLoadFac, compressLut, incCorner);
00334 
00335   if(!rank) {
00336     std::cout<<"nlevels final: "<<nlevels<<std::endl;
00337   }
00338 
00339   MPI_Barrier(MPI_COMM_WORLD);  
00340 
00341   if(!rank) {
00342     std::cout << "Created DA for all levels."<< std::endl;
00343   }
00344 
00345   MPI_Barrier(MPI_COMM_WORLD);
00346   ot::PrintDAMG(damg);
00347 
00348   unsigned int* allSizes = NULL;
00349 
00350   if(!rank) {
00351     allSizes = new unsigned int[6*size]; 
00352   }
00353 
00354   for(int lev = 0; lev < nlevels; lev++) {
00355     for(int auxCtr = 0; auxCtr < 2; auxCtr++) {
00356       ot::DA* currDa = NULL;
00357 
00358       if(auxCtr == 0) {
00359         currDa = damg[lev]->da;
00360       }else if(damg[lev]->da_aux) {
00361         currDa = damg[lev]->da_aux;
00362       } else {
00363         break;
00364       }
00365 
00366       unsigned int localSizes[6];
00367 
00368       //Pre-ghost size
00369       localSizes[0] = currDa->getIdxElementBegin();
00370 
00371       //Own element size
00372       localSizes[2] = currDa->getElementSize();
00373 
00374       unsigned int localTotalSize = currDa->getLocalBufferSize();
00375 
00376       //Post-ghost size
00377       localSizes[1] = (localTotalSize - (currDa->getIdxPostGhostBegin()));
00378 
00379       //Independent size
00380       localSizes[3] = currDa->getIndependentSize();
00381 
00382       //PreGhostElementSize
00383       localSizes[4] = currDa->getPreGhostElementSize();
00384 
00385       //Own Boundary size
00386       localSizes[5] = ((currDa->getIdxPostGhostBegin()) - (currDa->getIdxElementEnd()));
00387 
00388       MPI_Gather(localSizes, 6, MPI_UNSIGNED, allSizes, 6,
00389           MPI_UNSIGNED, 0, MPI_COMM_WORLD);
00390 
00391       MPI_Barrier(MPI_COMM_WORLD);
00392 
00393       if(!rank) {
00394         std::cout<<"Printing Info for lev: "<<lev<<" auxCtr: "<<auxCtr<<std::endl; 
00395         for(int i = 0; i < size; i++) {
00396           std::cout<<"Proc: "<<i<<" PreGh: "<<allSizes[6*i]
00397             <<" PostGh: "<<allSizes[(6*i)+1]<<" own: "
00398             <<allSizes[(6*i)+2]<<" Ind: "<<allSizes[(6*i)+3]
00399             <<" preGhElems: "<<allSizes[(6*i)+4]<<" LocalBnd: "
00400             <<allSizes[(6*i)+5]<<std::endl;
00401           fflush(stdout);
00402         }
00403       }
00404       fflush(stdout);
00405     }//end da or da_aux
00406   }//end lev
00407 
00408   if(!rank) {
00409     delete [] allSizes;
00410   }
00411 
00412   MPI_Barrier(MPI_COMM_WORLD);
00413 
00414   if(setMatPropsUsingPts) {
00415     PetscTruth setMatPropsAtCoarsest;
00416     PetscOptionsHasName(0,"-setMatPropsAtCoarsest",&setMatPropsAtCoarsest);
00417     if(setMatPropsAtCoarsest) {
00418       //Coarsest is only 1 processor. So to align with coarse blocks, everything
00419       //must be sent to p0
00420       if(!rank) {
00421         par::scatterValues<ot::TreeNode>(matPropNodes, linOct,
00422             totalNumPts, MPI_COMM_WORLD);
00423       } else {
00424         par::scatterValues<ot::TreeNode>(matPropNodes, linOct,
00425             0, MPI_COMM_WORLD);
00426       }
00427       matPropNodes.clear();
00428 
00429       std::vector<double> matPropPts(3*(linOct.size()));
00430       for(int i=0;i<linOct.size();i++) {
00431         matPropPts[3*i] =
00432           (((double)(linOct[i].getX())) + 0.5)/((double)(1u << maxDepth));
00433         matPropPts[(3*i)+1] =
00434           (((double)(linOct[i].getY())) +0.5)/((double)(1u << maxDepth));
00435         matPropPts[(3*i)+2] =
00436           (((double)(linOct[i].getZ())) +0.5)/((double)(1u << maxDepth));
00437       }//end for i
00438       linOct.clear();
00439 
00440       PetscReal lapFac = 0.0;
00441       PetscTruth optFound;
00442       PetscOptionsGetReal("lap","-MatPropFac",&lapFac,&optFound);
00443       std::vector<double> lapJumps((matPropPts.size()/3));
00444       for(int i=0;i<lapJumps.size();i++) {
00445         lapJumps[i] = lapFac;
00446       }
00447       SetCoarseToFineFromPts(damg, matPropPts, lapJumps);
00448     } else {
00449       if(damg[nlevels-1]->da->iAmActive()) {
00450         int npesActive;
00451         int rankActive;
00452 
00453         MPI_Comm_size(damg[nlevels-1]->da->getCommActive(), &npesActive);        
00454         MPI_Comm_rank(damg[nlevels-1]->da->getCommActive(), &rankActive);        
00455 
00456         unsigned int avgSize = (totalNumPts/npesActive);
00457         unsigned int remSize = (totalNumPts % npesActive);
00458 
00459         if(rankActive < remSize) {
00460           par::scatterValues<ot::TreeNode>(matPropNodes, linOct,
00461               (avgSize + 1), MPI_COMM_WORLD);
00462         } else {
00463           par::scatterValues<ot::TreeNode>(matPropNodes, linOct,
00464               avgSize, MPI_COMM_WORLD);
00465         }
00466 
00467         matPropNodes.clear();
00468         for(unsigned int i = 0; i < linOct.size(); i++) {
00469           unsigned int xpt, ypt, zpt;
00470           linOct[i].getAnchor(xpt,ypt,zpt);
00471           ot::TreeNode tmpNode(xpt,ypt,zpt,(maxDepth+1),3,(maxDepth+1));
00472           matPropNodes.push_back(tmpNode); 
00473         }
00474         linOct.clear();
00475 
00476         std::vector<ot::TreeNode> minsArray(npesActive);
00477 
00478         std::vector<ot::TreeNode> blocks = damg[nlevels-1]->da->getBlocks();
00479 
00480         assert(!blocks.empty());
00481 
00482         ot::TreeNode firstBlock = blocks[0];
00483 
00484         par::Mpi_Allgather<ot::TreeNode>(&firstBlock, &(*(minsArray.begin())), 1, 
00485             damg[nlevels-1]->da->getCommActive());
00486 
00487         int *sendCnt = new int[npesActive];
00488         int *recvCnt = new int[npesActive];
00489         int *sendOffsets = new int[npesActive];
00490         int *recvOffsets = new int[npesActive];
00491 
00492         //compute the total number of nodes being sent to each proc ...
00493         for (int i = 0; i < npesActive; i++) {
00494           sendCnt[i]=0;
00495           recvCnt[i]=0;
00496         }
00497 
00498         unsigned int blockCtr = 0;
00499         for(unsigned int tempCtr = 0; tempCtr < matPropNodes.size(); tempCtr++) {
00500           while( ((blockCtr + 1) < npesActive) &&
00501               (matPropNodes[tempCtr] >= minsArray[blockCtr+1]) ) {
00502             blockCtr++;
00503           }
00504           if ( (blockCtr + 1) == npesActive) {
00505             sendCnt[blockCtr] += (matPropNodes.size() - tempCtr);
00506             break;
00507           } else {
00508             sendCnt[blockCtr]++;
00509           }
00510         }
00511 
00512         minsArray.clear();
00513 
00514         // communicate with other procs how many you shall be sending and get how
00515         // many to recieve from whom.
00516         par::Mpi_Alltoall<int>(sendCnt, recvCnt, 1, damg[nlevels-1]->da->getCommActive());
00517 
00518         unsigned int totalRecv = 0;
00519         for (unsigned int i = 0; i < npesActive; i++) {
00520           totalRecv += recvCnt[i];
00521         }//end for i
00522 
00523         linOct.resize(totalRecv);
00524 
00525         sendOffsets[0] = 0;
00526         recvOffsets[0] = 0;
00527 
00528         // compute offsets ...
00529         for (int i=1; i < npesActive; i++) {
00530           sendOffsets[i] = sendOffsets[i-1] + sendCnt[i-1];
00531           recvOffsets[i] = recvOffsets[i-1] + recvCnt[i-1];
00532         }//end for i
00533 
00534         // perform All2Allv
00535         par::Mpi_Alltoallv_sparse<ot::TreeNode>(&(*(matPropNodes.begin())), sendCnt, sendOffsets,        
00536             &(*(linOct.begin())), recvCnt, recvOffsets, damg[nlevels-1]->da->getCommActive());
00537 
00538         matPropNodes.clear();
00539 
00540         // clean up ...
00541         delete [] sendCnt;
00542         sendCnt = NULL;
00543 
00544         delete [] recvCnt;
00545         recvCnt = NULL;
00546 
00547         delete [] sendOffsets;
00548         sendOffsets = NULL;
00549 
00550         delete [] recvOffsets;
00551         recvOffsets = NULL;
00552 
00553       } else {
00554         par::scatterValues<ot::TreeNode>(matPropNodes, linOct,
00555             0, MPI_COMM_WORLD);
00556         matPropNodes.clear();
00557       }
00558 
00559       std::vector<double> matPropPts(3*(linOct.size()));
00560       for(int i=0;i<linOct.size();i++) {
00561         matPropPts[3*i] =
00562           (((double)(linOct[i].getX())) + 0.5)/((double)(1u << maxDepth));
00563         matPropPts[(3*i)+1] =
00564           (((double)(linOct[i].getY())) +0.5)/((double)(1u << maxDepth));
00565         matPropPts[(3*i)+2] =
00566           (((double)(linOct[i].getZ())) +0.5)/((double)(1u << maxDepth));
00567       }//end for i
00568       linOct.clear();
00569 
00570       PetscReal lapFac = 0.0;
00571       PetscTruth optFound;
00572       PetscOptionsGetReal("lap","-MatPropFac",&lapFac,&optFound);
00573       std::vector<double> lapJumps((matPropPts.size()/3));
00574       for(int i=0;i<lapJumps.size();i++) {
00575         lapJumps[i] = lapFac;
00576       }
00577       SetUserContextsFromPts(damg, matPropPts, lapJumps); 
00578     }
00579   } else {
00580     SetUserContexts(damg);
00581   }
00582 
00583   MPI_Barrier(MPI_COMM_WORLD);
00584   if(!rank) {
00585     std::cout << "Set User Contexts all levels."<< std::endl;
00586   }
00587 
00588   MPI_Barrier(MPI_COMM_WORLD);
00589 
00590   PetscInt       jacType = 1;
00591   PetscOptionsGetInt(0,"-jacType",&jacType,0);
00592 
00593   MPI_Barrier(MPI_COMM_WORLD);
00594   if(!rank) {
00595     std::cout << "Creating stencils..."<< std::endl;
00596   }
00597   MPI_Barrier(MPI_COMM_WORLD);
00598 
00599   createLmatType2(LaplacianType2Stencil);
00600 
00601   createMmatType2(MassType2Stencil);
00602 
00603   if(jacType == 3) {
00604     createLmatType1(LaplacianType1Stencil);
00605     createMmatType1(MassType1Stencil);
00606   }
00607   createShFnMat(ShapeFnStencil);
00608 
00609   MPI_Barrier(MPI_COMM_WORLD);
00610   if(!rank) {
00611     std::cout << "Created all stencils."<< std::endl;
00612   }
00613   MPI_Barrier(MPI_COMM_WORLD);
00614 
00615   PetscInt rhsType = 2;
00616   PetscOptionsGetInt(0,"-rhsType",&rhsType,0);
00617 
00618   //Function handles
00619   PetscErrorCode (*ComputeRHSHandle)(ot::DAMG damg,Vec rhs) = NULL;
00620   PetscErrorCode (*CreateJacobianHandle)(ot::DAMG damg,Mat *B) = NULL;
00621   PetscErrorCode (*ComputeJacobianHandle)(ot::DAMG damg,Mat J, Mat B) = NULL;
00622 
00623   if(rhsType == 0) {
00624     ComputeRHSHandle = ComputeRHS0;
00625   } else if (rhsType == 1) {
00626     ComputeRHSHandle = ComputeRHS1;
00627   } else if (rhsType == 2) {
00628     ComputeRHSHandle = ComputeRHS2;
00629   } else if (rhsType == 3) {
00630     ComputeRHSHandle = ComputeRHS3;
00631   } else if (rhsType == 4) {
00632     ComputeRHSHandle = ComputeRHS4;
00633   } else if (rhsType == 5) {
00634     ComputeRHSHandle = ComputeRHS5;
00635   } else if (rhsType == 6) {
00636     ComputeRHSHandle = ComputeRHS6;
00637   } else if (rhsType == 7) {
00638     ComputeRHSHandle = ComputeRHS7;
00639   } else if (rhsType == 8) {
00640     ComputeRHSHandle = ComputeRHS8;
00641   } else {
00642     assert(false);
00643   }
00644 
00645   if(jacType == 1) {
00646     CreateJacobianHandle = CreateJacobian1;
00647     ComputeJacobianHandle = ComputeJacobian1;
00648   } else if (jacType == 2) {
00649     CreateJacobianHandle = CreateJacobian2;
00650     ComputeJacobianHandle = ComputeJacobian2;
00651   } else if (jacType == 3) {
00652     CreateJacobianHandle = CreateJacobian3;
00653     ComputeJacobianHandle = ComputeJacobian3;
00654     //Skip the finest and the coarsest levels. For the other levels, J and B
00655     //must be different
00656     for(int i = 1; i < (nlevels-1); i++) {
00657       ot::DAMGCreateJMatrix(damg[i], CreateJacobianHandle);
00658     }
00659   } else {
00660     assert(false);
00661   }
00662 
00663   //Global Function Handles for using KSP_Shell (will be used @ the coarsest grid if not all
00664   //processors are active on the coarsest grid)
00665   if (jacType == 1) {
00666     ot::getPrivateMatricesForKSP_Shell = getPrivateMatricesForKSP_Shell_Jac1;
00667   } else if (jacType == 2) {
00668     ot::getPrivateMatricesForKSP_Shell = getPrivateMatricesForKSP_Shell_Jac2;
00669   } else if (jacType == 3) {
00670     ot::getPrivateMatricesForKSP_Shell = getPrivateMatricesForKSP_Shell_Jac3;
00671   } else {
00672     assert(false);
00673   }
00674 
00675   ot::DAMGSetKSP(damg, CreateJacobianHandle, ComputeJacobianHandle, ComputeRHSHandle);
00676 
00677   MPI_Barrier(MPI_COMM_WORLD);
00678   if(!rank) {
00679     std::cout << "Called DAMGSetKSP"<< std::endl;
00680   }
00681   MPI_Barrier(MPI_COMM_WORLD);
00682 
00683   PetscReal norm2;
00684   PetscReal normInf;
00685 
00686   PetscTruth setRandomGuess = PETSC_FALSE;
00687   PetscOptionsHasName(0,"-setRandomGuess",&setRandomGuess);
00688 
00689   if(setRandomGuess) { 
00690     ot::DAMGSetInitialGuess(damg,ot::DAMGInitialGuessCurrent);
00691 
00692     PetscRandom rctx;  
00693     PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
00694     PetscRandomSetType(rctx,PETSCRAND48);
00695     PetscInt randomSeed = 12345;
00696     PetscOptionsGetInt(0,"-randomSeed",&randomSeed,0);
00697     if(!rank) {
00698       std::cout<<"Using Random Seed: "<<randomSeed<<std::endl;
00699     }
00700     PetscRandomSetSeed(rctx,randomSeed);
00701     PetscRandomSeed(rctx);
00702     PetscRandomSetFromOptions(rctx);
00703 
00704     VecSetRandom((DAMGGetx(damg)),rctx);
00705 
00706     PetscRandomDestroy(rctx);
00707 
00708     VecNorm((DAMGGetx(damg)),NORM_INFINITY,&normInf);
00709     VecNorm((DAMGGetx(damg)),NORM_2,&norm2);
00710 
00711     MPI_Barrier(MPI_COMM_WORLD);
00712     if(!rank) {
00713       std::cout << "Solving with random intial guess"<< std::endl;
00714       std::cout<<"Initial Norm2: "<<norm2<<" NormInf: "<<normInf<<std::endl;
00715     }
00716     MPI_Barrier(MPI_COMM_WORLD);
00717   }else {
00718     MPI_Barrier(MPI_COMM_WORLD);
00719     if(!rank) {
00720       std::cout << "Solving with 0-intial guess"<< std::endl;
00721     }
00722     MPI_Barrier(MPI_COMM_WORLD);
00723   }
00724 
00725   startTime = MPI_Wtime();
00726   ot::DAMGSolve(damg);
00727   endTime = MPI_Wtime();
00728 
00729   localTime = (endTime - startTime);
00730   par::Mpi_Reduce<double>(&localTime, &totalTime, 1, MPI_MAX, 0, MPI_COMM_WORLD);
00731 
00732   MPI_Barrier(MPI_COMM_WORLD);
00733 
00734   if (!rank) {
00735     std::cout << GRN << "Done Solve" << NRM << std::endl;
00736     std::cout << "Solve Time: "<<totalTime << std::endl;
00737   }
00738 
00739   MPI_Barrier(MPI_COMM_WORLD);
00740 
00741   double error = ComputeError5(DAMGGetDA(damg), DAMGGetx(damg));
00742 
00743   if(!rank) {
00744     std::cout<<"L2 Error: "<<error<<std::endl;
00745   }
00746 
00747   double testError = TestError5(DAMGGetDA(damg));
00748 
00749   if(!rank) {
00750     std::cout<<"Test Error: "<<testError<<std::endl;
00751   }
00752 
00753   destroyLmatType2(LaplacianType2Stencil);
00754   destroyMmatType2(MassType2Stencil);
00755   if(jacType == 3) {
00756     destroyLmatType1(LaplacianType1Stencil);
00757     destroyMmatType1(MassType1Stencil);
00758   }
00759   destroyShFnMat(ShapeFnStencil);
00760 
00761   if (!rank) {
00762     std::cout << GRN << "Destroyed Stencils" << NRM << std::endl;
00763   }
00764 
00765   DestroyUserContexts(damg);
00766 
00767   if (!rank) {
00768     std::cout << GRN << "Destroyed User Contexts." << NRM << std::endl;
00769   }
00770 
00771   MPI_Barrier(MPI_COMM_WORLD);
00772 
00773   DAMGDestroy(damg);
00774 
00775   if (!rank) {
00776     std::cout << GRN << "Destroyed DAMG" << NRM << std::endl;
00777   }
00778   MPI_Barrier(MPI_COMM_WORLD);
00779 
00780 #ifdef PETSC_USE_LOG
00781   PetscLogStagePop();
00782 #endif
00783   balOct.clear();
00784 
00785   ot::DAMG_Finalize();
00786   if (!rank) {
00787     std::cout << GRN << "Finalizing ..." << NRM << std::endl;
00788   }
00789   MPI_Barrier(MPI_COMM_WORLD);
00790   PetscFinalize();
00791 }//end function
00792 

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