src/expde/evpar/commbos.cc

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
//

// ------------------------------------------------------------
//
// communi.cc
//
// ------------------------------------------------------------


// Include level 0:
#ifdef COMP_GNUOLD
 #include <iostream.h>
 #include <fstream.h>
 #include <time.h>
 #include <math.h>
#else
 #include <iostream>
 #include <fstream>
 #include <ctime>
 #include <cmath>
#endif 



// Include level 0:
#include "../parser.h"

// Include level 1:
#include "../paramete.h"
#include "../abbrevi.h"
#include "../math_lib/math_lib.h"

// Include level 2:
#include "../basic/basic.h"

// Include level 3:
#include "../domain/domain.h"

// Include level 4:
#include "../formulas/boundy.h"
#include "../formulas/loc_sten.h"

// Include level 5
#include "../grid/gpar.h"
#include "../grid/parallel.h"
#include "../grid/mgcoeff.h"
#include "../grid/sto_man.h"
#include "../grid/gridbase.h"
#include "../grid/grid.h"

// 1. communication of boundary stencils
// 2. communication of stencils




// 1. communication of boundary stencils


bool Grid_base::Send_boundary_stencils_in_direction(Point_hashtable0* poi, 
                                                    int level,
                                                    Index3D my_lev_index,
                                                    Index3D next_index) {
  // here level is the level of the cell point  
  Index3D I;
  int co;

  /* old  
  if(poi->isCell()) { 
    if(poi->Give_Tiefe() == level) {
      if(I < my_lev_index) {
        if(poi->typ==bo_cell) {
  */

  I = poi->ind;
  // send one layer in right and left direction
  for(co=0;co<8;++co) {
    // study corners of cell
    if(I.neighbour((dir_sons)co) << next_index)
      return true;
  }

  /* old
    }
    }
    }
    }
  */

  return false;
}



void Grid_base::Prepare_communication_boundary_stencils() {
  int l,i;

  // -2. set zero 
  for(i=0;i<26;++i)
    for(l=0;l<Max_number_levels;++l) {
      number_send_bo_stencil[l][i] = 0;
      number_receive_bo_stencil[l][i] = 0;
    }

  // -1. iterate levels
  for(l=min_level+1;l<=max_level;++l)
  //for(l=1;l<=max_level+1;++l) 
    Prepare_communication_boundary_stencils(l);
}

