src/expde/grid/parallel.h

Go to the documentation of this file.

//    expde: expression templates for partial differential equations.
//    Copyright (C) 2001  Christoph Pflaum
//    This program is free software; you can redistribute it and/or modify
//    it under the terms of the GNU General Public License as published by
//    the Free Software Foundation; either version 2 of the License, or
//    (at your option) any later version.
//
//    This program is distributed in the hope that it will be useful,
//    but WITHOUT ANY WARRANTY; without even the implied warranty of
//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//    GNU General Public License for more details.
//
//                 SEE  Notice1.doc made by 
//                 LAWRENCE LIVERMORE NATIONAL LABORATORY
//

// ------------------------------------------------------------
// parallel.h
//
// ------------------------------------------------------------

#ifndef PARALLEL_H_
#define PARALLEL_H_
       
class Parallel_Info;

// Erklaerung:
// 3-dimensionale strukturierte Aufteilung der Prozessoren
// Hierzu gibt es verschiedene Moeglichkeiten; es wird eine optimale gesucht
// Eine Moeglichkeit der Aufteilung ist in einem Objekt vom Typ Partitioning 
// abgespeichert.
// Die unterschiedlichen Aufteilung werden durch
// a) durch unterschiedlich feine Aufteilungen der Tiefe p_level
//




//  Partitioning
// ---------------------------------

// hier gehoert noch folgendes rein:
// n_nec:   Anzahl der notwendigen Prozessoren
// p_level: Feinheit des Gitters zur Aufteilung der Prozessoren
// opt_I,opt_J,opt_K 

class Partitioning {
 private:
  bool *start;   // matrix of size^3=(2^p_level)^3
  int size;
 public:
  Partitioning(int level);
  ~Partitioning();
  void Put_proc(int i, int j, int k);
  bool Is_this_a_proc(Index3D Ind);
  void Print();
};

//  quality of decomposition
// ---------------------------------
class quality {
 public:
  quality(int n_p, int m_i, int a_i) : p_nec(n_p), minimal_intersec(m_i),
    maximal_intersec(a_i) {};
  quality() : p_nec(-1), minimal_intersec(0), maximal_intersec(0) {};
  bool operator<(quality other);
  bool positive_quality() { return p_nec>0; };
  int  Give_p_nec() { return p_nec; };
  void Print() {
    cout << " Quality: " << p_nec << " processors and itersec: "
         << minimal_intersec << endl;
  }
 private:
  int p_nec;            // number of necessary processors
  int minimal_intersec; // number of minimal intesection 
                        // for boundary processors
  int maximal_intersec; // number of minimal intesection 
};

inline bool quality::operator<(quality other) {
  /*
  if(other.p_nec>p_nec)  return true;
  if(other.p_nec==p_nec) return other.minimal_intersec>minimal_intersec;
  return false;
  */
  if(other.p_nec<=0) return false;
  if(p_nec<=0)       return true;
  if(other.maximal_intersec<maximal_intersec)  return true;
  if(other.maximal_intersec==maximal_intersec) 
    return other.minimal_intersec>minimal_intersec;
  return false;
}

//  rough domain
// ---------------------------------
//  rough_domain is constructed by looking at points
//  it stores all cells adjecent to an interior point to be
//  a interior cell, these cells must be covered by the processors
class rough_domain {
 private:
  int p_level;       // level procesors can be 0 <= p_level_min, ...
                     // n_parallel = p_level_min+1,...
  int refine_level;  // refine_level    can be 0,1, ...
  int n_interior;    // Zweipotenz(p_level+refine_level)
  int n_side;        // Zweipotenz(refine_level)
  bool   *start;     // rough domain informations
                     // size: n_interior-1 x n_interior-1 x n_interior-1
                     // entspricht: ]0,1[^3
  void set(int i, int j, int k);
  int l;
  int ll;
  // periodic
  bool is_periodic;

 public:
  rough_domain(int P_level, int Refine_level, All_Domains* dom,
               D3vector A_bounding, double H_bounding);
  ~rough_domain();

