omg.h

Go to the documentation of this file.

#ifndef __OCT_MG_H__
#define __OCT_MG_H__

#include "petscpc.h"
#include "petscksp.h"
#include <vector>

#ifdef PETSC_USE_LOG

#include "petscsys.h"

namespace ot {
  extern int pcKspShellSetupEvent;
  extern int pcKspShellDestroyEvent;
  extern int pcKspShellApplyEvent;
  extern int setDaEvent;
  extern int setUpEvent;
  extern int createRp1Event;
  extern int createRp2Event;
  extern int setKspEvent;
  extern int restrictEvent;
  extern int dummyRestrictEvent;
  extern int prolongEvent;
  extern int scatterEvent;
  extern int damgInitEvent;
  extern int damgFinalEvent;

  extern int setDAstage1Event;
  extern int setDAstage2Event;
  extern int setDAstage3Event;
  extern int setDAstage4Event;
  extern int setDAstage5Event;
  extern int setDAstage6Event;
}

#define PROF_SET_DA_STAGE1_BEGIN PetscLogEventBegin(setDAstage1Event,0,0,0,0);
#define PROF_SET_DA_STAGE1_END PetscLogEventEnd(setDAstage1Event,0,0,0,0); 

#define PROF_SET_DA_STAGE2_BEGIN PetscLogEventBegin(setDAstage2Event,0,0,0,0);
#define PROF_SET_DA_STAGE2_END PetscLogEventEnd(setDAstage2Event,0,0,0,0); 

#define PROF_SET_DA_STAGE3_BEGIN PetscLogEventBegin(setDAstage3Event,0,0,0,0);
#define PROF_SET_DA_STAGE3_END PetscLogEventEnd(setDAstage3Event,0,0,0,0); 

#define PROF_SET_DA_STAGE4_BEGIN PetscLogEventBegin(setDAstage4Event,0,0,0,0);
#define PROF_SET_DA_STAGE4_END PetscLogEventEnd(setDAstage4Event,0,0,0,0); 

#define PROF_SET_DA_STAGE5_BEGIN PetscLogEventBegin(setDAstage5Event,0,0,0,0);
#define PROF_SET_DA_STAGE5_END PetscLogEventEnd(setDAstage5Event,0,0,0,0); 

#define PROF_SET_DA_STAGE6_BEGIN PetscLogEventBegin(setDAstage6Event,0,0,0,0);
#define PROF_SET_DA_STAGE6_END PetscLogEventEnd(setDAstage6Event,0,0,0,0); 

#define PROF_PC_KSP_SHELL_SETUP_BEGIN  PetscFunctionBegin; \
  PetscLogEventBegin(pcKspShellSetupEvent,0,0,0,0);

#define PROF_PC_KSP_SHELL_SETUP_END PetscLogEventEnd(pcKspShellSetupEvent,0,0,0,0); \
  PetscFunctionReturn(0);

#define PROF_PC_KSP_SHELL_DESTROY_BEGIN  PetscFunctionBegin; \
  PetscLogEventBegin(pcKspShellDestroyEvent,0,0,0,0);

#define PROF_PC_KSP_SHELL_DESTROY_END PetscLogEventEnd(pcKspShellDestroyEvent,0,0,0,0); \
  PetscFunctionReturn(0);

#define PROF_PC_KSP_SHELL_APPLY_BEGIN  PetscFunctionBegin; \
  PetscLogEventBegin(pcKspShellApplyEvent,0,0,0,0);

#define PROF_PC_KSP_SHELL_APPLY_END PetscLogEventEnd(pcKspShellApplyEvent,0,0,0,0); \
  PetscFunctionReturn(0);

#define PROF_MG_INIT_BEGIN  PetscFunctionBegin; \
  PetscLogEventBegin(damgInitEvent,0,0,0,0);

#define PROF_MG_INIT_END PetscLogEventEnd(damgInitEvent,0,0,0,0); \
  PetscFunctionReturn(0);

#define PROF_MG_FINAL_BEGIN  PetscFunctionBegin; \
  PetscLogEventBegin(damgFinalEvent,0,0,0,0);

#define PROF_MG_FINAL_END PetscLogEventEnd(damgFinalEvent,0,0,0,0); \
  PetscFunctionReturn(0);

