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
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;
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 + 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
00071 ot::getPrivateMatricesForKSP_Shell = getPrivateMatricesForKSP_Shell_Jac2;
00072
00073 solve_neumann(
00074
00075 pts,
00076 CalcConstCoef,
00077 CalcVarRHS,
00078 30,
00079
00080
00081 sol,
00082 damg
00083 );
00084
00085
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
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
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
00194 DAMGDestroy(damg);
00195
00196 ot::DAMG_Finalize();
00197 PetscFinalize();
00198 return 0;
00199 }
00200