  int Give_n_side() { return n_side; }
  /*
  int p_necessary(int I, int J, int K,
                  int p);               // number of processors necessary for
                                        // I,J,K = 0,1, ..., n_side-1
  */
  quality calc_quality(int I, int J, int K,
                       int p);          // number of processors necessary for
                                        // I,J,K = 0,1, ..., n_side-1
  int Calc_optimal_partitioning(int p,  // p number of processors
                                int& opt_I, int& opt_J, int& opt_K,
                                Partitioning* parti);
  void Set_partitioning(int opt_I,int opt_J,int opt_K,Partitioning* parti);
  bool In_domain(int i, int j, int k);               // cell indomain ??
  int Give_n_interior() { return n_interior; };
};



class Grid_base;

//  Fuer Hashtable proc: Processes
// ---------------------------------
class Point_hashtable_proc {
  friend class Parallel_Info;
  friend class Grid_base;
 private:
  Index3D ind;      // midpoint of processes
  int num_proc;     // number   of process

  Point_hashtable_proc* next;
 public:
  bool used_on_coarser_grid;
  bool interior_processor;    // braucht man das ??????

  //  I Index des  Mittelpunktes der Zelle oder Kante
  Point_hashtable_proc(Index3D I) : ind(I) { 
    num_proc = -1;
    next = NULL;
    interior_processor   = true;
    used_on_coarser_grid = false;
  };
  Index3D Give_Index()    { return ind; }
  int     Give_num_proc() { return num_proc; }
  void    Set_num_proc(int p) {  num_proc=p; }

 ~Point_hashtable_proc() {
    if(next!=NULL) delete next;
  };
};

inline int hashtable_proc_function(int x, int y, int z, int length) {
  return (x + 101*y + 10007*z) % length;
}

//  Fuer Iteration durch Hashtable proc:
// -----------------------------------------
// folgendes define verwendet : point_proc:
#define iterate_hash_proc  for(p_iter=0;p_iter<hashtable_proc_leng;++p_iter) for(point_proc=hashtable_proc_start[p_iter];point_proc!=NULL;point_proc=point_proc->next) 

//  Fuer Hashtable proc_corner:
// ---------------------------------
class Point_hashtable_proc_corner {
  friend class Parallel_Info;
  friend class Grid_base;
 private:
  Index3D ind;      // corner point 
  int num_proc;     // number of responsible process

  Point_hashtable_proc_corner* next;
 public:
  int maximal_level; // with respect to points

  //  I Index des  Mittelpunktes der Zelle oder Kante
  Point_hashtable_proc_corner(Index3D I) : ind(I) { 
    num_proc = -1;
    next = NULL;
    maximal_level = 0;
  };
  Index3D Give_Index()    { return ind; }
  int     Give_num_proc() { return num_proc; }
  void    Set_num_proc(int p) {  num_proc=p; }

 ~Point_hashtable_proc_corner() {
    if(next!=NULL) delete next;
  };
};

inline int hashtable_proc_corner_function(int x, int y, int z, int length) {
  return (x + 101*y + 10007*z) % length;
}

//  Fuer Iteration durch Hashtable proc_corner: 
// -----------------------------------------
// folgendes define verwendet : point_proc:
#define iterate_hash_proc_corner  for(p_iter_corner=0;p_iter_corner<hashtable_proc_corner_leng;++p_iter_corner) for(point_proc_corner=hashtable_proc_corner_start[p_iter_corner];point_proc_corner!=NULL;point_proc_corner=point_proc_corner->next) 








//  Class describing processes
// ----------------------------

//class Parallel_Info : public Grid_gen_parameters {
class Parallel_Info  {
 public:
  // function for testing pre_Grid generation
  Parallel_Info(int n_max, All_Domains* dom, int number_processes, 
                Grid_gen_parameters& gpara);
  Parallel_Info(double meshsize, All_Domains* dom, int number_processes, 
                Grid_gen_parameters& gpara);

