Athena: restart.c Source File

Go to the documentation of this file.
00001 #include "copyright.h"
00002 /*============================================================================*/
00003 /*! \file restart.c
00004  *  \brief Functions for writing and reading restart files.
00005  *
00006  * PURPOSE: Functions for writing and reading restart files.  Restart files
00007  *   begin with the input parameter file (athinput.XX) used to start the run in
00008  *   text format, followed by binary format data for selected quantities (time,
00009  *   dt, etc) for each of the Grid structures being updated by this processor.
00010  *   Lastly, any problem-specific data in binary format is included, which must
00011  *   be written by a user function in the problem generator.  Since all of the
00012  *   mesh data (NLevels, number of Domains, size of Grids, etc.) can be
00013  *   recalculated from the athinput file, only the dependent variables (time,
00014  *   dt, nstep, U, B, etc.) are dumped in restart files.
00015  *
00016  *   The athinput file included at the start of the restart file is parsed by
00017  *   main() in Step 2, using the standard par_* functions, with the parameters
00018  *   used by init_mesh() and init_grid() to re-initialize the grid.  Parameter
00019  *   values read from the athinput file contained in the resfile can be
00020  *   superceded by input from the command line, or another input file.
00021  *
00022  * MPI parallel jobs must be restarted on the same number of processors as they
00023  * were run originally.
00024  *
00025  * With SMR, restart files contain ALL levels and domains being updated by each
00026  * processor in one file, written in the default directory for the process.
00027  *
00028  * CONTAINS PUBLIC FUNCTIONS: 
00029  * - restart_grids() - reads nstep,time,dt,ConsS and B from restart file 
00030  * - dump_restart()  - writes a restart file
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 /*! \fn void restart_grids(char *res_file, MeshS *pM)
00046  *  \brief Reads nstep, time, dt, and arrays of ConsS and interface B
00047  *   for each of the Grid structures in the restart file.  
00048  *
00049  *    By the time this
00050  *   function is called (in Step 6 of main()), the Mesh hierarchy has already
00051  *   been re-initialized by init_mesh() and init_grid() in Step 4 of main()
00052  *   using parameters in the athinput file at the start of this restart file,
00053  *   the command line, or from a new input file.
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 /* Open the restart file */
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 /* Skip over the parameter file at the start of the restart file */
00079 
00080   do{
00081     fgets(line,MAXLEN,fp);
00082   }while(strncmp(line,"<par_end>",9) != 0);
00083 
00084 /* read nstep */
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 /* read time */
00092 
00093   fgets(line,MAXLEN,fp);   /* Read the '\n' preceeding the next string */
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 /* read dt */
00100 
00101   fgets(line,MAXLEN,fp);    /* Read the '\n' preceeding the next string */
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 /* Now loop over all Domains containing a Grid on this processor */
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 /* propagate time and dt to all Grids */
00121 
00122       pG->time = pM->time;
00123       pG->dt   = pM->dt;
00124 
00125 /* Read the density */
00126 
00127       fgets(line,MAXLEN,fp); /* Read the '\n' preceeding the next string */
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 /* Read the x1-momentum */
00140 
00141       fgets(line,MAXLEN,fp); /* Read the '\n' preceeding the next string */
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 /* Read the x2-momentum */
00154 
00155       fgets(line,MAXLEN,fp); /* Read the '\n' preceeding the next string */
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 /* Read the x3-momentum */
00168 
00169       fgets(line,MAXLEN,fp); /* Read the '\n' preceeding the next string */
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 /* Read energy density */
00183 
00184       fgets(line,MAXLEN,fp); /* Read the '\n' preceeding the next string */
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 /* if there is more than one cell in each dimension, need to read one more
00199  * face-centered field component than the number of cells.  [ijk]b is
00200  * the number of extra cells to be read  */
00201 
00202       if (ie > is) ib=1;
00203       if (je > js) jb=1;
00204       if (ke > ks) kb=1;
00205 
00206 /* Read the face-centered x1 B-field */
00207 
00208       fgets(line,MAXLEN,fp); /* Read the '\n' preceeding the next string */
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 /* Read the face-centered x2 B-field */
00221 
00222       fgets(line,MAXLEN,fp); /* Read the '\n' preceeding the next string */
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 /* Read the face-centered x3 B-field */
00235 
00236       fgets(line,MAXLEN,fp); /* Read the '\n' preceeding the next string */
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 /* initialize the cell center magnetic fields as either the average of the face
00249  * centered field if there is more than one cell in that dimension, or just
00250  * the face centered field if not  */
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 /* Read any passively advected scalars */
00266 /* Following code only works if NSCALARS < 10 */
00267 
00268       for (n=0; n<NSCALARS; n++) {
00269         fgets(line,MAXLEN,fp); /* Read the '\n' preceeding the next string */
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 /* Read particle properties and the complete particle list */
00286 
00287       fgets(line,MAXLEN,fp); /* Read the '\n' preceeding the next string */
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++) {          /* particle property list */
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);  /* coef to calc Reynolds number */
00303 
00304       for (i=0; i<pG->partypes; i++)
00305         fread(&(pG->grproperty[i].integrator),sizeof(short),1,fp);
00306 
00307 /* Read the x1-positions */
00308 
00309       fgets(line,MAXLEN,fp); /* Read the '\n' preceeding the next string */
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 /* Read the x2-positions */
00318 
00319       fgets(line,MAXLEN,fp); /* Read the '\n' preceeding the next string */
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 /* Read the x3-positions */
00328 
00329       fgets(line,MAXLEN,fp); /* Read the '\n' preceeding the next string */
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 /* Read the v1 velocity */
00338 
00339       fgets(line,MAXLEN,fp); /* Read the '\n' preceeding the next string */
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 /* Read the v2 velocity */
00348 
00349       fgets(line,MAXLEN,fp); /* Read the '\n' preceeding the next string */
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 /* Read the v3 velocity */
00358 
00359       fgets(line,MAXLEN,fp); /* Read the '\n' preceeding the next string */
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 /* Read particle properties */
00368 
00369       fgets(line,MAXLEN,fp); /* Read the '\n' preceeding the next string */
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;        /* grid particle */
00376       }
00377 
00378 /* Read particle my_id */
00379 
00380       fgets(line,MAXLEN,fp); /* Read the '\n' preceeding the next string */
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 /* Read particle init_id */
00390 
00391       fgets(line,MAXLEN,fp); /* Read the '\n' preceeding the next string */
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 /* count the number of particles with different types */
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 /* PARTICLES */
00408 
00409     }
00410   }} /* End loop over all Domains --------------------------------------------*/
00411 
00412 /* Call a user function to read his/her problem-specific data! */
00413 
00414   fgets(line,MAXLEN,fp); /* Read the '\n' preceeding the next string */
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 /*! \fn void dump_restart(MeshS *pM, OutputS *pout)
00427  *  \brief Writes a restart file, including problem-specific data from
00428  *   a user defined function  */
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 /* Allocate memory for buffer */
00452   bufsize = 262144 / sizeof(Real);  /* 256 KB worth of Reals */
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;  /* Right now, we just don't write instead of aborting completely */
00456   }
00457 #ifdef PARTICLES
00458   ibufsize = 262144 / sizeof(int);  /* 256 KB worth of ints */
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;  /* Right now, we just don't write instead of aborting completely */
00462   }
00463   lbufsize = 262144 / sizeof(long);  /* 256 KB worth of longs */
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;  /* Right now, we just don't write instead of aborting completely */
00467   }
00468   sbufsize = 262144 / sizeof(short);  /* 256 KB worth of longs */
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;  /* Right now, we just don't write instead of aborting completely */
00472   }
00473 #endif
00474 
00475 /* Create filename and Open the output file */
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 /* Add the current time & nstep to the parameter file */
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 /* Write the current state of the parameter file */
00493 
00494   par_dump(2,fp);
00495 
00496 /* Write out the current simulation step number */
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 /* Write out the current simulation time */
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 /* Write out the current simulation time step */
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 /* Now loop over all Domains containing a Grid on this processor */
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 /* Write the density */
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 /* Write the x1-momentum */
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 /* Write the x2-momentum */
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 /* Write the x3-momentum */
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 /* Write energy density */
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 /* see comments in restart_grid_block() for use of [ijk]b */
00626 
00627   if (ie > is) ib = 1;
00628   if (je > js) jb = 1;
00629   if (ke > ks) kb = 1;
00630 
00631 /* Write the x1-field */
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 /* Write the x2-field */
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 /* Write the x3-field */
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 /* Write out passively advected scalars */
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 /* Write out the number of particles */
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 /* Write out the particle properties */
00722     
00723 #ifdef FEEDBACK
00724       nprop = 5;
00725 #else
00726       nprop = 4;
00727 #endif
00728       fwrite(&(pG->partypes),sizeof(int),1,fp); /* number of particle types */
00729       for (i=0; i<pG->partypes; i++) {          /* particle property list */
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);  /* coef for Reynolds number */
00747     
00748       for (i=0; i<pG->partypes; i++) {         /* particle integrator type */
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 /* Write x1 */
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 /* Write x2 */
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 /* Write x3 */
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 /* Write v1 */
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 /* Write v2 */
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 /* Write v3 */
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 /* Write properties */
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 /* Write my_id */
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 /* Write init_id */
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 /*PARTICLES*/
00906 
00907     }
00908   }}  /*---------- End loop over all Domains ---------------------------------*/
00909     
00910 /* call a user function to write his/her problem-specific data! */
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 }