00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
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
00026
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
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;
00040 double x1, x2, x3;
00041 #if (NSCALARS > 0)
00042 int n;
00043 #endif
00044
00045
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
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
00074
00075 ndata0 = iu-il+1;
00076 ndata1 = ju-jl+1;
00077 ndata2 = ku-kl+1;
00078
00079
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
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
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
00119
00120
00121 fprintf(pfile,"# vtk DataFile Version 2.0\n");
00122
00123
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
00134
00135 fprintf(pfile,"BINARY\n");
00136
00137
00138
00139
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
00159
00160 fprintf(pfile,"CELL_DATA %d \n", (iu-il+1)*(ju-jl+1)*(ku-kl+1));
00161
00162
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
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
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
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
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
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
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
00319
00320 fclose(pfile);
00321 free(data);
00322 if(strcmp(pOut->out,"prim") == 0) free_3d_array(W);
00323 }}
00324 }
00325 }
00326 return;
00327 }