void Grid_base::Prepare_communication_boundary_stencils(int level) {
  // here level is the level of the cell point
  int d;
  Index3D I;
  
  int number_bo_send, number_bo_rec, num_message, num;
  Index3D my_lev_index;
  Index3D next_index;

  int* send_info;
  int* receive_info;
  //  int counter_bo_stencil[Max_number_levels];
  int counter_bo_stencil;
  double* sp;
  Point_hashtable0* poi0;
  
  MPI_Request req[2];
  MPI_Status status[2];
  
  int rank_source, rank_destination;
  

  // -1. set my_lev_index 
  my_lev_index = Give_my_level_index(level);


  // construct extra array for invariant points which are needed.
  int lenght_invar_array = 0;
  iterate_hash0 {
    if(point0->isCell()) { 
      if(point0->Give_Tiefe() == level) {
        if(point0->ind < my_lev_index) {
          if(point0->typ==bo_cell) {
            lenght_invar_array++;
          }
        }
      }
    }
  }
  Point_hashtable0** invariant_array = new Point_hashtable0*[lenght_invar_array];
  int array_i = 0;
  iterate_hash0 {
    if(point0->isCell()) { 
      if(point0->Give_Tiefe() == level) {
        if(point0->ind < my_lev_index) {
          if(point0->typ==bo_cell) {
            invariant_array[array_i] = point0;
            array_i++;
          }
        }
      }
    }
  }

  for(d=0;d<26;++d) { 
    // 0. index of receiving next processor 
    next_index = Give_next_on_level(d,my_lev_index);

    // 1. part: basic infos about send and recieve
    if(d<6) {
      rank_source      = give_next_rank_source((dir_3D)d,level-1);
      rank_destination = give_next_rank_destination((dir_3D)d,level-1);  
    }
    else if(d<18) {
      rank_source      = give_next_rank_source((Edges_cell)(d-6),level-1);
      rank_destination = give_next_rank_destination((Edges_cell)(d-6),level-1);  
    }
    else {
      rank_source      = give_next_rank_source((dir_sons)(d-18),level-1);
      rank_destination = give_next_rank_destination((dir_sons)(d-18),level-1);
    }

    
    // 2. part: count number of boundary cells to be send
    //          and count numbers with respect to levels
    number_bo_send = 0;

    for(int i = 0;i < lenght_invar_array;i++) {
      if(Send_boundary_stencils_in_direction(invariant_array[i],level,
                                             my_lev_index,next_index)) {
        
        number_send_bo_stencil[level][d]++;
        ++number_bo_send;       
      }
    }
    
    // 3. construct send double array 
    if(number_send_bo_stencil[level][d]!=0)
      bo_stencil_send[level][d] = new double*[number_send_bo_stencil[level][d]];
    else bo_stencil_send[level][d] = NULL;
    
    // 4. part: send informations about number_bo_send
    number_bo_rec = 0;
    num_message = 0;
    if(rank_source != -1) {
      MPI_Irecv(&number_bo_rec,1,MPI_INT,rank_source,1,comm,
                &req[num_message]);
      ++num_message;
    }
    if(rank_destination != -1) {
      MPI_Isend(&number_bo_send,1,MPI_INT,rank_destination,1,comm,
                &req[num_message]);
      ++num_message;
    }
    MPI_Waitall(num_message,req,status);
    
    // 8. construct double array 
    number_receive_bo_stencil[level][d] = number_bo_rec;
    if(number_receive_bo_stencil[level][d]!=0)
      bo_stencil_receive[level][d] = 
        new double*[number_receive_bo_stencil[level][d]];
    else bo_stencil_receive[level][d] = NULL;
    
    // 5. part: build informations about boundary cells
    //          and set pointers of boundary stencils of send
    // How informations are stored:
    // ind_x, ind_y, ind_z,    ...              ,   ...    >   3 int
    // 5.1. set counter zero
    counter_bo_stencil = 0;
    
    // 5.2. built and set
    send_info = Give_buffer_send_info(3*number_bo_send);
    receive_info = Give_buffer_receive_info(3*number_bo_rec);
    
    num = 0;
    for(int i = 0;i < lenght_invar_array;i++) {
      point0 = invariant_array[i];
      if(Send_boundary_stencils_in_direction(point0,level,
                                             my_lev_index,next_index)) {

        // cell index
        I = point0->ind;
        // built
        send_info[num*3]   = I.I_x().get();
        send_info[num*3+1] = I.I_y().get();
        send_info[num*3+2] = I.I_z().get();
        num++;
        // set
        bo_stencil_send[level][d][counter_bo_stencil] =
          Give_pointer_stencil(I);
        counter_bo_stencil++;
      }
    }
    
    // 6. part: send cell index
    num_message = 0;
    if(rank_source != -1 && number_bo_rec>0) {
      MPI_Irecv(receive_info,3*number_bo_rec,MPI_INT,rank_source,2,comm,
                &req[num_message]);
      ++num_message;
    }
    if(rank_destination != -1 && number_bo_send>0) {
      MPI_Isend(send_info,3*number_bo_send,MPI_INT,rank_destination,2,comm,
                &req[num_message]);
      ++num_message;
    }
    // MPI_Barrier(comm);
    MPI_Waitall(num_message,req,status);
    
    // 9. part: set pointers of boundary stencils of receive
    // 9.1. set counter zero
    counter_bo_stencil = 0;
    
    // 9.2. set
    for(num=0;num<number_bo_rec;++num) {
      I = Index3D(receive_info[num*3],
                  receive_info[num*3+1],
                  receive_info[num*3+2]);
      
      if(developer_version) {
        if(I.Cell_index()==false) 
          cout << "\n Mistake I in Prepare_communication_boundary_stencils() "
               << endl;
        if(Exists_Cell(I)==false) {
          cout << "\n Mistake II in Prepare_communication_boundary_stencils() "
               << endl;
        }
        if(I.Tiefe()!=level) {
          cout << "\n Mistake III in Prepare_communication_boundary_stencils() "
               << endl;
        }
      }
      
      bo_stencil_receive[level][d][counter_bo_stencil] =
        Give_pointer_stencil(I);
      if(developer_version) {
        if(bo_stencil_receive[level][d][counter_bo_stencil]==NULL
           && Give_max_num_bo_cells()>0) {
          //       && Give_cell_typ(I)==bo_cell) 
          cout << "\n Mistake IIV in Prepare_communication_boundary_stencils() "
               << Give_cell_typ(I) 
               << " level: " << I.Tiefe()
               << " x: " << I.I_x().get()
               << " y: " << I.I_y().get()
               << " z: " << I.I_z().get()  
               << " my_rank: " << my_rank << endl;
          I.coordinate().Print();
          cout << endl;
        }
      }
      
      if(bo_stencil_receive[level][d][counter_bo_stencil]==NULL
         //      && Give_cell_typ(I)==ex_cell) {
         && ((I < my_index)==false)) {
        sp = newdouble(Give_max_num_bo_cells());
        poi0 = hashtable0_point(I);
        poi0->var = sp;
        poi0->typ = bo_cell;
        bo_stencil_receive[level][d][counter_bo_stencil] = sp;
      }
      
      counter_bo_stencil++;
    }
  }
  delete invariant_array;
}

