acoustic_solution()

void acoustic_solution

(

)

Acoustic Solution main function. The solution is evaluated and output files are printed.

                        {
    
   int i,j;
   int isig,iome,ifreq;
   DOUBLE delome,delsig,sigma,omega;
   int nfreq;
   int nphiu,npsiu;
   int mpi_ncntr;

   COMPLEX cfreq,hfreq;

   //DOUBLE norm, cond;


   int nelmb,ncntrc;
   int icounter;
   //KSP ksp;
   char fname[MAX_PATH];
   char fsolname[MAX_PATH];
   FILE *f;
#ifdef _ACOUSTO_DEBUG
   FILE *file;
   Vec vrhs;
   COMPLEX vrhsval;
   COMPLEX vrhsvalc;
   int istart,iend;
#endif

   delome = (modsolinfo->maxome-modsolinfo->minome) / (modsolinfo->nome-1);
   if(1   == modsolinfo->nome) delome = 0.0;
   if(0.0 == delome          ) delome = 0.0;
   delsig = (modsolinfo->maxsig-modsolinfo->minsig) / (modsolinfo->nsig-1);
   if(1   == modsolinfo->nsig) delsig = 0.0;
   if(0.0 == delsig          ) delsig = 0.0;

   nfreq = modsolinfo->nsig * modsolinfo->nome;

   nphiu  = modgeominfo->ncntr;  // MAKE THIS COHERENT 
   npsiu  = modgeominfo->ncntr;
   nelmb  = modgeominfo->nelmb;
   ncntrc = modgeominfo->ncntr + modgeominfo->nchief;

   // Allocation
   //solution =  calloc(1,sizeof(struct solution_struct));
   solution->sig     = calloc(modsolinfo->nsig,     sizeof(DOUBLE));
   solution->ome     = calloc(modsolinfo->nome,     sizeof(DOUBLE));
   solution->sources = calloc(modsolinfo->nsourc,   sizeof(struct acoustic_point));
   solution->dipoles = calloc(modsolinfo->ndipoles, sizeof(struct dipole));

   // Custom BC Allocation
    bc.lambda  = calloc(nphiu,sizeof(COMPLEX));
    for(j=0;j<nphiu;j++) bc.lambda[j] = 1.0;

    bc.gamma   = calloc(nphiu,sizeof(COMPLEX));
    for(j=0;j<nphiu;j++) bc.gamma[j]  = COMPLEX_ZERO;

    bc.func    = calloc(nphiu,sizeof(COMPLEX));
    for(j=0;j<nphiu;j++) bc.func[j] = COMPLEX_ZERO;

    bc.g.x    = 0.0;
    bc.g.y    = 0.0;
    bc.g.z    = 0.0;

   solution->planw   = calloc(modsolinfo->nplanw,  sizeof(struct acoustic_point));
                                                        
//   if(modsolinfo->nsig == 1){
//     solution->sig[0] = modsolinfo->minsig;
//   }else{
//    for(i=0;i<modsolinfo->nsig;i++){
//      solution->sig[i] = modsolinfo->minsig + i*delsig;
//    }
//   }
//   if(modsolinfo->nome == 1){
//     solution->ome[0] = modsolinfo->minome;
//   }else{
//     for(i=0;i<modsolinfo->nome;i++){
//       solution->ome[i] = modsolinfo->minome + i*delome;
//     }
//   }
  
   // Read data from files
   // Sources and plane wave.
   if(rank == 0){
    if(!read_acoustic_files()){
        logger(LOG_ERROR,"  **************************** \n");
        logger(LOG_ERROR,"  Error reading acoustic files \n");
        logger(LOG_ERROR,"  **************************** \n");
    }
   }


   // allocating indices on slave nodes for freq resp at microphones
   if(0 != rank){
     solution->i_freqresp_mic = calloc(modsolinfo->n_freqresp_mic,sizeof(int));
   }

   // Broadcasting files read on node 0 to all nodes
   MPI_Bcast(solution->sources,modsolinfo->nsourc   ,mpi_acoustic_point_type,0,MPI_COMM_WORLD);
   MPI_Bcast(solution->dipoles,modsolinfo->ndipoles ,mpi_dipole_type,0,MPI_COMM_WORLD);
   MPI_Bcast(solution->planw  ,modsolinfo->nplanw ,mpi_acoustic_point_type,0,MPI_COMM_WORLD);
   MPI_Bcast(solution->i_freqresp_mic  ,modsolinfo->n_freqresp_mic ,MPI_INT,0,MPI_COMM_WORLD);
   

#ifdef _ACOUSTO_DEBUG  
  sprintf(fname,"checkcast.%d",rank);
  logger(LOG_DEBUG,"filename [%s]\n",fname);
  file = fopen(fname,"w");
  fprintf(file,"SOURCES (%d)\n",modsolinfo->nsourc);
  for(i=0;i<modsolinfo->nsourc;i++){
    fprintf(file,f1I5F,
            i,
            solution->sources[i].pos.x,
            solution->sources[i].pos.y,
            solution->sources[i].pos.z,
            CREAL(solution->sources[i].amplitude),
            CIMAG(solution->sources[i].amplitude));
  }
  fprintf(file,"IMPEDANCES (%d)\n",modsolinfo->nimped);
  for(i=0;i<modsolinfo->nimped;i++){
    fprintf(file,f1I2F,solution->izimp[i],CREAL(solution->zimp[i]),CIMAG(solution->zimp[i]));
  }
  fprintf(file,"RADIANTS (%d)\n",modsolinfo->nradian);
  for(i=0;i<modsolinfo->nradian;i++){
    fprintf(file,f1I2F,solution->iradian[i],CREAL(solution->radian[i]),CIMAG(solution->radian[i]));
  }
  fclose(file);
#endif

   mpi_ncntr  = modgeominfo->ncntr + modgeominfo->nchief;
   CREATE_VEC(solution->sol,mpi_ncntr);

// Init frequency response files if needed   
   if( 0 == rank){
       if(0 != modsolinfo->n_freqresp_cnt){
         for(i=0;i<modsolinfo->n_freqresp_cnt;i++){
           get_filename(fname,"%s-fresp-cnt%d.out",rundetails->title,solution->i_freqresp_cnt[i]);  
           f = fopen(fname,"w+");       
           fprintf(f,"# Frequency response at control point %d\n",solution->i_freqresp_cnt[i]);
           fclose(f);
         }
       }
       if(0 != modsolinfo->n_freqresp_mic){
         for(i=0;i<modsolinfo->n_freqresp_mic;i++){
          get_filename(fname,"%s-fresp-mic%d.out",rundetails->title,solution->i_freqresp_mic[i]);  
          f = fopen(fname,"w+");        
          fprintf(f,"# Frequency response at microphone %d\n",solution->i_freqresp_mic[i]);
          fclose(f);
         }
       }
   }

   
   /***********************************
    *          Starting loop 
    ***********************************/
   icounter  =0;
   for(isig=0;isig<modsolinfo->nsig;isig++){
   sigma = modsolinfo->minsig + isig*delsig;
    for(iome=0;iome<modsolinfo->nome;iome++){
      omega = modsolinfo->minome + iome*delome;
//      cfreq = solution->sig[isig] + I*solution->ome[iome];
      cfreq = sigma + I*omega;
      hfreq = cfreq;

      ifreq = iome + modsolinfo->nome*isig;
      logger(LOG_INFO,"ifreq = %d, s=(%f +i %f), f=%f Hz\n",ifreq,(double)CREAL(cfreq),(double)CIMAG(cfreq),(double)(0.5*omega/(PI)));

      // Read and broadcast bc from files
      if(rank == 0){read_boundary_conditions(cfreq, ifreq);}

      logger(LOG_INFO,"Boundary conditions: transferring %d bytes to nodes\n",(3*nphiu*sizeof(COMPLEX)+sizeof(struct vector)));
      MPI_Bcast(&bc.nlambda  ,1 ,MPI_INT,0,MPI_COMM_WORLD);
      MPI_Bcast(&bc.ngamma   ,1 ,MPI_INT,0,MPI_COMM_WORLD);
      MPI_Bcast(&bc.nfunc    ,1 ,MPI_INT,0,MPI_COMM_WORLD);
      MPI_Bcast(&bc.ng       ,1 ,MPI_INT,0,MPI_COMM_WORLD);
      MPI_Bcast(bc.lambda    ,nphiu ,MPI_ACOUSTO_COMPLEX,0,MPI_COMM_WORLD);
      MPI_Bcast(bc.gamma     ,nphiu ,MPI_ACOUSTO_COMPLEX,0,MPI_COMM_WORLD);
      MPI_Bcast(bc.func      ,nphiu ,MPI_ACOUSTO_COMPLEX,0,MPI_COMM_WORLD);
      MPI_Bcast(&bc.g        ,1     ,mpi_vector_type   ,0,MPI_COMM_WORLD);

      logger(LOG_DEBUG,"   Evaluating Phi inc.\n");
      eval_phinc(cfreq,icounter);

      logger(LOG_DEBUG,"   Evaluating Normal Wash\n");
      eval_nw(cfreq,icounter);

      if(modsolinfo->onlyrep == 0){
        // building  matrices and RHS
        if(modsolinfo->pre_calculate_coefs){
          logger(LOG_DEBUG,"   Evaluating Y\n");
          matyfd(cfreq,icounter);
        }else{
          matyfdc(cfreq,icounter);
        }
  
        // Estimating condition number of matrix 
        if(PSEUDOINVERSE == modsolinfo->solver){
          solve_linear_system(icounter);
        }else if (GMRES == modsolinfo->solver){
          gmres(icounter);
        }
  
        if(modsolinfo->printout){
                printsol_on_master(icounter,cfreq);
        }
      }else{ // reading solution
        get_filename(fsolname,"%s-surf-%.4fHz.out",rundetails->title,(double)HZFREQ(cfreq));
        logger(LOG_INFO,"Skipping evaluation, reading solution from file %s. on master node\n",fsolname);
        if(0 == rank){
          read_solution(fsolname);
        }
        // Broadcast
        logger(LOG_INFO,"Solution read on master node. Broadcasting to nodes\n");
        MPI_Bcast(solution->sol,mpi_ncntr,mpi_complex,0,MPI_COMM_WORLD);

      }
      
      if(modsolinfo->pre_calculate_coefs){
        pre_mic(cfreq,icounter); /* evaluating pressure at microphones */
      }else{
        pre_mic_coef(cfreq,icounter); /* evaluating pressure at microphones */
      }

#ifdef _MYSQL_
      // Storing run on MySQL if enabled.
      if(mysqlinfo->active){
       mysql_save(cfreq);
      }
#endif

      // Frequency response analysis printout
      // Centers
      if( 0 == rank){
          if(0 != modsolinfo->n_freqresp_cnt){
            for(i=0;i<modsolinfo->n_freqresp_cnt;i++){
              get_filename(fname,"%s-fresp-cnt%d.out",rundetails->title,solution->i_freqresp_cnt[i]);  
              f = fopen(fname,"a");     
              fr_print_sol_cntr(f,solution->i_freqresp_cnt[i],icounter,cfreq);
              fclose(f);
            }
          }
      }

      // Microphones
      if(0 != modsolinfo->n_freqresp_mic){
        for(i=0;i<modsolinfo->n_freqresp_mic;i++){
          if(0 == rank){
           get_filename(fname,"%s-fresp-mic%d.out",rundetails->title,solution->i_freqresp_mic[i]);
           f = fopen(fname,"a");
          }
          fr_print_pre_mic(f,solution->i_freqresp_mic[i],icounter,cfreq);
          if(0 == rank){
            fclose(f);
          }
        }
      }
      
/* VTK files printout */      
      if(modsolinfo->printvtk){
          vtkoutb(icounter, cfreq);
          if(1==havemicsmesh){
                  vtkoutmmesh(icounter, cfreq);
          }else{
                  vtkoutmu2(icounter, cfreq);
          }
      }
      if(modsolinfo->printmsh){
          mshoutb(icounter, cfreq);
          mshoutm(icounter, cfreq);
      }

      icounter++;
    }
   }

}