00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
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;
00026 static double *data=NULL;
00027
00028 static int size_pdf=0;
00029 static int *pdf=NULL;
00030 #ifdef MPI_PARALLEL
00031 static int *cd_pdf=NULL;
00032 #endif
00033
00034 static char def_fmt[]="%21.15e";
00035
00036
00037
00038
00039
00040 void output_pdf(MeshS *pM, OutputS *pOut)
00041 {
00042 GridS *pG;
00043 FILE *pfile;
00044 char fmt[80];
00045 char fid[80];
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;
00052 double adev=0.0, sdev=0.0;
00053 double skew=0.0, kurt=0.0;
00054 double r, s, ep=0.0;
00055 #ifdef MPI_PARALLEL
00056 DomainS *pD;
00057 int ierr, cd_data_cnt, myID_Comm_Domain;
00058 #endif
00059
00060
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
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
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
00089
00090 #ifdef MPI_PARALLEL
00091 size_pdf = (int)sqrt((double)cd_data_cnt);
00092 #else
00093 size_pdf = (int)sqrt((double)size_dat);
00094 #endif
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){
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
00106 }
00107
00108
00109
00110 dmin = dmax = (*pOut->expr)(pG,is,js,ks);
00111
00112
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;
00138
00139 #else
00140
00141 mean /= (double)data_cnt;
00142
00143 #endif
00144
00145 if(data_cnt > 1){
00146
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
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
00185 }
00186
00187
00188 pOut->gmin = dmin < pOut->gmin ? dmin : pOut->gmin;
00189 pOut->gmax = dmax > pOut->gmax ? dmax : pOut->gmax;
00190
00191
00192
00193
00194
00195 if(dmax - dmin > 0.0){
00196
00197 for(n=0; n<size_pdf; n++) pdf[n] = 0;
00198
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
00210 ierr=MPI_Reduce(pdf,cd_pdf,size_pdf,MPI_INT,MPI_SUM,0,pD->Comm_Domain);
00211
00212 #endif
00213
00214 #ifdef MPI_PARALLEL
00215
00216 if(myID_Comm_Domain != 0) return;
00217 #endif
00218
00219
00220
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
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
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
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
00278 for(n=0; n<size_pdf; n++){
00279
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
00286 fprintf(pfile, fmt, dat, dpdf);
00287 }
00288 }
00289 else
00290 fprintf(pfile,fmt,dmax,1.0);
00291
00292 fclose(pfile);
00293
00294
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
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
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 }