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 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 + (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
00079 ot::getPrivateMatricesForKSP_Shell = getPrivateMatricesForKSP_Shell_Jac2;
00080
00081 solve_neumann(
00082
00083 pts,
00084 CalcVarCoef,
00085 CalcVarRHS,
00086 30,
00087
00088
00089 sol,
00090 damg
00091 );
00092
00093
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
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
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
00202 DAMGDestroy(damg);
00203
00204 ot::DAMG_Finalize();
00205 PetscFinalize();
00206 return 0;
00207 }
00208