#define PROF_MG_SET_KSP_BEGIN  PetscFunctionBegin; \
  PetscLogEventBegin(setKspEvent,0,0,0,0);

#define PROF_MG_SET_KSP_END PetscLogEventEnd(setKspEvent,0,0,0,0); \
  PetscFunctionReturn(0);

#define PROF_MG_SET_DA_BEGIN  PetscFunctionBegin; \
  PetscLogEventBegin(setDaEvent,0,0,0,0);

#define PROF_MG_SET_DA_END PetscLogEventEnd(setDaEvent,0,0,0,0); \
  PetscFunctionReturn(0);

#define PROF_MG_SETUP_BEGIN  PetscFunctionBegin; \
  PetscLogEventBegin(setUpEvent,0,0,0,0);

#define PROF_MG_SETUP_END PetscLogEventEnd(setUpEvent,0,0,0,0); \
  PetscFunctionReturn(0);

#define PROF_MG_CREATE_RP1_BEGIN  PetscFunctionBegin; \
  PetscLogEventBegin(createRp1Event,0,0,0,0);

#define PROF_MG_CREATE_RP1_END PetscLogEventEnd(createRp1Event,0,0,0,0); \
  PetscFunctionReturn(0);

#define PROF_MG_CREATE_RP2_BEGIN  PetscFunctionBegin; \
  PetscLogEventBegin(createRp2Event,0,0,0,0);

#define PROF_MG_CREATE_RP2_END PetscLogEventEnd(createRp2Event,0,0,0,0); \
  PetscFunctionReturn(0);

#define PROF_MG_RESTRICT_BEGIN PetscFunctionBegin; \
  PetscLogEventBegin(restrictEvent,0,0,0,0);

#define PROF_MG_RESTRICT_END PetscLogEventEnd(restrictEvent,0,0,0,0); \
  PetscFunctionReturn(0);

#define PROF_MG_RESTRICT_DUMMY_BEGIN PetscFunctionBegin; \
  PetscLogEventBegin(dummyRestrictEvent,0,0,0,0);

#define PROF_MG_RESTRICT_DUMMY_END PetscLogEventEnd(dummyRestrictEvent,0,0,0,0); \
  PetscFunctionReturn(0);

#define PROF_MG_PROLONG_BEGIN  PetscFunctionBegin; \
  PetscLogEventBegin(prolongEvent,0,0,0,0);

#define PROF_MG_PROLONG_END PetscLogEventEnd(prolongEvent,0,0,0,0); \
  PetscFunctionReturn(0);

#define PROF_MG_SCATTER_BEGIN PetscLogEventBegin(scatterEvent,0,0,0,0);

#define PROF_MG_SCATTER_END PetscLogEventEnd(scatterEvent,0,0,0,0); \
  return 0;

#else

#define PROF_SET_DA_STAGE1_BEGIN 
#define PROF_SET_DA_STAGE1_END 

#define PROF_SET_DA_STAGE2_BEGIN 
#define PROF_SET_DA_STAGE2_END 

#define PROF_SET_DA_STAGE3_BEGIN 
#define PROF_SET_DA_STAGE3_END 

#define PROF_SET_DA_STAGE4_BEGIN 
#define PROF_SET_DA_STAGE4_END 

#define PROF_SET_DA_STAGE5_BEGIN 
#define PROF_SET_DA_STAGE5_END 

#define PROF_SET_DA_STAGE6_BEGIN 
#define PROF_SET_DA_STAGE6_END 

#define PROF_PC_KSP_SHELL_SETUP_BEGIN PetscFunctionBegin; 

#define PROF_PC_KSP_SHELL_SETUP_END PetscFunctionReturn(0);

#define PROF_PC_KSP_SHELL_DESTROY_BEGIN PetscFunctionBegin; 

#define PROF_PC_KSP_SHELL_DESTROY_END PetscFunctionReturn(0);

#define PROF_PC_KSP_SHELL_APPLY_BEGIN PetscFunctionBegin; 

#define PROF_PC_KSP_SHELL_APPLY_END PetscFunctionReturn(0);

#define PROF_MG_INIT_BEGIN PetscFunctionBegin; 

#define PROF_MG_INIT_END PetscFunctionReturn(0);

