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

tstFBM2.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 <cstdlib>
00012 #include "colors.h"
00013 #include "externVars.h"
00014 #include "dendro.h"
00015 
00016 static char help[] = "FBM";
00017 
00018 //Don't want time to be synchronized. Need to check load imbalance.
00019 #ifdef MPI_WTIME_IS_GLOBAL
00020 #undef MPI_WTIME_IS_GLOBAL
00021 #endif
00022 
00023 #ifdef PETSC_USE_LOG
00024 //user-defined variables
00025 int Jac1DiagEvent;
00026 int Jac1MultEvent;
00027 int Jac1FinestDiagEvent;
00028 int Jac1FinestMultEvent;
00029 
00030 int Jac2DiagEvent;
00031 int Jac2MultEvent;
00032 int Jac2FinestDiagEvent;
00033 int Jac2FinestMultEvent;
00034 
00035 int Jac3DiagEvent;
00036 int Jac3MultEvent;
00037 int Jac3FinestDiagEvent;
00038 int Jac3FinestMultEvent;
00039 #endif
00040 
00041 double***** LaplacianType1Stencil; 
00042 double**** LaplacianType2Stencil; 
00043 double***** MassType1Stencil; 
00044 double**** MassType2Stencil; 
00045 double****** ShapeFnStencil;
00046 
00047 int main(int argc, char ** argv ) {     
00048   int npes, rank;
00049   unsigned int maxNumPts = 1;
00050   unsigned int dim = 3;
00051   bool compressLut = false;
00052   double mgLoadFac = 2.0;
00053   bool incCorner = 1;  
00054 
00055   PetscInitialize(&argc,&argv,"optionsFBM2",help);
00056   ot::RegisterEvents();
00057 
00058   ot::DAMG_Initialize(MPI_COMM_WORLD);
00059 
00060   MPI_Comm_size(MPI_COMM_WORLD, &npes);
00061   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
00062 
00063   PetscReal gamma = 0;
00064   PetscOptionsGetReal(0, "-gamma", &gamma, 0);
00065 
00066   PetscInt Nsample = 100;
00067   PetscOptionsGetInt(0, "-Nsample", &Nsample, 0);
00068 
00069   PetscInt maxDepth = 30;
00070   PetscOptionsGetInt(0, "-maxDepth", &maxDepth, 0);
00071 
00072   double hSample = 1.0/static_cast<double>(Nsample);
00073   std::vector<double> pts;
00074   if(!rank) {
00075     for(double z = 0; z < 1.0; z += hSample) {
00076       for(double y = 0; y < 1.0; y += hSample) {
00077         pts.push_back(gamma);
00078         pts.push_back(y);
00079         pts.push_back(z);
00080       }
00081     }
00082   }
00083 
00084   double gSize[3];
00085   gSize[0] = 1.0;
00086   gSize[1] = 1.0;
00087   gSize[2] = 1.0;
00088 
00089   //Construction
00090   std::vector<ot::TreeNode> linOct;
00091   ot::points2Octree(pts, gSize, linOct, dim, maxDepth, maxNumPts, MPI_COMM_WORLD);
00092   pts.clear();
00093 
00094   //Balancing...
00095   std::vector<ot::TreeNode> balOct;
00096   ot::balanceOctree(linOct, balOct, dim, maxDepth, incCorner, MPI_COMM_WORLD, NULL, NULL);
00097   linOct.clear();
00098 
00099   PetscInt       numRefinements = 0;
00100 
00101   PetscOptionsGetInt(0,"-numRefinements",&numRefinements,0);
00102   for(int i = 0; i < numRefinements; i++) {
00103     std::vector<ot::TreeNode> tmpOct = balOct;
00104     balOct.clear();
00105     ot::refineOctree(tmpOct, balOct); 
00106   }
00107 
00108   //Solve ...
00109 
00110   ot::DAMG        *damg;    
00111   int       nlevels = 1; //number of multigrid levels
00112   unsigned int       dof =1;// degrees of freedom per node  
00113 
00114   PetscInt nlevelsPetscInt = nlevels;
00115   PetscOptionsGetInt(0, "-nlevels", &nlevelsPetscInt, 0);
00116   nlevels = nlevelsPetscInt;
00117 
00118   if(!rank) {
00119     std::cout<<"nlevels initial: "<<nlevels<<std::endl;
00120   }
00121 
00122   MPI_Barrier(MPI_COMM_WORLD);
00123   double setupStartTime = MPI_Wtime();
00124 
00125   // Note: The user context for all levels will be set separately later.
00126   ot::DAMGCreateAndSetDA(PETSC_COMM_WORLD, nlevels, NULL, &damg,
00127       balOct, dof, mgLoadFac, compressLut, incCorner);
00128 
00129   MPI_Barrier(MPI_COMM_WORLD);
00130   double setupEndTime = MPI_Wtime();
00131 
00132   if(!rank) {
00133     std::cout<<"nlevels final: "<<nlevels<<std::endl;
00134   }
00135 
00136   if(!rank) {
00137     std::cout << "Created DA for all levels."<< std::endl;
00138   }
00139 
00140   ot::PrintDAMG(damg);
00141 
00142   createLmatType2(LaplacianType2Stencil);
00143   createMmatType2(MassType2Stencil);
00144   createShFnMat(ShapeFnStencil);
00145 
00146   ot::DAMGCreateSuppressedDOFs(damg);
00147 
00148   SetDirichletJacContexts(damg);
00149 
00150   //Global Function Handles for using KSP_Shell (will be used @ the coarsest grid if not all
00151   //processors are active on the coarsest grid)
00152   ot::getPrivateMatricesForKSP_Shell = getPrivateMatricesForKSP_Shell_DirichletJac;
00153 
00154   ot::DAMGSetKSP(damg, CreateDirichletLaplacian, ComputeDirichletLaplacian, ComputeFBM2_RHS);
00155 
00156   MPI_Barrier(MPI_COMM_WORLD);
00157   double solveStartTime = MPI_Wtime();
00158 
00159   ot::DAMGSolve(damg);
00160 
00161   MPI_Barrier(MPI_COMM_WORLD);
00162   double solveEndTime = MPI_Wtime();
00163 
00164   Vec solTrue;
00165   VecDuplicate(DAMGGetx(damg), &solTrue);
00166   SetSolutionFBM2(damg[nlevels - 1], solTrue);
00167   EnforceZeroFBM2(damg[nlevels - 1], DAMGGetx(damg));
00168 
00169   VecAXPY(solTrue, -1.0, DAMGGetx(damg));
00170 
00171   PetscReal maxNormErr;
00172   VecNorm(solTrue, NORM_INFINITY, &maxNormErr);
00173 
00174   VecDestroy(solTrue);
00175 
00176   if(!rank) {
00177     std::cout<<" Total Setup Time: "<<(setupEndTime - setupStartTime)<<std::endl;
00178     std::cout<<" Total Solve Time: "<<(solveEndTime - solveStartTime)<<std::endl;
00179     std::cout<<" maxNormErr (Pointwise): "<<maxNormErr<<std::endl;
00180   }
00181 
00182   destroyLmatType2(LaplacianType2Stencil);
00183   destroyMmatType2(MassType2Stencil);
00184   destroyShFnMat(ShapeFnStencil);
00185 
00186   DestroyDirichletJacContexts(damg);
00187 
00188   DAMGDestroy(damg);
00189 
00190   balOct.clear();
00191 
00192   ot::DAMG_Finalize();
00193 
00194   PetscFinalize();
00195 }//end function
00196 

Generated on Tue Mar 24 16:14:00 2009 for DENDRO by  doxygen 1.3.9.1