00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
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
00034
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
00051 if(pOut->dat_fmt == NULL){
00052 sprintf(fmt," %%12.8e");
00053 }
00054 else{
00055 sprintf(fmt," %s",pOut->dat_fmt);
00056 }
00057
00058
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
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
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
00086 if((pfile = fopen(fname,"w")) == NULL){
00087 ath_error("[dump_tab]: Unable to open ppm file %s\n",fname);
00088 }
00089
00090
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
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
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
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
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
00172 #ifndef BAROTROPIC
00173 fprintf(pfile," [%d]=E",col_cnt);
00174 col_cnt++;
00175 #endif
00176
00177
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
00186
00187
00188 #ifdef SELF_GRAVITY
00189 fprintf(pfile," [%d]=Phi",col_cnt);
00190 col_cnt++;
00191 #endif
00192
00193
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
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
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
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
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 }
00270 }
00271
00272 fclose(pfile);
00273
00274 return;
00275 }
00276
00277
00278
00279
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
00297 if(pOut->dat_fmt == NULL){
00298 sprintf(fmt," %%12.8e");
00299 }
00300 else{
00301 sprintf(fmt," %s",pOut->dat_fmt);
00302 }
00303
00304
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
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
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
00332 if((pfile = fopen(fname,"w")) == NULL){
00333 ath_error("[dump_tab]: Unable to open ppm file %s\n",fname);
00334 }
00335
00336
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
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
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
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
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
00417 #ifndef BAROTROPIC
00418 fprintf(pfile," [%d]=P",col_cnt);
00419 col_cnt++;
00420 #endif
00421
00422
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
00431
00432
00433 #ifdef SELF_GRAVITY
00434 fprintf(pfile," [%d]=Phi",col_cnt);
00435 col_cnt++;
00436 #endif
00437
00438
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
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
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
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
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 }
00519 }
00520
00521 fclose(pfile);
00522
00523 return;
00524 }