Athena: dump_vtk.c Source File

Go to the documentation of this file.
00001 #include "copyright.h"
00002 /*============================================================================*/
00003 /*! \file dump_vtk.c 
00004  *  \brief Function to write a dump in VTK "legacy" format.
00005  *
00006  * PURPOSE: Function to write a dump in VTK "legacy" format.  With SMR,
00007  *   dumps are made for all levels and domains, unless nlevel and ndomain are
00008  *   specified in <output> block.  Works for BOTH conserved and primitives.
00009  *
00010  * CONTAINS PUBLIC FUNCTIONS: 
00011  * - dump_vtk() - writes VTK dump (all variables).                            */
00012 /*============================================================================*/
00013 
00014 #include <stdio.h>
00015 #include <stdlib.h>
00016 #include <string.h>
00017 #include "defs.h"
00018 #include "athena.h"
00019 #include "prototypes.h"
00020 #ifdef PARTICLES
00021 #include "particles/particle.h"
00022 #endif
00023 
00024 /*----------------------------------------------------------------------------*/
00025 /*! \fn void dump_vtk(MeshS *pM, OutputS *pOut)
00026  *  \brief Writes VTK dump (all variables).                                   */
00027 
00028 void dump_vtk(MeshS *pM, OutputS *pOut)
00029 {
00030   GridS *pGrid;
00031   PrimS ***W;
00032   FILE *pfile;
00033   char *fname,*plev=NULL,*pdom=NULL;
00034   char levstr[8],domstr[8];
00035 /* Upper and Lower bounds on i,j,k for data dump */
00036   int i,j,k,il,iu,jl,ju,kl,ku,nl,nd;
00037   int big_end = ath_big_endian();
00038   int ndata0,ndata1,ndata2;
00039   float *data;   /* points to 3*ndata0 allocated floats */
00040   double x1, x2, x3;
00041 #if (NSCALARS > 0)
00042   int n;
00043 #endif
00044 
00045 /* Loop over all Domains in Mesh, and output Grid data */
00046 
00047   for (nl=0; nl<(pM->NLevels); nl++){
00048     for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
00049       if (pM->Domain[nl][nd].Grid != NULL){
00050 
00051 /* write files if domain and level match input, or are not specified (-1) */
00052       if ((pOut->nlevel == -1 || pOut->nlevel == nl) &&
00053           (pOut->ndomain == -1 || pOut->ndomain == nd)){
00054         pGrid = pM->Domain[nl][nd].Grid;
00055 
00056         il = pGrid->is, iu = pGrid->ie;
00057         jl = pGrid->js, ju = pGrid->je;
00058         kl = pGrid->ks, ku = pGrid->ke;
00059 
00060 #ifdef WRITE_GHOST_CELLS
00061         iu = pGrid->ie + nghost;
00062         il = pGrid->is - nghost;
00063 
00064         if(pGrid->Nx[1] > 1) {
00065           ju = pGrid->je + nghost;
00066           jl = pGrid->js - nghost;
00067         }
00068 
00069         if(pGrid->Nx[2] > 1) {
00070           ku = pGrid->ke + nghost;
00071           kl = pGrid->ks - nghost;
00072         }
00073 #endif /* WRITE_GHOST_CELLS */
00074 
00075         ndata0 = iu-il+1;
00076         ndata1 = ju-jl+1;
00077         ndata2 = ku-kl+1;
00078 
00079 /* calculate primitive variables, if needed */
00080 
00081         if(strcmp(pOut->out,"prim") == 0) {
00082           if((W = (PrimS***)calloc_3d_array(ndata2,ndata1,ndata0,sizeof(PrimS)))
00083              == NULL) ath_error("[dump_vtk]: 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,"vtk")) == NULL){
00103           ath_error("[dump_vtk]: Error constructing filename\n");
00104         }
00105 
00106         if((pfile = fopen(fname,"w")) == NULL){
00107           ath_error("[dump_vtk]: Unable to open vtk dump file\n");
00108           return;
00109         }
00110 
00111 /* Allocate memory for temporary array of floats */
00112 
00113         if((data = (float *)malloc(3*ndata0*sizeof(float))) == NULL){
00114           ath_error("[dump_vtk]: malloc failed for temporary array\n");
00115           return;
00116         }
00117 
00118 /* There are five basic parts to the VTK "legacy" file format.  */
00119 /*  1. Write file version and identifier */
00120 
00121         fprintf(pfile,"# vtk DataFile Version 2.0\n");
00122 
00123 /*  2. Header */
00124 
00125         if (strcmp(pOut->out,"cons") == 0){
00126           fprintf(pfile,"CONSERVED vars at time= %e, level= %i, domain= %i\n",
00127             pGrid->time,nl,nd);
00128         } else if(strcmp(pOut->out,"prim") == 0) {
00129           fprintf(pfile,"PRIMITIVE vars at time= %e, level= %i, domain= %i\n",
00130             pGrid->time,nl,nd);
00131         }
00132 
00133 /*  3. File format */
00134 
00135         fprintf(pfile,"BINARY\n");
00136 
00137 /*  4. Dataset structure */
00138 
00139 /* Set the Grid origin */
00140 
00141         x1 = pGrid->MinX[0];
00142         x2 = pGrid->MinX[1];
00143         x3 = pGrid->MinX[2];
00144 
00145         fprintf(pfile,"DATASET STRUCTURED_POINTS\n");
00146         if (pGrid->Nx[1] == 1) {
00147           fprintf(pfile,"DIMENSIONS %d %d %d\n",iu-il+2,1,1);
00148         } else {
00149           if (pGrid->Nx[2] == 1) {
00150             fprintf(pfile,"DIMENSIONS %d %d %d\n",iu-il+2,ju-jl+2,1);
00151           } else {
00152             fprintf(pfile,"DIMENSIONS %d %d %d\n",iu-il+2,ju-jl+2,ku-kl+2);
00153           }
00154         }
00155         fprintf(pfile,"ORIGIN %e %e %e \n",x1,x2,x3);
00156         fprintf(pfile,"SPACING %e %e %e \n",pGrid->dx1,pGrid->dx2,pGrid->dx3);
00157 
00158 /*  5. Data  */
00159 
00160         fprintf(pfile,"CELL_DATA %d \n", (iu-il+1)*(ju-jl+1)*(ku-kl+1));
00161 
00162 /* Write density */
00163 
00164         fprintf(pfile,"SCALARS density float\n");
00165         fprintf(pfile,"LOOKUP_TABLE default\n");
00166         for (k=kl; k<=ku; k++) {
00167           for (j=jl; j<=ju; j++) {
00168             for (i=il; i<=iu; i++) {
00169               if (strcmp(pOut->out,"cons") == 0){
00170                 data[i-il] = (float)pGrid->U[k][j][i].d;
00171               } else if(strcmp(pOut->out,"prim") == 0) {
00172                 data[i-il] = (float)W[k-kl][j-jl][i-il].d;
00173               }
00174             }
00175             if(!big_end) ath_bswap(data,sizeof(float),iu-il+1);
00176             fwrite(data,sizeof(float),(size_t)ndata0,pfile);
00177           }
00178         }
00179 
00180 /* Write momentum or velocity */
00181 
00182         if (strcmp(pOut->out,"cons") == 0){
00183           fprintf(pfile,"\nVECTORS momentum float\n");
00184         } else if(strcmp(pOut->out,"prim") == 0) {
00185           fprintf(pfile,"\nVECTORS velocity float\n");
00186         }
00187         for (k=kl; k<=ku; k++) {
00188           for (j=jl; j<=ju; j++) {
00189             for (i=il; i<=iu; i++) {
00190               if (strcmp(pOut->out,"cons") == 0){
00191                 data[3*(i-il)  ] = (float)pGrid->U[k][j][i].M1;
00192                 data[3*(i-il)+1] = (float)pGrid->U[k][j][i].M2;
00193                 data[3*(i-il)+2] = (float)pGrid->U[k][j][i].M3;
00194               } else if(strcmp(pOut->out,"prim") == 0) {
00195                 data[3*(i-il)  ] = (float)W[k-kl][j-jl][i-il].V1;
00196                 data[3*(i-il)+1] = (float)W[k-kl][j-jl][i-il].V2;
00197                 data[3*(i-il)+2] = (float)W[k-kl][j-jl][i-il].V3;
00198               }
00199             }
00200             if(!big_end) ath_bswap(data,sizeof(float),3*(iu-il+1));
00201             fwrite(data,sizeof(float),(size_t)(3*ndata0),pfile);
00202           }
00203         }
00204 
00205 /* Write total energy or pressure */
00206 
00207 #ifndef BAROTROPIC
00208         if (strcmp(pOut->out,"cons") == 0){
00209           fprintf(pfile,"\nSCALARS total_energy float\n");
00210         } else if(strcmp(pOut->out,"prim") == 0) {
00211           fprintf(pfile,"\nSCALARS pressure float\n");
00212         }
00213         fprintf(pfile,"LOOKUP_TABLE default\n");
00214         for (k=kl; k<=ku; k++) {
00215           for (j=jl; j<=ju; j++) {
00216             for (i=il; i<=iu; i++) {
00217               if (strcmp(pOut->out,"cons") == 0){
00218                 data[i-il] = (float)pGrid->U[k][j][i].E;
00219               } else if(strcmp(pOut->out,"prim") == 0) {
00220                 data[i-il] = (float)W[k-kl][j-jl][i-il].P;
00221               }
00222             }
00223             if(!big_end) ath_bswap(data,sizeof(float),iu-il+1);
00224             fwrite(data,sizeof(float),(size_t)ndata0,pfile);
00225           }
00226         }
00227 #endif
00228 
00229 /* Write cell centered B */
00230 
00231 #ifdef MHD
00232         fprintf(pfile,"\nVECTORS cell_centered_B float\n");
00233         for (k=kl; k<=ku; k++) {
00234           for (j=jl; j<=ju; j++) {
00235             for (i=il; i<=iu; i++) {
00236               data[3*(i-il)] = (float)pGrid->U[k][j][i].B1c;
00237               data[3*(i-il)+1] = (float)pGrid->U[k][j][i].B2c;
00238               data[3*(i-il)+2] = (float)pGrid->U[k][j][i].B3c;
00239             }
00240             if(!big_end) ath_bswap(data,sizeof(float),3*(iu-il+1));
00241             fwrite(data,sizeof(float),(size_t)(3*ndata0),pfile);
00242           }
00243         }
00244 #endif
00245 
00246 /* Write gravitational potential */
00247 
00248 #ifdef SELF_GRAVITY
00249         fprintf(pfile,"\nSCALARS gravitational_potential float\n");
00250         fprintf(pfile,"LOOKUP_TABLE default\n");
00251         for (k=kl; k<=ku; k++) {
00252           for (j=jl; j<=ju; j++) {
00253             for (i=il; i<=iu; i++) {
00254               data[i-il] = (float)pGrid->Phi[k][j][i];
00255             }
00256             if(!big_end) ath_bswap(data,sizeof(float),iu-il+1);
00257             fwrite(data,sizeof(float),(size_t)ndata0,pfile);
00258           }
00259         }
00260 #endif
00261 
00262 /* Write binned particle grid */
00263 
00264 #ifdef PARTICLES
00265         if (pOut->out_pargrid) {
00266           fprintf(pfile,"\nSCALARS particle_density float\n");
00267           fprintf(pfile,"LOOKUP_TABLE default\n");
00268           for (k=kl; k<=ku; k++) {
00269             for (j=jl; j<=ju; j++) {
00270               for (i=il; i<=iu; i++) {
00271                 data[i-il] = pGrid->Coup[k][j][i].grid_d;
00272               }
00273               if(!big_end) ath_bswap(data,sizeof(float),iu-il+1);
00274               fwrite(data,sizeof(float),(size_t)ndata0,pfile);
00275             }
00276           }
00277           fprintf(pfile,"\nVECTORS particle_momentum float\n");
00278           for (k=kl; k<=ku; k++) {
00279             for (j=jl; j<=ju; j++) {
00280               for (i=il; i<=iu; i++) {
00281                 data[3*(i-il)] = pGrid->Coup[k][j][i].grid_v1;
00282                 data[3*(i-il)+1] = pGrid->Coup[k][j][i].grid_v2;
00283                 data[3*(i-il)+2] = pGrid->Coup[k][j][i].grid_v3;
00284               }
00285               if(!big_end) ath_bswap(data,sizeof(float),3*(iu-il+1));
00286               fwrite(data,sizeof(float),(size_t)(3*ndata0),pfile);
00287             }
00288           }
00289         }
00290 #endif
00291 
00292 /* Write passive scalars */
00293 
00294 #if (NSCALARS > 0)
00295         for (n=0; n<NSCALARS; n++){
00296           if (strcmp(pOut->out,"cons") == 0){
00297             fprintf(pfile,"\nSCALARS scalar[%d] float\n",n);
00298           } else if(strcmp(pOut->out,"prim") == 0) {
00299             fprintf(pfile,"\nSCALARS specific_scalar[%d] float\n",n);
00300           }
00301           fprintf(pfile,"LOOKUP_TABLE default\n");
00302           for (k=kl; k<=ku; k++) {
00303             for (j=jl; j<=ju; j++) {
00304               for (i=il; i<=iu; i++) {
00305                 if (strcmp(pOut->out,"cons") == 0){
00306                   data[i-il] = (float)pGrid->U[k][j][i].s[n];
00307                 } else if(strcmp(pOut->out,"prim") == 0) {
00308                   data[i-il] = (float)W[k-kl][j-jl][i-il].r[n];
00309                 }
00310               }
00311               if(!big_end) ath_bswap(data,sizeof(float),iu-il+1);
00312               fwrite(data,sizeof(float),(size_t)ndata0,pfile);
00313             }
00314           }
00315         }
00316 #endif
00317 
00318 /* close file and free memory */
00319 
00320         fclose(pfile);
00321         free(data);
00322         if(strcmp(pOut->out,"prim") == 0) free_3d_array(W);
00323       }}
00324     }
00325   }
00326   return;
00327 }