Athena: dump_tab.c Source File

Go to the documentation of this file.
00001 #include "copyright.h"
00002 /*============================================================================*/
00003 /*! \file dump_tab.c
00004  *  \brief Functions to write a dump as a formatted table.
00005  *
00006  * PURPOSE: Functions to write a dump as a formatted table.  The resulting
00007  *   output files can be extremely large, so they are realy only useful for 1D
00008  *   calculations, some 2D calculations, and for very small 3D runs.  With SMR,
00009  *   dumps are made for all levels and domains, unless nlevel and ndomain are
00010  *   specified in <output> block.
00011  *
00012  * REMINDER: use the slicing option available in output_tab() to write selected
00013  *   variables as a formatted table along any arbitrary 1D slice, or in any
00014  *   sub-volume, of 2D or 3D calculations.
00015  *
00016  * CONTAINS PUBLIC FUNCTIONS:
00017  * - dump_tab_cons() - writes conserved variables as formatted table
00018  * - dump_tab_prim() - writes primitive variables as formatted table          */
00019 /*============================================================================*/
00020 
00021 #include <stdio.h>
00022 #include <string.h>
00023 #include <math.h>
00024 #include "defs.h"
00025 #include "athena.h"
00026 #include "globals.h"
00027 #include "prototypes.h"
00028 #ifdef PARTICLES
00029 #include "particles/particle.h"
00030 #endif
00031 
00032 /*----------------------------------------------------------------------------*/
00033 /*! \fn void dump_tab_cons(MeshS *pM, OutputS *pOut)
00034  *  \brief Output CONSERVED variables  */
00035 
00036 void dump_tab_cons(MeshS *pM, OutputS *pOut)
00037 {
00038   GridS *pG;
00039   int nl,nd,i,j,k,il,iu,jl,ju,kl,ku;
00040   FILE *pfile;
00041   char *fname,*plev=NULL,*pdom=NULL;
00042   char levstr[8],domstr[8];
00043   Real x1,x2,x3;
00044   char zone_fmt[20], fmt[80];
00045   int col_cnt, nmax;
00046 #if (NSCALARS > 0)
00047   int n;
00048 #endif
00049 
00050 /* Add a white space to the format, setup format for integer zone columns */
00051   if(pOut->dat_fmt == NULL){
00052      sprintf(fmt," %%12.8e"); /* Use a default format */
00053   }
00054   else{
00055     sprintf(fmt," %s",pOut->dat_fmt);
00056   }
00057 
00058 /* Loop over all Domains in Mesh, and output Grid data */
00059 
00060   for (nl=0; nl<(pM->NLevels); nl++){
00061     for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
00062       if (pM->Domain[nl][nd].Grid != NULL){
00063 
00064 /* write files if domain and level match input, or are not specified (-1) */
00065       if ((pOut->nlevel == -1 || pOut->nlevel == nl) &&
00066           (pOut->ndomain == -1 || pOut->ndomain == nd)){
00067         pG = pM->Domain[nl][nd].Grid;
00068         col_cnt = 1;
00069 
00070 /* construct output filename. */
00071         if (nl>0) {
00072           plev = &levstr[0];
00073           sprintf(plev,"lev%d",nl);
00074         }
00075         if (nd>0) {
00076           pdom = &domstr[0];
00077           sprintf(pdom,"dom%d",nd);
00078         }
00079 
00080         if((fname = ath_fname(plev,pM->outfilename,plev,pdom,num_digit,
00081             pOut->num,NULL,"tab")) == NULL){
00082           ath_error("[dump_tab]: Error constructing filename\n");
00083         }
00084 
00085 /* open output file */
00086         if((pfile = fopen(fname,"w")) == NULL){
00087           ath_error("[dump_tab]: Unable to open ppm file %s\n",fname);
00088         }
00089 
00090 /* Upper and Lower bounds on i,j,k for data dump */
00091         il = pG->is; iu = pG->ie;
00092         jl = pG->js; ju = pG->je;
00093         kl = pG->ks; ku = pG->ke;
00094 
00095         nmax =  pG->Nx[0] > pG->Nx[1]  ? pG->Nx[0] : pG->Nx[1];
00096         nmax = (pG->Nx[2] > nmax ? pG->Nx[2] : nmax);
00097 
00098 #ifdef WRITE_GHOST_CELLS
00099         iu = pG->ie + nghost;
00100         il = pG->is - nghost;
00101 
00102         if(pG->Nx[1] > 1) {
00103           ju = pG->je + nghost;
00104           jl = pG->js - nghost;
00105         }
00106 
00107         if(pG->Nx[2] > 1) {
00108           ku = pG->ke + nghost;
00109           kl = pG->ks - nghost;
00110         }
00111         nmax += 2*nghost;
00112 #endif
00113         sprintf(zone_fmt,"%%%dd", (int)(2+log10((double)(nmax))));
00114 
00115 /* Write out some header information */
00116 
00117         if (pG->Nx[0] > 1) {
00118           fprintf(pfile,"# Nx1 = %d\n",iu-il+1);
00119           fprintf(pfile,"# x1-size = %g\n",(iu-il+1)*pG->dx1);
00120         }
00121         if (pG->Nx[1] > 1) {
00122           fprintf(pfile,"# Nx2 = %d\n",ju-jl+1);
00123           fprintf(pfile,"# x2-size = %g\n",(ju-jl+1)*pG->dx2);
00124         }
00125         if (pG->Nx[2] > 1) {
00126           fprintf(pfile,"# Nx3 = %d\n",ku-kl+1);
00127           fprintf(pfile,"# x3-size = %g\n",(ku-kl+1)*pG->dx3);
00128         }
00129         fprintf(pfile,"# CONSERVED vars at Time= %g, level= %i, domain= %i\n",
00130           pM->time,nl,nd);
00131 
00132 /* write out i,j,k column headers.  Note column number is embedded in header */
00133 
00134         fprintf(pfile,"# [%d]=i-zone",col_cnt);
00135         col_cnt++;
00136         if (pG->Nx[1] > 2) {
00137           fprintf(pfile," [%d]=j-zone",col_cnt);
00138           col_cnt++;
00139         }
00140         if (pG->Nx[2] > 3) {
00141           fprintf(pfile," [%d]=k-zone",col_cnt);
00142           col_cnt++;
00143         }
00144 
00145 /* write out x1,x2,x3 column headers.  */
00146 
00147         if (pG->Nx[0] > 1) {
00148           fprintf(pfile," [%d]=x1",col_cnt);
00149           col_cnt++;
00150         }
00151         if (pG->Nx[1] > 2) {
00152           fprintf(pfile," [%d]=x2",col_cnt);
00153           col_cnt++;
00154         }
00155         if (pG->Nx[2] > 3) {
00156           fprintf(pfile," [%d]=x3",col_cnt);
00157           col_cnt++;
00158         }
00159 
00160 /* write out d,M1,M2,M3 column headers */
00161 
00162         fprintf(pfile," [%d]=d",col_cnt);
00163         col_cnt++;
00164         fprintf(pfile," [%d]=M1",col_cnt);
00165         col_cnt++;
00166         fprintf(pfile," [%d]=M2",col_cnt);
00167         col_cnt++;
00168         fprintf(pfile," [%d]=M3",col_cnt);
00169         col_cnt++;
00170 
00171 /* write out E column header, if not barotropic */
00172 #ifndef BAROTROPIC
00173         fprintf(pfile," [%d]=E",col_cnt);
00174         col_cnt++;
00175 #endif /* BAROTROPIC */
00176 
00177 /* write out magnetic field component column headers, if mhd */
00178 #ifdef MHD
00179         fprintf(pfile," [%d]=B1c",col_cnt);
00180         col_cnt++;
00181         fprintf(pfile," [%d]=B2c",col_cnt);
00182         col_cnt++;
00183         fprintf(pfile," [%d]=B3c",col_cnt);
00184         col_cnt++;
00185 #endif /* MHD */
00186 
00187 /* write out column header for gravitational potential (self-gravity) */
00188 #ifdef SELF_GRAVITY
00189         fprintf(pfile," [%d]=Phi",col_cnt);
00190         col_cnt++;
00191 #endif
00192 
00193 /* write out column headers for particles */
00194 #ifdef PARTICLES
00195         if (pOut->out_pargrid) {
00196           fprintf(pfile," [%d]=dpar",col_cnt);
00197           col_cnt++;
00198           fprintf(pfile," [%d]=M1par",col_cnt);
00199           col_cnt++;
00200           fprintf(pfile," [%d]=M2par",col_cnt);
00201           col_cnt++;
00202           fprintf(pfile," [%d]=M3par",col_cnt);
00203           col_cnt++;
00204         }
00205 #endif
00206 
00207 /* write out column headers for passive scalars */
00208 #if (NSCALARS > 0)
00209         for (n=0; n<NSCALARS; n++) {
00210           fprintf(pfile," [%d]=s%d",col_cnt,n);
00211           col_cnt++;
00212         }
00213 #endif
00214         fprintf(pfile,"\n");
00215 
00216 /* Write out data */
00217 
00218         for(k=kl; k<=ku; k++){
00219           for(j=jl; j<=ju; j++){
00220             for(i=il; i<=iu; i++){
00221               cc_pos(pG,i,j,k,&x1,&x2,&x3);
00222 
00223               if (pG->Nx[0] > 1) fprintf(pfile,zone_fmt,i);
00224               if (pG->Nx[1] > 1) fprintf(pfile,zone_fmt,j);
00225               if (pG->Nx[2] > 1) fprintf(pfile,zone_fmt,k);
00226               if (pG->Nx[0] > 1) fprintf(pfile,fmt,x1);
00227               if (pG->Nx[1] > 1) fprintf(pfile,fmt,x2);
00228               if (pG->Nx[2] > 1) fprintf(pfile,fmt,x3);
00229 
00230 /* Dump all variables */
00231 
00232               fprintf(pfile,fmt,pG->U[k][j][i].d);
00233               fprintf(pfile,fmt,pG->U[k][j][i].M1);
00234               fprintf(pfile,fmt,pG->U[k][j][i].M2);
00235               fprintf(pfile,fmt,pG->U[k][j][i].M3);
00236 
00237 #ifndef BAROTROPIC
00238               fprintf(pfile,fmt,pG->U[k][j][i].E);
00239 #endif /* BAROTROPIC */
00240 
00241 #ifdef MHD
00242               fprintf(pfile,fmt,pG->U[k][j][i].B1c);
00243               fprintf(pfile,fmt,pG->U[k][j][i].B2c);
00244               fprintf(pfile,fmt,pG->U[k][j][i].B3c);
00245 #endif
00246 
00247 #ifdef SELF_GRAVITY
00248               fprintf(pfile,fmt,pG->Phi[k][j][i]);
00249 #endif
00250 
00251 #ifdef PARTICLES
00252               if (pOut->out_pargrid) {
00253                 fprintf(pfile,fmt,pG->Coup[k][j][i].grid_d);
00254                 fprintf(pfile,fmt,pG->Coup[k][j][i].grid_v1);
00255                 fprintf(pfile,fmt,pG->Coup[k][j][i].grid_v2);
00256                 fprintf(pfile,fmt,pG->Coup[k][j][i].grid_v3);
00257               }
00258 #endif
00259 
00260 #if (NSCALARS > 0)
00261               for (n=0; n<NSCALARS; n++) fprintf(pfile,fmt,pG->U[k][j][i].s[n]);
00262 #endif
00263 
00264               fprintf(pfile,"\n");
00265             }
00266           }
00267         }
00268       }}
00269     } /* end loop over domains */
00270   } /* end loop over levels */
00271 
00272   fclose(pfile);
00273 
00274   return;
00275 }
00276 
00277 /*----------------------------------------------------------------------------*/
00278 /*! \fn void dump_tab_prim(MeshS *pM, OutputS *pOut)
00279  *  \brief Output PRIMITIVE variables.  */
00280 
00281 void dump_tab_prim(MeshS *pM, OutputS *pOut)
00282 {
00283   GridS *pG;
00284   int nl,nd,i,j,k,il,iu,jl,ju,kl,ku;
00285   FILE *pfile;
00286   char *fname,*plev=NULL,*pdom=NULL;
00287   char levstr[8],domstr[8];
00288   PrimS W;
00289   Real x1,x2,x3,d1;
00290   char zone_fmt[20], fmt[80];
00291   int col_cnt, nmax;
00292 #if (NSCALARS > 0)
00293   int n;
00294 #endif
00295 
00296 /* Add a white space to the format, setup format for integer zone columns */
00297   if(pOut->dat_fmt == NULL){
00298     sprintf(fmt," %%12.8e"); /* Use a default format */
00299   }
00300   else{
00301     sprintf(fmt," %s",pOut->dat_fmt);
00302   }
00303 
00304 /* Loop over all Domains in Mesh, and output Grid data */
00305 
00306   for (nl=0; nl<(pM->NLevels); nl++){
00307     for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
00308       if (pM->Domain[nl][nd].Grid != NULL){
00309 
00310 /* write files if domain and level match input, or are not specified (-1) */
00311       if ((pOut->nlevel == -1 || pOut->nlevel == nl) &&
00312           (pOut->ndomain == -1 || pOut->ndomain == nd)){
00313         pG = pM->Domain[nl][nd].Grid;
00314         col_cnt = 1;
00315 
00316 /* construct output filename. */
00317         if (nl>0) {
00318           plev = &levstr[0];
00319           sprintf(plev,"lev%d",nl);
00320         }
00321         if (nd>0) {
00322           pdom = &domstr[0];
00323           sprintf(pdom,"dom%d",nd);
00324         }
00325 
00326         if((fname = ath_fname(plev,pM->outfilename,plev,pdom,num_digit,
00327             pOut->num,NULL,"tab")) == NULL){
00328           ath_error("[dump_tab]: Error constructing filename\n");
00329         }
00330 
00331 /* open output file */
00332         if((pfile = fopen(fname,"w")) == NULL){
00333           ath_error("[dump_tab]: Unable to open ppm file %s\n",fname);
00334         }
00335 
00336 /* Upper and Lower bounds on i,j,k for data dump */
00337 
00338         il = pG->is; iu = pG->ie;
00339         jl = pG->js; ju = pG->je;
00340         kl = pG->ks; ku = pG->ke;
00341 
00342         nmax =  pG->Nx[0] > pG->Nx[1]  ? pG->Nx[0] : pG->Nx[1];
00343         nmax = (pG->Nx[2] > nmax ? pG->Nx[2] : nmax);
00344 
00345 #ifdef WRITE_GHOST_CELLS
00346         iu = pG->ie + nghost;
00347         il = pG->is - nghost;
00348 
00349         if(pG->Nx[1] > 1) {
00350           ju = pG->je + nghost;
00351           jl = pG->js - nghost;
00352         }
00353 
00354         if(pG->Nx[2] > 1) {
00355           ku = pG->ke + nghost;
00356           kl = pG->ks - nghost;
00357         }
00358         nmax += 2*nghost;
00359 #endif
00360         sprintf(zone_fmt,"%%%dd", (int)(2+log10((double)(nmax))));
00361 
00362 /* Write out some header information */
00363 
00364         if (pG->Nx[0] > 1) {
00365           fprintf(pfile,"# Nx1 = %d\n",iu-il+1);
00366           fprintf(pfile,"# x1-size = %g\n",(iu-il+1)*pG->dx1);
00367         }
00368         if (pG->Nx[1] > 1) {
00369           fprintf(pfile,"# Nx2 = %d\n",ju-jl+1);
00370           fprintf(pfile,"# x2-size = %g\n",(ju-jl+1)*pG->dx2);
00371         }
00372         if (pG->Nx[2] > 1) {
00373           fprintf(pfile,"# Nx3 = %d\n",ku-kl+1);
00374           fprintf(pfile,"# x3-size = %g\n",(ku-kl+1)*pG->dx3);
00375         }
00376         fprintf(pfile,"# PRIMITIVE vars at Time = %g, level= %i, domain= %i\n",
00377           pM->time,nl,nd);
00378 
00379 /* write out i,j,k column headers.  Note column number is embedded in header */
00380 
00381         fprintf(pfile,"# [%d]=i-zone",col_cnt);
00382         col_cnt++;
00383         if (pG->Nx[1] > 2) {
00384           fprintf(pfile," [%d]=j-zone",col_cnt);
00385           col_cnt++;
00386         }
00387         if (pG->Nx[2] > 3) {
00388           fprintf(pfile," [%d]=k-zone",col_cnt);
00389           col_cnt++;
00390         }
00391 
00392 /* write out x1,x2,x3 column headers.  */
00393 
00394         fprintf(pfile," [%d]=x1",col_cnt);
00395         col_cnt++;
00396         if (pG->Nx[1] > 2) {
00397           fprintf(pfile," [%d]=x2",col_cnt);
00398           col_cnt++;
00399         }
00400         if (pG->Nx[2] > 3) {
00401           fprintf(pfile," [%d]=x3",col_cnt);
00402           col_cnt++;
00403         }
00404 
00405 /* write out d,V1,V2,V3 column headers */
00406 
00407         fprintf(pfile," [%d]=d",col_cnt);
00408         col_cnt++;
00409         fprintf(pfile," [%d]=V1",col_cnt);
00410         col_cnt++;
00411         fprintf(pfile," [%d]=V2",col_cnt);
00412         col_cnt++;
00413         fprintf(pfile," [%d]=V3",col_cnt);
00414         col_cnt++;
00415 
00416 /* write out P column header, if not barotropic */
00417 #ifndef BAROTROPIC
00418         fprintf(pfile," [%d]=P",col_cnt);
00419         col_cnt++;
00420 #endif /* BAROTROPIC */
00421 
00422 /* write out magnetic field component column headers, if mhd */
00423 #ifdef MHD
00424         fprintf(pfile," [%d]=B1c",col_cnt);
00425         col_cnt++;
00426         fprintf(pfile," [%d]=B2c",col_cnt);
00427         col_cnt++;
00428         fprintf(pfile," [%d]=B3c",col_cnt);
00429         col_cnt++;
00430 #endif /* MHD */
00431 
00432 /* write out column header for gravitational potential (self-gravity) */
00433 #ifdef SELF_GRAVITY
00434         fprintf(pfile," [%d]=Phi",col_cnt);
00435         col_cnt++;
00436 #endif
00437 
00438 /* write out column headers for particles */
00439 #ifdef PARTICLES
00440         if (pOut->out_pargrid) {
00441           fprintf(pfile," [%d]=dpar",col_cnt);
00442           col_cnt++;
00443           fprintf(pfile," [%d]=V1par",col_cnt);
00444           col_cnt++;
00445           fprintf(pfile," [%d]=V2par",col_cnt);
00446           col_cnt++;
00447           fprintf(pfile," [%d]=V3par",col_cnt);
00448           col_cnt++;
00449         }
00450 #endif
00451 
00452 /* write out column headers for passive scalars */
00453 #if (NSCALARS > 0)
00454         for (n=0; n<NSCALARS; n++) {
00455           fprintf(pfile," [%d]=s%d",col_cnt,n);
00456           col_cnt++;
00457         }
00458 #endif
00459         fprintf(pfile,"\n");
00460 
00461 /* Write out data */
00462 
00463         for(k=kl; k<=ku; k++){
00464           for(j=jl; j<=ju; j++){
00465             for(i=il; i<=iu; i++){
00466               cc_pos(pG,i,j,k,&x1,&x2,&x3);
00467               W = Cons_to_Prim(&(pG->U[k][j][i])); 
00468 
00469               if (pG->Nx[0] > 1) fprintf(pfile,zone_fmt,i);
00470               if (pG->Nx[1] > 1) fprintf(pfile,zone_fmt,j);
00471               if (pG->Nx[2] > 1) fprintf(pfile,zone_fmt,k);
00472               if (pG->Nx[0] > 1) fprintf(pfile,fmt,x1);
00473               if (pG->Nx[1] > 1) fprintf(pfile,fmt,x2);
00474               if (pG->Nx[2] > 1) fprintf(pfile,fmt,x3);
00475 
00476 /* Dump all variables */
00477 
00478               fprintf(pfile,fmt,W.d);
00479               fprintf(pfile,fmt,W.V1);
00480               fprintf(pfile,fmt,W.V2);
00481               fprintf(pfile,fmt,W.V3);
00482 
00483 #ifndef BAROTROPIC
00484               fprintf(pfile,fmt,W.P);
00485 #endif /* BAROTROPIC */
00486 
00487 #ifdef MHD
00488               fprintf(pfile,fmt,W.B1c);
00489               fprintf(pfile,fmt,W.B2c);
00490               fprintf(pfile,fmt,W.B3c);
00491 #endif
00492 
00493 #ifdef SELF_GRAVITY
00494               fprintf(pfile,fmt,pG->Phi[k][j][i]);
00495 #endif
00496 
00497 #ifdef PARTICLES
00498               if (pOut->out_pargrid) {
00499                 fprintf(pfile,fmt,pG->Coup[k][j][i].grid_d);
00500                 if (pG->Coup[k][j][i].grid_d>0.0)
00501                   d1 = 1.0/pG->Coup[k][j][i].grid_d;
00502                 else
00503                   d1 = 0.0;
00504                 fprintf(pfile,fmt,pG->Coup[k][j][i].grid_v1*d1);
00505                 fprintf(pfile,fmt,pG->Coup[k][j][i].grid_v2*d1);
00506                 fprintf(pfile,fmt,pG->Coup[k][j][i].grid_v3*d1);
00507               }
00508 #endif
00509 
00510 #if (NSCALARS > 0)
00511               for (n=0; n<NSCALARS; n++) fprintf(pfile,fmt,W.r[n]);
00512 #endif
00513               fprintf(pfile,"\n");
00514             }
00515           }
00516         }
00517       }}
00518     } /* end loop over domains */
00519   } /* end loop over levels */
00520 
00521   fclose(pfile);
00522 
00523   return;
00524 }