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