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
00019 #ifdef MPI_WTIME_IS_GLOBAL
00020 #undef MPI_WTIME_IS_GLOBAL
00021 #endif
00022
00023 #ifdef PETSC_USE_LOG
00024
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 unsigned int maxDepth=30;
00052 bool compressLut=false;
00053 double mgLoadFac = 2.0;
00054 bool incCorner = 1;
00055
00056 PetscInitialize(&argc, &argv,"optionsFBM",help);
00057 ot::RegisterEvents();
00058
00059 ot::DAMG_Initialize(MPI_COMM_WORLD);
00060
00061 MPI_Comm_size(MPI_COMM_WORLD,&npes);
00062 MPI_Comm_rank(MPI_COMM_WORLD,&rank);
00063
00064 std::vector<ot::TreeNode> balOct;
00065 if(argc > 1) {
00066
00067 char pFile[256];
00068 sprintf(pFile, "%s%d_%d.pts", argv[1], rank, npes);
00069
00070 std::vector<double> pts;
00071 ot::readPtsFromFile(pFile, pts);
00072
00073 double gSize[3];
00074 gSize[0] = 1.0;
00075 gSize[1] = 1.0;
00076 gSize[2] = 1.0;
00077
00078 std::vector<ot::TreeNode> linOct;
00079 ot::points2Octree(pts, gSize, linOct, dim, maxDepth, maxNumPts, MPI_COMM_WORLD);
00080 pts.clear();
00081
00082 ot::balanceOctree(linOct, balOct, dim, maxDepth, incCorner, MPI_COMM_WORLD, NULL, NULL);
00083 linOct.clear();
00084
00085 } else {
00086
00087 PetscInt regLev = 4;
00088 PetscOptionsGetInt(0, "-regLev", ®Lev, 0);
00089 createRegularOctree(balOct, regLev, dim, maxDepth, MPI_COMM_WORLD);
00090
00091 }
00092
00093 PetscInt numRefinements = 0;
00094
00095 PetscOptionsGetInt(0,"-numRefinements",&numRefinements,0);
00096 for(int i = 0; i < numRefinements; i++) {
00097 std::vector<ot::TreeNode> tmpOct = balOct;
00098 balOct.clear();
00099 ot::refineOctree(tmpOct, balOct);
00100 }
00101
00102
00103
00104 ot::DAMG *damg;
00105 int nlevels = 1;
00106 unsigned int dof =1;
00107
00108 PetscInt nlevelsPetscInt = nlevels;
00109 PetscOptionsGetInt(0, "-nlevels", &nlevelsPetscInt, 0);
00110 nlevels = nlevelsPetscInt;
00111
00112 if(!rank) {
00113 std::cout<<"nlevels initial: "<<nlevels<<std::endl;
00114 }
00115
00116 MPI_Barrier(MPI_COMM_WORLD);
00117 double setupStartTime = MPI_Wtime();
00118
00119
00120 ot::DAMGCreateAndSetDA(PETSC_COMM_WORLD, nlevels, NULL, &damg,
00121 balOct, dof, mgLoadFac, compressLut, incCorner);
00122
00123 MPI_Barrier(MPI_COMM_WORLD);
00124 double setupEndTime = MPI_Wtime();
00125
00126 if(!rank) {
00127 std::cout<<"nlevels final: "<<nlevels<<std::endl;
00128 }
00129
00130 if(!rank) {
00131 std::cout << "Created DA for all levels."<< std::endl;
00132 }
00133
00134 ot::PrintDAMG(damg);
00135
00136 createLmatType2(LaplacianType2Stencil);
00137 createMmatType2(MassType2Stencil);
00138 createShFnMat(ShapeFnStencil);
00139
00140 ot::DAMGCreateSuppressedDOFs(damg);
00141
00142 SetDirichletJacContexts(damg);
00143
00144
00145
00146 ot::getPrivateMatricesForKSP_Shell = getPrivateMatricesForKSP_Shell_DirichletJac;
00147
00148 ot::DAMGSetKSP(damg, CreateDirichletJacobian, ComputeDirichletJacobian, ComputeTestFBM_RHS);
00149
00150 MPI_Barrier(MPI_COMM_WORLD);
00151 double solveStartTime = MPI_Wtime();
00152
00153 ot::DAMGSolve(damg);
00154
00155 MPI_Barrier(MPI_COMM_WORLD);
00156 double solveEndTime = MPI_Wtime();
00157
00158 Vec solTrue;
00159 VecDuplicate(DAMGGetx(damg), &solTrue);
00160 SetSolutionFBM(damg[nlevels - 1], solTrue);
00161 EnforceZeroFBM(damg[nlevels - 1], DAMGGetx(damg));
00162
00163 VecAXPY(solTrue, -1.0, DAMGGetx(damg));
00164
00165 PetscReal maxNormErr;
00166 VecNorm(solTrue, NORM_INFINITY, &maxNormErr);
00167
00168 VecDestroy(solTrue);
00169
00170 if(!rank) {
00171 std::cout<<" Total Setup Time: "<<(setupEndTime - setupStartTime)<<std::endl;
00172 std::cout<<" Total Solve Time: "<<(solveEndTime - solveStartTime)<<std::endl;
00173 std::cout<<" maxNormErr (Pointwise): "<<maxNormErr<<std::endl;
00174 }
00175
00176 destroyLmatType2(LaplacianType2Stencil);
00177 destroyMmatType2(MassType2Stencil);
00178 destroyShFnMat(ShapeFnStencil);
00179
00180 DestroyDirichletJacContexts(damg);
00181
00182 DAMGDestroy(damg);
00183
00184 balOct.clear();
00185
00186 ot::DAMG_Finalize();
00187
00188 PetscFinalize();
00189 }
00190