xref: /phasta/converterIO/converterO2N.cc (revision 16223cb9c3f88b34f2cb94151b5cf5ffc1aac5e2)
1 #include <stdio.h>
2 #include <iostream>
3 #include <string.h>
4 #include <stdlib.h>
5 //#define OMPI_SKIP_MPICXX 1 //Added in the CMakeList.txt file
6 #include <mpi.h>
7 #include <math.h>
8 #include "phastaIO.h"
9 #include <sys/stat.h>
10 #include <sys/types.h>
11 
12 enum {
13   DIR_FANOUT = 2048
14 };
15 
16 inline int
17 cscompare( const char teststring[],
18 	   const char targetstring[] )
19 {
20   char* s1 = const_cast<char*>(teststring);
21   char* s2 = const_cast<char*>(targetstring);
22 
23   while( *s1 == ' ') s1++;
24   while( *s2 == ' ') s2++;
25   while( ( *s1 )
26 	 && ( *s2 )
27 	 && ( *s2 != '?')
28 	 && ( tolower( *s1 )==tolower( *s2 ) ) ) {
29     s1++;
30     s2++;
31     while( *s1 == ' ') s1++;
32     while( *s2 == ' ') s2++;
33   }
34   if ( !( *s1 ) || ( *s1 == '?') ) return 1;
35   else return 0;
36 }
37 
38 inline int
39 computenitems(const int localpartid, const int fieldid, const int myrank, const char *fieldName, int ***para, const int intHeader, const int numVariables) {
40 // This routine computes the number of items in the data block based on
41 // - the name of the fields
42 // - the integers read from the header
43 
44   int nItems = -1;
45 
46   if (cscompare("nbc values",fieldName))
47     nItems = para[localpartid][fieldid][0] * (numVariables+1);
48   else if  (cscompare("nbc codes",fieldName))
49     nItems = para[localpartid][fieldid][0] * 2;
50   else if ( intHeader==1)
51     nItems = para[localpartid][fieldid][0];
52   else
53     nItems = para[localpartid][fieldid][0] * para[localpartid][fieldid][1];
54   return nItems;
55 }
56 
57 
58 int main(int argc, char *argv[]) {
59 
60   MPI_Init(&argc,&argv);
61 
62   int myrank, N_procs;
63   MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
64   MPI_Comm_size(MPI_COMM_WORLD, &N_procs);
65 
66   FILE * pFile,* AccessoryFileHandle;
67   char target[1024], pool[256], tempString[128],numpe[8],numstart[8];
68   char * temp, * token;
69   int fieldCompareMark;
70   bool readField;
71   int i, j, k, N_geombc_double, N_geombc_integer, N_restart_double;
72   int N_restart_integer, N_steps, N_parts, N_files;
73 
74   pFile = fopen("./IO.O2N.input","r");
75   if (pFile == NULL)
76     printf("Error openning\n");
77 
78   fgets( target, 1024, pFile );
79   token = strtok ( target, ";" );strcpy(pool,token);
80   temp = strtok ( pool, ":" );temp = strtok ( NULL, ":" );
81   N_geombc_double = atoi(temp);
82 
83   fgets( target, 1024, pFile );
84   token = strtok ( target, ";" );strcpy(pool,token);
85   temp = strtok ( pool, ":" );temp = strtok ( NULL, ":" );
86   N_geombc_integer = atoi(temp);
87 
88   fgets( target, 1024, pFile );
89   token = strtok ( target, ";" );strcpy(pool,token);
90   temp = strtok ( pool, ":" );temp = strtok ( NULL, ":" );
91   N_restart_double = atoi(temp);
92 
93   fgets( target, 1024, pFile );
94   token = strtok ( target, ";" );strcpy(pool,token);
95   temp = strtok ( pool, ":" );temp = strtok ( NULL, ":" );
96   N_restart_integer = atoi(temp);
97 
98   fgets( target, 1024, pFile );
99   token = strtok ( target, ";" );strcpy(pool,token);
100   temp = strtok ( pool, ":" );temp = strtok ( NULL, ":" );
101   strncpy(numstart,temp,1);
102   N_steps = atoi(temp);
103 
104   fgets( target, 1024, pFile );
105   token = strtok ( target, ";" );strcpy(pool,token);
106   temp = strtok ( pool, ":" );temp = strtok ( NULL, ":" );
107   strncpy(numpe,temp,1);
108   N_parts = atoi(temp);
109 
110   if(myrank==0){
111     printf("numpe is %d and start is %d\n",N_parts,N_steps);
112   }
113 
114   fgets( target, 1024, pFile );
115   token = strtok ( target, ";" );strcpy(pool,token);
116   temp = strtok ( pool, ":" );temp = strtok ( NULL, ":" );
117   N_files = atoi(temp);
118 
119   double ***Dfield; int ***Ifield;
120   int ***paraD, ***paraI, *expectD, *expectI;
121   char **fieldNameD, **fileTypeD, **dataTypeD, **headerTypeD;
122   char **fieldNameI, **fileTypeI, **dataTypeI, **headerTypeI;
123 
124   int WriteLockD[N_geombc_double];
125   int WriteLockI[N_geombc_integer];
126 
127   int nppp = N_parts/N_procs;
128   int startpart = myrank * nppp +1;
129   int endpart = startpart + nppp - 1;
130   char gfname[64], numTemp[128];
131   int iarray[10], igeom, isize,TempFileHandle;
132 
133   ////////////////////////////////////////
134   // Test if the user has given the right  parameters related to the number of parts, procs and SyncIO files
135   ////////////////////////////////////////
136 
137   if (N_parts%N_files!=0)
138     {
139       printf("Input error: number of parts should be a multiple of number of files!\n");
140       printf("Please modify the IO.O2N.input file!\n");
141       return 0;
142     }
143   if (N_procs%N_files!=0)
144     {
145       printf("Input error: number of procs should be a multiple of number of files!\n");
146       printf("Please modify the IO.O2N.input file!\n");
147       return 0;
148     }
149   if (N_parts%N_procs!=0)
150     {
151       printf("Input error: number of parts should be a multiple of number of procs!\n");
152       printf("Please modify the IO.O2N.input file!\n");
153       return 0;
154     }
155 
156   /////////////////////////////////////
157   // Create numpe.in and numstart.dat
158   ////////////////////////////////////////
159 
160   if (myrank==0)
161     {
162 
163 //MR CHANGE
164       bzero((void*)gfname,64);
165       sprintf(gfname,"./%d-procs_case-SyncIO-%d",N_parts,N_files);
166       if(0<mkdir(gfname,0777)) { printf("ERROR - Could not create procs_case-SyncIO directory\n"); return 1; }
167 //MR CHANGE END
168 
169       bzero((void*)gfname,64);
170       sprintf(gfname,"./%d-procs_case-SyncIO-%d/numstart.dat",N_parts,N_files);
171       AccessoryFileHandle=fopen(gfname,"w");
172       fprintf(AccessoryFileHandle,"%d",N_steps);
173       fclose(AccessoryFileHandle);
174 
175       bzero((void*)gfname,64);
176       sprintf(gfname,"./%d-procs_case-SyncIO-%d/numpe.in",N_parts,N_files);
177       AccessoryFileHandle=fopen(gfname,"w");
178       fprintf(AccessoryFileHandle,"%d",N_parts);
179       fclose(AccessoryFileHandle);
180     }
181 
182   ///////////////////////////////////////
183   Dfield = new double**[nppp];
184   Ifield = new int**[nppp];
185 
186   paraD = new int**[nppp];
187   paraI = new int**[nppp];
188   /////////////////////////////////////
189 
190   int interiorMark[nppp], boundaryMark[nppp], codesMark[nppp], valuesMark[nppp];
191   int numVariables[nppp];
192 //  int numBoundaryFields[nppp], numInteriorFields[nppp];
193 
194   if(myrank==0){
195     printf("Starting to read some headers in the restart.##.## and geombc.dat.## files\n");
196   }
197 
198   // The number of variables does not vary from one part to another. Do not read this info from every part.
199   // The ideal would be to ask only rank 0 to read only part 0 and broadcast the information to the other ranks. It saves a lot of time for very big meshes.
200   // Right now, every rank will read the number of variables from the solution field of one part only, which is already better than reading this info from every part.
201     int ithree = 3;
202     j = 0;
203 
204     // Test if file exist in the procs_case directory
205     bzero((void*)gfname,64);
206     sprintf(gfname,"./%d-procs_case/restart.%d.%d", N_parts, N_steps, startpart+j);
207     int subdir = -1;
208     bool existsubdir = false;
209     FILE *pfiletest;
210     pfiletest = fopen(gfname, "r");
211     if (pfiletest == NULL ) {
212       // Test if file exist in the procs_case/subdir directory
213       subdir = (startpart+j-1) / DIR_FANOUT;
214       bzero((void*)gfname,64);
215       sprintf(gfname,"./%d-procs_case/%d/restart.%d.%d",N_parts, subdir, N_steps, startpart+j);
216       pfiletest = fopen(gfname, "r");
217       if (pfiletest == NULL ) {
218         printf("[%d] File %s does not exit - abort\n",myrank, gfname);
219         abort();
220       }
221       else{
222         existsubdir = true;
223         fclose(pfiletest);
224       }
225     }
226     else {
227       fclose(pfiletest);
228     }
229 
230     // Debug
231     //printf("Rank: %d - subdir: %d - path: %s\n",myrank, subdir, gfname);
232 
233     openfile(gfname,"read",&TempFileHandle);
234     readheader( &TempFileHandle,
235 		   "solution",
236 		   (void*)iarray,
237 		   &ithree,
238 		   "double",
239 		   "binary" );
240     closefile(&TempFileHandle, "read");
241     for ( j = 0; j < nppp; j++ )
242       numVariables[j] = iarray[1]; //iarray[i] contains the number of variables from the header of the solution field
243 
244     MPI_Barrier(MPI_COMM_WORLD); //added by MR
245 
246 
247   /////////////////////////
248   for ( i = 0; i < nppp; i++ )
249     {
250       Dfield[i] = new double*[N_geombc_double]; // Space for the datablock of double format
251       paraD[i] = new int*[N_geombc_double];     // Integers in the header of each field of double format
252       Ifield[i] = new int*[N_geombc_integer];   // Space for the datablock of integer format
253       paraI[i] = new int*[N_geombc_integer];    // Integers in the header of each field of integer format
254 
255     }
256 
257   expectD = new int[N_geombc_double];          // Expected number of integers in the header for each field of double format
258   expectI = new int[N_geombc_integer];         // Expected number of integers in the header for each field of integer format
259 
260   fieldNameD = new char*[N_geombc_double];     // Name of the field in double format
261   fileTypeD = new char*[N_geombc_double];      // geombc or restart (useless if associated with geombc file but read)
262   dataTypeD = new char*[N_geombc_double];      // Integer or double (useless if associated with double data but read)
263   headerTypeD = new char*[N_geombc_double];    // block (means with data block) or header (just a header with no block)
264 
265   fieldNameI = new char*[N_geombc_integer];
266   fileTypeI = new char*[N_geombc_integer];
267   dataTypeI = new char*[N_geombc_integer];
268   headerTypeI = new char*[N_geombc_integer];
269 
270   /////////////////////////
271 
272   for ( i = 0; i < N_geombc_double; i++ )
273     {
274       WriteLockD[i]=0;
275 
276       fieldNameD[i] = new char[128];
277       fileTypeD[i] = new char[128];
278       dataTypeD[i] = new char[128];
279       headerTypeD[i] = new char[128];
280     }
281 
282   for ( i = 0; i < N_geombc_integer; i++ )
283     {
284       WriteLockI[i]=0; //This take value 1 if the field requested in IO.input is not found (typically for tpblocks)
285 
286       fieldNameI[i] = new char[128];
287       fileTypeI[i] = new char[128];
288       dataTypeI[i] = new char[128];
289       headerTypeI[i] = new char[128];
290     }
291 
292   ////////////////////////////////////////////////////////////////
293   // Reading IO.input
294   // temporary fix: in the new version the double and integer
295   //                can mix and match to avoid the order confusion
296   ////////////////////////////////////////////////////////////////
297 
298   char S1[128],S2[128],S3[128],S4[128],S5[128];
299   int double_counter=0,int_counter=0;
300 
301   for ( i = 0; i < N_geombc_double+N_geombc_integer; i++ )
302     {
303       fgets( target, 1024, pFile );
304       temp = strtok( target, ";" );
305       token = strtok( temp, "," );
306       strcpy( S1, token );
307       token = strtok ( NULL, "," );
308       strcpy( S2, token );
309       token = strtok ( NULL, "," );
310       strcpy( S3, token );
311       token = strtok ( NULL, "," );
312       strcpy( S4, token );
313       token = strtok ( NULL, "," );
314       strcpy( S5, token );
315 
316       //if (cscompare(S3,"double"))
317       if (cscompare("double",S3))
318 	{
319 	  strcpy( fileTypeD[double_counter], S1 );
320 	  strcpy( fieldNameD[double_counter], S2 );
321 	  strcpy( dataTypeD[double_counter], S3 );
322 	  strcpy( headerTypeD[double_counter], S4 );
323 	  strcpy( numTemp, S5 );
324 	  expectD[double_counter] = atoi (numTemp);
325 	  double_counter++;
326 	}
327 
328       //if (cscompare(S3,"integer"))
329       if (cscompare("integer",S3))
330 	{
331 	  strcpy( fileTypeI[int_counter], S1 );
332 	  strcpy( fieldNameI[int_counter], S2 );
333 	  strcpy( dataTypeI[int_counter], S3 );
334 	  strcpy( headerTypeI[int_counter], S4 );
335 	  strcpy( numTemp, S5 );
336 	  expectI[int_counter] = atoi (numTemp);
337 	  int_counter++;
338 	}
339     }
340 
341   //////////////////////////////////////////////////////////////
342 
343   //for ( i = 0; i < N_geombc_double; i++) {
344     //printf("%d %s %s %s %s %d\n", myrank, fileTypeD[i], fieldNameD[i], dataTypeD[i], headerTypeD[i], expectD[i]);
345   //}
346 
347 
348 //  printf("rank %d is waiting\n",myrank);
349   MPI_Barrier(MPI_COMM_WORLD); //already there
350 
351   if(myrank==0){
352     printf("Starting to read some blocks (doubles) in the geombc.dat.## files\n");
353   }
354 
355   for ( i = 0; i < nppp; i++ )
356     {
357       if(existsubdir) {
358         subdir = (startpart+i-1) / DIR_FANOUT;
359         sprintf(gfname,"./%d-procs_case/%d/geombc.dat.%d",N_parts, subdir, startpart+i);
360       }
361       else
362         sprintf(gfname,"./%d-procs_case/geombc.dat.%d",N_parts,startpart+i);
363 
364       openfile(gfname,"read",&igeom);
365 
366 //      MPI_Barrier(MPI_COMM_WORLD);
367 
368       for ( j = 0; j < N_geombc_double; j++ )
369 	{
370 
371 //	  printf("rank: %d - read double: %s\n",myrank,fieldNameD[j]);
372 	  paraD[i][j] = new int[expectD[j]];
373 
374 	  for ( k = 0; k < 10; k++ )
375 	    iarray[k]=0;
376 
377 	  iarray[0]=-1;
378 	  WriteLockD[j]=0;
379 	  readheader( &igeom,
380 		       fieldNameD[j],
381 		       (void*)iarray,
382 		       &expectD[j],
383 		       "double",
384 		       "binary" );
385 	  if ( iarray[0]==-1 )  // The field requested in IO.O2N.input has not been found (should be a tpblocks)
386 	      WriteLockD[j]=1;
387 
388 //          printf("rank: %d - part: %d - field: %s - iarray: %d\n",myrank,startpart+i,fieldNameD[j],iarray[0]);
389 
390 //	  MPI_Barrier(MPI_COMM_WORLD);
391 
392 	  // get the parameter list in data header ...
393 	  // different fields have different ways to get this list ...
394 	  for ( k = 0; k < expectD[j]; k++ )
395 	    paraD[i][j][k] = iarray[k];
396 
397 	  if ( WriteLockD[j]==1) // Put the value of the expected integers in the header to 0 when field not present
398 	    for ( k = 0; k < expectD[j]; k++ )
399 	      paraD[i][j][k] = 0;
400 
401 /*          int iproc;
402           for(iproc=0; iproc<N_procs; iproc++){
403             MPI_Barrier(MPI_COMM_WORLD);
404             if(iproc == myrank){
405               printf(" iproc: %d ", iproc);
406               printf("part: %d myrank: %d field: %s header: ",startpart+i,myrank,fieldNameD[j]);
407               for ( k = 0; k < expectD[j]; k++ )
408                 printf(" %d ",iarray[k]);
409               printf("\n");
410             }
411             usleep(100);
412             MPI_Barrier(MPI_COMM_WORLD);
413           }*/
414 
415 
416 	  if ( cscompare("block",headerTypeD[j]) )
417 	    {
418 /*	      if ( expectD[j]==1)
419 		isize = paraD[i][j][0];
420 	      else
421 		isize = paraD[i][j][0] * paraD[i][j][1];
422 
423 	      if (cscompare("nbc values",fieldNameD[j]))
424 		isize = paraD[i][j][0] * (numVariables[i]+1);
425 */
426 
427 //              int test;
428 //              test = computenitems(i,j,myrank,fieldNameD[j],paraD,expectD[j],numVariables[i]);
429 //              printf("irank: %d fieldname: %s ParaD: %d\n",myrank,fieldNameD[j], paraD[0][0][0], numVariables[i]);
430 //              if(test != isize)
431 //                printf("PROBLEM fieldname: %s part: %d isize: %d test: %d\n",fieldNameD[j],startpart+i,isize,test);
432 //              else
433 //                printf("fieldname: %s part: %d isize: %d test: %d\n",fieldNameD[j],startpart+i,isize,test);
434 
435               isize = computenitems(i,j,myrank,fieldNameD[j],paraD,expectD[j],numVariables[i]);
436               //printf("fieldname: %s part: %d isize: %d\n",fieldNameD[j],startpart+i,isize);
437 
438 	      Dfield[i][j] = new double[isize];
439 	      readdatablock( &igeom,
440 			      fieldNameD[j],
441 			      (void*)Dfield[i][j],
442 			      &isize,
443 			      "double",
444 			      "binary" );
445 	    }
446 	}
447 //      MPI_Barrier(MPI_COMM_WORLD);
448       closefile(&igeom, "read");
449     }
450 
451   MPI_Barrier(MPI_COMM_WORLD);
452   if(myrank==0){
453     printf("Starting to read some blocks (integers) in the geombc.dat.## files\n");
454   }
455 
456   // Count the number of interior and boundary tpblocks for 2 new headers named
457   // 'total number of different interior tpblocks' and
458   // 'total number of different boundary tpblocks'
459   int interiorCounter, boundaryCounter;
460   interiorCounter=0;
461   boundaryCounter=0;
462   for ( j = 0; j < N_geombc_integer; j++ )
463     {
464        if (cscompare("connectivity interior",fieldNameI[j]))
465 	{
466                   //printf("part: %d, fieldNameI[j]: %s\n",GPID,fieldNameI[j]);
467 		  interiorCounter++;
468         }
469         else if (cscompare("connectivity boundary",fieldNameI[j]))
470 	{
471                   //printf("part: %d, fieldNameI[j]: %s\n",GPID,fieldNameI[j]);
472 		  boundaryCounter++;
473         }
474     }
475 
476   MPI_Barrier(MPI_COMM_WORLD);
477   if(myrank==0){
478     printf("There are %d total connectivity interior and %d total connectivity boundary\n", interiorCounter, boundaryCounter);
479   }
480 
481 
482 
483   // Now, start to read the integer fields
484   for ( i = 0; i < nppp; i++ )
485     {
486       if(existsubdir) {
487         subdir = (startpart+i-1) / DIR_FANOUT;
488         sprintf(gfname,"./%d-procs_case/%d/geombc.dat.%d",N_parts, subdir, startpart+i);
489       }
490       else
491         sprintf(gfname,"./%d-procs_case/geombc.dat.%d", N_parts, startpart+i);
492 
493       openfile(gfname,"read",&igeom);
494 
495 //      MPI_Barrier(MPI_COMM_WORLD);
496 
497 //      printf("gfname is %s and nppp is %d myrank %d\n",gfname,i,myrank);
498 
499       for ( j = 0; j < N_geombc_integer; j++ )
500 	{
501 
502 //	  printf("Writing integer ... %s\n",fieldNameI[j]);
503 	  paraI[i][j] = new int[expectI[j]];
504 
505 	  for ( k = 0; k < 10; k++ )
506 	    iarray[k]=0;
507 
508 	  //	  printf("myrank %d and i %d j %d numBou is %d\n",myrank,i,j,numBoundaryFields[i]);
509 
510 	  WriteLockI[j]=0;
511 	  iarray[0]=-1;
512 
513           if ( cscompare("total number of interior tpblocks",fieldNameI[j] ) )
514           {
515              iarray[0] = interiorCounter; //New header that does not exist in the posix file
516           }
517           else if ( cscompare("total number of boundary tpblocks",fieldNameI[j] ) )
518           {
519              iarray[0] = boundaryCounter; //New header that does not exist in the posix file
520           }
521           else
522           {
523 	     readheader( &igeom,
524 		          fieldNameI[j],
525 		          (void*)iarray,
526 		          &expectI[j],
527 		          "integer",
528 		          "binary" );
529 	     if ( iarray[0]==-1)
530 	       WriteLockI[j]=1; // The field was not found in the posix geombc file
531           }
532 
533 	  //MPI_Barrier(MPI_COMM_WORLD);
534 
535 /*          int iproc;
536           for(iproc=0; iproc<N_procs; iproc++){
537             MPI_Barrier(MPI_COMM_WORLD);
538             if(iproc == myrank){
539               printf(" iproc: %d ", iproc);
540               printf("part: %d myrank: %d field: %s header: ",startpart+i,myrank,fieldNameI[j]);
541               for ( k = 0; k < expectI[j]; k++ )
542                 printf(" %d ",iarray[k]);
543               printf("\n");
544             }
545             usleep(100);
546             MPI_Barrier(MPI_COMM_WORLD);
547           }*/
548 
549 	  for ( k = 0; k < expectI[j]; k++ )
550 	    paraI[i][j][k] = iarray[k];
551 
552 	  if ( WriteLockI[j]==1) //The field is not present but SyncIO needs it to read collectively. Put 0.
553 	    for ( k = 0; k < expectI[j]; k++ )
554 	      paraI[i][j][k] = 0;
555 
556 	  if ( cscompare("block",headerTypeI[j]) )
557 	    {
558 /*	      if ( expectI[j]==1)
559 		isize = paraI[i][j][0];
560 	      else
561 		isize = paraI[i][j][0] * paraI[i][j][1];
562 
563 	      if (cscompare("nbc codes",fieldNameI[j]))
564 		isize = paraI[i][j][0] * 2;
565 */
566 //              int test;
567 //              test = computenitems(i,j,myrank,fieldNameI[j],paraI,expectI[j],numVariables[i]);
568 //              printf("irank: %d fieldname: %s ParaI: %d\n",myrank,fieldNameI[j], parapI[0][0][0], numVariables[i]);
569 //              if(test != isize)
570 //                printf("PROBLEM fieldname: %s part: %d isize: %d test: %d\n",fieldNameI[j],startpart+i,isize,test);
571 //              else
572 //                printf("fieldname: %s part: %d isize: %d test: %d\n",fieldNameI[j],startpart+i,isize,test);
573 
574               isize = computenitems(i,j,myrank,fieldNameI[j],paraI,expectI[j],numVariables[i]);
575               //printf("fieldname: %s part: %d isize: %d\n",fieldNameI[j],startpart+i,isize);
576 
577 	      Ifield[i][j] = new int[isize];
578 	      readdatablock( &igeom,
579 			      fieldNameI[j],
580 			      (void*)Ifield[i][j],
581 			      &isize,
582 			      "integer",
583 			      "binary" );
584 	    }
585 	}
586 //      MPI_Barrier(MPI_COMM_WORLD);
587       closefile(&igeom, "read");
588     }
589 
590   MPI_Barrier(MPI_COMM_WORLD); //added by MR
591 
592   ///////////////////// Writing ///////////////////////////////
593 
594   int nppf = N_parts/N_files;
595   int N_geombc = N_geombc_double + N_geombc_integer;
596   int writeHandle, GPID;
597   char fname[255],fieldtag[255];
598 
599   bzero((void*)fname,255);
600 //MR CHANGE
601 //  sprintf(fname,"./%d-procs_case-SyncIO",N_parts);
602 //  if(0<mkdir(fname,0777)) { printf("ERROR - Could not create procs_case-SyncIO directory\n"); return 1; }
603 //MR CHANGE END
604   sprintf(fname,"./%d-procs_case-SyncIO-%d/geombc-dat.%d",N_parts,N_files,((int)(myrank/(N_procs/N_files))+1));
605 
606   MPI_Barrier(MPI_COMM_WORLD);
607   if(myrank==0){
608     printf("Starting to write some blocks (doubles) in the geombc.dat-## files\n");
609   }
610 
611   initphmpiio(&N_geombc, &nppf, &N_files,&writeHandle, "write");
612 //  initphmpiio(&N_geombc, &nppf, &N_files,&writeHandle);
613   openfile(fname, "write", &writeHandle);
614 
615   for (  i = 0; i < nppp; i++  )
616     {
617       valuesMark[i]=0;
618     }
619 
620   for ( j = 0; j < N_geombc_double; j++ )
621     {
622       for (  i = 0; i < nppp; i++  )
623 	{
624 
625 	  //if ( WriteLockD[i] == 0 )
626 	    {
627 	      GPID = startpart + i;
628 	      bzero((void*)fieldtag,255);
629 
630 	      fieldCompareMark=0;
631 	      if (cscompare("nbc values",fieldNameD[j]))
632 		{
633 		  fieldCompareMark = 1;
634 		  valuesMark[i]++;
635 		  bzero((void*)fieldNameD[j],128);
636 		  sprintf(fieldNameD[j],"nbc values%d",valuesMark[i]);
637 
638 //		  if ( valuesMark[i]>numBoundaryFields[i] )
639 //		   for ( k = 0; k < expectD[j]; k++ )
640 //		      paraD[i][j][k] = 0;
641 		}
642 
643 	      sprintf(fieldtag,"%s@%d",fieldNameD[j],GPID);
644 
645 
646 /*	      if ( expectD[j]==1 )
647 		isize = paraD[i][j][0];
648 	      else
649 		isize = paraD[i][j][0] * paraD[i][j][1];
650 
651 	      //specially designed for nbc values fields
652 	      //check the size in presolver codes
653 	      //Yeah, you have to open restart to get the size
654 	      if ( fieldCompareMark==1 )
655 		isize = paraD[i][j][0] * (numVariables[i]+1);
656 */
657 	      if ( cscompare("header",headerTypeD[j]) )
658 		isize = 0;
659               else // block
660                 isize = computenitems(i,j,myrank,fieldNameD[j],paraD,expectD[j],numVariables[i]);
661 
662 	      for ( k = 0; k < expectD[j]; k++ )
663 		iarray[k] = paraD[i][j][k];
664 
665               //printf("write fieldname: %s part: %d isize: %d iarray: %d\n",fieldNameD[j],startpart+i,isize, iarray[0]);
666 
667 	      writeheader( &writeHandle,
668 			    fieldtag,
669 			    (void*)iarray,
670 			    &expectD[j],
671 			    &isize,
672 			    "double",
673 			    "binary");
674 	      writedatablock( &writeHandle,
675 			       fieldtag,
676 			       (void*)Dfield[i][j],
677 			       &isize,
678 			       "double",
679 			       "binary");
680 
681 	      if ( cscompare("block",headerTypeD[j]) )
682 		delete [] Dfield[i][j];
683 	    }
684 	  delete [] paraD[i][j];
685 	}
686     }
687 
688   for (  i = 0; i < nppp; i++  )
689     {
690       interiorMark[i]=0;
691       boundaryMark[i]=0;
692       codesMark[i]=0;
693     }
694 
695   MPI_Barrier(MPI_COMM_WORLD);
696   if(myrank==0){
697     printf("Starting to write some blocks (integers) in the geombc.dat-## files\n");
698   }
699 
700   // Now the other fields listed in IO.O2N.input
701   for ( j = 0; j < N_geombc_integer; j++ )
702     {
703       for (  i = 0; i < nppp; i++  )
704 	{
705 	  //if ( WriteLockI[i] == 0 )
706 	    {
707 	      GPID = startpart + i;
708 	      bzero((void*)fieldtag,255);
709 
710 	      if (cscompare("connectivity interior",fieldNameI[j]))
711 		{
712 		  interiorMark[i]++;
713 		  bzero((void*)fieldNameI[j],128);
714 		  sprintf(fieldNameI[j],"connectivity interior%d",interiorMark[i]);
715 
716 //		  if ( interiorMark[i]>numInteriorFields[i] )
717 //		    for ( k = 0; k < expectI[j]; k++ )
718 //		      paraI[i][j][k] = 0;
719 
720 		}
721 
722 	      if (cscompare("connectivity boundary",fieldNameI[j]))
723 		{
724 		  boundaryMark[i]++;
725 		  bzero((void*)fieldNameI[j],128);
726 		  sprintf(fieldNameI[j],"connectivity boundary%d",boundaryMark[i]);
727 
728 //		  if ( boundaryMark[i]>numBoundaryFields[i] )
729 //		    for ( k = 0; k < expectI[j]; k++ )
730 //		      paraI[i][j][k] = 0;
731 		}
732 
733 	      fieldCompareMark=0;
734 	      if (cscompare("nbc codes",fieldNameI[j]))
735 		{
736 		  fieldCompareMark=1;
737 		  codesMark[i]++;
738 		  bzero((void*)fieldNameI[j],128);
739 		  sprintf(fieldNameI[j],"nbc codes%d",codesMark[i]);
740 
741 //		  if ( codesMark[i]>numBoundaryFields[i] )
742 //		    for ( k = 0; k < expectI[j]; k++ )
743 //		      paraI[i][j][k] = 0;
744 		}
745 
746 	      sprintf(fieldtag,"%s@%d",fieldNameI[j],GPID);
747 
748 //	      if ( expectI[j]==1)
749 //		isize = paraI[i][j][0];
750 //	      else
751 //		isize = paraI[i][j][0] * paraI[i][j][1];
752 
753 //MR CHANGE
754 //              printf("rank,j,i,isize %d %d %d %d\n",myrank,j,i,isize);
755 
756 
757 	      //specially designed for nbc codes fields
758 	      //check the size in presolver codes
759 //	      if (fieldCompareMark==1)
760 //		isize = paraI[i][j][0] * 2;
761 
762 	      if ( cscompare("header",headerTypeI[j]) )
763 		isize = 0;
764               else
765                 isize = computenitems(i,j,myrank,fieldNameI[j],paraI,expectI[j],numVariables[i]);
766 
767 	      for ( k = 0; k < expectI[j]; k++ )
768 		iarray[k] = paraI[i][j][k];
769 
770 //              printf("write fieldname: %s part: %d isize: %d iarray: %d\n",fieldNameI[j],startpart+i,isize, iarray[0]);
771 
772 	      writeheader( &writeHandle,
773 			    fieldtag,
774 			    (void*)iarray,
775 			    &expectI[j],
776 			    &isize,
777 			    "integer",
778 			    "binary");
779 	      writedatablock( &writeHandle,
780 			       fieldtag,
781 			       (void*)Ifield[i][j],
782 			       &isize,
783 			       "integer",
784 			       "binary" );
785 
786 	      if ( cscompare("block",headerTypeI[j]) ){
787 //                printf("rank %d - deleting Ifield %d %d\n",myrank,i,j);
788 		delete [] Ifield[i][j];
789 //                printf("rank %d - Ifield deleted %d %d\n",myrank,i,j);
790               }
791 	    }
792 //          printf("rank %d - deleting paraI %d %d\n",myrank,i,j);
793 	  delete [] paraI[i][j];
794 //          printf("rank %d - paraI deleted %d %d\n",myrank,i,j);
795 	}
796 
797     }
798 
799   MPI_Barrier(MPI_COMM_WORLD);
800   if(myrank==0){
801     printf("Closing geombc-dat.##.## files\n");
802   }
803 
804   closefile(&writeHandle, "write");
805   finalizephmpiio(&writeHandle);
806 
807   if(myrank==0){
808     printf("Free memory related to geombc-dat.##.## files\n");
809   }
810 
811   for ( i = 0; i < nppp; i++ )
812     {
813       delete [] Dfield[i];
814       delete [] paraD[i];
815 
816       delete [] Ifield[i];
817       delete [] paraI[i];
818 
819      }
820 
821   for ( i = 0; i < N_geombc_double; i++ )
822     {
823       delete [] fieldNameD[i];
824       delete [] fileTypeD[i];
825       delete [] dataTypeD[i];
826       delete [] headerTypeD[i];
827     }
828 
829   for ( i = 0; i < N_geombc_integer; i++ )
830     {
831 
832       delete [] fieldNameI[i];
833       delete [] fileTypeI[i];
834       delete [] dataTypeI[i];
835       delete [] headerTypeI[i];
836     }
837 
838   delete [] Dfield;
839   delete [] Ifield;
840 
841   delete [] paraD;
842   delete [] paraI;
843 
844   delete [] expectD;
845   delete [] expectI;
846 
847   delete [] fieldNameD;
848   delete [] fileTypeD;
849   delete [] dataTypeD;
850   delete [] headerTypeD;
851 
852   delete [] fieldNameI;
853   delete [] fileTypeI;
854   delete [] dataTypeI;
855   delete [] headerTypeI;
856 
857   MPI_Barrier(MPI_COMM_WORLD);
858   if(myrank==0){
859     printf("Done with geombc-dat.##.## files\n");
860   }
861 
862   /////////////////////// restart data ////////////////////////////
863 
864   int irestart;
865 
866   Dfield = new double**[N_restart_double];
867   Ifield = new int**[N_restart_integer];
868 
869   paraD = new int**[N_restart_double];
870   paraI = new int**[N_restart_integer];
871 
872   expectD = new int[N_restart_double];
873   expectI = new int[N_restart_integer];
874 
875   fieldNameD = new char*[N_restart_double];
876   fileTypeD = new char*[N_restart_double];
877   dataTypeD = new char*[N_restart_double];
878   headerTypeD = new char*[N_restart_double];
879 
880   fieldNameI = new char*[N_restart_integer];
881   fileTypeI = new char*[N_restart_integer];
882   dataTypeI = new char*[N_restart_integer];
883   headerTypeI = new char*[N_restart_integer];
884 
885   if (N_restart_double>0)
886     for ( i = 0; i < N_restart_double; i++ )
887       {
888 	WriteLockD[i]=0;
889 	Dfield[i] = new double*[nppp];
890 
891 	paraD[i] = new int*[nppp];
892 
893 	fieldNameD[i] = new char[128];
894 	fileTypeD[i] = new char[128];
895 	dataTypeD[i] = new char[128];
896 	headerTypeD[i] = new char[128];
897       }
898 
899   if (N_restart_integer>0)
900     for ( i = 0; i < N_restart_integer; i++ )
901       {
902 	WriteLockI[i]=0;
903 	Ifield[i] = new int*[nppp];
904 
905 	paraI[i] = new int*[nppp];
906 
907 	fieldNameI[i] = new char[128];
908 	fileTypeI[i] = new char[128];
909 	dataTypeI[i] = new char[128];
910 	headerTypeI[i] = new char[128];
911       }
912 
913   ////////////////////////////////////////////////////////////////
914   // temporary fix: in the new version the double and integer
915   //                can mix and match to avoid the order confusion
916   ////////////////////////////////////////////////////////////////
917 
918   double_counter=0,int_counter=0;
919 
920   for ( i = 0; i < N_restart_double+N_restart_integer; i++ )
921     {
922       fgets( target, 1024, pFile );
923       temp = strtok( target, ";" );
924       token = strtok( temp, "," );
925       strcpy( S1, token );
926       token = strtok ( NULL, "," );
927       strcpy( S2, token );
928       token = strtok ( NULL, "," );
929       strcpy( S3, token );
930       token = strtok ( NULL, "," );
931       strcpy( S4, token );
932       token = strtok ( NULL, "," );
933       strcpy( S5, token );
934 
935       if (cscompare(S3,"double"))
936 	{
937 	  strcpy( fileTypeD[double_counter], S1 );
938 	  strcpy( fieldNameD[double_counter], S2 );
939 	  strcpy( dataTypeD[double_counter], S3 );
940 	  strcpy( headerTypeD[double_counter], S4 );
941 	  strcpy( numTemp, S5 );
942 	  expectD[double_counter] = atoi (numTemp);
943 	  double_counter++;
944 	}
945 
946       if (cscompare(S3,"integer"))
947 	{
948 	  strcpy( fileTypeI[int_counter], S1 );
949 	  strcpy( fieldNameI[int_counter], S2 );
950 	  strcpy( dataTypeI[int_counter], S3 );
951 	  strcpy( headerTypeI[int_counter], S4 );
952 	  strcpy( numTemp, S5 );
953 	  expectI[int_counter] = atoi (numTemp);
954 	  int_counter++;
955 	}
956     }
957 
958   MPI_Barrier(MPI_COMM_WORLD);
959   if(myrank==0){
960     printf("Starting to read some blocks (doubles) in the restart.dat.##.## files\n");
961   }
962 
963   for ( i = 0; i < N_restart_double; i++ )
964     {
965       for ( j = 0; j < nppp; j++ )
966 	{
967 
968           if(existsubdir) {
969             subdir = (startpart+j-1) / DIR_FANOUT;
970 	    sprintf(gfname,"./%d-procs_case/%d/restart.%d.%d",N_parts, subdir, N_steps,startpart+j);
971           }
972           else
973 	    sprintf(gfname,"./%d-procs_case/restart.%d.%d",N_parts, N_steps, startpart+j);
974 
975 	  openfile(gfname,"read",&irestart);
976 
977 	  for ( k = 0; k < 10; k++ )
978 	    iarray[k]=0;
979 
980 	  paraD[i][j] = new int[expectD[i]];
981 
982 	  iarray[0]=-1;
983 	  readheader( &irestart,
984 		       fieldNameD[i],
985 		       (void*)iarray,
986 		       &expectD[i],
987 		       "double",
988 		       "binary" );
989 
990 	  for ( k = 0; k < expectD[i]; k++ )
991 	    paraD[i][j][k] = iarray[k];
992 
993 	  if ( iarray[0]==-1 )
994 	      WriteLockD[i]=1;
995 	  if ( WriteLockD[i]==0 )
996 	    {
997 	      if ( cscompare("block",headerTypeD[i]) )
998 		{
999 		  if ( expectD[i]==1)
1000 		    isize = paraD[i][j][0];
1001 		  else
1002 		    isize = paraD[i][j][0] * paraD[i][j][1];
1003 
1004 		  Dfield[i][j] = new double[isize];
1005 		  readdatablock( &irestart,
1006 				  fieldNameD[i],
1007 				  (void*)Dfield[i][j],
1008 				  &isize,
1009 				  "double",
1010 				  "binary" );
1011 
1012 		}
1013 	    }
1014 	  closefile(&irestart, "read");
1015 	}
1016     }
1017 
1018   MPI_Barrier(MPI_COMM_WORLD);
1019   if(myrank==0){
1020     printf("Starting to read some blocks (integers) in the restart.dat.##.## files\n");
1021   }
1022 
1023 
1024   for ( i = 0; i < N_restart_integer; i++ )
1025     {
1026       for ( j = 0; j < nppp; j++ )
1027 	{
1028 
1029           if(existsubdir) {
1030             subdir = (startpart+j-1) / DIR_FANOUT;
1031 	    sprintf(gfname,"./%d-procs_case/%d/restart.%d.%d",N_parts, subdir, N_steps, startpart+j);
1032           }
1033           else
1034 	    sprintf(gfname,"./%d-procs_case/restart.%d.%d",N_parts, N_steps, startpart+j);
1035 
1036 	  openfile(gfname,"read",&irestart);
1037 
1038 	  for ( k = 0; k < 10; k++ )
1039 	    iarray[k]=0;
1040 
1041 	  paraI[i][j] = new int[expectI[i]];
1042 
1043 	  iarray[0]=-1;
1044 	  readheader( &irestart,
1045 		       fieldNameI[i],
1046 		       (void*)iarray,
1047 		       &expectI[i],
1048 		       "integer",
1049 		       "binary" );
1050 
1051 	  for ( k = 0; k < expectI[i]; k++ )
1052 	    paraI[i][j][k] = iarray[k];
1053 
1054 	  if ( iarray[0]==-1 )
1055 	      WriteLockI[i]=1;
1056 	  if ( WriteLockI[i]==0 )
1057 	    {
1058 
1059 	      if ( cscompare("block",headerTypeI[i]) )
1060 		{
1061 		  if ( expectI[i]==1)
1062 		    isize = paraI[i][j][0];
1063 		  else
1064 		    isize = paraI[i][j][0] * paraI[i][j][1];
1065 
1066 		  Ifield[i][j] = new int[isize];
1067 		  readdatablock( &irestart,
1068 				  fieldNameI[i],
1069 				  (void*)Ifield[i][j],
1070 				  &isize,
1071 				  "integer",
1072 				  "binary" );
1073 		}
1074 	    }
1075 	  closefile(&irestart, "read");
1076 	}
1077     }
1078 
1079   fclose(pFile);
1080 
1081   ///////////////////// Writing ///////////////////////////////
1082 
1083   int N_restart = N_restart_double + N_restart_integer;
1084 
1085   bzero((void*)fname,255);
1086 //MR CHANGE
1087 //  sprintf(fname,"./%d-procs_case/restart-dat.%d.%d",N_parts,N_steps,((int)(myrank/(N_procs/N_files))+1));
1088   sprintf(fname,"./%d-procs_case-SyncIO-%d/restart-dat.%d.%d",N_parts,N_files,N_steps,((int)(myrank/(N_procs/N_files))+1));
1089 //MR CHANGE END
1090   initphmpiio(&N_restart, &nppf, &N_files,&writeHandle, "write");
1091 //  initphmpiio(&N_restart, &nppf, &N_files,&writeHandle);
1092   openfile(fname, "write", &writeHandle);
1093 
1094   MPI_Barrier(MPI_COMM_WORLD);
1095   if(myrank==0){
1096     printf("Starting to write some blocks (doubles) in the restart-dat.##.## files\n");
1097   }
1098 
1099   for ( i = 0; i < N_restart_double; i++ )
1100     {
1101       for (  j = 0; j < nppp; j++  )
1102 	{
1103 
1104 	  if (WriteLockD[i]==0)
1105 	    {
1106 	      GPID = startpart + j;
1107 	      bzero((void*)fieldtag,255);
1108 	      sprintf(fieldtag,"%s@%d",fieldNameD[i],GPID);
1109 
1110 	      if ( expectD[i]==1)
1111 		isize = paraD[i][j][0];
1112 	      else
1113 		isize = paraD[i][j][0] * paraD[i][j][1];
1114 
1115 	      for ( k = 0; k < expectD[i]; k++ )
1116 		iarray[k] = paraD[i][j][k];
1117 
1118 	      if ( cscompare("header",headerTypeD[i]) )
1119 		isize = 0;
1120 
1121 	      writeheader( &writeHandle,
1122 			    fieldtag,
1123 			    (void*)iarray,
1124 			    &expectD[i],
1125 			    &isize,
1126 			    "double",
1127 			    "binary");
1128 
1129 	      writedatablock( &writeHandle,
1130 			       fieldtag,
1131 			       (void*)Dfield[i][j],
1132 			       &isize,
1133 			       "double",
1134 			       "binary" );
1135 	      if ( cscompare("block",headerTypeD[i]) )
1136 		delete [] Dfield[i][j];
1137 	    }
1138 	  delete [] paraD[i][j];
1139 	}
1140     }
1141 
1142   MPI_Barrier(MPI_COMM_WORLD);
1143   if(myrank==0){
1144     printf("Starting to write some blocks (integers) in the restart.dat.##.## files\n");
1145   }
1146 
1147   for ( i = 0; i < N_restart_integer; i++ )
1148     {
1149       for (  j = 0; j < nppp; j++  )
1150 	{
1151 
1152 	  if (WriteLockI[i]==0)
1153 	    {
1154 	      GPID = startpart + j;
1155 	      bzero((void*)fieldtag,255);
1156 	      sprintf(fieldtag,"%s@%d",fieldNameI[i],GPID);
1157 
1158 	      if ( expectI[i]==1)
1159 		isize = paraI[i][j][0];
1160 	      else
1161 		isize = paraI[i][j][0] * paraI[i][j][1];
1162 
1163 	      for ( k = 0; k < expectI[i]; k++ )
1164 		iarray[k] = paraI[i][j][k];
1165 
1166 	      if ( cscompare("header",headerTypeI[i]) )
1167 		isize = 0;
1168 
1169 	      writeheader( &writeHandle,
1170 			    fieldtag,
1171 			    (void*)iarray,
1172 			    &expectI[i],
1173 			    &isize,
1174 			    "integer",
1175 			    "binary");
1176 
1177 	      writedatablock( &writeHandle,
1178 			       fieldtag,
1179 			       (void*)Ifield[i][j],
1180 			       &isize,
1181 			       "integer",
1182 			       "binary" );
1183 
1184 	      if ( cscompare("block",headerTypeI[i]) )
1185 		delete [] Ifield[i][j];
1186 	    }
1187 	  delete [] paraI[i][j];
1188 	}
1189 
1190     }
1191 
1192   MPI_Barrier(MPI_COMM_WORLD);
1193   if(myrank==0){
1194     printf("Closing restart-dat.##.## files\n");
1195   }
1196 
1197   closefile(&writeHandle, "write");
1198   finalizephmpiio(&writeHandle);
1199 
1200   MPI_Barrier(MPI_COMM_WORLD);
1201   if(myrank==0){
1202     printf("Free memory related to restart-dat.##.## files\n");
1203   }
1204 
1205 
1206   for ( i = 0; i < N_restart_double; i++ )
1207     {
1208       delete [] Dfield[i];
1209       delete [] paraD[i];
1210 
1211       delete [] fieldNameD[i];
1212       delete [] fileTypeD[i];
1213       delete [] dataTypeD[i];
1214       delete [] headerTypeD[i];
1215     }
1216 
1217   for ( i = 0; i < N_restart_integer; i++ )
1218     {
1219       delete [] Ifield[i];
1220       delete [] paraI[i];
1221 
1222       delete [] fieldNameI[i];
1223       delete [] fileTypeI[i];
1224       delete [] dataTypeI[i];
1225       delete [] headerTypeI[i];
1226     }
1227 
1228   delete [] Dfield;
1229   delete [] Ifield;
1230 
1231   delete [] paraD;
1232   delete [] paraI;
1233 
1234   delete [] expectD;
1235   delete [] expectI;
1236 
1237   delete [] fieldNameD;
1238   delete [] fileTypeD;
1239   delete [] dataTypeD;
1240   delete [] headerTypeD;
1241 
1242   delete [] fieldNameI;
1243   delete [] fileTypeI;
1244   delete [] dataTypeI;
1245   delete [] headerTypeI;
1246 
1247 
1248   if (myrank==0)
1249     {
1250       printf("\nFinished transfer, please check data using:\n");
1251       printf(" grep -a ': <' filename \n\n");
1252       printf("Note that the size of the fields is computed based on previous geombc files\n");
1253       printf("Check the routine 'computenitems' if you have any reason to think it has changes for the fields you are interested in\n\n");
1254     }
1255 
1256   MPI_Finalize();
1257 
1258 }
1259 
1260 
1261