#define PROF_MG_FINAL_BEGIN PetscFunctionBegin; 

#define PROF_MG_FINAL_END PetscFunctionReturn(0);

#define PROF_MG_SET_KSP_BEGIN PetscFunctionBegin; 

#define PROF_MG_SET_KSP_END PetscFunctionReturn(0);

#define PROF_MG_SET_DA_BEGIN PetscFunctionBegin; 

#define PROF_MG_SET_DA_END PetscFunctionReturn(0);

#define PROF_MG_SETUP_BEGIN PetscFunctionBegin; 

#define PROF_MG_SETUP_END PetscFunctionReturn(0);

#define PROF_MG_CREATE_RP1_BEGIN PetscFunctionBegin; 

#define PROF_MG_CREATE_RP1_END PetscFunctionReturn(0);

#define PROF_MG_CREATE_RP2_BEGIN PetscFunctionBegin; 

#define PROF_MG_CREATE_RP2_END PetscFunctionReturn(0);

#define PROF_MG_SCATTER_BEGIN 

#define PROF_MG_SCATTER_END return 0;

#define PROF_MG_RESTRICT_BEGIN PetscFunctionBegin; 

#define PROF_MG_RESTRICT_END PetscFunctionReturn(0);

#define PROF_MG_RESTRICT_DUMMY_BEGIN PetscFunctionBegin; 

#define PROF_MG_RESTRICT_DUMMY_END PetscFunctionReturn(0);

#define PROF_MG_PROLONG_BEGIN PetscFunctionBegin; 

#define PROF_MG_PROLONG_END PetscFunctionReturn(0);

#endif

namespace ot {

  typedef struct {
    bool iAmActive; 
    MPI_Comm commActive; 
    PC pc; 
    KSP ksp_private; 
    Vec rhs_private; 
    Vec sol_private; 
  } PC_KSP_Shell;

  PetscErrorCode PC_KSP_Shell_SetUp(void* ctx);

  PetscErrorCode PC_KSP_Shell_Apply(void* ctx, Vec in, Vec out);

  PetscErrorCode PC_KSP_Shell_Destroy(void* ctx);

  class FineTouchedStatus {
    public:
      unsigned char flags[8];

      FineTouchedStatus() {
        for(int i=0;i<8;i++) {
          flags[i] = 0;
        }
      }    

      FineTouchedStatus(const FineTouchedStatus & other) {
        for(int i=0; i<8; i++) {
          this->flags[i] = other.flags[i];
        }
      }

      FineTouchedStatus & operator = (FineTouchedStatus const  & other) {
        if(this == (&other)) {return *this;}    
        for(int i=0; i<8; i++) {
          this->flags[i] = other.flags[i];
        }
        return *this;
      }//end fn.

      FineTouchedStatus & operator +=(FineTouchedStatus const & other) {
        //Self-Assignment not a problem!
        for(int i=0; i<8; i++) {
          this->flags[i] |= other.flags[i];
        }
        return *this;
      }//end fn.

      bool  operator == ( FineTouchedStatus const  &other) const {
        for(int i = 0; i < 8; i++ ) {
          if(this->flags[i] != other.flags[i]) {
            return false;
          }
        }
        return true;
      }

      bool  operator != ( FineTouchedStatus const  &other) const {
        return (!((*this) == other));
      }

      const FineTouchedStatus operator + (FineTouchedStatus const &other) const {
        return ((FineTouchedStatus(*this)) += other);
      }

  };

  class FineTouchedDummyStatus {
    public:
      unsigned char flags[16];

      FineTouchedDummyStatus() {
        for(int i=0;i<16;i++) {
          flags[i] = 0;
        }
      }    

      FineTouchedDummyStatus(const FineTouchedDummyStatus & other) {
        for(int i=0; i<16; i++) {
          this->flags[i] = other.flags[i];
        }
      }

      FineTouchedDummyStatus & operator = (FineTouchedDummyStatus const  & other) {
        if(this == (&other)) {return *this;}    
        for(int i=0; i<16; i++) {
          this->flags[i] = other.flags[i];
        }
        return *this;
      }//end fn.

      FineTouchedDummyStatus & operator +=(FineTouchedDummyStatus const & other) {
        //Self-Assignment not a problem!
        for(int i=0; i<16; i++) {
          this->flags[i] |= other.flags[i];
        }
        return *this;
      }//end fn.

