pre_mic()

void pre_mic	(	COMPLEX	cfreq,
		int	icounter
	)

Evaluates solution at microphones. Coefficients are read on disk. This function is called when pre_calculate_coefs == 1.

Parameters

[in]	cfreq	COMPLEX frequency
[in]	icounter	frequency index

                                        {
  int i;
  int iphiu;
  int nphiu;
  int nphib;

  int mpi_ncntr,mpi_nelmu;
  int bstride,cstride;

  int locRows,locCols,locCRows;
  int gRows,gCols,gCRows;
  int ielmb,ielmu;
  int j;
  int descpsiu[9];
  int descgpsiu[9];
  int mm,nn;
  int ione  = 1;
  int izero = 0;
  int info;
  struct vector g;
  struct vector un;
  struct panel4* mpi_elements;

  DOUBLE gdotn;


  COMPLEX cexpo;
  DOUBLE bb,cc,dd,theta;
  COMPLEX rho_s;

  Vec psincvec;
  Vec phincvec;

  COMPLEX phiscval  = COMPLEX_ZERO; // Cauchy Data
  COMPLEX phincval  = COMPLEX_ZERO;
  COMPLEX phitcval  = COMPLEX_ZERO;

  COMPLEX psincval  = COMPLEX_ZERO;
  COMPLEX psiscval  = COMPLEX_ZERO;
  COMPLEX psitcval  = COMPLEX_ZERO;

  COMPLEX phincmval = COMPLEX_ZERO; // Potential at microphones
  COMPLEX phiscmval = COMPLEX_ZERO;
  COMPLEX phitmval  = COMPLEX_ZERO;

  COMPLEX prescmval = COMPLEX_ZERO;
  COMPLEX pretmval  = COMPLEX_ZERO;

  DOUBLE* ctmp;
  int ic,is;
  int bsymmoffset,csymmoffset;


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

#ifdef _ACOUSTO_DEBUG
  int imics;
  char fname[MAX_PATH];
  FILE* f;
#endif

  COMPLEX lambda;
  COMPLEX gamma;
  COMPLEX func;
  COMPLEX gamma_over_lambda;
  COMPLEX lambda_over_gamma;

  mpi_elements = geometry->elements;
  mpi_ncntr  = modgeominfo->ncntr + modgeominfo->nchief;
  mpi_nelmu  = modgeominfo->ncntr;
  nphiu      = modgeominfo->ncntr;

  gRows      = modgeominfo->nmics;
  gCols      = modgeominfo->ncntr;
  gCRows     = nphiu;

  locRows  = numroc_(&gRows  ,&runinfo->mics_block_size,&runinfo->myrow,&izero,&runinfo->nprows);
  locCols  = numroc_(&gCols  ,&runinfo->col_block_size,&runinfo->mycol,&izero,&runinfo->npcols);
  locCRows = numroc_(&gCRows ,&runinfo->row_block_size,&runinfo->myrow,&izero,&runinfo->nprows);

  nrow  = runinfo->krow;
  if(runinfo->krow <=0 || runinfo->krow >= locRows) nrow = locRows;

  if(runinfo->krow > 0){
   micscoefs->B = calloc(2*nrow*locCols*runinfo->nsymm,sizeof(DOUBLE));
   micscoefs->C = calloc(3*nrow*locCols*runinfo->nsymm,sizeof(DOUBLE));
  }
  ctmp = calloc(locRows,sizeof(DOUBLE));

  nphib = modgeominfo->nelmb;

  if(0==icounter){
    CREATE_VEC(solution->vecPhium ,locRows);
    CREATE_VEC(solution->vecPrescm,locRows);
    CREATE_VEC(solution->vecPretm ,locRows);
  }
  VecZeroEntries(solution->vecPhium  ,locRows);
  VecZeroEntries(solution->vecPrescm ,locRows);
  VecZeroEntries(solution->vecPretm  ,locRows);

  psincvec = calloc(nphiu,sizeof(COMPLEX)); // allocating vector
  phincvec = calloc(nphiu,sizeof(COMPLEX)); // allocating vector

  // Gathering psiu on node 0 and broadcasting
  // Collecting solution on node zero and broadcasting to everyone
  mm = gCRows;// * runinfo->nprows;
  nn = runinfo->npcols;
  descinit_(descpsiu   , &gCRows  ,&ione, &runinfo->row_block_size  ,&ione,&izero,&izero,&runinfo->ctxt,&gCRows  ,&info);
  descinit_(descgpsiu  , &mm      ,&ione, &gCRows                   ,&ione,&izero,&izero,&runinfo->ctxt,&gCRows  ,&info);
  pzgemr2d_(&gCRows,&ione,solution->vecPsinc,&ione, &ione,descpsiu, psincvec,&ione,&ione,descgpsiu,&runinfo->ctxt);
  pzgemr2d_(&gCRows,&ione,solution->vecPhinc,&ione, &ione,descpsiu, phincvec,&ione,&ione,descgpsiu,&runinfo->ctxt);

  MPI_Bcast(psincvec,gCRows,MPI_ACOUSTO_COMPLEX,0,MPI_COMM_WORLD);
  MPI_Bcast(phincvec,gCRows,MPI_ACOUSTO_COMPLEX,0,MPI_COMM_WORLD);

  rho_s = cfreq*modsolinfo->rho;


  /* Since we are splitting also on columns Pressure on microphones  must be summed on the columns
   * of the grid.
   * e.g.:
   *  for pressure on microphones in row 0 of the process grid (PRE00, PRE01, PRE02) 
   *  the actual pressure on the microphones is the sum on the colums (PRE00+PRE01+PRE02)
   *
   * This is due to the choice of splitting the body elements along the column of the process
   * grid and to the fact that the evaluation of pressure is simply additive.
   *
   * This also applies to potential at microphones
   *
   *         0          1         2  
   *   --------------------------------
   *   |           |         |        |
   * 0 |   PRE00   |  PRE01  |  PRE02 |
   *   |           |         |        |
   *   --------------------------------
   *   |           |         |        |
   * 1 |   PRE11   |  PRE11  |  PRE12 |
   *   |           |         |        |
   *   --------------------------------
   * 
   **/

  if(runinfo->krow > 0){
    get_filename(fbname,"mbcoefs.bin.%d",rank); 
    get_filename(fcname,"mccoefs.bin.%d",rank);
    fbcoef = fopen(fbname,"rb"); fccoef = fopen(fcname,"rb");
  }

  j12(locRows,nrow,&nrlast,&nloop);
  for(iloop=0;iloop<nloop;iloop++){  // LOOP ON NLOOP
   ic = 0;
   nrowi = nrow;
   if(iloop == nloop-1) nrowi = nrlast;
   bstride = nrowi*locCols;
   cstride = nrowi*locCols;
   if(runinfo->krow > 0){
     nio = fread(micscoefs->B,sizeof(DOUBLE),2*nrowi*locCols*runinfo->nsymm,fbcoef);
     nio = fread(micscoefs->C,sizeof(DOUBLE),3*nrowi*locCols*runinfo->nsymm,fccoef);
   }
   for(irow=0;irow<nrowi;irow++){    // LOOP ON NROW

    i= j21(irow ,iloop,nrow);
    phiscmval = COMPLEX_ZERO;

    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 * 3*nrowi*locCols;
                        ielmb = ielmu + is*mpi_nelmu;
                        iphiu = ielmu;

                        ic = irow*locCols +j;

                        /* Coefficients at microphones */
                        bb    = micscoefs->B[ic             + bsymmoffset];
                        cc    = micscoefs->C[ic             + csymmoffset];
                        dd    = micscoefs->C[ic + cstride   + csymmoffset];
                        theta = micscoefs->C[ic + cstride*2 + csymmoffset];

                        cexpo  = CEXP(-cfreq*theta);

                        psincval = psincvec[iphiu];           // Psi_incident   on surface
                        phincval = phincvec[iphiu];           // Phi_incident   on surface

                        vec_copy(&un,mpi_elements[ielmb].n);
                        lambda = bc.lambda[ielmu];
                        gamma  = bc.gamma[ielmu];
                        func   = bc.func[ielmu];
                        gamma_over_lambda = gamma/lambda;
                        lambda_over_gamma = lambda/gamma;
                        g.x    = bc.g.x;
                        g.y    = bc.g.y;
                        g.z    = bc.g.z;
                        gdotn  = vec_dot(g,un);

                        if (1 == modsolinfo->dirneu){          // Dirichlet-type formulation

                                if(1 == modsolinfo->knw){              // Phi scatttering on surface
                                        psitcval = solution->sol[iphiu];
                                        psiscval = solution->sol[iphiu] - psincval;                                             // Form A
                                        phitcval = (func+gdotn)/gamma - lambda_over_gamma * psitcval;
                                        phiscmval += bb*cexpo*psitcval + (cc + cfreq*dd)*phitcval*cexpo;
                                }else if (2 == modsolinfo->knw){
                                        psiscval = solution->sol[iphiu];
                                        psitcval = psiscval + psincval;                                                                         // Form B
                                        phitcval = (func+gdotn)/gamma - lambda_over_gamma * psitcval;
                                        phiscval = phitcval - phincval;
                                        phiscmval += bb*cexpo*psiscval + (cc + cfreq*dd)*phiscval*cexpo;
                                }

                        }else if (2 == modsolinfo->dirneu){

                                if(1 == modsolinfo->knw){              // Phi scatttering on surface
                                        phitcval = solution->sol[iphiu];
                                        phiscval = solution->sol[iphiu] - phincval; // Form A
                                }else if (2 == modsolinfo->knw){
                                        phiscval = solution->sol[iphiu];
                                        phitcval = phiscval + phincval;  // Form B
                                }

                                if(1 == modsolinfo->knw){
                                  phiscmval += bb*cexpo*((func+gdotn)/lambda - gamma_over_lambda * phitcval) + (cc + cfreq*dd)*phitcval*cexpo;
                                }else if (2 == modsolinfo->knw){
                                  phiscmval += bb*cexpo*((func+gdotn)/lambda - gamma_over_lambda * phincval - gamma_over_lambda*phiscval - psincval) + (cc+cfreq*dd)*phiscval*cexpo;
                                }
                        }
        }
    }

    phincmval = solution->vecPhincm[i];
    phitmval  = phincmval + phiscmval;

    //    prescmval = -rho_s * phiscmval;
    //    pretmval  = -rho_s * phitmval;

    prescmval = phiscmval;
    pretmval  = phitmval;
    
    solution->vecPhium [i] = phiscmval;

    solution->vecPrescm[i] = prescmval;
    solution->vecPretm [i] = pretmval;

   } // END OF LOOP ON NROW
  } // END OF LOOP ON NLOOP
  if(runinfo->krow > 0){
    fclose(fbcoef); fclose(fccoef);
  }
  Czgsum2d(runinfo->ctxt,"R"," ",locRows,1,solution->vecPhium  ,locRows,-1,-1);
  Czgsum2d(runinfo->ctxt,"R"," ",locRows,1,solution->vecPrescm ,locRows,-1,-1);
  Czgsum2d(runinfo->ctxt,"R"," ",locRows,1,solution->vecPretm  ,locRows,-1,-1);


   // summing on process columns

#ifdef _ACOUSTO_DEBUG
   sprintf(fname,"pre.%d%d",runinfo->myrow,runinfo->mycol);
   f = fopen(fname,"w");
   for(i=0;i<locRows;i++){
     sc_l2g(i,runinfo->myrow,gRows,runinfo->nprows,runinfo->mics_block_size,&imics);
     fprintf(f," %d %f %f %f %f %f %f\n",imics,
                            CREAL(solution->vecPhium[i] ),CIMAG(solution->vecPhium[i]),
                            CREAL(solution->vecPretm[i] ),CIMAG(solution->vecPretm[i]),
                            CREAL(solution->vecPrescm[i]),CIMAG(solution->vecPrescm[i])
                            );
   }
   fclose(f);
#endif

  /* Printing pressures */

  if(modsolinfo->printout){
          print_premic_on_master(cfreq);
  }

  VecDestroy(psincvec);
  VecDestroy(phincvec);

  if(runinfo->krow > 0){
    free(micscoefs->B);
    free(micscoefs->C);
  }

}