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
00020 #ifdef MPI_WTIME_IS_GLOBAL
00021 #undef MPI_WTIME_IS_GLOBAL
00022 #endif
00023
00024 #ifdef PETSC_USE_LOG
00025
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
00135 MPI_Barrier(MPI_COMM_WORLD);
00136 if(!rank){
00137 std::cout << " reading "<<pFile<<std::endl;
00138 }
00139 ot::readPtsFromFile(pFile, pts);
00140 if(!rank){
00141 std::cout << " finished reading "<<pFile<<std::endl;
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
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 }
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",®Lev,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
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
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
00254
00255
00256
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
00273
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 }
00286 balOct = tmpBalOct;
00287 tmpBalOct.clear();
00288
00289 if(writeB) {
00290 ot::writeNodesToFile(bFile,balOct);
00291 }
00292
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
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;
00311 PetscInt numRefinements = 0;
00312 unsigned int dof =1;
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
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
00369 localSizes[0] = currDa->getIdxElementBegin();
00370
00371
00372 localSizes[2] = currDa->getElementSize();
00373
00374 unsigned int localTotalSize = currDa->getLocalBufferSize();
00375
00376
00377 localSizes[1] = (localTotalSize - (currDa->getIdxPostGhostBegin()));
00378
00379
00380 localSizes[3] = currDa->getIndependentSize();
00381
00382
00383 localSizes[4] = currDa->getPreGhostElementSize();
00384
00385
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 }
00406 }
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
00419
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 }
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
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
00515
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 }
00522
00523 linOct.resize(totalRecv);
00524
00525 sendOffsets[0] = 0;
00526 recvOffsets[0] = 0;
00527
00528
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 }
00533
00534
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
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 }
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
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
00655
00656 for(int i = 1; i < (nlevels-1); i++) {
00657 ot::DAMGCreateJMatrix(damg[i], CreateJacobianHandle);
00658 }
00659 } else {
00660 assert(false);
00661 }
00662
00663
00664
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 }
00792