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 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 tlescp2 = TMRC() 325 impistat=0 326 impistat2=0 327c call summary_stop() 328 329 tcorecp(1) = tcorecp(1) + telmcp2-telmcp1 ! elem. formation 330 tcorecp(2) = tcorecp(2) + tlescp2-tlescp1 ! linear alg. solution 331 ENDIF 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) 348c 349c---------------------------------------------------------------------- 350c 351c This is the 2nd interface routine to the Farzin's linear equation 352c solver library. 353c 354c input: 355c y (nshg,ndof) : Y-variables at n+alpha_f 356c ac (nshg,ndof) : Primvar. accel. variable n+alpha_m 357c yold (nshg,ndof) : Y-variables at beginning of step 358c x (numnp,nsd) : node coordinates 359c iBC (nshg) : BC codes 360c BC (nshg,ndofBC) : BC constraint parameters 361c iper (nshg) : periodic nodal information 362c 363c output: 364c y (nshg,ndof) : Y-variables at n+alpha_f 365c ac (nshg,ndof) : Primvar. accel. variable n+alpha_m 366c 367c 368c The followings are preliminary steps required to use Farzin's 369c solver library. New way of writing has to be used such as 370c 371c | E | | dS | = | RScal | 372c 373c---------------------------------------------------------------------- 374c 375 use pointer_data 376 377 include "common.h" 378 include "mpif.h" 379 include "auxmpi.h" 380c 381 real*8 y(nshg,ndof), ac(nshg,ndof), 382 & yold(nshg,ndof), acold(nshg,ndof), 383 & u(nshg,nsd), uold(nshg,nsd), 384 & x(numnp,nsd), BC(nshg,ndofBC), 385 & res(nshg,1), 386 & flowDiag(nshg,4), 387 & sclrDiag(nshg,1), lhsS(nnz_tot), 388 & apermS(nshg,nPermDimsS), atempS(nshg,nTmpDimsS) 389 390c 391 real*8 shp(MAXTOP,maxsh,MAXQPT), 392 & shgl(MAXTOP,nsd,maxsh,MAXQPT), 393 & shpb(MAXTOP,maxsh,MAXQPT), 394 & shglb(MAXTOP,nsd,maxsh,MAXQPT) 395c 396 integer usr(100), eqnType, 397 & rowp(nshg*nnz), colm(nshg+1), 398 & iBC(nshg), ilwork(nlwork), 399 & iper(nshg) 400c 401 real*8 yAlpha(nshg,ndof), acAlpha(nshg,ndof), 402 & uAlpha(nshg,nsd), 403 & lesP(nshg,1), lesQ(nshg,1), 404 & solinc(nshg,1), CGsol(nshg), 405 & tcorecpscal(2) 406 407c 408c.... *******************>> Element Data Formation <<****************** 409c 410c.... compute solution at n+alpha 411c 412 call itrYAlpha( uold, yold, acold, 413 & u, y, ac, 414 & uAlpha, yAlpha, acAlpha) 415c 416c.... form the LHS matrices, the residual vector (at alpha) 417c 418 impistat=2 419 impistat2=1 420 telmcp1 = TMRC() 421 call ElmGMRSclr(yAlpha,acAlpha, x, 422 & shp, shgl, iBC, 423 & BC, shpb, shglb, 424 & res, iper, ilwork, 425 & rowp, colm, lhsS ) 426 telmcp2 = TMRC() 427 impistat=0 428 impistat2=0 429c 430c.... lesSolve : main matrix solver 431c 432 lesId = numeqns(1+nsolt+isclr) 433 eqnType = 2 434c 435c.... setup the linear algebra solver 436c 437 438 impistat=2 439 impistat2=1 440 tlescp1 = TMRC() 441 call usrNew ( usr, eqnType, apermS, 442 & atempS, res, solinc, 443 & flowDiag, sclrDiag, lesP, 444 & lesQ, iBC, BC, 445 & iper, ilwork, numpe, 446 & nshg, nshl, nPermDimsS, 447 & nTmpDimsS, rowp, colm, 448 & rlhsK, rlhsP, lhsS, 449 & nnz_tot, CGsol ) 450c 451c.... solve linear system 452c 453 call myfLesSolve ( lesId, usr ) 454 call getSol ( usr, solinc ) 455 456 if (numpe > 1) then 457 call commu ( solinc, ilwork, 1, 'out') 458 endif 459 tlescp2 = TMRC() 460 impistat=0 461 impistat2=0 462 463 tcorecpscal(1) = tcorecpscal(1) + telmcp2-telmcp1 ! elem. formation 464 tcorecpscal(2) = tcorecpscal(2) + tlescp2-tlescp1 ! linear alg. solution 465 466 nsolsc=5+isclr 467 call rstaticSclr (res, y, solinc, nsolsc) ! output scalar stats 468c 469c.... end 470c 471 return 472 end 473 474 475 476 477 478