pre_mic_coef()

void pre_mic_coef	(	COMPLEX	cfreq,
		int	icounter
	)

Evaluates solution at microphones. Coefficients are calculated. This function is called when pre_calculate_coefs == 0.

Parameters

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

                                             {
  int i,icntr;
  int iphiu;
  int nphiu;
  int nphib;

  int mpi_ncntr,mpi_nelmu;

  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;
  DOUBLE gdotn;


  COMPLEX cexpo;
  DOUBLE bb,cc,dd,theta,thetad;
  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;

  int is;
  struct panel4* mpi_elements;
  struct vector* mpi_cntr;
  int kode;
  DOUBLE source,doublet,ratelet;

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

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

//  char fname2[25];
//  FILE* f2;
//  sprintf(fname2,"coemicdump_dipaolo.%d%d",runinfo->myrow,runinfo->mycol);
//  f2 = fopen(fname2,"w");

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

  gRows      = modgeominfo->nmics;
//  gCols      = modgeominfo->nelmb;
  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);

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


   for(i=0;i<locRows;i++){
    phiscmval = COMPLEX_ZERO;
    sc_l2g(i,runinfo->myrow,gRows,runinfo->nprows,runinfo->mics_block_size,&icntr);

    for(j=0;j<locCols;j++){
          // getting local indices
                sc_l2g(j,runinfo->mycol,gCols,runinfo->npcols,runinfo->col_block_size,&iphiu);

            for(is=0;is<runinfo->nsymm;is++){

                        ielmu = iphiu;
                        ielmb = ielmu + is*mpi_nelmu;

                        /* Coefficients at microphones */
                        kode = 0; // singular panels flag not applicable here
                        delaya(&thetad,ielmb,&mpi_elements[ielmb],mpi_cntr[icntr]);                              // evaluating acoustic delay
                        intgba(1, 0, kode, &mpi_cntr[icntr], &mpi_elements[ielmb], &source, &doublet, &ratelet); // integrating

                        bb    = source;
                        cc    = doublet;
                        dd    = ratelet/runinfo->vsound;
                        theta = thetad/runinfo->vsound;

                        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;
                                }
                        }

            }  // END OF SYMMETRY LOOP

    }

    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 LocRows

  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);

//  fclose(f2);

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

}