void Grid::Update_boundary_stencil(int num_stencil, int level) {
  int num, iv;
  int d;
  MPI_Request req[12];
  MPI_Status status[12];
  int num_message;
  double*  send_info;
  double*  receive_info;
  int number_send;
  int number_receive;
  int number_start;


  //  int number_bo_stencils;
  double** bo_stencil_sendd;
  double** bo_stencil_received;
  
  int rank_source, rank_destination;
  

  if(I_am_active()) {
    // 0. Step: number start:
    number_start = cell_number_of_bo_stencil(num_stencil);

    for(d=0;d<26;++d) { 
      // 1. part: basic infos about send and recieve
      if(d<6) {
        rank_source      = give_next_rank_source((dir_3D)d,level-1);
        rank_destination = give_next_rank_destination((dir_3D)d,level-1);  
      }
      else if(d<18) {
        rank_source      = give_next_rank_source((Edges_cell)(d-6),level-1);
        rank_destination = give_next_rank_destination((Edges_cell)(d-6),level-1);  
      }
      else {
        rank_source      = give_next_rank_source((dir_sons)(d-18),level-1);
        rank_destination = give_next_rank_destination((dir_sons)(d-18),level-1);
      }
      
      // 2. Step: calc number to be send and receive
      number_send = 
        number_send_bo_stencil[level][d] * 64;
      number_receive = 
        number_receive_bo_stencil[level][d] * 64;
      
      // 3. Step: build informations about bo_stencil(double) to be send  
      send_info    = Give_buffer_send(d,number_send); 
      receive_info = Give_buffer_receive(d,number_receive); 
      
      bo_stencil_sendd=bo_stencil_send[level][d];
      for(num=0;num<number_send_bo_stencil[level][d];++num) {
        for(iv=0;iv<64;++iv) {
          send_info[num*64+iv] = 
            bo_stencil_sendd[num][number_start+iv];
        }
      }
      
      // 4. Step: send doubles
      num_message = 0;
      if(rank_source != -1 && number_receive>0) {
        MPI_Irecv(receive_info ,number_receive,
                  MPI_DOUBLE,rank_source,10,comm,
                  &req[num_message]);
        ++num_message;
      }      
      if(rank_destination != -1 && number_send>0) {
        MPI_Isend(send_info,number_send,
                  MPI_DOUBLE,rank_destination,10,comm,
                  &req[num_message]);
        ++num_message;
      }
      MPI_Waitall(num_message,req,status);
      
      // 5. part: set informations about bo_stencil(double)
      bo_stencil_received=bo_stencil_receive[level][d];
      for(num=0;num<number_receive_bo_stencil[level][d];++num) {
        
        //  cout << " num: " << num << endl;
        
        for(iv=0;iv<64;++iv) {
          bo_stencil_received[num][number_start+iv] = 
            receive_info[num*64+iv];
        }
      }
    }
  }
}


