00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
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
00028
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
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
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
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
00073
00074 ndata[0] = iu-il+1;
00075 ndata[1] = ju-jl+1;
00076 ndata[2] = ku-kl+1;
00077
00078
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
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
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
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
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;
00146 #endif
00147 fwrite(dat,sizeof(float),2,p_binfile);
00148
00149
00150
00151 dat[0] = (float)pGrid->time;
00152 dat[1] = (float)pGrid->dt;
00153 fwrite(dat,sizeof(float),2,p_binfile);
00154
00155
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
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
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
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 }