matyfdc()

void matyfdc	(	COMPLEX	cfreq,
		int	icounter
	)

Evaluates System Matrix and RHS of the linear System. Valid for Neumann-type & Dirichlet-type Boundary Conditions. Coefficients are evaluated here and discarded after use.

Parameters

cfreq	COMPLEX frequency
icounter	index of frequency

                                        {

  int icntr;
  int ielmb,ielmu;
  int nphiu, mpi_nelmu;
  int nphib;

  int mpi_ncntr;
  COMPLEX cexpc;
  DOUBLE bb,cc,dd,thetac,estar;
  COMPLEX matval = 0.0 + I*0.0;;
  COMPLEX yy = 0.0 + I*0.0;;
  int cstride;
  int kode;

  //DOUBLE *C;

  int gRows,gCols;     // global dim of coefficient matrix
  int locRows,locCols; // local dim of coefficient matrix
  int izero = 0;
  int i,j;
  COMPLEX lambda;
  COMPLEX gamma;
  COMPLEX func;
  struct vector g;
  struct vector un;
  COMPLEX gamma_over_lambda;
  COMPLEX lambda_over_gamma;
  DOUBLE gdotn;

#ifdef _ACOUSTO_DEBUG
  FILE* f;
  FILE* cmatfile;
  char fname[15];
  int ipsiu;
#endif

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


 //  int nrow,nrowi,iloop,nloop,irow,nio,nrlast;

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

   MPI_Group* mpi_row_groups;
   MPI_Group mpi_group_world;
   MPI_Comm* mpi_comm;
   int* row_ranks;
   COMPLEX* ctmp;
   Vec psiuvec;
   Vec psincvec;
   Vec phincvec;
   int gCRows;
   int nn;
   int mm;
   int descpsiu[9];
   int descgpsiu[9];
   int ione = 1;
   int info;


  // number of control points,phib and phiu points
  mpi_elements = geometry->elements;
  mpi_cntr     = geometry->cntrc;
  mpi_ncntr  = modgeominfo->ncntr + modgeominfo->nchief;
  nphiu      = modgeominfo->ncntr;
  nphib      = modgeominfo->nelmb;
  mpi_nelmu  = nphiu;

  gRows  = mpi_ncntr;
  gCols  = nphiu;


#ifdef _ACOUSTO_DEBUG
  cmatfile = fopen("CD.matrix.out","w");
  fprintf(cmatfile,"CD=[");
#endif

  locRows  = numroc_(&gRows ,&runinfo->row_block_size,&runinfo->myrow,&izero,&runinfo->nprows);
  locCols  = numroc_(&gCols ,&runinfo->col_block_size,&runinfo->mycol,&izero,&runinfo->npcols);
  cstride = locRows * locCols; // C stride

  // nrow partition on node to limit memory usage (see manual)
  if (solution->matY == NULL) {
    solution->matY = calloc(locRows*locCols,sizeof(COMPLEX)); // allocation
  } else {
    for(i = 0; i < locRows*locCols; i++)
      solution->matY[i] = COMPLEX_ZERO;
  }

  gCRows   = nphiu;
  psiuvec  = calloc(nphiu,sizeof(COMPLEX)); // allocating vector
  psincvec = calloc(nphiu,sizeof(COMPLEX)); // allocating vector
  phincvec = calloc(nphiu,sizeof(COMPLEX)); // allocating vector
  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->vecPsiu ,&ione, &ione,descpsiu, psiuvec ,&ione,&ione,descgpsiu,&runinfo->ctxt);
  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(psiuvec ,gCRows,MPI_ACOUSTO_COMPLEX,0,MPI_COMM_WORLD);
  MPI_Bcast(psincvec,gCRows,MPI_ACOUSTO_COMPLEX,0,MPI_COMM_WORLD);
  MPI_Bcast(phincvec,gCRows,MPI_ACOUSTO_COMPLEX,0,MPI_COMM_WORLD);

  CREATE_VEC(solution->rhs,locRows); // allocation
  VecZeroEntries(solution->rhs,locRows); // reset

  CREATE_VEC(ctmp,locRows); // allocation

   for(i=0;i<locRows;i++){    // LOOP ON ICNTR
    sc_l2g(i,runinfo->myrow,gRows,runinfo->nprows,runinfo->row_block_size,&icntr);

    if(modsolinfo->knw == 1){
      if(runinfo->mycol == 0){
          if (modsolinfo->dirneu == 1) solution->rhs[i] =  -phincvec[icntr]; // DIRICHLET-TYPE
          if (modsolinfo->dirneu == 2) solution->rhs[i] =   phincvec[icntr]; // NEUMANN-TYPE
      }else{
        solution->rhs[i] = COMPLEX_ZERO;
      }
    }else{
      solution->rhs[i] = COMPLEX_ZERO;
    }

    for(j=0;j<locCols;j++){ // LOOP ON NELMB
        sc_l2g(j,runinfo->mycol,gCols,runinfo->npcols,runinfo->col_block_size,&ielmu);
                //ic = irow*locCols +j;

                estar = (icntr == ielmu ) ? 0.5 : 0.0;
                solution->matY[locRows*j + i ] += estar;

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

                ielmb = ielmu + is*mpi_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

                        bb     = source;
                        cc     = doublet;
                        dd     = ((1 == kode) ? (DOUBLE)0.0 : (DOUBLE)ratelet/runinfo->vsound );
                        thetac = theta/runinfo->vsound;
                        cexpc  = CEXP(-cfreq*thetac);

                        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);
                        matval =  (cc+cfreq*dd) * cexpc ;
                        //
                        // THIS IS CLASSICAL YMAT=E*I-(C+s*D). NEEDED TO BUILD NEUMANN & DIRICHLET RHSs.
                        // NOTE: matval includes delay
                        //
                        solution->matY[locRows*j + i ] -= matval;
                        yy     = solution->matY[locRows*j + i ];

                        if(modsolinfo->knw == 1){
                                if (modsolinfo->dirneu == 1) solution->rhs[i] += yy * (func+gdotn)/gamma;                       // DIRICHLET
                                if (modsolinfo->dirneu == 2) solution->rhs[i] += bb * cexpc * (func+gdotn)/lambda;      // NEUMANN
                        }else if(modsolinfo->knw == 2){
                                if (modsolinfo->dirneu == 1) solution->rhs[i] += yy * ((func+gdotn)/gamma - psincvec[ielmu]*lambda_over_gamma - phincvec[ielmu]);                       // DIRICHLET
                                if (modsolinfo->dirneu == 2) solution->rhs[i] += bb * cexpc*((func+gdotn)/lambda - phincvec[ielmu]*gamma_over_lambda - psincvec[ielmu]);        // NEUMANN
                        }
                        //
                        // HERE SYSTEM MATRIX IS BUILT ON YMAT.
                        //
                        if (modsolinfo->dirneu == 1) solution->matY[locRows*j + i ]  = bb*cexpc + yy * lambda_over_gamma; //DIRICHLET (matrix entry is rewritten)
                        if (modsolinfo->dirneu == 2) solution->matY[locRows*j + i ] += bb * gamma_over_lambda *cexpc; //NEUMANN   (matrix entry is updated)

#ifdef _ACOUSTO_DEBUG
        fprintf(cmatfile,"%f + %%i*%f ,",CREAL(solution->matY[locRows*j + i ]),CIMAG(solution->matY[locRows*j + i ]));;
#endif

      } // END OF LOOP ON NSYMM
