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
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
00070 MPI_Barrier(MPI_COMM_WORLD);
00071 if(!rank){
00072 std::cout << " reading "<<pFile<<std::endl;
00073 }
00074 ot::readPtsFromFile(pFile, pts);
00075 if(!rank){
00076 std::cout << " finished reading "<<pFile<<std::endl;
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
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 }
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
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
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
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;
00162
00163
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 }
00225