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 "phString.h"
15 #include "syncio.h"
16 #include "posixio.h"
17 #include "streamio.h"
18 #include "new_interface.h"
19 #include <FCMangle.h>
20
21 extern char phasta_iotype[80];
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];
usrNew(UsrHd usrHd,int * eqnType,double * aperm,double * atemp,double * resf,double * solinc,double * flowDiag,double * sclrDiag,double * lesP,double * lesQ,int * iBC,double * BC,int * iper,int * ilwork,int * numpe,int * nNodes,int * nenl,int * nPermDims,int * nTmpDims,int * rowp,int * colm,double * lhsK,double * lhsP,double * lhsS,int * nnz_tot,double * CGsol)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*
usrPointer(UsrHd usrHd,Integer id,Integer offset,Integer nDims)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 iarray[3];
250 int size, nitems;
251 double* projVec;
252 int i, j, count;
253
254 nPrjs = (Integer) lesGetPar( lesArray[ *lesId ], LES_ACT_PRJS );
255 PrjSrcId = (Integer) lesGetPar( lesArray[ *lesId ], LES_PRJ_VEC_ID );
256
257 if ( PrjSrcId < 0 ) PrjSrcId += *nPermDims;
258
259 projVec = (double*)malloc( nPrjs * ( *nshg ) * sizeof( double ) );
260
261 count = 0;
262 for( i = PrjSrcId; i < PrjSrcId+nPrjs; i ++ ) {
263 for( j = 0 ; j < *nshg; j++ ) {
264 projVec[ count++ ] = aperm[ (*nshg) * i + j ];
265 }
266 }
267
268 iarray[ 0 ] = *nshg;
269 iarray[ 1 ] = nPrjs;
270 nitems = 2;
271 size = (*nshg)*nPrjs;
272
273 int name_length;
274 name_length = 18;
275 Write_Field(myrank,"a","projection vectors",&name_length, (void *)projVec,"d", nshg, &nPrjs, lstep);
276
277 free(projVec);
278
279 nPresPrjs = (Integer) lesGetPar( lesArray[ *lesId ], LES_ACT_PRES_PRJS );
280 PresPrjSrcId =(Integer)lesGetPar( lesArray[ *lesId ], LES_PRES_PRJ_VEC_ID );
281 if ( PresPrjSrcId < 0 ) PresPrjSrcId += *nPermDims;
282
283 projVec = (double*)malloc( nPresPrjs * ( *nshg ) * sizeof( double ) );
284
285 count = 0;
286 for( i = PresPrjSrcId; i < (PresPrjSrcId + nPresPrjs) ; i ++ ) {
287 for( j = 0 ; j < *nshg; j++ ) {
288 projVec[ count++ ] = aperm[ (*nshg) * i + j ];
289 }
290 }
291
292 iarray[ 0 ] = *nshg;
293 iarray[ 1 ] = nPresPrjs;
294 nitems = 2;
295 size = (*nshg)*nPresPrjs;
296
297 name_length = 27;
298 Write_Field(myrank,"a","pressure projection vectors",&name_length, projVec,"d", nshg, &nPresPrjs, lstep);
299
300 free( projVec);
301 }
302
303 void
304 readlesrestart( Integer* lesId,
305 Real* aperm,
306 Integer* nshg,
307 Integer* myrank,
308 Integer* lstep ,
309 Integer* nPermDims ) {
310
311 int nPrjs, PrjSrcId;
312 int nPresPrjs, PresPrjSrcId;
313 char filename[255];
314 phio_fp fileHandle = NULL;
315 int iarray[3]={-1,-1,-1};
316 int size, nitems;
317 int itwo=2;
318 int lnshg;
319 double* projVec;
320 int i,j,count;
321 int nfields;
322 int numParts;
323 int nprocs;
324 int nppf;
325
326 numParts = workfc.numpe;
327 nprocs = workfc.numpe;
328 // Calculate number of parts each proc deal with and where it start and end ...
329 int nppp = numParts/nprocs; // nppp : Number of parts per proc ...
330 int startpart = *myrank * nppp +1; // Part id from which I (myrank) start ...
331 int endpart = startpart + nppp - 1; // Part id to which I (myrank) end ...
332
333 if( outpar.input_mode == -1 )
334 streamio_setup_read(&fileHandle, streamio_get_gr());
335 else if( outpar.input_mode == 0 )
336 posixio_setup(&fileHandle, 'r');
337 else if( outpar.input_mode > 0 )
338 syncio_setup_read(outpar.nsynciofiles, &fileHandle);
339 phio_constructName(fileHandle,"restart",filename);
340 phstr_appendInt(filename, *lstep);
341 phstr_appendStr(filename, ".");
342 phio_openfile(filename, fileHandle);
343
344 if ( !fileHandle ) return; // See phastaIO.cc for error fileHandle
345 phio_readheader(fileHandle, "projection vectors", (void*)iarray,
346 &itwo, "integer", phasta_iotype);
347
348 if ( iarray[0] != *nshg ) {
349 phio_closefile(fileHandle);
350 if(workfc.myrank==workfc.master)
351 printf("projection vectors are being initialized to zero (SAFE)\n");
352 return;
353 }
354
355 lnshg = iarray[ 0 ] ;
356 nPrjs = iarray[ 1 ] ;
357
358 size = (*nshg)*nPrjs;
359 projVec = (double*)malloc( size * sizeof( double ));
360
361 phio_readdatablock(fileHandle, "projection vectors", (void*)projVec,
362 &size, "double", phasta_iotype );
363
364 lesSetPar( lesArray[ *lesId ], LES_ACT_PRJS, (Real) nPrjs );
365 PrjSrcId = (Integer) lesGetPar( lesArray[ *lesId ], LES_PRJ_VEC_ID );
366 if ( PrjSrcId < 0 ) PrjSrcId += *nPermDims;
367
368 count = 0;
369 for( i = PrjSrcId; i < PrjSrcId+nPrjs; i ++ ) {
370 for( j = 0 ; j < *nshg; j++ ) {
371 aperm[ (*nshg) * i + j ] = projVec[ count++ ] ;
372 }
373 }
374
375 free( projVec );
376
377 iarray[0] = -1; iarray[1] = -1; iarray[2] = -1;
378
379 phio_readheader(fileHandle, "pressure projection vectors", (void*)iarray,
380 &itwo, "integer", phasta_iotype );
381
382 lnshg = iarray[ 0 ] ;
383 nPresPrjs = iarray[ 1 ] ;
384
385 if ( lnshg != *nshg ) {
386 phio_closefile(fileHandle);
387 if(workfc.myrank==workfc.master)
388 printf("pressure projection vectors are being initialized to zero (SAFE)\n");
389 return;
390 }
391
392 size = (*nshg)*nPresPrjs;
393 projVec = (double*)malloc( size * sizeof( double ));
394
395 phio_readdatablock(fileHandle, "pressure projection vectors", (void*)projVec,
396 &size, "double", phasta_iotype );
397
398 lesSetPar( lesArray[ *lesId ], LES_ACT_PRES_PRJS, (Real) nPresPrjs );
399 PresPrjSrcId=(Integer)lesGetPar( lesArray[ *lesId ], LES_PRES_PRJ_VEC_ID );
400 if ( PresPrjSrcId < 0 ) PresPrjSrcId += *nPermDims;
401
402 count = 0;
403 for( i = PresPrjSrcId; i < PresPrjSrcId+nPresPrjs; i ++ ) {
404 for( j = 0 ; j < *nshg; j++ ) {
405 aperm[ (*nshg) * i + j ] = projVec[ count++ ] ;
406 }
407 }
408
409 free( projVec );
410
411 phio_closefile(fileHandle);
412 }
413
414 void myflessolve( Integer* lesId,
415 UsrHd usrHd){
416 lesSolve( lesArray[ *lesId ], usrHd );
417 }
418
419
420 int solverlicenseserver(char key[]){
421 #ifdef intel
422 strcpy(key,"C:\\cygwin\\license.dat");
423 #else
424 char* env_server_name;
425 env_server_name = getenv("LES_LICENSE_SERVER");
426 if(env_server_name) strcpy(key, env_server_name);
427 else {
428 if(workfc.myrank==workfc.master) {
429 fprintf(stderr,"environment variable LES_LICENSE_SERVER not defined \n");
430 fprintf(stderr,"using wesley as default \n");
431 }
432 strcpy(key, "acusim.license.scorec.rpi.edu");
433 }
434 #endif
435 return 1;
436 }
437