  int Give_my_rank() const { return my_rank; };
  Index3D Give_my_index() const { return my_index; }
  Index3D Give_my_level_index(int level); // level here with respect to 
                                          // cellpoint
  Index3D Give_level_index_of(int rang, int level);
  int Give_coarsest_level_of(int rang);
   
  int Give_n_parallel() const { return n_parallel; } 
  int Give_my_coarsest_level()            // coarsest level for this process
    { return my_coarsest_level; }         // with respect to points

  bool exists_process(Index3D I);
  bool I_am_active();

  inline int give_number_of_processes();
  inline int give_number_of_active_processes();
  void process_info();

  void Print_processors();
  void Print_region_processes_UCD_normalized(ofstream *Datei);

  MPI_Comm Give_comm() { return comm; };
  bool        I_am_responsible_for(Index3D I) { return I<my_index; }
  inline bool I_am_responsible_for(Index3D I, int level); // level w.r. points
  bool        I_am_responsible_for(IndexBo Ib);

  inline All_Domains* Give_domain() { return domain; }
  
  ~Parallel_Info();

 protected:
  Parallel_Info(All_Domains* dom, MPI_Comm comm);

  // ^^^^^^^^^^^^^^^^^^^^  Part A: Processes ^^^^^^^^^^^^^^^^^^^
  void Generate_processes();
  void Generate_seriel_process();

  // informations for all processes  ///////////////////////////////
  All_Domains* domain;  // domain
  MPI_Comm comm;
  int p;          // number of processes
  int n_parallel; 

  void Add_proc_cell(Index3D I);               // Add cell for process
  void Add_proc_corner(Index3D I, int l);      // Add corner process

  // transform computational coordinate in physical coordinate
  inline D3vector coord_in_domain(const Index3D  I) const;

  // for bounding box
  double   H_bounding;   
  D3vector A_bounding;

  // for making grid through one grid point
  void make_grid_with_grid_point(D3vector gp, double fms);

  // for old bounding box and local box
  bool contained_in_old_bounding_box(Index3D I) const;
  bool contained_in_boxes(Index3D I) const;  // local and old bounding

  // maximal and minimal level
  int max_level;
  int min_level;

  // informations for all processes  ///////////////////////////////
  int my_rank;             // my rank of process
  Index3D my_index;        // my index
  int my_coarsest_level;   // coarsest level for this process
                           // with respect to points
  Index3D *my_level_index; // my index in coarser grids than n_parallel
  Index3D Give_next_on_level(int d,          // 0 <= d < 26
                             Index3D me);    // next proc index
                                             // if it exists or not  

  // void Set_coarser_index_and_level();
  void Calc_coarse_processors();
  void Recursion_calc_coarse_processors(Index3D I);
  

  // ^^^^^^^^^^^^^^^^^^^^  Part  B:  Communication ^^^^^^^^^^^^^^^^
  // buffers used for several communications
  int* Give_buffer_send_info(int l);
  int* Give_buffer_receive_info(int l);

  // communication general in direction d (dir_3D)
  void Set_rank_3D();
  int give_next_rank_source(dir_3D d,int level);       // -1 if not exists
  int give_next_rank_destination(dir_3D d,int level);  // -1 if not exists

  // communication general in edge direction d (Edges_cell)
  int give_next_rank_source(Edges_cell ed);       // -1 if not exists
  int give_next_rank_destination(Edges_cell ed);  // -1 if not exists
  int give_next_rank_source(Edges_cell ed, int level);      // -1 if not exists
  int give_next_rank_destination(Edges_cell ed, int level); // -1 if not exists

  // communication general in corner direction d (dir_sons)
  int give_next_rank_source(dir_sons ed);       // -1 if not exists
  int give_next_rank_destination(dir_sons ed);  // -1 if not exists
  int give_next_rank_source(dir_sons ed, int level);       // -1 if not exists
  int give_next_rank_destination(dir_sons ed, int level);  // -1 if not exists

