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 timedataC 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 timedataC 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 nshl = lcblk(10,iblk) 360 mattyp = lcblk(7,iblk) 361 ndofl = lcblk(8,iblk) 362 nsymdl = lcblk(9,iblk) 363 npro = lcblk(1,iblk+1) - iel 364 ngauss = nint(lcsyst) 365c 366c.... compute and assemble diffusive flux vector residual, qres, 367c and lumped mass matrix, rmass 368 369 allocate (tmpshp(nshl,MAXQPT)) 370 allocate (tmpshgl(nsd,nshl,MAXQPT)) 371 372 tmpshp(1:nshl,:) = shp(lcsyst,1:nshl,:) 373 tmpshgl(:,1:nshl,:) = shgl(lcsyst,:,1:nshl,:) 374 375 call AsIq (y, x, 376 & tmpshp, 377 & tmpshgl, 378 & mien(iblk)%p, mxmudmi(iblk)%p, 379 & qres, 380 & rmass) 381 382 deallocate ( tmpshp ) 383 deallocate ( tmpshgl ) 384 enddo 385 386c 387c.... take care of periodic boundary conditions 388c 389 390 call qpbc( rmass, qres, iBC, iper, ilwork ) 391c 392 endif ! computation of global diffusive flux 393c 394c.... loop over element blocks to compute element residuals 395c 396c 397c.... initialize the arrays 398c 399 res = zero 400 rmes = zero ! to avoid trap_uninitialized 401 if (lhs. eq. 1) lhsK = zero 402 if (iprec .ne. 0) BDiag = zero 403 flxID = zero 404c 405c.... loop over the element-blocks 406c 407 do iblk = 1, nelblk 408c 409c.... set up the parameters 410c 411 iblkts = iblk ! used in timeseries 412 nenl = lcblk(5,iblk) ! no. of vertices per element 413 iel = lcblk(1,iblk) 414 lelCat = lcblk(2,iblk) 415 lcsyst = lcblk(3,iblk) 416 iorder = lcblk(4,iblk) 417 nenl = lcblk(5,iblk) ! no. of vertices per element 418 nshl = lcblk(10,iblk) 419 mattyp = lcblk(7,iblk) 420 ndofl = lcblk(8,iblk) 421 nsymdl = lcblk(9,iblk) 422 npro = lcblk(1,iblk+1) - iel 423 inum = iel + npro - 1 424 ngauss = nint(lcsyst) 425c 426c.... compute and assemble the residual and tangent matrix 427c 428 429 if(lhs.eq.1) then 430 allocate (EGmass(npro,nedof,nedof)) 431 EGmass = zero 432 else 433 allocate (EGmass(1,1,1)) 434 endif 435 436 allocate (tmpshp(nshl,MAXQPT)) 437 allocate (tmpshgl(nsd,nshl,MAXQPT)) 438 tmpshp(1:nshl,:) = shp(lcsyst,1:nshl,:) 439 tmpshgl(:,1:nshl,:) = shgl(lcsyst,:,1:nshl,:) 440 441 call AsIGMR (y, ac, 442 & x, mxmudmi(iblk)%p, 443 & tmpshp, 444 & tmpshgl, mien(iblk)%p, 445 & mmat(iblk)%p, res, 446 & rmes, BDiag, 447 & qres, EGmass, 448 & rerr ) 449 if(lhs.eq.1) then 450c 451c.... satisfy the BC's on the implicit LHS 452c 453 call bc3LHS (iBC, BC, mien(iblk)%p, 454 & EGmass ) 455 456c 457c.... Fill-up the global sparse LHS mass matrix 458c 459 call fillsparseC( mien(iblk)%p, EGmass, 460 1 lhsK, row, col) 461 endif 462c 463 deallocate ( EGmass ) 464 deallocate ( tmpshp ) 465 deallocate ( tmpshgl ) 466c 467c.... end of interior element loop 468c 469 enddo 470!ifdef DEBUG !Nicholas Mati 471! call write_debug(myrank, 'res-afterAsIGMR'//char(0), 472! & 'res'//char(0), res, 473! & 'd'//char(0), nshg, nflow, lstep) 474! call write_debug(myrank, 'y-afterAsIGMR'//char(0), 475! & 'y'//char(0), y, 476! & 'd'//char(0), nshg, ndof, lstep) 477!endif //DEBUG 478c 479c.... --------------------> boundary elements <-------------------- 480c 481c.... loop over the boundary elements 482c 483 do iblk = 1, nelblb 484c 485c.... set up the parameters 486c 487 iel = lcblkb(1,iblk) 488 lelCat = lcblkb(2,iblk) 489 lcsyst = lcblkb(3,iblk) 490 iorder = lcblkb(4,iblk) 491 nenl = lcblkb(5,iblk) ! no. of vertices per element 492 nenbl = lcblkb(6,iblk) ! no. of vertices per bdry. face 493 mattyp = lcblkb(7,iblk) 494 ndofl = lcblkb(8,iblk) 495 nshl = lcblkb(9,iblk) 496 nshlb = lcblkb(10,iblk) 497 npro = lcblkb(1,iblk+1) - iel 498 if(lcsyst.eq.3) lcsyst=nenbl 499 ngaussb = nintb(lcsyst) 500 501c 502c.... compute and assemble the residuals corresponding to the 503c boundary integral 504c 505 506 allocate (tmpshpb(nshl,MAXQPT)) 507 allocate (tmpshglb(nsd,nshl,MAXQPT)) 508 if(lhs.eq.1 .and. iLHScond >= 1) then 509 allocate (EGmass(npro,nshl,nshl)) 510 EGmass = zero 511 else 512 allocate (EGmass(1,1,1)) 513 endif 514 515 tmpshpb(1:nshl,:) = shpb(lcsyst,1:nshl,:) 516 tmpshglb(:,1:nshl,:) = shglb(lcsyst,:,1:nshl,:) 517 518 call AsBMFG (y, x, 519 & tmpshpb, tmpshglb, 520 & mienb(iblk)%p, mmatb(iblk)%p, 521 & miBCB(iblk)%p, mBCB(iblk)%p, 522 & res, rmes, 523 & EGmass) 524 if(lhs == 1 .and. iLHScond > 0) then 525 call fillSparseC_BC(mienb(iblk)%p, EGmass, 526 & lhsk, row, col) 527 endif 528 529 deallocate (EGmass) 530 deallocate (tmpshpb) 531 deallocate (tmpshglb) 532 enddo !end of boundary element loop 533 534!ifdef DEBUG !Nicholas Mati 535! call write_debug(myrank, 'res-afterAsBMFG'//char(0), 536! & 'res'//char(0), res, 537! & 'd'//char(0), nshg, nflow, lstep) 538! call MPI_ABORT(MPI_COMM_WORLD) 539!endif //DEBUG 540 541 542c 543 ttim(80) = ttim(80) + secs(0.0) 544c 545c before the commu we need to rotate the residual vector for axisymmetric 546c boundary conditions (so that off processor periodicity is a dof add instead 547c of a dof combination). Take care of all nodes now so periodicity, like 548c commu is a simple dof add. 549c 550 if(iabc==1) then !are there any axisym bc's 551 call rotabc(res(1,2), iBC, 'in ') 552c Bdiagvec(:,1)=BDiag(:,1,1) 553c Bdiagvec(:,2)=BDiag(:,2,2) 554c Bdiagvec(:,3)=BDiag(:,3,3) 555c Bdiagvec(:,4)=BDiag(:,4,4) 556c Bdiagvec(:,5)=BDiag(:,5,5) 557c call rotabc(Bdiagvec(1,2), iBC, 'in ') 558c BDiag(:,:,:)=zero 559c BDiag(:,1,1)=Bdiagvec(:,1) 560c BDiag(:,2,2)=Bdiagvec(:,2) 561c BDiag(:,3,3)=Bdiagvec(:,3) 562c BDiag(:,4,4)=Bdiagvec(:,4) 563c BDiag(:,5,5)=Bdiagvec(:,5) 564 endif 565 566c.... --------------------> communications <------------------------- 567c 568 569 if (numpe > 1) then 570 call commu (res , ilwork, nflow , 'in ') 571 572 call MPI_BARRIER (MPI_COMM_WORLD,ierr) 573 574 if(iprec .ne. 0) call commu (BDiag, ilwork, nflow*nflow, 'in ') 575 endif 576 577c 578c.... ----------------------> post processing <---------------------- 579c 580c.... satisfy the BCs on the residual 581c 582 call bc3Res (y, iBC, BC, res, iper, ilwork) 583c 584c.... satisfy the BCs on the block-diagonal preconditioner 585c 586c This code fragment would extract the "on processor diagonal 587c blocks". lhsK alread has the BC's applied to it (using BC3lhs), 588c though it was on an ebe basis. For now, we forgo this and still 589c form BDiag before BC3lhs, leaving the need to still apply BC's 590c below. Keep in mind that if we used the code fragment below we 591c would still need to make BDiag periodic since BC3lhs did not do 592c that part. 593c 594 if (iprec .ne. 0) then 595c$$$ do iaa=1,nshg 596c$$$ k = sparseloc( row(col(iaa)), col(iaa+1)-colm(iaa), iaa ) 597c$$$ & + colm(iaa)-1 598c$$$ do idof=1,nflow 599c$$$ do jdof=1,nflow 600c$$$ idx=idof+(jdof-1)*nflow 601c$$$ BDiag(iaa,idof,jdof)=lhsK(idx,k) 602c$$$ enddo 603c$$$ enddo 604c$$$ enddo 605 call bc3BDg (y, iBC, BC, BDiag, iper, ilwork) 606 endif 607c 608c.... return 609c 610 call timer ('Back ') 611 return 612 end 613c 614c 615 616c 617 subroutine ElmGMRSclr(y, ac, 618 & x, elDw, 619 & shp, shgl, iBC, 620 & BC, shpb, shglb, 621 & rest, rmest, Diag, 622 & iper, ilwork, EGmasst) 623c 624c---------------------------------------------------------------------- 625c 626c This routine computes the LHS mass matrix, the RHS residual 627c vector, and the preconditioning matrix, for use with the GMRES 628c solver. 629c 630c Zdenek Johan, Winter 1991. (Fortran 90) 631c Chris Whiting, Winter 1998. (Matrix EBE-GMRES) 632c---------------------------------------------------------------------- 633c 634 use pointer_data 635c 636 include "common.h" 637 include "mpif.h" 638c 639 dimension y(nshg,ndof), ac(nshg,ndof), 640 & x(numnp,nsd), 641 & iBC(nshg), 642 & BC(nshg,ndofBC), 643 & rest(nshg), Diag(nshg), 644 & rmest(nshg), BDiag(nshg,nflow,nflow), 645 & iper(nshg), EGmasst(numel,nshape,nshape) 646c 647 dimension shp(MAXTOP,maxsh,MAXQPT), 648 & shgl(MAXTOP,nsd,maxsh,MAXQPT), 649 & shpb(MAXTOP,maxsh,MAXQPT), 650 & shglb(MAXTOP,nsd,maxsh,MAXQPT) 651c 652 dimension qrest(nshg), rmasst(nshg) 653c 654 dimension ilwork(nlwork) 655 dimension elDw(numel) 656c 657 real*8, allocatable :: tmpshp(:,:), tmpshgl(:,:,:) 658 real*8, allocatable :: tmpshpb(:,:), tmpshglb(:,:,:) 659 real*8, allocatable :: elDwl(:) 660c 661 ttim(80) = ttim(80) - tmr() 662c 663c.... set up the timer 664c 665 666 call timer ('Elm_Form') 667c 668c.... --------------------> interior elements <-------------------- 669c 670c.... set up parameters 671c 672 intrul = intg (1,itseq) 673 intind = intpt (intrul) 674c 675 ires = 1 676 677c if (idiff>=1) then ! global reconstruction of qdiff 678c 679c loop over element blocks for the global reconstruction 680c of the diffusive flux vector, q, and lumped mass matrix, rmass 681c 682c qrest = zero 683c rmasst = zero 684c 685c do iblk = 1, nelblk 686c 687c.... set up the parameters 688c 689c iel = lcblk(1,iblk) 690c lelCat = lcblk(2,iblk) 691c lcsyst = lcblk(3,iblk) 692c iorder = lcblk(4,iblk) 693c nenl = lcblk(5,iblk) ! no. of vertices per element 694c mattyp = lcblk(7,iblk) 695c ndofl = lcblk(8,iblk) 696c nsymdl = lcblk(9,iblk) 697c npro = lcblk(1,iblk+1) - iel 698c 699c nintg = numQpt (nsd,intrul,lcsyst) 700c 701c 702c.... compute and assemble diffusive flux vector residual, qres, 703c and lumped mass matrix, rmass 704c 705c call AsIq (y, x, 706c & shp(1,intind,lelCat), 707c & shgl(1,intind,lelCat), 708c & mien(iblk)%p, mxmudmi(iblk)%p, 709c & qres, rmass ) 710c 711c enddo 712c 713c 714c.... compute qi for node A, i.e., qres <-- qres/rmass 715c 716c if (numpe > 1) then 717c call commu (qres , ilwork, (ndof-1)*nsd , 'in ') 718c call commu (rmass , ilwork, 1 , 'in ') 719c endif 720c 721c.... take care of periodic boundary conditions 722c 723c call qpbc( rmass, qres, iBC, iper ) 724c 725c rmass = one/rmass 726c 727c do i=1, (nflow-1)*nsd 728c qres(:,i) = rmass*qres(:,i) 729c enddo 730c 731c if(numpe > 1) then 732c call commu (qres, ilwork, (nflow-1)*nsd, 'out') 733c endif 734c 735c endif ! computation of global diffusive flux 736c 737c.... loop over element blocks to compute element residuals 738c 739c 740c.... initialize the arrays 741c 742 rest = zero 743 rmest = zero ! to avoid trap_uninitialized 744 if (lhs .eq. 1) EGmasst = zero 745 if (iprec. ne. 0) Diag = zero 746c 747c.... loop over the element-blocks 748c 749 do iblk = 1, nelblk 750c 751c 752 nenl = lcblk(5,iblk) ! no. of vertices per element 753 iel = lcblk(1,iblk) 754 lelCat = lcblk(2,iblk) 755 lcsyst = lcblk(3,iblk) 756 iorder = lcblk(4,iblk) 757 nenl = lcblk(5,iblk) ! no. of vertices per element 758 nshl = lcblk(10,iblk) 759 mattyp = lcblk(7,iblk) 760 ndofl = lcblk(8,iblk) 761 nsymdl = lcblk(9,iblk) 762 npro = lcblk(1,iblk+1) - iel 763 inum = iel + npro - 1 764 ngauss = nint(lcsyst) 765c 766c.... compute and assemble the residual and tangent matrix 767c 768 allocate (tmpshp(nshl,MAXQPT)) 769 allocate (tmpshgl(nsd,nshl,MAXQPT)) 770 771 tmpshp(1:nshl,:) = shp(lcsyst,1:nshl,:) 772 tmpshgl(:,1:nshl,:) = shgl(lcsyst,:,1:nshl,:) 773 774 allocate (elDwl(npro)) 775c 776 call AsIGMRSclr(y, 777 & ac, 778 & x, elDwl, 779 & tmpshp, tmpshgl, 780 & mien(iblk)%p, 781 & mmat(iblk)%p, rest, 782 & rmest, 783 & qrest, EGmasst(iel:inum,:,:), 784 & Diag ) 785c 786 787c.... satisfy the BC's on the implicit LHS 788c 789 call bc3LHSSclr (iBC, mien(iblk)%p, EGmasst(iel:inum,:,:) ) 790c 791 elDw(iel:inum)=elDwl(1:npro) 792 deallocate ( elDwl ) 793 deallocate ( tmpshp ) 794 deallocate ( tmpshgl ) 795c.... end of interior element loop 796c 797 enddo 798c 799c.... --------------------> boundary elements <-------------------- 800c 801c.... set up parameters 802c 803 intrul = intg (2,itseq) 804 intind = intptb (intrul) 805c 806c.... loop over the boundary elements 807c 808 do iblk = 1, nelblb 809c 810c.... set up the parameters 811c 812 iel = lcblkb(1,iblk) 813 lelCat = lcblkb(2,iblk) 814 lcsyst = lcblkb(3,iblk) 815 iorder = lcblkb(4,iblk) 816 nenl = lcblkb(5,iblk) ! no. of vertices per element 817 nenbl = lcblkb(6,iblk) ! no. of vertices per bdry. face 818 mattyp = lcblkb(7,iblk) 819 ndofl = lcblkb(8,iblk) 820 nshl = lcblkb(9,iblk) 821 nshlb = lcblkb(10,iblk) 822 npro = lcblkb(1,iblk+1) - iel 823 if(lcsyst.eq.3) lcsyst=nenbl 824 ngaussb = nintb(lcsyst) 825c 826c.... compute and assemble the residuals corresponding to the 827c boundary integral 828c 829 830 allocate (tmpshpb(nshl,MAXQPT)) 831 allocate (tmpshglb(nsd,nshl,MAXQPT)) 832 833 tmpshpb(1:nshl,:) = shpb(lcsyst,1:nshl,:) 834 tmpshglb(:,1:nshl,:) = shglb(lcsyst,:,1:nshl,:) 835c 836 call AsBMFGSclr (y, x, 837 & tmpshpb, 838 & tmpshglb, 839 & mienb(iblk)%p, mmatb(iblk)%p, 840 & miBCB(iblk)%p, mBCB(iblk)%p, 841 & rest, rmest) 842c 843 deallocate ( tmpshpb ) 844 deallocate ( tmpshglb ) 845 846c.... end of boundary element loop 847c 848 enddo 849 850 851 ttim(80) = ttim(80) + tmr() 852c 853c.... --------------------> communications <------------------------- 854c 855 856 if (numpe > 1) then 857 call commu (rest , ilwork, 1 , 'in ') 858 859 call MPI_BARRIER (MPI_COMM_WORLD,ierr) 860 861 if(iprec .ne. 0) call commu (Diag, ilwork, 1, 'in ') 862 endif 863 864c 865c.... ----------------------> post processing <---------------------- 866c 867c.... satisfy the BCs on the residual 868c 869 call bc3ResSclr (y, iBC, BC, rest, iper, ilwork) 870c 871c.... satisfy the BCs on the preconditioner 872c 873 call bc3BDgSclr (iBC, Diag, iper, ilwork) 874c 875c.... return 876c 877 call timer ('Back ') 878 return 879 end 880 881