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 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
00090 std::vector<ot::TreeNode> linOct;
00091 ot::points2Octree(pts, gSize, linOct, dim, maxDepth, maxNumPts, MPI_COMM_WORLD);
00092 pts.clear();
00093
00094
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
00109
00110 ot::DAMG *damg;
00111 int nlevels = 1;
00112 unsigned int dof =1;
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
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
00151
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 }
00196