fmm3d_setup_mpi.cpp

Go to the documentation of this file.
00001 /* Kernel Independent Fast Multipole Method
00002    Copyright (C) 2004 Lexing Ying, New York University
00003 
00004 This program is free software; you can redistribute it and/or modify
00005 it under the terms of the GNU General Public License as published by
00006 the Free Software Foundation; either version 2, or (at your option)
00007 any later version.
00008 
00009 This program is distributed in the hope that it will be useful, but WITHOUT
00010 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00011 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00012 for more details.
00013 
00014 You should have received a copy of the GNU General Public License
00015 along with this program; see the file COPYING.  If not, write to the Free
00016 Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
00017 02111-1307, USA.  */
00018 #include "fmm3d_mpi.hpp"
00019 
00020 using std::cerr;
00021 using std::endl;
00022 
00023 // ----------------------------------------------------------------------
00024 /* FMM setup sets up all of the essential information needed for
00025  * running the fast multipole 3d mpi implementation.  A major
00026  * componenet is the setup of the local essential tree, which
00027  * builds the octree, source data requirements and target data
00028  * requirements
00029  */
00030 #undef __FUNCT__
00031 #define __FUNCT__ "FMM3d_MPI::setup"
00032 int FMM3d_MPI::setup()
00033 {
00034   //begin
00035   //ebiLogInfo( "setup............."); 
00036   //-----------------------------------------------------
00037   PetscTruth flg = PETSC_FALSE;
00038   pC( PetscOptionsGetInt(prefix().c_str(), "-np",     &_np,     &flg) );  pA(flg==true);
00039   
00040   //-----------------------------------------------------
00041   pA(_srcPos!=NULL && _srcNor!=NULL && _trgPos!=NULL);
00042   /* Build a new Communication Object for the local essential tree (LET)
00043         * and run setup on it.  See let3d_mpi.hpp and let3d_mpi.cpp for
00044         * more information */
00045   _let = new Let3d_MPI(prefix()+"let3d_");
00046   _let->srcPos()=_srcPos;
00047   _let->trgPos()=_trgPos;
00048   _let->ctr()=_ctr;
00049   _let->rootLevel()=_rootLevel; 
00050 
00051   clock_t ck0, ck1;
00052   ck0 = clock();
00053   pC( _let->setup() );  
00054   
00055   //2. decide _eq_mm and _mul_mm, and get matmgnt based on that
00056   switch(_knl.kernelType()) {
00057          //laplace kernels
00058   case KNL_LAP_S_U: _knl_mm = Kernel3d_MPI(KNL_LAP_S_U, _knl.coefs()); break;
00059   case KNL_LAP_D_U: _knl_mm = Kernel3d_MPI(KNL_LAP_S_U, _knl.coefs()); break;
00060          //stokes kernels
00061   case KNL_STK_S_U: _knl_mm = Kernel3d_MPI(KNL_STK_F_U, _knl.coefs()); break;
00062   case KNL_STK_S_P: _knl_mm = Kernel3d_MPI(KNL_LAP_S_U, vector<double>()); break;
00063   case KNL_STK_D_U: _knl_mm = Kernel3d_MPI(KNL_STK_F_U, _knl.coefs()); break;
00064   case KNL_STK_D_P: _knl_mm = Kernel3d_MPI(KNL_LAP_S_U, vector<double>()); break;
00065          //navier kernels
00066   case KNL_NAV_S_U: _knl_mm = Kernel3d_MPI(KNL_NAV_S_U, _knl.coefs()); break;
00067   case KNL_NAV_D_U: _knl_mm = Kernel3d_MPI(KNL_NAV_S_U, _knl.coefs()); break;
00068   default: pA(0);
00069   }
00070   
00071   _matmgnt  = MatMgnt3d_MPI::getmmptr(_knl_mm, _np);
00072 
00073   //3. self setup
00074   _nodeVec.resize( _let->nodeVec().size() );
00075     // ----------------------------------------------------------------------
00076   /* In the following scope, use the LET to build the global data vectors
00077         *_glbSrcExaPos, _glbSrcExaNor, _glbSrcExaDen, _glbSrcUpwEquDen
00078         * These vectors are globally available, so we use VecCreateMPI.
00079         * See http://www-unix.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-current/docs/manualpages/Vec/VecCreateMPI.html
00080         * for more information.
00081         */
00082   {
00083          //begin
00084          //ebiLogInfo( "setup...gdata");
00085          //1. create cvecs
00086          int lclGlbSrcNodeCnt = _let->lclGlbSrcNodeCnt();
00087          int lclGlbSrcExaCnt = _let->lclGlbSrcExaCnt();
00088          pC( VecCreateMPI(mpiComm(), lclGlbSrcExaCnt*dim(),         PETSC_DETERMINE, &_glbSrcExaPos) );
00089          pC( VecCreateMPI(mpiComm(), lclGlbSrcExaCnt*dim(),         PETSC_DETERMINE, &_glbSrcExaNor) );  
00090          pC( VecCreateMPI(mpiComm(), lclGlbSrcExaCnt*srcDOF(),        PETSC_DETERMINE, &_glbSrcExaDen) );
00091          pC( VecCreateMPI(mpiComm(), lclGlbSrcNodeCnt*datSze(UE),  PETSC_DETERMINE, &_glbSrcUpwEquDen) );
00092   }
00093   // ----------------------------------------------------------------------
00094   /* Build contributor data for a specific processor */
00095   {
00096          //begin
00097          //ebiLogInfo("setup...cdata");
00098          //1. create vecs
00099          int ctbSrcNodeCnt = _let->ctbSrcNodeCnt();
00100          int ctbSrcExaCnt = _let->ctbSrcExaCnt();
00101          /* Create standard sequential array.  See the following for more information:
00102           * http://www-unix.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-current/docs/manualpages/Vec/VecCreateSeq.html
00103           * The following creates sequential vectors on each processor for the contributor positions, normals, etc.
00104           * based on the sizes of contributor data found when building the LET and figuring out how many nodes this
00105           * processor contributes to.
00106           */
00107          pC( VecCreateSeq(PETSC_COMM_SELF, ctbSrcExaCnt*dim(),  &_ctbSrcExaPos) );
00108          pC( VecCreateSeq(PETSC_COMM_SELF, ctbSrcExaCnt*dim(),  &_ctbSrcExaNor) );
00109          pC( VecCreateSeq(PETSC_COMM_SELF, ctbSrcExaCnt*srcDOF(), &_ctbSrcExaDen) );
00110          pC( VecCreateSeq(PETSC_COMM_SELF, ctbSrcNodeCnt*datSze(UE), &_ctbSrcUpwEquDen) );
00111          pC( VecCreateSeq(PETSC_COMM_SELF, ctbSrcNodeCnt*datSze(UC), &_ctbSrcUpwChkVal) );
00112   
00113          /* 2. create scatters
00114           * Index Sets are used to setup the vector scatters
00115           * See http://www-unix.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-current/docs/manualpages/IS/index.html
00116           * for more information */
00117          vector<int>& ctb2GlbSrcNodeMap = _let->ctb2GlbSrcNodeMap();
00118          vector<int>& ctb2GlbSrcExaMap = _let->ctb2GlbSrcExaMap();
00119          /* self index set */
00120          IS selfis;
00121          /* Combine index set */
00122          IS combis;
00123          /* Creates a data structure for an index set containing a list of evenly spaced integers.
00124           * Starts at 0 and increments by 1 until size ctbSrcExaCnt*dim().  Store in selfis.  See:
00125           * http://www-unix.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-current/docs/manualpages/IS/ISCreateStride.html
00126           */
00127          pC( ISCreateStride(PETSC_COMM_SELF, ctbSrcExaCnt*dim(), 0, 1, &selfis) );
00128          for(int i=0; i<ctb2GlbSrcExaMap.size(); i++){
00129                 ctb2GlbSrcExaMap[i]*=dim();
00130          }
00131          /* Creates a data structure for an index set containing a list of integers. The indices are relative to entries, not blocks
00132           * ctb2GlbSrcExaMap.size() = length of index set
00133           * dim() = number of elements in each block
00134           * ctb2GlbSrcExaMap = list of integers
00135           * combis = new index set
00136           * See for more info: 
00137           * http://www-unix.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-current/docs/manualpages/IS/ISCreateBlock.html
00138           */
00139          pC( ISCreateBlock( PETSC_COMM_SELF, dim(), ctb2GlbSrcExaMap.size(), &(ctb2GlbSrcExaMap[0]), &combis) );
00140          for(int i=0; i<ctb2GlbSrcExaMap.size(); i++){
00141                 ctb2GlbSrcExaMap[i]/=dim();
00142          }
00143          /* http://www-unix.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-current/docs/manualpages/Vec/VecScatterCreate.html */
00144          pC( VecScatterCreate(_ctbSrcExaPos, selfis, _glbSrcExaPos, combis, &_ctb2GlbSrcExaPos) );
00145          pC( ISDestroy(selfis) );
00146          pC( ISDestroy(combis) );
00147          
00148          pC( ISCreateStride(PETSC_COMM_SELF, ctbSrcExaCnt*srcDOF(), 0, 1, &selfis) );
00149          for(int i=0; i<ctb2GlbSrcExaMap.size(); i++) ctb2GlbSrcExaMap[i]*=srcDOF();
00150          pC( ISCreateBlock( PETSC_COMM_SELF, srcDOF(),ctb2GlbSrcExaMap.size(), &(ctb2GlbSrcExaMap[0]), &combis) );
00151          for(int i=0; i<ctb2GlbSrcExaMap.size(); i++) ctb2GlbSrcExaMap[i]/=srcDOF();
00152          pC( VecScatterCreate(_ctbSrcExaDen, selfis, _glbSrcExaDen, combis, &_ctb2GlbSrcExaDen) );
00153          pC( ISDestroy(selfis) );
00154          pC( ISDestroy(combis) );
00155          
00156          pC( ISCreateStride(PETSC_COMM_SELF, ctbSrcNodeCnt*datSze(UE),0,1,&selfis) );
00157          for(int i=0; i<ctb2GlbSrcNodeMap.size(); i++) ctb2GlbSrcNodeMap[i]*=datSze(UE);
00158          pC( ISCreateBlock( PETSC_COMM_SELF, datSze(UE), ctb2GlbSrcNodeMap.size(), &(ctb2GlbSrcNodeMap[0]), &combis) );
00159          for(int i=0; i<ctb2GlbSrcNodeMap.size(); i++) ctb2GlbSrcNodeMap[i]/=datSze(UE);
00160          pC( VecScatterCreate(_ctbSrcUpwEquDen, selfis, _glbSrcUpwEquDen, combis, &_ctb2GlbSrcUpwEquDen) );
00161          pC( ISDestroy(selfis) );
00162          pC( ISDestroy(combis) );
00163          
00164          //3. gather the position using the pos scatter
00165          PetscScalar zero=0.0;  pC( VecSet(_ctbSrcExaPos, zero) );
00166          int procLclstart, procLclend; _let->procLclRan(_srcPos, procLclstart, procLclend);
00167          double* parr;      pC( VecGetArray(    _srcPos, &parr) );
00168          double* narr;      pC( VecGetArray(    _srcNor, &narr) );
00169          vector<int> ordvec;  pC( _let->upwOrderCollect(ordvec) );
00170          for(int i=0; i<ordvec.size(); i++) {
00171                 int gNodeIdx = ordvec[i];
00172                 if(_let->node(gNodeIdx).tag() & LET_CBTRNODE) { //contributor
00173                   if(_let->terminal(gNodeIdx)==true) { //terminal cbtr
00174                          DblNumMat ctbSrcExaPos(this->ctbSrcExaPos(gNodeIdx));
00175                          DblNumMat ctbSrcExaNor(this->ctbSrcExaNor(gNodeIdx));
00176                          vector<int>& curVecIdxs = _let->node(gNodeIdx).ctbSrcOwnVecIdxs();
00177                          for(int k=0; k<curVecIdxs.size(); k++) {
00178                                 int poff = curVecIdxs[k] - procLclstart;
00179                                 for(int d=0; d<dim(); d++) {
00180                                   ctbSrcExaPos(d,k) = parr[poff*dim()+d]; //Pos
00181                                   ctbSrcExaNor(d,k) = narr[poff*dim()+d]; //Nor
00182                                 }
00183                          }
00184                   }
00185                 }
00186          }
00187          pC( VecRestoreArray(_srcPos, &parr) );
00188          pC( VecRestoreArray(_srcNor, &narr) );
00189          
00190          pC( VecScatterBegin(_ctbSrcExaPos, _glbSrcExaPos, INSERT_VALUES, SCATTER_FORWARD, _ctb2GlbSrcExaPos) );
00191          pC( VecScatterEnd(  _ctbSrcExaPos, _glbSrcExaPos, INSERT_VALUES, SCATTER_FORWARD, _ctb2GlbSrcExaPos) );
00192          
00193          pC( VecScatterBegin(_ctbSrcExaNor, _glbSrcExaNor, INSERT_VALUES, SCATTER_FORWARD, _ctb2GlbSrcExaPos) );
00194          pC( VecScatterEnd(  _ctbSrcExaNor, _glbSrcExaNor, INSERT_VALUES, SCATTER_FORWARD, _ctb2GlbSrcExaPos) );
00195   }
00196   // ---------------------------------------------------------------------- 
00197   {
00198          //begin
00199          //ebiLogInfo("setup...udata");
00200          //1. create vecs
00201          int usrSrcNodeCnt = _let->usrSrcNodeCnt();
00202          int usrSrcExaCnt = _let->usrSrcExaCnt();
00203          pC( VecCreateSeq(PETSC_COMM_SELF, usrSrcExaCnt*dim(),  &_usrSrcExaPos) );
00204          pC( VecCreateSeq(PETSC_COMM_SELF, usrSrcExaCnt*dim(),  &_usrSrcExaNor) );
00205          pC( VecCreateSeq(PETSC_COMM_SELF, usrSrcExaCnt*srcDOF(), &_usrSrcExaDen) );
00206          pC( VecCreateSeq(PETSC_COMM_SELF, usrSrcNodeCnt*datSze(UE), &_usrSrcUpwEquDen) );
00207          
00208          //2. scatter user -> global
00209          vector<int>& usr2GlbSrcNodeMap = _let->usr2GlbSrcNodeMap();
00210          vector<int>& usr2GlbSrcExaMap = _let->usr2GlbSrcExaMap();
00211          IS selfis;
00212          IS combis;
00213          pC( ISCreateStride(PETSC_COMM_SELF, usrSrcExaCnt*dim(), 0, 1, &selfis) );
00214          for(int i=0; i<usr2GlbSrcExaMap.size(); i++) usr2GlbSrcExaMap[i]*=dim();
00215          pC( ISCreateBlock( PETSC_COMM_SELF, dim(), usr2GlbSrcExaMap.size(), &(usr2GlbSrcExaMap[0]), &combis) );
00216          for(int i=0; i<usr2GlbSrcExaMap.size(); i++) usr2GlbSrcExaMap[i]/=dim();
00217          pC( VecScatterCreate(_usrSrcExaPos, selfis, _glbSrcExaPos, combis, &_usr2GlbSrcExaPos) );
00218          pC( ISDestroy(selfis) );
00219          pC( ISDestroy(combis) );
00220          
00221          pC( ISCreateStride(PETSC_COMM_SELF, usrSrcExaCnt*srcDOF(), 0, 1, &selfis) );
00222          for(int i=0; i<usr2GlbSrcExaMap.size(); i++) usr2GlbSrcExaMap[i]*=srcDOF();
00223          pC( ISCreateBlock( PETSC_COMM_SELF, srcDOF(),usr2GlbSrcExaMap.size(), &(usr2GlbSrcExaMap[0]), &combis) );
00224          for(int i=0; i<usr2GlbSrcExaMap.size(); i++) usr2GlbSrcExaMap[i]/=srcDOF();
00225          pC( VecScatterCreate(_usrSrcExaDen, selfis, _glbSrcExaDen, combis, &_usr2GlbSrcExaDen) );
00226          pC( ISDestroy(selfis) );
00227          pC( ISDestroy(combis) );
00228          
00229          pC( ISCreateStride(PETSC_COMM_SELF, usrSrcNodeCnt*datSze(UE),0,1,&selfis) );
00230          for(int i=0; i<usr2GlbSrcNodeMap.size(); i++) usr2GlbSrcNodeMap[i]*=datSze(UE);
00231          pC( ISCreateBlock( PETSC_COMM_SELF, datSze(UE), usr2GlbSrcNodeMap.size(), &(usr2GlbSrcNodeMap[0]), &combis) );
00232          for(int i=0; i<usr2GlbSrcNodeMap.size(); i++) usr2GlbSrcNodeMap[i]/=datSze(UE);
00233          pC( VecScatterCreate(_usrSrcUpwEquDen, selfis, _glbSrcUpwEquDen, combis, &_usr2GlbSrcUpwEquDen) );
00234          pC( ISDestroy(selfis) );
00235          pC( ISDestroy(combis) );
00236          
00237          //3. do vecscatter
00238          pC( VecScatterBegin(_glbSrcExaPos, _usrSrcExaPos, INSERT_VALUES, SCATTER_REVERSE, _usr2GlbSrcExaPos) );
00239          pC( VecScatterEnd(  _glbSrcExaPos, _usrSrcExaPos, INSERT_VALUES, SCATTER_REVERSE, _usr2GlbSrcExaPos) );
00240          
00241          pC( VecScatterBegin(_glbSrcExaNor, _usrSrcExaNor, INSERT_VALUES, SCATTER_REVERSE, _usr2GlbSrcExaPos) );
00242          pC( VecScatterEnd(  _glbSrcExaNor, _usrSrcExaNor, INSERT_VALUES, SCATTER_REVERSE, _usr2GlbSrcExaPos) );
00243          
00244          vector<int> ordvec; pC( _let->upwOrderCollect(ordvec) );
00245          //4. allocate UNExt* and write data
00246          for(int i=0; i<ordvec.size(); i++) {
00247                 int gNodeIdx = ordvec[i];
00248                 if(_let->node(gNodeIdx).tag() & LET_EVTRNODE) {
00249                   for(vector<int>::iterator vi=_let->node(gNodeIdx).Vnodes().begin(); vi!=_let->node(gNodeIdx).Vnodes().end(); vi++) {
00250                          node(*vi).vLstOthNum() ++;
00251                   }
00252                 }
00253          }
00254   }
00255   // ---------------------------------------------------------------------- 
00256   {
00257          //begin
00258          //ebiLogInfo("setup...edata");
00259          //1. create vecs
00260          int evaTrgNodeCnt = _let->evaTrgNodeCnt();
00261          int evaTrgExaCnt = _let->evaTrgExaCnt();
00262          pC( VecCreateSeq(PETSC_COMM_SELF, evaTrgExaCnt*dim(),  &_evaTrgExaPos) );
00263          pC( VecCreateSeq(PETSC_COMM_SELF, evaTrgExaCnt*trgDOF(), &_evaTrgExaVal) );
00264          pC( VecCreateSeq(PETSC_COMM_SELF, evaTrgNodeCnt*datSze(DE), &_evaTrgDwnEquDen) );
00265          pC( VecCreateSeq(PETSC_COMM_SELF, evaTrgNodeCnt*datSze(DC), &_evaTrgDwnChkVal) );
00266          
00267          //2. gather data from _trgPos
00268          int procLclStart, procLclEnd; _let->procLclRan(_trgPos, procLclStart, procLclEnd);
00269          double* parr; pC( VecGetArray(_trgPos, &parr) );
00270          vector<int> ordvec; pC( _let->upwOrderCollect(ordvec) );
00271          for(int i=0; i<ordvec.size(); i++) {
00272                 int gNodeIdx = ordvec[i];
00273                 if(_let->node(gNodeIdx).tag() & LET_EVTRNODE) {
00274                   if(_let->terminal(gNodeIdx)==true) {
00275                          DblNumMat evaTrgExaPos(this->evaTrgExaPos(gNodeIdx));
00276                          vector<int>& curVecIdxs = _let->node(gNodeIdx).evaTrgOwnVecIdxs();
00277                          for(int k=0; k<curVecIdxs.size(); k++) {
00278                                 int poff = curVecIdxs[k]-procLclStart;
00279                                 for(int d=0; d<dim(); d++)
00280                                   evaTrgExaPos(d,k) = parr[poff*dim()+d];
00281                          }
00282                   }
00283                 }
00284          }
00285          pC( VecRestoreArray(_trgPos, &parr) );
00286          //3. allocate ENExt
00287          for(int i=0; i<ordvec.size(); i++) {
00288                 int gNodeIdx = ordvec[i];
00289                 if( _let->node(gNodeIdx).tag() & LET_EVTRNODE) {
00290                   node(gNodeIdx).vLstInNum() = _let->node(gNodeIdx).Vnodes().size();
00291                 }
00292          }
00293   }
00294     
00295   return(0);
00296 }
00297 
00298 

Generated on Sun Dec 4 21:12:40 2005 for fmm3d_mpi by  doxygen 1.4.5