//          fprintf(f2,"%d %d %f %f %f %f\n",icntr,ielmu,CREAL(bbtmp),CIMAG(bbtmp),CREAL(ytmp),CIMAG(ytmp));
    }




#ifdef _ACOUSTO_DEBUG
        fprintf(cmatfile,";\n");
#endif
   } // END OF LOOP ON NROW
  /*} // END OF LOOP ON NLOOP */

//  fclose(f2);

#ifdef _ACOUSTO_DEBUG
        fprintf(cmatfile,"];\n");
#endif

   
  // Reducing RHS if cols in proc grid are > 1
  if(runinfo->npcols > 1){
   logger(LOG_INFO," Reducing RHS on columns\n");
   mpi_row_groups = calloc(sizeof(MPI_Group) , runinfo->nprows);
   mpi_comm = calloc(sizeof(MPI_Comm) , runinfo->nprows);
   row_ranks = calloc(sizeof(int) , runinfo->npcols);
   MPI_Comm_group(MPI_COMM_WORLD, &mpi_group_world);
   for(i=0;i<runinfo->nprows;i++){
      for(j=0;j<runinfo->npcols;j++){
        row_ranks[j] = j + i*runinfo->npcols;
      }
      MPI_Group_incl(mpi_group_world,runinfo->npcols,row_ranks,&mpi_row_groups[i]);
      MPI_Comm_create(MPI_COMM_WORLD, mpi_row_groups[i], &mpi_comm[i]);
   }
 
   for(i=0;i<runinfo->nprows;i++){
 
     if(MPI_COMM_NULL != mpi_comm[i]){
      MPI_Allreduce((void*)solution->rhs,
                    (void*)&ctmp[0],
                    locRows,
                    mpi_complex,
                    mpi_op_complex_sum,
                    mpi_comm[i]);
     }
   }
 
 
   for(i=0;i<locRows;i++){
     /*
     printf("(%d) %d %f %f - %f %f\n",rank,i,
            CREAL(solution->rhs[i]),CIMAG(solution->rhs[i]),
            CREAL(ctmp[i]),CIMAG(ctmp[i]));
     */
 
     solution->rhs[i] = ctmp[i];
   }
  }


  if(psiuvec)  free(psiuvec);
  if(psincvec) free(psincvec);
  if(phincvec) free(phincvec);
  if(ctmp)     free(ctmp);