      bool  operator == ( FineTouchedDummyStatus const  &other) const {
        for(int i = 0; i < 16; i++ ) {
          if(this->flags[i] != other.flags[i]) {
            return false;
          }
        }
        return true;
      }

      bool  operator != ( FineTouchedDummyStatus const  &other) const {
        return (!((*this) == other));
      }

      const FineTouchedDummyStatus operator + (FineTouchedDummyStatus const &other) const {
        return ((FineTouchedDummyStatus(*this)) += other);
      }

  };

  //Forward Declarations
  class DA;
  class TreeNode;

  typedef struct {
    ot::DA* dac; 
    ot::DA* daf; 
    unsigned int minIndependentSize; 
    unsigned char* suppressedDOFc; 
    unsigned char* suppressedDOFf; 
    std::vector<ot::FineTouchedStatus >* fineTouchedFlags; 
    unsigned int dof; 
    MPI_Comm comm;
    Vec tmp; //For R/P-type2 scatter
    Vec addRtmp;
    Vec addPtmp;
    int * sendSzP; 
    int * sendOffP; 
    int * recvSzP; 
    int * recvOffP;
    int * sendSzR; 
    int * sendOffR; 
    int * recvSzR; 
    int * recvOffR;
  } TransferOpData;

  struct _p_DAMG {

    ot::DA* da; 
    ot::DA* da_aux;  
    unsigned int dof; 
    unsigned char* suppressedDOF; 
    unsigned char* suppressedDOFaux; 
    Vec            x,b,r;                
    Mat            J;                    
    Mat            B;                   
    Mat            R;                    
    int       nlevels;              
    int       totalLevels; 
    MPI_Comm       comm; 
    PetscErrorCode (*solve)(_p_DAMG**, int);

    void           *user;        
    KSP            ksp;  
    PetscErrorCode (*initialguess)(_p_DAMG*, Vec); 
    PetscErrorCode (*rhs)(_p_DAMG*,Vec); 
  };//end struct definition

  typedef struct _p_DAMG* DAMG; 

  //Public Functions

  PetscErrorCode DAMGCreateAndSetDA(MPI_Comm comm, int& nlevels,
      void* user, DAMG** damg, std::vector<ot::TreeNode>& finestOctree,
      unsigned int dof = 1, double loadFac = 1.5,
      bool compressLut = true, bool incCorner = true);

  PetscErrorCode DAMGCreateJMatrix(DAMG damg, PetscErrorCode (*crJ)(DAMG, Mat*J));

  PetscErrorCode DAMGSetKSP(DAMG* damg, PetscErrorCode (*crjac)(DAMG, Mat* B),
      PetscErrorCode (*computeJac)(DAMG, Mat J, Mat B), PetscErrorCode (*rhs)(DAMG, Vec));

  PetscErrorCode DAMGSolve(DAMG* damg);

  PetscErrorCode DAMGSetInitialGuess(DAMG*, PetscErrorCode (*)(DAMG, Vec));

  PetscErrorCode DAMGInitialGuessCurrent(DAMG, Vec);

  void PrintDAMG(DAMG*);

  int DAMGCreateSuppressedDOFs(DAMG* damg);

  void DAMG_InitPrivateType1(MPI_Comm comm);

  void DAMG_InitPrivateType2(MPI_Comm comm);

  void DAMG_InitPrivateType3(MPI_Comm comm);

  PetscErrorCode DAMG_Initialize(MPI_Comm comm);

  PetscErrorCode DAMG_Finalize();

