xref: /phasta/phSolver/common/new_interface.c (revision 32eed1412d9a3e748fe578866226c6021c4817b2)
1 /* This file provides interface functions for 'partial ' random
2    access into the PHASTA input files
3 
4    Anil Karanam March 2001 */
5 
6 #include <stdio.h>
7 #include <string.h>
8 #include <ctype.h>
9 #include <stdlib.h>
10 #include <time.h>
11 #include <math.h>
12 #include "mpi.h"
13 #include "phastaIO.h"
14 #include "rdtsc.h"
15 #include <FCMangle.h>
16 #include "new_interface.h"
17 
18 //MR CHANGE
19 #include "common_c.h"
20 //MR CHANGE END
21 
22 #ifdef intel
23 #include <winsock2.h>
24 #else
25 #include <unistd.h>
26 #include <strings.h>
27 #endif
28 
29 //extern double cpu_speed = 2600000000.0; //for Jaguar XT5
30 //extern double cpu_speed = 2100000000.0;  //for Jaguar xt4
31 //extern double cpu_speed =  850000000.0;  //for Intrepid
32 
33 void igetMinMaxAvg(int *ivalue, double *stats, int *statRanks) {
34   int isThisRank;
35 
36   double *value = (double*)malloc(sizeof(double));
37   *value = 1.0*(*ivalue);
38 
39   rgetMinMaxAvg(value,stats,statRanks);
40 
41   /* MPI_Allreduce(value,&stats[0],1,MPI_DOUBLE,MPI_MIN,MPI_COMM_WORLD);
42   isThisRank=workfc.numpe+1;
43   if(*value==stats[0])
44     isThisRank=workfc.myrank;
45   MPI_Allreduce(&isThisRank,&statRanks[0],1,MPI_INT,MPI_MIN,MPI_COMM_WORLD);
46 
47   MPI_Allreduce(value,&stats[1],1,MPI_DOUBLE,MPI_MAX,MPI_COMM_WORLD);
48   isThisRank=workfc.numpe+1;
49   if(*value==stats[1])
50     isThisRank=workfc.myrank;
51   MPI_Allreduce(&isThisRank,&statRanks[1],1,MPI_INT,MPI_MIN,MPI_COMM_WORLD);
52 
53   MPI_Allreduce(value,&stats[2],1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
54   stats[2] /= workfc.numpe; */
55 
56   free(value);
57 }
58 
59 void rgetMinMaxAvg(double *value, double *stats, int *statRanks) {
60   int isThisRank;
61 
62   MPI_Allreduce(value,&stats[0],1,MPI_DOUBLE,MPI_MIN,MPI_COMM_WORLD);
63   isThisRank=workfc.numpe+1;
64   if(*value==stats[0])
65     isThisRank=workfc.myrank;
66   MPI_Allreduce(&isThisRank,&statRanks[0],1,MPI_INT,MPI_MIN,MPI_COMM_WORLD);
67 
68   MPI_Allreduce(value,&stats[1],1,MPI_DOUBLE,MPI_MAX,MPI_COMM_WORLD);
69   isThisRank=workfc.numpe+1;
70   if(*value==stats[1])
71     isThisRank=workfc.myrank;
72   MPI_Allreduce(&isThisRank,&statRanks[1],1,MPI_INT,MPI_MIN,MPI_COMM_WORLD);
73 
74   MPI_Allreduce(value,&stats[2],1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
75   stats[2] /= workfc.numpe;
76 
77   double sqValue = (*value)*(*value), sqValueAvg = 0.;
78   MPI_Allreduce(&sqValue,&sqValueAvg,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
79   sqValueAvg /= workfc.numpe;
80   // stats[3] = sqValueAvg;
81 
82   stats[3] = sqrt(sqValueAvg-stats[2]*stats[2]);
83 }
84 
85 void print_mesh_stats(void) {
86   int statRanks[2];
87   double iStats[4], rStats[4];
88 
89   igetMinMaxAvg(&conpar.nshg,iStats,statRanks);
90   if(workfc.myrank==workfc.master)
91     printf("nshg    : min [%d,%d], max[%d,%d] and avg[.,%d] (rms=%d)\n",statRanks[0],(int)iStats[0],statRanks[1],(int)iStats[1],(int)iStats[2],(int)iStats[3]);
92   igetMinMaxAvg(&conpar.numel,iStats,statRanks);
93   if(workfc.myrank==workfc.master)
94     printf("numel   : min [%d,%d], max[%d,%d] and avg[.,%d] (rms=%d)\n",statRanks[0],(int)iStats[0],statRanks[1],(int)iStats[1],(int)iStats[2],(int)iStats[3]);
95   igetMinMaxAvg(&conpar.numelb,iStats,statRanks);
96   if(workfc.myrank==workfc.master)
97     printf("numelb  : min [%d,%d], max[%d,%d] and avg[.,%d] (rms=%d)\n",statRanks[0],(int)iStats[0],statRanks[1],(int)iStats[1],(int)iStats[2],(int)iStats[3]);
98   igetMinMaxAvg(&conpar.nnz_tot,iStats,statRanks);
99   if(workfc.myrank==workfc.master) {
100     printf("nnz_tot : min [%d,%d], max[%d,%d] and avg[.,%d] (rms=%d)\n",statRanks[0],(int)iStats[0],statRanks[1],(int)iStats[1],(int)iStats[2],(int)iStats[3]);
101     printf("\n");
102   }
103 }
104 
105 void print_mpi_stats(void) {
106   int statRanks[2];
107   double iStats[4], rStats[4];
108 
109 // NS equations
110   igetMinMaxAvg(&mpistats.iISend,iStats,statRanks);
111   if(workfc.myrank==workfc.master)
112     printf("iISend : min [%d,%d], max[%d,%d] and avg[.,%d] (rms=%d)\n",statRanks[0],(int)iStats[0],statRanks[1],(int)iStats[1],(int)iStats[2],(int)iStats[3]);
113   igetMinMaxAvg(&mpistats.iIRecv,iStats,statRanks);
114   if(workfc.myrank==workfc.master)
115     printf("iIRecv : min [%d,%d], max[%d,%d] and avg[.,%d] (rms=%d)\n",statRanks[0],(int)iStats[0],statRanks[1],(int)iStats[1],(int)iStats[2],(int)iStats[3]);
116   igetMinMaxAvg(&mpistats.iWaitAll,iStats,statRanks);
117   if(workfc.myrank==workfc.master)
118     printf("iWtAll : min [%d,%d], max[%d,%d] and avg[.,%d] (rms=%d)\n",statRanks[0],(int)iStats[0],statRanks[1],(int)iStats[1],(int)iStats[2],(int)iStats[3]);
119   igetMinMaxAvg(&mpistats.iAllR,iStats,statRanks);
120   if(workfc.myrank==workfc.master)
121     printf("iAllR  : min [%d,%d], max[%d,%d] and avg[.,%d] (rms=%d)\n",statRanks[0],(int)iStats[0],statRanks[1],(int)iStats[1],(int)iStats[2],(int)iStats[3]);
122 
123   rgetMinMaxAvg(&mpistats.rISend,rStats,statRanks);
124   if(workfc.myrank==workfc.master)
125     printf("rISend : min [%d,%2.5f], max[%d,%2.5f] and avg[.,%2.5f] (rms=%2.5f)\n",statRanks[0],rStats[0],statRanks[1],rStats[1],rStats[2],rStats[3]);
126   rgetMinMaxAvg(&mpistats.rIRecv,rStats,statRanks);
127   if(workfc.myrank==workfc.master)
128     printf("rIRecv : min [%d,%2.5f], max[%d,%2.5f] and avg[.,%2.5f] (rms=%2.5f)\n",statRanks[0],rStats[0],statRanks[1],rStats[1],rStats[2],rStats[3]);
129   rgetMinMaxAvg(&mpistats.rWaitAll,rStats,statRanks);
130   if(workfc.myrank==workfc.master)
131     printf("rWtAll : min [%d,%2.5f], max[%d,%2.5f] and avg[.,%2.5f] (rms=%2.5f)\n",statRanks[0],rStats[0],statRanks[1],rStats[1],rStats[2],rStats[3]);
132   rgetMinMaxAvg(&mpistats.rCommu,rStats,statRanks);
133   if(workfc.myrank==workfc.master)
134     printf("rCommu : min [%d,%2.5f], max[%d,%2.5f] and avg[.,%2.5f] (rms=%2.5f)\n",statRanks[0],rStats[0],statRanks[1],rStats[1],rStats[2],rStats[3]);
135   rgetMinMaxAvg(&mpistats.rAllR,rStats,statRanks);
136   if(workfc.myrank==workfc.master) {
137     printf("rAllR  : min [%d,%2.5f], max[%d,%2.5f] and avg[.,%2.5f] (rms=%2.5f)\n",statRanks[0],rStats[0],statRanks[1],rStats[1],rStats[2],rStats[3]);
138     printf("\n");
139   }
140 // Scalars
141   igetMinMaxAvg(&mpistats.iISendScal,iStats,statRanks);
142   if(workfc.myrank==workfc.master)
143     printf("iISendScal : min [%d,%d], max[%d,%d] and avg[.,%d] (rms=%d)\n",statRanks[0],(int)iStats[0],statRanks[1],(int)iStats[1],(int)iStats[2],(int)iStats[3]);
144   igetMinMaxAvg(&mpistats.iIRecvScal,iStats,statRanks);
145   if(workfc.myrank==workfc.master)
146     printf("iIRecvScal : min [%d,%d], max[%d,%d] and avg[.,%d] (rms=%d)\n",statRanks[0],(int)iStats[0],statRanks[1],(int)iStats[1],(int)iStats[2],(int)iStats[3]);
147   igetMinMaxAvg(&mpistats.iWaitAllScal,iStats,statRanks);
148   if(workfc.myrank==workfc.master)
149     printf("iWtAllScal : min [%d,%d], max[%d,%d] and avg[.,%d] (rms=%d)\n",statRanks[0],(int)iStats[0],statRanks[1],(int)iStats[1],(int)iStats[2],(int)iStats[3]);
150   igetMinMaxAvg(&mpistats.iAllRScal,iStats,statRanks);
151   if(workfc.myrank==workfc.master)
152     printf("iAllRScal : min [%d,%d], max[%d,%d] and avg[.,%d] (rms=%d)\n",statRanks[0],(int)iStats[0],statRanks[1],(int)iStats[1],(int)iStats[2],(int)iStats[3]);
153 
154   rgetMinMaxAvg(&mpistats.rISendScal,rStats,statRanks);
155   if(workfc.myrank==workfc.master)
156     printf("rISendScal : min [%d,%2.5f], max[%d,%2.5f] and avg[.,%2.5f] (rms=%2.5f)\n",statRanks[0],rStats[0],statRanks[1],rStats[1],rStats[2],rStats[3]);
157   rgetMinMaxAvg(&mpistats.rIRecvScal,rStats,statRanks);
158   if(workfc.myrank==workfc.master)
159     printf("rIRecvScal : min [%d,%2.5f], max[%d,%2.5f] and avg[.,%2.5f] (rms=%2.5f)\n",statRanks[0],rStats[0],statRanks[1],rStats[1],rStats[2],rStats[3]);
160   rgetMinMaxAvg(&mpistats.rWaitAllScal,rStats,statRanks);
161   if(workfc.myrank==workfc.master)
162     printf("rWtAllScal : min [%d,%2.5f], max[%d,%2.5f] and avg[.,%2.5f] (rms=%2.5f)\n",statRanks[0],rStats[0],statRanks[1],rStats[1],rStats[2],rStats[3]);
163   rgetMinMaxAvg(&mpistats.rCommuScal,rStats,statRanks);
164   if(workfc.myrank==workfc.master)
165     printf("rCommuScal : min [%d,%2.5f], max[%d,%2.5f] and avg[.,%2.5f] (rms=%2.5f)\n",statRanks[0],rStats[0],statRanks[1],rStats[1],rStats[2],rStats[3]);
166   rgetMinMaxAvg(&mpistats.rAllRScal,rStats,statRanks);
167   if(workfc.myrank==workfc.master)
168     printf("rAllRScal  : min [%d,%2.5f], max[%d,%2.5f] and avg[.,%2.5f] (rms=%2.5f)\n",statRanks[0],rStats[0],statRanks[1],rStats[1],rStats[2],rStats[3]);
169 
170 
171 }
172 
173 //void print_system_stats(double tcorecp[2]) {
174 void print_system_stats(double *tcorecp, double *tcorecpscal) {
175   int statRanks[2];
176   double iStats[4], rStats[4];
177   double syst_assembly, syst_solve;
178 
179 // NS equations
180   syst_assembly = tcorecp[0];
181   syst_solve = tcorecp[1];
182 
183   rgetMinMaxAvg(&syst_assembly,rStats,statRanks);
184   if(workfc.myrank==workfc.master)
185     printf("Elm. form. : min [%d,%2.5f], max[%d,%2.5f] and avg[.,%2.5f] (rms=%2.5f)\n",statRanks[0],rStats[0],statRanks[1],rStats[1],rStats[2],rStats[3]);
186 
187   rgetMinMaxAvg(&syst_solve,rStats,statRanks);
188   if(workfc.myrank==workfc.master)
189     printf("Lin. alg. sol : min [%d,%2.5f], max[%d,%2.5f] and avg[.,%2.5f] (rms=%2.5f)\n",statRanks[0],rStats[0],statRanks[1],rStats[1],rStats[2],rStats[3]);
190 
191 // Scalars
192   syst_assembly = tcorecpscal[0];
193   syst_solve = tcorecpscal[1];
194 
195   rgetMinMaxAvg(&syst_assembly,rStats,statRanks);
196   if(workfc.myrank==workfc.master)
197     printf("Elm. form. Scal. : min [%d,%2.5f], max[%d,%2.5f] and avg[.,%2.5f] (rms=%2.5f)\n",statRanks[0],rStats[0],statRanks[1],rStats[1],rStats[2],rStats[3]);
198 
199   rgetMinMaxAvg(&syst_solve,rStats,statRanks);
200   if(workfc.myrank==workfc.master) {
201     printf("Lin. alg. sol Scal. : min [%d,%2.5f], max[%d,%2.5f] and avg[.,%2.5f] (rms=%2.5f)\n",statRanks[0],rStats[0],statRanks[1],rStats[1],rStats[2],rStats[3]);
202     printf("\n");
203   }
204   //printf("rank %d - syst_assembly %f - syst_solve %f\n",workfc.myrank,syst_assembly,syst_solve);
205 }
206 
207 
208 
209 void countfieldstowriterestart()
210 {
211   int nfields;
212 
213 //     printf("TEST: %d %d %d %d %d\n",timdat.istep,timdat.itseq,inpdat.nstep[0],inpdat.nstep[1],timdat.lstep);
214 
215   nfields = 2; //solution, time derivatives
216 
217   if(outpar.ivort == 1){
218     nfields++; //vorticity
219   }
220 
221   if(abs(turbvar.itwmod) != 1 && outpar.iowflux == 1) {
222     nfields++; //instantaneous wss in bflux.f
223   }
224 
225 //   if(ideformwall.eq.1) not handled yet
226 
227 // KEJ says save all with this version
228 // BEFORE we do animations, we will want to have a "saveall" frequency flag
229 // if(timdat.istep == inpdat.nstep[timdat.itseq-1]){ //Last time step of the computation
230 
231     //projection vectors and pressure projection vectors (call saveLesRestart in itrdrv)
232     nfields = nfields +2;
233 
234     //if Print Error Indicators = true (call write_error in itrdrv)
235     if(turbvar.ierrcalc == 1){
236       nfields++;
237     }
238 
239     //if Print ybar = True (call write_field(myrank,'a','ybar',4,... in itrdrv)
240     if(outpar.ioybar == 1){
241       nfields++;  //ybar
242 
243       //phase average fields
244       if(outpar.nphasesincycle >0) {
245         nfields = nfields + outpar.nphasesincycle;
246       }
247 
248       if(abs(turbvar.itwmod) != 1 && outpar.iowflux == 1) {
249         nfields++; //wssbar
250       }
251 
252     }
253 
254     if(turbvari.irans < 0) {
255       nfields++; //dwal
256     }
257 
258 //  }
259 
260   outpar.nsynciofieldswriterestart = nfields;
261 
262   if(workfc.myrank == 0) {
263     printf("Number of fields to write in restart files: %d\n", nfields);
264   }
265 }
266 
267 
268 void
269 Write_Restart(  int* pid,
270                 int* stepno,
271                 int* nshg,
272                 int* numVars,
273                 double* array1,
274                 double* array2 ) {
275 
276     char fname[255];
277     char rfile[60];
278     char existingfile[30], linkfile[30];
279     int irstou;
280     int magic_number = 362436;
281     int* mptr = &magic_number;
282 //    time_t timenow = time ( &timenow);
283     double version=0.0;
284     int isize, nitems;
285     int iarray[10];
286 
287     /*sprintf(rfile,"restart.%d.%d",*stepno,*pid+1);
288     openfile_(rfile,"write", &irstou);
289 
290     // writing the top ascii header for the restart file
291 
292     writestring_( &irstou,"# PHASTA Input File Version 2.0\n");
293     writestring_( &irstou,
294                   "# format \"keyphrase : sizeofnextblock usual headers\"\n");
295 
296     bzero( (void*)fname, 255 );
297     sprintf(fname,"# Output generated by phasta version (NOT YET CURRENT): %lf \n", version);
298     writestring_( &irstou, fname );
299 
300     bzero( (void*)fname, 255 );
301     gethostname(fname,255);
302     writestring_( &irstou,"# This result was produced on: ");
303     writestring_( &irstou, fname );
304     writestring_( &irstou,"\n");
305 
306     bzero( (void*)fname, 255 );
307     sprintf(fname,"# %s\n", ctime( &timenow ));
308     writestring_( &irstou, fname );
309 
310     isize = 1;
311     nitems = 1;
312     iarray[ 0 ] = 1;
313     writeheader_( &irstou, "byteorder magic number ",
314                   (void*)iarray, &nitems, &isize, "integer", phasta_iotype );
315 
316     nitems = 1;
317     writedatablock_( &irstou, "byteorder magic number ",
318                      (void*)mptr, &nitems, "integer", phasta_iotype );
319 
320 
321     bzero( (void*)fname, 255 );
322     sprintf(fname,"number of modes : < 0 > %d\n", *nshg);
323     writestring_( &irstou, fname );
324 
325     bzero( (void*)fname, 255 );
326     sprintf(fname,"number of variables : < 0 > %d\n", *numVars);
327     writestring_( &irstou, fname );
328 
329 
330     isize = (*nshg)*(*numVars);
331     nitems = 3;
332     iarray[ 0 ] = (*nshg);
333     iarray[ 1 ] = (*numVars);
334     iarray[ 2 ] = (*stepno);
335     writeheader_( &irstou, "solution ",
336                   (void*)iarray, &nitems, &isize, "double", phasta_iotype );
337 
338 
339     nitems = (*nshg)*(*numVars);
340     writedatablock_( &irstou, "solution ",
341                      (void*)(array1), &nitems, "double", phasta_iotype );
342 
343 
344 
345     nitems = 3;
346     writeheader_( &irstou, "time derivative of solution ",
347                   (void*)iarray, &nitems, &isize, "double", phasta_iotype );
348 
349 
350     nitems = (*nshg)*(*numVars);
351     writedatablock_( &irstou, "time derivative of solution ",
352                      (void*)(array2), &nitems, "double", phasta_iotype );
353 
354 
355     closefile_( &irstou, "write" );
356     */
357     //MPI_Barrier(MPI_COMM_WORLD);
358 
359     /////////////////////////////// Start of writing using new-lib ////////////////////////////
360 
361 //MR CHANGE
362     int nfiles;
363     int nfields;
364     int numparts;
365     int irank;
366     int nprocs;
367 
368     //  First, count the number of fields to write and store the result in
369     countfieldstowriterestart();
370 
371     //  Retrieve and compute the parameters required for SyncIO
372     nfiles = outpar.nsynciofiles;
373     nfields = outpar.nsynciofieldswriterestart;
374     numparts = workfc.numpe;
375     irank = *pid; //workfc.myrank;
376     nprocs = workfc.numpe;
377 //MR CHANGE END
378     int nppf = numparts/nfiles;
379     int GPID;
380 
381     // Calculate number of parts each proc deal with and where it start and end ...
382     int nppp = numparts/nprocs;// nppp : Number of parts per proc ...
383     int startpart = irank * nppp +1;// Part id from which I (myrank) start ...
384     int endpart = startpart + nppp - 1;// Part id to which I (myrank) end ...
385 
386     int descriptor;
387     char filename[255],path[255],fieldtag_s[255];
388     bzero((void*)filename,255);
389     bzero((void*)fieldtag_s,255);
390 
391     sprintf(filename,"restart-dat.%d.%d",*stepno,((int)(irank/(nprocs/nfiles))+1));
392 
393 //    unsigned long long timer_start;
394 //    unsigned long long timer_end;
395 //    double time_span;
396 
397 //MR CHANGE
398 //  Measure the time - Start the timer
399 //    MPI_Barrier(MPI_COMM_WORLD);
400 //    timer_start = rdtsc();
401 //MR CHANGE END
402 
403     initphmpiio(&nfields, &nppf, &nfiles, &f_descriptor, "write");
404 
405 //MR CHANGE
406 //  Measure the time - End of timer
407 //    timer_end = rdtsc();
408 //    time_span=(double)((timer_end-timer_start)/cpu_speed);
409     if (*pid==0) {
410 //      printf("\n*****************************\n");
411 //      printf("Time: 'initphmpiio' of %s with %d fields and %d files is:    %f s\n",filename,nfields,nfiles,time_span);
412       printf("Filename is %s \n",filename);
413     }
414 //MR CHANGE END
415 
416 
417 //MR CHANGE
418 //  Measure the time - Start the timer
419 //    MPI_Barrier(MPI_COMM_WORLD);
420 //    timer_start = rdtsc();
421 //MR CHANGE END
422 
423     openfile(filename, "write", &f_descriptor);
424 
425 //MR CHANGE
426 //  Measure the time - End of timer
427 //    MPI_Barrier(MPI_COMM_WORLD);
428 //    timer_end = rdtsc();
429 //    time_span=(double)((timer_end-timer_start)/cpu_speed);
430 //    if (*pid==0) {
431 //      printf("Time: 'openfile' of %s with %d fields and %d files is:    %f s\n",filename,nfields,nfiles,time_span);
432 //      printf("*****************************\n");
433 //    }
434 //MR CHANGE END
435 
436     field_flag=0;
437 
438      int i;
439      for ( i = 0; i < nppp; i++) { //This loop is useful only if several parts per processor
440      // GPID : global part id, corresponds to rank ...
441         // e.g : (in this example)
442         // proc 0 : 1--4
443         // proc 1 : 5--8 ...
444         GPID = startpart + i;
445 
446         // Write solution field ...
447         sprintf(fieldtag_s,"solution@%d",GPID);
448 
449         isize = (*nshg)*(*numVars);
450         nitems = 3;
451         iarray[ 0 ] = (*nshg);
452         iarray[ 1 ] = (*numVars);
453         iarray[ 2 ] = (*stepno);
454 
455 //MR CHANGE
456 //  Measure the time - Start the timer
457 //        MPI_Barrier(MPI_COMM_WORLD);
458 //        timer_start = rdtsc();
459 //MR CHANGE END
460 
461         writeheader( &f_descriptor, fieldtag_s, (void*)iarray, &nitems, &isize, "double", phasta_iotype);
462 
463 //MR CHANGE
464 //  Measure the time - End of timer
465 //        MPI_Barrier(MPI_COMM_WORLD);
466 //        timer_end = rdtsc();
467 //        time_span=(double)((timer_end-timer_start)/cpu_speed);
468 //        if (*pid==0) {
469 //          printf("\n*****************************\n");
470 //          printf("Time: header for 'Solution':    %f s\n",time_span);
471 //        }
472 //MR CHANGE END
473 
474         nitems = (*nshg)*(*numVars);
475 
476 //MR CHANGE
477 //  Measure the time - Start the timer
478 //        MPI_Barrier(MPI_COMM_WORLD);
479 //        timer_start = rdtsc();
480 //MR CHANGE END
481 
482         writedatablock( &f_descriptor, fieldtag_s, (void*)(array1), &isize, "double", phasta_iotype );
483 
484 //MR CHANGE
485 //  Measure the time - End of timer
486 //        MPI_Barrier(MPI_COMM_WORLD);
487 //        timer_end = rdtsc();
488 //        time_span=(double)((timer_end-timer_start)/cpu_speed);
489 
490 //        int isizemin,isizemax,isizetot;
491 //        double sizemin,sizemax,sizeavg,sizetot,rate;
492 
493 //        MPI_Allreduce(&isize,&isizemin,1,MPI_INT,MPI_MIN,MPI_COMM_WORLD);
494 //        MPI_Allreduce(&isize,&isizemax,1,MPI_INT,MPI_MAX,MPI_COMM_WORLD);
495 //        MPI_Allreduce(&isize,&isizetot,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
496 
497 //        sizemin=(double)(8.0*isizemin/1024.0/1024.0);
498 //        sizemax=(double)(8.0*isizemax/1024.0/1024.0);
499 //        sizetot=(double)(8.0*isizetot/1024.0/1024.0);
500 //        sizeavg=(double)(1.0*sizetot/workfc.numpe);
501 //        rate=(double)(1.0*sizetot/time_span);
502 
503 //        if (*pid==0) {
504 //          printf("Time: block for 'Solution':    %f s\n",time_span);
505 //          printf("Time: block:   Min= %f MB; Max= %f MB; Avg= %f MB; Tot= %f MB; Rate= %f MB/s; \n",sizemin,sizemax,sizeavg,sizetot,rate);
506 
507 //        }
508 //MR CHANGE END
509 
510     }
511     field_flag++;
512 
513     for ( i = 0; i < nppp; i++) {
514 
515         // GPID : global part id, corresponds to rank ...
516         // e.g : (in this example)
517         // proc 0 : 1--4
518         // proc 1 : 5--8 ...
519         GPID = startpart + i;
520 
521         // Write solution field ...
522         sprintf(fieldtag_s,"time derivative of solution@%d",GPID);
523 
524         isize = (*nshg)*(*numVars);
525         nitems = 3;
526         iarray[ 0 ] = (*nshg);
527         iarray[ 1 ] = (*numVars);
528         iarray[ 2 ] = (*stepno);
529 
530 //MR CHANGE
531 //  Measure the time - Start the timer
532 //        MPI_Barrier(MPI_COMM_WORLD);
533 //        timer_start = rdtsc();
534 //MR CHANGE END
535 
536         writeheader( &f_descriptor, fieldtag_s, (void*)iarray, &nitems, &isize, "double", phasta_iotype);
537 
538 //MR CHANGE
539 //  Measure the time - End of timer
540 //        MPI_Barrier(MPI_COMM_WORLD);
541 //        timer_end = rdtsc();
542 //        time_span=(double)((timer_end-timer_start)/cpu_speed);
543 //        if (*pid==0) {
544 //          printf("Time: header for 'Time Derivative of solution':    %f s\n",time_span);
545 //        }
546 //MR CHANGE END
547 
548         nitems = (*nshg)*(*numVars);
549 
550 //MR CHANGE
551 //  Measure the time - Start the timer
552 //        MPI_Barrier(MPI_COMM_WORLD);
553 //        timer_start = rdtsc();
554 //MR CHANGE END
555 
556         writedatablock( &f_descriptor, fieldtag_s, (void*)(array2), &isize, "double", phasta_iotype );
557 
558 //MR CHANGE
559 //  Measure the time - End of timer
560 //        MPI_Barrier(MPI_COMM_WORLD);
561 //        timer_end = rdtsc();
562 //        time_span=(double)((timer_end-timer_start)/cpu_speed);
563 
564 //        int isizemin,isizemax,isizetot;
565 //        double sizemin,sizemax,sizeavg,sizetot,rate;
566 
567 //        MPI_Allreduce(&isize,&isizemin,1,MPI_INT,MPI_MIN,MPI_COMM_WORLD);
568 //        MPI_Allreduce(&isize,&isizemax,1,MPI_INT,MPI_MAX,MPI_COMM_WORLD);
569 //        MPI_Allreduce(&isize,&isizetot,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
570 
571 //        sizemin=(double)(8.0*isizemin/1024.0/1024.0);
572 //        sizemax=(double)(8.0*isizemax/1024.0/1024.0);
573 //        sizetot=(double)(8.0*isizetot/1024.0/1024.0);
574 //        sizeavg=sizetot/workfc.numpe;
575 //        rate=sizetot/time_span;
576 
577 //        if (*pid==0) {
578 //          printf("Time: block for 'Time Derivative of Solution':    %f s\n",time_span);
579 //          printf("Time: block:   Min= %f MB; Max= %f MB; Avg= %f MB; Tot= %f MB; Rate= %f MB/s; \n",sizemin,sizemax,sizeavg,sizetot,rate);
580 //          printf("*****************************\n");
581 
582 //        }
583 //MR CHANGE END
584 
585     }
586     field_flag++;
587 
588     if (field_flag==nfields){
589 
590 //MR CHANGE
591 //  Measure the time - Start the timer
592 //      MPI_Barrier(MPI_COMM_WORLD);
593 //      timer_start = rdtsc();
594 //MR CHANGE END
595 
596       closefile(&f_descriptor, "write");
597 
598 //MR CHANGE
599 //    Measure the time - End of timer
600 //      MPI_Barrier(MPI_COMM_WORLD);
601 //      timer_end = rdtsc();
602 //      time_span=(double)((timer_end-timer_start)/cpu_speed);
603 //      if (*pid==0) {
604 //        printf("\n*****************************\n");
605 //        printf("Time: 'closefile' is:    %f s\n",time_span);
606 //      }
607 //MR CHANGE END
608 
609 //MR CHANGE
610 //  Measure the time - Start the timer
611 //      MPI_Barrier(MPI_COMM_WORLD);
612 //      timer_start = rdtsc();
613 //MR CHANGE END
614 
615       finalizephmpiio(&f_descriptor);
616 
617 //MR CHANGE
618 //    Measure the time - End of timer
619 //      MPI_Barrier(MPI_COMM_WORLD);
620 //      timer_end = rdtsc();
621 //      time_span=(double)((timer_end-timer_start)/cpu_speed);
622       if (*pid==0) {
623 //        printf("Time: 'finalizephmpiio' is:    %f s\n",time_span);
624 //        printf("Last field %d '%s' finished! \n",nfields, fieldtag_s);
625         printf("\n");
626 //        printf("*****************************\n");
627       }
628     }
629 //MR CHANGE END
630 
631 
632 
633     ///////////////////////////////////////////////////////////////////////////////////////////
634 
635     /* create a soft link of the restart we just wrote to restart.latest
636      this is the file the next run will always try to start from */
637 
638 /*    sprintf( linkfile, "restart.latest.%d", *pid+1 );
639     unlink( linkfile );
640     sprintf( existingfile, "restart.%d.%d", *stepno, *pid+1 );
641     link( existingfile, linkfile );
642 */
643 }
644 
645 void
646 Write_Error(  int* pid,
647               int* stepno,
648               int* nshg,
649               int* numVars,
650               double* array1 ) {
651 
652 
653     char fname[255];
654     char rfile[60];
655     int irstou;
656     int magic_number = 362436;
657     int* mptr = &magic_number;
658     //printf("Time is commented\n");
659     //time_t timenow = time ( &timenow);
660     //printf("Yes\n");
661     double version=0.0;
662     int isize, nitems;
663     int iarray[10];
664 
665     /*sprintf(rfile,"restart.%d.%d",*stepno,*pid+1);
666     openfile_(rfile,"append", &irstou);
667 
668     isize = (*nshg)*(*numVars);
669     nitems = 3;
670     iarray[ 0 ] = (*nshg);
671     iarray[ 1 ] = (*numVars);
672     iarray[ 2 ] = (*stepno);
673     writeheader_( &irstou, "errors", (void*)iarray, &nitems, &isize, "double", phasta_iotype );
674 
675 
676     nitems = (*nshg)*(*numVars);
677     writedatablock_( &irstou, "errors ", (void*)(array1), &nitems, "double", phasta_iotype );
678 
679     closefile_( &irstou, "append" );*/
680 
681     /////////////////////////////// Start of writing using new-lib ////////////////////////////
682 
683     int nfiles;
684     int nfields;
685     int numparts;
686     int irank;
687     int nprocs;
688 
689 //    unsigned long long timer_start;
690 //    unsigned long long timer_end;
691 //    double time_span;
692 
693     nfiles = outpar.nsynciofiles;
694     nfields = outpar.nsynciofieldswriterestart;
695     numparts = workfc.numpe;
696     irank = *pid; //workfc.myrank;
697     nprocs = workfc.numpe;
698 
699     int nppf = numparts/nfiles;
700     int GPID;
701 
702     // Calculate number of parts each  proc deal with and where it start and end ...
703     int nppp = numparts/nprocs;// nppp : Number of parts per proc ...
704     int startpart = irank * nppp +1;// Part id from which I (myrank) start ...
705     int endpart = startpart + nppp - 1;// Part id to which I (myrank) end ...
706 
707     field_flag++;
708 
709     char fieldtag[255];
710 
711     int i;
712     for ( i = 0; i < nppp; i++  ) {
713         GPID = startpart + i;
714         sprintf(fieldtag,"errors@%d",GPID);
715 
716         if(*pid==0) {
717 //          printf("\n*****************************\n");
718           printf("\n");
719           printf("The %d/%d th field to be written is '%s'\n",field_flag,nfields,fieldtag);
720         }
721 
722         isize = (*nshg)*(*numVars);
723         nitems = 3;
724         iarray[ 0 ] = (*nshg);
725         iarray[ 1 ] = (*numVars);
726         iarray[ 2 ] = (*stepno);
727 
728 //MR CHANGE
729 //  Measure the time - Start the timer
730 //        MPI_Barrier(MPI_COMM_WORLD);
731 //        timer_start = rdtsc();
732 //MR CHANGE END
733 
734         writeheader( &f_descriptor, fieldtag, (void*)iarray, &nitems, &isize, "double", phasta_iotype);
735 
736 //MR CHANGE
737 //  Measure the time - End of timer
738 //        MPI_Barrier(MPI_COMM_WORLD);
739 //        timer_end = rdtsc();
740 //        time_span=(double)((timer_end-timer_start)/cpu_speed);
741 //        if (*pid==0) {
742 //          printf("Time: header for 'error':    %f s\n",time_span);
743 //        }
744 //MR CHANGE END
745 
746 //MR CHANGE
747 //  Measure the time - Start the timer
748 //        MPI_Barrier(MPI_COMM_WORLD);
749 //        timer_start = rdtsc();
750 //MR CHANGE END
751 
752         writedatablock( &f_descriptor, fieldtag, (void*)array1, &isize, "double", phasta_iotype );
753 
754 //MR CHANGE
755 //  Measure the time - End of timer
756 //        MPI_Barrier(MPI_COMM_WORLD);
757 //        timer_end = rdtsc();
758 //        time_span=(double)((timer_end-timer_start)/cpu_speed);
759 
760 //        int isizemin,isizemax,isizetot;
761 //        double sizemin,sizemax,sizeavg,sizetot,rate;
762 
763 //        MPI_Allreduce(&isize,&isizemin,1,MPI_INT,MPI_MIN,MPI_COMM_WORLD);
764 //        MPI_Allreduce(&isize,&isizemax,1,MPI_INT,MPI_MAX,MPI_COMM_WORLD);
765 //        MPI_Allreduce(&isize,&isizetot,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
766 
767 //        sizemin=(double)(8.0*isizemin/1024.0/1024.0);
768 //        sizemax=(double)(8.0*isizemax/1024.0/1024.0);
769 //        sizetot=(double)(8.0*isizetot/1024.0/1024.0);
770 //        sizeavg=sizetot/workfc.numpe;
771 //        rate=sizetot/time_span;
772 
773 //        if (*pid==0) {
774 //          printf("Time: block for 'error':    %f s\n",time_span);
775 //          printf("Time: block:   Min= %f MB; Max= %f MB; Avg= %f MB; Tot= %f MB; Rate= %f MB/s; \n",sizemin,sizemax,sizeavg,sizetot,rate);
776 //          printf("*****************************\n");
777 //        }
778 //MR CHANGE END
779 
780     }
781 
782 //     MPI_Barrier(MPI_COMM_WORLD);
783 //     timer_end = rdtsc();
784 //     time_span=(double)(timer_end-timer_start)/cpu_speed;
785 
786 //     if (*pid==0) {
787 //         printf("Field 'error' written in:     %f s\n",time_span);
788 //         printf("Write field '%s' finished! \n",fieldtag);
789 //     }
790 
791     if (field_flag==nfields){
792 
793 //MR CHANGE
794 //  Measure the time - Start the timer
795 //      MPI_Barrier(MPI_COMM_WORLD);
796 //      timer_start = rdtsc();
797 //MR CHANGE END
798 
799       closefile(&f_descriptor, "write");
800 
801 //MR CHANGE
802 //    Measure the time - End of timer
803 //      MPI_Barrier(MPI_COMM_WORLD);
804 //      timer_end = rdtsc();
805 //      time_span=(double)((timer_end-timer_start)/cpu_speed);
806 //      if (*pid==0) {
807 //        printf("\n*****************************\n");
808 //        printf("Time: 'closefile' is:    %f s\n",time_span);
809 //      }
810 //MR CHANGE END
811 
812 //MR CHANGE
813 //  Measure the time - Start the timer
814 //      MPI_Barrier(MPI_COMM_WORLD);
815 //      timer_start = rdtsc();
816 //MR CHANGE END
817 
818       finalizephmpiio(&f_descriptor);
819 
820 //MR CHANGE
821 //    Measure the time - End of timer
822 //      MPI_Barrier(MPI_COMM_WORLD);
823 //      timer_end = rdtsc();
824 //      time_span=(double)((timer_end-timer_start)/cpu_speed);
825       if (*pid==0) {
826 //        printf("Time: 'finalizephmpiio' is:    %f s\n",time_span);
827         printf("Last field %d '%s' finished! \n",nfields, fieldtag);
828         printf("\n");
829 //        printf("*****************************\n");
830       }
831     }
832 //MR CHANGE END
833 
834     ///////////////////////////////////////////////////////////////////////////////////////////
835 
836 
837 }
838 
839 
840 void
841 Write_Displ(  int* pid,
842               int* stepno,
843               int* nshg,
844               int* numVars,
845               double* array1 ) {
846 
847 
848     char fname[255];
849     char rfile[60];
850     int irstou;
851     int magic_number = 362436;
852     int* mptr = &magic_number;
853     time_t timenow = time ( &timenow);
854     double version=0.0;
855     int isize, nitems;
856     int iarray[10];
857 
858     sprintf(rfile,"restart.%d.%d",*stepno,*pid+1);
859     openfile(rfile,"append", &irstou);
860 
861     isize = (*nshg)*(*numVars);
862     nitems = 3;
863     iarray[ 0 ] = (*nshg);
864     iarray[ 1 ] = (*numVars);
865     iarray[ 2 ] = (*stepno);
866     writeheader( &irstou, "displacement", (void*)iarray, &nitems, &isize, "double", phasta_iotype );
867 
868     nitems = (*nshg)*(*numVars);
869     writedatablock( &irstou, "displacement", (void*)(array1), &nitems, "double", phasta_iotype );
870 
871     closefile( &irstou, "append" );
872 }
873 
874 void
875 Write_Field(  int *pid,
876               char* filemode,
877               char* fieldtag,
878               int* tagsize,
879               void* array,
880               char* arraytype,
881               int* nshg,
882               int* numvars,
883               int* stepno) {
884 
885     //printf("Rank is %d, field is %s, tagsize is %d, nshg is %d, numvars is %d\n",*pid,fieldtag,*tagsize,*nshg,*numvars);
886 
887 //     char rfile[32];
888     // assuming restart.sn.(pid+1)
889 //     sprintf(rfile,"restart.%d.%d",*stepno,*pid+1);
890 
891     char *fieldlabel = (char *)malloc((*tagsize+1)*sizeof(char));
892     strncpy(fieldlabel, fieldtag, *tagsize);
893     fieldlabel[*tagsize] = '\0';
894 
895     int irstou;
896     int magic_number = 362436;
897     int* mptr = &magic_number;
898     double version=0.0;
899     int isize, nitems;
900     int iarray[10];
901 
902     char fmode[10];
903     if(!strncmp(filemode,"w",1))
904       strcpy(fmode,"write");
905     else // default is append
906       strcpy(fmode,"append");
907 
908     char datatype[10];
909     if(!strncmp(arraytype,"i",1))
910       strcpy(datatype,"int");
911     else // default is double
912       strcpy(datatype,"double");
913 
914 /*     openfile_(rfile, fmode, &irstou);
915 
916      nitems = 3; // assuming field will write 3 items in iarray
917      iarray[ 0 ] = (*nshg);
918      iarray[ 1 ] = (*numvars);
919      iarray[ 2 ] = (*stepno);
920 
921      isize = (*nshg)*(*numvars);
922      writeheader_( &irstou, fieldlabel, (void*)iarray, &nitems, &isize, datatype, phasta_iotype );
923 
924      nitems = (*nshg)*(*numvars);
925      writedatablock_( &irstou, fieldlabel, array, &nitems, datatype, phasta_iotype );
926      closefile_( &irstou, fmode);
927 */
928     /////////////////////////////// Start of writing using new-lib ////////////////////////////
929 
930     int nfiles;
931     int nfields;
932     int numparts;
933     int irank;
934     int nprocs;
935 
936 //    unsigned long long timer_start;
937 //    unsigned long long timer_end;
938 //    double time_span;
939 
940     nfiles = outpar.nsynciofiles;
941     nfields = outpar.nsynciofieldswriterestart;
942     numparts = workfc.numpe;
943     irank = *pid; //workfc.myrank;
944     nprocs = workfc.numpe;
945 
946     int nppf = numparts/nfiles;
947     int GPID;
948 
949     // Calculate number of parts each  proc deal with and where it start and end ...
950     int nppp = numparts/nprocs;// nppp : Number of parts per proc ...
951     int startpart = irank * nppp +1;// Part id from which I (myrank) start ...
952     int endpart = startpart + nppp - 1;// Part id to which I (myrank) end ...
953 
954     char filename[255],path[255],fieldtag_s[255];
955     bzero((void*)filename,255);
956     bzero((void*)fieldtag_s,255);
957 
958     strncpy(fieldlabel, fieldtag, *tagsize);
959 
960     field_flag++;
961     if(*pid==0) {
962 //      printf("\n*****************************\n");
963       printf("\n");
964       printf("The %d/%d th field to be written is '%s'\n",field_flag,nfields,fieldlabel);
965     }
966 
967     sprintf(filename,"restart-dat.%d.%d",*stepno,((int)(irank/(nprocs/nfiles))+1));
968 
969 //     MPI_Barrier(MPI_COMM_WORLD);
970 //     timer_start = rdtsc();
971 
972     int i;
973     for ( i = 0; i < nppp; i++  ) {
974         GPID = startpart + i;
975 
976         // Write solution field ...
977         sprintf(fieldtag_s,"%s@%d",fieldlabel,GPID);
978 
979         isize = (*nshg)*(*numvars);
980         nitems = 3;
981         iarray[ 0 ] = (*nshg);
982         iarray[ 1 ] = (*numvars);
983         iarray[ 2 ] = (*stepno);
984 
985 //MR CHANGE
986 //  Measure the time - Start the timer
987 //        MPI_Barrier(MPI_COMM_WORLD);
988 //        timer_start = rdtsc();
989 //MR CHANGE END
990 
991         writeheader( &f_descriptor, fieldtag_s, (void*)iarray, &nitems, &isize, datatype, phasta_iotype);
992 
993 //MR CHANGE
994 //  Measure the time - End of timer
995 //        MPI_Barrier(MPI_COMM_WORLD);
996 //        timer_end = rdtsc();
997 //        time_span=(double)((timer_end-timer_start)/cpu_speed);
998 //        if (*pid==0) {
999 //          printf("Time: header for '%s':    %f s\n",fieldtag_s,time_span);
1000 //        }
1001 //MR CHANGE END
1002 
1003         nitems = (*nshg)*(*numvars);
1004 
1005 //MR CHANGE
1006 //  Measure the time - Start the timer
1007 //        MPI_Barrier(MPI_COMM_WORLD);
1008 //        timer_start = rdtsc();
1009 //MR CHANGE END
1010 
1011         writedatablock( &f_descriptor, fieldtag_s, array, &isize, datatype, phasta_iotype );
1012 
1013 //MR CHANGE
1014 //  Measure the time - End of timer
1015 //        MPI_Barrier(MPI_COMM_WORLD);
1016 //        timer_end = rdtsc();
1017 //        time_span=(double)((timer_end-timer_start)/cpu_speed);
1018 
1019 //        int isizemin,isizemax,isizetot;
1020 //        double sizemin,sizemax,sizeavg,sizetot,rate;
1021 
1022 //        MPI_Allreduce(&isize,&isizemin,1,MPI_INT,MPI_MIN,MPI_COMM_WORLD);
1023 //        MPI_Allreduce(&isize,&isizemax,1,MPI_INT,MPI_MAX,MPI_COMM_WORLD);
1024 //        MPI_Allreduce(&isize,&isizetot,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
1025 
1026 //        sizemin=(double)(8.0*isizemin/1024.0/1024.0);
1027 //        sizemax=(double)(8.0*isizemax/1024.0/1024.0);
1028 //        sizetot=(double)(8.0*isizetot/1024.0/1024.0);
1029 //        sizeavg=sizetot/workfc.numpe;
1030 //        rate=sizetot/time_span;
1031 
1032 //        if (*pid==0) {
1033 //          printf("Time: block for '%s':    %f s\n",fieldtag_s,time_span);
1034 //          printf("Time: block:   Min= %f MB; Max= %f MB; Avg= %f MB; Tot= %f MB; Rate= %f MB/s; \n",sizemin,sizemax,sizeavg,sizetot,rate);
1035 //          printf("*****************************\n");
1036 //        }
1037 //MR CHANGE END
1038 
1039     }
1040 
1041 //     MPI_Barrier(MPI_COMM_WORLD);
1042 //     timer_end = rdtsc();
1043 //     time_span=(double)(timer_end-timer_start)/cpu_speed;
1044 
1045 //     if (*pid==0) {
1046 //         printf("Field '%s' written in:     %f s\n",fieldtag,time_span);
1047 //         printf("Write field '%s' finished! \n",fieldtag_s);
1048 //     }
1049 
1050 //     if (field_flag==nfields){
1051 //       closefile_(&f_descriptor, "write");
1052 //       finalizephmpiio_(&f_descriptor);
1053 //       if(*pid==0) {
1054 //         printf("Last field %d '%s' finished! \n",nfields, fieldtag_s);
1055 //         printf("\n*****************************\n");
1056 //       }
1057 //     }
1058 
1059     if (field_flag==nfields){
1060 
1061 //MR CHANGE
1062 //  Measure the time - Start the timer
1063 //      MPI_Barrier(MPI_COMM_WORLD);
1064 //      timer_start = rdtsc();
1065 //MR CHANGE END
1066 
1067       closefile(&f_descriptor, "write");
1068 
1069 //MR CHANGE
1070 //    Measure the time - End of timer
1071 //      MPI_Barrier(MPI_COMM_WORLD);
1072 //      timer_end = rdtsc();
1073 //      time_span=(double)((timer_end-timer_start)/cpu_speed);
1074 //      if (*pid==0) {
1075 //        printf("\n*****************************\n");
1076 //        printf("Time: 'closefile' is:    %f s\n",time_span);
1077 //      }
1078 //MR CHANGE END
1079 
1080 //MR CHANGE
1081 //  Measure the time - Start the timer
1082 //      MPI_Barrier(MPI_COMM_WORLD);
1083 //      timer_start = rdtsc();
1084 //MR CHANGE END
1085 
1086       finalizephmpiio(&f_descriptor);
1087 
1088 //MR CHANGE
1089 //    Measure the time - End of timer
1090 //      MPI_Barrier(MPI_COMM_WORLD);
1091 //      timer_end = rdtsc();
1092 //      time_span=(double)((timer_end-timer_start)/cpu_speed);
1093       if (*pid==0) {
1094 //        printf("Time: 'finalizephmpiio' is:    %f s\n",time_span);
1095         printf("Last field %d '%s' finished! \n",nfields, fieldtag);
1096         printf("\n");
1097 //        printf("*****************************\n");
1098       }
1099     }
1100 //MR CHANGE END
1101 
1102     ///////////////////////////////////////////////////////////////////////////////////////////
1103 
1104     free(fieldlabel);
1105 }
1106 
1107 //MR CHANGE
1108 void
1109 Write_PhAvg(  int* pid,
1110               char* filemode,
1111               char* fieldtag,
1112               int* tagsize,
1113               int* iphase,
1114               void* array,
1115               char* arraytype,
1116               int* nshg,
1117               int* numvars,
1118               int* stepno) {
1119 
1120     char rfile[32];
1121     // assuming restart_phase_avg_<sn>.<iphase>.<pid+1>
1122     sprintf(rfile,"restart_phase_avg_%d.%d.%d",*stepno,*iphase,*pid+1);
1123 
1124     char *fieldlabel = (char *)malloc((*tagsize+1)*sizeof(char));
1125     strncpy(fieldlabel, fieldtag, *tagsize);
1126     fieldlabel[*tagsize] = '\0';
1127 
1128     int irstou;
1129     int isize, nitems;
1130     int iarray[10];
1131 
1132     char fmode[10];
1133     if(!strncmp(filemode,"w",1))
1134       strcpy(fmode,"write");
1135     else // default is append
1136       strcpy(fmode,"append");
1137 
1138     char datatype[10];
1139     if(!strncmp(arraytype,"i",1))
1140       strcpy(datatype,"int");
1141     else // default is double
1142       strcpy(datatype,"double");
1143 
1144     openfile(rfile, fmode, &irstou);
1145 
1146     if(!strcmp(fmode,"write")) {
1147       // may be create a routine for 'top' portion under write mode
1148       int magic_number = 362436;
1149       int* mptr = &magic_number;
1150       time_t timenow = time ( &timenow);
1151       double version=0.0;
1152 
1153       /* writing the top ascii header for the restart file */
1154 
1155       writestring( &irstou,"# PHASTA Input File Version 2.0\n");
1156       writestring( &irstou,
1157                     "# format \"keyphrase : sizeofnextblock usual headers\"\n");
1158 
1159       char fname[255];
1160       bzero( (void*)fname, 255 );
1161       sprintf(fname,"# Output generated by phasta version (NOT YET CURRENT): %lf \n", version);
1162       writestring( &irstou, fname );
1163 
1164       bzero( (void*)fname, 255 );
1165       gethostname(fname,255);
1166       writestring( &irstou,"# This result was produced on: ");
1167       writestring( &irstou, fname );
1168       writestring( &irstou,"\n");
1169 
1170       bzero( (void*)fname, 255 );
1171       sprintf(fname,"# %s\n", ctime( &timenow ));
1172       writestring( &irstou, fname );
1173 
1174       isize = 1;
1175       nitems = 1;
1176       iarray[ 0 ] = 1;
1177       writeheader( &irstou, "byteorder magic number ",
1178                     (void*)iarray, &nitems, &isize, "integer", phasta_iotype );
1179       nitems = 1;
1180       writedatablock( &irstou, "byteorder magic number ",
1181                        (void*)mptr, &nitems, "integer", phasta_iotype );
1182     }
1183 
1184     nitems = 3; // assuming field will write 3 items in iarray
1185     iarray[ 0 ] = (*nshg);
1186     iarray[ 1 ] = (*numvars);
1187     iarray[ 2 ] = (*stepno);
1188 
1189     isize = (*nshg)*(*numvars);
1190     writeheader( &irstou, fieldlabel, (void*)iarray, &nitems, &isize, datatype, phasta_iotype );
1191 
1192     nitems = (*nshg)*(*numvars);
1193     writedatablock( &irstou, fieldlabel, array, &nitems, datatype, phasta_iotype );
1194 
1195     closefile( &irstou, fmode);
1196 
1197     free(fieldlabel);
1198 }
1199 
1200 //MR CHANGE
1201 void
1202 Write_PhAvg2( int* pid,
1203               char* filemode,
1204               char* fieldtag,
1205               int* tagsize,
1206               int* iphase,
1207               int* nphasesincycle,
1208               void* array,
1209               char* arraytype,
1210               int* nshg,
1211               int* numvars,
1212               int* stepno) {
1213 
1214 //     char rfile[32];
1215     // assuming restart.sn.(pid+1)
1216 //     sprintf(rfile,"restart.%d.%d",*stepno,*pid+1);
1217 
1218     int addtagsize; // phase number is added to the name of the field
1219     if(*iphase<10)
1220       addtagsize=1;
1221     else if(*iphase<100)
1222       addtagsize=2;
1223     else if(*iphase<1000)
1224       addtagsize=3;
1225 
1226     int tagsize2;
1227     tagsize2=*tagsize+addtagsize;
1228 
1229 //     char *fieldlabel = (char *)malloc((*tagsize+1)*sizeof(char));
1230 //     strncpy(fieldlabel, fieldtag, *tagsize);
1231 //     fieldlabel[*tagsize] = '\0';
1232 
1233     char *fieldlabel = (char *)malloc((tagsize2+1)*sizeof(char));
1234     strncpy(fieldlabel, fieldtag, *tagsize);
1235     fieldlabel[tagsize2] = '\0';
1236 
1237     char straddtagsize[10];
1238     sprintf(straddtagsize,"%d",*iphase);
1239 
1240     if(*iphase<10) {
1241       fieldlabel[tagsize2-1]=straddtagsize[0];
1242     }
1243     else if(*iphase<100) {
1244       fieldlabel[tagsize2-2]=straddtagsize[0];
1245       fieldlabel[tagsize2-1]=straddtagsize[1];
1246     }
1247     else if(*iphase<1000) {
1248       fieldlabel[tagsize2-3]=straddtagsize[0];
1249       fieldlabel[tagsize2-2]=straddtagsize[1];
1250       fieldlabel[tagsize2-1]=straddtagsize[2];
1251     }
1252 
1253     int irstou;
1254     int magic_number = 362436;
1255     int* mptr = &magic_number;
1256     double version=0.0;
1257     int isize, nitems;
1258     int iarray[10];
1259 
1260     char fmode[10];
1261     if(!strncmp(filemode,"w",1))
1262       strcpy(fmode,"write");
1263     else // default is append
1264       strcpy(fmode,"append");
1265 
1266     char datatype[10];
1267     if(!strncmp(arraytype,"i",1))
1268       strcpy(datatype,"int");
1269     else // default is double
1270       strcpy(datatype,"double");
1271 
1272 //
1273 // //     if(*iphase==1) //open the file but then keep it open for the remaining cycles
1274 //     openfile_(rfile, fmode, &irstou);
1275 //
1276 // //     printf("iphase: %d - pid: %d - irstou %d\n",*iphase,*pid,irstou);
1277 //
1278 //
1279 //     nitems = 3; // assuming field will write 3 items in iarray
1280 //     iarray[ 0 ] = (*nshg);
1281 //     iarray[ 1 ] = (*numvars);
1282 //     iarray[ 2 ] = (*stepno);
1283 //
1284 //     isize = (*nshg)*(*numvars);
1285 //     writeheader_( &irstou, fieldlabel, (void*)iarray, &nitems, &isize, datatype, phasta_iotype );
1286 //
1287 //     nitems = (*nshg)*(*numvars);
1288 //     writedatablock_( &irstou, fieldlabel, array, &nitems, datatype, phasta_iotype );
1289 //
1290 // //     if(*iphase==*nphasesincycle) //close the file after nphasesincycle
1291 //       closefile_( &irstou, fmode);
1292 //
1293 
1294     /////////////////////////////// Start of writing using new-lib ////////////////////////////
1295 
1296     int nfiles;
1297     int nfields;
1298     int numparts;
1299     int irank;
1300     int nprocs;
1301 //    unsigned long long timer_start;
1302 //    unsigned long long timer_end;
1303 //    double time_span;
1304 
1305     nfiles = outpar.nsynciofiles;
1306     nfields = outpar.nsynciofieldswriterestart;
1307     numparts = workfc.numpe;
1308     irank = *pid; //workfc.myrank;
1309     nprocs = workfc.numpe;
1310 
1311     int nppf = numparts/nfiles;
1312     int GPID;
1313 
1314     // Calculate number of parts each  proc deal with and where it start and end ...
1315     int nppp = numparts/nprocs;// nppp : Number of parts per proc ...
1316     int startpart = irank * nppp +1;// Part id from which I (myrank) start ...
1317     int endpart = startpart + nppp - 1;// Part id to which I (myrank) end ...
1318 
1319     //int descriptor;
1320     char filename[255],path[255],fieldtag_s[255];
1321     bzero((void*)filename,255);
1322     bzero((void*)fieldtag_s,255);
1323 
1324 //     char * namer;
1325 //     namer = strtok(fieldlabel," ");
1326 //     strncpy(fieldlabel, fieldtag, *tagsize);
1327 
1328     field_flag++;
1329     if(*pid==0) {
1330 //      printf("\n*****************************\n");
1331       printf("\n");
1332       printf("The %d/%d th field to be written is '%s'\n",field_flag,nfields,fieldlabel);
1333     }
1334 
1335     sprintf(filename,"restart-dat.%d.%d",*stepno,((int)(irank/(nprocs/nfiles))+1));
1336 
1337     int i;
1338     for ( i = 0; i < nppp; i++  ) {
1339         GPID = startpart + i;
1340 
1341         // Write solution field ...
1342         sprintf(fieldtag_s,"%s@%d",fieldlabel,GPID);
1343 
1344         //printf("This is %d and fieldtag_s is %s \n",myrank,fieldtag_s);
1345 
1346         isize = (*nshg)*(*numvars);
1347         nitems = 3;
1348         iarray[ 0 ] = (*nshg);
1349         iarray[ 1 ] = (*numvars);
1350         iarray[ 2 ] = (*stepno);
1351 
1352 //MR CHANGE
1353 //  Measure the time - Start the timer
1354 //        MPI_Barrier(MPI_COMM_WORLD);
1355 //        timer_start = rdtsc();
1356 //MR CHANGE END
1357 
1358         writeheader( &f_descriptor, fieldtag_s, (void*)iarray, &nitems, &isize, "double", phasta_iotype);
1359 
1360 //MR CHANGE
1361 //  Measure the time - End of timer
1362 //        MPI_Barrier(MPI_COMM_WORLD);
1363 //        timer_end = rdtsc();
1364 //        time_span=(double)((timer_end-timer_start)/cpu_speed);
1365 //        if (*pid==0) {
1366 //          printf("Time: header for '%s':    %f s\n",fieldtag_s,time_span);
1367 //        }
1368 //MR CHANGE END
1369 
1370         nitems = (*nshg)*(*numvars);
1371 
1372 //MR CHANGE
1373 //  Measure the time - Start the timer
1374 //        MPI_Barrier(MPI_COMM_WORLD);
1375 //        timer_start = rdtsc();
1376 //MR CHANGE END
1377 
1378         writedatablock( &f_descriptor, fieldtag_s, array, &isize, "double", phasta_iotype );
1379 
1380 //MR CHANGE
1381 //  Measure the time - End of timer
1382 //        MPI_Barrier(MPI_COMM_WORLD);
1383 //        timer_end = rdtsc();
1384 //        time_span=(double)((timer_end-timer_start)/cpu_speed);
1385 
1386 //        int isizemin,isizemax,isizetot;
1387 //        double sizemin,sizemax,sizeavg,sizetot,rate;
1388 
1389 //        MPI_Allreduce(&isize,&isizemin,1,MPI_INT,MPI_MIN,MPI_COMM_WORLD);
1390 //        MPI_Allreduce(&isize,&isizemax,1,MPI_INT,MPI_MAX,MPI_COMM_WORLD);
1391 //        MPI_Allreduce(&isize,&isizetot,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);
1392 
1393 //        sizemin=(double)(8.0*isizemin/1024.0/1024.0);
1394 //        sizemax=(double)(8.0*isizemax/1024.0/1024.0);
1395 //        sizetot=(double)(8.0*isizetot/1024.0/1024.0);
1396 //        sizeavg=sizetot/workfc.numpe;
1397 //        rate=sizetot/time_span;
1398 
1399 //        if (*pid==0) {
1400 //          printf("Time: block for '%s':    %f s\n",fieldtag_s,time_span);
1401 //          printf("Time: block:   Min= %f MB; Max= %f MB; Avg= %f MB; Tot= %f MB; Rate= %f MB/s; \n",sizemin,sizemax,sizeavg,sizetot,rate);
1402 //          printf("*****************************\n");
1403 //        }
1404 //MR CHANGE END
1405 
1406     }
1407 
1408 //     if (*pid==0) {
1409 //         printf("Field '%s' written in:     %f s\n",fieldtag,time_span);
1410 //         printf("Write field '%s' finished! \n",fieldtag_s);
1411 //     }
1412 
1413 //
1414 //     if (field_flag==nfields){
1415 //       closefile_(&f_descriptor, "write");
1416 //       finalizephmpiio_(&f_descriptor);
1417 //       if(*pid==0) {
1418 //         printf("Last field %d '%s' finished! \n",nfields, fieldtag_s);
1419 //         printf("\n*****************************\n");
1420 //       }
1421 //     }
1422 
1423     if (field_flag==nfields){
1424 
1425 //MR CHANGE
1426 //  Measure the time - Start the timer
1427 //      MPI_Barrier(MPI_COMM_WORLD);
1428 //      timer_start = rdtsc();
1429 //MR CHANGE END
1430 
1431       closefile(&f_descriptor, "write");
1432 
1433 //MR CHANGE
1434 //    Measure the time - End of timer
1435 //      MPI_Barrier(MPI_COMM_WORLD);
1436 //      timer_end = rdtsc();
1437 //      time_span=(double)((timer_end-timer_start)/cpu_speed);
1438 //     if (*pid==0) {
1439 //        printf("\n*****************************\n");
1440 //        printf("Time: 'closefile' is:    %f s\n",time_span);
1441 //      }
1442 //MR CHANGE END
1443 
1444 //MR CHANGE
1445 //  Measure the time - Start the timer
1446 //      MPI_Barrier(MPI_COMM_WORLD);
1447 //      timer_start = rdtsc();
1448 //MR CHANGE END
1449 
1450       finalizephmpiio(&f_descriptor);
1451 
1452 //MR CHANGE
1453 //    Measure the time - End of timer
1454 //      MPI_Barrier(MPI_COMM_WORLD);
1455 //      timer_end = rdtsc();
1456 //      time_span=(double)((timer_end-timer_start)/cpu_speed);
1457       if (*pid==0) {
1458 //        printf("Time: 'finalizephmpiio' is:    %f s\n",time_span);
1459 //        printf("Last field %d '%s' finished! \n",nfields, fieldtag);
1460         printf("\n");
1461 //        printf("*****************************\n");
1462       }
1463     }
1464 //MR CHANGE END
1465 
1466     ///////////////////////////////////////////////////////////////////////////////////////////
1467 
1468     free(fieldlabel);
1469 }
1470 
1471 
1472 void
1473 Write_d2wall(   int* pid,
1474                 int* numnp,
1475                 double* array1 ) {
1476 
1477 //    time_t timenow = time ( &timenow);
1478     int isize, nitems;
1479     int iarray[10];
1480 
1481 //    MPI_Barrier(MPI_COMM_WORLD);
1482 
1483     /////////////////////////////// Start of writing using new-lib ////////////////////////////
1484 
1485     int nfiles;
1486     int nfields;
1487     int numparts;
1488     int irank;
1489     int nprocs;
1490 
1491     //  First, count the number of fields to write and store the result in
1492     //countfieldstowriterestart();
1493 
1494     //  Retrieve and compute the parameters required for SyncIO
1495     nfiles = outpar.nsynciofiles;
1496     nfields = 1; //outpar.nsynciofieldswriterestart;  // Only the distance to the walls in d2wall
1497     numparts = workfc.numpe;
1498     irank = *pid; //workfc.myrank;
1499     nprocs = workfc.numpe;
1500     int nppf = numparts/nfiles;
1501     int GPID;
1502 
1503     // Calculate number of parts each proc deal with and where it start and end ...
1504     int nppp = numparts/nprocs;// nppp : Number of parts per proc ...
1505     int startpart = irank * nppp +1;// Part id from which I (myrank) start ...
1506     int endpart = startpart + nppp - 1;// Part id to which I (myrank) end ...
1507 
1508     int descriptor;
1509     char filename[255],path[255],fieldtag_s[255];
1510     bzero((void*)filename,255);
1511     bzero((void*)fieldtag_s,255);
1512 
1513     sprintf(filename,"d2wall.%d",((int)(irank/(nprocs/nfiles))+1));
1514 
1515     if (irank==0) {
1516       printf("Filename is %s \n",filename);
1517     }
1518 
1519     initphmpiio(&nfields, &nppf, &nfiles, &f_descriptor, "write");
1520 
1521     openfile(filename, "write", &f_descriptor);
1522 
1523     field_flag=0;
1524 
1525      int i;
1526      for ( i = 0; i < nppp; i++) { //This loop is useful only if several parts per processor
1527      // GPID : global part id, corresponds to rank ...
1528         // e.g : (in this example)
1529         // proc 0 : 1--4
1530         // proc 1 : 5--8 ...
1531         GPID = startpart + i;
1532 
1533         // Write solution field ...
1534         sprintf(fieldtag_s,"d2wall@%d",GPID);
1535 
1536         isize = (*numnp);
1537         nitems = 2;
1538         iarray[ 0 ] = (*numnp);
1539         iarray[ 1 ] = 1; //numVars = 1
1540 
1541         writeheader( &f_descriptor, fieldtag_s, (void*)iarray, &nitems, &isize, "double", phasta_iotype);
1542 
1543         //nitems = (*nshg)*(*numVars);
1544         //nitems = (*numnp);
1545 
1546         writedatablock( &f_descriptor, fieldtag_s, (void*)(array1), &isize, "double", phasta_iotype );
1547 
1548 
1549     }
1550     field_flag++;
1551 
1552     if (field_flag==nfields){
1553 
1554       closefile(&f_descriptor, "write");
1555 
1556       finalizephmpiio(&f_descriptor);
1557 
1558       if (irank==0) {
1559         printf("\n");
1560       }
1561     }
1562 }
1563 
1564