1 subroutine SolFlow(y, ac, u, 2 & yold, acold, uold, 3 & x, iBC, 4 & BC, res, 5 & nPermDims, nTmpDims, aperm, 6 & atemp, iper, 7 & ilwork, shp, shgl, 8 & shpb, shglb, rowp, 9 & colm, lhsK, lhsP, 10 & solinc, rerr, tcorecp, 11 & GradV, sumtime 12#ifdef HAVE_SVLS 13 & ,svLS_lhs, svLS_ls, svLS_nFaces) 14#else 15 & ) 16#endif 17c 18c---------------------------------------------------------------------- 19c 20c This is the 2nd interface routine to the linear equation 21c solver library that uses the CGP and GMRES methods. 22c 23c input: 24c y (nshg,ndof) : Y-variables at n+alpha_f 25c ac (nshg,ndof) : Primvar. accel. variable n+alpha_m 26c yold (nshg,ndof) : Y-variables at beginning of step 27c acold (nshg,ndof) : Primvar. accel. at beginning of step 28c x (numnp,nsd) : node coordinates 29c iBC (nshg) : BC codes 30c BC (nshg,ndofBC) : BC constraint parameters 31c iper (nshg) : periodic nodal information 32c 33c output: 34c res (nshg,nflow) : preconditioned residual 35c y (nshg,ndof) : Y-variables at n+alpha_f 36c ac (nshg,ndof) : Primvar. accel. variable n+alpha_m 37c 38c 39c The followings are preliminary steps required to use Farzin's 40c solver library. New way of writing has to be used such as 41c 42c | K G | | du | | Rmom | 43c | | | | = | | 44c | G^t C | | dp | | Rcon | 45c 46c | E | | dT | = | Rtemp | 47c 48c where 49c 50c xKebe : K_ab = dRmom_a/du_b xTe : E_ab = dRtemp_a/dT_b 51c 52c G_ab = dRmom_a/dp_b 53c xGoC : 54c C_ab = dRcon_a/dp_b 55c 56c resf = Rmon Rcon rest = Rtemp 57c 58c 59c Zdenek Johan, Winter 1991. (Fortran 90) 60c Juin Kim, Summer 1998. (Incompressible flow solver interface) 61c Alberto Figueroa. CMM-FSI 62c---------------------------------------------------------------------- 63c 64 use pointer_data 65#ifdef AMG 66 use ramg_data 67#endif 68 69 include "common.h" 70 include "mpif.h" 71 include "auxmpi.h" 72#ifdef HAVE_SVLS 73 include "svLS.h" 74#endif 75c 76C 77C Argument variables 78C 79 INTEGER npermdims 80 INTEGER ntmpdims 81C 82C Local variables 83C 84 INTEGER lesid 85C 86 REAL*8 rdtmp 87C 88#ifdef HAVE_SVLS 89 TYPE(svLS_lhsType), INTENT(INOUT) :: svLS_lhs 90 TYPE(svLS_lsType), INTENT(INOUT) :: svLS_ls 91#endif 92 93 real*8 y(nshg,ndof), ac(nshg,ndof), 94 & yold(nshg,ndof), acold(nshg,ndof), 95 & u(nshg,nsd), uold(nshg,nsd), 96 & x(numnp,nsd), BC(nshg,ndofBC), 97 & res(nshg,nflow), tmpres(nshg,nflow), 98 & flowDiag(nshg,4), 99 & aperm(nshg,nPermDims), atemp(nshg,nTmpDims), 100 & sclrDiag(nshg,1), 101 & lhsK(9,nnz_tot), lhsP(4,nnz_tot), 102 & GradV(nshg,nsdsq) 103c 104 real*8 shp(MAXTOP,maxsh,MAXQPT), 105 & shgl(MAXTOP,nsd,maxsh,MAXQPT), 106 & shpb(MAXTOP,maxsh,MAXQPT), 107 & shglb(MAXTOP,nsd,maxsh,MAXQPT) 108c 109 integer usr(100), eqnType,temp, 110 & rowp(nshg*nnz), colm(nshg+1), 111 & iBC(nshg), ilwork(nlwork), 112 & iper(nshg) 113c 114 real*8 yAlpha(nshg,ndof), acAlpha(nshg,ndof), 115 & uAlpha(nshg,nsd), 116 & lesP(nshg,4), lesQ(nshg,4), 117 & solinc(nshg,ndof), CGsol(nshg) 118 119 real*8 tcorecp(2) 120 121 real*8 rerr(nshg,10), rtmp(nshg,4),rtemp 122 123 real*8 msum(4),mval(4),cpusec(10) 124 REAL*8 sumtime 125#ifdef HAVE_SVLS 126 INTEGER svLS_nFaces 127#endif 128 INTEGER dof, i, j, k, l 129 INTEGER, ALLOCATABLE :: incL(:) 130 REAL*8, ALLOCATABLE :: faceRes(:), Res4(:,:), Val4(:,:) 131 132c 133c.... *******************>> Element Data Formation <<****************** 134c 135c 136c.... set the parameters for flux and surface tension calculations 137c 138c 139 temp = npro 140 141 142 idflx = 0 143 if(idiff >= 1 ) idflx= (nflow-1) * nsd 144 if (isurf == 1) idflx=nflow*nsd 145c.... compute solution at n+alpha 146c 147 call itrYAlpha( uold, yold, acold, 148 & u, y, ac, 149 & uAlpha, yAlpha, acAlpha) 150 151c 152c.... form the LHS matrices, the residual vector (at alpha) 153c 154c call summary_start() 155 impistat=1 156 impistat2=1 157 telmcp1 = TMRC() 158 call ElmGMR (uAlpha, yAlpha, acAlpha, x, 159 & shp, shgl, iBC, 160 & BC, shpb, shglb, 161 & res, iper, ilwork, 162 & rowp, colm, lhsK, 163 & lhsP, rerr, GradV ) 164 telmcp2 = TMRC() 165 impistat=0 166 impistat2=0 167c call summary_stop() 168 169 170 tmpres(:,:) = res(:,:) 171 iblk = 1 172#ifdef HAVE_SVLS 173 IF (svLSFlag .EQ. 1) THEN 174 175c#################################################################### 176! Here calling svLS 177 178 ALLOCATE(faceRes(svLS_nFaces), incL(svLS_nFaces)) 179 faceRes=zero ! function to compute this left out at this time but would be needed to enable adnvanced p vs. Q BC's 180 incL = 1 181 dof = 4 182 IF (.NOT.ALLOCATED(Res4)) THEN 183 ALLOCATE (Res4(dof,nshg), Val4(dof*dof,nnz_tot)) 184 END IF 185 186 DO i=1, nshg 187 Res4(1:dof,i) = res(i,1:dof) 188 END DO 189 190 DO i=1, nnz_tot 191 Val4(1:3,i) = lhsK(1:3,i) 192 Val4(5:7,i) = lhsK(4:6,i) 193 Val4(9:11,i) = lhsK(7:9,i) 194 Val4(13:15,i) = lhsP(1:3,i) 195 Val4(16,i) = lhsP(4,i) 196 END DO 197 198 !Val4(4:12:4,:) = -lhsP(1:3,:)^t 199 DO i=1, nshg 200 Do j=colm(i), colm(i+1) - 1 201 k = rowp(j) 202 DO l=colm(k), colm(k+1) - 1 203 IF (rowp(l) .EQ. i) THEN 204 Val4(4:12:4,l) = -lhsP(1:3,j) 205 EXIT 206 END IF 207 END DO 208 END DO 209 END DO 210 CALL svLS_SOLVE(svLS_lhs, svLS_ls, dof, Res4, Val4, incL, 211 2 faceRes) 212 213 if(myrank.eq.master) write(*,*) 'svLS outer iterations', svLS_ls%RI%itr 214 statsflow(1)=1.0*svLS_ls%GM%itr 215 statsflow(4)=1.0*svLS_ls%CG%itr 216 DO i=1, nshg 217 solinc(i,1:dof) = Res4(1:dof,i) 218 END DO 219 ENDIF 220#endif 221 222#ifdef HAVE_LESLIB 223 if(leslib.eq.1) then 224c 225c.... lesSolve : main matrix solver 226c 227 lesId = numeqns(1) 228 eqnType = 1 229 230c call summary_start() 231 impistat=1 232 impistat2=1 233 tlescp1 = TMRC() 234#ifdef AMG 235 ! Initial Time Profiling 236 call cpu_time(cpusec(1)) 237 if (irun_amg_prec.gt.0) then 238 call ramg_control(colm,rowp,lhsK,lhsP, 239 & ilwork,BC,iBC,iper) 240 end if 241 242 call cpu_time(cpusec(6)) 243 if (irun_amg_prec.gt.0) then 244 ramg_flag = 1 245 if (irun_amg_prec.eq.2) then ! Some setup work (mode a) 246 ramg_window = 1.0 247 ramg_redo = 0 248 endif 249 do while (ramg_flag.le.irun_amg_prec) 250 ! if smart solve, possible run solve twice 251 ! restart only if meets plateau 252#endif 253 254c 255c.... setup the linear algebra solver 256c 257 rtmp = res(:,1:4) 258 call usrNew ( usr, eqnType, aperm, 259 & atemp, rtmp, solinc, 260 & flowDiag, sclrDiag, lesP, 261 & lesQ, iBC, BC, 262 & iper, ilwork, numpe, 263 & nshg, nshl, nPermDims, 264 & nTmpDims, rowp, colm, 265 & lhsK, lhsP, rdtmp, 266 & nnz_tot, CGsol ) 267c 268c.... solve linear system 269c 270 call myfLesSolve ( lesId, usr ) 271#ifdef AMG 272 ramg_flag = ramg_flag + 2 ! Default no second run in mode a 273 if (irun_amg_prec.eq.3) then 274 if (maxIters.gt.int(statsflow(4))) then 275 ramg_flag = ramg_flag + 1 ! Default no second run in mode b 276 endif 277 endif 278 enddo 279 else 280c 281c.... setup the linear algebra solver 282c 283 rtmp = res(:,1:4) 284 call usrNew ( usr, eqnType, aperm, 285 & atemp, rtmp, solinc, 286 & flowDiag, sclrDiag, lesP, 287 & lesQ, iBC, BC, 288 & iper, ilwork, numpe, 289 & nshg, nshl, nPermDims, 290 & nTmpDims, rowp, colm, 291 & lhsK, lhsP, rdtmp, 292 & nnz_tot, CGsol ) 293 294 call myfLesSolve( lesId, usr ) 295 endif 296 297 call cpu_time(cpusec(3)) 298 299 ramg_time(1) = ramg_time(1) + cpusec(3)-cpusec(1) 300 ramg_time(7) = ramg_time(7) + cpusec(6)-cpusec(1) 301 302 ! ramg_time: 1 : local total 303 ! 4 : local VG-cycle 304 ! 7 : local setup time 305 ! 11 : Ap-product level 1 306 ! 12 : Ap-product level >1 307 ! 13 : Prolongation/Restriction 308 ! 20 : local boundary MLS time 309 310 if (myrank.eq.master) then 311 write(*,*) 312 endif 313 call ramg_print_time(" == AMG == Total ACUSIM Solver:", 314 & ramg_time(1)) 315 call ramg_print_time(" == AMG == Setup: ",ramg_time(7)) 316 call ramg_print_time(" == AMG == Prec'd cycle: ",ramg_time(4)) 317 call ramg_print_time(" == AMG == Ap product(level=1): ", 318 & ramg_time(11)) 319 call ramg_print_time(" == AMG == Ap product(level>=2): ", 320 & ramg_time(12)) 321 call ramg_print_time(" == AMG == Pro/Restr ", 322 & ramg_time(13)) 323 call ramg_print_time(" == AMG == Boundary Ap (GS only)", 324 & ramg_time(20)) 325 if (myrank.eq.master) then 326 write(*,*) 327 endif 328 329#endif 330 331 ! End Time profiling output 332 333 call getSol ( usr, solinc ) 334 335 if (numpe > 1) then 336 call commu ( solinc, ilwork, nflow, 'out') 337 endif 338 ENDIF ! end of leslib flow solve 339#endif 340 tlescp2 = TMRC() 341 impistat=0 342 impistat2=0 343c call summary_stop() 344 345 tcorecp(1) = tcorecp(1) + telmcp2-telmcp1 ! elem. formation 346 tcorecp(2) = tcorecp(2) + tlescp2-tlescp1 ! linear alg. solution 347 call rstatic (res, y, solinc) ! output flow stats 348c 349c.... end 350c 351 return 352 end 353 354 subroutine SolSclr(y, ac, u, 355 & yold, acold, uold, 356 & x, iBC, 357 & BC, nPermDimsS, nTmpDimsS, 358 & apermS, atempS, iper, 359 & ilwork, shp, shgl, 360 & shpb, shglb, rowp, 361 & colm, lhsS, solinc, 362 & tcorecpscal 363#ifdef HAVE_SVLS 364 & ,svLS_lhs, svLS_ls, svLS_nFaces) 365#else 366 & ) 367#endif 368c 369c---------------------------------------------------------------------- 370c 371c This is the 2nd interface routine to the Farzin's linear equation 372c solver library. 373c 374c input: 375c y (nshg,ndof) : Y-variables at n+alpha_f 376c ac (nshg,ndof) : Primvar. accel. variable n+alpha_m 377c yold (nshg,ndof) : Y-variables at beginning of step 378c x (numnp,nsd) : node coordinates 379c iBC (nshg) : BC codes 380c BC (nshg,ndofBC) : BC constraint parameters 381c iper (nshg) : periodic nodal information 382c 383c output: 384c y (nshg,ndof) : Y-variables at n+alpha_f 385c ac (nshg,ndof) : Primvar. accel. variable n+alpha_m 386c 387c 388c The followings are preliminary steps required to use Farzin's 389c solver library. New way of writing has to be used such as 390c 391c | E | | dS | = | RScal | 392c 393c---------------------------------------------------------------------- 394c 395 use pointer_data 396 397 include "common.h" 398 include "mpif.h" 399 include "auxmpi.h" 400#ifdef HAVE_SVLS 401 include "svLS.h" 402#endif 403c 404 real*8 y(nshg,ndof), ac(nshg,ndof), 405 & yold(nshg,ndof), acold(nshg,ndof), 406 & u(nshg,nsd), uold(nshg,nsd), 407 & x(numnp,nsd), BC(nshg,ndofBC), 408 & res(nshg,1), 409 & flowDiag(nshg,4), 410 & sclrDiag(nshg,1), lhsS(nnz_tot), 411 & apermS(nshg,nPermDimsS), atempS(nshg,nTmpDimsS) 412 413c 414 real*8 shp(MAXTOP,maxsh,MAXQPT), 415 & shgl(MAXTOP,nsd,maxsh,MAXQPT), 416 & shpb(MAXTOP,maxsh,MAXQPT), 417 & shglb(MAXTOP,nsd,maxsh,MAXQPT) 418c 419 integer usr(100), eqnType, 420 & rowp(nshg*nnz), colm(nshg+1), 421 & iBC(nshg), ilwork(nlwork), 422 & iper(nshg) 423c 424 real*8 yAlpha(nshg,ndof), acAlpha(nshg,ndof), 425 & uAlpha(nshg,nsd), 426 & lesP(nshg,1), lesQ(nshg,1), 427 & solinc(nshg,1), CGsol(nshg), 428 & tcorecpscal(2) 429#ifdef HAVE_SVLS 430 TYPE(svLS_lhsType), INTENT(INOUT) :: svLS_lhs 431 TYPE(svLS_lsType), INTENT(INOUT) :: svLS_ls 432 INTEGER svLS_nFaces 433#endif 434 REAL*8 sumtime 435 INTEGER dof, i, j, k, l 436 INTEGER, ALLOCATABLE :: incL(:) 437 REAL*8, ALLOCATABLE :: faceRes(:), Res1(:,:), Val1(:,:) 438 439c 440c.... *******************>> Element Data Formation <<****************** 441c 442c.... compute solution at n+alpha 443c 444 call itrYAlpha( uold, yold, acold, 445 & u, y, ac, 446 & uAlpha, yAlpha, acAlpha) 447c 448c.... form the LHS matrices, the residual vector (at alpha) 449c 450 impistat=2 451 impistat2=1 452 telmcp1 = TMRC() 453 call ElmGMRSclr(yAlpha,acAlpha, x, 454 & shp, shgl, iBC, 455 & BC, shpb, shglb, 456 & res, iper, ilwork, 457 & rowp, colm, lhsS ) 458 telmcp2 = TMRC() 459 impistat=0 460 impistat2=0 461 statssclr(1)=0 462#ifdef HAVE_SVLS 463 IF (svLSFlag .EQ. 1) THEN 464 465c#################################################################### 466! Here calling svLS 467 468 ALLOCATE(faceRes(svLS_nFaces), incL(svLS_nFaces)) 469 faceRes=zero 470 incL = 1 471 dof = 1 472 IF (.NOT.ALLOCATED(Res1)) THEN 473 ALLOCATE (Res1(dof,nshg), Val1(dof*dof,nnz_tot)) 474 END IF 475 476 DO i=1, nshg 477 Res1(1,i) = res(i,1) 478 END DO 479 480 DO i=1, nnz_tot 481 Val1(1,i) = lhsS(i) 482 END DO 483 484 CALL svLS_SOLVE(svLS_lhs, svLS_ls, dof, Res1, Val1, incL, 485 2 faceRes) 486 statssclr(1)=1.0*svLS_ls%RI%itr 487 DO i=1, nshg 488 solinc(i,1) = Res1(1,i) 489 END DO 490 ENDIF 491#endif 492#ifdef HAVE_LESLIB 493 if(leslib.eq.1) then 494c 495c.... lesSolve : main matrix solver 496c 497 lesId = numeqns(1+nsolt+isclr) 498 eqnType = 2 499c 500c.... setup the linear algebra solver 501c 502 503 impistat=2 504 impistat2=1 505 tlescp1 = TMRC() 506 call usrNew ( usr, eqnType, apermS, 507 & atempS, res, solinc, 508 & flowDiag, sclrDiag, lesP, 509 & lesQ, iBC, BC, 510 & iper, ilwork, numpe, 511 & nshg, nshl, nPermDimsS, 512 & nTmpDimsS, rowp, colm, 513 & rlhsK, rlhsP, lhsS, 514 & nnz_tot, CGsol ) 515c 516c.... solve linear system 517c 518 call myfLesSolve ( lesId, usr ) 519 call getSol ( usr, solinc ) 520 521 if (numpe > 1) then 522 call commu ( solinc, ilwork, 1, 'out') 523 endif 524 ENDIF ! leslib conditional 525#endif 526 tlescp2 = TMRC() 527 impistat=0 528 impistat2=0 529 530 tcorecpscal(1) = tcorecpscal(1) + telmcp2-telmcp1 ! elem. formation 531 tcorecpscal(2) = tcorecpscal(2) + tlescp2-tlescp1 ! linear alg. solution 532 533 nsolsc=5+isclr 534 call rstaticSclr (res, y, solinc, nsolsc) ! output scalar stats 535c 536c.... end 537c 538 return 539 end 540 541 542 543 544 545