  PetscErrorCode DAMGDestroy(DAMG* damg);

#define DAMGGetRHS(ctx)              (ctx)[(ctx)[0]->nlevels-1]->b

#define DAMGGetr(ctx)              (ctx)[(ctx)[0]->nlevels-1]->r

#define DAMGGetx(ctx)              (ctx)[(ctx)[0]->nlevels-1]->x

#define DAMGGetComm(ctx)           (ctx)[(ctx)[0]->nlevels-1]->comm

#define DAMGGetJ(ctx)              (ctx)[(ctx)[0]->nlevels-1]->J

#define DAMGGetB(ctx)              (ctx)[(ctx)[0]->nlevels-1]->B

#define DAMGGetFine(ctx)           (ctx)[(ctx)[0]->nlevels-1]

#define DAMGGetKSP(ctx)            (ctx)[(ctx)[0]->nlevels-1]->ksp

#define DAMGGetDA(ctx)             ((ctx)[(ctx)[0]->nlevels-1]->da)

#define DAMGGetUser(ctx,level)     ((ctx)[level]->user)

#define DAMGSetUser(ctx,level,usr) ((ctx)[level]->user = usr,0)

#define DAMGGetLevels(ctx)         (ctx)[0]->nlevels

#define DAMGGetDAMG(ctx)              (ctx)[(ctx)[0]->nlevels-1]

  PetscErrorCode DAMGSetNullSpace(DAMG* damg, PetscTruth, int, PetscErrorCode (*)(DAMG,Vec[]));

  int scatterValues(Vec in, Vec out, PetscInt inSz, PetscInt outSz,
      int *& sendSz, int *& sendOff, int *& recvSz, int *& recvOff, MPI_Comm comm);

  PetscErrorCode DAMGSetUp(DAMG*);

  PetscErrorCode DAMGSolveKSP(DAMG *damg, int level);

  PetscErrorCode DAMGSetUpLevel(DAMG* damg, KSP ksp, int nlevels);

  /*Matrix-Free Intergrid Transfer Operators */
  int destroyRmatType1Stencil(double *****&lut);
  int destroyRmatType2Stencil(double ****&lut);
  int destroyVtxMaps(unsigned short ****&map1, unsigned short *****&map2,
      unsigned short *****&map3, unsigned short ******&map4);

  int readRmatType1Stencil(double *****&lut);
  int readRmatType2Stencil(double ****&lut);
  int readVtxMaps(unsigned short ****&map1, unsigned short *****&map2,
      unsigned short *****&map3, unsigned short ******&map4);

  int IreadRmatType1Stencil(double *****&lut, int rank);
  int IreadRmatType2Stencil(double ****&lut, int rank);
  int IreadVtxMaps(unsigned short ****&map1, unsigned short *****&map2,
      unsigned short *****&map3, unsigned short ******&map4, int rank);


  //Coarse and fine are aligned. No need to use da_aux
  PetscErrorCode createInterpolationType1(DAMG , DAMG, Mat *);
  PetscErrorCode prolongMatVecType1(Mat, Vec, Vec);
  PetscErrorCode restrictMatVecType1(Mat, Vec, Vec);
  PetscErrorCode dummyRestrictMatVecType1(TransferOpData* data);

  //Coarse and fine are NOT aligned. Need to use da_aux.
  PetscErrorCode  createInterpolationType2(DAMG , DAMG, Mat *);
  PetscErrorCode prolongMatVecType2(Mat, Vec, Vec);
  PetscErrorCode   restrictMatVecType2(Mat, Vec, Vec);

  PetscErrorCode addProlongMatVec(Mat, Vec, Vec, Vec);
  PetscErrorCode  addRestrictMatVec(Mat, Vec, Vec, Vec);

  PetscErrorCode  rpMatDestroy(Mat);


}//end namespace


namespace par {

  //Forward Declaration
  template <typename T>
    class Mpi_datatype;

  //Used only with ReadFromGhosts
  template <>
    class Mpi_datatype< ot::FineTouchedStatus > {
      public: 
        static MPI_Datatype value()
        {
          static bool  first = true;
          static MPI_Datatype datatype;

          if (first)
          {
            first = false;
            MPI_Type_contiguous(sizeof(ot::FineTouchedStatus), MPI_BYTE, &datatype);
            MPI_Type_commit(&datatype);
          }

          return datatype;
        }
    };

  //Used only with WriteToGhosts
  template <>
    class Mpi_datatype< ot::FineTouchedDummyStatus > {
      public: 
        static MPI_Datatype value()
        {
          static bool  first = true;
          static MPI_Datatype datatype;

          if (first)
          {
            first = false;
            MPI_Type_contiguous(sizeof(ot::FineTouchedDummyStatus), MPI_BYTE, &datatype);
            MPI_Type_commit(&datatype);
          }

          return datatype;
        }
    };

}//end namespace par


#endif

---