acoustic_body_coefs()

int acoustic_body_coefs

(

)

Evaluates B and C coefficients (please refer to the user manual for details on the theory) The actual evaluation is done in parallel on all the nodes of the domain (SPMD). The matrices are partitioned according to scalapack 2D block cyclic distribution (locRows x locCols). In each loop over the nodes, krow lines of the coefficient matrices can be loaded into memory to minimize memory usage and then written on two different files:

bbcoefs.bin.<nodeid>  B coefficients
bccoefs.bin.<nodeid>  C coefficients

where <nodeid> is the rank of the node (to avoid name-clashing on multiple process running on the same node). The same files are then read in the Matrices assembly procedure (matrices.c), basically we dictate here the overall partitioning method for the run.

See also

ac_coef_mics.c

matrices.c

integrals.c

                         {

   int npsib,nphib,nelmu;
   int icntr,ielmb;
   int ielmu,is;
   int mpi_nelmb;
   struct panel4* mpi_elements;
   int mpi_ncntr;
   struct vector* mpi_cntr;

   int kode;
   int bstride;
   int cstride;

   DOUBLE theta;
   DOUBLE source,doublet,ratelet;

   int gRows,gCols;     // global dim of coefficient matrix
   int locRows,locCols; // local dim of coefficient matrix

   int izero=0;

//#define _ACOUSTO_DEBUG
#ifdef _ACOUSTO_DEBUG
   FILE* f;
   char fname[15];
#endif

   FILE* fbcoef = NULL;
   FILE* fccoef = NULL;
   int nrow,nrowi,iloop,nloop,irow,nio,nrlast;
   char fbname[MAX_PATH];
   char fcname[MAX_PATH];

   DOUBLE* ctmp;
   int i,j;
   int ic;
   int bsymmoffset,csymmoffset;

   long long int perc,ioldcurrent,icurrent,ntotal;


   logger(LOG_INFO,MSG_BODYCOEF_START);

   // assigning number of psi and phi points (coincident with number of elements)
   npsib = modgeominfo->nelmb;
   nphib = modgeominfo->nelmb;
   nelmu = modgeominfo->ncntr;


   /* copying variables into local storage (just for ease of use) */
   mpi_nelmb    = modgeominfo->nelmb;
   mpi_elements = geometry->elements;
   mpi_ncntr    = modgeominfo->ncntr + modgeominfo->nchief;
   mpi_cntr     = geometry->cntrc;


#ifdef _ACOUSTO_DEBUG
  if(rank == 0){
    sprintf(fname,"elements.%d",rank);
    logger(LOG_DEBUG,"filename [%s]\n",fname);
    f = fopen(fname,"w");
    for(i=0;i<modgeominfo->nelmb;i++){
      fprintf(f,"# ielmb=%d\n",i);
      fprintf(f,VECTOR_PRINTF,mpi_elements[i].x[0].x, mpi_elements[i].x[0].y, mpi_elements[i].x[0].z);
      fprintf(f,VECTOR_PRINTF,mpi_elements[i].x[1].x, mpi_elements[i].x[1].y, mpi_elements[i].x[1].z);
      fprintf(f,VECTOR_PRINTF,mpi_elements[i].x[2].x, mpi_elements[i].x[2].y, mpi_elements[i].x[2].z);
      fprintf(f,VECTOR_PRINTF,mpi_elements[i].x[3].x, mpi_elements[i].x[3].y, mpi_elements[i].x[3].z);
      fprintf(f,"\n\n");
    }
    fclose(f);
    sprintf(fname,"cntr.%d",rank);
    logger(LOG_DEBUG,"filename [%s]\n",fname);
    f = fopen(fname,"w");
    for(i=0;i<mpi_ncntr;i++){
      fprintf(f,VECTOR_PRINTF,mpi_cntr[i].x, mpi_cntr[i].y, mpi_cntr[i].z);
    }
    fclose(f);
    sprintf(fname,"normals.%d",rank);
    logger(LOG_DEBUG,"filename [%s]\n",fname);
    f = fopen(fname,"w");
    for(i=0;i<mpi_ncntr;i++){
      fprintf(f,VECTOR_PRINTF,mpi_cntr[i].x, mpi_cntr[i].y, mpi_cntr[i].z);
      fprintf(f,VECTOR_PRINTF,mpi_cntr[i].x+mpi_elements[i].n.x, 
                              mpi_cntr[i].y+mpi_elements[i].n.y, 
                              mpi_cntr[i].z+mpi_elements[i].n.z);
      fprintf(f,"\n\n");
    }
    fclose(f);
  }
#endif
   gRows = mpi_ncntr;
   //gCols = mpi_nelmb;
   gCols = modgeominfo->ncntr;

   bodycoefs = calloc(1,sizeof(struct bodycoefs_struct)); // allocating main body coefs struct

   // Building array descriptor
   locRows = numroc_(&gRows,&runinfo->row_block_size,&runinfo->myrow,&izero,&runinfo->nprows);
   locCols = numroc_(&gCols,&runinfo->col_block_size,&runinfo->mycol,&izero,&runinfo->npcols);

   // nrow lines each loop
   nrow  = runinfo->krow;
   if(runinfo->krow <= 0 || runinfo->krow >= locRows) nrow = locRows;

   // ALLOCATING COEFS
   bodycoefs->B =calloc(2*nrow*locCols*runinfo->nsymm,sizeof(DOUBLE));  // 2 B coefficients, nrow*partition size*symmetry factor
   bodycoefs->C =calloc(4*nrow*locCols*runinfo->nsymm,sizeof(DOUBLE));  // 4 C coefficients, nrow*partition size*symmetry factor
   logger(LOG_DEBUG,"Total allocated memory for body coefficients=%d bytes\n", 
                     (2*nrow*locCols*runinfo->nsymm*sizeof(DOUBLE) + 4*nrow*locCols*runinfo->nsymm*sizeof(DOUBLE)));

   ctmp = calloc(locRows,sizeof(DOUBLE));


  if(runinfo->krow > 0){ // only if krow is used
   // setting names of files (different for each node in case of multi-process run)
   get_filename(fbname,"bbcoefs.bin.%d",rank); 
   get_filename(fcname,"bccoefs.bin.%d",rank);
   logger(LOG_DEBUG,"body coefficients will be in files \"%s\" and \"%s\"\n", fbname,fcname);
   // opening files in binary-write mode
   fbcoef = fopen(fbname,"wb"); 
   fccoef = fopen(fcname,"wb");
  }

  // evaluating local block of coefficient
  j12(locRows,nrow,&nrlast,&nloop); // calculating start row and number of loops

  ntotal = (nloop-1) * nrow * locCols * runinfo->nsymm;
  ioldcurrent = 0;
  icurrent = 0;
  ioldcurrent = -1;

  for(iloop=0;iloop<nloop;iloop++){  // LOOP ON NLOOP
   ic = 0;
   nrowi = nrow;
   if(iloop == nloop-1) nrowi = nrlast; // setting number of rows=number of last rows on final loop
   //calculating stride for arrays
   bstride = nrowi*locCols;
   cstride = nrowi*locCols;

   for(irow=0;irow<nrowi;irow++){    // LOOP ON NROW
    i= j21(irow,iloop,nrow); // calculating control point index
    sc_l2g(i,runinfo->myrow,gRows,runinfo->nprows,runinfo->row_block_size,&icntr);
    for(j=0;j<locCols;j++){
      // getting local indices

      sc_l2g(j,runinfo->mycol,gCols,runinfo->npcols,runinfo->col_block_size,&ielmu);
      for(is=0;is<runinfo->nsymm;is++){
       
       bsymmoffset = is * 2*nrowi*locCols;
       csymmoffset = is * 4*nrowi*locCols;
       ielmb = ielmu + is*nelmu; 

       kode = (ielmb == icntr) ? 1 : 0; // singular panels flag
       
       delaya(&theta,ielmb,&mpi_elements[ielmb],mpi_cntr[icntr]);                               // evaluating acoustic delay
       intgba(1, 0, kode, &mpi_cntr[icntr], &mpi_elements[ielmb], &source, &doublet, &ratelet); // integrating

       ic = irow*locCols +j;
        
       bodycoefs->B[ic             + bsymmoffset]  = source;
       bodycoefs->B[ic + bstride   + bsymmoffset]  = theta/runinfo->vsound;
       
       bodycoefs->C[ic             + csymmoffset]  = doublet;
       bodycoefs->C[ic + cstride   + csymmoffset]  = ((1 == kode) ? (DOUBLE)0.0 : (DOUBLE)ratelet/runinfo->vsound );
       bodycoefs->C[ic + cstride*2 + csymmoffset]  = theta/runinfo->vsound;

       
       bodycoefs->C[ic + cstride*3 + csymmoffset] = (icntr == ielmb) ? 0.5 : 0.0;
//       ic++;

       icurrent++;
       if( (icurrent*10/ntotal) != (ioldcurrent*10/ntotal)){
         perc =icurrent*100/ntotal;
         logger(LOG_INFO,"Completed %lld \%% \n",perc);
       }
       ioldcurrent = icurrent;
      }
    
    
    } // LOOP ON locCOLS END
   } // LOOP ON NROW END
   if(runinfo->krow > 0){
     // append on files
     nio = fwrite(bodycoefs->B,sizeof(DOUBLE),2*nrowi*locCols*runinfo->nsymm,fbcoef); 
     nio = fwrite(bodycoefs->C,sizeof(DOUBLE),4*nrowi*locCols*runinfo->nsymm,fccoef); 
   }
  } // LOOP ON NLOOP END

   if(runinfo->krow > 0){
    fclose(fbcoef);
    fclose(fccoef);
   }
  //zgsum2d_(&runinfo->ctxt,"R"," ",&locRows,&runinfo->npcols,ctmp,&locRows,&ione,&ione,&imone,&imone,&imone);
  /*
  Cdgsum2d(runinfo->ctxt,"R"," ",locRows,1,ctmp,locRows,-1,-1);


  for(icntr=0;icntr<gRows;icntr++){
   sc_g2l(icntr,gRows,runinfo->nprows,runinfo->row_block_size,&myi,&i);
   sc_g2l(icntr,gCols,runinfo->npcols,runinfo->col_block_size,&myj,&j);
   if(myi == runinfo->myrow && 
      myj == runinfo->mycol){
     bodycoefs->C[(i*locCols+j) + cstride*3]  = 1.0 + ctmp[i]; // setting C[3] value
   } 
  }
  */

  // out of main loop for diagonal element (must optimize this)
  /*
  for(icntr=0;icntr<gRows;icntr++){
   dsum = 0.0;
   for(ielmb=0;ielmb<gCols;ielmb++){ // summing on col
     kode = (ielmb == icntr) ? 1 : 0; // singular panels flag
     intgba(1, 0, kode, &mpi_cntr[icntr], &mpi_elements[ielmb], &source, &doublet, &ratelet); // integrating
     dsum += doublet;
   }
   sc_g2l(icntr,gRows,runinfo->nprows,runinfo->row_block_size,&myi,&i);
   sc_g2l(icntr,gCols,runinfo->npcols,runinfo->col_block_size,&myj,&j);
   if(myi == runinfo->myrow && 
      myj == runinfo->mycol){
     bodycoefs->C[(i*locCols+j) + cstride*3]  = 1.0 + dsum; // setting C[3] value
   } 
  }
  */

#ifdef _ACOUSTO_DEBUG
 if(runinfo->krow > 0){ // only if krow is used
   get_filename(fbname,"bbcoefs.bin.%d",rank); get_filename(fcname,"bccoefs.bin.%d",rank);
   fbcoef = fopen(fbname,"rb"); fccoef = fopen(fcname,"rb");
 }

 sprintf(fname,"bcoef.%d%d",runinfo->myrow,runinfo->mycol);
 f = fopen(fname,"w");
  j12(locRows,nrow,&nrlast,&nloop);
  for(iloop=0;iloop<nloop;iloop++){  // LOOP ON NLOOP
   nrowi = nrow;
   if(iloop == nloop-1) nrowi = nrlast;
   bstride = nrowi*locCols;
   cstride = nrowi*locCols;
   if(runinfo->krow > 0){ // only if krow is used
     nio = fread(bodycoefs->B,sizeof(DOUBLE),2*nrowi*locCols*runinfo->nsymm,fbcoef); 
     nio = fread(bodycoefs->C,sizeof(DOUBLE),4*nrowi*locCols*runinfo->nsymm,fccoef); 
   }
   for(irow=0;irow<nrowi;irow++){    // LOOP ON NROW
   i= j21(irow,iloop,nrow);
   sc_l2g(i,runinfo->myrow,gRows,runinfo->nprows,runinfo->row_block_size,&icntr);
   for(j=0;j<locCols;j++){
      sc_l2g(j,runinfo->mycol,gCols,runinfo->npcols,runinfo->col_block_size,&ielmu);
      for(is=0;is<runinfo->nsymm;is++){
       
       ielmb = ielmu + is*nelmu;
       bsymmoffset = is * 2*nrowi*locCols;
       csymmoffset = is * 4*nrow*locCols;
       fprintf(f,"%d %d %f %f %f %f %f %f\n",icntr,ielmb,
                bodycoefs->B[(irow*locCols+j)             + bsymmoffset],
                bodycoefs->B[(irow*locCols+j) + bstride   + bsymmoffset],
                bodycoefs->C[(irow*locCols+j)             + csymmoffset],
                bodycoefs->C[(irow*locCols+j) + cstride   + csymmoffset],
                bodycoefs->C[(irow*locCols+j) + cstride*2 + csymmoffset],
                bodycoefs->C[(irow*locCols+j) + cstride*3 + csymmoffset]);
      }

   }
  }
  }
  if(runinfo->krow > 0){ // only if krow is used
    fclose(fbcoef);
    fclose(fccoef);
  }


 fclose(f);
#endif
#undef _ACOUSTO_DEBUG
 if(runinfo->krow > 0){
  free(bodycoefs->B);
  free(bodycoefs->C);
 }

 free(ctmp);
 logger(LOG_INFO,MSG_BODYCOEF_END);
 return 1;

}