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