1 subroutine ElmMFG (y, ac, x, 2 & shp, shgl, iBC, 3 & BC, shpb, shglb, 4 & res, rmes, BDiag, 5 & iper, ilwork, rerr) 6c 7c---------------------------------------------------------------------- 8c 9c This routine calculates the preconditioning matrix and the 10c R.H.S. residual vector for the Matrix-free Iterative solver. 11c 12c 13c Zdenek Johan, Winter 1991. (Fortran 90) 14c---------------------------------------------------------------------- 15c 16 use pointer_data 17c 18 include "common.h" 19 include "mpif.h" 20c 21 dimension y(nshg,ndof), ac(nshg,ndof), 22 & x(numnp,nsd), 23 & iBC(nshg), 24 & BC(nshg,ndofBC), 25 & res(nshg,nflow), 26 & rmes(nshg,nflow), BDiag(nshg,nflow,nflow), 27 & iper(nshg) 28c 29 dimension shp(MAXTOP,maxsh,MAXQPT), 30 & shgl(MAXTOP,nsd,maxsh,MAXQPT), 31 & shpb(MAXTOP,maxsh,MAXQPT), 32 & shglb(MAXTOP,nsd,maxsh,MAXQPT) 33c 34 dimension qres(nshg,idflx), rmass(nshg), 35 & tmp(nshg) 36c 37 dimension ilwork(nlwork) 38 39 real*8 rerr(nshg,10) 40 real*8, allocatable :: tmpshp(:,:), tmpshgl(:,:,:) 41 real*8, allocatable :: tmpshpb(:,:), tmpshglb(:,:,:) 42 43 ttim(80) = ttim(80) - secs(0.0) 44c 45c.... set up the timer 46c 47 48 call timer ('Elm_Form') 49c 50c.... --------------------> interior elements <-------------------- 51c 52c.... set up parameters 53c 54 ires = 3 55 56 if (idiff==1 .or. idiff==3 .or. isurf==1) then ! global reconstruction 57 ! of qdiff 58c 59c loop over element blocks for the global reconstruction 60c of the diffusive flux vector, q, and lumped mass matrix, rmass 61c 62 qres = zero 63 rmass = zero 64 65 do iblk = 1, nelblk 66c 67c.... set up the parameters 68c 69c 70 nenl = lcblk(5,iblk) ! no. of vertices per element 71 iel = lcblk(1,iblk) 72 lelCat = lcblk(2,iblk) 73 lcsyst = lcblk(3,iblk) 74 iorder = lcblk(4,iblk) 75 nenl = lcblk(5,iblk) ! no. of vertices per element 76 nshl = lcblk(10,iblk) 77 mattyp = lcblk(7,iblk) 78 ndofl = lcblk(8,iblk) 79 nsymdl = lcblk(9,iblk) 80 npro = lcblk(1,iblk+1) - iel 81 ngauss = nint(lcsyst) 82c 83c.... compute and assemble diffusive flux vector residual, qres, 84c and lumped mass matrix, rmass 85 86 allocate (tmpshp(nshl,MAXQPT)) 87 allocate (tmpshgl(nsd,nshl,MAXQPT)) 88 89 tmpshp(1:nshl,:) = shp(lcsyst,1:nshl,:) 90 tmpshgl(:,1:nshl,:) = shgl(lcsyst,:,1:nshl,:) 91 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 103 enddo 104c 105c.... take care of periodic boundary conditions 106c 107 108 call qpbc( rmass, qres, iBC, iper, ilwork ) 109 110c 111 endif ! computation of global diffusive flux 112c 113c.... loop over element blocks to compute element residuals 114c 115c 116c.... initialize the arrays 117c 118 res = zero 119 rmes = zero 120c 121 if (iprec .ne. 0) BDiag = zero 122c 123c.... loop over the element-blocks 124c 125 do iblk = 1, nelblk 126c 127c.... set up the parameters 128c 129 nenl = lcblk(5,iblk) ! no. of vertices per element 130 iel = lcblk(1,iblk) 131 lelCat = lcblk(2,iblk) 132 lcsyst = lcblk(3,iblk) 133 iorder = lcblk(4,iblk) 134 nenl = lcblk(5,iblk) ! no. of vertices per element 135 nshl = lcblk(10,iblk) 136 mattyp = lcblk(7,iblk) 137 ndofl = lcblk(8,iblk) 138 nsymdl = lcblk(9,iblk) 139 npro = lcblk(1,iblk+1) - iel 140 inum = iel + npro - 1 141 ngauss = nint(lcsyst) 142c 143c.... compute and assemble the residuals and the preconditioner 144c 145 allocate (tmpshp(nshl,MAXQPT)) 146 allocate (tmpshgl(nsd,nshl,MAXQPT)) 147 148 tmpshp(1:nshl,:) = shp(lcsyst,1:nshl,:) 149 tmpshgl(:,1:nshl,:) = shgl(lcsyst,:,1:nshl,:) 150 151 call AsIMFG (y, ac, 152 & x, mxmudmi(iblk)%p, 153 & tmpshp, tmpshgl, 154 & mien(iblk)%p, 155 & mmat(iblk)%p, res, 156 & rmes, BDiag, 157 & qres, rerr ) 158 159 deallocate (tmpshp) 160 deallocate (tmpshgl) 161 162c 163c.... end of interior element loop 164c 165 enddo 166c 167c.... --------------------> boundary elements <-------------------- 168c 169c.... loop over the elements 170c 171 do iblk = 1, nelblb 172c 173c.... set up the parameters 174c 175c 176 iel = lcblkb(1,iblk) 177 lelCat = lcblkb(2,iblk) 178 lcsyst = lcblkb(3,iblk) 179 iorder = lcblkb(4,iblk) 180 nenl = lcblkb(5,iblk) ! no. of vertices per element 181 nenbl = lcblkb(6,iblk) ! no. of vertices per bdry. face 182 mattyp = lcblkb(7,iblk) 183 ndofl = lcblkb(8,iblk) 184 nshl = lcblkb(9,iblk) 185 nshlb = lcblkb(10,iblk) 186 npro = lcblkb(1,iblk+1) - iel 187 if(lcsyst.eq.3) lcsyst=nenbl 188 ngaussb = nintb(lcsyst) 189c 190c.... compute and assemble the residuals corresponding to the 191c boundary integral 192c 193 194 allocate (tmpshpb(nshl,MAXQPT)) 195 allocate (tmpshglb(nsd,nshl,MAXQPT)) 196 197 tmpshpb(1:nshl,:) = shpb(lcsyst,1:nshl,:) 198 tmpshglb(:,1:nshl,:) = shglb(lcsyst,:,1:nshl,:) 199 200 call AsBMFG (y, x, 201 & tmpshpb, tmpshglb, 202 & mienb(iblk)%p, mmatb(iblk)%p, 203 & miBCB(iblk)%p, mBCB(iblk)%p, 204 & res, rmes) 205 206 deallocate (tmpshpb) 207 deallocate (tmpshglb) 208c 209c.... end of boundary element loop 210c 211 enddo 212c 213 ttim(80) = ttim(80) + secs(0.0) 214c 215c.... --------------------> communications <------------------------- 216c 217 if((iabc==1)) then !are there any axisym bc's 218 call rotabc(res(1,2), iBC, 'in ') 219 call rotabc(rmes(1,2), iBC, 'in ') 220 endif 221 222 if (numpe > 1) then 223c 224 225 call commu (res , ilwork, nflow , 'in ') 226 227 call MPI_BARRIER (MPI_COMM_WORLD,ierr) 228 229 call commu (rmes , ilwork, nflow , 'in ') 230 if(iprec.ne.0) call commu(BDiag,ilwork, nflow*nflow, 'in ') 231 endif 232 233c 234c.... ----------------------> post processing <---------------------- 235c 236c.... satisfy the BCs on the residual and the modified residual 237c 238 call bc3Res (y, iBC, BC, res, iper, ilwork) 239 call bc3Res (y, iBC, BC, rmes, iper, ilwork) 240 241c 242c.... satisfy the BCs on the block-diagonal preconditioner 243c 244 if (iprec .ne. 0) then 245 call bc3BDg (y, iBC, BC, BDiag, iper, ilwork) 246 endif 247 248 249c 250c.... return 251c 252 call timer ('Back ') 253 return 254 end 255 256 257