Athena: dump_history.c Source File

Go to the documentation of this file.
00001 #include "copyright.h"
00002 /*============================================================================*/
00003 /*! \file dump_history.c
00004  *  \brief Functions to write dumps of scalar "history" variables in a
00005  *  formatted table.
00006  *
00007  * PURPOSE: Functions to write dumps of scalar "history" variables in a
00008  *   formatted table.  "History" variables are mostly VOLUME AVERAGES, e.g.
00009  *   - S = \Sum_{ijk} q[k][j][i]*dx1[i]*dx2[j]*dx3[k] /
00010  *             \Sum_{ijk} dx1[i]*dx2[j]*dx3[k],
00011  *
00012  *   although other quantities like time and dt are also output.  To compute
00013  *   TOTAL, VOLUME INTEGRATED quantities, just multiply by the Domain volume,
00014  *   which is output in the first line of the history file.  History data is
00015  *   written periodically, giving the time-evolution of that quantitity.
00016  *   Default (hardwired) values are:
00017  *   - scal[0] = time
00018  *   - scal[1] = dt
00019  *   - scal[2] = mass
00020  *   - scal[3] = total energy
00021  *   - scal[4] = d*v1
00022  *   - scal[5] = d*v2
00023  *   - scal[6] = d*v3
00024  *   - scal[7] = 0.5*d*v1**2
00025  *   - scal[8] = 0.5*d*v2**2
00026  *   - scal[9] = 0.5*d*v3**2
00027  *   - scal[10] = 0.5*b1**2
00028  *   - scal[11] = 0.5*b2**2
00029  *   - scal[12] = 0.5*b3**2
00030  *   - scal[13] = d*Phi
00031  *   - scal[14+NSCALARS] = passively advected scalars
00032  *   - scal[15+NSCALARS] = angular momentum
00033  *
00034  * More variables can be hardwired by increasing NSCAL=number of variables, and
00035  * adding calculation of desired quantities below.
00036  *
00037  * Alternatively, up to MAX_USR_H_COUNT new history variables can be added using
00038  * dump_history_enroll() in the problem generator.
00039  *
00040  * With SMR, data is averaged over each Domain separately, and dumped to
00041  * separate files with the level and domain number encoded in the filename.
00042  * Dumps are always made for all levels and domains, and are written in lev#
00043  * directories of the root (rank=0) process FOR THAT DOMAIN COMMUNICATOR.
00044  *
00045  * With SR, the default (hardwired) variables are:
00046  *   - scal[0] = time
00047  *   - scal[1] = dt
00048  *   - scal[2] = mass
00049  *   - scal[3] = total energy
00050  *   - scal[4] = x1-mom
00051  *   - scal[5] = x2-mom
00052  *   - scal[6] = x3-mom
00053  *   - scal[7] = (U^t)**2
00054  *   - scal[8] = (U^x)**2
00055  *   - scal[9] = (U^y)**2
00056  *   - scal[10] = (U^z)**2
00057  *   - scal[11] = (b^t)**2
00058  *   - scal[12] = (b^x)**2
00059  *   - scal[13] = (b^y)**2
00060  *   - scal[14] = (b^z)**2
00061  *   - scal[15] = |b|**2
00062  *   - scal[16] = T^00_EM
00063  *
00064  * CONTAINS PUBLIC FUNCTIONS: 
00065  * - dump_history()        - Writes variables as formatted table
00066  * - dump_history_enroll() - Adds new user-defined history variables          */
00067 /*============================================================================*/
00068 
00069 #include <stdio.h>
00070 #include "defs.h"
00071 #include "athena.h"
00072 #include "globals.h"
00073 #include "prototypes.h"
00074 
00075 /* Maximum Number of default history dump columns. */
00076 #define NSCAL 14
00077 
00078 /* Maximum number of history dump columns that the user routine can add. */
00079 #define MAX_USR_H_COUNT 30
00080 
00081 /* Array of strings / labels for the user added history column. */
00082 static char *usr_label[MAX_USR_H_COUNT];
00083 
00084 /* Array of history dump function pointers for user added history columns. */
00085 static ConsFun_t phst_fun[MAX_USR_H_COUNT];
00086 
00087 static int usr_hst_cnt = 0; /* User History Counter <= MAX_USR_H_COUNT */
00088 
00089 /*----------------------------------------------------------------------------*/
00090 /*! \fn void dump_history(MeshS *pM, OutputS *pOut)
00091  *  \brief Function to write dumps of scalar "history" variables in a
00092  *  formatted table. */
00093 
00094 void dump_history(MeshS *pM, OutputS *pOut)
00095 {
00096   GridS *pG;
00097   DomainS *pD;
00098   int i,j,k,is,ie,js,je,ks,ke,nl,nd;
00099   double dVol, scal[NSCAL + NSCALARS + MAX_USR_H_COUNT], d1;
00100   FILE *pfile;
00101   char *fname,*plev=NULL,*pdom=NULL,*pdir=NULL,fmt[80];
00102   char levstr[8],domstr[8],dirstr[20];
00103   int n, total_hst_cnt, mhst, myID_Comm_Domain=1;
00104 #ifdef MPI_PARALLEL
00105   double my_scal[NSCAL + NSCALARS + MAX_USR_H_COUNT]; /* My Volume averages */
00106   int ierr;
00107 #endif
00108 #ifdef CYLINDRICAL
00109   Real x1,x2,x3;
00110 #endif
00111 #ifdef SPECIAL_RELATIVITY
00112   PrimS W;
00113   Real g, g2, g_2;
00114   Real bx, by, bz, vB, b2, Bmag2;
00115 #endif
00116 
00117 
00118   total_hst_cnt = 9 + NSCALARS + usr_hst_cnt;
00119 #ifdef ADIABATIC
00120   total_hst_cnt++;
00121 #endif
00122 #ifdef MHD
00123   total_hst_cnt += 3;
00124 #endif
00125 #ifdef SELF_GRAVITY
00126   total_hst_cnt += 1;
00127 #endif
00128 #ifdef CYLINDRICAL
00129   total_hst_cnt++;  /* for angular momentum */
00130 #endif
00131 #ifdef SPECIAL_RELATIVITY
00132   total_hst_cnt = 12 + usr_hst_cnt;
00133 #ifdef MHD
00134    total_hst_cnt += 6;
00135 #endif
00136 #endif
00137 
00138 
00139 /* Add a white space to the format */
00140   if(pOut->dat_fmt == NULL){
00141     sprintf(fmt," %%14.6e"); /* Use a default format */
00142   }
00143   else{
00144     sprintf(fmt," %s",pOut->dat_fmt);
00145   }
00146 
00147 /* store time and dt in first two elements of output vector */
00148 
00149   scal[0] = pM->time;
00150   scal[1] = pM->dt;
00151 
00152 /* Loop over all Domains in Mesh, and output Grid data */
00153 
00154   for (nl=0; nl<(pM->NLevels); nl++){
00155 
00156     for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
00157 
00158       if (pM->Domain[nl][nd].Grid != NULL){
00159         pG = pM->Domain[nl][nd].Grid;
00160         pD = (DomainS*)&(pM->Domain[nl][nd]);
00161         is = pG->is, ie = pG->ie;
00162         js = pG->js, je = pG->je;
00163         ks = pG->ks, ke = pG->ke;
00164 
00165         for (i=2; i<total_hst_cnt; i++) {
00166           scal[i] = 0.0;
00167         }
00168  
00169 /* Compute history variables */
00170 
00171         for (k=ks; k<=ke; k++) {
00172           for (j=js; j<=je; j++) {
00173             for (i=is; i<=ie; i++) {
00174               dVol = 1.0; 
00175               if (pG->dx1 > 0.0) dVol *= pG->dx1;
00176               if (pG->dx2 > 0.0) dVol *= pG->dx2;
00177               if (pG->dx3 > 0.0) dVol *= pG->dx3;
00178 #ifndef SPECIAL_RELATIVITY
00179 #ifdef CYLINDRICAL
00180               cc_pos(pG,i,j,k,&x1,&x2,&x3);
00181               dVol *= x1;
00182 #endif
00183 
00184               mhst = 2;
00185               scal[mhst] += dVol*pG->U[k][j][i].d;
00186               d1 = 1.0/pG->U[k][j][i].d;
00187 #ifndef BAROTROPIC
00188               mhst++;
00189               scal[mhst] += dVol*pG->U[k][j][i].E;
00190 #endif
00191               mhst++;
00192               scal[mhst] += dVol*pG->U[k][j][i].M1;
00193               mhst++;
00194               scal[mhst] += dVol*pG->U[k][j][i].M2;
00195               mhst++;
00196               scal[mhst] += dVol*pG->U[k][j][i].M3;
00197               mhst++;
00198               scal[mhst] += dVol*0.5*SQR(pG->U[k][j][i].M1)*d1;
00199               mhst++;
00200               scal[mhst] += dVol*0.5*SQR(pG->U[k][j][i].M2)*d1;
00201               mhst++;
00202               scal[mhst] += dVol*0.5*SQR(pG->U[k][j][i].M3)*d1;
00203 #ifdef MHD
00204               mhst++;
00205               scal[mhst] += dVol*0.5*SQR(pG->U[k][j][i].B1c);
00206               mhst++;
00207               scal[mhst] += dVol*0.5*SQR(pG->U[k][j][i].B2c);
00208               mhst++;
00209               scal[mhst] += dVol*0.5*SQR(pG->U[k][j][i].B3c);
00210 #endif
00211 #ifdef SELF_GRAVITY
00212               mhst++;
00213               scal[mhst] += dVol*pG->U[k][j][i].d*pG->Phi[k][j][i];
00214 #endif
00215 #if (NSCALARS > 0)
00216               for(n=0; n<NSCALARS; n++){
00217                 mhst++;
00218                 scal[mhst] += dVol*pG->U[k][j][i].s[n];
00219               }
00220 #endif
00221 
00222 #ifdef CYLINDRICAL
00223               mhst++;
00224               scal[mhst] += dVol*(x1*pG->U[k][j][i].M2);
00225 #endif
00226 
00227 #else /* SPECIAL_RELATIVITY */
00228 
00229               W = Cons_to_Prim (&(pG->U[k][j][i]));
00230         
00231               /* calculate gamma */
00232               g   = pG->U[k][j][i].d/W.d;
00233               g2  = SQR(g);
00234               g_2 = 1.0/g2;
00235 
00236               mhst = 2;
00237               scal[mhst] += dVol*pG->U[k][j][i].d;
00238               mhst++;
00239               scal[mhst] += dVol*pG->U[k][j][i].E;
00240               mhst++;
00241               scal[mhst] += dVol*pG->U[k][j][i].M1;
00242               mhst++;
00243               scal[mhst] += dVol*pG->U[k][j][i].M2;
00244               mhst++;
00245               scal[mhst] += dVol*pG->U[k][j][i].M3;
00246 
00247               mhst++;
00248               scal[mhst] += dVol*SQR(g);
00249               mhst++;
00250               scal[mhst] += dVol*SQR(g*W.V1);
00251               mhst++;
00252               scal[mhst] += dVol*SQR(g*W.V2);
00253               mhst++;
00254               scal[mhst] += dVol*SQR(g*W.V3);
00255 
00256               mhst++;
00257               scal[mhst] += dVol*W.P;
00258 
00259 #ifdef MHD
00260 
00261               vB = W.V1*pG->U[k][j][i].B1c + W.V2*W.B2c + W.V3*W.B3c;
00262               Bmag2 = SQR(pG->U[k][j][i].B1c) + SQR(W.B2c) + SQR(W.B3c);
00263         
00264               bx = g*(pG->U[k][j][i].B1c*g_2 + vB*W.V1);
00265               by = g*(W.B2c*g_2 + vB*W.V2);
00266               bz = g*(W.B3c*g_2 + vB*W.V3);
00267         
00268               b2 = Bmag2*g_2 + vB*vB;
00269 
00270               mhst++;
00271               scal[mhst] += dVol*(g*vB*g*vB);
00272               mhst++;
00273               scal[mhst] += dVol*bx*bx;
00274               mhst++;
00275               scal[mhst] += dVol*by*by;
00276               mhst++;
00277               scal[mhst] += dVol*bz*bz;
00278               mhst++;
00279               scal[mhst] += dVol*b2;
00280               mhst++;
00281               scal[mhst] += dVol*(Bmag2*(1.0 - 0.5*g_2) - SQR(vB) / 2.0);
00282 
00283 #endif /* MHD */
00284 
00285 #endif  /* SPECIAL_RELATIVITY */
00286 
00287 /* Calculate the user defined history variables */
00288               for(n=0; n<usr_hst_cnt; n++){
00289                 mhst++;
00290                 scal[mhst] += dVol*(*phst_fun[n])(pG, i, j, k);
00291               }
00292             }
00293           }
00294         }
00295 
00296 /* Compute the sum over all Grids in Domain */
00297 
00298 #ifdef MPI_PARALLEL 
00299         for(i=2; i<total_hst_cnt; i++){
00300           my_scal[i] = scal[i];
00301         }
00302         ierr = MPI_Reduce(&(my_scal[2]), &(scal[2]), (total_hst_cnt - 2),
00303           MPI_DOUBLE, MPI_SUM, 0, pD->Comm_Domain);
00304 #endif
00305 
00306 /* Only the parent (rank=0) process computes the average and writes output.
00307  * For single-processor jobs, myID_Comm_world is always zero. */
00308 
00309 #ifdef MPI_PARALLEL
00310         ierr = MPI_Comm_rank(pD->Comm_Domain, &myID_Comm_Domain);
00311 #endif
00312         if((myID_Comm_Domain==0) || (myID_Comm_world==0)){  /* I'm the parent */
00313 
00314 /* Compute volume averages */
00315 
00316           dVol = pD->MaxX[0] - pD->MinX[0];
00317 #ifdef CYLINDRICAL
00318           dVol = (pD->MaxX[0]*pD->MaxX[0]) - (pD->MinX[0]*pD->MinX[0]);
00319 #endif
00320           if (pD->Nx[1] > 1) dVol *= (pD->MaxX[1] - pD->MinX[1]);
00321           if (pD->Nx[2] > 1) dVol *= (pD->MaxX[2] - pD->MinX[2]);
00322           for(i=2; i<total_hst_cnt; i++){
00323             scal[i] /= dVol;
00324           }
00325 
00326 /* Create filename and open file.  History files are always written in lev#
00327  * directories of root process (rank=0 in MPI_COMM_WORLD) */
00328 #ifdef MPI_PARALLEL
00329           if (nl>0) {
00330             plev = &levstr[0];
00331             sprintf(plev,"lev%d",nl);
00332             pdir = &dirstr[0];
00333             sprintf(pdir,"../id0/lev%d",nl);
00334           }
00335 #else
00336           if (nl>0) {
00337             plev = &levstr[0];
00338             sprintf(plev,"lev%d",nl);
00339             pdir = &dirstr[0];
00340             sprintf(pdir,"lev%d",nl);
00341           }
00342 #endif
00343 
00344           if (nd>0) {
00345             pdom = &domstr[0];
00346             sprintf(pdom,"dom%d",nd);
00347           }
00348 
00349           fname = ath_fname(pdir,pM->outfilename,plev,pdom,0,0,NULL,"hst");
00350           if(fname == NULL){
00351             ath_perr(-1,"[dump_history]: Unable to create history filename\n");
00352           }
00353           pfile = fopen(fname,"a");
00354           if(pfile == NULL){
00355             ath_perr(-1,"[dump_history]: Unable to open the history file\n");
00356           }
00357 
00358 /* Write out column headers, but only for first dump */
00359 
00360           mhst = 0;
00361           if(pOut->num == 0){
00362             fprintf(pfile,
00363          "# Athena history dump for level=%i domain=%i volume=%e\n",nl,nd,dVol);
00364             mhst++;
00365             fprintf(pfile,"#   [%i]=time   ",mhst);
00366             mhst++;
00367             fprintf(pfile,"   [%i]=dt      ",mhst);
00368 #ifndef SPECIAL_RELATIVITY
00369             mhst++;
00370             fprintf(pfile,"   [%i]=mass    ",mhst);
00371 #ifdef ADIABATIC
00372             mhst++;
00373             fprintf(pfile,"   [%i]=total E ",mhst);
00374 #endif
00375             mhst++;
00376             fprintf(pfile,"   [%i]=x1 Mom. ",mhst);
00377             mhst++;
00378             fprintf(pfile,"   [%i]=x2 Mom. ",mhst);
00379             mhst++;
00380             fprintf(pfile,"   [%i]=x3 Mom. ",mhst);
00381             mhst++;
00382             fprintf(pfile,"   [%i]=x1-KE   ",mhst);
00383             mhst++;
00384             fprintf(pfile,"   [%i]=x2-KE   ",mhst);
00385             mhst++;
00386             fprintf(pfile,"   [%i]=x3-KE   ",mhst);
00387 #ifdef MHD
00388             mhst++;
00389             fprintf(pfile,"   [%i]=x1-ME   ",mhst);
00390             mhst++;
00391             fprintf(pfile,"   [%i]=x2-ME   ",mhst);
00392             mhst++;
00393             fprintf(pfile,"   [%i]=x3-ME   ",mhst);
00394 #endif
00395 #ifdef SELF_GRAVITY
00396             mhst++;
00397             fprintf(pfile,"   [%i]=grav PE ",mhst);
00398 #endif
00399 #if (NSCALARS > 0)
00400             for(n=0; n<NSCALARS; n++){
00401               mhst++;
00402               fprintf(pfile,"  [%i]=scalar %i",mhst,n);
00403             }
00404 #endif
00405 
00406 #ifdef CYLINDRICAL
00407             mhst++;
00408             fprintf(pfile,"   [%i]=Ang.Mom.",mhst);
00409 #endif
00410 
00411 #else /* SPECIAL_RELATIVITY */
00412             mhst++;
00413             fprintf(pfile,"   [%i]=mass    ",mhst);
00414             mhst++;
00415             fprintf(pfile,"   [%i]=total E ",mhst);
00416             mhst++;
00417             fprintf(pfile,"   [%i]=x1 Mom. ",mhst);
00418             mhst++;
00419             fprintf(pfile,"   [%i]=x2 Mom. ",mhst);
00420             mhst++;
00421             fprintf(pfile,"   [%i]=x3 Mom." ,mhst);
00422             mhst++;
00423             fprintf(pfile,"   [%i]=Gamma   ",mhst);
00424             mhst++;
00425             fprintf(pfile,"   [%i]=x1-KE   ",mhst);
00426             mhst++;
00427             fprintf(pfile,"   [%i]=x2-KE   ",mhst);
00428             mhst++;
00429             fprintf(pfile,"   [%i]=x3-KE  " ,mhst);
00430             mhst++;
00431             fprintf(pfile,"   [%i]=Press  " ,mhst);
00432 #ifdef MHD
00433             mhst++;
00434             fprintf(pfile,"   [%i]=x0-ME  " ,mhst);
00435             mhst++;
00436             fprintf(pfile,"   [%i]=x1-ME  " ,mhst);
00437             mhst++;
00438             fprintf(pfile,"   [%i]=x2-ME  " ,mhst);
00439             mhst++;
00440             fprintf(pfile,"   [%i]=x3-ME  " ,mhst);
00441             mhst++;
00442             fprintf(pfile,"   [%i]=bsq    " ,mhst);
00443             mhst++;
00444             fprintf(pfile,"   [%i]=T^00_EM" ,mhst);
00445 #endif
00446 #endif /* SPECIAL_RELATIVITY */
00447 
00448             for(n=0; n<usr_hst_cnt; n++){
00449               mhst++;
00450               fprintf(pfile,"  [%i]=%s",mhst,usr_label[n]);
00451             }
00452             fprintf(pfile,"\n#\n");
00453           }
00454 
00455 /* Write out data, and close file */
00456 
00457           for (i=0; i<total_hst_cnt; i++) {
00458             fprintf(pfile,fmt,scal[i]);
00459           }
00460           fprintf(pfile,"\n");
00461           fclose(pfile);
00462   
00463         }
00464       }
00465     }
00466   }
00467 
00468   return;
00469 }
00470 
00471 /*----------------------------------------------------------------------------*/
00472 /*! \fn void dump_history_enroll(const ConsFun_t pfun, const char *label)
00473  *  \brief Adds new user-defined history variables            */
00474 
00475 void dump_history_enroll(const ConsFun_t pfun, const char *label){
00476 
00477   if(usr_hst_cnt >= MAX_USR_H_COUNT)
00478     ath_error("[dump_history_enroll]: MAX_USR_H_COUNT = %d exceeded\n",
00479               MAX_USR_H_COUNT);
00480 
00481 /* Copy the label string */
00482   if((usr_label[usr_hst_cnt] = ath_strdup(label)) == NULL)
00483     ath_error("[dump_history_enroll]: Error on sim_strdup(\"%s\")\n",label);
00484 
00485 /* Store the function pointer */
00486   phst_fun[usr_hst_cnt] = pfun;
00487 
00488   usr_hst_cnt++;
00489 
00490   return;
00491 
00492 }
00493 
00494 #undef NSCAL
00495 #undef MAX_USR_H_COUNT