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