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