#ifdef _ACOUSTO_DEBUG
FILE* f;
char fname[64];
int ipsiu;
  sprintf(fname,"maty.%d%d",runinfo->myrow,runinfo->mycol);
  f = fopen(fname,"w");
  for(i=0;i<locRows;i++){    // LOOP ON ROWS
    sc_l2g(i,runinfo->myrow,gRows,runinfo->nprows,runinfo->row_block_size,&icntr);
    for(j=0;j<locCols;j++){ // LOOP on COLS
      sc_l2g(j,runinfo->mycol,gCols,runinfo->npcols,runinfo->col_block_size,&ipsiu);
      fprintf(f,"%d %d %f %f\n",icntr,ipsiu,CREAL(solution->matY[i+locRows*j]), CIMAG(solution->matY[i+locRows*j]));
    }
  }
  fclose(f);




  sprintf(fname,"rhs.%d.out",modgeominfo->nelmb);  
  f = fopen(fname,"w");
  for(i=0;i<locRows;i++){    // LOOP ON ROWS
    sc_l2g(i,runinfo->myrow,gRows,runinfo->nprows,runinfo->row_block_size,&icntr);
    fprintf(f,"%f %f %f %f %f %f\n",
                       geometry->cntr[icntr].x,
                       geometry->cntr[icntr].y,
                       geometry->cntr[icntr].z,
                       CREAL(solution->rhs[i]),
                       CIMAG(solution->rhs[i]),
                       CABS(solution->rhs[i]));
  }
  fclose(f);

#endif
#undef _ACOUSTO_DEBUG

#ifdef _ACOUSTO_DEBUG
FILE* f;
char fname[64];
  sprintf(fname,"rhs.%d.out",modgeominfo->nelmb);  
  f = fopen(fname,"w");
  for(i=0;i<locRows;i++){    // LOOP ON ROWS
    sc_l2g(i,runinfo->myrow,gRows,runinfo->nprows,runinfo->row_block_size,&icntr);
    fprintf(f,"%f %f %f %f %f %f\n",
                       geometry->cntr[icntr].x,
                       geometry->cntr[icntr].y,
                       geometry->cntr[icntr].z,
                       CREAL(solution->rhs[i]),
                       CIMAG(solution->rhs[i]),
                       CABS(solution->rhs[i]));
  }
  fclose(f);

#endif
#undef _ACOUSTO_DEBUG
}