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