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 #include <math.h>
00035 #include <stdio.h>
00036 #include <stdlib.h>
00037 #include <string.h>
00038 #include "defs.h"
00039 #include "athena.h"
00040 #include "globals.h"
00041 #include "prototypes.h"
00042 #include "particles/particle.h"
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056 void restart_grids(char *res_file, MeshS *pM)
00057 {
00058 GridS *pG;
00059 FILE *fp;
00060 char line[MAXLEN];
00061 int i,j,k,is,ie,js,je,ks,ke,nl,nd;
00062 #ifdef MHD
00063 int ib=0,jb=0,kb=0;
00064 #endif
00065 #if (NSCALARS > 0)
00066 int n;
00067 char scalarstr[16];
00068 #endif
00069 #ifdef PARTICLES
00070 long p;
00071 #endif
00072
00073
00074
00075 if((fp = fopen(res_file,"r")) == NULL)
00076 ath_error("[restart_grids]: Error opening the restart file\nIf this is a MPI job, make sure each file from each processor is in the same directory.\n");
00077
00078
00079
00080 do{
00081 fgets(line,MAXLEN,fp);
00082 }while(strncmp(line,"<par_end>",9) != 0);
00083
00084
00085
00086 fgets(line,MAXLEN,fp);
00087 if(strncmp(line,"N_STEP",6) != 0)
00088 ath_error("[restart_grids]: Expected N_STEP, found %s",line);
00089 fread(&(pM->nstep),sizeof(int),1,fp);
00090
00091
00092
00093 fgets(line,MAXLEN,fp);
00094 fgets(line,MAXLEN,fp);
00095 if(strncmp(line,"TIME",4) != 0)
00096 ath_error("[restart_grids]: Expected TIME, found %s",line);
00097 fread(&(pM->time),sizeof(Real),1,fp);
00098
00099
00100
00101 fgets(line,MAXLEN,fp);
00102 fgets(line,MAXLEN,fp);
00103 if(strncmp(line,"TIME_STEP",9) != 0)
00104 ath_error("[restart_grids]: Expected TIME_STEP, found %s",line);
00105 fread(&(pM->dt),sizeof(Real),1,fp);
00106
00107
00108
00109 for (nl=0; nl<=(pM->NLevels)-1; nl++){
00110 for (nd=0; nd<=(pM->DomainsPerLevel[nl])-1; nd++){
00111 if (pM->Domain[nl][nd].Grid != NULL) {
00112 pG=pM->Domain[nl][nd].Grid;
00113 is = pG->is;
00114 ie = pG->ie;
00115 js = pG->js;
00116 je = pG->je;
00117 ks = pG->ks;
00118 ke = pG->ke;
00119
00120
00121
00122 pG->time = pM->time;
00123 pG->dt = pM->dt;
00124
00125
00126
00127 fgets(line,MAXLEN,fp);
00128 fgets(line,MAXLEN,fp);
00129 if(strncmp(line,"DENSITY",7) != 0)
00130 ath_error("[restart_grids]: Expected DENSITY, found %s",line);
00131 for (k=ks; k<=ke; k++) {
00132 for (j=js; j<=je; j++) {
00133 for (i=is; i<=ie; i++) {
00134 fread(&(pG->U[k][j][i].d),sizeof(Real),1,fp);
00135 }
00136 }
00137 }
00138
00139
00140
00141 fgets(line,MAXLEN,fp);
00142 fgets(line,MAXLEN,fp);
00143 if(strncmp(line,"1-MOMENTUM",10) != 0)
00144 ath_error("[restart_grids]: Expected 1-MOMENTUM, found %s",line);
00145 for (k=ks; k<=ke; k++) {
00146 for (j=js; j<=je; j++) {
00147 for (i=is; i<=ie; i++) {
00148 fread(&(pG->U[k][j][i].M1),sizeof(Real),1,fp);
00149 }
00150 }
00151 }
00152
00153
00154
00155 fgets(line,MAXLEN,fp);
00156 fgets(line,MAXLEN,fp);
00157 if(strncmp(line,"2-MOMENTUM",10) != 0)
00158 ath_error("[restart_grids]: Expected 2-MOMENTUM, found %s",line);
00159 for (k=ks; k<=ke; k++) {
00160 for (j=js; j<=je; j++) {
00161 for (i=is; i<=ie; i++) {
00162 fread(&(pG->U[k][j][i].M2),sizeof(Real),1,fp);
00163 }
00164 }
00165 }
00166
00167
00168
00169 fgets(line,MAXLEN,fp);
00170 fgets(line,MAXLEN,fp);
00171 if(strncmp(line,"3-MOMENTUM",10) != 0)
00172 ath_error("[restart_grids]: Expected 3-MOMENTUM, found %s",line);
00173 for (k=ks; k<=ke; k++) {
00174 for (j=js; j<=je; j++) {
00175 for (i=is; i<=ie; i++) {
00176 fread(&(pG->U[k][j][i].M3),sizeof(Real),1,fp);
00177 }
00178 }
00179 }
00180
00181 #ifndef BAROTROPIC
00182
00183
00184 fgets(line,MAXLEN,fp);
00185 fgets(line,MAXLEN,fp);
00186 if(strncmp(line,"ENERGY",6) != 0)
00187 ath_error("[restart_grids]: Expected ENERGY, found %s",line);
00188 for (k=ks; k<=ke; k++) {
00189 for (j=js; j<=je; j++) {
00190 for (i=is; i<=ie; i++) {
00191 fread(&(pG->U[k][j][i].E),sizeof(Real),1,fp);
00192 }
00193 }
00194 }
00195 #endif
00196
00197 #ifdef MHD
00198
00199
00200
00201
00202 if (ie > is) ib=1;
00203 if (je > js) jb=1;
00204 if (ke > ks) kb=1;
00205
00206
00207
00208 fgets(line,MAXLEN,fp);
00209 fgets(line,MAXLEN,fp);
00210 if(strncmp(line,"1-FIELD",7) != 0)
00211 ath_error("[restart_grids]: Expected 1-FIELD, found %s",line);
00212 for (k=ks; k<=ke; k++) {
00213 for (j=js; j<=je; j++) {
00214 for (i=is; i<=ie+ib; i++) {
00215 fread(&(pG->B1i[k][j][i]),sizeof(Real),1,fp);
00216 }
00217 }
00218 }
00219
00220
00221
00222 fgets(line,MAXLEN,fp);
00223 fgets(line,MAXLEN,fp);
00224 if(strncmp(line,"2-FIELD",7) != 0)
00225 ath_error("[restart_grids]: Expected 2-FIELD, found %s",line);
00226 for (k=ks; k<=ke; k++) {
00227 for (j=js; j<=je+jb; j++) {
00228 for (i=is; i<=ie; i++) {
00229 fread(&(pG->B2i[k][j][i]),sizeof(Real),1,fp);
00230 }
00231 }
00232 }
00233
00234
00235
00236 fgets(line,MAXLEN,fp);
00237 fgets(line,MAXLEN,fp);
00238 if(strncmp(line,"3-FIELD",7) != 0)
00239 ath_error("[restart_grids]: Expected 3-FIELD, found %s",line);
00240 for (k=ks; k<=ke+kb; k++) {
00241 for (j=js; j<=je; j++) {
00242 for (i=is; i<=ie; i++) {
00243 fread(&(pG->B3i[k][j][i]),sizeof(Real),1,fp);
00244 }
00245 }
00246 }
00247
00248
00249
00250
00251
00252 for (k=ks; k<=ke; k++) {
00253 for (j=js; j<=je; j++) {
00254 for (i=is; i<=ie; i++) {
00255 pG->U[k][j][i].B1c = pG->B1i[k][j][i];
00256 pG->U[k][j][i].B2c = pG->B2i[k][j][i];
00257 pG->U[k][j][i].B3c = pG->B3i[k][j][i];
00258 if(ib==1) pG->U[k][j][i].B1c=0.5*(pG->B1i[k][j][i] +pG->B1i[k][j][i+1]);
00259 if(jb==1) pG->U[k][j][i].B2c=0.5*(pG->B2i[k][j][i] +pG->B2i[k][j+1][i]);
00260 if(kb==1) pG->U[k][j][i].B3c=0.5*(pG->B3i[k][j][i] +pG->B3i[k+1][j][i]);
00261 }}}
00262 #endif
00263
00264 #if (NSCALARS > 0)
00265
00266
00267
00268 for (n=0; n<NSCALARS; n++) {
00269 fgets(line,MAXLEN,fp);
00270 fgets(line,MAXLEN,fp);
00271 sprintf(scalarstr, "SCALAR %d", n);
00272 if(strncmp(line,scalarstr,8) != 0)
00273 ath_error("[restart_grids]: Expected %s, found %s",scalarstr,line);
00274 for (k=ks; k<=ke; k++) {
00275 for (j=js; j<=je; j++) {
00276 for (i=is; i<=ie; i++) {
00277 fread(&(pG->U[k][j][i].s[n]),sizeof(Real),1,fp);
00278 }
00279 }
00280 }
00281 }
00282 #endif
00283
00284 #ifdef PARTICLES
00285
00286
00287 fgets(line,MAXLEN,fp);
00288 fgets(line,MAXLEN,fp);
00289 if(strncmp(line,"PARTICLE LIST",13) != 0)
00290 ath_error("[restart_grids]: Expected PARTICLE LIST, found %s",line);
00291 fread(&(pG->nparticle),sizeof(long),1,fp);
00292 fread(&(pG->partypes),sizeof(int),1,fp);
00293 for (i=0; i<pG->partypes; i++) {
00294 #ifdef FEEDBACK
00295 fread(&(pG->grproperty[i].m),sizeof(Real),1,fp);
00296 #endif
00297 fread(&(pG->grproperty[i].rad),sizeof(Real),1,fp);
00298 fread(&(pG->grproperty[i].rho),sizeof(Real),1,fp);
00299 fread(&(tstop0[i]),sizeof(Real),1,fp);
00300 fread(&(grrhoa[i]),sizeof(Real),1,fp);
00301 }
00302 fread(&(alamcoeff),sizeof(Real),1,fp);
00303
00304 for (i=0; i<pG->partypes; i++)
00305 fread(&(pG->grproperty[i].integrator),sizeof(short),1,fp);
00306
00307
00308
00309 fgets(line,MAXLEN,fp);
00310 fgets(line,MAXLEN,fp);
00311 if(strncmp(line,"PARTICLE X1",11) != 0)
00312 ath_error("[restart_grids]: Expected PARTICLE X1, found %s",line);
00313 for (p=0; p<pG->nparticle; p++) {
00314 fread(&(pG->particle[p].x1),sizeof(Real),1,fp);
00315 }
00316
00317
00318
00319 fgets(line,MAXLEN,fp);
00320 fgets(line,MAXLEN,fp);
00321 if(strncmp(line,"PARTICLE X2",11) != 0)
00322 ath_error("[restart_grids]: Expected PARTICLE X2, found %s",line);
00323 for (p=0; p<pG->nparticle; p++) {
00324 fread(&(pG->particle[p].x2),sizeof(Real),1,fp);
00325 }
00326
00327
00328
00329 fgets(line,MAXLEN,fp);
00330 fgets(line,MAXLEN,fp);
00331 if(strncmp(line,"PARTICLE X3",11) != 0)
00332 ath_error("[restart_grids]: Expected PARTICLE X3, found %s",line);
00333 for (p=0; p<pG->nparticle; p++) {
00334 fread(&(pG->particle[p].x3),sizeof(Real),1,fp);
00335 }
00336
00337
00338
00339 fgets(line,MAXLEN,fp);
00340 fgets(line,MAXLEN,fp);
00341 if(strncmp(line,"PARTICLE V1",11) != 0)
00342 ath_error("[restart_grids]: Expected PARTICLE V1, found %s",line);
00343 for (p=0; p<pG->nparticle; p++) {
00344 fread(&(pG->particle[p].v1),sizeof(Real),1,fp);
00345 }
00346
00347
00348
00349 fgets(line,MAXLEN,fp);
00350 fgets(line,MAXLEN,fp);
00351 if(strncmp(line,"PARTICLE V2",11) != 0)
00352 ath_error("[restart_grids]: Expected PARTICLE V2, found %s",line);
00353 for (p=0; p<pG->nparticle; p++) {
00354 fread(&(pG->particle[p].v2),sizeof(Real),1,fp);
00355 }
00356
00357
00358
00359 fgets(line,MAXLEN,fp);
00360 fgets(line,MAXLEN,fp);
00361 if(strncmp(line,"PARTICLE V3",11) != 0)
00362 ath_error("[restart_grids]: Expected PARTICLE V3, found %s",line);
00363 for (p=0; p<pG->nparticle; p++) {
00364 fread(&(pG->particle[p].v3),sizeof(Real),1,fp);
00365 }
00366
00367
00368
00369 fgets(line,MAXLEN,fp);
00370 fgets(line,MAXLEN,fp);
00371 if(strncmp(line,"PARTICLE PROPERTY",17) != 0)
00372 ath_error("[restart_grids]: Expected PARTICLE PROPERTY, found %s",line);
00373 for (p=0; p<pG->nparticle; p++) {
00374 fread(&(pG->particle[p].property),sizeof(int),1,fp);
00375 pG->particle[p].pos = 1;
00376 }
00377
00378
00379
00380 fgets(line,MAXLEN,fp);
00381 fgets(line,MAXLEN,fp);
00382 if(strncmp(line,"PARTICLE MY_ID",14) != 0)
00383 ath_error("[restart_grids]: Expected PARTICLE MY_ID, found %s",line);
00384 for (p=0; p<pG->nparticle; p++) {
00385 fread(&(pG->particle[p].my_id),sizeof(long),1,fp);
00386 }
00387
00388 #ifdef MPI_PARALLEL
00389
00390
00391 fgets(line,MAXLEN,fp);
00392 fgets(line,MAXLEN,fp);
00393 if(strncmp(line,"PARTICLE INIT_ID",16) != 0)
00394 ath_error("[restart_grids]: Expected PARTICLE INIT_ID, found %s",line);
00395 for (p=0; p<pG->nparticle; p++) {
00396 fread(&(pG->particle[p].init_id),sizeof(int),1,fp);
00397 }
00398 #endif
00399
00400
00401
00402 for (i=0; i<pG->partypes; i++)
00403 pG->grproperty[i].num = 0;
00404 for (p=0; p<pG->nparticle; p++)
00405 pG->grproperty[pG->particle[p].property].num += 1;
00406
00407 #endif
00408
00409 }
00410 }}
00411
00412
00413
00414 fgets(line,MAXLEN,fp);
00415 fgets(line,MAXLEN,fp);
00416 if(strncmp(line,"USER_DATA",9) != 0)
00417 ath_error("[restart_grids]: Expected USER_DATA, found %s",line);
00418 problem_read_restart(pM, fp);
00419
00420 fclose(fp);
00421
00422 return;
00423 }
00424
00425
00426
00427
00428
00429
00430 void dump_restart(MeshS *pM, OutputS *pout)
00431 {
00432 GridS *pG;
00433 FILE *fp;
00434 char *fname;
00435 int i,j,k,is,ie,js,je,ks,ke,nl,nd;
00436 #ifdef MHD
00437 int ib=0,jb=0,kb=0;
00438 #endif
00439 #if (NSCALARS > 0)
00440 int n;
00441 #endif
00442 #ifdef PARTICLES
00443 int nprop, *ibuf = NULL, ibufsize, lbufsize, sbufsize, nibuf, nlbuf, nsbuf;
00444 long np, p, *lbuf = NULL;
00445 short *sbuf = NULL;
00446 nibuf = 0; nsbuf = 0; nlbuf = 0;
00447 #endif
00448 int bufsize, nbuf = 0;
00449 Real *buf = NULL;
00450
00451
00452 bufsize = 262144 / sizeof(Real);
00453 if ((buf = (Real*)calloc_1d_array(bufsize, sizeof(Real))) == NULL) {
00454 ath_perr(-1,"[dump_restart]: Error allocating memory for buffer\n");
00455 return;
00456 }
00457 #ifdef PARTICLES
00458 ibufsize = 262144 / sizeof(int);
00459 if ((ibuf = (int*)calloc_1d_array(ibufsize, sizeof(int))) == NULL) {
00460 ath_perr(-1,"[dump_restart]: Error allocating memory for buffer\n");
00461 return;
00462 }
00463 lbufsize = 262144 / sizeof(long);
00464 if ((lbuf = (long*)calloc_1d_array(lbufsize, sizeof(long))) == NULL) {
00465 ath_perr(-1,"[dump_restart]: Error allocating memory for buffer\n");
00466 return;
00467 }
00468 sbufsize = 262144 / sizeof(short);
00469 if ((sbuf = (short*)calloc_1d_array(sbufsize, sizeof(short))) == NULL) {
00470 ath_perr(-1,"[dump_restart]: Error allocating memory for buffer\n");
00471 return;
00472 }
00473 #endif
00474
00475
00476
00477 if((fname = ath_fname(NULL,pM->outfilename,NULL,NULL,num_digit,
00478 pout->num,NULL,"rst")) == NULL){
00479 ath_error("[dump_restart]: Error constructing filename\n");
00480 }
00481
00482 if((fp = fopen(fname,"wb")) == NULL){
00483 ath_error("[dump_restart]: Unable to open restart file\n");
00484 return;
00485 }
00486
00487
00488
00489 par_setd("time","time","%e",pM->time,"Current Simulation Time");
00490 par_seti("time","nstep","%d",pM->nstep,"Current Simulation Time Step");
00491
00492
00493
00494 par_dump(2,fp);
00495
00496
00497
00498 fprintf(fp,"N_STEP\n");
00499 if(fwrite(&(pM->nstep),sizeof(int),1,fp) != 1)
00500 ath_error("[dump_restart]: fwrite() error\n");
00501
00502
00503
00504 fprintf(fp,"\nTIME\n");
00505 if(fwrite(&(pM->time),sizeof(Real),1,fp) != 1)
00506 ath_error("[dump_restart]: fwrite() error\n");
00507
00508
00509
00510 fprintf(fp,"\nTIME_STEP\n");
00511 if(fwrite(&(pM->dt),sizeof(Real),1,fp) != 1)
00512 ath_error("[dump_restart]: fwrite() error\n");
00513
00514
00515
00516 for (nl=0; nl<=(pM->NLevels)-1; nl++){
00517 for (nd=0; nd<=(pM->DomainsPerLevel[nl])-1; nd++){
00518 if (pM->Domain[nl][nd].Grid != NULL) {
00519 pG=pM->Domain[nl][nd].Grid;
00520 is = pG->is;
00521 ie = pG->ie;
00522 js = pG->js;
00523 je = pG->je;
00524 ks = pG->ks;
00525 ke = pG->ke;
00526
00527
00528
00529 fprintf(fp,"\nDENSITY\n");
00530 for (k=ks; k<=ke; k++) {
00531 for (j=js; j<=je; j++) {
00532 for (i=is; i<=ie; i++) {
00533 buf[nbuf++] = pG->U[k][j][i].d;
00534 if ((nbuf+1) > bufsize) {
00535 fwrite(buf,sizeof(Real),nbuf,fp);
00536 nbuf = 0;
00537 }
00538 }
00539 }
00540 }
00541 if (nbuf > 0) {
00542 fwrite(buf,sizeof(Real),nbuf,fp);
00543 nbuf = 0;
00544 }
00545
00546
00547
00548 fprintf(fp,"\n1-MOMENTUM\n");
00549 for (k=ks; k<=ke; k++) {
00550 for (j=js; j<=je; j++) {
00551 for (i=is; i<=ie; i++) {
00552 buf[nbuf++] = pG->U[k][j][i].M1;
00553 if ((nbuf+1) > bufsize) {
00554 fwrite(buf,sizeof(Real),nbuf,fp);
00555 nbuf = 0;
00556 }
00557 }
00558 }
00559 }
00560 if (nbuf > 0) {
00561 fwrite(buf,sizeof(Real),nbuf,fp);
00562 nbuf = 0;
00563 }
00564
00565
00566
00567 fprintf(fp,"\n2-MOMENTUM\n");
00568 for (k=ks; k<=ke; k++) {
00569 for (j=js; j<=je; j++) {
00570 for (i=is; i<=ie; i++) {
00571 buf[nbuf++] = pG->U[k][j][i].M2;
00572 if ((nbuf+1) > bufsize) {
00573 fwrite(buf,sizeof(Real),nbuf,fp);
00574 nbuf = 0;
00575 }
00576 }
00577 }
00578 }
00579 if (nbuf > 0) {
00580 fwrite(buf,sizeof(Real),nbuf,fp);
00581 nbuf = 0;
00582 }
00583
00584
00585
00586 fprintf(fp,"\n3-MOMENTUM\n");
00587 for (k=ks; k<=ke; k++) {
00588 for (j=js; j<=je; j++) {
00589 for (i=is; i<=ie; i++) {
00590 buf[nbuf++] = pG->U[k][j][i].M3;
00591 if ((nbuf+1) > bufsize) {
00592 fwrite(buf,sizeof(Real),nbuf,fp);
00593 nbuf = 0;
00594 }
00595 }
00596 }
00597 }
00598 if (nbuf > 0) {
00599 fwrite(buf,sizeof(Real),nbuf,fp);
00600 nbuf = 0;
00601 }
00602
00603 #ifndef BAROTROPIC
00604
00605
00606 fprintf(fp,"\nENERGY\n");
00607 for (k=ks; k<=ke; k++) {
00608 for (j=js; j<=je; j++) {
00609 for (i=is; i<=ie; i++) {
00610 buf[nbuf++] = pG->U[k][j][i].E;
00611 if ((nbuf+1) > bufsize) {
00612 fwrite(buf,sizeof(Real),nbuf,fp);
00613 nbuf = 0;
00614 }
00615 }
00616 }
00617 }
00618 if (nbuf > 0) {
00619 fwrite(buf,sizeof(Real),nbuf,fp);
00620 nbuf = 0;
00621 }
00622 #endif
00623
00624 #ifdef MHD
00625
00626
00627 if (ie > is) ib = 1;
00628 if (je > js) jb = 1;
00629 if (ke > ks) kb = 1;
00630
00631
00632
00633 fprintf(fp,"\n1-FIELD\n");
00634 for (k=ks; k<=ke; k++) {
00635 for (j=js; j<=je; j++) {
00636 for (i=is; i<=ie+ib; i++) {
00637 buf[nbuf++] = pG->B1i[k][j][i];
00638 if ((nbuf+1) > bufsize) {
00639 fwrite(buf,sizeof(Real),nbuf,fp);
00640 nbuf = 0;
00641 }
00642 }
00643 }
00644 }
00645 if (nbuf > 0) {
00646 fwrite(buf,sizeof(Real),nbuf,fp);
00647 nbuf = 0;
00648 }
00649
00650
00651
00652 fprintf(fp,"\n2-FIELD\n");
00653 for (k=ks; k<=ke; k++) {
00654 for (j=js; j<=je+jb; j++) {
00655 for (i=is; i<=ie; i++) {
00656 buf[nbuf++] = pG->B2i[k][j][i];
00657 if ((nbuf+1) > bufsize) {
00658 fwrite(buf,sizeof(Real),nbuf,fp);
00659 nbuf = 0;
00660 }
00661 }
00662 }
00663 }
00664 if (nbuf > 0) {
00665 fwrite(buf,sizeof(Real),nbuf,fp);
00666 nbuf = 0;
00667 }
00668
00669
00670
00671 fprintf(fp,"\n3-FIELD\n");
00672 for (k=ks; k<=ke+kb; k++) {
00673 for (j=js; j<=je; j++) {
00674 for (i=is; i<=ie; i++) {
00675 buf[nbuf++] = pG->B3i[k][j][i];
00676 if ((nbuf+1) > bufsize) {
00677 fwrite(buf,sizeof(Real),nbuf,fp);
00678 nbuf = 0;
00679 }
00680 }
00681 }
00682 }
00683 if (nbuf > 0) {
00684 fwrite(buf,sizeof(Real),nbuf,fp);
00685 nbuf = 0;
00686 }
00687 #endif
00688
00689
00690
00691 #if (NSCALARS > 0)
00692 for (n=0; n<NSCALARS; n++) {
00693 fprintf(fp,"\nSCALAR %d\n", n);
00694 for (k=ks; k<=ke; k++) {
00695 for (j=js; j<=je; j++) {
00696 for (i=is; i<=ie; i++) {
00697 buf[nbuf++] = pG->U[k][j][i].s[n];
00698 if ((nbuf+1) > bufsize) {
00699 fwrite(buf,sizeof(Real),nbuf,fp);
00700 nbuf = 0;
00701 }
00702 }
00703 }
00704 }
00705 if (nbuf > 0) {
00706 fwrite(buf,sizeof(Real),nbuf,fp);
00707 nbuf = 0;
00708 }
00709 }
00710 #endif
00711
00712 #ifdef PARTICLES
00713
00714
00715 fprintf(fp,"\nPARTICLE LIST\n");
00716 np = 0;
00717 for (p=0; p<pG->nparticle; p++)
00718 if (pG->particle[p].pos == 1) np += 1;
00719 fwrite(&(np),sizeof(long),1,fp);
00720
00721
00722
00723 #ifdef FEEDBACK
00724 nprop = 5;
00725 #else
00726 nprop = 4;
00727 #endif
00728 fwrite(&(pG->partypes),sizeof(int),1,fp);
00729 for (i=0; i<pG->partypes; i++) {
00730 #ifdef FEEDBACK
00731 buf[nbuf++] = pG->grproperty[i].m;
00732 #endif
00733 buf[nbuf++] = pG->grproperty[i].rad;
00734 buf[nbuf++] = pG->grproperty[i].rho;
00735 buf[nbuf++] = tstop0[i];
00736 buf[nbuf++] = grrhoa[i];
00737 if ((nbuf+nprop) > bufsize) {
00738 fwrite(buf,sizeof(Real),nbuf,fp);
00739 nbuf = 0;
00740 }
00741 }
00742 if (nbuf > 0) {
00743 fwrite(buf,sizeof(Real),nbuf,fp);
00744 nbuf = 0;
00745 }
00746 fwrite(&(alamcoeff),sizeof(Real),1,fp);
00747
00748 for (i=0; i<pG->partypes; i++) {
00749 sbuf[nsbuf++] = pG->grproperty[i].integrator;
00750 if ((nsbuf+1) > sbufsize) {
00751 fwrite(sbuf,sizeof(short),nsbuf,fp);
00752 nsbuf = 0;
00753 }
00754 }
00755 if (nsbuf > 0) {
00756 fwrite(sbuf,sizeof(short),nsbuf,fp);
00757 nsbuf = 0;
00758 }
00759
00760
00761
00762 fprintf(fp,"\nPARTICLE X1\n");
00763 for (p=0;p<pG->nparticle;p++)
00764 if (pG->particle[p].pos == 1){
00765 buf[nbuf++] = pG->particle[p].x1;
00766 if ((nbuf+1) > bufsize) {
00767 fwrite(buf,sizeof(Real),nbuf,fp);
00768 nbuf = 0;
00769 }
00770 }
00771 if (nbuf > 0) {
00772 fwrite(buf,sizeof(Real),nbuf,fp);
00773 nbuf = 0;
00774 }
00775
00776
00777
00778 fprintf(fp,"\nPARTICLE X2\n");
00779 for (p=0;p<pG->nparticle;p++)
00780 if (pG->particle[p].pos == 1){
00781 buf[nbuf++] = pG->particle[p].x2;
00782 if ((nbuf+1) > bufsize) {
00783 fwrite(buf,sizeof(Real),nbuf,fp);
00784 nbuf = 0;
00785 }
00786 }
00787 if (nbuf > 0) {
00788 fwrite(buf,sizeof(Real),nbuf,fp);
00789 nbuf = 0;
00790 }
00791
00792
00793
00794 fprintf(fp,"\nPARTICLE X3\n");
00795 for (p=0;p<pG->nparticle;p++)
00796 if (pG->particle[p].pos == 1){
00797 buf[nbuf++] = pG->particle[p].x3;
00798 if ((nbuf+1) > bufsize) {
00799 fwrite(buf,sizeof(Real),nbuf,fp);
00800 nbuf = 0;
00801 }
00802 }
00803 if (nbuf > 0) {
00804 fwrite(buf,sizeof(Real),nbuf,fp);
00805 nbuf = 0;
00806 }
00807
00808
00809
00810 fprintf(fp,"\nPARTICLE V1\n");
00811 for (p=0;p<pG->nparticle;p++)
00812 if (pG->particle[p].pos == 1){
00813 buf[nbuf++] = pG->particle[p].v1;
00814 if ((nbuf+1) > bufsize) {
00815 fwrite(buf,sizeof(Real),nbuf,fp);
00816 nbuf = 0;
00817 }
00818 }
00819 if (nbuf > 0) {
00820 fwrite(buf,sizeof(Real),nbuf,fp);
00821 nbuf = 0;
00822 }
00823
00824
00825
00826 fprintf(fp,"\nPARTICLE V2\n");
00827 for (p=0;p<pG->nparticle;p++)
00828 if (pG->particle[p].pos == 1){
00829 buf[nbuf++] = pG->particle[p].v2;
00830 if ((nbuf+1) > bufsize) {
00831 fwrite(buf,sizeof(Real),nbuf,fp);
00832 nbuf = 0;
00833 }
00834 }
00835 if (nbuf > 0) {
00836 fwrite(buf,sizeof(Real),nbuf,fp);
00837 nbuf = 0;
00838 }
00839
00840
00841
00842 fprintf(fp,"\nPARTICLE V3\n");
00843 for (p=0;p<pG->nparticle;p++)
00844 if (pG->particle[p].pos == 1){
00845 buf[nbuf++] = pG->particle[p].v3;
00846 if ((nbuf+1) > bufsize) {
00847 fwrite(buf,sizeof(Real),nbuf,fp);
00848 nbuf = 0;
00849 }
00850 }
00851 if (nbuf > 0) {
00852 fwrite(buf,sizeof(Real),nbuf,fp);
00853 nbuf = 0;
00854 }
00855
00856
00857
00858 fprintf(fp,"\nPARTICLE PROPERTY\n");
00859 for (p=0;p<pG->nparticle;p++)
00860 if (pG->particle[p].pos == 1){
00861 ibuf[nibuf++] = pG->particle[p].property;
00862 if ((nibuf+1) > ibufsize) {
00863 fwrite(ibuf,sizeof(int),nibuf,fp);
00864 nibuf = 0;
00865 }
00866 }
00867 if (nibuf > 0) {
00868 fwrite(ibuf,sizeof(int),nibuf,fp);
00869 nibuf = 0;
00870 }
00871
00872
00873
00874 fprintf(fp,"\nPARTICLE MY_ID\n");
00875 for (p=0;p<pG->nparticle;p++)
00876 if (pG->particle[p].pos == 1){
00877 lbuf[nlbuf++] = pG->particle[p].my_id;
00878 if ((nlbuf+1) > lbufsize) {
00879 fwrite(lbuf,sizeof(long),nlbuf,fp);
00880 nlbuf = 0;
00881 }
00882 }
00883 if (nlbuf > 0) {
00884 fwrite(lbuf,sizeof(long),nlbuf,fp);
00885 nlbuf = 0;
00886 }
00887
00888 #ifdef MPI_PARALLEL
00889
00890
00891 fprintf(fp,"\nPARTICLE INIT_ID\n");
00892 for (p=0;p<pG->nparticle;p++)
00893 if (pG->particle[p].pos == 1){
00894 ibuf[nibuf++] = pG->particle[p].init_id;
00895 if ((nibuf+1) > ibufsize) {
00896 fwrite(ibuf,sizeof(int),nibuf,fp);
00897 nibuf = 0;
00898 }
00899 }
00900 if (nibuf > 0) {
00901 fwrite(ibuf,sizeof(int),nibuf,fp);
00902 nibuf = 0;
00903 }
00904 #endif
00905 #endif
00906
00907 }
00908 }}
00909
00910
00911
00912 fprintf(fp,"\nUSER_DATA\n");
00913 problem_write_restart(pM, fp);
00914
00915 fclose(fp);
00916
00917 free_1d_array(buf);
00918
00919 return;
00920 }