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