xref: /phasta/phSolver/common/new_interface.c (revision 20fccb6efa214427fa44d1eb8bf102ff5b377fb3)
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_Displ(  int* pid,
408               int* stepno,
409               int* nshg,
410               int* numVars,
411               double* array1 ) {
412   fprintf(stderr, "This function is dead...exiting\n");
413   exit(1);
414 }
415 
416 void
417 Write_Field(  int *pid,
418               char* filemode,
419               char* fieldtag,
420               int* tagsize,
421               void* array,
422               char* arraytype,
423               int* nshg,
424               int* numvars,
425               int* stepno) {
426     char *fieldlabel = (char *)malloc((*tagsize+1)*sizeof(char));
427     strncpy(fieldlabel, fieldtag, *tagsize);
428     fieldlabel[*tagsize] = '\0';
429 
430     int irstou;
431     double version=0.0;
432     int isize, nitems;
433     int iarray[10];
434 
435     char fmode[10];
436     if(!strncmp(filemode,"w",1))
437       strcpy(fmode,"write");
438     else // default is append
439       strcpy(fmode,"append");
440 
441     char datatype[10];
442     if(!strncmp(arraytype,"i",1))
443       strcpy(datatype,"int");
444     else // default is double
445       strcpy(datatype,"double");
446 
447     int nfiles;
448     int nfields;
449     int numparts;
450     int irank;
451     int nprocs;
452 
453     nfiles = outpar.nsynciofiles;
454     nfields = outpar.nsynciofieldswriterestart;
455     numparts = workfc.numpe;
456     irank = *pid; //workfc.myrank;
457     nprocs = workfc.numpe;
458 
459     // Calculate number of parts each  proc deal with and where it start and end ...
460     int nppp = numparts/nprocs;// nppp : Number of parts per proc ...
461     int startpart = irank * nppp +1;// Part id from which I (myrank) start ...
462     int endpart = startpart + nppp - 1;// Part id to which I (myrank) end ...
463 
464     strncpy(fieldlabel, fieldtag, *tagsize);
465 
466     field_flag++;
467     if(*pid==0) {
468       printf("\n");
469       printf("The %d/%d th field to be written is '%s'\n",field_flag,nfields,fieldlabel);
470     }
471 
472     int i;
473     for ( i = 0; i < nppp; i++  ) {
474         // Write solution field ...
475         isize = (*nshg)*(*numvars);
476         nitems = 3;
477         iarray[ 0 ] = (*nshg);
478         iarray[ 1 ] = (*numvars);
479         iarray[ 2 ] = (*stepno);
480 
481         phio_writeheader(f_descriptor, fieldlabel, (void*)iarray, &nitems,
482             &isize, datatype, phasta_iotype);
483         nitems = (*nshg)*(*numvars);
484         phio_writedatablock(f_descriptor, fieldlabel, array, &isize,
485             datatype, phasta_iotype );
486     }
487     if (field_flag==nfields){
488       phio_closefile(f_descriptor);
489       if (*pid==0) {
490         printf("Last field %d '%s' finished! \n",nfields, fieldtag);
491         printf("\n");
492       }
493     }
494     free(fieldlabel);
495 }
496 
497 void
498 Write_PhAvg2( int* pid,
499               char* filemode,
500               char* fieldtag,
501               int* tagsize,
502               int* iphase,
503               int* nphasesincycle,
504               void* array,
505               char* arraytype,
506               int* nshg,
507               int* numvars,
508               int* stepno) {
509     int addtagsize=0; // phase number is added to the name of the field
510     if(*iphase<10)
511       addtagsize=1;
512     else if(*iphase<100)
513       addtagsize=2;
514     else if(*iphase<1000)
515       addtagsize=3;
516 
517     int tagsize2;
518     tagsize2=*tagsize+addtagsize;
519 
520     char *fieldlabel = (char *)malloc((tagsize2+1)*sizeof(char));
521     strncpy(fieldlabel, fieldtag, *tagsize);
522     fieldlabel[tagsize2] = '\0';
523 
524     char straddtagsize[10];
525     sprintf(straddtagsize,"%d",*iphase);
526 
527     if(*iphase<10) {
528       fieldlabel[tagsize2-1]=straddtagsize[0];
529     }
530     else if(*iphase<100) {
531       fieldlabel[tagsize2-2]=straddtagsize[0];
532       fieldlabel[tagsize2-1]=straddtagsize[1];
533     }
534     else if(*iphase<1000) {
535       fieldlabel[tagsize2-3]=straddtagsize[0];
536       fieldlabel[tagsize2-2]=straddtagsize[1];
537       fieldlabel[tagsize2-1]=straddtagsize[2];
538     }
539 
540     int irstou;
541     double version=0.0;
542     int isize, nitems;
543     int iarray[10];
544 
545     char fmode[10];
546     if(!strncmp(filemode,"w",1))
547       strcpy(fmode,"write");
548     else // default is append
549       strcpy(fmode,"append");
550 
551     char datatype[10];
552     if(!strncmp(arraytype,"i",1))
553       strcpy(datatype,"int");
554     else // default is double
555       strcpy(datatype,"double");
556 
557     int nfiles;
558     int nfields;
559     int numparts;
560     int irank;
561     int nprocs;
562 
563     nfiles = outpar.nsynciofiles;
564     nfields = outpar.nsynciofieldswriterestart;
565     numparts = workfc.numpe;
566     irank = *pid; //workfc.myrank;
567     nprocs = workfc.numpe;
568 
569     // Calculate number of parts each  proc deal with and where it start and end ...
570     int nppp = numparts/nprocs;// nppp : Number of parts per proc ...
571     int startpart = irank * nppp +1;// Part id from which I (myrank) start ...
572     int endpart = startpart + nppp - 1;// Part id to which I (myrank) end ...
573 
574     field_flag++;
575     if(*pid==0) {
576       printf("\n");
577       printf("The %d/%d th field to be written is '%s'\n",field_flag,nfields,fieldlabel);
578     }
579 
580     int i;
581     for ( i = 0; i < nppp; i++  ) {
582         // Write solution field ...
583         isize = (*nshg)*(*numvars);
584         nitems = 3;
585         iarray[ 0 ] = (*nshg);
586         iarray[ 1 ] = (*numvars);
587         iarray[ 2 ] = (*stepno);
588         phio_writeheader(f_descriptor, fieldlabel, (void*)iarray, &nitems,
589             &isize, "double", phasta_iotype);
590         nitems = (*nshg)*(*numvars);
591         phio_writedatablock(f_descriptor, fieldlabel, array, &isize,
592             "double", phasta_iotype );
593     }
594     if (field_flag==nfields){
595       phio_closefile(f_descriptor);
596       if (*pid==0) {
597         printf("\n");
598       }
599     }
600     free(fieldlabel);
601 }
602