Athena: dump_binary.c Source File

Go to the documentation of this file.
00001 #include "copyright.h"
00002 /*============================================================================*/
00003 /*! \file dump_binary.c
00004  *  \brief Function to write an unformatted dump of the field variables.
00005  *
00006  * PURPOSE: Function to write an unformatted dump of the field variables that
00007  *   can be read, e.g., by IDL scripts.  With SMR, dumps are made for all levels
00008  *   and domains, unless nlevel and ndomain are specified in <output> block.
00009  *
00010  * CONTAINS PUBLIC FUNCTIONS: 
00011  * - dump_binary() - writes either conserved or primitive variables depending
00012  *                 on value of pOut->out read from input block.               */
00013 /*============================================================================*/
00014 
00015 #include <stdio.h>
00016 #include <stdlib.h>
00017 #include <string.h>
00018 #include "defs.h"
00019 #include "athena.h"
00020 #include "globals.h"
00021 #include "prototypes.h"
00022 #ifdef PARTICLES
00023 #include "particles/particle.h"
00024 #endif
00025 
00026 /*----------------------------------------------------------------------------*/
00027 /*! \fn void dump_binary(MeshS *pM, OutputS *pOut)
00028  *  \brief Function to write an unformatted dump of the field variables. */
00029 
00030 void dump_binary(MeshS *pM, OutputS *pOut)
00031 {
00032   GridS *pGrid;
00033   PrimS ***W;
00034   FILE *p_binfile;
00035   char *fname,*plev=NULL,*pdom=NULL;
00036   char levstr[8],domstr[8];
00037   int n,ndata[7];
00038 /* Upper and Lower bounds on i,j,k for data dump */
00039   int i,j,k,il,iu,jl,ju,kl,ku,nl,nd;
00040   float dat[2],*datax,*datay,*dataz;
00041   Real *pData,x1,x2,x3;
00042   int coordsys = -1;
00043 
00044 /* Loop over all Domains in Mesh, and output Grid data */
00045 
00046   for (nl=0; nl<(pM->NLevels); nl++){
00047     for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
00048       if (pM->Domain[nl][nd].Grid != NULL){
00049 
00050 /* write files if domain and level match input, or are not specified (-1) */ 
00051       if ((pOut->nlevel == -1 || pOut->nlevel == nl) &&
00052           (pOut->ndomain == -1 || pOut->ndomain == nd)){
00053         pGrid = pM->Domain[nl][nd].Grid;
00054 
00055         il = pGrid->is, iu = pGrid->ie;
00056         jl = pGrid->js, ju = pGrid->je;
00057         kl = pGrid->ks, ku = pGrid->ke;
00058 
00059 #ifdef WRITE_GHOST_CELLS
00060         il = pGrid->is - nghost;
00061         iu = pGrid->ie + nghost;
00062 
00063         if(pGrid->Nx[1] > 1){
00064           jl = pGrid->js - nghost;
00065           ju = pGrid->je + nghost;
00066         }
00067 
00068         if(pGrid->Nx[2] > 1){
00069           kl = pGrid->ks - nghost;
00070           ku = pGrid->ke + nghost;
00071         }
00072 #endif /* WRITE_GHOST_CELLS */
00073 
00074         ndata[0] = iu-il+1;
00075         ndata[1] = ju-jl+1;
00076         ndata[2] = ku-kl+1;
00077 
00078 /* calculate primitive variables, if needed */
00079 
00080         if(strcmp(pOut->out,"prim") == 0) {
00081           if((W = (PrimS***)
00082             calloc_3d_array(ndata[2],ndata[1],ndata[0],sizeof(PrimS))) == NULL)
00083             ath_error("[dump_bin]: failed to allocate Prim array\n");
00084 
00085           for (k=kl; k<=ku; k++) {
00086           for (j=jl; j<=ju; j++) {
00087           for (i=il; i<=iu; i++) {
00088             W[k-kl][j-jl][i-il] = Cons_to_Prim(&(pGrid->U[k][j][i]));
00089           }}}
00090         }
00091 
00092 /* construct filename, open file */
00093         if (nl>0) {
00094           plev = &levstr[0];
00095           sprintf(plev,"lev%d",nl);
00096         }
00097         if (nd>0) {
00098           pdom = &domstr[0];
00099           sprintf(pdom,"dom%d",nd);
00100         }
00101         if((fname = ath_fname(plev,pM->outfilename,plev,pdom,num_digit,
00102             pOut->num,NULL,"bin")) == NULL){
00103           ath_error("[dump_binary]: Error constructing filename\n");
00104         }
00105 
00106         if((p_binfile = fopen(fname,"wb")) == NULL){
00107           ath_error("[dump_binary]: Unable to open binary dump file\n");
00108           return;
00109         }
00110 
00111 /* Write the coordinate system information */
00112 #if defined CARTESIAN
00113         coordsys = -1;
00114 #elif defined CYLINDRICAL
00115         coordsys = -2;
00116 #elif defined SPHERICAL
00117         coordsys = -3;
00118 #endif
00119         fwrite(&coordsys,sizeof(int),1,p_binfile);
00120 
00121 /* Write number of zones and variables */
00122         ndata[3] = NVAR;
00123         ndata[4] = NSCALARS;
00124 #ifdef SELF_GRAVITY
00125         ndata[5] = 1;
00126 #else
00127         ndata[5] = 0;
00128 #endif
00129 #ifdef PARTICLES
00130         ndata[6] = 1;
00131 #else
00132         ndata[6] = 0;
00133 #endif
00134         fwrite(ndata,sizeof(int),7,p_binfile);
00135 
00136 /* Write (gamma-1) and isothermal sound speed */
00137 
00138 #ifdef ISOTHERMAL
00139         dat[0] = (float)0.0;
00140         dat[1] = (float)Iso_csound;
00141 #elif defined ADIABATIC
00142         dat[0] = (float)Gamma_1 ;
00143         dat[1] = (float)0.0;
00144 #else
00145         dat[0] = dat[1] = 0.0; /* Anything better to put here? */
00146 #endif
00147         fwrite(dat,sizeof(float),2,p_binfile);
00148 
00149 /* Write time, dt */
00150 
00151         dat[0] = (float)pGrid->time;
00152         dat[1] = (float)pGrid->dt;
00153         fwrite(dat,sizeof(float),2,p_binfile);
00154  
00155 /* Allocate Memory */
00156 
00157         if((datax = (float *)malloc(ndata[0]*sizeof(float))) == NULL){
00158           ath_error("[dump_binary]: malloc failed for temporary array\n");
00159           return;
00160         }
00161         if((datay = (float *)malloc(ndata[1]*sizeof(float))) == NULL){
00162           ath_error("[dump_binary]: malloc failed for temporary array\n");
00163           return;
00164         }
00165         if((dataz = (float *)malloc(ndata[2]*sizeof(float))) == NULL){
00166           ath_error("[dump_binary]: malloc failed for temporary array\n");
00167           return;
00168         }
00169 
00170 /* compute x,y,z coordinates of cell centers, and write out */
00171 
00172         for (i=il; i<=iu; i++) {
00173           cc_pos(pGrid,i,jl,kl,&x1,&x2,&x3);
00174           pData = ((Real *) &(x1));
00175           datax[i-il] = (float)(*pData);
00176         }
00177         fwrite(datax,sizeof(float),(size_t)ndata[0],p_binfile);
00178 
00179         for (j=jl; j<=ju; j++) {
00180           cc_pos(pGrid,il,j,kl,&x1,&x2,&x3);
00181           pData = ((Real *) &(x2));
00182           datay[j-jl] = (float)(*pData);
00183         }
00184         fwrite(datay,sizeof(float),(size_t)ndata[1],p_binfile);
00185 
00186         for (k=kl; k<=ku; k++) {
00187           cc_pos(pGrid,il,jl,k,&x1,&x2,&x3);
00188           pData = ((Real *) &(x3));
00189           dataz[k-kl] = (float)(*pData);
00190         }
00191         fwrite(dataz,sizeof(float),(size_t)ndata[2],p_binfile);
00192 
00193 /* Write cell-centered data (either conserved or primitives) */
00194 
00195         for (n=0;n<NVAR; n++) {
00196           for (k=0; k<ndata[2]; k++) {
00197           for (j=0; j<ndata[1]; j++) {
00198             for (i=0; i<ndata[0]; i++) {
00199 
00200               if (strcmp(pOut->out,"cons") == 0){
00201                 pData = ((Real*)&(pGrid->U[k+kl][j+jl][i+il])) + n;
00202               } else if(strcmp(pOut->out,"prim") == 0) {
00203                 pData = ((Real*)&(W[k][j][i])) + n;
00204               }
00205               datax[i] = (float)(*pData);
00206 
00207             }
00208             fwrite(datax,sizeof(float),(size_t)ndata[0],p_binfile);
00209 
00210           }}
00211         }
00212 
00213 #ifdef SELF_GRAVITY
00214         for (k=0; k<ndata[2]; k++) {
00215         for (j=0; j<ndata[1]; j++) {
00216           for (i=0; i<ndata[0]; i++) {
00217             pData = &(pGrid->Phi[k+kl][j+jl][i+il]);
00218             datax[i] = (float)(*pData);
00219           }
00220           fwrite(datax,sizeof(float),(size_t)ndata[0],p_binfile);
00221         }}
00222 #endif
00223 
00224 #ifdef PARTICLES
00225         if (pOut->out_pargrid) {
00226           for (k=0; k<ndata[2]; k++) {
00227           for (j=0; j<ndata[1]; j++) {
00228             for (i=0; i<ndata[0]; i++) {
00229               datax[i] = pGrid->Coup[k+kl][j+jl][i+il].grid_d;
00230             }
00231             fwrite(datax,sizeof(float),(size_t)ndata[0],p_binfile);
00232           }}
00233           for (k=0; k<ndata[2]; k++) {
00234           for (j=0; j<ndata[1]; j++) {
00235             for (i=0; i<ndata[0]; i++) {
00236               datax[i] = pGrid->Coup[k+kl][j+jl][i+il].grid_v1;
00237             }
00238             fwrite(datax,sizeof(float),(size_t)ndata[0],p_binfile);
00239           }}
00240           for (k=0; k<ndata[2]; k++) {
00241           for (j=0; j<ndata[1]; j++) {
00242             for (i=0; i<ndata[0]; i++) {
00243               datax[i] = pGrid->Coup[k+kl][j+jl][i+il].grid_v2;
00244             }
00245             fwrite(datax,sizeof(float),(size_t)ndata[0],p_binfile);
00246           }}
00247           for (k=0; k<ndata[2]; k++) {
00248           for (j=0; j<ndata[1]; j++) {
00249             for (i=0; i<ndata[0]; i++) {
00250               datax[i] = pGrid->Coup[k+kl][j+jl][i+il].grid_v3;
00251             }
00252             fwrite(datax,sizeof(float),(size_t)ndata[0],p_binfile);
00253           }}
00254         }
00255 #endif
00256 
00257 /* close file and free memory */
00258         fclose(p_binfile); 
00259         free(datax); 
00260         free(datay); 
00261         free(dataz); 
00262         if(strcmp(pOut->out,"prim") == 0) free_3d_array(W);
00263       }}
00264     }
00265   }
00266 }