  // for communication during expression template evaluation
  // -------------------------------------------------------
  // I. storage for faces:
  // a) number of points for each face
  int** number_receive_face;
  int** number_send_face;
     // number of points for each face for prolongation
  int** number_receive_face_prolongation;
  int** number_send_face_prolongation;
  // b) where the start of "var" is
  double**** var_receive;
  double**** var_send;
     // where the start of "var" for prolongation is
  double**** var_receive_prolongation;
  double**** var_send_prolongation;
  // c) buffer
  double* Give_buffer_send(int d, int l); 
  double* Give_buffer_receive(int d, int l);


  // Explaination:
  // double** xxxx[num][level][direction];

  // interior points at faces: contained in var_receive[...]
  //                                 and in    var_send[...]
  // Zellfreiheitsgrade and Randpunkte:
  //                           contained in var_receive[E]  [N]  [T]
  //                                 and in    var_send[E]  [N]  [T]
  // II. storage for boundary stencil:
  // a) number of points for each direction
  int number_receive_bo_stencil[Max_number_levels][26];
  int number_send_bo_stencil[Max_number_levels][26];
  // b) where the start of "stencil" is
  double** bo_stencil_receive[Max_number_levels][26];
  double** bo_stencil_send[Max_number_levels][26];
  // c) buffer
  // same as above

  // II. storage for stencil:
  // a) number of points for each direction
  int number_receive_stencil[Max_number_levels][26];
  int number_send_stencil[Max_number_levels][26];
  // b) where the start of "stencil" is
  double** stencil_receive[Max_number_levels][26];
  double** stencil_send[Max_number_levels][26];
  // c) buffer
  // same as above




// private:
  void Recursion_study_cells(Index3D I, Partitioning* parti);
  // for partitioning
  int refine_level;       // 0,1,2,... optimzation parameter for partitioning
  int p_level;            // level procesors can be 0 <= p_level_min, ...
                          // n_parallel = p_level_min+1,...
  int number_of_active_procs; // number of active processors

  // for old bounding box
  Index1D Ixl, Ixr;
  Index1D Iyl, Iyr;
  Index1D Izl, Izr;

  double  pxl, pxr;
  double  pyl, pyr;
  double  pzl, pzr;

  // values for a local box
  IndexSurface coarse_faces[6];
  bool         info_coarse_faces[6];
  IndexEdge    coarse_edges[12];
  bool         info_coarse_edges[12];
  Index3D      coarse_corners[8];
  bool         info_coarse_corners[8];

  // buffers used for several communications
  int* buffer_send_info;
  int* buffer_receive_info;
  int length_send_info;
  int length_receive_info;

  // for communication
  // ------------------
  // next processor for level >= n_parallel-1
  int next_rank_source[6];
  int next_rank_destination[6];

  int next_rank_source_edges[12];
  int next_rank_destination_edges[12];

  int next_rank_source_corners[8];
  int next_rank_destination_corners[8];

  // next processor for level < n_parallel-1
  int* next_rank_source_coarse[6];
  int* next_rank_destination_coarse[6];

  int* next_rank_source_coarse_edges[12];
  int* next_rank_destination_coarse_edges[12];

  int* next_rank_source_coarse_corners[8];
  int* next_rank_destination_coarse_corners[8];

  // buffers ...
  double* buffer_send[26];
  double* buffer_receive[26];
  int length_send[26];
  int length_receive[26];

  // Hashtable
  // for iteration
  int p_iter;
  int p_iter_corner;
  // Funktionen zur Anwendung:
  int  Give_proc_number(Index3D I);
  bool Exists_proc(Index3D I);
  int  Give_proc_corner_number(Index3D I);
  bool Exists_proc_corner(Index3D I);


  Point_hashtable_proc* hashtable_proc_point(Index3D I) const;
  Point_hashtable_proc_corner* hashtable_proc_corner_point(Index3D I) const;

  void Initialize_hash_proc(int lenght);
  void Initialize_hash_proc_corner(int lenght);

  // Speicher des hashtable_proc:
  Point_hashtable_proc* point_proc;          //fuer Iterationen
  int hashtable_proc_leng;               // Laenge des hashtable
  int hashtable_proc_occ;                // Besetzungszahl des hashtable
  Point_hashtable_proc** hashtable_proc_start;  //Zeiger auf Tabelle

