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