lbmpi.c

基于格子Boltzmann方法开源可视化软件的源代码 具有很高的实用价值。对学习LBM方法及其软件开发非常游泳

//##############################################################################
//
// Copyright (C), 2005, Michael Sukop and Danny Thorne
//
// lbmpi.c
//
//  - Note that this is adapted from a previous implementation in an
//    old mcmp code from last summer (2004).
//

void lbmpi_construct( lbmpi_ptr lbmpi, lattice_ptr lattice, int argc, char **argv)
{
  int i, j;
  int ierr;

  // Routine inititialization calls.
  MPI_Init( &argc, &argv);
  MPI_Comm_size( MPI_COMM_WORLD, &lbmpi->NumProcs);
  MPI_Comm_rank( MPI_COMM_WORLD, &lbmpi->ProcID);

  // Determine coordinates (myX,myY) of subdomain in
  // PX-by-PY array of subdomains.
  lbmpi->myX = lbmpi->ProcID%PX;
  lbmpi->myY = (int)floor((double)lbmpi->ProcID/(double)PX);

  // Determine ID of processes with adjacent subdomains.
  lbmpi->NorthID = PX*((lbmpi->myY+1)%PY) + (lbmpi->myX);
  lbmpi->SouthID = PX*((lbmpi->myY-1+PY)%PY) + (lbmpi->myX);
  lbmpi->EastID  = PX*((lbmpi->myY)) + (lbmpi->myX+1)%PX;
  lbmpi->WestID  = PX*((lbmpi->myY)) + (lbmpi->myX-1+PX)%PX;

  // Say "hi".
  printf(
   "Proc %d (%d,%d) of %d (%d-by-%d) says, \"Hi!\" "
   "// (N,S,E,W) = (%d,%d,%d,%d)\n", 
    lbmpi->ProcID,   lbmpi->myX, lbmpi->myY, 
    lbmpi->NumProcs, PX,  PY,
    lbmpi->NorthID, lbmpi->SouthID, lbmpi->EastID, lbmpi->WestID );

  // Allocate space for datatypes.
  lbmpi_allocate_datatypes( lbmpi, lattice);

  // Create data type structure MPI_South2North.
  //
  //  - This will be used in the communication step after
  //    streaming to send/recv elements from north subdomains
  //    to south subdomains.
  //

  printf("\n%d: Making MPI_South2North... \n", lbmpi_get_ProcID(lbmpi));
  MPI_Address( 
    get_ftemp_ptr( lattice,
                   /*subs*/0, 
                   /*j*/get_sj(lattice), 
                   /*i*/get_si(lattice), 
                   /*a*/0),
    (MPI_Aint *)lbmpi_get_Index0_ptr(lbmpi));

  // Acquire memory addresses of all the elements.  With these
  // we will compute the indices needed by MPI_Type_struct 
  // below.
  for(i=get_si(lattice); i<=get_ei(lattice); i++)
  {
    MPI_Address( get_ftemp_ptr(lattice,0,get_sj(lattice),i,2), 
                 &( lbmpi->AddrsNS[ 3*(i-0)+0]) );
    MPI_Address( get_ftemp_ptr(lattice,0,get_sj(lattice),i,5),
                 &( lbmpi->AddrsNS[ 3*(i-0)+1]) );
    MPI_Address( get_ftemp_ptr(lattice,0,get_sj(lattice),i,6),
                 &( lbmpi->AddrsNS[ 3*(i-0)+2]) );
  
  } /* for(i=1; i<=get_LX(lattice); i++) */
  
  // Stuff needed by MPI_Type_struct.
  for(i=1; i<=get_LX(lattice); i++)
  {
    // All the block lengths are one (1).
    lbmpi->BlockLengthsNS[ 3*(i-1)+0] = 1;
    lbmpi->BlockLengthsNS[ 3*(i-1)+1] = 1;
    lbmpi->BlockLengthsNS[ 3*(i-1)+2] = 1;

    // Compute offsets from the first element.
    lbmpi->IndicesNS[ 3*(i-1)+0] = lbmpi->AddrsNS[ 3*(i-1)+0]-lbmpi->AddrsNS[0];
    lbmpi->IndicesNS[ 3*(i-1)+1] = lbmpi->AddrsNS[ 3*(i-1)+1]-lbmpi->AddrsNS[0];
    lbmpi->IndicesNS[ 3*(i-1)+2] = lbmpi->AddrsNS[ 3*(i-1)+2]-lbmpi->AddrsNS[0];

    // All the types are doubles.
    lbmpi->TypesNS[ 3*(i-1)+0] = MPI_DOUBLE;
    lbmpi->TypesNS[ 3*(i-1)+1] = MPI_DOUBLE;
    lbmpi->TypesNS[ 3*(i-1)+2] = MPI_DOUBLE;

  } /* for(i=1; i<=get_LX(lattice); i++) */


  ierr = MPI_Type_struct(
           /* int count */                3*get_LX(lattice),
           /* int blocklens[] */          lbmpi->BlockLengthsNS,
           /* MPI_Aint indices[] */       lbmpi->IndicesNS,
           /* MPI_Datatype old_types[] */ lbmpi->TypesNS,
           /* MPI_Datatype *newtype */    &lbmpi->MPI_South2North    );

  ierr = MPI_Type_commit( 
           /* MPI_Datatype *datatype */   &lbmpi->MPI_South2North    );

  printf("\n%d: Done making MPI_South2North\n", lbmpi_get_ProcID(lbmpi));

#if 1
  // Output the indices for inspection...
  sprintf( lbmpi->iobuf, " ");
  for(i=1; i<=get_LX(lattice); i++)
  {
    sprintf( lbmpi->iobuf, "%s%d %d %d ",
      lbmpi->iobuf,
      lbmpi->IndicesNS[ 3*(i-1)+0],
      lbmpi->IndicesNS[ 3*(i-1)+1],
      lbmpi->IndicesNS[ 3*(i-1)+2] );
  }
  printf("\n%d: MPI_South2North { %s }\n", 
    lbmpi_get_ProcID(lbmpi), 
    lbmpi->iobuf);
#endif

  // Create data type structure MPI_North2South.
  //
  //  - This will be used in the communication step after
  //    streaming to send/recv elements from south subdomains
  //    to north subdomains.
  //

  printf("\n%d: Making MPI_North2South... \n", lbmpi_get_ProcID(lbmpi));
  MPI_Address( 
    get_ftemp_ptr( lattice,
                   /*subs*/0, 
                   /*j*/get_ej(lattice), 
                   /*i*/get_si(lattice), 
                   /*a*/4),
    (MPI_Aint *)lbmpi_get_Index0_ptr(lbmpi));

  // Acquire memory addresses of all the elements.  With these
  // we will compute the indices needed by MPI_Type_struct 
  // below.
  for(i=get_si(lattice); i<=get_ei(lattice); i++)
  {
    MPI_Address( 
      get_ftemp_ptr(lattice,0,get_ej(lattice),i,4),
      &( lbmpi->AddrsNS[ 3*(i-0)+0]));
    MPI_Address( 
      get_ftemp_ptr(lattice,0,get_ej(lattice),i,7),
      &( lbmpi->AddrsNS[ 3*(i-0)+1]));
    MPI_Address( 
      get_ftemp_ptr(lattice,0,get_ej(lattice),i,8),
      &( lbmpi->AddrsNS[ 3*(i-0)+2]));
  
  } /* for(i=1; i<=get_LX(lattice); i++) */

  // Stuff needed by MPI_Type_struct.
  for(i=1; i<=get_LX(lattice); i++)
  {
    // All the block lengths are one (1).
    lbmpi->BlockLengthsNS[ 3*(i-1)+0] = 1;
    lbmpi->BlockLengthsNS[ 3*(i-1)+1] = 1;
    lbmpi->BlockLengthsNS[ 3*(i-1)+2] = 1;

    // Compute offsets from the first element.
    lbmpi->IndicesNS[ 3*(i-1)+0] = lbmpi->AddrsNS[ 3*(i-1)+0]-lbmpi->AddrsNS[0];
    lbmpi->IndicesNS[ 3*(i-1)+1] = lbmpi->AddrsNS[ 3*(i-1)+1]-lbmpi->AddrsNS[0];
    lbmpi->IndicesNS[ 3*(i-1)+2] = lbmpi->AddrsNS[ 3*(i-1)+2]-lbmpi->AddrsNS[0];

    // All the types are doubles.
    lbmpi->TypesNS[ 3*(i-1)+0] = MPI_DOUBLE;
    lbmpi->TypesNS[ 3*(i-1)+1] = MPI_DOUBLE;
    lbmpi->TypesNS[ 3*(i-1)+2] = MPI_DOUBLE;

  } /* for(i=1; i<=get_LX(lattice); i++) */

  ierr = MPI_Type_struct(
           /* int count */                3*get_LX(lattice),
           /* int blocklens[] */          lbmpi->BlockLengthsNS,
           /* MPI_Aint indices[] */       lbmpi->IndicesNS,
           /* MPI_Datatype old_types[] */ lbmpi->TypesNS,
           /* MPI_Datatype *newtype */    &lbmpi->MPI_North2South    );

  ierr = MPI_Type_commit( 
           /* MPI_Datatype *datatype */   &lbmpi->MPI_North2South    );

  printf("\n%d: Done making MPI_North2South\n", lbmpi->ProcID);

#if 1
  // Output the indices for inspection...
  sprintf( lbmpi->iobuf, " ");
  for(i=1; i<=get_LX(lattice); i++)
  {
    sprintf( lbmpi->iobuf, "%s%d %d %d ",
      lbmpi->iobuf,
      lbmpi->IndicesNS[ 3*(i-1)+0],
      lbmpi->IndicesNS[ 3*(i-1)+1],
      lbmpi->IndicesNS[ 3*(i-1)+2] );
  }
  printf("\n%d: MPI_North2South { %s }\n", 
    lbmpi_get_ProcID(lbmpi),
    lbmpi->iobuf);
#endif

  // Create data type structure MPI_East2West.
  //
  //  - This will be used in the communication step after
  //    streaming to send/recv elements from east subdomains
  //    to west subdomains.
  //

  printf("\n%d: Making MPI_East2West... \n", lbmpi_get_ProcID(lbmpi));
  //MPI_Address( &( ftemp[1][1][get_LX(lattice)][3]), &Index0);
  MPI_Address( 
    get_ftemp_ptr( lattice,
                   /*subs*/0, 
                   /*j*/get_sj(lattice), 
                   /*i*/get_ei(lattice), 
                   /*a*/3),
    (MPI_Aint *)lbmpi_get_Index0_ptr(lbmpi));

  // Acquire memory addresses of all the elements.  With these
  // we will compute the indices needed by MPI_Type_struct 
  // below.
  for(j=get_sj(lattice); j<=get_ej(lattice); j++)
  {
    MPI_Address(
      get_ftemp_ptr(lattice,0,j,get_ei(lattice),3),
      &( lbmpi->AddrsEW[ 3*(j-0)+0]));
    MPI_Address( 
      get_ftemp_ptr(lattice,0,j,get_ei(lattice),6),
      &( lbmpi->AddrsEW[ 3*(j-0)+1]));
    MPI_Address( 
      get_ftemp_ptr(lattice,0,j,get_ei(lattice),7),
      &( lbmpi->AddrsEW[ 3*(j-0)+2]));

  } /* for(j=1; j<=get_LY(lattice); j++) */

  // Stuff needed by MPI_Type_struct.
  for(j=1; j<=get_LY(lattice); j++)
  {
    // All the block lengths are one (1).
    lbmpi->BlockLengthsEW[ 3*(j-1)+0] = 1;
    lbmpi->BlockLengthsEW[ 3*(j-1)+1] = 1;
    lbmpi->BlockLengthsEW[ 3*(j-1)+2] = 1;

    // Compute offsets from the first element.
    lbmpi->IndicesEW[ 3*(j-1)+0] = lbmpi->AddrsEW[ 3*(j-1)+0]-lbmpi->AddrsEW[0];
    lbmpi->IndicesEW[ 3*(j-1)+1] = lbmpi->AddrsEW[ 3*(j-1)+1]-lbmpi->AddrsEW[0];
    lbmpi->IndicesEW[ 3*(j-1)+2] = lbmpi->AddrsEW[ 3*(j-1)+2]-lbmpi->AddrsEW[0];

    // All the types are doubles.
    lbmpi->TypesEW[ 3*(j-1)+0] = MPI_DOUBLE;
    lbmpi->TypesEW[ 3*(j-1)+1] = MPI_DOUBLE;
    lbmpi->TypesEW[ 3*(j-1)+2] = MPI_DOUBLE;

  } /* for(j=1; j<=get_LY(lattice); j++) */

  ierr = MPI_Type_struct(
           /* int count */                3*get_LY(lattice),
           /* int blocklens[] */          lbmpi->BlockLengthsEW,
           /* MPI_Aint indices[] */       lbmpi->IndicesEW,
           /* MPI_Datatype old_types[] */ lbmpi->TypesEW,
           /* MPI_Datatype *newtype */    &lbmpi->MPI_East2West    );

  ierr = MPI_Type_commit( 
           /* MPI_Datatype *datatype */   &lbmpi->MPI_East2West    );

  printf("\n%d: Done making MPI_East2West\n", lbmpi_get_ProcID(lbmpi));

#if 1
  // Output the indices for inspection...
  sprintf( lbmpi->iobuf, " ");
  for(j=1; j<=get_LY(lattice); j++)
  {
    sprintf( lbmpi->iobuf, "%s%d %d %d ",
      lbmpi->iobuf,
      lbmpi->IndicesEW[ 3*(j-1)+0],
      lbmpi->IndicesEW[ 3*(j-1)+1],
      lbmpi->IndicesEW[ 3*(j-1)+2] );
  }
  printf("\n%d: MPI_East2West { %s }\n", lbmpi_get_ProcID(lbmpi), lbmpi->iobuf);
#endif

#if 0

// Create data type structure MPI_West2East.
//
//  - This will be used in the communication step after
//    streaming to send/recv elements from west subdomains
//    to east subdomains.
//

printf("\n%d: Making MPI_West2East... \n", ProcID);
MPI_Address( &( ftemp[1][1][1][1]), &Index0);

// Acquire memory addresses of all the elements.  With these
// we will compute the indices needed by MPI_Type_struct 
// below.
for(j=1; j<=get_LY(lattice); j++)
{
  MPI_Address( &( ftemp[1][j][1][1]), &( AddrsEW[ 3*(j-1)+0]));
  MPI_Address( &( ftemp[1][j][1][5]), &( AddrsEW[ 3*(j-1)+1]));
  MPI_Address( &( ftemp[1][j][1][8]), &( AddrsEW[ 3*(j-1)+2]));

} /* for(j=1; j<=get_LY(lattice); j++) */

// Stuff needed by MPI_Type_struct.
for(j=1; j<=get_LY(lattice); j++)
{
  // All the block lengths are one (1).
  BlockLengthsEW[ 3*(j-1)+0] = 1;
  BlockLengthsEW[ 3*(j-1)+1] = 1;
  BlockLengthsEW[ 3*(j-1)+2] = 1;

  // Compute offsets from the first element.
  IndicesEW[ 3*(j-1)+0] = AddrsEW[ 3*(j-1)+0]-AddrsEW[0];
  IndicesEW[ 3*(j-1)+1] = AddrsEW[ 3*(j-1)+1]-AddrsEW[0];
  IndicesEW[ 3*(j-1)+2] = AddrsEW[ 3*(j-1)+2]-AddrsEW[0];

  // All the types are doubles.
  TypesEW[ 3*(j-1)+0] = MPI_DOUBLE;
  TypesEW[ 3*(j-1)+1] = MPI_DOUBLE;
  TypesEW[ 3*(j-1)+2] = MPI_DOUBLE;

} /* for(j=1; j<=get_LY(lattice); j++) */

ierr = MPI_Type_struct(
         /* int count */                3*get_LY(lattice),
         /* int blocklens[] */          BlockLengthsEW,
         /* MPI_Aint indices[] */       IndicesEW,
         /* MPI_Datatype old_types[] */ TypesEW,
         /* MPI_Datatype *newtype */    &MPI_West2East    );

ierr = MPI_Type_commit( 
         /* MPI_Datatype *datatype */   &MPI_West2East    );

printf("\n%d: Done making MPI_West2East\n", ProcID);

#if 0
// Output the indices for inspection...
sprintf( iobuf, " ");
for(j=1; j<=get_LY(lattice); j++)
{
  sprintf( iobuf, "%s%d %d %d ",
    iobuf,
    IndicesEW[ 3*(j-1)+0],
    IndicesEW[ 3*(j-1)+1],
    IndicesEW[ 3*(j-1)+2] );
}
printf("\n%d: MPI_West2East { %s }\n", ProcID, iobuf);
#endif

#endif
} /* void lbmpi_construct( lbmpi_ptr lbmpi, lattice_ptr lattice, argc, argv) */