  // Speicher des hashtable_proc_corner:
  Point_hashtable_proc_corner* point_proc_corner;  //fuer Iterationen
  int hashtable_proc_corner_leng;                  // Laenge des hashtable
  int hashtable_proc_corner_occ;                // Besetzungszahl des hashtable
  Point_hashtable_proc_corner** hashtable_proc_corner_start; 
                                                //Zeiger auf Tabelle
};



// Implementierung einiger Memberfunktionen

// transform computational coordinate in physical coordinate
inline D3vector Parallel_Info::coord_in_domain(const Index3D  I) const {
  D3vector V;
  V = I.coordinate();
  return D3vector(A_bounding.x + V.x*H_bounding,
                  A_bounding.y + V.y*H_bounding,
                  A_bounding.z + V.z*H_bounding);
};

inline int Parallel_Info::give_number_of_processes() {
  return p;
};

inline int Parallel_Info::give_number_of_active_processes() {
  return number_of_active_procs;
}

inline bool Parallel_Info::I_am_active() {
  if(parallel_version)
    return my_rank<number_of_active_procs;
  else 
    return true;
}

// hash proc

inline Point_hashtable_proc* Parallel_Info::hashtable_proc_point(Index3D I)
     const {
  Point_hashtable_proc* point;
  point = hashtable_proc_start
    [hashtable_proc_function(I.ind_x.index,I.ind_y.index,I.ind_z.index,
                         hashtable_proc_leng)];
  while(point!=NULL) {
    if(point->ind==I) return point;
    point = point->next;
  }
  return NULL;
}

inline Point_hashtable_proc_corner* 
Parallel_Info::hashtable_proc_corner_point(Index3D I)
     const {
  static Point_hashtable_proc_corner* point;
  point = hashtable_proc_corner_start
    [hashtable_proc_corner_function(I.ind_x.index,I.ind_y.index,I.ind_z.index,
                                    hashtable_proc_corner_leng)];
  while(point!=NULL) {
    if(point->ind==I) return point;
    point = point->next;
  }
  return NULL;
}

inline bool Parallel_Info::exists_process(Index3D I) {
  return hashtable_proc_point(I)!=NULL;
}


inline int Parallel_Info::give_next_rank_source(dir_3D d,int level) {
  if(level>=n_parallel-1)
    return next_rank_source[d];
  return next_rank_source_coarse[d][level];
}

inline int Parallel_Info::give_next_rank_destination(dir_3D d,int level) {
  if(level>=n_parallel-1)
    return next_rank_destination[d];
  return next_rank_destination_coarse[d][level];
}

inline int Parallel_Info::give_next_rank_source(Edges_cell d,
                                                int level) {
  if(level>=n_parallel-1)
    return next_rank_source_edges[d];
  return next_rank_source_coarse_edges[d][level];
}

inline int Parallel_Info::give_next_rank_destination(Edges_cell d,
                                                     int level) {
  if(level>=n_parallel-1)
    return next_rank_destination_edges[d];
  return next_rank_destination_coarse_edges[d][level];
}

inline int Parallel_Info::give_next_rank_source(dir_sons d,
                                                int level) {
  if(level>=n_parallel-1)
    return next_rank_source_corners[d];
  return next_rank_source_coarse_corners[d][level];
}

inline int Parallel_Info::give_next_rank_destination(dir_sons d,
                                                     int level) {
  if(level>=n_parallel-1)
    return next_rank_destination_corners[d];
  return next_rank_destination_coarse_corners[d][level];
}



inline bool Parallel_Info::I_am_responsible_for(IndexBo Ib) {
  return (Ib.I.next(Ib.d,MAX(n_parallel,Ib.I.Tiefe())+1)) < my_index;
}

inline bool Parallel_Info::I_am_responsible_for(Index3D I, int level) {
  // level w.r. points
  return I < Give_my_level_index(level+1);
}

#endif
     

---