1 subroutine ElmGMRe (y, ac, x, 2 & shp, shgl, iBC, 3 & BC, shpb, shglb, 4 & res, rmes, BDiag, 5 & iper, ilwork, EGmass, rerr) 6c 7c---------------------------------------------------------------------- 8c 9c This routine computes the LHS mass matrix, the RHS residual 10c vector, and the preconditioning matrix, for use with the GMRES 11c solver. 12c 13c Zdenek Johan, Winter 1991. (Fortran 90) 14c Chris Whiting, Winter 1998. (Matrix EBE-GMRES) 15c---------------------------------------------------------------------- 16c 17 use pointer_data 18 use timedata 19c 20 include "common.h" 21 include "mpif.h" 22c 23 dimension y(nshg,ndof), ac(nshg,ndof), 24 & x(numnp,nsd), 25 & iBC(nshg), 26 & BC(nshg,ndofBC), 27 & res(nshg,nflow), 28 & rmes(nshg,nflow), BDiag(nshg,nflow,nflow), 29 & iper(nshg), EGmass(numel,nedof,nedof) 30c 31 dimension shp(MAXTOP,maxsh,MAXQPT), 32 & shgl(MAXTOP,nsd,maxsh,MAXQPT), 33 & shpb(MAXTOP,maxsh,MAXQPT), 34 & shglb(MAXTOP,nsd,maxsh,MAXQPT) 35c 36 dimension qres(nshg, idflx), rmass(nshg) 37c 38 dimension ilwork(nlwork) 39 40 real*8 Bdiagvec(nshg,nflow), rerr(nshg,10) 41 42 real*8, allocatable :: tmpshp(:,:), tmpshgl(:,:,:) 43 real*8, allocatable :: tmpshpb(:,:), tmpshglb(:,:,:) 44 45 ttim(80) = ttim(80) - secs(0.0) 46c 47c.... set up the timer 48c 49 50 call timer ('Elm_Form') 51c 52c.... --------------------> interior elements <-------------------- 53c 54c.... set up parameters 55c 56 ires = 1 57c 58 if (idiff==1 .or. idiff==3 .or. isurf==1) then ! global reconstruction 59 ! of qdiff 60c 61c loop over element blocks for the global reconstruction 62c of the diffusive flux vector, q, and lumped mass matrix, rmass 63c 64 qres = zero 65 rmass = zero 66 67 do iblk = 1, nelblk 68c 69c.... set up the parameters 70c 71 nenl = lcblk(5,iblk) ! no. of vertices per element 72 iel = lcblk(1,iblk) 73 lelCat = lcblk(2,iblk) 74 lcsyst = lcblk(3,iblk) 75 iorder = lcblk(4,iblk) 76 nenl = lcblk(5,iblk) ! no. of vertices per element 77 nshl = lcblk(10,iblk) 78 mattyp = lcblk(7,iblk) 79 ndofl = lcblk(8,iblk) 80 nsymdl = lcblk(9,iblk) 81 npro = lcblk(1,iblk+1) - iel 82 ngauss = nint(lcsyst) 83c 84c.... compute and assemble diffusive flux vector residual, qres, 85c and lumped mass matrix, rmass 86 87 allocate (tmpshp(nshl,MAXQPT)) 88 allocate (tmpshgl(nsd,nshl,MAXQPT)) 89 90 tmpshp(1:nshl,:) = shp(lcsyst,1:nshl,:) 91 tmpshgl(:,1:nshl,:) = shgl(lcsyst,:,1:nshl,:) 92 93 call AsIq (y, x, 94 & tmpshp, 95 & tmpshgl, 96 & mien(iblk)%p, mxmudmi(iblk)%p, 97 & qres, 98 & rmass) 99 100 deallocate ( tmpshp ) 101 deallocate ( tmpshgl ) 102 enddo 103 104c 105c.... take care of periodic boundary conditions 106c 107 108 call qpbc( rmass, qres, iBC, iper, ilwork ) 109c 110 endif ! computation of global diffusive flux 111c 112c.... loop over element blocks to compute element residuals 113c 114c 115c.... initialize the arrays 116c 117 res = zero 118 rmes = zero ! to avoid trap_uninitialized 119 if (lhs. eq. 1) EGmass = zero 120 if (iprec .ne. 0) BDiag = zero 121 flxID = zero 122 123c 124c.... loop over the element-blocks 125c 126 do iblk = 1, nelblk 127c 128c.... set up the parameters 129c 130 iblkts = iblk ! used in timeseries 131 nenl = lcblk(5,iblk) ! no. of vertices per element 132 iel = lcblk(1,iblk) 133 lelCat = lcblk(2,iblk) 134 lcsyst = lcblk(3,iblk) 135 iorder = lcblk(4,iblk) 136 nenl = lcblk(5,iblk) ! no. of vertices per element 137 nshl = lcblk(10,iblk) 138 mattyp = lcblk(7,iblk) 139 ndofl = lcblk(8,iblk) 140 nsymdl = lcblk(9,iblk) 141 npro = lcblk(1,iblk+1) - iel 142 inum = iel + npro - 1 143 ngauss = nint(lcsyst) 144c 145c.... compute and assemble the residual and tangent matrix 146c 147 148 allocate (tmpshp(nshl,MAXQPT)) 149 allocate (tmpshgl(nsd,nshl,MAXQPT)) 150 151 tmpshp(1:nshl,:) = shp(lcsyst,1:nshl,:) 152 tmpshgl(:,1:nshl,:) = shgl(lcsyst,:,1:nshl,:) 153 154 call AsIGMR (y, ac, 155 & x, mxmudmi(iblk)%p, 156 & tmpshp, 157 & tmpshgl, mien(iblk)%p, 158 & mmat(iblk)%p, res, 159 & rmes, BDiag, 160 & qres, EGmass(iel:inum,:,:), 161 & rerr) 162c 163c.... satisfy the BC's on the implicit LHS 164c 165 call bc3LHS (iBC, BC, mien(iblk)%p, 166 & EGmass(iel:inum,:,:) ) 167 168 deallocate ( tmpshp ) 169 deallocate ( tmpshgl ) 170c 171c.... end of interior element loop 172c 173 enddo 174c 175c.... --------------------> boundary elements <-------------------- 176c 177c.... loop over the boundary elements 178c 179 do iblk = 1, nelblb 180c 181c.... set up the parameters 182c 183 iel = lcblkb(1,iblk) 184 lelCat = lcblkb(2,iblk) 185 lcsyst = lcblkb(3,iblk) 186 iorder = lcblkb(4,iblk) 187 nenl = lcblkb(5,iblk) ! no. of vertices per element 188 nenbl = lcblkb(6,iblk) ! no. of vertices per bdry. face 189 mattyp = lcblkb(7,iblk) 190 ndofl = lcblkb(8,iblk) 191 nshl = lcblkb(9,iblk) 192 nshlb = lcblkb(10,iblk) 193 npro = lcblkb(1,iblk+1) - iel 194 if(lcsyst.eq.3) lcsyst=nenbl 195 ngaussb = nintb(lcsyst) 196 197c 198c.... compute and assemble the residuals corresponding to the 199c boundary integral 200c 201 202 allocate (tmpshpb(nshl,MAXQPT)) 203 allocate (tmpshglb(nsd,nshl,MAXQPT)) 204 205 tmpshpb(1:nshl,:) = shpb(lcsyst,1:nshl,:) 206 tmpshglb(:,1:nshl,:) = shglb(lcsyst,:,1:nshl,:) 207 208 call AsBMFG (y, x, 209 & tmpshpb, tmpshglb, 210 & mienb(iblk)%p, mmatb(iblk)%p, 211 & miBCB(iblk)%p, mBCB(iblk)%p, 212 & res, rmes) 213 214 deallocate (tmpshpb) 215 deallocate (tmpshglb) 216c 217c.... end of boundary element loop 218c 219 enddo 220c 221 ttim(80) = ttim(80) + secs(0.0) 222c 223c before the commu we need to rotate the residual vector for axisymmetric 224c boundary conditions (so that off processor periodicity is a dof add instead 225c of a dof combination). Take care of all nodes now so periodicity, like 226c commu is a simple dof add. 227c 228c if(iabc==1) !are there any axisym bc's 229c & call rotabc(res(1,2), iBC, BC, nflow, 'in ') 230 if(iabc==1) then !are there any axisym bc's 231 call rotabc(res(1,2), iBC, 'in ') 232c Bdiagvec(:,1)=BDiag(:,1,1) 233c Bdiagvec(:,2)=BDiag(:,2,2) 234c Bdiagvec(:,3)=BDiag(:,3,3) 235c Bdiagvec(:,4)=BDiag(:,4,4) 236c Bdiagvec(:,5)=BDiag(:,5,5) 237c call rotabc(Bdiagvec(1,2), iBC, BC, 2, 'in ') 238c BDiag(:,:,:)=zero 239c BDiag(:,1,1)=Bdiagvec(:,1) 240c BDiag(:,2,2)=Bdiagvec(:,2) 241c BDiag(:,3,3)=Bdiagvec(:,3) 242c BDiag(:,4,4)=Bdiagvec(:,4) 243c BDiag(:,5,5)=Bdiagvec(:,5) 244 endif 245 246c.... --------------------> communications <------------------------- 247c 248 249 if (numpe > 1) then 250 call commu (res , ilwork, nflow , 'in ') 251 252 call MPI_BARRIER (MPI_COMM_WORLD,ierr) 253 254 if(iprec .ne. 0) call commu (BDiag, ilwork, nflow*nflow, 'in ') 255 endif 256 257c 258c.... ----------------------> post processing <---------------------- 259c 260c.... satisfy the BCs on the residual 261c 262 call bc3Res (y, iBC, BC, res, iper, ilwork) 263c 264c.... satisfy the BCs on the block-diagonal preconditioner 265c 266 if (iprec .ne. 0) then 267 call bc3BDg (y, iBC, BC, BDiag, iper, ilwork) 268 endif 269c 270c.... return 271c 272 call timer ('Back ') 273 return 274 end 275c 276cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 277ccccccccccccc SPARSE 278c_______________________________________________________________ 279 280 subroutine ElmGMRs (y, ac, x, 281 & shp, shgl, iBC, 282 & BC, shpb, shglb, 283 & res, rmes, BDiag, 284 & iper, ilwork, lhsK, 285 & col, row, rerr) 286c 287c---------------------------------------------------------------------- 288c 289c This routine computes the LHS mass matrix, the RHS residual 290c vector, and the preconditioning matrix, for use with the GMRES 291c solver. 292c 293c Zdenek Johan, Winter 1991. (Fortran 90) 294c Chris Whiting, Winter 1998. (Matrix EBE-GMRES) 295c---------------------------------------------------------------------- 296c 297 use pointer_data 298 use timedata 299c 300 include "common.h" 301 include "mpif.h" 302c 303 integer col(nshg+1), row(nnz*nshg) 304 real*8 lhsK(nflow*nflow,nnz_tot) 305 306 dimension y(nshg,ndof), ac(nshg,ndof), 307 & x(numnp,nsd), 308 & iBC(nshg), 309 & BC(nshg,ndofBC), 310 & res(nshg,nflow), 311 & rmes(nshg,nflow), BDiag(nshg,nflow,nflow), 312 & iper(nshg) 313c 314 dimension shp(MAXTOP,maxsh,MAXQPT), 315 & shgl(MAXTOP,nsd,maxsh,MAXQPT), 316 & shpb(MAXTOP,maxsh,MAXQPT), 317 & shglb(MAXTOP,nsd,maxsh,MAXQPT) 318c 319 dimension qres(nshg, idflx), rmass(nshg) 320c 321 dimension ilwork(nlwork) 322 323 real*8 Bdiagvec(nshg,nflow), rerr(nshg,10) 324 325 real*8, allocatable :: tmpshp(:,:), tmpshgl(:,:,:) 326 real*8, allocatable :: tmpshpb(:,:), tmpshglb(:,:,:) 327 real*8, allocatable :: EGmass(:,:,:) 328 ttim(80) = ttim(80) - secs(0.0) 329c 330c.... set up the timer 331c 332 333 call timer ('Elm_Form') 334c 335c.... --------------------> interior elements <-------------------- 336c 337c.... set up parameters 338c 339 ires = 1 340c 341 if (idiff==1 .or. idiff==3 .or. isurf==1) then ! global reconstruction 342 ! of qdiff 343c 344c loop over element blocks for the global reconstruction 345c of the diffusive flux vector, q, and lumped mass matrix, rmass 346c 347 qres = zero 348 rmass = zero 349 350 do iblk = 1, nelblk 351c 352c.... set up the parameters 353c 354 nenl = lcblk(5,iblk) ! no. of vertices per element 355 iel = lcblk(1,iblk) 356 lelCat = lcblk(2,iblk) 357 lcsyst = lcblk(3,iblk) 358 iorder = lcblk(4,iblk) 359 nenl = lcblk(5,iblk) ! no. of vertices per element 360 nshl = lcblk(10,iblk) 361 mattyp = lcblk(7,iblk) 362 ndofl = lcblk(8,iblk) 363 nsymdl = lcblk(9,iblk) 364 npro = lcblk(1,iblk+1) - iel 365 ngauss = nint(lcsyst) 366c 367c.... compute and assemble diffusive flux vector residual, qres, 368c and lumped mass matrix, rmass 369 370 allocate (tmpshp(nshl,MAXQPT)) 371 allocate (tmpshgl(nsd,nshl,MAXQPT)) 372 373 tmpshp(1:nshl,:) = shp(lcsyst,1:nshl,:) 374 tmpshgl(:,1:nshl,:) = shgl(lcsyst,:,1:nshl,:) 375 376 call AsIq (y, x, 377 & tmpshp, 378 & tmpshgl, 379 & mien(iblk)%p, mxmudmi(iblk)%p, 380 & qres, 381 & rmass) 382 383 deallocate ( tmpshp ) 384 deallocate ( tmpshgl ) 385 enddo 386 387c 388c.... take care of periodic boundary conditions 389c 390 391 call qpbc( rmass, qres, iBC, iper, ilwork ) 392c 393 endif ! computation of global diffusive flux 394c 395c.... loop over element blocks to compute element residuals 396c 397c 398c.... initialize the arrays 399c 400 res = zero 401 rmes = zero ! to avoid trap_uninitialized 402 if (lhs. eq. 1) lhsK = zero 403 if (iprec .ne. 0) BDiag = zero 404 flxID = zero 405c 406c.... loop over the element-blocks 407c 408 do iblk = 1, nelblk 409c 410c.... set up the parameters 411c 412 iblkts = iblk ! used in timeseries 413 nenl = lcblk(5,iblk) ! no. of vertices per element 414 iel = lcblk(1,iblk) 415 lelCat = lcblk(2,iblk) 416 lcsyst = lcblk(3,iblk) 417 iorder = lcblk(4,iblk) 418 nenl = lcblk(5,iblk) ! no. of vertices per element 419 nshl = lcblk(10,iblk) 420 mattyp = lcblk(7,iblk) 421 ndofl = lcblk(8,iblk) 422 nsymdl = lcblk(9,iblk) 423 npro = lcblk(1,iblk+1) - iel 424 inum = iel + npro - 1 425 ngauss = nint(lcsyst) 426c 427c.... compute and assemble the residual and tangent matrix 428c 429 430 if(lhs.eq.1) then 431 allocate (EGmass(npro,nedof,nedof)) 432 EGmass = zero 433 else 434 allocate (EGmass(1,1,1)) 435 endif 436 437 allocate (tmpshp(nshl,MAXQPT)) 438 allocate (tmpshgl(nsd,nshl,MAXQPT)) 439 tmpshp(1:nshl,:) = shp(lcsyst,1:nshl,:) 440 tmpshgl(:,1:nshl,:) = shgl(lcsyst,:,1:nshl,:) 441 442 call AsIGMR (y, ac, 443 & x, mxmudmi(iblk)%p, 444 & tmpshp, 445 & tmpshgl, mien(iblk)%p, 446 & mmat(iblk)%p, res, 447 & rmes, BDiag, 448 & qres, EGmass, 449 & rerr ) 450 if(lhs.eq.1) then 451c 452c.... satisfy the BC's on the implicit LHS 453c 454 call bc3LHS (iBC, BC, mien(iblk)%p, 455 & EGmass ) 456 457c 458c.... Fill-up the global sparse LHS mass matrix 459c 460 call fillsparseC( mien(iblk)%p, EGmass, 461 1 lhsK, row, col) 462 endif 463c 464 deallocate ( EGmass ) 465 deallocate ( tmpshp ) 466 deallocate ( tmpshgl ) 467c 468c.... end of interior element loop 469c 470 enddo 471c 472c.... --------------------> boundary elements <-------------------- 473c 474c.... loop over the boundary elements 475c 476 do iblk = 1, nelblb 477c 478c.... set up the parameters 479c 480 iel = lcblkb(1,iblk) 481 lelCat = lcblkb(2,iblk) 482 lcsyst = lcblkb(3,iblk) 483 iorder = lcblkb(4,iblk) 484 nenl = lcblkb(5,iblk) ! no. of vertices per element 485 nenbl = lcblkb(6,iblk) ! no. of vertices per bdry. face 486 mattyp = lcblkb(7,iblk) 487 ndofl = lcblkb(8,iblk) 488 nshl = lcblkb(9,iblk) 489 nshlb = lcblkb(10,iblk) 490 npro = lcblkb(1,iblk+1) - iel 491 if(lcsyst.eq.3) lcsyst=nenbl 492 ngaussb = nintb(lcsyst) 493 494c 495c.... compute and assemble the residuals corresponding to the 496c boundary integral 497c 498 499 allocate (tmpshpb(nshl,MAXQPT)) 500 allocate (tmpshglb(nsd,nshl,MAXQPT)) 501 502 tmpshpb(1:nshl,:) = shpb(lcsyst,1:nshl,:) 503 tmpshglb(:,1:nshl,:) = shglb(lcsyst,:,1:nshl,:) 504 505 call AsBMFG (y, x, 506 & tmpshpb, tmpshglb, 507 & mienb(iblk)%p, mmatb(iblk)%p, 508 & miBCB(iblk)%p, mBCB(iblk)%p, 509 & res, rmes) 510 511 deallocate (tmpshpb) 512 deallocate (tmpshglb) 513c 514c.... end of boundary element loop 515c 516 enddo 517c 518 ttim(80) = ttim(80) + secs(0.0) 519c 520c before the commu we need to rotate the residual vector for axisymmetric 521c boundary conditions (so that off processor periodicity is a dof add instead 522c of a dof combination). Take care of all nodes now so periodicity, like 523c commu is a simple dof add. 524c 525 if(iabc==1) then !are there any axisym bc's 526 call rotabc(res(1,2), iBC, 'in ') 527c Bdiagvec(:,1)=BDiag(:,1,1) 528c Bdiagvec(:,2)=BDiag(:,2,2) 529c Bdiagvec(:,3)=BDiag(:,3,3) 530c Bdiagvec(:,4)=BDiag(:,4,4) 531c Bdiagvec(:,5)=BDiag(:,5,5) 532c call rotabc(Bdiagvec(1,2), iBC, 'in ') 533c BDiag(:,:,:)=zero 534c BDiag(:,1,1)=Bdiagvec(:,1) 535c BDiag(:,2,2)=Bdiagvec(:,2) 536c BDiag(:,3,3)=Bdiagvec(:,3) 537c BDiag(:,4,4)=Bdiagvec(:,4) 538c BDiag(:,5,5)=Bdiagvec(:,5) 539 endif 540 541c.... --------------------> communications <------------------------- 542c 543 544 if (numpe > 1) then 545 call commu (res , ilwork, nflow , 'in ') 546 547 call MPI_BARRIER (MPI_COMM_WORLD,ierr) 548 549 if(iprec .ne. 0) call commu (BDiag, ilwork, nflow*nflow, 'in ') 550 endif 551 552c 553c.... ----------------------> post processing <---------------------- 554c 555c.... satisfy the BCs on the residual 556c 557 call bc3Res (y, iBC, BC, res, iper, ilwork) 558c 559c.... satisfy the BCs on the block-diagonal preconditioner 560c 561c This code fragment would extract the "on processor diagonal 562c blocks". lhsK alread has the BC's applied to it (using BC3lhs), 563c though it was on an ebe basis. For now, we forgo this and still 564c form BDiag before BC3lhs, leaving the need to still apply BC's 565c below. Keep in mind that if we used the code fragment below we 566c would still need to make BDiag periodic since BC3lhs did not do 567c that part. 568c 569 if (iprec .ne. 0) then 570c$$$ do iaa=1,nshg 571c$$$ k = sparseloc( row(col(iaa)), col(iaa+1)-colm(iaa), iaa ) 572c$$$ & + colm(iaa)-1 573c$$$ do idof=1,nflow 574c$$$ do jdof=1,nflow 575c$$$ idx=idof+(jdof-1)*nflow 576c$$$ BDiag(iaa,idof,jdof)=lhsK(idx,k) 577c$$$ enddo 578c$$$ enddo 579c$$$ enddo 580 call bc3BDg (y, iBC, BC, BDiag, iper, ilwork) 581 endif 582c 583c.... return 584c 585 call timer ('Back ') 586 return 587 end 588c 589c 590 591c 592 subroutine ElmGMRSclr(y, ac, 593 & x, elDw, 594 & shp, shgl, iBC, 595 & BC, shpb, shglb, 596 & rest, rmest, Diag, 597 & iper, ilwork, EGmasst) 598c 599c---------------------------------------------------------------------- 600c 601c This routine computes the LHS mass matrix, the RHS residual 602c vector, and the preconditioning matrix, for use with the GMRES 603c solver. 604c 605c Zdenek Johan, Winter 1991. (Fortran 90) 606c Chris Whiting, Winter 1998. (Matrix EBE-GMRES) 607c---------------------------------------------------------------------- 608c 609 use pointer_data 610c 611 include "common.h" 612 include "mpif.h" 613c 614 dimension y(nshg,ndof), ac(nshg,ndof), 615 & x(numnp,nsd), 616 & iBC(nshg), 617 & BC(nshg,ndofBC), 618 & rest(nshg), Diag(nshg), 619 & rmest(nshg), BDiag(nshg,nflow,nflow), 620 & iper(nshg), EGmasst(numel,nshape,nshape) 621c 622 dimension shp(MAXTOP,maxsh,MAXQPT), 623 & shgl(MAXTOP,nsd,maxsh,MAXQPT), 624 & shpb(MAXTOP,maxsh,MAXQPT), 625 & shglb(MAXTOP,nsd,maxsh,MAXQPT) 626c 627 dimension qrest(nshg), rmasst(nshg) 628c 629 dimension ilwork(nlwork) 630 dimension elDw(numel) 631c 632 real*8, allocatable :: tmpshp(:,:), tmpshgl(:,:,:) 633 real*8, allocatable :: tmpshpb(:,:), tmpshglb(:,:,:) 634 real*8, allocatable :: elDwl(:) 635c 636 ttim(80) = ttim(80) - tmr() 637c 638c.... set up the timer 639c 640 641 call timer ('Elm_Form') 642c 643c.... --------------------> interior elements <-------------------- 644c 645c.... set up parameters 646c 647 intrul = intg (1,itseq) 648 intind = intpt (intrul) 649c 650 ires = 1 651 652c if (idiff>=1) then ! global reconstruction of qdiff 653c 654c loop over element blocks for the global reconstruction 655c of the diffusive flux vector, q, and lumped mass matrix, rmass 656c 657c qrest = zero 658c rmasst = zero 659c 660c do iblk = 1, nelblk 661c 662c.... set up the parameters 663c 664c iel = lcblk(1,iblk) 665c lelCat = lcblk(2,iblk) 666c lcsyst = lcblk(3,iblk) 667c iorder = lcblk(4,iblk) 668c nenl = lcblk(5,iblk) ! no. of vertices per element 669c mattyp = lcblk(7,iblk) 670c ndofl = lcblk(8,iblk) 671c nsymdl = lcblk(9,iblk) 672c npro = lcblk(1,iblk+1) - iel 673c 674c nintg = numQpt (nsd,intrul,lcsyst) 675c 676c 677c.... compute and assemble diffusive flux vector residual, qres, 678c and lumped mass matrix, rmass 679c 680c call AsIq (y, x, 681c & shp(1,intind,lelCat), 682c & shgl(1,intind,lelCat), 683c & mien(iblk)%p, mxmudmi(iblk)%p, 684c & qres, rmass ) 685c 686c enddo 687c 688c 689c.... compute qi for node A, i.e., qres <-- qres/rmass 690c 691c if (numpe > 1) then 692c call commu (qres , ilwork, (ndof-1)*nsd , 'in ') 693c call commu (rmass , ilwork, 1 , 'in ') 694c endif 695c 696c.... take care of periodic boundary conditions 697c 698c call qpbc( rmass, qres, iBC, iper ) 699c 700c rmass = one/rmass 701c 702c do i=1, (nflow-1)*nsd 703c qres(:,i) = rmass*qres(:,i) 704c enddo 705c 706c if(numpe > 1) then 707c call commu (qres, ilwork, (nflow-1)*nsd, 'out') 708c endif 709c 710c endif ! computation of global diffusive flux 711c 712c.... loop over element blocks to compute element residuals 713c 714c 715c.... initialize the arrays 716c 717 rest = zero 718 rmest = zero ! to avoid trap_uninitialized 719 if (lhs .eq. 1) EGmasst = zero 720 if (iprec. ne. 0) Diag = zero 721c 722c.... loop over the element-blocks 723c 724 do iblk = 1, nelblk 725c 726c 727 nenl = lcblk(5,iblk) ! no. of vertices per element 728 iel = lcblk(1,iblk) 729 lelCat = lcblk(2,iblk) 730 lcsyst = lcblk(3,iblk) 731 iorder = lcblk(4,iblk) 732 nenl = lcblk(5,iblk) ! no. of vertices per element 733 nshl = lcblk(10,iblk) 734 mattyp = lcblk(7,iblk) 735 ndofl = lcblk(8,iblk) 736 nsymdl = lcblk(9,iblk) 737 npro = lcblk(1,iblk+1) - iel 738 inum = iel + npro - 1 739 ngauss = nint(lcsyst) 740c 741c.... compute and assemble the residual and tangent matrix 742c 743 allocate (tmpshp(nshl,MAXQPT)) 744 allocate (tmpshgl(nsd,nshl,MAXQPT)) 745 746 tmpshp(1:nshl,:) = shp(lcsyst,1:nshl,:) 747 tmpshgl(:,1:nshl,:) = shgl(lcsyst,:,1:nshl,:) 748 749 allocate (elDwl(npro)) 750c 751 call AsIGMRSclr(y, 752 & ac, 753 & x, elDwl, 754 & tmpshp, tmpshgl, 755 & mien(iblk)%p, 756 & mmat(iblk)%p, rest, 757 & rmest, 758 & qrest, EGmasst(iel:inum,:,:), 759 & Diag ) 760c 761 762c.... satisfy the BC's on the implicit LHS 763c 764 call bc3LHSSclr (iBC, mien(iblk)%p, EGmasst(iel:inum,:,:) ) 765c 766 elDw(iel:iel+npro)=elDwl(1:npro) 767 deallocate ( elDwl ) 768 deallocate ( tmpshp ) 769 deallocate ( tmpshgl ) 770c.... end of interior element loop 771c 772 enddo 773c 774c.... --------------------> boundary elements <-------------------- 775c 776c.... set up parameters 777c 778 intrul = intg (2,itseq) 779 intind = intptb (intrul) 780c 781c.... loop over the boundary elements 782c 783 do iblk = 1, nelblb 784c 785c.... set up the parameters 786c 787 iel = lcblkb(1,iblk) 788 lelCat = lcblkb(2,iblk) 789 lcsyst = lcblkb(3,iblk) 790 iorder = lcblkb(4,iblk) 791 nenl = lcblkb(5,iblk) ! no. of vertices per element 792 nenbl = lcblkb(6,iblk) ! no. of vertices per bdry. face 793 mattyp = lcblkb(7,iblk) 794 ndofl = lcblkb(8,iblk) 795 nshl = lcblkb(9,iblk) 796 nshlb = lcblkb(10,iblk) 797 npro = lcblkb(1,iblk+1) - iel 798 if(lcsyst.eq.3) lcsyst=nenbl 799 ngaussb = nintb(lcsyst) 800c 801c.... compute and assemble the residuals corresponding to the 802c boundary integral 803c 804 805 allocate (tmpshpb(nshl,MAXQPT)) 806 allocate (tmpshglb(nsd,nshl,MAXQPT)) 807 808 tmpshpb(1:nshl,:) = shpb(lcsyst,1:nshl,:) 809 tmpshglb(:,1:nshl,:) = shglb(lcsyst,:,1:nshl,:) 810c 811 call AsBMFGSclr (y, x, 812 & tmpshpb, 813 & tmpshglb, 814 & mienb(iblk)%p, mmatb(iblk)%p, 815 & miBCB(iblk)%p, mBCB(iblk)%p, 816 & rest, rmest) 817c 818 deallocate ( tmpshpb ) 819 deallocate ( tmpshglb ) 820 821c.... end of boundary element loop 822c 823 enddo 824 825 826 ttim(80) = ttim(80) + tmr() 827c 828c.... --------------------> communications <------------------------- 829c 830 831 if (numpe > 1) then 832 call commu (rest , ilwork, 1 , 'in ') 833 834 call MPI_BARRIER (MPI_COMM_WORLD,ierr) 835 836 if(iprec .ne. 0) call commu (Diag, ilwork, 1, 'in ') 837 endif 838 839c 840c.... ----------------------> post processing <---------------------- 841c 842c.... satisfy the BCs on the residual 843c 844 call bc3ResSclr (y, iBC, BC, rest, iper, ilwork) 845c 846c.... satisfy the BCs on the preconditioner 847c 848 call bc3BDgSclr (iBC, Diag, iper, ilwork) 849c 850c.... return 851c 852 call timer ('Back ') 853 return 854 end 855 856