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 if (numpe < 1) then 273 call commu (res , ilwork, nflow , 'in ') 274 275 call MPI_BARRIER (MPI_COMM_WORLD,ierr) 276 endif 277 278c 279c.... ----------------------> post processing <---------------------- 280c 281c.... satisfy the BCs on the residual 282c 283 call bc3Res (y, iBC, BC, res, iper, ilwork) 284c 285c 286c.... return 287c 288! if (numpe > 1) then 289 if (numpe < 1) then 290 call commu (res , ilwork, nflow , 'out') 291 call MPI_BARRIER (MPI_COMM_WORLD,ierr) 292 endif 293 call timer ('Back ') 294 return 295 end 296c 297c 298 299c 300 301 subroutine ElmGMRPETScSclr(y, ac, 302 & x, elDw, 303 & shp, shgl, iBC, 304 & BC, shpb, shglb, 305 & rest, rmest, 306 & iper, ilwork, lhsPs) 307c 308c---------------------------------------------------------------------- 309c 310c This routine computes the LHS mass matrix, the RHS residual 311c vector, and the preconditioning matrix, for use with the GMRES 312c solver. 313c 314c Zdenek Johan, Winter 1991. (Fortran 90) 315c Chris Whiting, Winter 1998. (Matrix EBE-GMRES) 316c---------------------------------------------------------------------- 317c 318 use pointer_data 319c 320 include "common.h" 321 include "mpif.h" 322c 323 dimension y(nshg,ndof), ac(nshg,ndof), 324 & x(numnp,nsd), 325 & iBC(nshg), 326 & BC(nshg,ndofBC), 327 & rest(nshg), Diag(1), 328 & rmest(nshg), 329 & iper(nshg) 330c 331 dimension shp(MAXTOP,maxsh,MAXQPT), 332 & shgl(MAXTOP,nsd,maxsh,MAXQPT), 333 & shpb(MAXTOP,maxsh,MAXQPT), 334 & shglb(MAXTOP,nsd,maxsh,MAXQPT) 335c 336 dimension qrest(nshg), rmasst(nshg) 337 real*8 elDw(numel) 338c 339 dimension ilwork(nlwork) 340c 341 real*8, allocatable :: tmpshp(:,:), tmpshgl(:,:,:) 342 real*8, allocatable :: tmpshpb(:,:), tmpshglb(:,:,:) 343 real*8, allocatable :: EGMasst(:,:,:) 344 real*8, allocatable :: ElDwl(:) 345c 346 ttim(80) = ttim(80) - tmr() 347 iprec=0 348c 349c.... set up the timer 350c 351 352 call timer ('Elm_Form') 353c 354c.... --------------------> interior elements <-------------------- 355c 356c.... set up parameters 357c 358 intrul = intg (1,itseq) 359 intind = intpt (intrul) 360c 361 ires = 1 362 363c 364c.... initialize the arrays 365c 366 rest = zero 367 rmest = zero ! to avoid trap_uninitialized 368c 369c.... loop over the element-blocks 370c 371 do iblk = 1, nelblk 372c 373c 374 nenl = lcblk(5,iblk) ! no. of vertices per element 375 iel = lcblk(1,iblk) 376 lelCat = lcblk(2,iblk) 377 lcsyst = lcblk(3,iblk) 378 iorder = lcblk(4,iblk) 379 nshl = lcblk(10,iblk) 380 mattyp = lcblk(7,iblk) 381 ndofl = lcblk(8,iblk) 382 nsymdl = lcblk(9,iblk) 383 npro = lcblk(1,iblk+1) - iel 384 inum = iel + npro - 1 385 ngauss = nint(lcsyst) 386c 387c.... compute and assemble the residual and tangent matrix 388c 389 allocate (tmpshp(nshl,MAXQPT)) 390 allocate (tmpshgl(nsd,nshl,MAXQPT)) 391 if(lhs.eq.1) then 392 allocate (EGMasst(npro,nshl,nshl)) 393 EGMasst=zero 394 endif 395 396 tmpshp(1:nshl,:) = shp(lcsyst,1:nshl,:) 397 tmpshgl(:,1:nshl,:) = shgl(lcsyst,:,1:nshl,:) 398 399 allocate (elDwl(npro)) 400c 401 call AsIGMRSclr(y, 402 & ac, 403 & x, elDwl, 404 & tmpshp, tmpshgl, 405 & mien(iblk)%p, 406 & mmat(iblk)%p, rest, 407 & rmest, 408 & qrest, EGmasst, 409 & Diag ) 410c 411 elDw(iel:inum) = elDwl(1:npro) 412 deallocate ( elDwl ) 413 414 if(lhs.eq.1) then 415c.... satisfy the BCs on the implicit LHS 416c 417 call bc3LHSSclr (iBC, mien(iblk)%p, EGmasst ) 418c 419c 420c.... Fill-up the global sparse LHS mass matrix 421c 422 call fillsparsecpetscs( mieng(iblk)%p, EGmasst, lhsPs) 423 endif 424 if(lhs.eq.1) deallocate ( EGmasst ) 425 deallocate ( tmpshp ) 426 deallocate ( tmpshgl ) 427c.... end of interior element loop 428c 429 enddo 430c 431c.... --------------------> boundary elements <-------------------- 432c 433c.... set up parameters 434c 435 intrul = intg (2,itseq) 436 intind = intptb (intrul) 437c 438c.... loop over the boundary elements 439c 440 do iblk = 1, nelblb 441c 442c.... set up the parameters 443c 444 iel = lcblkb(1,iblk) 445 lelCat = lcblkb(2,iblk) 446 lcsyst = lcblkb(3,iblk) 447 iorder = lcblkb(4,iblk) 448 nenl = lcblkb(5,iblk) ! no. of vertices per element 449 nenbl = lcblkb(6,iblk) ! no. of vertices per bdry. face 450 mattyp = lcblkb(7,iblk) 451 ndofl = lcblkb(8,iblk) 452 nshl = lcblkb(9,iblk) 453 nshlb = lcblkb(10,iblk) 454 npro = lcblkb(1,iblk+1) - iel 455 if(lcsyst.eq.3) lcsyst=nenbl 456 ngaussb = nintb(lcsyst) 457c 458c.... compute and assemble the residuals corresponding to the 459c boundary integral 460c 461 462 allocate (tmpshpb(nshl,MAXQPT)) 463 allocate (tmpshglb(nsd,nshl,MAXQPT)) 464 465 tmpshpb(1:nshl,:) = shpb(lcsyst,1:nshl,:) 466 tmpshglb(:,1:nshl,:) = shglb(lcsyst,:,1:nshl,:) 467c 468 call AsBMFGSclr (y, x, 469 & tmpshpb, 470 & tmpshglb, 471 & mienb(iblk)%p, mmatb(iblk)%p, 472 & miBCB(iblk)%p, mBCB(iblk)%p, 473 & rest, rmest) 474c 475 deallocate ( tmpshpb ) 476 deallocate ( tmpshglb ) 477 478c.... end of boundary element loop 479c 480 enddo 481 482 483 ttim(80) = ttim(80) + tmr() 484c 485c.... --------------------> communications <------------------------- 486c 487 488 if (numpe < 1) then 489 call commu (rest , ilwork, 1 , 'in ') 490 491 call MPI_BARRIER (MPI_COMM_WORLD,ierr) 492 493 endif 494 495c 496c.... ----------------------> post processing <---------------------- 497c 498c.... satisfy the BCs on the residual 499c 500 call bc3ResSclr (y, iBC, BC, rest, iper, ilwork) 501 if (numpe < 1) then 502 call commu (rest , ilwork, 1 , 'out') 503 call MPI_BARRIER (MPI_COMM_WORLD,ierr) 504 endif 505c 506c.... return 507c 508 call timer ('Back ') 509 return 510 end 511