// 2. communication of interior stencils

bool Grid_base::Send_interior_stencils_in_direction(Point_hashtable0* poi, 
                                           int level,
                                           Index3D my_lev_index,
                                           Index3D next_index) {
  // here level is the level of the cell point  
  Index3D I;
  int co;
  
  I = poi->ind;
  if(poi->isCell()) {
    if(poi->Give_Tiefe() == level) {
      if(I < my_lev_index) {
        if(poi->typ==int_cell) {
          // send one layer in right and left direction
          for(co=0;co<8;++co) {
            // study corners of cell
            if(I.neighbour((dir_sons)co) << next_index)
              return true;
          }
        }
      }
    }
  }
  return false;
}



void Grid_base::Prepare_communication_interior_stencils() {
  int l,i;

  // -2. set zero 
  for(i=0;i<26;++i)
    for(l=0;l<Max_number_levels;++l) {
      number_send_stencil[l][i] = 0;
      number_receive_stencil[l][i] = 0;
    }

  // -1. iterate levels
  if(Give_max_num_coarse_int_cells()>0)
    for(l=min_level+1;l<=max_level;++l)
      //for(l=1;l<=max_level+1;++l) 
      Prepare_communication_interior_stencils(l);
}

void Grid_base::Prepare_communication_interior_stencils(int level) {
  // here level is the level of the cell point
  int d;
  Index3D I;
  
  int number_send, number_rec, num_message, num;
  Index3D my_lev_index;
  Index3D next_index;

  int* send_info;
  int* receive_info;
  //  int counter_stencil[Max_number_levels];
  int counter_stencil;
  double* sp;
  Point_hashtable0* poi0;
  
  MPI_Request req[2];
  MPI_Status status[2];
  
  int rank_source, rank_destination;
  

  // -1. set my_lev_index 
  my_lev_index = Give_my_level_index(level);

  for(d=0;d<26;++d) { 
    // 0. index of receiving next processor 
    next_index = Give_next_on_level(d,my_lev_index);

    // 1. part: basic infos about send and recieve
    if(d<6) {
      rank_source      = give_next_rank_source((dir_3D)d,level-1);
      rank_destination = give_next_rank_destination((dir_3D)d,level-1);  
    }
    else if(d<18) {
      rank_source      = give_next_rank_source((Edges_cell)(d-6),level-1);
      rank_destination = give_next_rank_destination((Edges_cell)(d-6),level-1);  
    }
    else {
      rank_source      = give_next_rank_source((dir_sons)(d-18),level-1);
      rank_destination = give_next_rank_destination((dir_sons)(d-18),level-1);
    }

    
    // 2. part: count number of boundary cells to be send
    //          and count numbers with respect to levels
    number_send = 0;
    iterate_hash0 {
      if(Send_interior_stencils_in_direction(point0,level,
                                             my_lev_index,next_index)) {
        
        number_send_stencil[level][d]++;
        ++number_send;  
      }
    }
    
    // 3. construct send double array 
    if(number_send_stencil[level][d]!=0)
      stencil_send[level][d] = new double*[number_send_stencil[level][d]];
    else stencil_send[level][d] = NULL;
    
    // 4. part: send informations about number_send
    number_rec = 0;
    num_message = 0;
    if(rank_source != -1) {
      MPI_Irecv(&number_rec,1,MPI_INT,rank_source,1,comm,
                &req[num_message]);
      ++num_message;
    }
    if(rank_destination != -1) {
      MPI_Isend(&number_send,1,MPI_INT,rank_destination,1,comm,
                &req[num_message]);
      ++num_message;
    }
    MPI_Waitall(num_message,req,status);
    
    // 8. construct double array 
    number_receive_stencil[level][d] = number_rec;
    if(number_receive_stencil[level][d]!=0)
      stencil_receive[level][d] = 
        new double*[number_receive_stencil[level][d]];
    else stencil_receive[level][d] = NULL;
    
    // 5. part: build informations about boundary cells
    //          and set pointers of boundary interior_stencils of send
    // How informations are stored:
    // ind_x, ind_y, ind_z,    ...              ,   ...    >   3 int
    // 5.1. set counter zero
    counter_stencil = 0;
    
    // 5.2. built and set
    send_info = Give_buffer_send_info(3*number_send);
    receive_info = Give_buffer_receive_info(3*number_rec);
    
    num = 0;
    iterate_hash0 {
      if(Send_interior_stencils_in_direction(point0,level,
                                             my_lev_index,next_index)) {

        // cell index
        I = point0->ind;
        // built
        send_info[num*3]   = I.I_x().get();
        send_info[num*3+1] = I.I_y().get();
        send_info[num*3+2] = I.I_z().get();
        num++;
        // set
        stencil_send[level][d][counter_stencil] =
          Give_pointer_stencil(I);
        counter_stencil++;
      }
    }
    
    // 6. part: send cell index
    num_message = 0;
    if(rank_source != -1 && number_rec>0) {
      MPI_Irecv(receive_info,3*number_rec,MPI_INT,rank_source,2,comm,
                &req[num_message]);
      ++num_message;
    }
    if(rank_destination != -1 && number_send>0) {
      MPI_Isend(send_info,3*number_send,MPI_INT,rank_destination,2,comm,
                &req[num_message]);
      ++num_message;
    }
    // MPI_Barrier(comm);
    MPI_Waitall(num_message,req,status);
    
    // 9. part: set pointers of boundary interior_stencils of receive
    // 9.1. set counter zero
    counter_stencil = 0;
    
    // 9.2. set
    for(num=0;num<number_rec;++num) {
      I = Index3D(receive_info[num*3],
                  receive_info[num*3+1],
                  receive_info[num*3+2]);
      
      if(developer_version) {
        if(I.Cell_index()==false) 
          cout << "\n Mistake I in Prepare_communication_interior_stencils() "
               << endl;
        if(Exists_Cell(I)==false) {
          cout << "\n Mistake II in Prepare_communication_interior_stencils() "
               << endl;
        }
        if(I.Tiefe()!=level) {
          cout << "\n Mistake III in Prepare_communication_interior_stencils() "
               << endl;
        }
      }
      
      stencil_receive[level][d][counter_stencil] =
        Give_pointer_stencil(I);
      if(developer_version) {
        if(stencil_receive[level][d][counter_stencil]==NULL
           && Give_max_num_coarse_int_cells()>0) {
          //       && Give_cell_typ(I)==bo_cell) 
          cout << "\n Mistake IIV in Prepare_communication_interior_stencils() "
               << Give_cell_typ(I) 
               << " level: " << I.Tiefe()
               << " x: " << I.I_x().get()
               << " y: " << I.I_y().get()
               << " z: " << I.I_z().get()  
               << " my_rank: " << my_rank << endl;
          I.coordinate().Print();
          cout << endl;
        }
      }
      
      if(stencil_receive[level][d][counter_stencil]==NULL
         //      && Give_cell_typ(I)==ex_cell) {
         && ((I < my_index)==false)) {
        sp = newdouble(Give_max_num_coarse_int_cells());
        poi0 = hashtable0_point(I);
        poi0->var = sp;
        poi0->typ = bo_cell;
        stencil_receive[level][d][counter_stencil] = sp;
      }
      
      counter_stencil++;
    }
  }
}

