1c 2cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 3ccccccccccccc SPARSE 4c_______________________________________________________________ 5 6 subroutine ElmGMRPETSc (y, ac, x, 7 & shp, shgl, iBC, 8 & BC, shpb, shglb, 9 & res, rmes, 10 & iper, ilwork, 11 & rerr, lhsP) 12c 13c---------------------------------------------------------------------- 14c 15c This routine computes the LHS mass matrix, the RHS residual 16c vector, for use with the GMRES 17c solver. 18c 19c Zdenek Johan, Winter 1991. (Fortran 90) 20c Chris Whiting, Winter 1998. (Matrix EBE-GMRES) 21c---------------------------------------------------------------------- 22c 23 use pointer_data 24 use timedataC 25c 26 include "common.h" 27 include "mpif.h" 28c 29! integer col(nshg+1), row(nnz*nshg) 30! real*8 lhsK(nflow*nflow,nnz_tot) 31 real*8 BDiag(1,1,1) 32 33 dimension y(nshg,ndof), ac(nshg,ndof), 34 & x(numnp,nsd), 35 & iBC(nshg), 36 & BC(nshg,ndofBC), 37 & res(nshg,nflow), 38 & rmes(nshg,nflow), 39 & iper(nshg) 40c 41 dimension shp(MAXTOP,maxsh,MAXQPT), 42 & shgl(MAXTOP,nsd,maxsh,MAXQPT), 43 & shpb(MAXTOP,maxsh,MAXQPT), 44 & shglb(MAXTOP,nsd,maxsh,MAXQPT) 45c 46 dimension qres(nshg, idflx), rmass(nshg) 47c 48 dimension ilwork(nlwork) 49 50 real*8 Bdiagvec(nshg,nflow), rerr(nshg,10) 51 52 real*8, allocatable :: tmpshp(:,:), tmpshgl(:,:,:) 53 real*8, allocatable :: tmpshpb(:,:), tmpshglb(:,:,:) 54 real*8, allocatable :: EGmass(:,:,:) 55 integer gnode(numnp) 56 ttim(80) = ttim(80) - secs(0.0) 57 iprec=0 58c 59c.... set up the timer 60c 61! call MPI_BARRIER (MPI_COMM_WORLD,ierr) 62! if(myrank.eq.0) write (*,*) 'top of elmgmrpetsc ' 63 64 call timer ('Elm_Form') 65c 66c.... --------------------> interior elements <-------------------- 67c 68c.... set up parameters 69c 70 ires = 1 71c 72 if (idiff==1 .or. idiff==3 .or. isurf==1) then ! global reconstruction 73 ! of qdiff 74c 75c loop over element blocks for the global reconstruction 76c of the diffusive flux vector, q, and lumped mass matrix, rmass 77c 78 qres = zero 79 rmass = zero 80 81 do iblk = 1, nelblk 82c 83c.... set up the parameters 84c 85 nenl = lcblk(5,iblk) ! no. of vertices per element 86 iel = lcblk(1,iblk) 87 lelCat = lcblk(2,iblk) 88 lcsyst = lcblk(3,iblk) 89 iorder = lcblk(4,iblk) 90 nenl = lcblk(5,iblk) ! no. of vertices per element 91 nshl = lcblk(10,iblk) 92 mattyp = lcblk(7,iblk) 93 ndofl = lcblk(8,iblk) 94 nsymdl = lcblk(9,iblk) 95 npro = lcblk(1,iblk+1) - iel 96 ngauss = nint(lcsyst) 97c 98c.... compute and assemble diffusive flux vector residual, qres, 99c and lumped mass matrix, rmass 100 101 allocate (tmpshp(nshl,MAXQPT)) 102 allocate (tmpshgl(nsd,nshl,MAXQPT)) 103 104 tmpshp(1:nshl,:) = shp(lcsyst,1:nshl,:) 105 tmpshgl(:,1:nshl,:) = shgl(lcsyst,:,1:nshl,:) 106 107 call AsIq (y, x, 108 & tmpshp, 109 & tmpshgl, 110 & mien(iblk)%p, mxmudmi(iblk)%p, 111 & qres, 112 & rmass) 113 114 deallocate ( tmpshp ) 115 deallocate ( tmpshgl ) 116 enddo 117 118c 119c.... take care of periodic boundary conditions 120c 121 122 call qpbc( rmass, qres, iBC, iper, ilwork ) 123c 124 endif ! computation of global diffusive flux 125c 126c.... loop over element blocks to compute element residuals 127c 128c 129c.... initialize the arrays 130c 131 res = zero 132 rmes = zero ! to avoid trap_uninitialized 133 !if (lhs. eq. 1) lhsK = zero 134 flxID = zero 135c 136c.... loop over the element-blocks 137c 138! call commu (res , ilwork, nflow , 'in ') !FOR TEST 139! call MPI_BARRIER (MPI_COMM_WORLD,ierr) 140! if(myrank.eq.0) write (*,*) 'after res zeroed ' 141 do iblk = 1, nelblk 142c 143c.... set up the parameters 144c 145 iblkts = iblk ! used in timeseries 146 nenl = lcblk(5,iblk) ! no. of vertices per element 147 iel = lcblk(1,iblk) 148 lelCat = lcblk(2,iblk) 149 lcsyst = lcblk(3,iblk) 150 iorder = lcblk(4,iblk) 151 nenl = lcblk(5,iblk) ! no. of vertices per element 152 nshl = lcblk(10,iblk) 153 mattyp = lcblk(7,iblk) 154 ndofl = lcblk(8,iblk) 155 nsymdl = lcblk(9,iblk) 156 npro = lcblk(1,iblk+1) - iel 157 inum = iel + npro - 1 158 ngauss = nint(lcsyst) 159c 160c.... compute and assemble the residual and tangent matrix 161c 162 163 if(lhs.eq.1) then 164 allocate (EGmass(npro,nedof,nedof)) 165 EGmass = zero 166 else 167 allocate (EGmass(1,1,1)) 168 endif 169 170 allocate (tmpshp(nshl,MAXQPT)) 171 allocate (tmpshgl(nsd,nshl,MAXQPT)) 172 tmpshp(1:nshl,:) = shp(lcsyst,1:nshl,:) 173 tmpshgl(:,1:nshl,:) = shgl(lcsyst,:,1:nshl,:) 174 175 call AsIGMR (y, ac, 176 & x, mxmudmi(iblk)%p, 177 & tmpshp, 178 & tmpshgl, mien(iblk)%p, 179 & mmat(iblk)%p, res, 180 & rmes, BDiag, 181 & qres, EGmass, 182 & rerr ) 183 if(lhs.eq.1) then 184c 185c.... satisfy the BCs on the implicit LHS 186c 187 call bc3LHS (iBC, BC, mien(iblk)%p, 188 & EGmass ) 189 190c 191c.... Fill-up the global sparse LHS mass matrix 192c 193! if(myrank.eq.0) write (*,*) 'before fillsparsepetscc ',iblk 194 call cycle_count_start() 195 call fillsparsecpetscc( mieng(iblk)%p, EGmass, lhsP) 196 call cycle_count_stop() 197 endif 198c 199 deallocate ( EGmass ) 200 deallocate ( tmpshp ) 201 deallocate ( tmpshgl ) 202c 203c.... end of interior element loop 204c 205 enddo 206 if(lhs.eq.1) call cycle_count_print() 207! call commu (res , ilwork, nflow , 'in ') !FOR TEST 208! call MPI_BARRIER (MPI_COMM_WORLD,ierr) 209! if(myrank.eq.0) write (*,*) 'after fillsparsepetscc ' 210c 211c.... --------------------> boundary elements <-------------------- 212c 213c.... loop over the boundary elements 214c 215 do iblk = 1, nelblb 216c 217c.... set up the parameters 218c 219 iel = lcblkb(1,iblk) 220 lelCat = lcblkb(2,iblk) 221 lcsyst = lcblkb(3,iblk) 222 iorder = lcblkb(4,iblk) 223 nenl = lcblkb(5,iblk) ! no. of vertices per element 224 nenbl = lcblkb(6,iblk) ! no. of vertices per bdry. face 225 mattyp = lcblkb(7,iblk) 226 ndofl = lcblkb(8,iblk) 227 nshl = lcblkb(9,iblk) 228 nshlb = lcblkb(10,iblk) 229 npro = lcblkb(1,iblk+1) - iel 230 if(lcsyst.eq.3) lcsyst=nenbl 231 ngaussb = nintb(lcsyst) 232 233c 234c.... compute and assemble the residuals corresponding to the 235c boundary integral 236c 237 238 allocate (tmpshpb(nshl,MAXQPT)) 239 allocate (tmpshglb(nsd,nshl,MAXQPT)) 240 241 tmpshpb(1:nshl,:) = shpb(lcsyst,1:nshl,:) 242 tmpshglb(:,1:nshl,:) = shglb(lcsyst,:,1:nshl,:) 243 244 call AsBMFG (y, x, 245 & tmpshpb, tmpshglb, 246 & mienb(iblk)%p, mmatb(iblk)%p, 247 & miBCB(iblk)%p, mBCB(iblk)%p, 248 & res, rmes) 249 250 deallocate (tmpshpb) 251 deallocate (tmpshglb) 252c 253c.... end of boundary element loop 254c 255 enddo 256c 257 ttim(80) = ttim(80) + secs(0.0) 258c 259c before the commu we need to rotate the residual vector for axisymmetric 260c boundary conditions (so that off processor periodicity is a dof add instead 261c of a dof combination). Take care of all nodes now so periodicity, like 262c commu is a simple dof add. 263c 264 if(iabc==1) then !are there any axisym bc's 265 call rotabc(res(1,2), iBC, 'in ') 266 endif 267 268c.... --------------------> communications <------------------------- 269c 270 271 if (numpe > 1) then 272 call commu (res , ilwork, nflow , 'in ') 273 274 call MPI_BARRIER (MPI_COMM_WORLD,ierr) 275 endif 276 277c 278c.... ----------------------> post processing <---------------------- 279c 280c.... satisfy the BCs on the residual 281c 282 call bc3Res (y, iBC, BC, res, iper, ilwork) 283c 284c 285c.... return 286c 287 if (numpe > 1) then 288 call commu (res , ilwork, nflow , 'out') 289 call MPI_BARRIER (MPI_COMM_WORLD,ierr) 290 endif 291 call timer ('Back ') 292 return 293 end 294c 295c 296 297c 298 299 subroutine ElmGMRPETScSclr(y, ac, 300 & x, elDw, 301 & shp, shgl, iBC, 302 & BC, shpb, shglb, 303 & rest, rmest, 304 & iper, ilwork, lhsPs) 305c 306c---------------------------------------------------------------------- 307c 308c This routine computes the LHS mass matrix, the RHS residual 309c vector, and the preconditioning matrix, for use with the GMRES 310c solver. 311c 312c Zdenek Johan, Winter 1991. (Fortran 90) 313c Chris Whiting, Winter 1998. (Matrix EBE-GMRES) 314c---------------------------------------------------------------------- 315c 316 use pointer_data 317c 318 include "common.h" 319 include "mpif.h" 320c 321 dimension y(nshg,ndof), ac(nshg,ndof), 322 & x(numnp,nsd), 323 & iBC(nshg), 324 & BC(nshg,ndofBC), 325 & rest(nshg), Diag(1), 326 & rmest(nshg), 327 & iper(nshg) 328c 329 dimension shp(MAXTOP,maxsh,MAXQPT), 330 & shgl(MAXTOP,nsd,maxsh,MAXQPT), 331 & shpb(MAXTOP,maxsh,MAXQPT), 332 & shglb(MAXTOP,nsd,maxsh,MAXQPT) 333c 334 dimension qrest(nshg), rmasst(nshg) 335 real*8 elDw(numel) 336c 337 dimension ilwork(nlwork) 338c 339 real*8, allocatable :: tmpshp(:,:), tmpshgl(:,:,:) 340 real*8, allocatable :: tmpshpb(:,:), tmpshglb(:,:,:) 341 real*8, allocatable :: EGMasst(:,:,:) 342 real*8, allocatable :: ElDwl(:) 343c 344 ttim(80) = ttim(80) - tmr() 345 iprec=0 346c 347c.... set up the timer 348c 349 350 call timer ('Elm_Form') 351c 352c.... --------------------> interior elements <-------------------- 353c 354c.... set up parameters 355c 356 intrul = intg (1,itseq) 357 intind = intpt (intrul) 358c 359 ires = 1 360 361c 362c.... initialize the arrays 363c 364 rest = zero 365 rmest = zero ! to avoid trap_uninitialized 366c 367c.... loop over the element-blocks 368c 369 do iblk = 1, nelblk 370c 371c 372 nenl = lcblk(5,iblk) ! no. of vertices per element 373 iel = lcblk(1,iblk) 374 lelCat = lcblk(2,iblk) 375 lcsyst = lcblk(3,iblk) 376 iorder = lcblk(4,iblk) 377 nshl = lcblk(10,iblk) 378 mattyp = lcblk(7,iblk) 379 ndofl = lcblk(8,iblk) 380 nsymdl = lcblk(9,iblk) 381 npro = lcblk(1,iblk+1) - iel 382 inum = iel + npro - 1 383 ngauss = nint(lcsyst) 384c 385c.... compute and assemble the residual and tangent matrix 386c 387 allocate (tmpshp(nshl,MAXQPT)) 388 allocate (tmpshgl(nsd,nshl,MAXQPT)) 389 if(lhs.eq.1) then 390 allocate (EGMasst(npro,nshl,nshl)) 391 EGMasst=zero 392 endif 393 394 tmpshp(1:nshl,:) = shp(lcsyst,1:nshl,:) 395 tmpshgl(:,1:nshl,:) = shgl(lcsyst,:,1:nshl,:) 396 397 allocate (elDwl(npro)) 398c 399 call AsIGMRSclr(y, 400 & ac, 401 & x, elDwl, 402 & tmpshp, tmpshgl, 403 & mien(iblk)%p, 404 & mmat(iblk)%p, rest, 405 & rmest, 406 & qrest, EGmasst, 407 & Diag ) 408c 409 elDw(iel:inum) = elDwl(1:npro) 410 deallocate ( elDwl ) 411 412 if(lhs.eq.1) then 413c.... satisfy the BCs on the implicit LHS 414c 415 call bc3LHSSclr (iBC, mien(iblk)%p, EGmasst ) 416c 417c 418c.... Fill-up the global sparse LHS mass matrix 419c 420 call fillsparsecpetscs( mieng(iblk)%p, EGmasst, lhsPs) 421 endif 422 if(lhs.eq.1) deallocate ( EGmasst ) 423 deallocate ( tmpshp ) 424 deallocate ( tmpshgl ) 425c.... end of interior element loop 426c 427 enddo 428c 429c.... --------------------> boundary elements <-------------------- 430c 431c.... set up parameters 432c 433 intrul = intg (2,itseq) 434 intind = intptb (intrul) 435c 436c.... loop over the boundary elements 437c 438 do iblk = 1, nelblb 439c 440c.... set up the parameters 441c 442 iel = lcblkb(1,iblk) 443 lelCat = lcblkb(2,iblk) 444 lcsyst = lcblkb(3,iblk) 445 iorder = lcblkb(4,iblk) 446 nenl = lcblkb(5,iblk) ! no. of vertices per element 447 nenbl = lcblkb(6,iblk) ! no. of vertices per bdry. face 448 mattyp = lcblkb(7,iblk) 449 ndofl = lcblkb(8,iblk) 450 nshl = lcblkb(9,iblk) 451 nshlb = lcblkb(10,iblk) 452 npro = lcblkb(1,iblk+1) - iel 453 if(lcsyst.eq.3) lcsyst=nenbl 454 ngaussb = nintb(lcsyst) 455c 456c.... compute and assemble the residuals corresponding to the 457c boundary integral 458c 459 460 allocate (tmpshpb(nshl,MAXQPT)) 461 allocate (tmpshglb(nsd,nshl,MAXQPT)) 462 463 tmpshpb(1:nshl,:) = shpb(lcsyst,1:nshl,:) 464 tmpshglb(:,1:nshl,:) = shglb(lcsyst,:,1:nshl,:) 465c 466 call AsBMFGSclr (y, x, 467 & tmpshpb, 468 & tmpshglb, 469 & mienb(iblk)%p, mmatb(iblk)%p, 470 & miBCB(iblk)%p, mBCB(iblk)%p, 471 & rest, rmest) 472c 473 deallocate ( tmpshpb ) 474 deallocate ( tmpshglb ) 475 476c.... end of boundary element loop 477c 478 enddo 479 480 481 ttim(80) = ttim(80) + tmr() 482c 483c.... --------------------> communications <------------------------- 484c 485 486 if (numpe > 1) then 487 call commu (rest , ilwork, 1 , 'in ') 488 489 call MPI_BARRIER (MPI_COMM_WORLD,ierr) 490 491 endif 492 493c 494c.... ----------------------> post processing <---------------------- 495c 496c.... satisfy the BCs on the residual 497c 498 call bc3ResSclr (y, iBC, BC, rest, iper, ilwork) 499 if (numpe > 1) then 500 call commu (rest , ilwork, 1 , 'out') 501 call MPI_BARRIER (MPI_COMM_WORLD,ierr) 502 endif 503c 504c.... return 505c 506 call timer ('Back ') 507 return 508 end 509