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

omgNeumann_ex3.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 CalcVarCoef(const std::vector<double> & pts, std::vector<double> & coef)
00034 {
00035   for(int i=0,j=0; i<pts.size(); i+=3,j+=2)
00036   {
00037     double s = pts[i]*pts[i] + pts[i+1]*pts[i+1] + pts[i+2]*pts[i+2];
00038     coef[j] = s+1;
00039     coef[j+1] = s ;
00040   }
00041 }
00042 
00043 void CalcVarRHS(const std::vector<double> & pts, std::vector<double> & rhs)
00044 {
00045   for(int i=0,j=0; i<pts.size(); i+=3,j++)
00046   {
00047     double s = pts[i]*pts[i] + pts[i+1]*pts[i+1] + pts[i+2]*pts[i+2];
00048     rhs[j] = (s/*mass_factor*/ + (s+1)*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]) + 2*pts[i]*w*M_PI*sin(w*M_PI*pts[i])*cos(w*M_PI*pts[i+1])*cos(w*M_PI*pts[i+2]) + 2*pts[i+1]*w*M_PI*cos(w*M_PI*pts[i])*sin(w*M_PI*pts[i+1])*cos(w*M_PI*pts[i+2]) + 2*pts[i+2]*w*M_PI*cos(w*M_PI*pts[i])*cos(w*M_PI*pts[i+1])*sin(w*M_PI*pts[i+2]);
00049   }
00050 }
00051 
00052 int main(int argc, char ** argv )
00053 {
00054   PetscInitialize(&argc,&argv,"options","DENDRO example\n");
00055   ot::RegisterEvents();
00056 
00057   ot::DAMG_Initialize(MPI_COMM_WORLD);
00058 
00059   int size, rank;
00060   MPI_Comm_size(MPI_COMM_WORLD,&size);
00061   MPI_Comm_rank(MPI_COMM_WORLD,&rank);
00062 
00063   if (argc<2) 
00064   {
00065     std::cout<<"Usage: mpirun -np <numProcs> "<<argv[0]<<" <filePrefix>"<<std::endl;
00066     return(-1);
00067   }
00068 
00069   std::ostringstream fName;
00070   fName<<argv[1]<<rank<<"_"<<size<<".pts";
00071 
00072   std::vector<double> pts;
00073   ot::readPtsFromFile(const_cast<char*>(fName.str().c_str()), pts);
00074   
00075   Vec sol;
00076   ot::DAMG * damg;
00077   
00078   // 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 
00079   ot::getPrivateMatricesForKSP_Shell = getPrivateMatricesForKSP_Shell_Jac2;
00080 
00081   solve_neumann(
00082     /* input parameters: */
00083      pts,
00084      CalcVarCoef,
00085      CalcVarRHS,
00086      30, /* upper bound */
00087 
00088      /* output parameters */
00089      sol,
00090      damg
00091     );
00092  
00093   // compare solution with the exact one
00094   unsigned numMultigridLevels = damg[0]->totalLevels;
00095   ot::DA* da=damg[numMultigridLevels-1]->da;
00096   PetscScalar *solArray;
00097   da->vecGetBuffer(sol,solArray,false,false,true,1);
00098   double maxDiff = 0;
00099   unsigned maxDepth=da->getMaxDepth()-1;
00100 
00101   if(da->iAmActive())
00102   { 
00103     // first loop over those octants which do not contain ghost nodes.
00104     for(da->init<ot::DA_FLAGS::INDEPENDENT>(); da->curr() < da->end<ot::DA_FLAGS::INDEPENDENT>(); da->next<ot::DA_FLAGS::INDEPENDENT>())  
00105     {
00106       Point pt = da->getCurrentOffset();
00107       unsigned levelhere = da->getLevel(da->curr()) - 1;
00108       double hxOct = ldexp(1.0,-levelhere);
00109       double x = ldexp(pt.x(),-maxDepth );
00110       double y = ldexp(pt.y(),-maxDepth );
00111       double z = ldexp(pt.z(),-maxDepth );
00112 
00113       unsigned int indices[8];
00114       da->getNodeIndices(indices); 
00115 
00116       double coord[8][3] = {
00117         {0.0,0.0,0.0},
00118         {1.0,0.0,0.0},
00119         {0.0,1.0,0.0},
00120         {1.0,1.0,0.0},
00121         {0.0,0.0,1.0},
00122         {1.0,0.0,1.0},
00123         {0.0,1.0,1.0},
00124         {1.0,1.0,1.0}
00125       };
00126 
00127       unsigned char hn = da->getHangingNodeIndex(da->curr());
00128 
00129       for(int i = 0; i < 8; i++)
00130         if (!(hn & (1 << i)))
00131         {
00132           double xhere, yhere, zhere;
00133           xhere = x + coord[i][0]*hxOct; 
00134           yhere = y + coord[i][1]*hxOct; 
00135           zhere = z + coord[i][2]*hxOct; 
00136 
00137           double currDiff = fabs ( solArray[indices[i]] - cos(w*M_PI*xhere)*cos(w*M_PI*yhere)*cos(w*M_PI*zhere) );      
00138           if (maxDiff < currDiff)
00139             maxDiff=currDiff;
00140         }
00141     }  
00142     
00143     // 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
00144     size_t myFirstNode = da->getIdxElementBegin();
00145     size_t postFirstNode = da->getIdxPostGhostBegin();
00146     for(da->init<ot::DA_FLAGS::DEPENDENT>(); da->curr() < da->end<ot::DA_FLAGS::DEPENDENT>(); da->next<ot::DA_FLAGS::DEPENDENT>())  
00147     {
00148       Point pt = da->getCurrentOffset();
00149       unsigned levelhere = da->getLevel(da->curr()) - 1;
00150       double hxOct = ldexp(1.0,-levelhere);
00151       double x = ldexp(pt.x(),-maxDepth );
00152       double y = ldexp(pt.y(),-maxDepth );
00153       double z = ldexp(pt.z(),-maxDepth );
00154 
00155       unsigned int indices[8];
00156       da->getNodeIndices(indices); 
00157 
00158       double coord[8][3] = {
00159         {0.0,0.0,0.0},
00160         {1.0,0.0,0.0},
00161         {0.0,1.0,0.0},
00162         {1.0,1.0,0.0},
00163         {0.0,0.0,1.0},
00164         {1.0,0.0,1.0},
00165         {0.0,1.0,1.0},
00166         {1.0,1.0,1.0}
00167       };
00168 
00169       unsigned char hn = da->getHangingNodeIndex(da->curr());
00170 
00171       for(int i = 0; i < 8; i++)
00172         if ( indices[i]>=myFirstNode && indices[i]<postFirstNode && !(hn & (1 << i)) )
00173         {
00174           double xhere, yhere, zhere;
00175           xhere = x + coord[i][0]*hxOct; 
00176           yhere = y + coord[i][1]*hxOct; 
00177           zhere = z + coord[i][2]*hxOct; 
00178 
00179           double currDiff = fabs ( solArray[indices[i]] - cos(w*M_PI*xhere)*cos(w*M_PI*yhere)*cos(w*M_PI*zhere) );      
00180           if (maxDiff < currDiff)
00181              maxDiff=currDiff;
00182         }
00183     }  
00184   }
00185   da->vecRestoreBuffer(sol,solArray,false,false,true,1); 
00186 
00187   double gMaxDiff;
00188   par::Mpi_Reduce<double>(&maxDiff, &gMaxDiff, 1, MPI_MAX, 0, MPI_COMM_WORLD);  
00189   if (!rank)
00190   {
00191     std::cout<<"Maximum difference: "<<gMaxDiff<<std::endl;
00192   }
00193     
00194   double solmax,solmin;
00195   VecMax(sol,NULL,&solmax);
00196   VecMin(sol,NULL,&solmin);
00197 
00198   if (!rank)
00199     std::cout<<"solution min: " << solmin << "; solution max: " << solmax << std::endl;
00200 
00201   // destroy multigrid context; this destroys the solution as well 
00202   DAMGDestroy(damg);
00203   
00204   ot::DAMG_Finalize();
00205   PetscFinalize();
00206   return 0;
00207 }
00208 

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