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 #include <stdio.h>
00070 #include "defs.h"
00071 #include "athena.h"
00072 #include "globals.h"
00073 #include "prototypes.h"
00074
00075
00076 #define NSCAL 14
00077
00078
00079 #define MAX_USR_H_COUNT 30
00080
00081
00082 static char *usr_label[MAX_USR_H_COUNT];
00083
00084
00085 static ConsFun_t phst_fun[MAX_USR_H_COUNT];
00086
00087 static int usr_hst_cnt = 0;
00088
00089
00090
00091
00092
00093
00094 void dump_history(MeshS *pM, OutputS *pOut)
00095 {
00096 GridS *pG;
00097 DomainS *pD;
00098 int i,j,k,is,ie,js,je,ks,ke,nl,nd;
00099 double dVol, scal[NSCAL + NSCALARS + MAX_USR_H_COUNT], d1;
00100 FILE *pfile;
00101 char *fname,*plev=NULL,*pdom=NULL,*pdir=NULL,fmt[80];
00102 char levstr[8],domstr[8],dirstr[20];
00103 int n, total_hst_cnt, mhst, myID_Comm_Domain=1;
00104 #ifdef MPI_PARALLEL
00105 double my_scal[NSCAL + NSCALARS + MAX_USR_H_COUNT];
00106 int ierr;
00107 #endif
00108 #ifdef CYLINDRICAL
00109 Real x1,x2,x3;
00110 #endif
00111 #ifdef SPECIAL_RELATIVITY
00112 PrimS W;
00113 Real g, g2, g_2;
00114 Real bx, by, bz, vB, b2, Bmag2;
00115 #endif
00116
00117
00118 total_hst_cnt = 9 + NSCALARS + usr_hst_cnt;
00119 #ifdef ADIABATIC
00120 total_hst_cnt++;
00121 #endif
00122 #ifdef MHD
00123 total_hst_cnt += 3;
00124 #endif
00125 #ifdef SELF_GRAVITY
00126 total_hst_cnt += 1;
00127 #endif
00128 #ifdef CYLINDRICAL
00129 total_hst_cnt++;
00130 #endif
00131 #ifdef SPECIAL_RELATIVITY
00132 total_hst_cnt = 12 + usr_hst_cnt;
00133 #ifdef MHD
00134 total_hst_cnt += 6;
00135 #endif
00136 #endif
00137
00138
00139
00140 if(pOut->dat_fmt == NULL){
00141 sprintf(fmt," %%14.6e");
00142 }
00143 else{
00144 sprintf(fmt," %s",pOut->dat_fmt);
00145 }
00146
00147
00148
00149 scal[0] = pM->time;
00150 scal[1] = pM->dt;
00151
00152
00153
00154 for (nl=0; nl<(pM->NLevels); nl++){
00155
00156 for (nd=0; nd<(pM->DomainsPerLevel[nl]); nd++){
00157
00158 if (pM->Domain[nl][nd].Grid != NULL){
00159 pG = pM->Domain[nl][nd].Grid;
00160 pD = (DomainS*)&(pM->Domain[nl][nd]);
00161 is = pG->is, ie = pG->ie;
00162 js = pG->js, je = pG->je;
00163 ks = pG->ks, ke = pG->ke;
00164
00165 for (i=2; i<total_hst_cnt; i++) {
00166 scal[i] = 0.0;
00167 }
00168
00169
00170
00171 for (k=ks; k<=ke; k++) {
00172 for (j=js; j<=je; j++) {
00173 for (i=is; i<=ie; i++) {
00174 dVol = 1.0;
00175 if (pG->dx1 > 0.0) dVol *= pG->dx1;
00176 if (pG->dx2 > 0.0) dVol *= pG->dx2;
00177 if (pG->dx3 > 0.0) dVol *= pG->dx3;
00178 #ifndef SPECIAL_RELATIVITY
00179 #ifdef CYLINDRICAL
00180 cc_pos(pG,i,j,k,&x1,&x2,&x3);
00181 dVol *= x1;
00182 #endif
00183
00184 mhst = 2;
00185 scal[mhst] += dVol*pG->U[k][j][i].d;
00186 d1 = 1.0/pG->U[k][j][i].d;
00187 #ifndef BAROTROPIC
00188 mhst++;
00189 scal[mhst] += dVol*pG->U[k][j][i].E;
00190 #endif
00191 mhst++;
00192 scal[mhst] += dVol*pG->U[k][j][i].M1;
00193 mhst++;
00194 scal[mhst] += dVol*pG->U[k][j][i].M2;
00195 mhst++;
00196 scal[mhst] += dVol*pG->U[k][j][i].M3;
00197 mhst++;
00198 scal[mhst] += dVol*0.5*SQR(pG->U[k][j][i].M1)*d1;
00199 mhst++;
00200 scal[mhst] += dVol*0.5*SQR(pG->U[k][j][i].M2)*d1;
00201 mhst++;
00202 scal[mhst] += dVol*0.5*SQR(pG->U[k][j][i].M3)*d1;
00203 #ifdef MHD
00204 mhst++;
00205 scal[mhst] += dVol*0.5*SQR(pG->U[k][j][i].B1c);
00206 mhst++;
00207 scal[mhst] += dVol*0.5*SQR(pG->U[k][j][i].B2c);
00208 mhst++;
00209 scal[mhst] += dVol*0.5*SQR(pG->U[k][j][i].B3c);
00210 #endif
00211 #ifdef SELF_GRAVITY
00212 mhst++;
00213 scal[mhst] += dVol*pG->U[k][j][i].d*pG->Phi[k][j][i];
00214 #endif
00215 #if (NSCALARS > 0)
00216 for(n=0; n<NSCALARS; n++){
00217 mhst++;
00218 scal[mhst] += dVol*pG->U[k][j][i].s[n];
00219 }
00220 #endif
00221
00222 #ifdef CYLINDRICAL
00223 mhst++;
00224 scal[mhst] += dVol*(x1*pG->U[k][j][i].M2);
00225 #endif
00226
00227 #else
00228
00229 W = Cons_to_Prim (&(pG->U[k][j][i]));
00230
00231
00232 g = pG->U[k][j][i].d/W.d;
00233 g2 = SQR(g);
00234 g_2 = 1.0/g2;
00235
00236 mhst = 2;
00237 scal[mhst] += dVol*pG->U[k][j][i].d;
00238 mhst++;
00239 scal[mhst] += dVol*pG->U[k][j][i].E;
00240 mhst++;
00241 scal[mhst] += dVol*pG->U[k][j][i].M1;
00242 mhst++;
00243 scal[mhst] += dVol*pG->U[k][j][i].M2;
00244 mhst++;
00245 scal[mhst] += dVol*pG->U[k][j][i].M3;
00246
00247 mhst++;
00248 scal[mhst] += dVol*SQR(g);
00249 mhst++;
00250 scal[mhst] += dVol*SQR(g*W.V1);
00251 mhst++;
00252 scal[mhst] += dVol*SQR(g*W.V2);
00253 mhst++;
00254 scal[mhst] += dVol*SQR(g*W.V3);
00255
00256 mhst++;
00257 scal[mhst] += dVol*W.P;
00258
00259 #ifdef MHD
00260
00261 vB = W.V1*pG->U[k][j][i].B1c + W.V2*W.B2c + W.V3*W.B3c;
00262 Bmag2 = SQR(pG->U[k][j][i].B1c) + SQR(W.B2c) + SQR(W.B3c);
00263
00264 bx = g*(pG->U[k][j][i].B1c*g_2 + vB*W.V1);
00265 by = g*(W.B2c*g_2 + vB*W.V2);
00266 bz = g*(W.B3c*g_2 + vB*W.V3);
00267
00268 b2 = Bmag2*g_2 + vB*vB;
00269
00270 mhst++;
00271 scal[mhst] += dVol*(g*vB*g*vB);
00272 mhst++;
00273 scal[mhst] += dVol*bx*bx;
00274 mhst++;
00275 scal[mhst] += dVol*by*by;
00276 mhst++;
00277 scal[mhst] += dVol*bz*bz;
00278 mhst++;
00279 scal[mhst] += dVol*b2;
00280 mhst++;
00281 scal[mhst] += dVol*(Bmag2*(1.0 - 0.5*g_2) - SQR(vB) / 2.0);
00282
00283 #endif
00284
00285 #endif
00286
00287
00288 for(n=0; n<usr_hst_cnt; n++){
00289 mhst++;
00290 scal[mhst] += dVol*(*phst_fun[n])(pG, i, j, k);
00291 }
00292 }
00293 }
00294 }
00295
00296
00297
00298 #ifdef MPI_PARALLEL
00299 for(i=2; i<total_hst_cnt; i++){
00300 my_scal[i] = scal[i];
00301 }
00302 ierr = MPI_Reduce(&(my_scal[2]), &(scal[2]), (total_hst_cnt - 2),
00303 MPI_DOUBLE, MPI_SUM, 0, pD->Comm_Domain);
00304 #endif
00305
00306
00307
00308
00309 #ifdef MPI_PARALLEL
00310 ierr = MPI_Comm_rank(pD->Comm_Domain, &myID_Comm_Domain);
00311 #endif
00312 if((myID_Comm_Domain==0) || (myID_Comm_world==0)){
00313
00314
00315
00316 dVol = pD->MaxX[0] - pD->MinX[0];
00317 #ifdef CYLINDRICAL
00318 dVol = (pD->MaxX[0]*pD->MaxX[0]) - (pD->MinX[0]*pD->MinX[0]);
00319 #endif
00320 if (pD->Nx[1] > 1) dVol *= (pD->MaxX[1] - pD->MinX[1]);
00321 if (pD->Nx[2] > 1) dVol *= (pD->MaxX[2] - pD->MinX[2]);
00322 for(i=2; i<total_hst_cnt; i++){
00323 scal[i] /= dVol;
00324 }
00325
00326
00327
00328 #ifdef MPI_PARALLEL
00329 if (nl>0) {
00330 plev = &levstr[0];
00331 sprintf(plev,"lev%d",nl);
00332 pdir = &dirstr[0];
00333 sprintf(pdir,"../id0/lev%d",nl);
00334 }
00335 #else
00336 if (nl>0) {
00337 plev = &levstr[0];
00338 sprintf(plev,"lev%d",nl);
00339 pdir = &dirstr[0];
00340 sprintf(pdir,"lev%d",nl);
00341 }
00342 #endif
00343
00344 if (nd>0) {
00345 pdom = &domstr[0];
00346 sprintf(pdom,"dom%d",nd);
00347 }
00348
00349 fname = ath_fname(pdir,pM->outfilename,plev,pdom,0,0,NULL,"hst");
00350 if(fname == NULL){
00351 ath_perr(-1,"[dump_history]: Unable to create history filename\n");
00352 }
00353 pfile = fopen(fname,"a");
00354 if(pfile == NULL){
00355 ath_perr(-1,"[dump_history]: Unable to open the history file\n");
00356 }
00357
00358
00359
00360 mhst = 0;
00361 if(pOut->num == 0){
00362 fprintf(pfile,
00363 "# Athena history dump for level=%i domain=%i volume=%e\n",nl,nd,dVol);
00364 mhst++;
00365 fprintf(pfile,"# [%i]=time ",mhst);
00366 mhst++;
00367 fprintf(pfile," [%i]=dt ",mhst);
00368 #ifndef SPECIAL_RELATIVITY
00369 mhst++;
00370 fprintf(pfile," [%i]=mass ",mhst);
00371 #ifdef ADIABATIC
00372 mhst++;
00373 fprintf(pfile," [%i]=total E ",mhst);
00374 #endif
00375 mhst++;
00376 fprintf(pfile," [%i]=x1 Mom. ",mhst);
00377 mhst++;
00378 fprintf(pfile," [%i]=x2 Mom. ",mhst);
00379 mhst++;
00380 fprintf(pfile," [%i]=x3 Mom. ",mhst);
00381 mhst++;
00382 fprintf(pfile," [%i]=x1-KE ",mhst);
00383 mhst++;
00384 fprintf(pfile," [%i]=x2-KE ",mhst);
00385 mhst++;
00386 fprintf(pfile," [%i]=x3-KE ",mhst);
00387 #ifdef MHD
00388 mhst++;
00389 fprintf(pfile," [%i]=x1-ME ",mhst);
00390 mhst++;
00391 fprintf(pfile," [%i]=x2-ME ",mhst);
00392 mhst++;
00393 fprintf(pfile," [%i]=x3-ME ",mhst);
00394 #endif
00395 #ifdef SELF_GRAVITY
00396 mhst++;
00397 fprintf(pfile," [%i]=grav PE ",mhst);
00398 #endif
00399 #if (NSCALARS > 0)
00400 for(n=0; n<NSCALARS; n++){
00401 mhst++;
00402 fprintf(pfile," [%i]=scalar %i",mhst,n);
00403 }
00404 #endif
00405
00406 #ifdef CYLINDRICAL
00407 mhst++;
00408 fprintf(pfile," [%i]=Ang.Mom.",mhst);
00409 #endif
00410
00411 #else
00412 mhst++;
00413 fprintf(pfile," [%i]=mass ",mhst);
00414 mhst++;
00415 fprintf(pfile," [%i]=total E ",mhst);
00416 mhst++;
00417 fprintf(pfile," [%i]=x1 Mom. ",mhst);
00418 mhst++;
00419 fprintf(pfile," [%i]=x2 Mom. ",mhst);
00420 mhst++;
00421 fprintf(pfile," [%i]=x3 Mom." ,mhst);
00422 mhst++;
00423 fprintf(pfile," [%i]=Gamma ",mhst);
00424 mhst++;
00425 fprintf(pfile," [%i]=x1-KE ",mhst);
00426 mhst++;
00427 fprintf(pfile," [%i]=x2-KE ",mhst);
00428 mhst++;
00429 fprintf(pfile," [%i]=x3-KE " ,mhst);
00430 mhst++;
00431 fprintf(pfile," [%i]=Press " ,mhst);
00432 #ifdef MHD
00433 mhst++;
00434 fprintf(pfile," [%i]=x0-ME " ,mhst);
00435 mhst++;
00436 fprintf(pfile," [%i]=x1-ME " ,mhst);
00437 mhst++;
00438 fprintf(pfile," [%i]=x2-ME " ,mhst);
00439 mhst++;
00440 fprintf(pfile," [%i]=x3-ME " ,mhst);
00441 mhst++;
00442 fprintf(pfile," [%i]=bsq " ,mhst);
00443 mhst++;
00444 fprintf(pfile," [%i]=T^00_EM" ,mhst);
00445 #endif
00446 #endif
00447
00448 for(n=0; n<usr_hst_cnt; n++){
00449 mhst++;
00450 fprintf(pfile," [%i]=%s",mhst,usr_label[n]);
00451 }
00452 fprintf(pfile,"\n#\n");
00453 }
00454
00455
00456
00457 for (i=0; i<total_hst_cnt; i++) {
00458 fprintf(pfile,fmt,scal[i]);
00459 }
00460 fprintf(pfile,"\n");
00461 fclose(pfile);
00462
00463 }
00464 }
00465 }
00466 }
00467
00468 return;
00469 }
00470
00471
00472
00473
00474
00475 void dump_history_enroll(const ConsFun_t pfun, const char *label){
00476
00477 if(usr_hst_cnt >= MAX_USR_H_COUNT)
00478 ath_error("[dump_history_enroll]: MAX_USR_H_COUNT = %d exceeded\n",
00479 MAX_USR_H_COUNT);
00480
00481
00482 if((usr_label[usr_hst_cnt] = ath_strdup(label)) == NULL)
00483 ath_error("[dump_history_enroll]: Error on sim_strdup(\"%s\")\n",label);
00484
00485
00486 phst_fun[usr_hst_cnt] = pfun;
00487
00488 usr_hst_cnt++;
00489
00490 return;
00491
00492 }
00493
00494 #undef NSCAL
00495 #undef MAX_USR_H_COUNT