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