00001 #include "copyright.h"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083 #ifndef __POWERPC__
00084 #define HAVE_DLFCN
00085 #endif
00086
00087 #include <ctype.h>
00088 #include <math.h>
00089 #include <stdio.h>
00090 #include <stdlib.h>
00091 #include <string.h>
00092 #ifdef HAVE_DLFCN
00093 #include <dlfcn.h>
00094 #endif
00095
00096 #include "defs.h"
00097 #include "athena.h"
00098 #include "globals.h"
00099 #include "palette.h"
00100 #include "prototypes.h"
00101 #include "particles/prototypes.h"
00102
00103 #define MAXOUT_DEFAULT 10
00104
00105 static int out_count = 0;
00106 static OutputS *OutArray = NULL;
00107 static OutputS rst_out;
00108 static int rst_flag = 0;
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119 Real expr_d (const GridS *pG, const int i, const int j, const int k);
00120 Real expr_M1 (const GridS *pG, const int i, const int j, const int k);
00121 Real expr_M2 (const GridS *pG, const int i, const int j, const int k);
00122 Real expr_M3 (const GridS *pG, const int i, const int j, const int k);
00123 Real expr_E (const GridS *pG, const int i, const int j, const int k);
00124 Real expr_B1c(const GridS *pG, const int i, const int j, const int k);
00125 Real expr_B2c(const GridS *pG, const int i, const int j, const int k);
00126 Real expr_B3c(const GridS *pG, const int i, const int j, const int k);
00127 Real expr_ME (const GridS *pG, const int i, const int j, const int k);
00128 Real expr_V1 (const GridS *pG, const int i, const int j, const int k);
00129 Real expr_V2 (const GridS *pG, const int i, const int j, const int k);
00130 Real expr_V3 (const GridS *pG, const int i, const int j, const int k);
00131 Real expr_P (const GridS *pG, const int i, const int j, const int k);
00132 Real expr_cs2(const GridS *pG, const int i, const int j, const int k);
00133 Real expr_S (const GridS *pG, const int i, const int j, const int k);
00134 #ifdef SPECIAL_RELATIVITY
00135 Real expr_G (const GridS *pG, const int i, const int j, const int k);
00136 #endif
00137 #ifdef PARTICLES
00138
00139
00140 extern Real expr_dpar (const GridS *pG, const int i, const int j, const int k);
00141
00142
00143 extern Real expr_M1par(const GridS *pG, const int i, const int j, const int k);
00144
00145
00146 extern Real expr_M2par(const GridS *pG, const int i, const int j, const int k);
00147
00148
00149 extern Real expr_M3par(const GridS *pG, const int i, const int j, const int k);
00150
00151
00152 extern Real expr_V1par(const GridS *pG, const int i, const int j, const int k);
00153
00154
00155 extern Real expr_V2par(const GridS *pG, const int i, const int j, const int k);
00156
00157
00158 extern Real expr_V3par(const GridS *pG, const int i, const int j, const int k);
00159 int check_particle_binning(char *out);
00160 #endif
00161 static ConsFun_t getexpr(const int n, const char *expr);
00162 static void free_output(OutputS *pout);
00163 static void parse_slice(char *block, char *axname, Real *l, Real *u, int *flag);
00164 float *getRGB(char *name);
00165
00166
00167
00168
00169
00170
00171 void init_output(MeshS *pM)
00172 {
00173 int i,j,outn,maxout,nl,nd;
00174 char block[80], *fmt, defid[10];
00175 OutputS new_out;
00176 int usr_expr_flag;
00177
00178 maxout = par_geti_def("job","maxout",MAXOUT_DEFAULT);
00179
00180
00181
00182 if((OutArray = malloc(maxout*sizeof(OutputS))) == NULL){
00183 ath_error("[init_output]: Error allocating output array\n");
00184 }
00185
00186
00187
00188
00189 for (outn=1; outn<=maxout; outn++) {
00190
00191 sprintf(block,"output%d",outn);
00192
00193
00194
00195 if((par_exist(block,"out_fmt") == 0) && (par_exist(block,"name") == 0)){
00196 ath_perr(-1,"[init_output]: neither %s/out_fmt, nor %s/name exist\n",
00197 block, block);
00198 continue;
00199 }
00200
00201
00202 memset(&new_out,0,sizeof(OutputS));
00203
00204
00205 new_out.t = par_getd_def(block,"time",pM->time);
00206 new_out.num = par_geti_def(block,"num",0);
00207
00208 new_out.dt = par_getd(block,"dt");
00209 new_out.n = outn;
00210
00211
00212 nl = new_out.nlevel = par_geti_def(block,"level",-1);
00213 nd = new_out.ndomain = par_geti_def(block,"domain",-1);
00214
00215 if (par_exist(block,"dat_fmt")) new_out.dat_fmt = par_gets(block,"dat_fmt");
00216
00217
00218
00219 sprintf(defid,"out%d",outn);
00220 new_out.id = par_gets_def(block,"id",defid);
00221
00222 if(par_exist(block,"out_fmt"))
00223 fmt = new_out.out_fmt = par_gets(block,"out_fmt");
00224
00225
00226
00227
00228
00229 new_out.out = par_gets_def(block,"out","cons");
00230
00231 #ifdef PARTICLES
00232
00233 new_out.out_pargrid = par_geti_def(block,"pargrid",
00234 check_particle_binning(new_out.out));
00235 if ((new_out.out_pargrid < 0) || (new_out.out_pargrid >1)) {
00236 ath_perr(-1,"[init_output]: %s/pargrid must be 0 or 1\n",
00237 block, block);
00238 continue;
00239 }
00240
00241
00242
00243 if(par_exist(block,"par_prop")) {
00244 new_out.par_prop = get_usr_par_prop(par_gets(block,"par_prop"));
00245 if (new_out.par_prop == NULL) {
00246 ath_pout(0,"[init_output]: Particle selection function not found! \
00247 Now use the default one.\n");
00248 new_out.par_prop = property_all;
00249 }
00250 }
00251 else
00252 new_out.par_prop = property_all;
00253 #endif
00254
00255
00256
00257 if(strcmp(new_out.out,"cons") == 0){
00258
00259 if(par_exist(block,"name")){
00260
00261 char *name = par_gets(block,"name");
00262
00263 new_out.out_fun = get_usr_out_fun(name);
00264 if(new_out.out_fun == NULL){
00265 free_output(&new_out);
00266 ath_error("Unsupported output named %s in %s/out_fmt=%s\n",
00267 name,block,fmt);
00268 }
00269 free(name); name = NULL;
00270 goto add_it;
00271 }
00272 else if (strcmp(fmt,"bin")==0){
00273 new_out.out_fun = dump_binary;
00274 #ifdef PARTICLES
00275 new_out.out_pargrid = 1;
00276 #endif
00277 #ifdef PARTICLES
00278 new_out.out_pargrid = 1;
00279 #endif
00280 goto add_it;
00281 }
00282 else if (strcmp(fmt,"hst")==0){
00283 new_out.out_fun = dump_history;
00284 goto add_it;
00285 }
00286 #ifdef PARTICLES
00287 else if (strcmp(fmt,"phst")==0){
00288 new_out.fun = dump_particle_history;
00289 goto add_it;
00290 }
00291 #endif
00292 else if (strcmp(fmt,"tab")==0){
00293 new_out.out_fun = dump_tab_cons;
00294 #ifdef PARTICLES
00295 new_out.out_pargrid = 1;
00296 #endif
00297 goto add_it;
00298 }
00299 else if (strcmp(fmt,"rst")==0){
00300 new_out.res_fun = dump_restart;
00301 rst_flag = 1;
00302 rst_out = new_out;
00303 ath_pout(0,"Added out%d\n",outn);
00304 continue;
00305 }
00306 else if (strcmp(fmt,"vtk")==0){
00307 new_out.out_fun = dump_vtk;
00308 #ifdef PARTICLES
00309 new_out.out_pargrid = 1;
00310 #endif
00311 goto add_it;
00312 }
00313 #ifdef PARTICLES
00314 else if (strcmp(fmt,"lis")==0){
00315 new_out.fun = dump_particle_binary;
00316 goto add_it;
00317 }
00318 #endif
00319 else{
00320 ath_error("Unsupported dump mode for %s/out_fmt=%s for out=cons\n",
00321 block,fmt);
00322 }
00323 }
00324
00325
00326
00327 if(strcmp(new_out.out,"prim") == 0){
00328
00329 if(par_exist(block,"name")){
00330
00331 char *name = par_gets(block,"name");
00332
00333 new_out.out_fun = get_usr_out_fun(name);
00334 if(new_out.out_fun == NULL){
00335 free_output(&new_out);
00336 ath_error("Unsupported output named %s in %s/out_fmt=%s\n",
00337 name,block,fmt);
00338 }
00339 free(name); name = NULL;
00340 goto add_it;
00341 }
00342 else if (strcmp(fmt,"bin")==0){
00343 new_out.out_fun = dump_binary;
00344 goto add_it;
00345 }
00346 else if (strcmp(fmt,"tab")==0){
00347 new_out.out_fun = dump_tab_prim;
00348 #ifdef PARTICLES
00349 new_out.out_pargrid = 1;
00350 #endif
00351 goto add_it;
00352 }
00353 else if (strcmp(fmt,"vtk")==0){
00354 new_out.out_fun = dump_vtk;
00355 goto add_it;
00356 }
00357 else{
00358 ath_error("Unsupported dump mode for %s/out_fmt=%s for out=prim\n",
00359 block,fmt);
00360 }
00361 }
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371 if(par_exist(block,"usr_expr_flag"))
00372 usr_expr_flag = par_geti(block,"usr_expr_flag");
00373 else
00374 usr_expr_flag = 0;
00375
00376
00377 if(usr_expr_flag)
00378 new_out.expr = get_usr_expr(new_out.out);
00379 else
00380 new_out.expr = getexpr(outn, new_out.out);
00381
00382 if (new_out.expr == NULL) {
00383 ath_perr(-1,"Could not parse expression %s, skipping it\n",
00384 new_out.out);
00385 free_output(&new_out);
00386 continue;
00387 }
00388
00389
00390 new_out.ndim = 1;
00391 for (i=1; i<3; i++) if (pM->Nx[i]>1) new_out.ndim++;
00392
00393 new_out.x1l = pM->RootMinX[0];
00394 new_out.x1u = pM->RootMaxX[0];
00395 new_out.reduce_x1 = 0;
00396 parse_slice(block,"x1",&new_out.x1l,&new_out.x1u,&new_out.reduce_x1);
00397 if (new_out.reduce_x1 != 0) new_out.ndim--;
00398
00399 new_out.x2l = pM->RootMinX[1];
00400 new_out.x2u = pM->RootMaxX[1];
00401 new_out.reduce_x2 = 0;
00402 parse_slice(block,"x2",&new_out.x2l,&new_out.x2u,&new_out.reduce_x2);
00403 if (pM->Nx[1] > 1 && new_out.reduce_x2 != 0) new_out.ndim--;
00404
00405 new_out.x3l = pM->RootMinX[2];
00406 new_out.x3u = pM->RootMaxX[2];
00407 new_out.reduce_x3 = 0;
00408 parse_slice(block,"x3",&new_out.x3l,&new_out.x3u,&new_out.reduce_x3);
00409 if (pM->Nx[2] > 1 && new_out.reduce_x3 != 0) new_out.ndim--;
00410
00411 if (new_out.ndim <= 0) ath_error("Too many slices specified in %s\n",block);
00412
00413
00414 if(par_exist(block,"dmin") != 0){
00415 new_out.sdmin = 1;
00416 new_out.dmin = par_getd(block,"dmin");
00417 }
00418 new_out.gmin = (HUGE_NUMBER);
00419
00420 if(par_exist(block,"dmax") != 0){
00421 new_out.sdmax = 1;
00422 new_out.dmax = par_getd(block,"dmax");
00423 }
00424 new_out.gmax = -1.0*(HUGE_NUMBER);
00425
00426
00427 if (strcmp(fmt,"ppm") == 0) {
00428 new_out.palette = par_gets_def(block,"palette","rainbow");
00429
00430 new_out.rgb = getRGB(new_out.palette);
00431 if ( (new_out.der = (float *) malloc(3*256*sizeof(float))) == NULL) {
00432 free_output(&new_out);
00433 ath_error("[init_output]: malloc returned a NULL pointer\n");
00434 }
00435 for(j=0; j<3; j++)
00436 for (i=0; i<255; i++)
00437 new_out.der[3*i+j] = new_out.rgb[3*(i+1)+j] - new_out.rgb[3*i+j];
00438 }
00439
00440
00441
00442
00443
00444
00445 if(par_exist(block,"name")){
00446
00447 char *name = par_gets(block,"name");
00448
00449 new_out.out_fun = get_usr_out_fun(name);
00450 if(new_out.out_fun == NULL){
00451 free_output(&new_out);
00452 ath_error("Unsupported output named %s in %s/out_fmt=%s\n",
00453 name,block,fmt);
00454 }
00455 free(name); name = NULL;
00456 }
00457 else if (strcmp(fmt,"pdf")==0)
00458 new_out.out_fun = output_pdf;
00459 else if (strcmp(fmt,"pgm")==0)
00460 new_out.out_fun = output_pgm;
00461 else if (strcmp(fmt,"ppm")==0)
00462 new_out.out_fun = output_ppm;
00463 else if (strcmp(fmt,"vtk")==0)
00464 new_out.out_fun = output_vtk;
00465 else if (strcmp(fmt,"tab")==0)
00466 new_out.out_fun = output_tab;
00467 else {
00468
00469 free_output(&new_out);
00470 ath_error("Unsupported %s/out_fmt=%s\n",block,fmt);
00471 }
00472
00473 add_it:
00474
00475
00476
00477 ath_pout(1,"OUTPUT: %d %d %s %s [%g : %g]\n",
00478 new_out.n, new_out.ndim, new_out.out_fmt,
00479 new_out.out, new_out.dmin, new_out.dmax);
00480
00481 OutArray[out_count] = new_out;
00482 out_count++;
00483 ath_pout(0,"Added out%d\n",outn);
00484
00485 }
00486
00487 }
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498 void data_output(MeshS *pM, const int flag)
00499 {
00500 int n;
00501 int dump_flag[MAXOUT_DEFAULT+1];
00502 char block[80];
00503
00504
00505
00506
00507 for (n=0; n<out_count; n++) {
00508 dump_flag[n] = flag;
00509 if (pM->time >= OutArray[n].t) {
00510 OutArray[n].t += OutArray[n].dt;
00511 dump_flag[n] = 1;
00512 }
00513 }
00514
00515
00516
00517 if(rst_flag){
00518 dump_flag[out_count] = flag;
00519 if(pM->time >= rst_out.t){
00520 rst_out.t += rst_out.dt;
00521 dump_flag[out_count] = 1;
00522 }
00523
00524 if(dump_flag[out_count] != 0){
00525
00526 for(n=0; n<out_count; n++){
00527
00528 if(OutArray[n].n > 0){
00529 sprintf(block,"output%d",OutArray[n].n);
00530 if (dump_flag[n] != 0) {
00531
00532
00533 par_seti(block,"num","%d",OutArray[n].num+1,"Next Output Number");
00534 } else {
00535 par_seti(block,"num","%d",OutArray[n].num,"Next Output Number");
00536 }
00537 par_setd(block,"time","%.15e",OutArray[n].t,"Next Output Time");
00538 }
00539 }
00540
00541 sprintf(block,"output%d",rst_out.n);
00542 par_seti(block,"num","%d",rst_out.num+1,"Next Output Number");
00543 par_setd(block,"time","%.15e",rst_out.t,"Next Output Time");
00544
00545
00546 (*(rst_out.res_fun))(pM,&(rst_out));
00547
00548 rst_out.num++;
00549 }
00550 }
00551
00552
00553
00554 for (n=0; n<out_count; n++) {
00555 if(dump_flag[n] != 0) {
00556
00557 #ifdef PARTICLES
00558 if (OutArray[n].out_pargrid == 1)
00559 particle_to_grid(pM, OutArray[n].par_prop);
00560 #endif
00561 (*OutArray[n].out_fun)(pM,&(OutArray[n]));
00562
00563 OutArray[n].num++;
00564
00565 }
00566 }
00567
00568 return;
00569 }
00570
00571
00572
00573
00574
00575
00576 void data_output_destruct(void)
00577 {
00578 int i;
00579 double global_min, global_max;
00580 #ifdef MPI_PARALLEL
00581 int ierr;
00582 #endif
00583
00584 for (i=0; i<out_count; i++) {
00585
00586
00587
00588 if (OutArray[i].out != NULL){
00589 if((strcmp(OutArray[i].out,"cons") != 0) &&
00590 (strcmp(OutArray[i].out,"prim") != 0)){
00591
00592 #ifdef MPI_PARALLEL
00593 ierr = MPI_Allreduce(&OutArray[i].gmin, &global_min, 1, MPI_DOUBLE,
00594 MPI_MIN, MPI_COMM_WORLD);
00595 ierr = MPI_Allreduce(&OutArray[i].gmax, &global_max, 1, MPI_DOUBLE,
00596 MPI_MAX, MPI_COMM_WORLD);
00597 #else
00598 global_min = OutArray[i].gmin;
00599 global_max = OutArray[i].gmax;
00600 #endif
00601 ath_pout(0,"Global min/max for %s: %g %g\n",OutArray[i].out,
00602 global_min, global_max);
00603 }
00604
00605 free(OutArray[i].out);
00606 }
00607 if (OutArray[i].out_fmt != NULL) free(OutArray[i].out_fmt);
00608 if (OutArray[i].dat_fmt != NULL) free(OutArray[i].dat_fmt);
00609 if (OutArray[i].id != NULL) free(OutArray[i].id);
00610 }
00611
00612 if(rst_flag){
00613 if (rst_out.out != NULL) free(rst_out.out);
00614 if (rst_out.out_fmt != NULL) free(rst_out.out_fmt);
00615 if (rst_out.dat_fmt != NULL) free(rst_out.dat_fmt);
00616 if (rst_out.id != NULL) free(rst_out.id);
00617 }
00618
00619 if (OutArray != NULL) {
00620 free(OutArray);
00621 OutArray = NULL;
00622 out_count = 0;
00623 }
00624
00625 return;
00626 }
00627
00628
00629
00630
00631
00632
00633
00634
00635
00636 Real ***OutData3(GridS *pgrid, OutputS *pout, int *Nx1, int *Nx2, int *Nx3)
00637 {
00638 Real ***data;
00639 int i,j,k,il,jl,kl,iu,ju,ku;
00640
00641 if (pout->ndim != 3) ath_error("[OutData3] <output%d> %s is %d-D, not 3-D\n",
00642 pout->n,pout->out, pout->ndim);
00643
00644 #ifdef WRITE_GHOST_CELLS
00645 if(pgrid->Nx[0] > 1){
00646 il = pgrid->is - nghost;
00647 iu = pgrid->ie + nghost;
00648 } else{
00649 il = pgrid->is;
00650 iu = pgrid->ie;
00651 }
00652
00653 if(pgrid->Nx[1] > 1){
00654 jl = pgrid->js - nghost;
00655 ju = pgrid->je + nghost;
00656 } else{
00657 jl = pgrid->js;
00658 ju = pgrid->je;
00659 }
00660
00661 if(pgrid->Nx[2] > 1){
00662 kl = pgrid->ks - nghost;
00663 ku = pgrid->ke + nghost;
00664 } else{
00665 kl = pgrid->ks;
00666 ku = pgrid->ke;
00667 }
00668 #else
00669 il = pgrid->is;
00670 iu = pgrid->ie;
00671 jl = pgrid->js;
00672 ju = pgrid->je;
00673 kl = pgrid->ks;
00674 ku = pgrid->ke;
00675 #endif
00676 *Nx1 = iu-il+1;
00677 *Nx2 = ju-jl+1;
00678 *Nx3 = ku-kl+1;
00679
00680 data = (Real***) calloc_3d_array(*Nx3,*Nx2,*Nx1,sizeof(Real));
00681 if (data == NULL) ath_error("[OutData3] Error creating 3D data array\n");
00682 for (k=0; k<*Nx3; k++)
00683 for (j=0; j<*Nx2; j++)
00684 for (i=0; i<*Nx1; i++)
00685 data[k][j][i] = (*pout->expr)(pgrid,i+il,j+jl,k+kl);
00686 return data;
00687 }
00688
00689
00690
00691
00692
00693
00694
00695
00696
00697
00698
00699 Real **OutData2(GridS *pgrid, OutputS *pout, int *Nx1, int *Nx2)
00700 {
00701 Real **data;
00702 Real factor,x1fc,x2fc,x3fc;
00703 int Nx3;
00704 int i,j,k,il,jl,kl,iu,ju,ku;
00705 int istart,iend,jstart,jend,kstart,kend;
00706
00707 if (pout->ndim != 2) ath_error("[OutData2] <output%d> %s is %d-D, not 2-D\n",
00708 pout->n,pout->out, pout->ndim);
00709
00710 #ifdef WRITE_GHOST_CELLS
00711 if(pgrid->Nx[0] > 1){
00712 il = pgrid->is - nghost;
00713 iu = pgrid->ie + nghost;
00714 } else{
00715 il = pgrid->is;
00716 iu = pgrid->ie;
00717 }
00718
00719 if(pgrid->Nx[1] > 1){
00720 jl = pgrid->js - nghost;
00721 ju = pgrid->je + nghost;
00722 } else{
00723 jl = pgrid->js;
00724 ju = pgrid->je;
00725 }
00726
00727 if(pgrid->Nx[2] > 1){
00728 kl = pgrid->ks - nghost;
00729 ku = pgrid->ke + nghost;
00730 } else{
00731 kl = pgrid->ks;
00732 ku = pgrid->ke;
00733 }
00734 #else
00735 il = pgrid->is;
00736 iu = pgrid->ie;
00737 jl = pgrid->js;
00738 ju = pgrid->je;
00739 kl = pgrid->ks;
00740 ku = pgrid->ke;
00741 #endif
00742 *Nx1 = iu-il+1;
00743 *Nx2 = ju-jl+1;
00744 Nx3 = ku-kl+1;
00745
00746
00747
00748 if (pgrid->Nx[2] == 1) {
00749 data = (Real**) calloc_2d_array(*Nx2,*Nx1,sizeof(Real));
00750 if (data == NULL) ath_error("[OutData2] Error creating 2D data array\n");
00751 for (j=0; j<*Nx2; j++) {
00752 for (i=0; i<*Nx1; i++) {
00753 data[j][i] = (*pout->expr)(pgrid,i+il,j+jl,kl);
00754 }
00755 }
00756 return data;
00757 }
00758
00759
00760
00761
00762 if (pout->reduce_x3 != 0) {
00763 if (pout->x3u < pgrid->MinX[2] || pout->x3l >= pgrid->MaxX[2]) return NULL;
00764
00765
00766 k=kl+1;
00767 fc_pos(pgrid,il,jl,k,&x1fc,&x2fc,&x3fc);
00768 while (pout->x3l >= x3fc) {
00769 k++;
00770 fc_pos(pgrid,il,jl,k,&x1fc,&x2fc,&x3fc);
00771 }
00772 kstart = k-1;
00773
00774 k=ku;
00775 fc_pos(pgrid,il,jl,k,&x1fc,&x2fc,&x3fc);
00776 while (pout->x3u < x3fc) {
00777 k--;
00778 fc_pos(pgrid,il,jl,k,&x1fc,&x2fc,&x3fc);
00779 }
00780 kend = k;
00781
00782
00783 data = (Real**) calloc_2d_array(*Nx2,*Nx1,sizeof(Real));
00784 if (data == NULL) ath_error("[OutData2] Error creating 2D data array\n");
00785 factor = 1.0/(kend - kstart + 1);
00786 for (j=0; j<*Nx2; j++) {
00787 for (i=0; i<*Nx1; i++) {
00788 data[j][i] = 0.0;
00789 for (k=kstart; k<=kend; k++)
00790 data[j][i] += (*pout->expr)(pgrid,i+il,j+jl,k+kl);
00791 data[j][i] *= factor;
00792 }
00793 }
00794
00795
00796 } else if (pout->reduce_x2 != 0) {
00797 if (pout->x2u < pgrid->MinX[1] || pout->x2l >= pgrid->MaxX[1]) return NULL;
00798
00799
00800 j=jl+1;
00801 fc_pos(pgrid,il,j,kl,&x1fc,&x2fc,&x3fc);
00802 while (pout->x2l >= x2fc) {
00803 j++;
00804 fc_pos(pgrid,il,j,kl,&x1fc,&x2fc,&x3fc);
00805 }
00806 jstart = j-1;
00807
00808 j=ju;
00809 fc_pos(pgrid,il,j,kl,&x1fc,&x2fc,&x3fc);
00810 while (pout->x2u < x2fc) {
00811 j--;
00812 fc_pos(pgrid,il,j,kl,&x1fc,&x2fc,&x3fc);
00813 }
00814 jend = j;
00815
00816
00817 data = (Real**) calloc_2d_array(Nx3,*Nx1,sizeof(Real));
00818 if (data == NULL) ath_error("[OutData2] Error creating 2D data array\n");
00819 factor = 1.0/(jend - jstart + 1);
00820 for (k=0; k<Nx3; k++) {
00821 for (i=0; i<*Nx1; i++) {
00822 data[k][i] = 0.0;
00823 for (j=jstart; j<=jend; j++)
00824 data[k][i] += (*pout->expr)(pgrid,i+il,j+jl,k+kl);
00825 data[k][i] *= factor;
00826 }
00827 }
00828 *Nx2 = Nx3;
00829
00830
00831 } else if (pout->reduce_x1 != 0) {
00832 if (pout->x1u < pgrid->MinX[0] || pout->x1l >= pgrid->MaxX[0]) return NULL;
00833
00834
00835 i=il+1;
00836 fc_pos(pgrid,i,jl,kl,&x1fc,&x2fc,&x3fc);
00837 while (pout->x1l >= x1fc) {
00838 i++;
00839 fc_pos(pgrid,i,jl,kl,&x1fc,&x2fc,&x3fc);
00840 }
00841 istart = i-1;
00842
00843 i=iu;
00844 fc_pos(pgrid,i,jl,kl,&x1fc,&x2fc,&x3fc);
00845 while (pout->x1u < x1fc) {
00846 i--;
00847 fc_pos(pgrid,i,jl,kl,&x1fc,&x2fc,&x3fc);
00848 }
00849 iend = i;
00850
00851
00852
00853 data = (Real**) calloc_2d_array(Nx3,*Nx2,sizeof(Real));
00854 if (data == NULL) ath_error("[OutData2] Error creating 2D data array\n");
00855 factor = 1.0/(iend - istart + 1);
00856 for (k=0; k<Nx3; k++) {
00857 for (j=0; j<*Nx2; j++) {
00858 data[k][j] = 0.0;
00859 for (i=istart; i<=iend; i++)
00860 data[k][j] += (*pout->expr)(pgrid,i+il,j+jl,k+kl);
00861 data[k][j] *= factor;
00862 }
00863 }
00864 *Nx1 = *Nx2;
00865 *Nx2 = Nx3;
00866 } else
00867 ath_perr(-1,"[OutData2]: Should not reach here\n");
00868 return data;
00869 }
00870
00871
00872
00873
00874
00875
00876
00877
00878
00879
00880
00881 Real *OutData1(GridS *pgrid, OutputS *pout, int *Nx1)
00882 {
00883 Real *data;
00884 Real factor,x1fc,x2fc,x3fc;
00885 int Nx2, Nx3;
00886 int i,j,k,il,jl,kl,iu,ju,ku;
00887 int istart,iend,jstart,jend,kstart,kend;
00888
00889 if (pout->ndim != 1) ath_error("[OutData1] <output%d> %s is %d-D, not 1-D\n",
00890 pout->n,pout->out, pout->ndim);
00891
00892 #ifdef WRITE_GHOST_CELLS
00893 if(pgrid->Nx[0] > 1){
00894 il = pgrid->is - nghost;
00895 iu = pgrid->ie + nghost;
00896 } else{
00897 il = pgrid->is;
00898 iu = pgrid->ie;
00899 }
00900
00901 if(pgrid->Nx[1] > 1){
00902 jl = pgrid->js - nghost;
00903 ju = pgrid->je + nghost;
00904 } else{
00905 jl = pgrid->js;
00906 ju = pgrid->je;
00907 }
00908
00909 if(pgrid->Nx[2] > 1){
00910 kl = pgrid->ks - nghost;
00911 ku = pgrid->ke + nghost;
00912 } else{
00913 kl = pgrid->ks;
00914 ku = pgrid->ke;
00915 }
00916 #else
00917 il = pgrid->is;
00918 iu = pgrid->ie;
00919 jl = pgrid->js;
00920 ju = pgrid->je;
00921 kl = pgrid->ks;
00922 ku = pgrid->ke;
00923 #endif
00924 *Nx1 = iu-il+1;
00925 Nx2 = ju-jl+1;
00926 Nx3 = ku-kl+1;
00927
00928
00929
00930 if (pgrid->Nx[1] == 1) {
00931 data = (Real*) calloc_1d_array(*Nx1,sizeof(Real));
00932 if (data == NULL) ath_error("[OutData1] Error creating 1D data array\n");
00933 for (i=0; i<*Nx1; i++) {
00934 data[i] = (*pout->expr)(pgrid,i+il,jl,kl);
00935 }
00936 return data;
00937 }
00938
00939
00940
00941
00942 if (pout->reduce_x1 == 0) {
00943 if (pout->x3u < pgrid->MinX[2] || pout->x3l >= pgrid->MaxX[2] ||
00944 pout->x2u < pgrid->MinX[1] || pout->x2l >= pgrid->MaxX[1]) return NULL;
00945
00946
00947 if (pgrid->Nx[2] == 1){
00948 kstart=kl;
00949 kend=ku;
00950 } else {
00951 k=kl+1;
00952 fc_pos(pgrid,il,jl,k,&x1fc,&x2fc,&x3fc);
00953 while (pout->x3l >= x3fc) {
00954 k++;
00955 fc_pos(pgrid,il,jl,k,&x1fc,&x2fc,&x3fc);
00956 }
00957 kstart = k-1;
00958
00959 k=ku;
00960 fc_pos(pgrid,il,jl,k,&x1fc,&x2fc,&x3fc);
00961 while (pout->x3u < x3fc) {
00962 k--;
00963 fc_pos(pgrid,il,jl,k,&x1fc,&x2fc,&x3fc);
00964 }
00965 kend = k;
00966 }
00967
00968
00969 j=jl+1;
00970 fc_pos(pgrid,il,j,kl,&x1fc,&x2fc,&x3fc);
00971 while (pout->x2l >= x2fc) {
00972 j++;
00973 fc_pos(pgrid,il,j,kl,&x1fc,&x2fc,&x3fc);
00974 }
00975 jstart = j-1;
00976
00977 j=ju;
00978 fc_pos(pgrid,il,j,kl,&x1fc,&x2fc,&x3fc);
00979 while (pout->x2u < x2fc) {
00980 j--;
00981 fc_pos(pgrid,il,j,kl,&x1fc,&x2fc,&x3fc);
00982 }
00983 jend = j;
00984
00985
00986 data = (Real*) calloc_1d_array(*Nx1,sizeof(Real));
00987 factor = 1.0/(kend - kstart + 1)/(jend - jstart + 1);
00988 for (i=0; i<*Nx1; i++) {
00989 data[i] = 0.0;
00990 for (k=kstart; k<=kend; k++)
00991 for (j=jstart; j<=jend; j++)
00992 data[i] += (*pout->expr)(pgrid,i+il,j+jl,k+kl);
00993 data[i] *= factor;
00994 }
00995
00996
00997 } else if (pout->reduce_x2 == 0) {
00998 if (pout->x3u < pgrid->MinX[2] || pout->x3l >= pgrid->MaxX[2] ||
00999 pout->x1u < pgrid->MinX[0] || pout->x1l >= pgrid->MaxX[0]) return NULL;
01000
01001
01002 if (pgrid->Nx[2] == 1){
01003 kstart=kl;
01004 kend=ku;
01005 } else {
01006 k=kl+1;
01007 fc_pos(pgrid,il,jl,k,&x1fc,&x2fc,&x3fc);
01008 while (pout->x3l >= x3fc) {
01009 k++;
01010 fc_pos(pgrid,il,jl,k,&x1fc,&x2fc,&x3fc);
01011 }
01012 kstart = k-1;
01013
01014 k=ku;
01015 fc_pos(pgrid,il,jl,k,&x1fc,&x2fc,&x3fc);
01016 while (pout->x3u < x3fc) {
01017 k--;
01018 fc_pos(pgrid,il,jl,k,&x1fc,&x2fc,&x3fc);
01019 }
01020 kend = k;
01021 }
01022
01023
01024 i=il+1;
01025 fc_pos(pgrid,i,jl,kl,&x1fc,&x2fc,&x3fc);
01026 while (pout->x1l >= x1fc) {
01027 i++;
01028 fc_pos(pgrid,i,jl,kl,&x1fc,&x2fc,&x3fc);
01029 }
01030 istart = i-1;
01031
01032 i=iu;
01033 fc_pos(pgrid,i,jl,kl,&x1fc,&x2fc,&x3fc);
01034 while (pout->x1u < x1fc) {
01035 i--;
01036 fc_pos(pgrid,i,jl,kl,&x1fc,&x2fc,&x3fc);
01037 }
01038 iend = i;
01039
01040
01041 data = (Real*) calloc_1d_array(Nx2,sizeof(Real));
01042 factor = 1.0/(kend - kstart + 1)/(iend - istart + 1);
01043 for (j=0; j<Nx2; j++) {
01044 data[j] = 0.0;
01045 for (k=kstart; k<=kend; k++)
01046 for (i=istart; i<=iend; i++)
01047 data[j] += (*pout->expr)(pgrid,i+il,j+jl,k+kl);
01048 data[j] *= factor;
01049 }
01050 *Nx1 = Nx2;
01051
01052
01053 } else if (pout->reduce_x3 == 0) {
01054 if (pout->x2u < pgrid->MinX[1] || pout->x2l >= pgrid->MaxX[1] ||
01055 pout->x1u < pgrid->MinX[0] || pout->x1l >= pgrid->MaxX[0]) return NULL;
01056
01057
01058 j=jl+1;
01059 fc_pos(pgrid,il,j,kl,&x1fc,&x2fc,&x3fc);
01060 while (pout->x2l >= x2fc) {
01061 j++;
01062 fc_pos(pgrid,il,j,kl,&x1fc,&x2fc,&x3fc);
01063 }
01064 jstart = j-1;
01065
01066 j=ju;
01067 fc_pos(pgrid,il,j,kl,&x1fc,&x2fc,&x3fc);
01068 while (pout->x2u < x2fc) {
01069 j--;
01070 fc_pos(pgrid,il,j,kl,&x1fc,&x2fc,&x3fc);
01071 }
01072 jend = j;
01073
01074
01075 i=il+1;
01076 fc_pos(pgrid,i,jl,kl,&x1fc,&x2fc,&x3fc);
01077 while (pout->x1l >= x1fc) {
01078 i++;
01079 fc_pos(pgrid,i,jl,kl,&x1fc,&x2fc,&x3fc);
01080 }
01081 istart = i-1;
01082
01083 i=iu;
01084 fc_pos(pgrid,i,jl,kl,&x1fc,&x2fc,&x3fc);
01085 while (pout->x1u < x1fc) {
01086 i--;
01087 fc_pos(pgrid,i,jl,kl,&x1fc,&x2fc,&x3fc);
01088 }
01089 iend = i;
01090
01091
01092 data = (Real*) calloc_1d_array(Nx3,sizeof(Real));
01093 factor = 1.0/(jend - jstart + 1)/(iend - istart + 1);
01094 for (k=0; k<Nx3; k++) {
01095 data[k] = 0.0;
01096 for (j=jstart; j<=jend; j++)
01097 for (i=istart; i<=iend; i++)
01098 data[k] += (*pout->expr)(pgrid,i+il,j+jl,k+kl);
01099 data[k] *= factor;
01100 }
01101 *Nx1 = Nx3;
01102 } else {
01103 ath_perr(-1,"[OutData1]: Should not reach here\n");
01104 }
01105
01106 return data;
01107 }
01108
01109
01110
01111
01112
01113
01114 Real expr_d(const GridS *pG, const int i, const int j, const int k) {
01115 return pG->U[k][j][i].d;
01116 }
01117
01118
01119 Real expr_M1(const GridS *pG, const int i, const int j, const int k) {
01120 return pG->U[k][j][i].M1;
01121 }
01122
01123
01124 Real expr_M2(const GridS *pG, const int i, const int j, const int k) {
01125 return pG->U[k][j][i].M2;
01126 }
01127
01128
01129 Real expr_M3(const GridS *pG, const int i, const int j, const int k) {
01130 return pG->U[k][j][i].M3;
01131 }
01132 #ifndef BAROTROPIC
01133
01134
01135 Real expr_E(const GridS *pG, const int i, const int j, const int k) {
01136 return pG->U[k][j][i].E;
01137 }
01138 #endif
01139
01140
01141
01142
01143 #ifdef MHD
01144
01145
01146 Real expr_B1c(const GridS *pG, const int i, const int j, const int k) {
01147 return pG->U[k][j][i].B1c;
01148 }
01149
01150
01151 Real expr_B2c(const GridS *pG, const int i, const int j, const int k) {
01152 return pG->U[k][j][i].B2c;
01153 }
01154
01155
01156 Real expr_B3c(const GridS *pG, const int i, const int j, const int k) {
01157 return pG->U[k][j][i].B3c;
01158 }
01159
01160
01161 Real expr_ME(const GridS *pG, const int i, const int j, const int k) {
01162 return 0.5*(pG->U[k][j][i].B1c*pG->U[k][j][i].B1c +
01163 pG->U[k][j][i].B2c*pG->U[k][j][i].B2c +
01164 pG->U[k][j][i].B3c*pG->U[k][j][i].B3c);
01165 }
01166 #endif
01167
01168
01169
01170
01171
01172
01173 Real expr_V1(const GridS *pG, const int i, const int j, const int k) {
01174 return pG->U[k][j][i].M1/pG->U[k][j][i].d;
01175 }
01176
01177
01178 Real expr_V2(const GridS *pG, const int i, const int j, const int k) {
01179 return pG->U[k][j][i].M2/pG->U[k][j][i].d;
01180 }
01181
01182
01183 Real expr_V3(const GridS *pG, const int i, const int j, const int k) {
01184 return pG->U[k][j][i].M3/pG->U[k][j][i].d;
01185 }
01186
01187
01188
01189 Real expr_P(const GridS *pG, const int i, const int j, const int k) {
01190 #ifdef ISOTHERMAL
01191 return pG->U[k][j][i].d*Iso_csound2;
01192 #else
01193 ConsS *gp = &(pG->U[k][j][i]);
01194 return Gamma_1*(gp->E
01195 #ifdef MHD
01196 - 0.5*(gp->B1c*gp->B1c + gp->B2c*gp->B2c + gp->B3c*gp->B3c)
01197 #endif
01198 - 0.5*(gp->M1*gp->M1 + gp->M2*gp->M2 + gp->M3*gp->M3)/gp->d);
01199 #endif
01200 }
01201
01202
01203
01204
01205
01206 #ifdef ADIABATIC
01207 Real expr_cs2(const GridS *pG, const int i, const int j, const int k)
01208 {
01209 ConsS *gp = &(pG->U[k][j][i]);
01210 return (Gamma*Gamma_1*(gp->E
01211 #ifdef MHD
01212 - 0.5*(gp->B1c*gp->B1c + gp->B2c*gp->B2c + gp->B3c*gp->B3c)
01213 #endif
01214 - 0.5*(gp->M1*gp->M1 + gp->M2*gp->M2 + gp->M3*gp->M3)/gp->d)/gp->d);
01215 }
01216 #endif
01217
01218
01219
01220
01221
01222 #ifdef ADIABATIC
01223 Real expr_S(const GridS *pG, const int i, const int j, const int k)
01224 {
01225 ConsS *gp = &(pG->U[k][j][i]);
01226 Real P = Gamma_1*(gp->E
01227 #ifdef MHD
01228 - 0.5*(gp->B1c*gp->B1c + gp->B2c*gp->B2c + gp->B3c*gp->B3c)
01229 #endif
01230 - 0.5*(gp->M1*gp->M1 + gp->M2*gp->M2 + gp->M3*gp->M3)/gp->d);
01231 return P/pow((double)gp->d, (double)Gamma);
01232 }
01233 #endif
01234
01235
01236
01237
01238
01239 #ifdef SPECIAL_RELATIVITY
01240 Real expr_G(const GridS *pG, const int i, const int j, const int k)
01241 {
01242 PrimS W;
01243 W = Cons_to_Prim(&(pG->U[k][j][i]));
01244 return 1.0/sqrt(1.0 - (SQR(W.V1)+SQR(W.V2)+SQR(W.V3)));
01245 }
01246 #endif
01247
01248
01249
01250
01251 #ifdef PARTICLES
01252 int check_particle_binning(char *out)
01253 {
01254 if (strcmp(out,"dpar")==0)
01255 return 1;
01256 else if (strcmp(out,"M1par")==0)
01257 return 1;
01258 else if (strcmp(out,"M2par")==0)
01259 return 1;
01260 else if (strcmp(out,"M3par")==0)
01261 return 1;
01262 else if (strcmp(out,"V1par")==0)
01263 return 1;
01264 else if (strcmp(out,"V2par")==0)
01265 return 1;
01266 else if (strcmp(out,"V3par")==0)
01267 return 1;
01268 else
01269 return 0;
01270 }
01271 #endif
01272
01273
01274
01275
01276
01277
01278
01279 static ConsFun_t getexpr(const int n, const char *expr)
01280 {
01281 char ename[32];
01282
01283 sprintf(ename,"expr_out%d",n);
01284
01285 if (strcmp(expr,"d")==0)
01286 return expr_d;
01287 else if (strcmp(expr,"M1")==0)
01288 return expr_M1;
01289 else if (strcmp(expr,"M2")==0)
01290 return expr_M2;
01291 else if (strcmp(expr,"M3")==0)
01292 return expr_M3;
01293 #ifndef BAROTROPIC
01294 else if (strcmp(expr,"E")==0)
01295 return expr_E;
01296 #endif
01297 #ifdef MHD
01298 else if (strcmp(expr,"B1c")==0)
01299 return expr_B1c;
01300 else if (strcmp(expr,"B2c")==0)
01301 return expr_B2c;
01302 else if (strcmp(expr,"B3c")==0)
01303 return expr_B3c;
01304 else if (strcmp(expr,"ME")==0)
01305 return expr_ME;
01306 #endif
01307 else if (strcmp(expr,"V1")==0)
01308 return expr_V1;
01309 else if (strcmp(expr,"V2")==0)
01310 return expr_V2;
01311 else if (strcmp(expr,"V3")==0)
01312 return expr_V3;
01313 else if (strcmp(expr,"P")==0)
01314 return expr_P;
01315 #ifdef ADIABATIC
01316 else if (strcmp(expr,"cs2")==0)
01317 return expr_cs2;
01318 #endif
01319 #ifdef ADIABATIC
01320 else if (strcmp(expr,"S")==0)
01321 return expr_S;
01322 #endif
01323 #ifdef SPECIAL_RELATIVITY
01324 else if (strcmp(expr,"G")==0)
01325 return expr_G;
01326 #endif
01327 #ifdef PARTICLES
01328 else if (strcmp(expr,"dpar")==0)
01329 return expr_dpar;
01330 else if (strcmp(expr,"M1par")==0)
01331 return expr_M1par;
01332 else if (strcmp(expr,"M2par")==0)
01333 return expr_M2par;
01334 else if (strcmp(expr,"M3par")==0)
01335 return expr_M3par;
01336 else if (strcmp(expr,"V1par")==0)
01337 return expr_V1par;
01338 else if (strcmp(expr,"V2par")==0)
01339 return expr_V2par;
01340 else if (strcmp(expr,"V3par")==0)
01341 return expr_V3par;
01342 #endif
01343 else {
01344 ath_perr(-1,"Unknown data expression\n");
01345 return NULL;
01346 }
01347 }
01348
01349
01350
01351
01352
01353
01354
01355
01356
01357 static void free_output(OutputS *pOut)
01358 {
01359 if(pOut->out != NULL) free(pOut->out);
01360 if(pOut->out_fmt != NULL) free(pOut->out_fmt);
01361 if(pOut->dat_fmt != NULL) free(pOut->dat_fmt);
01362 if(pOut->id != NULL) free(pOut->id);
01363 return;
01364 }
01365
01366
01367
01368
01369
01370
01371
01372
01373
01374
01375
01376
01377
01378
01379
01380
01381
01382
01383
01384
01385
01386
01387
01388
01389
01390
01391 static void parse_slice(char *block, char *axname, Real *l, Real *u, int *flag)
01392 {
01393 char *expr, *cp;
01394
01395 if (par_exist(block,axname)) {
01396 expr = par_gets(block,axname);
01397 cp = strchr(expr, ':');
01398 if (cp) {
01399 *cp++ = 0;
01400 while (*cp && isspace(*cp)) cp++;
01401 if (*cp)
01402 *u = atof(cp);
01403 cp = expr;
01404 while (*cp && isspace(*cp)) cp++;
01405 if (*cp)
01406 *l = atof(cp);
01407 } else {
01408 *l = *u = atof(expr);
01409 }
01410 if (*l > *u) {
01411 ath_error("[parse_slice]: lower slice limit %d > upper %d in %s\n",
01412 *l,*u,expr);
01413 }
01414 free(expr);
01415 *flag = 1;
01416 }
01417
01418 }
01419
01420
01421
01422
01423
01424
01425
01426
01427 float *getRGB(char *name)
01428 {
01429 int i;
01430
01431 for (i=0; rgb[i].name && rgb[i].rgb; i++) {
01432 if (strcmp(name,rgb[i].name) == 0)
01433 return rgb[i].rgb;
01434 }
01435
01436
01437
01438 ath_perr(-1,"Fatal error: could not find palette=%s, valid names are:\n",
01439 name);
01440 for (i=0; rgb[i].name && rgb[i].rgb; i++)
01441 ath_perr(-1,"%s ",rgb[i].name);
01442 ath_perr(-1,"\n");
01443 exit(EXIT_FAILURE);
01444
01445 return NULL;
01446 }