1 /*=========================================================================== 2 * 3 * "usr.c": user's function 4 * 5 *=========================================================================== 6 */ 7 #include <stdio.h> 8 #include <stdlib.h> 9 #include "les.h" 10 #include "usr.h" 11 //mr change 12 // #include "common_c.h" //not needed here any more because added in new_interface.c 13 //mr change end 14 #include "common_c.h" 15 #include "phastaIO.h" 16 #include "new_interface.h" 17 #include <FCMangle.h> 18 19 extern char phasta_iotype[80]; 20 extern int field_flag; //new_interface.c //TODO: figure out where these 21 extern int f_descriptor; //new_interface.c // should be 22 23 /*=========================================================================== 24 * 25 * "usrNew": Put all the values in usrHd 26 * 27 * From FORTRAN 28 * 29 * integer usr(100) 30 * dimension aperm(numnp,nperm) 31 * ... 32 * call usrnew( usr, aperm, ..., numnp, ...) 33 * 34 * 35 *=========================================================================== 36 */ 37 #include "mpi.h" 38 static int lmNum = 0; 39 static LesHd lesArray[8]; 40 void usrNew( UsrHd usrHd, 41 int* eqnType, 42 double* aperm, 43 double* atemp, 44 double* resf, 45 double* solinc, 46 double* flowDiag, 47 double* sclrDiag, 48 double* lesP, 49 double* lesQ, 50 int* iBC, 51 double* BC, 52 int* iper, 53 int* ilwork, 54 int* numpe, 55 int* nNodes, 56 int* nenl, 57 int* nPermDims, 58 int* nTmpDims, 59 int* rowp, 60 int* colm, 61 double* lhsK, 62 double* lhsP, 63 double* lhsS, 64 int* nnz_tot, 65 double* CGsol 66 ) 67 { 68 char* funcName = "usrNew" ; /* function name */ 69 70 /*--------------------------------------------------------------------------- 71 * Stick the parameters 72 *--------------------------------------------------------------------------- 73 */ 74 usrHd->eqnType = *eqnType ; 75 usrHd->aperm = aperm ; 76 usrHd->atemp = atemp ; 77 usrHd->resf = resf ; 78 usrHd->solinc = solinc ; 79 usrHd->flowDiag = flowDiag ; 80 usrHd->sclrDiag = sclrDiag ; 81 usrHd->lesP = lesP ; 82 usrHd->lesQ = lesQ ; 83 usrHd->iBC = iBC ; 84 usrHd->BC = BC ; 85 usrHd->iper = iper ; 86 usrHd->ilwork = ilwork ; 87 usrHd->numpe = *numpe ; 88 usrHd->nNodes = *nNodes ; 89 usrHd->nenl = *nenl ; 90 usrHd->nPermDims = *nPermDims ; 91 usrHd->nTmpDims = *nTmpDims ; 92 usrHd->rowp = rowp ; 93 usrHd->colm = colm ; 94 usrHd->lhsK = lhsK ; 95 usrHd->lhsP = lhsP ; 96 usrHd->lhsS = lhsS ; 97 usrHd->nnz_tot = nnz_tot ; 98 usrHd->CGsol = CGsol; 99 } /* end of usrNew() */ 100 101 /*=========================================================================== 102 * 103 * "usrPointer": Get the pointer 104 * 105 *=========================================================================== 106 */ 107 Real* 108 usrPointer( UsrHd usrHd, 109 Integer id, 110 Integer offset, 111 Integer nDims ) 112 { 113 char* funcName = "usrPointer";/* function name */ 114 Real* pnt ; /* pointer */ 115 116 /*--------------------------------------------------------------------------- 117 * Get the head of the memory 118 *--------------------------------------------------------------------------- 119 */ 120 if ( id == LES_RES_PNT ) { 121 122 pnt = usrHd->resf ; 123 id = 0 ; 124 125 } else if ( id == LES_SOL_PNT ) { 126 127 pnt = usrHd->solinc ; 128 id = 0 ; 129 130 } else if ( id < 0 ) { 131 132 pnt = usrHd->aperm ; 133 id = id + usrHd->nPermDims ; 134 135 } else { 136 137 pnt = usrHd->atemp ; 138 id = id ; 139 140 } 141 /*--------------------------------------------------------------------------- 142 * Get the offset 143 *--------------------------------------------------------------------------- 144 */ 145 pnt = pnt + (id + offset) * usrHd->nNodes ; 146 147 /*--------------------------------------------------------------------------- 148 * Return the pointer 149 *--------------------------------------------------------------------------- 150 */ 151 return( pnt ) ; 152 153 } /* end of usrPointer() */ 154 155 #define myflesnew FortranCInterface_GLOBAL_(myflesnew,MYFLESNEW) 156 #define myflessolve FortranCInterface_GLOBAL_(myflessolve,MYFLESSOLVE) 157 #define savelesrestart FortranCInterface_GLOBAL_(savelesrestart,SAVELESRESTART) 158 #define readlesrestart FortranCInterface_GLOBAL_(readlesrestart,READLESRESTART) 159 #define solverlicenseserver FortranCInterface_GLOBAL_(solverlicenseserver,SOLVERLICENSESERVER) 160 161 162 163 #ifdef intel 164 lesArray[ *lesId ] = lesNew( fileName, *lmport, &lmNum, *eqnType, 165 *nDofs, *minIters, *maxIters, *nKvecs, 166 *prjFlag, *nPrjs, *presPrjFlag, *nPresPrjs,presPrec, 167 *tol, *presTol, *verbose, stats, nPermDims, 168 nTmpDims ); 169 return ;} 170 /* the following is a fake function that was required when we moved to 171 a C++ main on in the MS Visual Studio environment. It fails to 172 link because it is looking for this function 173 */ 174 void _CrtDbgReport() { 175 return ;} 176 177 double __vcos_(double fg) { fflush(stdout); printf(" vcos got called \n"); fflush(stdout);} 178 double __vlog_(double fg) { fflush(stdout); printf(" vlog got called \n"); fflush(stdout);} 179 180 181 #endif /* we are in unix land... whew. secretly we have equivalenced fileName and */ 182 183 /* #ifdef LINUX*/ 184 /* void flush_(int* junk ){ return; }*/ 185 /* #endif*/ 186 void myflesnew( Integer* lesId, 187 Integer* lmport, 188 Integer* eqnType, 189 Integer* nDofs, 190 Integer* minIters, 191 Integer* maxIters, 192 Integer* nKvecs, 193 Integer* prjFlag, 194 Integer* nPrjs, 195 Integer* presPrjFlag, 196 Integer* nPresPrjs, 197 Real* tol, 198 Real* presTol, 199 Integer* verbose, 200 Real* stats, 201 Integer* nPermDims, 202 Integer* nTmpDims, 203 char* lmhost ) { 204 int procId; 205 #ifdef AMG 206 int presPrec=1; 207 #else 208 int presPrec=0; 209 #endif 210 MPI_Comm_rank( MPI_COMM_WORLD, &procId ) ; 211 if(lmNum==0){ 212 if(procId==0){ 213 lesArray[ *lesId ] = lesNew( lmhost, *lmport, &lmNum, *eqnType, 214 *nDofs, *minIters, *maxIters, *nKvecs, 215 *prjFlag, *nPrjs, *presPrjFlag, *nPresPrjs,presPrec, 216 *tol, *presTol, *verbose, stats, nPermDims, 217 nTmpDims ); 218 MPI_Bcast( &lmNum, 1, MPI_INT, 0, MPI_COMM_WORLD ) ; 219 } else { 220 MPI_Bcast( &lmNum, 1, MPI_INT, 0, MPI_COMM_WORLD ) ; 221 lesArray[ *lesId ] = lesNew( lmhost, *lmport, &lmNum, *eqnType, 222 *nDofs, *minIters, *maxIters, *nKvecs, 223 *prjFlag, *nPrjs, *presPrjFlag, *nPresPrjs,presPrec, 224 *tol, *presTol, *verbose, stats, nPermDims, 225 nTmpDims ); 226 } 227 } else { 228 lesArray[ *lesId ] = lesNew( lmhost, *lmport, &lmNum, *eqnType, 229 *nDofs, *minIters, *maxIters, *nKvecs, 230 *prjFlag, *nPrjs, *presPrjFlag, *nPresPrjs,presPrec, 231 *tol, *presTol, *verbose, stats, nPermDims, 232 nTmpDims ); 233 } 234 return ; 235 } 236 237 238 void 239 savelesrestart( Integer* lesId, 240 Real* aperm, 241 Integer* nshg, 242 Integer* myrank, 243 Integer* lstep, 244 Integer* nPermDims ) { 245 246 int nPrjs, PrjSrcId; 247 int nPresPrjs, PresPrjSrcId; 248 char filename[255]; 249 int fileHandle=0; 250 int iarray[3]; 251 int size, nitems; 252 double* projVec; 253 int i, j, count; 254 255 // sprintf( filename,"restart.%d.%d", *lstep, *myrank+1 ); 256 // openfile_( filename, "append", &fileHandle ); 257 258 nPrjs = (Integer) lesGetPar( lesArray[ *lesId ], LES_ACT_PRJS ); 259 PrjSrcId = (Integer) lesGetPar( lesArray[ *lesId ], LES_PRJ_VEC_ID ); 260 261 if ( PrjSrcId < 0 ) PrjSrcId += *nPermDims; 262 263 projVec = (double*)malloc( nPrjs * ( *nshg ) * sizeof( double ) ); 264 265 count = 0; 266 for( i = PrjSrcId; i < PrjSrcId+nPrjs; i ++ ) { 267 for( j = 0 ; j < *nshg; j++ ) { 268 projVec[ count++ ] = aperm[ (*nshg) * i + j ]; 269 } 270 } 271 272 //printf("Savelessrestart, nPrjs is %d\n",nPrjs); 273 274 iarray[ 0 ] = *nshg; 275 iarray[ 1 ] = nPrjs; 276 nitems = 2; 277 size = (*nshg)*nPrjs; 278 279 int name_length; 280 name_length = 18; 281 Write_Field(myrank,"a","projection vectors",&name_length, (void *)projVec,"d", nshg, &nPrjs, lstep); 282 283 //writeheader_( &fileHandle, "projection vectors ", (void*)iarray, 284 // &nitems, &size, "double", phasta_iotype ); 285 //nitems = size; 286 //writedatablock_( &fileHandle, "projection vectors ", (void*)projVec, 287 // &nitems, "double", phasta_iotype ); 288 free(projVec); 289 290 /*************************************************************************/ 291 292 nPresPrjs = (Integer) lesGetPar( lesArray[ *lesId ], LES_ACT_PRES_PRJS ); 293 PresPrjSrcId =(Integer)lesGetPar( lesArray[ *lesId ], LES_PRES_PRJ_VEC_ID ); 294 if ( PresPrjSrcId < 0 ) PresPrjSrcId += *nPermDims; 295 296 projVec = (double*)malloc( nPresPrjs * ( *nshg ) * sizeof( double ) ); 297 298 count = 0; 299 for( i = PresPrjSrcId; i < (PresPrjSrcId + nPresPrjs) ; i ++ ) { 300 for( j = 0 ; j < *nshg; j++ ) { 301 projVec[ count++ ] = aperm[ (*nshg) * i + j ]; 302 } 303 } 304 305 iarray[ 0 ] = *nshg; 306 iarray[ 1 ] = nPresPrjs; 307 nitems = 2; 308 size = (*nshg)*nPresPrjs; 309 310 name_length = 27; 311 Write_Field(myrank,"a","pressure projection vectors",&name_length, projVec,"d", nshg, &nPresPrjs, lstep); 312 313 //writeheader_( &fileHandle, "pressure projection vectors ", (void*)iarray, 314 // &nitems, &size, "double", phasta_iotype ); 315 // nitems = size; 316 317 //writedatablock_( &fileHandle, "pressure projection vectors ", 318 // (void*)projVec, &nitems, "double", phasta_iotype ); 319 free( projVec); 320 321 //closefile_( &fileHandle, "append" ); 322 } 323 324 void 325 readlesrestart( Integer* lesId, 326 Real* aperm, 327 Integer* nshg, 328 Integer* myrank, 329 Integer* lstep , 330 Integer* nPermDims ) { 331 332 int nPrjs, PrjSrcId; 333 int nPresPrjs, PresPrjSrcId; 334 char filename[255]; 335 int fileHandle=0; 336 int iarray[3]={-1,-1,-1}; 337 int size, nitems; 338 int itwo=2; 339 int lnshg; 340 double* projVec; 341 int i,j,count; 342 343 //MR CHANGE 344 // sprintf( filename,"restart.%d.%d", *lstep, *myrank+1 ); 345 346 // int nfiles=2; 347 // int numParts=8; 348 // int nfields=6; 349 int nfiles; 350 int nfields; 351 int numParts; 352 int nprocs; 353 int nppf; 354 355 nfiles = outpar.nsynciofiles; 356 // nfields = outpar.nsynciofieldsreadrestart; 357 numParts = workfc.numpe; 358 nprocs = workfc.numpe; 359 //MR CHANGE END 360 361 // int nppf = numParts/nfiles; 362 char fieldname[255]; 363 364 // MPI_Comm_size(MPI_COMM_WORLD, &numprocs); 365 // Calculate number of parts each proc deal with and where it start and end ... 366 int nppp = numParts/nprocs; // nppp : Number of parts per proc ... 367 int startpart = *myrank * nppp +1; // Part id from which I (myrank) start ... 368 int endpart = startpart + nppp - 1; // Part id to which I (myrank) end ... 369 370 sprintf( filename,"restart-dat.%d.%d",*lstep,((int)(*myrank/(numParts/nfiles))+1)); 371 queryphmpiio(filename, &nfields, &nppf); 372 initphmpiio(&nfields, &nppf, &nfiles,&fileHandle,"read"); 373 //MR CHANGE END 374 375 openfile( filename, "read", &fileHandle ); 376 377 378 // if ( fileHandle == 0 ) return; 379 if ( fileHandle < 0 ) return; // See phastaIO.cc for error fileHandle 380 // for ( i = 0; i < nppp; i++ ) // If one day several parts per file 381 // { 382 bzero((void*)fieldname,255); 383 sprintf(fieldname,"projection vectors@%d",startpart); 384 // readheader_( &fileHandle, "projection vectors", (void*)iarray, 385 // &itwo, "integer", phasta_iotype ); 386 readheader( &fileHandle, fieldname, (void*)iarray, 387 &itwo, "integer", phasta_iotype ); 388 // } 389 //MR CHANGE END 390 391 if ( iarray[0] != *nshg ) { 392 closefile( &fileHandle, "read" ); 393 if(workfc.myrank==workfc.master) 394 printf("projection vectors are being initialized to zero (SAFE)\n"); 395 return; 396 } 397 398 lnshg = iarray[ 0 ] ; 399 nPrjs = iarray[ 1 ] ; 400 401 size = (*nshg)*nPrjs; 402 projVec = (double*)malloc( size * sizeof( double )); 403 404 //MR CHANGE 405 // readdatablock_( &fileHandle, "projection vectors", (void*)projVec, 406 // &size, "double", phasta_iotype ); 407 readdatablock( &fileHandle, fieldname, (void*)projVec, 408 &size, "double", phasta_iotype ); 409 //MR CHANGE END 410 411 lesSetPar( lesArray[ *lesId ], LES_ACT_PRJS, (Real) nPrjs ); 412 PrjSrcId = (Integer) lesGetPar( lesArray[ *lesId ], LES_PRJ_VEC_ID ); 413 if ( PrjSrcId < 0 ) PrjSrcId += *nPermDims; 414 415 count = 0; 416 for( i = PrjSrcId; i < PrjSrcId+nPrjs; i ++ ) { 417 for( j = 0 ; j < *nshg; j++ ) { 418 aperm[ (*nshg) * i + j ] = projVec[ count++ ] ; 419 } 420 } 421 422 free( projVec ); 423 424 /************************************************************************/ 425 426 iarray[0] = -1; iarray[1] = -1; iarray[2] = -1; 427 428 //MR CHANGE 429 430 bzero((void*)fieldname,255); 431 sprintf(fieldname,"pressure projection vectors@%d",startpart); 432 433 //MR CHANGE END 434 435 436 readheader( &fileHandle, fieldname, (void*)iarray, 437 &itwo, "integer", phasta_iotype ); 438 439 lnshg = iarray[ 0 ] ; 440 nPresPrjs = iarray[ 1 ] ; 441 442 if ( lnshg != *nshg ) { 443 closefile( &fileHandle, "read" ); 444 if(workfc.myrank==workfc.master) 445 printf("pressure projection vectors are being initialized to zero (SAFE)\n"); 446 return; 447 } 448 449 size = (*nshg)*nPresPrjs; 450 projVec = (double*)malloc( size * sizeof( double )); 451 452 readdatablock( &fileHandle, fieldname, (void*)projVec, 453 &size, "double", phasta_iotype ); 454 455 lesSetPar( lesArray[ *lesId ], LES_ACT_PRES_PRJS, (Real) nPresPrjs ); 456 PresPrjSrcId=(Integer)lesGetPar( lesArray[ *lesId ], LES_PRES_PRJ_VEC_ID ); 457 if ( PresPrjSrcId < 0 ) PresPrjSrcId += *nPermDims; 458 459 count = 0; 460 for( i = PresPrjSrcId; i < PresPrjSrcId+nPresPrjs; i ++ ) { 461 for( j = 0 ; j < *nshg; j++ ) { 462 aperm[ (*nshg) * i + j ] = projVec[ count++ ] ; 463 } 464 } 465 466 free( projVec ); 467 468 closefile( &fileHandle, "read" ); 469 470 //MR CHANGE 471 finalizephmpiio(&fileHandle); 472 //MR CHANGE END 473 474 } 475 476 void myflessolve( Integer* lesId, 477 UsrHd usrHd){ 478 lesSolve( lesArray[ *lesId ], usrHd ); 479 } 480 481 482 int solverlicenseserver(char key[]){ 483 #ifdef intel 484 strcpy(key,"C:\\cygwin\\license.dat"); 485 #else 486 char* env_server_name; 487 env_server_name = getenv("LES_LICENSE_SERVER"); 488 if(env_server_name) strcpy(key, env_server_name); 489 else { 490 if(workfc.myrank==workfc.master) { 491 fprintf(stderr,"environment variable LES_LICENSE_SERVER not defined \n"); 492 fprintf(stderr,"using wesley as default \n"); 493 } 494 //MR CHANGE 495 // strcpy(key, "wesley.scorec.rpi.edu"); 496 strcpy(key, "acusim.license.scorec.rpi.edu"); 497 //MR CHANGE END 498 } 499 #endif 500 return 1; 501 } 502