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

omgNeumann_ex1.C

Go to the documentation of this file.
00001 #include <mpi.h>
00002 #include <cstdio>
00003 #include "oda.h"
00004 #include "omg.h"
00005 #include "Point.h"
00006 #include "parUtils.h"
00007 #include "octUtils.h"
00008 #include "TreeNode.h"
00009 #include "handleStencils.h"
00010 #include <cstdlib>
00011 #include "sys.h"
00012 #include "externVars.h"
00013 #include <sstream>
00014 #include "omgNeumann.h"
00015 
00016 //Don't want time to be synchronized. Need to check load imbalance.
00017 #ifdef MPI_WTIME_IS_GLOBAL
00018 #undef MPI_WTIME_IS_GLOBAL
00019 #endif
00020 
00021 #ifdef PETSC_USE_LOG
00022 
00023 int Jac2DiagEvent;
00024 int Jac2MultEvent;
00025 int Jac2FinestDiagEvent;
00026 int Jac2FinestMultEvent;
00027 
00028 #endif
00029 
00030 const double w = 1;  // frequency for test purposes
00031 
00032 
00033 void CalcConstCoef(const std::vector<double> & pts, std::vector<double> & coef)
00034 {
00035   coef.assign(coef.size(),1.0);
00036 }
00037 
00038 void CalcVarRHS(const std::vector<double> & pts, std::vector<double> & rhs)
00039 {
00040   for(int i=0,j=0; i<pts.size(); i+=3,j++)
00041    rhs[j] = (1 /*mass_factor*/ + 3*w*w*M_PI*M_PI) * cos(w*M_PI*pts[i])* cos(w*M_PI*pts[i+1])*cos(w*M_PI*pts[i+2]);      
00042 }
00043 
00044 int main(int argc, char ** argv )
00045 {
00046   PetscInitialize(&argc,&argv,"options","DENDRO example\n");
00047   ot::RegisterEvents();
00048 
00049   ot::DAMG_Initialize(MPI_COMM_WORLD);
00050 
00051   int size, rank;
00052   MPI_Comm_size(MPI_COMM_WORLD,&size);
00053   MPI_Comm_rank(MPI_COMM_WORLD,&rank);
00054 
00055   if (argc<2) 
00056   {
00057     std::cout<<"Usage: mpirun -np <numProcs> "<<argv[0]<<" <filePrefix>"<<std::endl;
00058     return(-1);
00059   }
00060 
00061   std::ostringstream fName;
00062   fName<<argv[1]<<rank<<"_"<<size<<".pts";
00063 
00064   std::vector<double> pts;
00065   ot::readPtsFromFile(const_cast<char*>(fName.str().c_str()), pts);
00066   
00067   Vec sol;
00068   ot::DAMG * damg;
00069   
00070   // This function pointer will be called while setting up the solver on the coarsest grid if not all processes are active on the coarse grid 
00071   ot::getPrivateMatricesForKSP_Shell = getPrivateMatricesForKSP_Shell_Jac2;
00072 
00073   solve_neumann(
00074     /* input parameters: */
00075      pts,
00076      CalcConstCoef,
00077      CalcVarRHS,
00078      30, /* upper bound */
00079 
00080      /* output parameters */
00081      sol,
00082      damg
00083     );
00084  
00085   // compare solution with the exact one
00086   unsigned numMultigridLevels = damg[0]->totalLevels;
00087   ot::DA* da=damg[numMultigridLevels-1]->da;
00088   PetscScalar *solArray;
00089   da->vecGetBuffer(sol,solArray,false,false,true,1);
00090   double maxDiff = 0;
00091   unsigned maxDepth=da->getMaxDepth()-1;
00092 
00093   if(da->iAmActive())
00094   { 
00095     // first loop over those octants which do not contain ghost nodes.
00096     for(da->init<ot::DA_FLAGS::INDEPENDENT>(); da->curr() < da->end<ot::DA_FLAGS::INDEPENDENT>(); da->next<ot::DA_FLAGS::INDEPENDENT>())  
00097     {
00098       Point pt = da->getCurrentOffset();
00099       unsigned levelhere = da->getLevel(da->curr()) - 1;
00100       double hxOct = ldexp(1.0,-levelhere);
00101       double x = ldexp(pt.x(),-maxDepth );
00102       double y = ldexp(pt.y(),-maxDepth );
00103       double z = ldexp(pt.z(),-maxDepth );
00104 
00105       unsigned int indices[8];
00106       da->getNodeIndices(indices); 
00107 
00108       double coord[8][3] = {
00109         {0.0,0.0,0.0},
00110         {1.0,0.0,0.0},
00111         {0.0,1.0,0.0},
00112         {1.0,1.0,0.0},
00113         {0.0,0.0,1.0},
00114         {1.0,0.0,1.0},
00115         {0.0,1.0,1.0},
00116         {1.0,1.0,1.0}
00117       };
00118 
00119       unsigned char hn = da->getHangingNodeIndex(da->curr());
00120 
00121       for(int i = 0; i < 8; i++)
00122         if (!(hn & (1 << i)))
00123         {
00124           double xhere, yhere, zhere;
00125           xhere = x + coord[i][0]*hxOct; 
00126           yhere = y + coord[i][1]*hxOct; 
00127           zhere = z + coord[i][2]*hxOct; 
00128 
00129           double currDiff = fabs ( solArray[indices[i]] - cos(w*M_PI*xhere)*cos(w*M_PI*yhere)*cos(w*M_PI*zhere) );      
00130           if (maxDiff < currDiff)
00131             maxDiff=currDiff;
00132         }
00133     }  
00134     
00135     // now loop over those octants which DO contain ghost nodes; now we need to  check if the node is ghost, and skip it if it is
00136     size_t myFirstNode = da->getIdxElementBegin();
00137     size_t postFirstNode = da->getIdxPostGhostBegin();
00138     for(da->init<ot::DA_FLAGS::DEPENDENT>(); da->curr() < da->end<ot::DA_FLAGS::DEPENDENT>(); da->next<ot::DA_FLAGS::DEPENDENT>())  
00139     {
00140       Point pt = da->getCurrentOffset();
00141       unsigned levelhere = da->getLevel(da->curr()) - 1;
00142       double hxOct = ldexp(1.0,-levelhere);
00143       double x = ldexp(pt.x(),-maxDepth );
00144       double y = ldexp(pt.y(),-maxDepth );
00145       double z = ldexp(pt.z(),-maxDepth );
00146 
00147       unsigned int indices[8];
00148       da->getNodeIndices(indices); 
00149 
00150       double coord[8][3] = {
00151         {0.0,0.0,0.0},
00152         {1.0,0.0,0.0},
00153         {0.0,1.0,0.0},
00154         {1.0,1.0,0.0},
00155         {0.0,0.0,1.0},
00156         {1.0,0.0,1.0},
00157         {0.0,1.0,1.0},
00158         {1.0,1.0,1.0}
00159       };
00160 
00161       unsigned char hn = da->getHangingNodeIndex(da->curr());
00162 
00163       for(int i = 0; i < 8; i++)
00164         if ( indices[i]>=myFirstNode && indices[i]<postFirstNode && !(hn & (1 << i)) )
00165         {
00166           double xhere, yhere, zhere;
00167           xhere = x + coord[i][0]*hxOct; 
00168           yhere = y + coord[i][1]*hxOct; 
00169           zhere = z + coord[i][2]*hxOct; 
00170 
00171           double currDiff = fabs ( solArray[indices[i]] - cos(w*M_PI*xhere)*cos(w*M_PI*yhere)*cos(w*M_PI*zhere) );      
00172           if (maxDiff < currDiff)
00173              maxDiff=currDiff;
00174         }
00175     }  
00176   }
00177   da->vecRestoreBuffer(sol,solArray,false,false,true,1); 
00178 
00179     double gMaxDiff;
00180     par::Mpi_Reduce<double>(&maxDiff, &gMaxDiff, 1, MPI_MAX, 0, MPI_COMM_WORLD);  
00181   if (!rank)
00182   {
00183     std::cout<<"Maximum difference: "<<gMaxDiff<<std::endl;
00184   }
00185     
00186   double solmax,solmin;
00187   VecMax(sol,NULL,&solmax);
00188   VecMin(sol,NULL,&solmin);
00189 
00190   if (!rank)
00191     std::cout<<"solution min: " << solmin << "; solution max: " << solmax << std::endl;
00192 
00193   // destroy multigrid context; this destroys the solution as well 
00194   DAMGDestroy(damg);
00195   
00196   ot::DAMG_Finalize();
00197   PetscFinalize();
00198   return 0;
00199 }
00200 

Generated on Tue Mar 24 16:13:59 2009 for DENDRO by  doxygen 1.3.9.1