xref: /phasta/phSolver/common/new_interface.c (revision 4c3261e2aac47bbd3436177e758cbd569440c5ce)
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 <stdlib.h>
7 #include <stdio.h>
8 #include <string.h>
9 #include <ctype.h>
10 #include <stdlib.h>
11 #include <time.h>
12 #include <math.h>
13 #include <assert.h>
14 #include "mpi.h"
15 #include "phastaIO.h"
16 #include "rdtsc.h"
17 #include <FCMangle.h>
18 #include "new_interface.h"
19 #include "phIO.h"
20 #include "syncio.h"
21 #include "posixio.h"
22 #include "streamio.h"
23 #include "common_c.h"
24 #include "tmrc.h"
25 
26 #ifdef intel
27 #include <winsock2.h>
28 #else
29 #include <unistd.h>
30 #include <strings.h>
31 #endif
32 
33 void igetMinMaxAvg(int *ivalue, double *stats, int *statRanks) {
34   int isThisRank;
35   double *value = (double*)malloc(sizeof(double));
36   *value = 1.0*(*ivalue);
37   rgetMinMaxAvg(value,stats,statRanks);
38   free(value);
39 }
40 
41 void rgetMinMaxAvg(double *value, double *stats, int *statRanks) {
42   int isThisRank;
43 
44   MPI_Allreduce(value,&stats[0],1,MPI_DOUBLE,MPI_MIN,MPI_COMM_WORLD);
45   isThisRank=workfc.numpe+1;
46   if(*value==stats[0])
47     isThisRank=workfc.myrank;
48   MPI_Allreduce(&isThisRank,&statRanks[0],1,MPI_INT,MPI_MIN,MPI_COMM_WORLD);
49 
50   MPI_Allreduce(value,&stats[1],1,MPI_DOUBLE,MPI_MAX,MPI_COMM_WORLD);
51   isThisRank=workfc.numpe+1;
52   if(*value==stats[1])
53     isThisRank=workfc.myrank;
54   MPI_Allreduce(&isThisRank,&statRanks[1],1,MPI_INT,MPI_MIN,MPI_COMM_WORLD);
55 
56   MPI_Allreduce(value,&stats[2],1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
57   stats[2] /= workfc.numpe;
58 
59   double sqValue = (*value)*(*value), sqValueAvg = 0.;
60   MPI_Allreduce(&sqValue,&sqValueAvg,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
61   sqValueAvg /= workfc.numpe;
62 
63   stats[3] = sqrt(sqValueAvg-stats[2]*stats[2]);
64 }
65 
66 void print_mesh_stats(void) {
67   int statRanks[2];
68   double iStats[4], rStats[4];
69 
70   igetMinMaxAvg(&conpar.nshg,iStats,statRanks);
71   if(workfc.myrank==workfc.master)
72     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]);
73   igetMinMaxAvg(&conpar.numel,iStats,statRanks);
74   if(workfc.myrank==workfc.master)
75     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]);
76   igetMinMaxAvg(&conpar.numelb,iStats,statRanks);
77   if(workfc.myrank==workfc.master)
78     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]);
79   igetMinMaxAvg(&conpar.nnz_tot,iStats,statRanks);
80   if(workfc.myrank==workfc.master) {
81     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]);
82     printf("\n");
83   }
84 }
85 
86 void print_mpi_stats(void) {
87   int statRanks[2];
88   double iStats[4], rStats[4];
89 
90 /* NS equations*/
91   igetMinMaxAvg(&mpistats.iISend,iStats,statRanks);
92   if(workfc.myrank==workfc.master)
93     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]);
94   igetMinMaxAvg(&mpistats.iIRecv,iStats,statRanks);
95   if(workfc.myrank==workfc.master)
96     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]);
97   igetMinMaxAvg(&mpistats.iWaitAll,iStats,statRanks);
98   if(workfc.myrank==workfc.master)
99     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]);
100   igetMinMaxAvg(&mpistats.iAllR,iStats,statRanks);
101   if(workfc.myrank==workfc.master)
102     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]);
103 
104   rgetMinMaxAvg(&mpistats.rISend,rStats,statRanks);
105   if(workfc.myrank==workfc.master)
106     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]);
107   rgetMinMaxAvg(&mpistats.rIRecv,rStats,statRanks);
108   if(workfc.myrank==workfc.master)
109     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]);
110   rgetMinMaxAvg(&mpistats.rWaitAll,rStats,statRanks);
111   if(workfc.myrank==workfc.master)
112     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]);
113   rgetMinMaxAvg(&mpistats.rCommu,rStats,statRanks);
114   if(workfc.myrank==workfc.master)
115     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]);
116   rgetMinMaxAvg(&mpistats.rAllR,rStats,statRanks);
117   if(workfc.myrank==workfc.master) {
118     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]);
119     printf("\n");
120   }
121 /* Scalars*/
122   igetMinMaxAvg(&mpistats.iISendScal,iStats,statRanks);
123   if(workfc.myrank==workfc.master)
124     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]);
125   igetMinMaxAvg(&mpistats.iIRecvScal,iStats,statRanks);
126   if(workfc.myrank==workfc.master)
127     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]);
128   igetMinMaxAvg(&mpistats.iWaitAllScal,iStats,statRanks);
129   if(workfc.myrank==workfc.master)
130     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]);
131   igetMinMaxAvg(&mpistats.iAllRScal,iStats,statRanks);
132   if(workfc.myrank==workfc.master)
133     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]);
134 
135   rgetMinMaxAvg(&mpistats.rISendScal,rStats,statRanks);
136   if(workfc.myrank==workfc.master)
137     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]);
138   rgetMinMaxAvg(&mpistats.rIRecvScal,rStats,statRanks);
139   if(workfc.myrank==workfc.master)
140     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]);
141   rgetMinMaxAvg(&mpistats.rWaitAllScal,rStats,statRanks);
142   if(workfc.myrank==workfc.master)
143     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]);
144   rgetMinMaxAvg(&mpistats.rCommuScal,rStats,statRanks);
145   if(workfc.myrank==workfc.master)
146     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]);
147   rgetMinMaxAvg(&mpistats.rAllRScal,rStats,statRanks);
148   if(workfc.myrank==workfc.master)
149     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]);
150 
151 
152 }
153 
154 void print_system_stats(double *tcorecp, double *tcorecpscal) {
155   int statRanks[2];
156   double iStats[4], rStats[4];
157   double syst_assembly, syst_solve;
158 
159 /* NS equations */
160   syst_assembly = tcorecp[0];
161   syst_solve = tcorecp[1];
162 
163   rgetMinMaxAvg(&syst_assembly,rStats,statRanks);
164   if(workfc.myrank==workfc.master)
165     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]);
166 
167   rgetMinMaxAvg(&syst_solve,rStats,statRanks);
168   if(workfc.myrank==workfc.master)
169     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]);
170 
171 /* Scalars */
172   syst_assembly = tcorecpscal[0];
173   syst_solve = tcorecpscal[1];
174 
175   rgetMinMaxAvg(&syst_assembly,rStats,statRanks);
176   if(workfc.myrank==workfc.master)
177     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]);
178 
179   rgetMinMaxAvg(&syst_solve,rStats,statRanks);
180   if(workfc.myrank==workfc.master) {
181     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]);
182     printf("\n");
183   }
184 }
185 
186 
187 
188 void countfieldstowriterestart()
189 {
190   int nfields = 3; /*magic number, solution, time derivatives*/
191 
192   if(outpar.ivort == 1){
193     nfields++; /*vorticity*/
194   }
195 
196   if(abs(turbvar.itwmod) != 1 && outpar.iowflux == 1) {
197     nfields++; /*instantaneous wss in bflux.f */
198   }
199 
200 /* Save every step  if(timdat.istep == inpdat.nstep[timdat.itseq-1]){ */ /*Last time step of the computation*/
201 
202     /*projection vectors and pressure projection vectors (call saveLesRestart in itrdrv)*/
203     if(matdat.matflg[0][0] ==-1) {
204         nfields = nfields +2;
205     }
206 
207     /*if Print Error Indicators = true (call write_error in itrdrv)*/
208     if(turbvar.ierrcalc == 1){
209       nfields++;
210     }
211 
212     /*if Print ybar = True (call write_field(myrank,'a','ybar',4,... in itrdrv)*/
213     if(outpar.ioybar == 1){
214       nfields++;  /*ybar*/
215 
216       /*phase average fields*/
217       if(outpar.nphasesincycle >0) {
218         nfields = nfields + outpar.nphasesincycle;
219       }
220 
221       if(abs(turbvar.itwmod) != 1 && outpar.iowflux == 1) {
222         nfields++; /*wssbar*/
223       }
224 
225     }
226 
227     if(turbvari.irans < 0) {
228       nfields++; /*dwal*/
229     }
230 
231 /*  }*/
232 
233   outpar.nsynciofieldswriterestart = nfields;
234 
235   if(workfc.myrank == 0) {
236     printf("Number of fields to write in restart files: %d\n", nfields);
237   }
238 }
239 
240 
241 void
242 Write_Restart(  int* pid,
243                 int* stepno,
244                 int* nshg,
245                 int* numVars,
246                 double* array1,
247                 double* array2 ) {
248 
249     char fname[255];
250     char rfile[60];
251     char existingfile[30], linkfile[30];
252     int irstou;
253     const char* magic_name = "byteorder magic number";
254     int magic_number = 362436;
255     double version=0.0;
256     int isize, nitems;
257     int iarray[10];
258     int nfiles;
259     int nfields;
260     int numparts;
261     int irank;
262     int nprocs;
263     int ione = 1;
264     double iotime = 0;
265 
266     /*  First, count the number of fields to write and store the result in*/
267     countfieldstowriterestart();
268 
269     /*  Retrieve and compute the parameters required for SyncIO*/
270     nfiles = outpar.nsynciofiles;
271     nfields = outpar.nsynciofieldswriterestart;
272     numparts = workfc.numpe;
273     irank = *pid; /*workfc.myrank;*/
274     nprocs = workfc.numpe;
275     int nppp = numparts/nprocs;   /* always 1 for PHASTA*/
276     int descriptor;
277     char filename[255];
278     bzero((void*)filename,255);
279 
280     iotime = TMRC();
281     if(outpar.output_mode == -1 )
282       streamio_setup_write(&f_descriptor, streamio_get_r());
283     else if(outpar.output_mode == 0 )
284       posixio_setup(&f_descriptor, 'w');
285     else if(outpar.output_mode > 0 )
286       syncio_setup_write(nfiles, nfields, numparts/nfiles, &f_descriptor);
287     else
288       exit(EXIT_FAILURE);
289     phio_constructName(f_descriptor,"restart",filename);
290     phstr_appendInt(filename, *stepno);
291     phstr_appendStr(filename, ".");
292     phio_openfile(filename, f_descriptor);
293 
294     field_flag=0;
295 
296     /* write the magic number*/
297     phio_writeheader(f_descriptor, magic_name, (void*)&magic_number, &ione,
298         &ione, "integer", phasta_iotype);
299     phio_writedatablock(f_descriptor, magic_name, (void*)&magic_number,
300         &ione, "integer", phasta_iotype );
301     field_flag++;
302 
303      int i;
304      for ( i = 0; i < nppp; i++) { /*This loop is useful only if several parts per processor*/
305         /* Write solution field ...*/
306         isize = (*nshg)*(*numVars);
307         nitems = 3;
308         iarray[ 0 ] = (*nshg);
309         iarray[ 1 ] = (*numVars);
310         iarray[ 2 ] = (*stepno);
311 
312         phio_writeheader(f_descriptor, "solution", (void*)iarray, &nitems,
313             &isize, "double", phasta_iotype);
314         nitems = (*nshg)*(*numVars);
315         phio_writedatablock(f_descriptor, "solution", (void*)(array1),
316             &isize, "double", phasta_iotype );
317     }
318     field_flag++;
319 
320     for ( i = 0; i < nppp; i++) {
321         /* Write solution field ...*/
322         isize = (*nshg)*(*numVars);
323         nitems = 3;
324         iarray[ 0 ] = (*nshg);
325         iarray[ 1 ] = (*numVars);
326         iarray[ 2 ] = (*stepno);
327         phio_writeheader(f_descriptor, "time derivative of solution",
328             (void*)iarray, &nitems, &isize, "double", phasta_iotype);
329         nitems = (*nshg)*(*numVars);
330         phio_writedatablock(f_descriptor, "time derivative of solution",
331             (void*)(array2), &isize, "double", phasta_iotype );
332     }
333     field_flag++;
334 
335     if (field_flag==nfields){
336       phio_closefile(f_descriptor);
337       if (*pid==0) {
338         printf("\n");
339       }
340     }
341     iotime = TMRC() - iotime;
342     if (workfc.master == workfc.myrank)
343       printf("time to write restart (seconds) %f\n",iotime);
344 }
345 
346 void
347 Write_Error(  int* pid,
348               int* stepno,
349               int* nshg,
350               int* numVars,
351               double* array1 ) {
352     char fname[255];
353     char rfile[60];
354     int irstou;
355     double version=0.0;
356     int isize, nitems;
357     int iarray[10];
358     int nfiles;
359     int nfields;
360     int numparts;
361     int irank;
362     int nprocs;
363 
364     nfiles = outpar.nsynciofiles;
365     nfields = outpar.nsynciofieldswriterestart;
366     numparts = workfc.numpe;
367     irank = *pid; /*workfc.myrank;*/
368     nprocs = workfc.numpe;
369 
370     /* Calculate number of parts each  proc deal with and where it start and end ...*/
371     int nppp = numparts/nprocs;/* nppp : Number of parts per proc ...*/
372     int startpart = irank * nppp +1;/* Part id from which I (myrank) start ...*/
373     int endpart = startpart + nppp - 1;/* Part id to which I (myrank) end ...*/
374 
375     field_flag++;
376 
377     int i;
378     for ( i = 0; i < nppp; i++  ) {
379 
380         if(*pid==0) {
381           printf("\n");
382           printf("The %d/%d th field to be written is 'errors'\n",field_flag,nfields);
383         }
384 
385         isize = (*nshg)*(*numVars);
386         nitems = 3;
387         iarray[ 0 ] = (*nshg);
388         iarray[ 1 ] = (*numVars);
389         iarray[ 2 ] = (*stepno);
390 
391         phio_writeheader(f_descriptor, "errors", (void*)iarray, &nitems,
392             &isize, "double", phasta_iotype);
393 
394         phio_writedatablock(f_descriptor, "errors", (void*)array1, &isize,
395             "double", phasta_iotype );
396     }
397     if (field_flag==nfields){
398       phio_closefile(f_descriptor);
399       if (*pid==0) {
400         printf("Last field %d 'errors' finished! \n",nfields);
401         printf("\n");
402       }
403     }
404 }
405 
406 void
407 Write_Field(  int *pid,
408               char* filemode,
409               char* fieldtag,
410               int* tagsize,
411               void* array,
412               char* arraytype,
413               int* nshg,
414               int* numvars,
415               int* stepno) {
416     char *fieldlabel = (char *)malloc((*tagsize+1)*sizeof(char));
417     strncpy(fieldlabel, fieldtag, *tagsize);
418     fieldlabel[*tagsize] = '\0';
419 
420     int irstou;
421     double version=0.0;
422     int isize, nitems;
423     int iarray[10];
424 
425     char fmode[10];
426     if(!strncmp(filemode,"w",1))
427       strcpy(fmode,"write");
428     else /* default is append*/
429       strcpy(fmode,"append");
430 
431     char datatype[10];
432     if(!strncmp(arraytype,"i",1))
433       strcpy(datatype,"int");
434     else /* default is double*/
435       strcpy(datatype,"double");
436 
437     int nfiles;
438     int nfields;
439     int numparts;
440     int irank;
441     int nprocs;
442 
443     nfiles = outpar.nsynciofiles;
444     nfields = outpar.nsynciofieldswriterestart;
445     numparts = workfc.numpe;
446     irank = *pid; /*workfc.myrank;*/
447     nprocs = workfc.numpe;
448 
449     /* Calculate number of parts each  proc deal with and where it start and end ...*/
450     int nppp = numparts/nprocs;/* nppp : Number of parts per proc ...*/
451     int startpart = irank * nppp +1;/* Part id from which I (myrank) start ...*/
452     int endpart = startpart + nppp - 1;/* Part id to which I (myrank) end ...*/
453 
454     strncpy(fieldlabel, fieldtag, *tagsize);
455 
456     field_flag++;
457     if(*pid==0) {
458       printf("\n");
459       printf("The %d/%d th field to be written is '%s'\n",field_flag,nfields,fieldlabel);
460     }
461 
462     int i;
463     for ( i = 0; i < nppp; i++  ) {
464         /* Write solution field ...*/
465         isize = (*nshg)*(*numvars);
466         nitems = 3;
467         iarray[ 0 ] = (*nshg);
468         iarray[ 1 ] = (*numvars);
469         iarray[ 2 ] = (*stepno);
470 
471         phio_writeheader(f_descriptor, fieldlabel, (void*)iarray, &nitems,
472             &isize, datatype, phasta_iotype);
473         nitems = (*nshg)*(*numvars);
474         phio_writedatablock(f_descriptor, fieldlabel, array, &isize,
475             datatype, phasta_iotype );
476     }
477     if (field_flag==nfields){
478       phio_closefile(f_descriptor);
479       if (*pid==0) {
480         printf("Last field %d '%s' finished! \n",nfields, fieldtag);
481         printf("\n");
482       }
483     }
484     free(fieldlabel);
485 }
486 
487 void
488 Write_PhAvg2( int* pid,
489               char* filemode,
490               char* fieldtag,
491               int* tagsize,
492               int* iphase,
493               int* nphasesincycle,
494               void* array,
495               char* arraytype,
496               int* nshg,
497               int* numvars,
498               int* stepno) {
499     int addtagsize=0; /* phase number is added to the name of the field*/
500     if(*iphase<10)
501       addtagsize=1;
502     else if(*iphase<100)
503       addtagsize=2;
504     else if(*iphase<1000)
505       addtagsize=3;
506 
507     int tagsize2;
508     tagsize2=*tagsize+addtagsize;
509 
510     char *fieldlabel = (char *)malloc((tagsize2+1)*sizeof(char));
511     strncpy(fieldlabel, fieldtag, *tagsize);
512     fieldlabel[tagsize2] = '\0';
513 
514     char straddtagsize[10];
515     sprintf(straddtagsize,"%d",*iphase);
516 
517     if(*iphase<10) {
518       fieldlabel[tagsize2-1]=straddtagsize[0];
519     }
520     else if(*iphase<100) {
521       fieldlabel[tagsize2-2]=straddtagsize[0];
522       fieldlabel[tagsize2-1]=straddtagsize[1];
523     }
524     else if(*iphase<1000) {
525       fieldlabel[tagsize2-3]=straddtagsize[0];
526       fieldlabel[tagsize2-2]=straddtagsize[1];
527       fieldlabel[tagsize2-1]=straddtagsize[2];
528     }
529 
530     int irstou;
531     double version=0.0;
532     int isize, nitems;
533     int iarray[10];
534 
535     char fmode[10];
536     if(!strncmp(filemode,"w",1))
537       strcpy(fmode,"write");
538     else /* default is append*/
539       strcpy(fmode,"append");
540 
541     char datatype[10];
542     if(!strncmp(arraytype,"i",1))
543       strcpy(datatype,"int");
544     else /* default is double*/
545       strcpy(datatype,"double");
546 
547     int nfiles;
548     int nfields;
549     int numparts;
550     int irank;
551     int nprocs;
552 
553     nfiles = outpar.nsynciofiles;
554     nfields = outpar.nsynciofieldswriterestart;
555     numparts = workfc.numpe;
556     irank = *pid; /*workfc.myrank;*/
557     nprocs = workfc.numpe;
558 
559     /* Calculate number of parts each  proc deal with and where it start and end ...*/
560     int nppp = numparts/nprocs;/* nppp : Number of parts per proc ...*/
561     int startpart = irank * nppp +1;/* Part id from which I (myrank) start ...*/
562     int endpart = startpart + nppp - 1;/* Part id to which I (myrank) end ...*/
563 
564     field_flag++;
565     if(*pid==0) {
566       printf("\n");
567       printf("The %d/%d th field to be written is '%s'\n",field_flag,nfields,fieldlabel);
568     }
569 
570     int i;
571     for ( i = 0; i < nppp; i++  ) {
572         /* Write solution field ...*/
573         isize = (*nshg)*(*numvars);
574         nitems = 3;
575         iarray[ 0 ] = (*nshg);
576         iarray[ 1 ] = (*numvars);
577         iarray[ 2 ] = (*stepno);
578         phio_writeheader(f_descriptor, fieldlabel, (void*)iarray, &nitems,
579             &isize, "double", phasta_iotype);
580         nitems = (*nshg)*(*numvars);
581         phio_writedatablock(f_descriptor, fieldlabel, array, &isize,
582             "double", phasta_iotype );
583     }
584     if (field_flag==nfields){
585       phio_closefile(f_descriptor);
586       if (*pid==0) {
587         printf("\n");
588       }
589     }
590     free(fieldlabel);
591 }
592