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