1 /* This file provides interface functions for 'partial ' random 2 access into the PHASTA input files 3 4 Anil Karanam March 2001 */ 5 6 #include <stdlib.h> 7 #include <stdio.h> 8 #include <string.h> 9 #include <ctype.h> 10 #include <stdlib.h> 11 #include <time.h> 12 #include <math.h> 13 #include <assert.h> 14 #include "mpi.h" 15 #include "phastaIO.h" 16 #include "rdtsc.h" 17 #include <FCMangle.h> 18 #include "new_interface.h" 19 #include "phIO.h" 20 #include "syncio.h" 21 #include "posixio.h" 22 #include "streamio.h" 23 #include "common_c.h" 24 #include "tmrc.h" 25 #include "phString.h" 26 27 #ifdef intel 28 #include <winsock2.h> 29 #else 30 #include <unistd.h> 31 #include <strings.h> 32 #endif 33 34 void igetMinMaxAvg(int *ivalue, double *stats, int *statRanks) { 35 double *value = (double*)malloc(sizeof(double)); 36 *value = 1.0*(*ivalue); 37 rgetMinMaxAvg(value,stats,statRanks); 38 free(value); 39 } 40 41 void rgetMinMaxAvg(double *value, double *stats, int *statRanks) { 42 int isThisRank; 43 double sqValue = 0., sqValueAvg = 0.; 44 45 MPI_Allreduce(value,&stats[0],1,MPI_DOUBLE,MPI_MIN,MPI_COMM_WORLD); 46 isThisRank=workfc.numpe+1; 47 if(*value==stats[0]) 48 isThisRank=workfc.myrank; 49 MPI_Allreduce(&isThisRank,&statRanks[0],1,MPI_INT,MPI_MIN,MPI_COMM_WORLD); 50 51 MPI_Allreduce(value,&stats[1],1,MPI_DOUBLE,MPI_MAX,MPI_COMM_WORLD); 52 isThisRank=workfc.numpe+1; 53 if(*value==stats[1]) 54 isThisRank=workfc.myrank; 55 MPI_Allreduce(&isThisRank,&statRanks[1],1,MPI_INT,MPI_MIN,MPI_COMM_WORLD); 56 57 MPI_Allreduce(value,&stats[2],1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); 58 stats[2] /= workfc.numpe; 59 60 sqValue = (*value)*(*value); 61 MPI_Allreduce(&sqValue,&sqValueAvg,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD); 62 sqValueAvg /= workfc.numpe; 63 64 stats[3] = sqrt(sqValueAvg-stats[2]*stats[2]); 65 } 66 67 void print_mesh_stats(void) { 68 int statRanks[2]; 69 double iStats[4]; 70 71 igetMinMaxAvg(&conpar.nshg,iStats,statRanks); 72 if(workfc.myrank==workfc.master) 73 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]); 74 igetMinMaxAvg(&conpar.numel,iStats,statRanks); 75 if(workfc.myrank==workfc.master) 76 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]); 77 igetMinMaxAvg(&conpar.numelb,iStats,statRanks); 78 if(workfc.myrank==workfc.master) 79 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]); 80 igetMinMaxAvg(&conpar.nnz_tot,iStats,statRanks); 81 if(workfc.myrank==workfc.master) { 82 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]); 83 printf("\n"); 84 } 85 } 86 87 void print_mpi_stats(void) { 88 int statRanks[2]; 89 double iStats[4], rStats[4]; 90 91 /* NS equations*/ 92 igetMinMaxAvg(&mpistats.iISend,iStats,statRanks); 93 if(workfc.myrank==workfc.master) 94 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]); 95 igetMinMaxAvg(&mpistats.iIRecv,iStats,statRanks); 96 if(workfc.myrank==workfc.master) 97 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]); 98 igetMinMaxAvg(&mpistats.iWaitAll,iStats,statRanks); 99 if(workfc.myrank==workfc.master) 100 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]); 101 igetMinMaxAvg(&mpistats.iAllR,iStats,statRanks); 102 if(workfc.myrank==workfc.master) 103 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]); 104 105 rgetMinMaxAvg(&mpistats.rISend,rStats,statRanks); 106 if(workfc.myrank==workfc.master) 107 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]); 108 rgetMinMaxAvg(&mpistats.rIRecv,rStats,statRanks); 109 if(workfc.myrank==workfc.master) 110 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]); 111 rgetMinMaxAvg(&mpistats.rWaitAll,rStats,statRanks); 112 if(workfc.myrank==workfc.master) 113 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]); 114 rgetMinMaxAvg(&mpistats.rCommu,rStats,statRanks); 115 if(workfc.myrank==workfc.master) 116 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]); 117 rgetMinMaxAvg(&mpistats.rAllR,rStats,statRanks); 118 if(workfc.myrank==workfc.master) { 119 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]); 120 printf("\n"); 121 } 122 /* Scalars*/ 123 igetMinMaxAvg(&mpistats.iISendScal,iStats,statRanks); 124 if(workfc.myrank==workfc.master) 125 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]); 126 igetMinMaxAvg(&mpistats.iIRecvScal,iStats,statRanks); 127 if(workfc.myrank==workfc.master) 128 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]); 129 igetMinMaxAvg(&mpistats.iWaitAllScal,iStats,statRanks); 130 if(workfc.myrank==workfc.master) 131 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]); 132 igetMinMaxAvg(&mpistats.iAllRScal,iStats,statRanks); 133 if(workfc.myrank==workfc.master) 134 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]); 135 136 rgetMinMaxAvg(&mpistats.rISendScal,rStats,statRanks); 137 if(workfc.myrank==workfc.master) 138 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]); 139 rgetMinMaxAvg(&mpistats.rIRecvScal,rStats,statRanks); 140 if(workfc.myrank==workfc.master) 141 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]); 142 rgetMinMaxAvg(&mpistats.rWaitAllScal,rStats,statRanks); 143 if(workfc.myrank==workfc.master) 144 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]); 145 rgetMinMaxAvg(&mpistats.rCommuScal,rStats,statRanks); 146 if(workfc.myrank==workfc.master) 147 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]); 148 rgetMinMaxAvg(&mpistats.rAllRScal,rStats,statRanks); 149 if(workfc.myrank==workfc.master) 150 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]); 151 152 153 } 154 155 void print_system_stats(double *tcorecp, double *tcorecpscal) { 156 int statRanks[2]; 157 double rStats[4]; 158 double syst_assembly, syst_solve; 159 160 /* NS equations */ 161 syst_assembly = tcorecp[0]; 162 syst_solve = tcorecp[1]; 163 164 rgetMinMaxAvg(&syst_assembly,rStats,statRanks); 165 if(workfc.myrank==workfc.master) 166 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]); 167 168 rgetMinMaxAvg(&syst_solve,rStats,statRanks); 169 if(workfc.myrank==workfc.master) 170 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]); 171 172 /* Scalars */ 173 syst_assembly = tcorecpscal[0]; 174 syst_solve = tcorecpscal[1]; 175 176 rgetMinMaxAvg(&syst_assembly,rStats,statRanks); 177 if(workfc.myrank==workfc.master) 178 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]); 179 180 rgetMinMaxAvg(&syst_solve,rStats,statRanks); 181 if(workfc.myrank==workfc.master) { 182 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]); 183 printf("\n"); 184 } 185 } 186 187 188 189 void countfieldstowriterestart() 190 { 191 int nfields = 3; /*magic number, solution, time derivatives*/ 192 193 if(outpar.ivort == 1){ 194 nfields++; /*vorticity*/ 195 } 196 197 if(abs(turbvar.itwmod) != 1 && outpar.iowflux == 1) { 198 nfields++; /*instantaneous wss in bflux.f */ 199 } 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 outpar.nsynciofieldswriterestart = nfields; 231 232 if(workfc.myrank == 0) { 233 printf("Number of fields to write in restart files: %d\n", nfields); 234 } 235 } 236 237 238 void 239 Write_Restart( int* pid, 240 int* stepno, 241 int* nshg, 242 int* numVars, 243 double* array1, 244 double* array2 ) { 245 246 const char* magic_name = "byteorder magic number"; 247 int magic_number = 362436; 248 int isize, nitems; 249 int iarray[10]; 250 int nfiles; 251 int nfields; 252 int numparts; 253 int nprocs; 254 int ione = 1; 255 double iotime = 0; 256 char filename[255]; 257 258 /* First, count the number of fields to write and store the result in*/ 259 countfieldstowriterestart(); 260 261 /* Retrieve and compute the parameters required for SyncIO*/ 262 nfiles = outpar.nsynciofiles; 263 nfields = outpar.nsynciofieldswriterestart; 264 numparts = workfc.numpe; 265 nprocs = workfc.numpe; 266 assert(numparts/nprocs == 1);/* Number of parts per proc ...*/ 267 268 bzero((void*)filename,255); 269 270 iotime = TMRC(); 271 if(outpar.output_mode == -1 ) 272 streamio_setup_write(&f_descriptor, streamio_get_r()); 273 else if(outpar.output_mode == 0 ) 274 posixio_setup(&f_descriptor, 'w'); 275 else if(outpar.output_mode > 0 ) 276 syncio_setup_write(nfiles, nfields, numparts/nfiles, &f_descriptor); 277 else 278 exit(EXIT_FAILURE); 279 phio_constructName(f_descriptor,"restart",filename); 280 phstr_appendInt(filename, *stepno); 281 phstr_appendStr(filename, "."); 282 phio_openfile(filename, f_descriptor); 283 284 field_flag=0; 285 286 /* write the magic number*/ 287 phio_writeheader(f_descriptor, magic_name, (void*)&magic_number, &ione, 288 &ione, "integer", phasta_iotype); 289 phio_writedatablock(f_descriptor, magic_name, (void*)&magic_number, 290 &ione, "integer", phasta_iotype ); 291 field_flag++; 292 293 /* Write solution field ...*/ 294 isize = (*nshg)*(*numVars); 295 nitems = 3; 296 iarray[ 0 ] = (*nshg); 297 iarray[ 1 ] = (*numVars); 298 iarray[ 2 ] = (*stepno); 299 300 phio_writeheader(f_descriptor, "solution", (void*)iarray, &nitems, 301 &isize, "double", phasta_iotype); 302 nitems = (*nshg)*(*numVars); 303 phio_writedatablock(f_descriptor, "solution", (void*)(array1), 304 &isize, "double", phasta_iotype ); 305 field_flag++; 306 307 /* Write solution field ...*/ 308 isize = (*nshg)*(*numVars); 309 nitems = 3; 310 iarray[ 0 ] = (*nshg); 311 iarray[ 1 ] = (*numVars); 312 iarray[ 2 ] = (*stepno); 313 phio_writeheader(f_descriptor, "time derivative of solution", 314 (void*)iarray, &nitems, &isize, "double", phasta_iotype); 315 nitems = (*nshg)*(*numVars); 316 phio_writedatablock(f_descriptor, "time derivative of solution", 317 (void*)(array2), &isize, "double", phasta_iotype ); 318 field_flag++; 319 320 if (field_flag==nfields){ 321 phio_closefile(f_descriptor); 322 if (*pid==0) { 323 printf("\n"); 324 } 325 } 326 iotime = TMRC() - iotime; 327 if (workfc.master == workfc.myrank) 328 printf("time to write restart (seconds) %f\n",iotime); 329 } 330 331 void 332 Write_Error( int* pid, 333 int* stepno, 334 int* nshg, 335 int* numVars, 336 double* array1 ) { 337 int isize, nitems; 338 int iarray[10]; 339 int nfields; 340 int numparts; 341 int nprocs; 342 343 (void)*pid; /*silence compiler warning*/ 344 345 nfields = outpar.nsynciofieldswriterestart; 346 numparts = workfc.numpe; 347 nprocs = workfc.numpe; 348 349 assert(numparts/nprocs == 1);/* Number of parts per proc ...*/ 350 351 field_flag++; 352 353 if(*pid==0) { 354 printf("\n"); 355 printf("The %d/%d th field to be written is 'errors'\n",field_flag,nfields); 356 } 357 358 isize = (*nshg)*(*numVars); 359 nitems = 3; 360 iarray[ 0 ] = (*nshg); 361 iarray[ 1 ] = (*numVars); 362 iarray[ 2 ] = (*stepno); 363 364 phio_writeheader(f_descriptor, "errors", (void*)iarray, &nitems, 365 &isize, "double", phasta_iotype); 366 367 phio_writedatablock(f_descriptor, "errors", (void*)array1, &isize, 368 "double", phasta_iotype ); 369 370 if (field_flag==nfields){ 371 phio_closefile(f_descriptor); 372 if (*pid==0) { 373 printf("Last field %d 'errors' finished! \n",nfields); 374 printf("\n"); 375 } 376 } 377 } 378 379 void 380 Write_Field( int *pid, 381 char* filemode, 382 char* fieldtag, 383 int* tagsize, 384 void* array, 385 char* arraytype, 386 int* nshg, 387 int* numvars, 388 int* stepno) { 389 char *fieldlabel = NULL; 390 391 int isize, nitems; 392 int iarray[10]; 393 char datatype[10]; 394 395 int nfields; 396 int numparts; 397 int nprocs; 398 399 (void)*filemode; /*silence compiler warning*/ 400 401 if(!strncmp(arraytype,"i",1)) 402 strcpy(datatype,"int"); 403 else /* default is double*/ 404 strcpy(datatype,"double"); 405 406 nfields = outpar.nsynciofieldswriterestart; 407 numparts = workfc.numpe; 408 nprocs = workfc.numpe; 409 410 assert(numparts/nprocs == 1);/* Number of parts per proc ...*/ 411 412 fieldlabel = (char *)malloc((*tagsize+1)*sizeof(char)); 413 strncpy(fieldlabel, fieldtag, *tagsize); 414 fieldlabel[*tagsize] = '\0'; 415 416 field_flag++; 417 if(*pid==0) { 418 printf("\n"); 419 printf("The %d/%d th field to be written is '%s'\n",field_flag,nfields,fieldlabel); 420 } 421 422 /* Write solution field ...*/ 423 isize = (*nshg)*(*numvars); 424 nitems = 3; 425 iarray[ 0 ] = (*nshg); 426 iarray[ 1 ] = (*numvars); 427 iarray[ 2 ] = (*stepno); 428 429 phio_writeheader(f_descriptor, fieldlabel, (void*)iarray, &nitems, 430 &isize, datatype, phasta_iotype); 431 nitems = (*nshg)*(*numvars); 432 phio_writedatablock(f_descriptor, fieldlabel, array, &isize, 433 datatype, phasta_iotype ); 434 435 if (field_flag==nfields){ 436 phio_closefile(f_descriptor); 437 if (*pid==0) { 438 printf("Last field %d '%s' finished! \n",nfields, fieldtag); 439 printf("\n"); 440 } 441 } 442 free(fieldlabel); 443 } 444 445 void 446 Write_PhAvg2( int* pid, 447 char* filemode, 448 char* fieldtag, 449 int* tagsize, 450 int* iphase, 451 int* nphasesincycle, 452 void* array, 453 char* arraytype, 454 int* nshg, 455 int* numvars, 456 int* stepno) { 457 int addtagsize=0; /* phase number is added to the name of the field*/ 458 int tagsize2 = 0; 459 char *fieldlabel = NULL; 460 char straddtagsize[10] = "\0"; 461 int isize, nitems; 462 int iarray[10]; 463 char datatype[10]; 464 int nfields; 465 int numparts; 466 int nprocs; 467 468 (void)*filemode; /*silence compiler warning*/ 469 (void)*nphasesincycle; /*silence compiler warning*/ 470 471 if(*iphase<10) 472 addtagsize=1; 473 else if(*iphase<100) 474 addtagsize=2; 475 else if(*iphase<1000) 476 addtagsize=3; 477 478 tagsize2=*tagsize+addtagsize; 479 480 fieldlabel = (char *)malloc((tagsize2+1)*sizeof(char)); 481 strncpy(fieldlabel, fieldtag, *tagsize); 482 fieldlabel[tagsize2] = '\0'; 483 484 sprintf(straddtagsize,"%d",*iphase); 485 486 if(*iphase<10) { 487 fieldlabel[tagsize2-1]=straddtagsize[0]; 488 } 489 else if(*iphase<100) { 490 fieldlabel[tagsize2-2]=straddtagsize[0]; 491 fieldlabel[tagsize2-1]=straddtagsize[1]; 492 } 493 else if(*iphase<1000) { 494 fieldlabel[tagsize2-3]=straddtagsize[0]; 495 fieldlabel[tagsize2-2]=straddtagsize[1]; 496 fieldlabel[tagsize2-1]=straddtagsize[2]; 497 } 498 499 if(!strncmp(arraytype,"i",1)) 500 strcpy(datatype,"int"); 501 else /* default is double*/ 502 strcpy(datatype,"double"); 503 504 nfields = outpar.nsynciofieldswriterestart; 505 numparts = workfc.numpe; 506 nprocs = workfc.numpe; 507 508 assert(numparts/nprocs == 1);/* Number of parts per proc ...*/ 509 510 field_flag++; 511 if(*pid==0) { 512 printf("\n"); 513 printf("The %d/%d th field to be written is '%s'\n",field_flag,nfields,fieldlabel); 514 } 515 516 /* Write solution field ...*/ 517 isize = (*nshg)*(*numvars); 518 nitems = 3; 519 iarray[ 0 ] = (*nshg); 520 iarray[ 1 ] = (*numvars); 521 iarray[ 2 ] = (*stepno); 522 phio_writeheader(f_descriptor, fieldlabel, (void*)iarray, &nitems, 523 &isize, "double", phasta_iotype); 524 nitems = (*nshg)*(*numvars); 525 phio_writedatablock(f_descriptor, fieldlabel, array, &isize, 526 "double", phasta_iotype ); 527 if (field_flag==nfields){ 528 phio_closefile(f_descriptor); 529 if (*pid==0) { 530 printf("\n"); 531 } 532 } 533 free(fieldlabel); 534 } 535