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