Athena: output_pdf.c Source File

Go to the documentation of this file.
00001 #include "copyright.h"
00002 /*============================================================================*/
00003 /*! \file output_pdf.c
00004  *  \brief Outputs Probability Distribution Functions of selected variables
00005  *   in formatted tabular form.  
00006  *
00007  * PURPOSE: Outputs Probability Distribution Functions of selected variables
00008  *   in formatted tabular form.  Fully MPI enabled, which requires passing
00009  *   lots of global sums and means (only the parent process produces output).
00010  *   With SMR, dumps are made for all levels and domains, unless nlevel and
00011  *   ndomain are specified in <output> block.
00012 
00013  *
00014  * CONTAINS PUBLIC FUNCTIONS: 
00015  * - output_pdf() - output PDFs
00016  *============================================================================*/
00017 
00018 #include <math.h>
00019 #include <stdio.h>
00020 #include <stdlib.h>
00021 #include "defs.h"
00022 #include "athena.h"
00023 #include "prototypes.h"
00024 
00025 static int size_dat=0; /* Number of elements in the data[] array */
00026 static double *data=NULL; /* Computed data array: data[size_dat] */
00027 
00028 static int size_pdf=0; /* Number of elements in the pdf[] array */
00029 static int *pdf=NULL; /* (non-normalized) PDF */
00030 #ifdef MPI_PARALLEL
00031 static int *cd_pdf=NULL; /* (non-normalized) complete Domain PDF */
00032 #endif /* MPI_PARALLEL */
00033 
00034 static char def_fmt[]="%21.15e"; /* A default tabular dump data format */
00035 
00036 /*----------------------------------------------------------------------------*/
00037 /*! \fn void output_pdf(MeshS *pM, OutputS *pOut)
00038  *  \brief Outputs PDFs. */
00039 
00040 void output_pdf(MeshS *pM, OutputS *pOut)
00041 {
00042   GridS *pG;
00043   FILE *pfile;
00044   char fmt[80];
00045   char fid[80]; /* File "id" for the statistics table */
00046   char *fname,*plev=NULL,*pdom=NULL,*pdir=NULL;
00047   char levstr[8],domstr[8],dirstr[8];
00048   int nl,nd,i,j,k,is,ie,js,je,ks,ke;
00049   int n, data_cnt;
00050   double dmin, dmax, delta, dpdf, dat, scl;
00051   double mean=0.0, var=0.0; /* mean and variance of the distribution */
00052   double adev=0.0, sdev=0.0; /* average & standard deviation */
00053   double skew=0.0, kurt=0.0; /* skewness and kurtosis of the distribution */
00054   double r, s, ep=0.0; /* Temp. variables for calculating the variance, etc. */
00055 #ifdef MPI_PARALLEL
00056   DomainS *pD;
00057   int ierr, cd_data_cnt, myID_Comm_Domain;
00058 #endif /* MPI_PARALLEL */
00059 
00060 /* Loop over all Domains in Mesh, and output Grid data */
00061 
00062   for (nl=0; nl<(pM->NLevels); nl++){
00063     for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
00064       if (pM->Domain[nl][nd].Grid != NULL){
00065 
00066 /* write files if domain and level match input, or are not specified (-1) */
00067       if ((pOut->nlevel == -1 || pOut->nlevel == nl) &&
00068           (pOut->ndomain == -1 || pOut->ndomain == nd)){
00069         pG = pM->Domain[nl][nd].Grid;
00070 
00071 #ifdef MPI_PARALLEL
00072         pD = (DomainS*)&(pM->Domain[nl][nd]);
00073         ierr = MPI_Comm_rank(pD->Comm_Domain, &myID_Comm_Domain);
00074         cd_data_cnt = (pD->Nx[0])*(pD->Nx[1])*(pD->Nx[2]);
00075 #endif
00076         is = pG->is, ie = pG->ie;
00077         js = pG->js, je = pG->je;
00078         ks = pG->ks, ke = pG->ke;
00079         data_cnt = (ie - is + 1)*(je - js + 1)*(ke - ks + 1);
00080 
00081 /* Are the requisite arrays allocated? */
00082         if(data == NULL){
00083           size_dat = data_cnt;
00084           data = (double *)calloc(size_dat,sizeof(double));
00085           if(data == NULL)
00086             ath_error("[output_pdf]: Failed to allocate data array\n");
00087 
00088 /* This choice for size_pdf represents a balance between
00089  * resolution in the PDF and "shot noise" in the data binning. */
00090 #ifdef MPI_PARALLEL
00091           size_pdf = (int)sqrt((double)cd_data_cnt);
00092 #else /* MPI_PARALLEL */
00093           size_pdf = (int)sqrt((double)size_dat);
00094 #endif /* MPI_PARALLEL */
00095           pdf = (int *)calloc(size_pdf,sizeof(int));
00096           if(pdf == NULL)
00097             ath_error("[output_pdf]: Failed to allocate pdf array\n");
00098 
00099 #ifdef MPI_PARALLEL
00100           if(myID_Comm_Domain == 0){ /* I'm the parent */
00101             cd_pdf = (int *)calloc(size_pdf,sizeof(int));
00102             if(cd_pdf == NULL)
00103               ath_error("[output_pdf]: Failed to allocate cd_pdf array\n");
00104           }
00105 #endif /* MPI_PARALLEL */
00106         }
00107 
00108 
00109 /* Initialize dmin, dmax */
00110         dmin = dmax = (*pOut->expr)(pG,is,js,ks);
00111 
00112 /* Fill the data array */
00113         n=0;
00114         for(k = ks; k<=ke; k++){
00115           for(j = js; j<=je; j++){
00116             for(i = is; i<=ie; i++){
00117               data[n] = (double)(*pOut->expr)(pG,i,j,k);
00118               dmin = data[n] < dmin ? data[n] : dmin;
00119               dmax = data[n] > dmax ? data[n] : dmax;
00120               mean += data[n];
00121               n++;
00122             }
00123           }
00124         }
00125 
00126 #ifdef MPI_PARALLEL
00127 
00128         dat = dmin;
00129         ierr = MPI_Allreduce(&dat,&dmin,1,MPI_DOUBLE,MPI_MIN,pD->Comm_Domain);
00130 
00131         dat = dmax;
00132         ierr = MPI_Allreduce(&dat,&dmax,1,MPI_DOUBLE,MPI_MAX,pD->Comm_Domain);
00133 
00134         dat = mean;
00135         ierr = MPI_Allreduce(&dat,&mean,1,MPI_DOUBLE,MPI_SUM,pD->Comm_Domain);
00136 
00137         mean /= (double)cd_data_cnt; /* Complete the calc. of the mean */
00138 
00139 #else /* MPI_PARALLEL */
00140 
00141         mean /= (double)data_cnt; /* Complete the calc. of the mean */
00142 
00143 #endif /* MPI_PARALLEL */
00144 
00145         if(data_cnt > 1){
00146 /* Calculate the variance, etc. with the corrected 2-pass formula */
00147           for(n=0; n<data_cnt; n++){
00148             s = data[n] - mean;
00149             adev += fabs(s);
00150             ep += s;
00151             var += (r = s*s);
00152             skew += (r *= s);
00153             kurt += (r *= s);
00154           }
00155 
00156 #ifdef MPI_PARALLEL
00157           dat = ep;
00158           ierr = MPI_Allreduce(&dat,&ep,1,MPI_DOUBLE,MPI_SUM,pD->Comm_Domain);
00159 
00160           dat = var;
00161           ierr = MPI_Allreduce(&dat,&var,1,MPI_DOUBLE,MPI_SUM,pD->Comm_Domain);
00162 
00163           dat = skew;
00164           ierr = MPI_Allreduce(&dat,&skew,1,MPI_DOUBLE,MPI_SUM,pD->Comm_Domain);
00165 
00166           dat = kurt;
00167           ierr = MPI_Allreduce(&dat,&kurt,1,MPI_DOUBLE,MPI_SUM,pD->Comm_Domain);
00168 
00169           adev /= (double)cd_data_cnt;
00170           var = (var - ep*ep/(double)cd_data_cnt)/(double)(cd_data_cnt-1);
00171           sdev = sqrt(var);
00172           if(sdev > 0.0){
00173             skew /= var*sdev*cd_data_cnt;
00174             kurt = kurt/(var*var*cd_data_cnt) - 3.0;
00175           }
00176 #else /* MPI_PARALLEL */
00177           adev /= (double)data_cnt;
00178           var = (var - ep*ep/(double)data_cnt)/(double)(data_cnt-1);
00179           sdev = sqrt(var);
00180           if(sdev > 0.0){
00181             skew /= var*sdev*data_cnt;
00182             kurt = kurt/(var*var*data_cnt) - 3.0;
00183           }
00184 #endif /* MPI_PARALLEL */
00185         }
00186 
00187 /* Store the global maximum and minimum of the quantity */
00188         pOut->gmin = dmin < pOut->gmin ? dmin : pOut->gmin;
00189         pOut->gmax = dmax > pOut->gmax ? dmax : pOut->gmax;
00190 
00191 /* Compute the pdf directly using sampling. Define size_pdf bins, each of equal
00192  * size, and fill them with the number of cells whose data value falls in the
00193  * range spanned by the bin. */
00194 
00195         if(dmax - dmin > 0.0){
00196 /* Initialize pdf[] to zero */
00197           for(n=0; n<size_pdf; n++) pdf[n] = 0;
00198 /* Calculate the number of cells whose data falls in each bin */
00199           scl = (double)size_pdf/(dmax - dmin);
00200           for(n=0; n<data_cnt; n++){
00201             i = (int)(scl*(data[n] - dmin));
00202             i = i < size_pdf ? i : size_pdf - 1;
00203             pdf[i]++;
00204           }
00205         }
00206 
00207 #ifdef MPI_PARALLEL
00208 
00209 /* Sum up the pdf in the array cd_pdf */
00210         ierr=MPI_Reduce(pdf,cd_pdf,size_pdf,MPI_INT,MPI_SUM,0,pD->Comm_Domain);
00211 
00212 #endif /* MPI_PARALLEL */
00213 
00214 #ifdef MPI_PARALLEL
00215 /* For parallel calculations, only the parent writes the output. */
00216         if(myID_Comm_Domain != 0) return;
00217 #endif /* MPI_PARALLEL */
00218 
00219 /* Create filename and open file.  pdf files are always written in lev#
00220  * directories of root (rank=0) process. */
00221 #ifdef MPI_PARALLEL
00222         if (nl>0) {
00223           plev = &levstr[0];
00224           sprintf(plev,"lev%d",nl);
00225           pdir = &dirstr[0];
00226           sprintf(pdir,"../lev%d",nl);
00227         }
00228 #else
00229         if (nl>0) {
00230           plev = &levstr[0];
00231           sprintf(plev,"lev%d",nl);
00232           pdir = &dirstr[0];
00233           sprintf(pdir,"lev%d",nl);
00234         }
00235 #endif
00236         if (nd>0) {
00237           pdom = &domstr[0];
00238           sprintf(pdom,"dom%d",nd);
00239         }
00240 
00241         fname = ath_fname(pdir,pM->outfilename,plev,pdom,num_digit,
00242           pOut->num,pOut->id,"prb");
00243         if(fname == NULL){
00244           ath_perr(-1,"[output_pdf]: Unable to create filename\n");
00245         }
00246         pfile = fopen(fname,"w");
00247         if(pfile == NULL){
00248           ath_perr(-1,"[output_pdf]: Unable to open pdf file\n");
00249         }
00250 
00251 /* Write out some extra information in a header */
00252         fprintf(pfile,"# Time = %21.15e\n",pG->time);
00253         fprintf(pfile,"# expr = \"%s\"\n",pOut->out);
00254         fprintf(pfile,"# Nbin = %d\n",((dmax - dmin) > 0.0 ? size_pdf : 1));
00255         fprintf(pfile,"# dmin = %21.15e\n",dmin); 
00256         fprintf(pfile,"# dmax = %21.15e\n",dmax); 
00257         fprintf(pfile,"# mean = %21.15e\n",mean); 
00258         fprintf(pfile,"# variance = %21.15e\n",var); 
00259         fprintf(pfile,"# std. dev. = %21.15e\n",sdev); 
00260         fprintf(pfile,"# avg. dev. = %21.15e\n",adev); 
00261         fprintf(pfile,"# skewness = %21.15e\n",skew); 
00262         fprintf(pfile,"# kurtosis = %21.15e\n#\n",kurt); 
00263 
00264 /* Add a white space to the format */
00265         if(pOut->dat_fmt == NULL)
00266           sprintf(fmt,"%s  %s\n",def_fmt, def_fmt);
00267         else
00268           sprintf(fmt,"%s  %s\n",pOut->dat_fmt,pOut->dat_fmt);
00269 
00270 /* write out the normalized Proabability Distribution Function */
00271         if(dmax - dmin > 0.0){
00272           delta = (dmax - dmin)/(double)(size_pdf);
00273 #ifdef MPI_PARALLEL
00274           scl = (double)size_pdf/(double)(cd_data_cnt*(dmax - dmin));
00275 #else
00276           scl = (double)size_pdf/(double)(data_cnt*(dmax - dmin));
00277 #endif /* MPI_PARALLEL */
00278           for(n=0; n<size_pdf; n++){
00279 /* Calculate the normalized Prob. Dist. Fun. */
00280             dat = dmin + (n + 0.5)*delta;
00281 #ifdef MPI_PARALLEL
00282             dpdf = (double)(cd_pdf[n])*scl;
00283 #else
00284             dpdf = (double)(pdf[n])*scl;
00285 #endif /* MPI_PARALLEL */
00286             fprintf(pfile, fmt, dat, dpdf);
00287           }
00288         }
00289         else
00290           fprintf(pfile,fmt,dmax,1.0);
00291 
00292         fclose(pfile);
00293 
00294 /* Also write a history type file on the statistics */
00295         sprintf(fid,"prb_stat.%s",pOut->id);
00296 
00297         fname = ath_fname(pdir,pM->outfilename,plev,pdom,0,0,fid,"tab");
00298         if(fname == NULL){
00299           ath_perr(-1,"[output_pdf]: Unable to create stats filename\n");
00300         }
00301         pfile = fopen(fname,"a");
00302         if(pfile == NULL){
00303           ath_perr(-1,"[output_pdf]: Unable to open stats file\n");
00304         }
00305 
00306         if(pOut->num == 0){
00307           fprintf(pfile,"# expr = \"%s\"\n#\n",pOut->out);
00308           fprintf(pfile,"# time  dmin  dmax  mean  variance  \"std. dev.\"  ");
00309           fprintf(pfile,"\"avg. dev.\"  skewness  kurtosis\n#\n");
00310         }
00311 
00312 /* Add a white space to the format */
00313         if(pOut->dat_fmt == NULL) sprintf(fmt," %s",def_fmt);
00314         else                      sprintf(fmt," %s",pOut->dat_fmt);
00315 
00316         fprintf(pfile,"%21.15e",pG->time);
00317 
00318 /* write out the table of statistics */
00319         fprintf(pfile,fmt,dmin);
00320         fprintf(pfile,fmt,dmax);
00321         fprintf(pfile,fmt,mean);
00322         fprintf(pfile,fmt,var);
00323         fprintf(pfile,fmt,sdev);
00324         fprintf(pfile,fmt,adev);
00325         fprintf(pfile,fmt,skew);
00326         fprintf(pfile,fmt,kurt);
00327         fprintf(pfile,"\n");
00328 
00329         fclose(pfile);
00330       }}
00331     }
00332   }
00333 
00334   return;
00335 }