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 #include <stdio.h>
00051 #include <math.h>
00052 #include "../defs.h"
00053 #include "../athena.h"
00054 #include "../prototypes.h"
00055 #include "prototypes.h"
00056 #include "particle.h"
00057
00058 #ifdef PARTICLES
00059
00060
00061 #define NSCAL 11
00062 #define NARAY 12
00063
00064
00065 #define MAX_USR_SCAL 10
00066 #define MAX_USR_ARAY 10
00067
00068
00069 static char *usr_label_scal[MAX_USR_SCAL];
00070 static char *usr_label_aray[MAX_USR_ARAY];
00071
00072
00073 static Parfun_t phst_scalfun[MAX_USR_SCAL];
00074 static Parfun_t phst_arayfun[MAX_USR_ARAY];
00075
00076 static int usr_scal_cnt = 0;
00077 static int usr_aray_cnt = 0;
00078
00079 extern Real expr_dpar(const Grid *pG, const int i, const int j, const int k);
00080
00081
00082
00083
00084
00085
00086
00087 void dump_particle_history(Grid *pGrid, Domain *pD, Output *pOut)
00088 {
00089 FILE *fid;
00090 int i,j,k,n,prp,mhst;
00091 long p,vol_rat,*npar;
00092 int tot_scal_cnt,tot_aray_cnt;
00093 Real scal[NSCAL+MAX_USR_SCAL],**array,rho,dvol;
00094 char fmt[20];
00095 Grain *gr;
00096 #ifdef MPI_PARALLEL
00097 Real my_scal[NSCAL+MAX_USR_SCAL],*sendbuf,*recvbuf;
00098 int err;
00099 long *my_npar;
00100
00101 my_npar = (long*)calloc_1d_array(pGrid->partypes, sizeof(long));
00102 sendbuf = (Real*)calloc_1d_array( (NARAY+MAX_USR_ARAY)*pGrid->partypes,
00103 sizeof(Real));
00104 recvbuf = (Real*)calloc_1d_array( (NARAY+MAX_USR_ARAY)*pGrid->partypes,
00105 sizeof(Real));
00106 #endif
00107
00108 npar = (long*)calloc_1d_array(pGrid->partypes, sizeof(long));
00109 array = (Real**)calloc_2d_array(NARAY+MAX_USR_ARAY, pGrid->partypes,
00110 sizeof(Real));
00111
00112 tot_scal_cnt = 11 + usr_scal_cnt;
00113 tot_aray_cnt = 12 + usr_aray_cnt;
00114
00115 particle_to_grid(pGrid, pD, property_all);
00116
00117
00118 if(pOut->dat_fmt == NULL){
00119 sprintf(fmt," %%13.5e");
00120 }
00121 else{
00122 sprintf(fmt," %s",pOut->dat_fmt);
00123 }
00124
00125 for (i=1; i<tot_scal_cnt; i++) {
00126 scal[i] = 0.0;
00127 }
00128
00129
00130 scal[0] = pGrid->time;
00131
00132
00133 scal[1] = 0.0;
00134 scal[2] = 0.0;
00135 scal[3] = 0.0;
00136 for (k=pGrid->ks; k<=pGrid->ke; k++)
00137 for (j=pGrid->js; j<=pGrid->je; j++)
00138 for (i=pGrid->is; i<=pGrid->ie; i++)
00139 {
00140 scal[1] = MAX(scal[1],expr_dpar(pGrid,i,j,k));
00141 #ifdef FEEDBACK
00142 scal[2] = MAX(scal[2],pGrid->Coup[k][j][i].FBstiff);
00143 scal[3]+= pGrid->Coup[k][j][i].Eloss;
00144 #endif
00145 }
00146
00147
00148 for(p=0; p<pGrid->nparticle; p++) {
00149 gr = &(pGrid->particle[p]);
00150 if (gr->pos == 1)
00151 {
00152 #ifdef FEEDBACK
00153 rho = pGrid->grproperty[gr->property].m;
00154 #else
00155 rho = 1.0;
00156 #endif
00157 mhst = 4;
00158 scal[mhst] += rho;
00159 mhst++;
00160 scal[mhst] += rho*gr->v1;
00161 mhst++;
00162 scal[mhst] += rho*gr->v2;
00163 mhst++;
00164 scal[mhst] += rho*gr->v3;
00165 mhst++;
00166 scal[mhst] += 0.5*rho*SQR(gr->v1);
00167 mhst++;
00168 scal[mhst] += 0.5*rho*SQR(gr->v2);
00169 mhst++;
00170 scal[mhst] += 0.5*rho*SQR(gr->v3);
00171
00172
00173 for(n=0; n<usr_scal_cnt; n++){
00174 mhst++;
00175 scal[mhst] += (*phst_scalfun[n])(pGrid, gr);
00176 }
00177 }
00178 }
00179
00180 #ifdef MPI_PARALLEL
00181 for (i=1; i<tot_scal_cnt; i++)
00182 my_scal[i] = scal[i];
00183
00184 err = MPI_Reduce(&(my_scal[1]),&(scal[1]),2,
00185 MPI_DOUBLE,MPI_MAX,0,MPI_COMM_WORLD);
00186 err = MPI_Reduce(&(my_scal[3]),&(scal[3]),8+usr_scal_cnt,
00187 MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);
00188 #endif
00189
00190
00191 vol_rat = (pD->ide - pD->ids + 1)*(pD->jde - pD->jds + 1)*
00192 (pD->kde - pD->kds + 1);
00193
00194 if(pGrid->my_id == 0){
00195
00196 dvol = 1.0/(double)vol_rat;
00197
00198 for (i=3; i<tot_scal_cnt; i++)
00199 scal[i] *= dvol;
00200 }
00201
00202
00203
00204 mhst = 0;
00205
00206 for (i=0; i<pGrid->partypes; i++)
00207 {
00208 npar[i] = 0;
00209
00210 for (j=0; j<tot_aray_cnt; j++)
00211 array[j][i] = 0.0;
00212 }
00213
00214
00215 for (p=0; p<pGrid->nparticle; p++)
00216 {
00217 gr = &(pGrid->particle[p]);
00218 if (gr->pos == 1)
00219 {
00220 prp = gr->property;
00221 npar[prp] += 1;
00222 array[0][prp] += gr->x1;
00223 array[1][prp] += gr->x2;
00224 array[2][prp] += gr->x3;
00225 array[3][prp] += gr->v1;
00226 array[4][prp] += gr->v2;
00227 array[5][prp] += gr->v3;
00228 }
00229 }
00230
00231 #ifdef MPI_PARALLEL
00232 for (i=0; i<pGrid->partypes; i++)
00233 my_npar[i] = npar[i];
00234
00235 for (j=0; j<6; j++)
00236 {
00237 n = j*pGrid->partypes;
00238 for (i=0; i<pGrid->partypes; i++)
00239 sendbuf[n+i] = array[j][i];
00240 }
00241
00242 err = MPI_Allreduce(my_npar, npar, pGrid->partypes,
00243 MPI_LONG, MPI_SUM,MPI_COMM_WORLD);
00244 err = MPI_Allreduce(sendbuf,recvbuf,6*pGrid->partypes,
00245 MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
00246
00247 for (j=0; j<6; j++)
00248 {
00249 n = j*pGrid->partypes;
00250 for (i=0; i<pGrid->partypes; i++)
00251 array[j][i] = recvbuf[n+i];
00252 }
00253 #endif
00254
00255 for (j=0; j<6; j++) {
00256 for (i=0; i<pGrid->partypes; i++) {
00257 array[j][i] = array[j][i]/MAX(1,npar[i]);
00258 }}
00259
00260
00261 for (p=0; p<pGrid->nparticle; p++)
00262 {
00263 gr = &(pGrid->particle[p]);
00264 if (gr->pos == 1)
00265 {
00266 prp = gr->property;
00267 array[6][prp] += SQR(gr->x1-array[0][prp]);
00268 array[7][prp] += SQR(gr->x2-array[1][prp]);
00269 array[8][prp] += SQR(gr->x3-array[2][prp]);
00270 array[9][prp] += SQR(gr->v1-array[3][prp]);
00271 array[10][prp]+= SQR(gr->v2-array[4][prp]);
00272 array[11][prp]+= SQR(gr->v3-array[5][prp]);
00273 mhst = 11;
00274 for (n=0; n<usr_aray_cnt; n++){
00275 mhst++;
00276 array[mhst][prp] += (*phst_arayfun[n])(pGrid, gr);
00277 }
00278 }
00279 }
00280
00281 #ifdef MPI_PARALLEL
00282 for (j=6; j<tot_aray_cnt; j++)
00283 {
00284 n = (j-6)*pGrid->partypes;
00285 for (i=0; i<pGrid->partypes; i++)
00286 sendbuf[n+i] = array[j][i];
00287 }
00288
00289 err = MPI_Reduce(sendbuf,recvbuf,(tot_aray_cnt-6)*pGrid->partypes,
00290 MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);
00291
00292 if (pGrid->my_id == 0)
00293 {
00294 for (j=6; j<tot_aray_cnt; j++)
00295 {
00296 n = (j-6)*pGrid->partypes;
00297 for (i=0; i<pGrid->partypes; i++)
00298 array[j][i] = recvbuf[n+i];
00299 }
00300
00301 }
00302 #endif
00303
00304 if (pGrid->my_id == 0) {
00305
00306 for (i=0; i<pGrid->partypes; i++) {
00307 for (j=6; j<12; j++)
00308 array[j][i] = sqrt(array[j][i]/MAX(1,npar[i]-1));
00309
00310 for (j=12; j<tot_aray_cnt; j++)
00311 array[j][i] = array[j][i]/MAX(1,npar[i]);
00312 }
00313
00314 }
00315
00316
00317
00318 if (pGrid->my_id == 0) {
00319
00320 #ifdef MPI_PARALLEL
00321 fid = ath_fopen("../",pGrid->outfilename,0,0,NULL,"phst","a");
00322 #else
00323 fid = ath_fopen(NULL,pGrid->outfilename,0,0,NULL,"phst","a");
00324 #endif
00325 if(fid == NULL){
00326 ath_perr(-1,"[dump_particle_history]: Unable to open the history file\n");
00327 return;
00328 }
00329
00330
00331 if(pOut->num == 0){
00332
00333 fprintf(fid,"#Global scalars:\n");
00334 mhst = 1;
00335 fprintf(fid,"# [%i]=time ",mhst);
00336 mhst++;
00337 fprintf(fid," [%i]=d_max ",mhst);
00338 mhst++;
00339 fprintf(fid," [%i]=stiffmax",mhst);
00340 mhst++;
00341 fprintf(fid," [%i]=Edot ",mhst);
00342 mhst++;
00343 fprintf(fid," [%i]=mass ",mhst);
00344 mhst++;
00345 fprintf(fid," [%i]=x1 Mom. ",mhst);
00346 mhst++;
00347 fprintf(fid," [%i]=x2 Mom. ",mhst);
00348 mhst++;
00349 fprintf(fid," [%i]=x3 Mom. ",mhst);
00350 mhst++;
00351 fprintf(fid," [%i]=x1-KE ",mhst);
00352 mhst++;
00353 fprintf(fid," [%i]=x2-KE ",mhst);
00354 mhst++;
00355 fprintf(fid," [%i]=x3-KE ",mhst);
00356 for(n=0; n<usr_scal_cnt; n++){
00357 mhst++;
00358 fprintf(fid," [%i]=%s",mhst,usr_label_scal[n]);
00359 }
00360
00361 fprintf(fid,"\n#Particle type dependent scalars:\n");
00362 mhst = 1;
00363 fprintf(fid,"# [%i]=x1 avg ",mhst);
00364 mhst++;
00365 fprintf(fid," [%i]=x2 avg ",mhst);
00366 mhst++;
00367 fprintf(fid," [%i]=x3 avg ",mhst);
00368 mhst++;
00369 fprintf(fid," [%i]=v1 avg ",mhst);
00370 mhst++;
00371 fprintf(fid," [%i]=v2 avg ",mhst);
00372 mhst++;
00373 fprintf(fid," [%i]=v3 avg ",mhst);
00374 mhst++;
00375 fprintf(fid," [%i]=x1 var ",mhst);
00376 mhst++;
00377 fprintf(fid," [%i]=x2 var ",mhst);
00378 mhst++;
00379 fprintf(fid," [%i]=x3 var ",mhst);
00380 mhst++;
00381 fprintf(fid," [%i]=v1 var ",mhst);
00382 mhst++;
00383 fprintf(fid," [%i]=v2 var ",mhst);
00384 mhst++;
00385 fprintf(fid," [%i]=v3 var ",mhst);
00386 for(n=0; n<usr_aray_cnt; n++){
00387 mhst++;
00388 fprintf(fid," [%i]=%s",mhst,usr_label_aray[n]);
00389 }
00390 fprintf(fid,"\n#\n");
00391 }
00392
00393
00394 for (i=0; i<tot_scal_cnt; i++) {
00395 fprintf(fid,fmt,scal[i]);
00396 }
00397 fprintf(fid,"\n");
00398
00399 for (j=0; j<pGrid->partypes; j++) {
00400
00401 for (i=0; i<tot_aray_cnt; i++)
00402 fprintf(fid,fmt,array[i][j]);
00403 fprintf(fid,"\n");
00404 }
00405 fprintf(fid,"\n");
00406
00407 fclose(fid);
00408 }
00409
00410 free(npar);
00411 free_2d_array(array);
00412 #ifdef MPI_PARALLEL
00413 free(my_npar);
00414 free(sendbuf);
00415 free(recvbuf);
00416 #endif
00417
00418 return;
00419 }
00420
00421
00422
00423
00424
00425
00426 void dump_parhistory_enroll(const Parfun_t pfun, const char *label,
00427 const int set)
00428 {
00429 if (set == 0)
00430 {
00431 if(usr_scal_cnt >= MAX_USR_SCAL)
00432 ath_error("[par_history_enroll]: MAX_USR_SCAL = %d exceeded\n",
00433 MAX_USR_SCAL);
00434
00435
00436 if((usr_label_scal[usr_scal_cnt] = ath_strdup(label)) == NULL)
00437 ath_error("[par_history_enroll]: Error on sim_strdup(\"%s\")\n",label);
00438
00439
00440 phst_scalfun[usr_scal_cnt] = pfun;
00441 usr_scal_cnt++;
00442 }
00443 else
00444 {
00445 if(usr_aray_cnt >= MAX_USR_ARAY)
00446 ath_error("[par_history_enroll]: MAX_USR_ARAY = %d exceeded\n",
00447 MAX_USR_ARAY);
00448
00449
00450 if((usr_label_aray[usr_aray_cnt] = ath_strdup(label)) == NULL)
00451 ath_error("[par_history_enroll]: Error on sim_strdup(\"%s\")\n",label);
00452
00453
00454 phst_arayfun[usr_aray_cnt] = pfun;
00455 usr_aray_cnt++;
00456 }
00457
00458 return;
00459 }
00460
00461 #undef NSCAL
00462 #undef NARAY
00463
00464 #undef MAX_USR_SCAL
00465 #undef MAX_USR_ARAY
00466
00467 #endif
00468