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#if defined(HAVE_SVLS) && defined(HAVE_ACUSOLVE) 221c#################################################################### 222 ELSE 223#elif defined(HAVE_SVLS) 224 ENDIF 225#endif 226#ifdef HAVE_ACUSOLVE 227c 228c.... lesSolve : main matrix solver 229c 230 lesId = numeqns(1) 231 eqnType = 1 232 233c call summary_start() 234 impistat=1 235 impistat2=1 236 tlescp1 = TMRC() 237#ifdef AMG 238 ! Initial Time Profiling 239 call cpu_time(cpusec(1)) 240 if (irun_amg_prec.gt.0) then 241 call ramg_control(colm,rowp,lhsK,lhsP, 242 & ilwork,BC,iBC,iper) 243 end if 244 245 call cpu_time(cpusec(6)) 246 if (irun_amg_prec.gt.0) then 247 ramg_flag = 1 248 if (irun_amg_prec.eq.2) then ! Some setup work (mode a) 249 ramg_window = 1.0 250 ramg_redo = 0 251 endif 252 do while (ramg_flag.le.irun_amg_prec) 253 ! if smart solve, possible run solve twice 254 ! restart only if meets plateau 255#endif 256 257c 258c.... setup the linear algebra solver 259c 260 rtmp = res(:,1:4) 261 call usrNew ( usr, eqnType, aperm, 262 & atemp, rtmp, solinc, 263 & flowDiag, sclrDiag, lesP, 264 & lesQ, iBC, BC, 265 & iper, ilwork, numpe, 266 & nshg, nshl, nPermDims, 267 & nTmpDims, rowp, colm, 268 & lhsK, lhsP, rdtmp, 269 & nnz_tot, CGsol ) 270c 271c.... solve linear system 272c 273 call myfLesSolve ( lesId, usr ) 274#ifdef AMG 275 ramg_flag = ramg_flag + 2 ! Default no second run in mode a 276 if (irun_amg_prec.eq.3) then 277 if (maxIters.gt.int(statsflow(4))) then 278 ramg_flag = ramg_flag + 1 ! Default no second run in mode b 279 endif 280 endif 281 enddo 282 else 283c 284c.... setup the linear algebra solver 285c 286 rtmp = res(:,1:4) 287 call usrNew ( usr, eqnType, aperm, 288 & atemp, rtmp, solinc, 289 & flowDiag, sclrDiag, lesP, 290 & lesQ, iBC, BC, 291 & iper, ilwork, numpe, 292 & nshg, nshl, nPermDims, 293 & nTmpDims, rowp, colm, 294 & lhsK, lhsP, rdtmp, 295 & nnz_tot, CGsol ) 296 297 call myfLesSolve( lesId, usr ) 298 endif 299 300 call cpu_time(cpusec(3)) 301 302 ramg_time(1) = ramg_time(1) + cpusec(3)-cpusec(1) 303 ramg_time(7) = ramg_time(7) + cpusec(6)-cpusec(1) 304 305 ! ramg_time: 1 : local total 306 ! 4 : local VG-cycle 307 ! 7 : local setup time 308 ! 11 : Ap-product level 1 309 ! 12 : Ap-product level >1 310 ! 13 : Prolongation/Restriction 311 ! 20 : local boundary MLS time 312 313 if (myrank.eq.master) then 314 write(*,*) 315 endif 316 call ramg_print_time(" == AMG == Total ACUSIM Solver:", 317 & ramg_time(1)) 318 call ramg_print_time(" == AMG == Setup: ",ramg_time(7)) 319 call ramg_print_time(" == AMG == Prec'd cycle: ",ramg_time(4)) 320 call ramg_print_time(" == AMG == Ap product(level=1): ", 321 & ramg_time(11)) 322 call ramg_print_time(" == AMG == Ap product(level>=2): ", 323 & ramg_time(12)) 324 call ramg_print_time(" == AMG == Pro/Restr ", 325 & ramg_time(13)) 326 call ramg_print_time(" == AMG == Boundary Ap (GS only)", 327 & ramg_time(20)) 328 if (myrank.eq.master) then 329 write(*,*) 330 endif 331 332#endif 333 334 ! End Time profiling output 335 336 call getSol ( usr, solinc ) 337 338 if (numpe > 1) then 339 call commu ( solinc, ilwork, nflow, 'out') 340 endif 341#endif 342#if defined(HAVE_SVLS) && defined(HAVE_ACUSOLVE) 343 ENDIF ! end of selection between solvers. 344#endif 345 tlescp2 = TMRC() 346 impistat=0 347 impistat2=0 348c call summary_stop() 349 350 tcorecp(1) = tcorecp(1) + telmcp2-telmcp1 ! elem. formation 351 tcorecp(2) = tcorecp(2) + tlescp2-tlescp1 ! linear alg. solution 352 call rstatic (res, y, solinc) ! output flow stats 353c 354c.... end 355c 356 return 357 end 358 359 subroutine SolSclr(y, ac, u, 360 & yold, acold, uold, 361 & x, iBC, 362 & BC, nPermDimsS, nTmpDimsS, 363 & apermS, atempS, iper, 364 & ilwork, shp, shgl, 365 & shpb, shglb, rowp, 366 & colm, lhsS, solinc, 367 & tcorecpscal 368#ifdef HAVE_SVLS 369 & ,svLS_lhs, svLS_ls, svLS_nFaces) 370#else 371 & ) 372#endif 373c 374c---------------------------------------------------------------------- 375c 376c This is the 2nd interface routine to the Farzin's linear equation 377c solver library. 378c 379c input: 380c y (nshg,ndof) : Y-variables at n+alpha_f 381c ac (nshg,ndof) : Primvar. accel. variable n+alpha_m 382c yold (nshg,ndof) : Y-variables at beginning of step 383c x (numnp,nsd) : node coordinates 384c iBC (nshg) : BC codes 385c BC (nshg,ndofBC) : BC constraint parameters 386c iper (nshg) : periodic nodal information 387c 388c output: 389c y (nshg,ndof) : Y-variables at n+alpha_f 390c ac (nshg,ndof) : Primvar. accel. variable n+alpha_m 391c 392c 393c The followings are preliminary steps required to use Farzin's 394c solver library. New way of writing has to be used such as 395c 396c | E | | dS | = | RScal | 397c 398c---------------------------------------------------------------------- 399c 400 use pointer_data 401 402 include "common.h" 403 include "mpif.h" 404 include "auxmpi.h" 405#ifdef HAVE_SVLS 406 include "svLS.h" 407#endif 408c 409 real*8 y(nshg,ndof), ac(nshg,ndof), 410 & yold(nshg,ndof), acold(nshg,ndof), 411 & u(nshg,nsd), uold(nshg,nsd), 412 & x(numnp,nsd), BC(nshg,ndofBC), 413 & res(nshg,1), 414 & flowDiag(nshg,4), 415 & sclrDiag(nshg,1), lhsS(nnz_tot), 416 & apermS(nshg,nPermDimsS), atempS(nshg,nTmpDimsS) 417 418c 419 real*8 shp(MAXTOP,maxsh,MAXQPT), 420 & shgl(MAXTOP,nsd,maxsh,MAXQPT), 421 & shpb(MAXTOP,maxsh,MAXQPT), 422 & shglb(MAXTOP,nsd,maxsh,MAXQPT) 423c 424 integer usr(100), eqnType, 425 & rowp(nshg*nnz), colm(nshg+1), 426 & iBC(nshg), ilwork(nlwork), 427 & iper(nshg) 428c 429 real*8 yAlpha(nshg,ndof), acAlpha(nshg,ndof), 430 & uAlpha(nshg,nsd), 431 & lesP(nshg,1), lesQ(nshg,1), 432 & solinc(nshg,1), CGsol(nshg), 433 & tcorecpscal(2) 434#ifdef HAVE_SVLS 435 TYPE(svLS_lhsType), INTENT(INOUT) :: svLS_lhs 436 TYPE(svLS_lsType), INTENT(INOUT) :: svLS_ls 437 INTEGER svLS_nFaces 438#endif 439 REAL*8 sumtime 440 INTEGER dof, i, j, k, l 441 INTEGER, ALLOCATABLE :: incL(:) 442 REAL*8, ALLOCATABLE :: faceRes(:), Res1(:,:), Val1(:,:) 443 444c 445c.... *******************>> Element Data Formation <<****************** 446c 447c.... compute solution at n+alpha 448c 449 call itrYAlpha( uold, yold, acold, 450 & u, y, ac, 451 & uAlpha, yAlpha, acAlpha) 452c 453c.... form the LHS matrices, the residual vector (at alpha) 454c 455 impistat=2 456 impistat2=1 457 telmcp1 = TMRC() 458 call ElmGMRSclr(yAlpha,acAlpha, x, 459 & shp, shgl, iBC, 460 & BC, shpb, shglb, 461 & res, iper, ilwork, 462 & rowp, colm, lhsS ) 463 telmcp2 = TMRC() 464 impistat=0 465 impistat2=0 466 statssclr(1)=0 467#ifdef HAVE_SVLS 468 IF (svLSFlag .EQ. 1) THEN 469 470c#################################################################### 471! Here calling svLS 472 473 ALLOCATE(faceRes(svLS_nFaces), incL(svLS_nFaces)) 474 faceRes=zero 475 incL = 1 476 dof = 1 477 IF (.NOT.ALLOCATED(Res1)) THEN 478 ALLOCATE (Res1(dof,nshg), Val1(dof*dof,nnz_tot)) 479 END IF 480 481 DO i=1, nshg 482 Res1(1,i) = res(i,1) 483 END DO 484 485 DO i=1, nnz_tot 486 Val1(1,i) = lhsS(i) 487 END DO 488 489 CALL svLS_SOLVE(svLS_lhs, svLS_ls, dof, Res1, Val1, incL, 490 2 faceRes) 491 statssclr(1)=1.0*svLS_ls%RI%itr 492 DO i=1, nshg 493 solinc(i,1) = Res1(1,i) 494 END DO 495#endif 496#if defined(HAVE_SVLS) && defined(HAVE_ACUSOLVE) 497c#################################################################### 498 ELSE 499#elif defined(HAVE_SVLS) 500 ENDIF 501#endif 502#ifdef HAVE_ACUSOLVE 503c 504c.... lesSolve : main matrix solver 505c 506 lesId = numeqns(1+nsolt+isclr) 507 eqnType = 2 508c 509c.... setup the linear algebra solver 510c 511 512 impistat=2 513 impistat2=1 514 tlescp1 = TMRC() 515 call usrNew ( usr, eqnType, apermS, 516 & atempS, res, solinc, 517 & flowDiag, sclrDiag, lesP, 518 & lesQ, iBC, BC, 519 & iper, ilwork, numpe, 520 & nshg, nshl, nPermDimsS, 521 & nTmpDimsS, rowp, colm, 522 & rlhsK, rlhsP, lhsS, 523 & nnz_tot, CGsol ) 524c 525c.... solve linear system 526c 527 call myfLesSolve ( lesId, usr ) 528 call getSol ( usr, solinc ) 529 530 if (numpe > 1) then 531 call commu ( solinc, ilwork, 1, 'out') 532 endif 533#endif 534#if defined(HAVE_SVLS) && defined(HAVE_ACUSOLVE) 535 ENDIF ! decision between solvers 536#endif 537 tlescp2 = TMRC() 538 impistat=0 539 impistat2=0 540 541 tcorecpscal(1) = tcorecpscal(1) + telmcp2-telmcp1 ! elem. formation 542 tcorecpscal(2) = tcorecpscal(2) + tlescp2-tlescp1 ! linear alg. solution 543 544 nsolsc=5+isclr 545 call rstaticSclr (res, y, solinc, nsolsc) ! output scalar stats 546c 547c.... end 548c 549 return 550 end 551 552 553 554 555 556