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