xref: /phasta/phSolver/common/new_interface.c (revision 3de22adcea938dcbea5f2b4ef2369694350775bd)
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 = 3; //magic number, 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 // Save every step  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     const char* magic_name = "byteorder magic number";
253     int magic_number = 362436;
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 ione = 1;
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     if(outpar.output_mode == -1 )
279       streamio_setup_write(&f_descriptor, streamio_get_r());
280     else if(outpar.output_mode == 0 )
281       posixio_setup(&f_descriptor, 'w');
282     else if(outpar.output_mode > 0 )
283       syncio_setup_write(nfiles, nfields, numparts/nfiles, &f_descriptor);
284     else
285       exit(EXIT_FAILURE);
286     phio_constructName(f_descriptor,"restart",filename);
287     phstr_appendInt(filename, *stepno);
288     phstr_appendStr(filename, ".");
289     phio_openfile(filename, f_descriptor);
290 
291     field_flag=0;
292 
293     // write the magic number
294     phio_writeheader(f_descriptor, magic_name, (void*)&magic_number, &ione,
295         &ione, "integer", phasta_iotype);
296     phio_writedatablock(f_descriptor, magic_name, (void*)&magic_number,
297         &ione, "integer", phasta_iotype );
298     field_flag++;
299 
300      int i;
301      for ( i = 0; i < nppp; i++) { //This loop is useful only if several parts per processor
302         // Write solution field ...
303         isize = (*nshg)*(*numVars);
304         nitems = 3;
305         iarray[ 0 ] = (*nshg);
306         iarray[ 1 ] = (*numVars);
307         iarray[ 2 ] = (*stepno);
308 
309         phio_writeheader(f_descriptor, "solution", (void*)iarray, &nitems,
310             &isize, "double", phasta_iotype);
311         nitems = (*nshg)*(*numVars);
312         phio_writedatablock(f_descriptor, "solution", (void*)(array1),
313             &isize, "double", phasta_iotype );
314     }
315     field_flag++;
316 
317     for ( i = 0; i < nppp; i++) {
318         // Write solution field ...
319         isize = (*nshg)*(*numVars);
320         nitems = 3;
321         iarray[ 0 ] = (*nshg);
322         iarray[ 1 ] = (*numVars);
323         iarray[ 2 ] = (*stepno);
324         phio_writeheader(f_descriptor, "time derivative of solution",
325             (void*)iarray, &nitems, &isize, "double", phasta_iotype);
326         nitems = (*nshg)*(*numVars);
327         phio_writedatablock(f_descriptor, "time derivative of solution",
328             (void*)(array2), &isize, "double", phasta_iotype );
329     }
330     field_flag++;
331 
332     if (field_flag==nfields){
333       phio_closefile(f_descriptor);
334       if (*pid==0) {
335         printf("\n");
336       }
337     }
338 }
339 
340 void
341 Write_Error(  int* pid,
342               int* stepno,
343               int* nshg,
344               int* numVars,
345               double* array1 ) {
346     char fname[255];
347     char rfile[60];
348     int irstou;
349     double version=0.0;
350     int isize, nitems;
351     int iarray[10];
352     int nfiles;
353     int nfields;
354     int numparts;
355     int irank;
356     int nprocs;
357 
358     nfiles = outpar.nsynciofiles;
359     nfields = outpar.nsynciofieldswriterestart;
360     numparts = workfc.numpe;
361     irank = *pid; //workfc.myrank;
362     nprocs = workfc.numpe;
363 
364     // Calculate number of parts each  proc deal with and where it start and end ...
365     int nppp = numparts/nprocs;// nppp : Number of parts per proc ...
366     int startpart = irank * nppp +1;// Part id from which I (myrank) start ...
367     int endpart = startpart + nppp - 1;// Part id to which I (myrank) end ...
368 
369     field_flag++;
370 
371     int i;
372     for ( i = 0; i < nppp; i++  ) {
373 
374         if(*pid==0) {
375           printf("\n");
376           printf("The %d/%d th field to be written is 'errors'\n",field_flag,nfields);
377         }
378 
379         isize = (*nshg)*(*numVars);
380         nitems = 3;
381         iarray[ 0 ] = (*nshg);
382         iarray[ 1 ] = (*numVars);
383         iarray[ 2 ] = (*stepno);
384 
385         phio_writeheader(f_descriptor, "errors", (void*)iarray, &nitems,
386             &isize, "double", phasta_iotype);
387 
388         phio_writedatablock(f_descriptor, "errors", (void*)array1, &isize,
389             "double", phasta_iotype );
390     }
391     if (field_flag==nfields){
392       phio_closefile(f_descriptor);
393       if (*pid==0) {
394         printf("Last field %d 'errors' finished! \n",nfields);
395         printf("\n");
396       }
397     }
398 }
399 
400 void
401 Write_Displ(  int* pid,
402               int* stepno,
403               int* nshg,
404               int* numVars,
405               double* array1 ) {
406   fprintf(stderr, "This function is dead...exiting\n");
407   exit(1);
408 }
409 
410 void
411 Write_Field(  int *pid,
412               char* filemode,
413               char* fieldtag,
414               int* tagsize,
415               void* array,
416               char* arraytype,
417               int* nshg,
418               int* numvars,
419               int* stepno) {
420     char *fieldlabel = (char *)malloc((*tagsize+1)*sizeof(char));
421     strncpy(fieldlabel, fieldtag, *tagsize);
422     fieldlabel[*tagsize] = '\0';
423 
424     int irstou;
425     double version=0.0;
426     int isize, nitems;
427     int iarray[10];
428 
429     char fmode[10];
430     if(!strncmp(filemode,"w",1))
431       strcpy(fmode,"write");
432     else // default is append
433       strcpy(fmode,"append");
434 
435     char datatype[10];
436     if(!strncmp(arraytype,"i",1))
437       strcpy(datatype,"int");
438     else // default is double
439       strcpy(datatype,"double");
440 
441     int nfiles;
442     int nfields;
443     int numparts;
444     int irank;
445     int nprocs;
446 
447     nfiles = outpar.nsynciofiles;
448     nfields = outpar.nsynciofieldswriterestart;
449     numparts = workfc.numpe;
450     irank = *pid; //workfc.myrank;
451     nprocs = workfc.numpe;
452 
453     // Calculate number of parts each  proc deal with and where it start and end ...
454     int nppp = numparts/nprocs;// nppp : Number of parts per proc ...
455     int startpart = irank * nppp +1;// Part id from which I (myrank) start ...
456     int endpart = startpart + nppp - 1;// Part id to which I (myrank) end ...
457 
458     strncpy(fieldlabel, fieldtag, *tagsize);
459 
460     field_flag++;
461     if(*pid==0) {
462       printf("\n");
463       printf("The %d/%d th field to be written is '%s'\n",field_flag,nfields,fieldlabel);
464     }
465 
466     int i;
467     for ( i = 0; i < nppp; i++  ) {
468         // Write solution field ...
469         isize = (*nshg)*(*numvars);
470         nitems = 3;
471         iarray[ 0 ] = (*nshg);
472         iarray[ 1 ] = (*numvars);
473         iarray[ 2 ] = (*stepno);
474 
475         phio_writeheader(f_descriptor, fieldlabel, (void*)iarray, &nitems,
476             &isize, datatype, phasta_iotype);
477         nitems = (*nshg)*(*numvars);
478         phio_writedatablock(f_descriptor, fieldlabel, array, &isize,
479             datatype, phasta_iotype );
480     }
481     if (field_flag==nfields){
482       phio_closefile(f_descriptor);
483       if (*pid==0) {
484         printf("Last field %d '%s' finished! \n",nfields, fieldtag);
485         printf("\n");
486       }
487     }
488     free(fieldlabel);
489 }
490 
491 void
492 Write_PhAvg2( int* pid,
493               char* filemode,
494               char* fieldtag,
495               int* tagsize,
496               int* iphase,
497               int* nphasesincycle,
498               void* array,
499               char* arraytype,
500               int* nshg,
501               int* numvars,
502               int* stepno) {
503     int addtagsize=0; // phase number is added to the name of the field
504     if(*iphase<10)
505       addtagsize=1;
506     else if(*iphase<100)
507       addtagsize=2;
508     else if(*iphase<1000)
509       addtagsize=3;
510 
511     int tagsize2;
512     tagsize2=*tagsize+addtagsize;
513 
514     char *fieldlabel = (char *)malloc((tagsize2+1)*sizeof(char));
515     strncpy(fieldlabel, fieldtag, *tagsize);
516     fieldlabel[tagsize2] = '\0';
517 
518     char straddtagsize[10];
519     sprintf(straddtagsize,"%d",*iphase);
520 
521     if(*iphase<10) {
522       fieldlabel[tagsize2-1]=straddtagsize[0];
523     }
524     else if(*iphase<100) {
525       fieldlabel[tagsize2-2]=straddtagsize[0];
526       fieldlabel[tagsize2-1]=straddtagsize[1];
527     }
528     else if(*iphase<1000) {
529       fieldlabel[tagsize2-3]=straddtagsize[0];
530       fieldlabel[tagsize2-2]=straddtagsize[1];
531       fieldlabel[tagsize2-1]=straddtagsize[2];
532     }
533 
534     int irstou;
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