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) svLS_lhs 83 TYPE(svLS_lsType) 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 DO i=1, nshg 202 solinc(i,1:dof) = Res4(1:dof,i) 203 END DO 204 205c#################################################################### 206 ELSE 207c 208c.... lesSolve : main matrix solver 209c 210 lesId = numeqns(1) 211 eqnType = 1 212 213c call summary_start() 214 impistat=1 215 impistat2=1 216 tlescp1 = TMRC() 217#ifdef AMG 218 ! Initial Time Profiling 219 call cpu_time(cpusec(1)) 220 if (irun_amg_prec.gt.0) then 221 call ramg_control(colm,rowp,lhsK,lhsP, 222 & ilwork,BC,iBC,iper) 223 end if 224 225 call cpu_time(cpusec(6)) 226 if (irun_amg_prec.gt.0) then 227 ramg_flag = 1 228 if (irun_amg_prec.eq.2) then ! Some setup work (mode a) 229 ramg_window = 1.0 230 ramg_redo = 0 231 endif 232 do while (ramg_flag.le.irun_amg_prec) 233 ! if smart solve, possible run solve twice 234 ! restart only if meets plateau 235#endif 236 237c 238c.... setup the linear algebra solver 239c 240 rtmp = res(:,1:4) 241 call usrNew ( usr, eqnType, aperm, 242 & atemp, rtmp, solinc, 243 & flowDiag, sclrDiag, lesP, 244 & lesQ, iBC, BC, 245 & iper, ilwork, numpe, 246 & nshg, nshl, nPermDims, 247 & nTmpDims, rowp, colm, 248 & lhsK, lhsP, rdtmp, 249 & nnz_tot, CGsol ) 250c 251c.... solve linear system 252c 253 call myfLesSolve ( lesId, usr ) 254#ifdef AMG 255 ramg_flag = ramg_flag + 2 ! Default no second run in mode a 256 if (irun_amg_prec.eq.3) then 257 if (maxIters.gt.int(statsflow(4))) then 258 ramg_flag = ramg_flag + 1 ! Default no second run in mode b 259 endif 260 endif 261 enddo 262 else 263c 264c.... setup the linear algebra solver 265c 266 rtmp = res(:,1:4) 267 call usrNew ( usr, eqnType, aperm, 268 & atemp, rtmp, solinc, 269 & flowDiag, sclrDiag, lesP, 270 & lesQ, iBC, BC, 271 & iper, ilwork, numpe, 272 & nshg, nshl, nPermDims, 273 & nTmpDims, rowp, colm, 274 & lhsK, lhsP, rdtmp, 275 & nnz_tot, CGsol ) 276 277 call myfLesSolve( lesId, usr ) 278 endif 279 280 call cpu_time(cpusec(3)) 281 282 ramg_time(1) = ramg_time(1) + cpusec(3)-cpusec(1) 283 ramg_time(7) = ramg_time(7) + cpusec(6)-cpusec(1) 284 285 ! ramg_time: 1 : local total 286 ! 4 : local VG-cycle 287 ! 7 : local setup time 288 ! 11 : Ap-product level 1 289 ! 12 : Ap-product level >1 290 ! 13 : Prolongation/Restriction 291 ! 20 : local boundary MLS time 292 293 if (myrank.eq.master) then 294 write(*,*) 295 endif 296 call ramg_print_time(" == AMG == Total ACUSIM Solver:", 297 & ramg_time(1)) 298 call ramg_print_time(" == AMG == Setup: ",ramg_time(7)) 299 call ramg_print_time(" == AMG == Prec'd cycle: ",ramg_time(4)) 300 call ramg_print_time(" == AMG == Ap product(level=1): ", 301 & ramg_time(11)) 302 call ramg_print_time(" == AMG == Ap product(level>=2): ", 303 & ramg_time(12)) 304 call ramg_print_time(" == AMG == Pro/Restr ", 305 & ramg_time(13)) 306 call ramg_print_time(" == AMG == Boundary Ap (GS only)", 307 & ramg_time(20)) 308 if (myrank.eq.master) then 309 write(*,*) 310 endif 311 312#endif 313 314 ! End Time profiling output 315 316 call getSol ( usr, solinc ) 317 318 if (numpe > 1) then 319 call commu ( solinc, ilwork, nflow, 'out') 320 endif 321 tlescp2 = TMRC() 322 impistat=0 323 impistat2=0 324c call summary_stop() 325 326 tcorecp(1) = tcorecp(1) + telmcp2-telmcp1 ! elem. formation 327 tcorecp(2) = tcorecp(2) + tlescp2-tlescp1 ! linear alg. solution 328 ENDIF 329 call rstatic (res, y, solinc) ! output flow stats 330c 331c.... end 332c 333 return 334 end 335 336 subroutine SolSclr(y, ac, u, 337 & yold, acold, uold, 338 & x, iBC, 339 & BC, nPermDimsS, nTmpDimsS, 340 & apermS, atempS, iper, 341 & ilwork, shp, shgl, 342 & shpb, shglb, rowp, 343 & colm, lhsS, solinc, 344 & tcorecpscal) 345c 346c---------------------------------------------------------------------- 347c 348c This is the 2nd interface routine to the Farzin's linear equation 349c solver library. 350c 351c input: 352c y (nshg,ndof) : Y-variables at n+alpha_f 353c ac (nshg,ndof) : Primvar. accel. variable n+alpha_m 354c yold (nshg,ndof) : Y-variables at beginning of step 355c x (numnp,nsd) : node coordinates 356c iBC (nshg) : BC codes 357c BC (nshg,ndofBC) : BC constraint parameters 358c iper (nshg) : periodic nodal information 359c 360c output: 361c y (nshg,ndof) : Y-variables at n+alpha_f 362c ac (nshg,ndof) : Primvar. accel. variable n+alpha_m 363c 364c 365c The followings are preliminary steps required to use Farzin's 366c solver library. New way of writing has to be used such as 367c 368c | E | | dS | = | RScal | 369c 370c---------------------------------------------------------------------- 371c 372 use pointer_data 373 374 include "common.h" 375 include "mpif.h" 376 include "auxmpi.h" 377c 378 real*8 y(nshg,ndof), ac(nshg,ndof), 379 & yold(nshg,ndof), acold(nshg,ndof), 380 & u(nshg,nsd), uold(nshg,nsd), 381 & x(numnp,nsd), BC(nshg,ndofBC), 382 & res(nshg,1), 383 & flowDiag(nshg,4), 384 & sclrDiag(nshg,1), lhsS(nnz_tot), 385 & apermS(nshg,nPermDimsS), atempS(nshg,nTmpDimsS) 386 387c 388 real*8 shp(MAXTOP,maxsh,MAXQPT), 389 & shgl(MAXTOP,nsd,maxsh,MAXQPT), 390 & shpb(MAXTOP,maxsh,MAXQPT), 391 & shglb(MAXTOP,nsd,maxsh,MAXQPT) 392c 393 integer usr(100), eqnType, 394 & rowp(nshg*nnz), colm(nshg+1), 395 & iBC(nshg), ilwork(nlwork), 396 & iper(nshg) 397c 398 real*8 yAlpha(nshg,ndof), acAlpha(nshg,ndof), 399 & uAlpha(nshg,nsd), 400 & lesP(nshg,1), lesQ(nshg,1), 401 & solinc(nshg,1), CGsol(nshg), 402 & tcorecpscal(2) 403 404c 405c.... *******************>> Element Data Formation <<****************** 406c 407c.... compute solution at n+alpha 408c 409 call itrYAlpha( uold, yold, acold, 410 & u, y, ac, 411 & uAlpha, yAlpha, acAlpha) 412c 413c.... form the LHS matrices, the residual vector (at alpha) 414c 415 impistat=2 416 impistat2=1 417 telmcp1 = TMRC() 418 call ElmGMRSclr(yAlpha,acAlpha, x, 419 & shp, shgl, iBC, 420 & BC, shpb, shglb, 421 & res, iper, ilwork, 422 & rowp, colm, lhsS ) 423 telmcp2 = TMRC() 424 impistat=0 425 impistat2=0 426c 427c.... lesSolve : main matrix solver 428c 429 lesId = numeqns(1+nsolt+isclr) 430 eqnType = 2 431c 432c.... setup the linear algebra solver 433c 434 435 impistat=2 436 impistat2=1 437 tlescp1 = TMRC() 438 call usrNew ( usr, eqnType, apermS, 439 & atempS, res, solinc, 440 & flowDiag, sclrDiag, lesP, 441 & lesQ, iBC, BC, 442 & iper, ilwork, numpe, 443 & nshg, nshl, nPermDimsS, 444 & nTmpDimsS, rowp, colm, 445 & rlhsK, rlhsP, lhsS, 446 & nnz_tot, CGsol ) 447c 448c.... solve linear system 449c 450 call myfLesSolve ( lesId, usr ) 451 call getSol ( usr, solinc ) 452 453 if (numpe > 1) then 454 call commu ( solinc, ilwork, 1, 'out') 455 endif 456 tlescp2 = TMRC() 457 impistat=0 458 impistat2=0 459 460 tcorecpscal(1) = tcorecpscal(1) + telmcp2-telmcp1 ! elem. formation 461 tcorecpscal(2) = tcorecpscal(2) + tlescp2-tlescp1 ! linear alg. solution 462 463 nsolsc=5+isclr 464 call rstaticSclr (res, y, solinc, nsolsc) ! output scalar stats 465c 466c.... end 467c 468 return 469 end 470 471 472 473 474 475