void Grid::Update_global_stencil(int num_stencil, int level) {
  int num, iv;
  int d;
  MPI_Request req[12];
  MPI_Status status[12];
  int num_message;
  double*  send_info;
  double*  receive_info;
  int number_send;
  int number_receive;
  int number_start;


  //  int number_interior_stencils;
  double** stencil_sendd;
  double** stencil_received;
  
  int rank_source, rank_destination;

  // 0. Step: number start:
  number_start = cell_number_of_stencil(num_stencil);  

  if(I_am_active()) {
    // I. Part: interior part
    for(d=0;d<26;++d) { 
      // 1. part: basic infos about send and recieve
      if(d<6) {
        rank_source      = give_next_rank_source((dir_3D)d,level-1);
        rank_destination = give_next_rank_destination((dir_3D)d,level-1);  
      }
      else if(d<18) {
        rank_source      = give_next_rank_source((Edges_cell)(d-6),level-1);
        rank_destination = give_next_rank_destination((Edges_cell)(d-6),level-1);  
      }
      else {
        rank_source      = give_next_rank_source((dir_sons)(d-18),level-1);
        rank_destination = give_next_rank_destination((dir_sons)(d-18),level-1);
      }
      
      // 2. Step: calc number to be send and receive
      number_send = 
        number_send_stencil[level][d] * 64;
      number_receive = 
        number_receive_stencil[level][d] * 64;
      
      // 3. Step: build informations about stencil(double) to be send  
      send_info    = Give_buffer_send(d,number_send); 
      receive_info = Give_buffer_receive(d,number_receive); 
      
      stencil_sendd=stencil_send[level][d];
      for(num=0;num<number_send_stencil[level][d];++num) {
        for(iv=0;iv<64;++iv) {
          send_info[num*64+iv] = 
            stencil_sendd[num][number_start+iv];
        }
      }
      
      // 4. Step: send doubles
      num_message = 0;
      if(rank_source != -1 && number_receive>0) {
        MPI_Irecv(receive_info ,number_receive,
                  MPI_DOUBLE,rank_source,10,comm,
                  &req[num_message]);
        ++num_message;
      }      
      if(rank_destination != -1 && number_send>0) {
        MPI_Isend(send_info,number_send,
                  MPI_DOUBLE,rank_destination,10,comm,
                  &req[num_message]);
        ++num_message;
      }
      MPI_Waitall(num_message,req,status);
      
      // 5. part: set informations about stencil(double)
      stencil_received=stencil_receive[level][d];
      for(num=0;num<number_receive_stencil[level][d];++num) {
        
        //  cout << " num: " << num << endl;
        
        for(iv=0;iv<64;++iv) {
          stencil_received[num][number_start+iv] = 
            receive_info[num*64+iv];
        }
      }
    }

    // II. Part: boundary part
    for(d=0;d<26;++d) { 
      // 1. part: basic infos about send and recieve
      if(d<6) {
        rank_source      = give_next_rank_source((dir_3D)d,level-1);
        rank_destination = give_next_rank_destination((dir_3D)d,level-1);  
      }
      else if(d<18) {
        rank_source      = give_next_rank_source((Edges_cell)(d-6),level-1);
        rank_destination = give_next_rank_destination((Edges_cell)(d-6),level-1);  
      }
      else {
        rank_source      = give_next_rank_source((dir_sons)(d-18),level-1);
        rank_destination = give_next_rank_destination((dir_sons)(d-18),level-1);
      }
      
      // 2. Step: calc number to be send and receive
      number_send = 
        number_send_bo_stencil[level][d] * 64;
      number_receive = 
        number_receive_bo_stencil[level][d] * 64;
      
      // 3. Step: build informations about bo_stencil(double) to be send  
      send_info    = Give_buffer_send(d,number_send); 
      receive_info = Give_buffer_receive(d,number_receive); 
      
      stencil_sendd=bo_stencil_send[level][d];
      for(num=0;num<number_send_bo_stencil[level][d];++num) {
        for(iv=0;iv<64;++iv) {
          send_info[num*64+iv] = 
            stencil_sendd[num][number_start+iv];
        }
      }
      
      // 4. Step: send doubles
      num_message = 0;
      if(rank_source != -1 && number_receive>0) {
        MPI_Irecv(receive_info ,number_receive,
                  MPI_DOUBLE,rank_source,10,comm,
                  &req[num_message]);
        ++num_message;
      }      
      if(rank_destination != -1 && number_send>0) {
        MPI_Isend(send_info,number_send,
                  MPI_DOUBLE,rank_destination,10,comm,
                  &req[num_message]);
        ++num_message;
      }
      MPI_Waitall(num_message,req,status);
      
      // 5. part: set informations about bo_stencil(double)
      stencil_received=bo_stencil_receive[level][d];
      for(num=0;num<number_receive_bo_stencil[level][d];++num) {
        
        //  cout << " num: " << num << endl;
        
        for(iv=0;iv<64;++iv) {
          stencil_received[num][number_start+iv] = 
            receive_info[num*64+iv];
        }
      }
    }
  }
}

---