// void lbmpi_allocate_data_structures( lbmpi, lattice)
//
// Allocate datatypes for 
void lbmpi_allocate_datatypes( lbmpi_ptr lbmpi, lattice_ptr lattice)
{
  //int BlockLengthsEW[3*LY];
  lbmpi->BlockLengthsEW = (int*)malloc( 3*get_LY(lattice)*sizeof(int));

  //MPI_Aint AddrsEW[3*LY];
  lbmpi->AddrsEW = (MPI_Aint*)malloc( 3*get_LY(lattice)*sizeof(MPI_Aint));

  //MPI_Aint IndicesEW[3*LY];
  lbmpi->IndicesEW = (MPI_Aint*)malloc( 3*get_LY(lattice)*sizeof(MPI_Aint));

  //MPI_Datatype TypesEW[3*LY];
  lbmpi->TypesEW = (MPI_Datatype*)malloc( 3*get_LY(lattice)*sizeof(MPI_Datatype));



  //int BlockLengthsNS[3*LX];
  lbmpi->BlockLengthsNS = (int*)malloc( 3*get_LX(lattice)*sizeof(int));

  //MPI_Aint AddrsNS[3*LX];
  lbmpi->AddrsNS = (MPI_Aint*)malloc( 3*get_LX(lattice)*sizeof(MPI_Aint));

  //MPI_Aint IndicesNS[3*LX];
  lbmpi->IndicesNS = (MPI_Aint*)malloc( 3*get_LX(lattice)*sizeof(MPI_Aint));

  //MPI_Datatype TypesNS[3*LX];
  lbmpi->TypesNS = (MPI_Datatype*)malloc( 3*get_LX(lattice)*sizeof(MPI_Datatype));


#if DUMP_AFTER_COMM || DUMP_BEFORE_COMM
  //double fmat[3*LY][3*LX];
  fmat = (double**)malloc(3*get_LY(lattice)*sizeof(double*))
  for(j=1;j<=get_LY(lattice);j++)
  {
    lbmpi->fmat[j] = (double*)malloc(3*get_LX(lattice)*sizeof(double));
  }

#endif /* DUMP_AFTER_COMM || DUMP_BEFORE_COMM */

} /* void lbmpi_allocate_datatypes( lbmpi, lattice) */

int lbmpi_get_ProcID( lbmpi_ptr lbmpi)
{
  return lbmpi->ProcID;
}

MPI_Aint *lbmpi_get_Index0_ptr( lbmpi_ptr lbmpi)
{
  return (MPI_Aint *)(&(lbmpi->Index0));
}