1 subroutine itrFDI (ypre, y, ac, x, 2 & rmes, uBrg, BDiag, 3 & iBC, BC, iper, 4 & ilwork, shp, shgl, 5 & shpb, shglb) 6c 7c---------------------------------------------------------------------- 8c 9c This subroutine computes the "optimum" finite difference 10c interval eps for the forward difference scheme 11c 12c Rmod(y + eps u) - Rmod(y) 13c --------------------------- 14c eps 15c 16c where u is the step and Rmod is the modified residual. 17c 18c Note: A good theoretical reference is 'Practical Optimization' by 19c P.E. Gill, W. Murray and M.H. Wright [1981]. 20c 21c input: 22c y (nshg,ndof) : Y-variables 23c ypre (nshg,ndof) : preconditioned Y-variables 24c x (numnp,nsd) : node coordinates 25c rmes (nshg,nflow) : modified residual 26c uBrg (nshg,nflow) : step 27c BDiag (nshg,nflow,nflow) : block-diagonal preconditioner 28c iBC (nshg) : BC codes 29c BC (nshg,ndofBC) : BC constraint parameters 30c 31c 32c 33c Zdenek Johan, Winter 1989. 34c Zdenek Johan, Winter 1991. (Fortran 90) 35c---------------------------------------------------------------------- 36c 37 include "common.h" 38c 39 dimension y(nshg,ndof), ypre(nshg,nflow), 40 & x(numnp,nsd), ac(nshg,ndof), 41 & rmes(nshg,nflow), 42 & uBrg(nshg,nflow), BDiag(nshg,nflow,nflow), 43 & iBC(nshg), BC(nshg,ndofBC), 44 & ilwork(nlwork), 45 & iper(nshg) 46c 47 dimension ytmp(nshg,nflow), rtmp(nshg,nflow) 48 dimension tmpy(nshg,ndof) 49c 50 dimension shp(MAXTOP,maxsh,MAXQPT), 51 & shgl(MAXTOP,nsd,maxsh,MAXQPT), 52 & shpb(MAXTOP,maxsh,MAXQPT), 53 & shglb(MAXTOP,nsd,maxsh,MAXQPT) 54c 55c.... compute the accuracy (cancellation error) -> epsA 56c 57 rtmp = zero 58c 59 ytmp = ypre 60c 61c call yshuffle(ytmp, 'new2old ') 62c 63 call i3LU (BDiag, ytmp, 'backward') 64c 65 call yshuffle(ytmp, 'old2new ') 66c 67 iabres = 1 68c 69 call itrRes (ytmp, y, 70 & x, shp, 71 & shgl, iBC, 72 & BC, shpb, 73 & shglb, rtmp, 74 & iper, ilwork, 75 & ac) 76c 77 iabres = 0 78c 79 call i3LU (BDiag, rtmp, 'forward ') 80c 81 rtmp = rtmp**2 82 call sumgat (rtmp, nflow, summed, ilwork) 83 epsA = (epsM**2) * sqrt(summed) 84c 85c.... compute the norm of the second derivative (truncation error) 86c 87c.... set interval 88c 89 epsSD = sqrt(epsM) 90c 91c.... compute the first residual 92c 93 rtmp = zero 94c 95c call yshuffle(ypre, 'new2old ') 96c 97 ytmp = ypre + epsSD * uBrg 98c 99 call i3LU (BDiag, ytmp, 'backward') 100c 101 call yshuffle(ytmp, 'old2new ') 102c 103 call itrRes (ytmp, y, 104 & x, shp, 105 & shgl, iBC, 106 & BC, shpb, 107 & shglb, rtmp, 108 & iper, ilwork) 109c 110c.... compute the second residual and add it to the first one 111c 112 ytmp = ypre - epsSD * uBrg 113c 114c call yshuffle(ypre, 'old2new ') 115c 116 call i3LU (BDiag, ytmp, 'backward') 117c 118 call yshuffle(ytmp, 'old2new ') 119 call itrRes (ytmp, y, 120 & x, shp, 121 & shgl, iBC, 122 & BC, shpb, 123 & shglb, rtmp, 124 & iper, ilwork) 125c 126 call i3LU (BDiag, rtmp, 'forward ') 127c 128c.... compute the second derivative and its norm 129c 130 rtmp = (( rtmp - two * rmes ) / epsM)**2 131c 132 call sumgat (rtmp, nflow, summed, ilwork) 133 SDnrm = sqrt(summed) 134c 135c.... compute the 'optimum' interval 136c 137 eGMRES = two * sqrt( epsA / SDnrm ) 138c 139c.... flop count 140c 141 flops = flops + 10*nflow*nshg+3*nshg 142c 143c.... end 144c 145 return 146 end 147 148 149 subroutine itrFDISclr (y, ypre, x, 150 & rmes, uBrg, BDiag, 151 & iBC, BC, engBC, iper, 152 & ilwork) 153c 154c---------------------------------------------------------------------- 155c 156c This subroutine computes the "optimum" finite difference 157c interval eps for the forward difference scheme 158c 159c Rmod(y + eps u) - Rmod(y) 160c --------------------------- 161c eps 162c 163c where u is the step and Rmod is the modified residual. 164c 165c Note: A good theoretical reference is 'Practical Optimization' by 166c P.E. Gill, W. Murray and M.H. Wright [1981]. 167c 168c input: 169c y (nshg,ndof) : Y-variables 170c ypre (nshg,ndof) : preconditioned Y-variables 171c x (numnp,nsd) : node coordinates 172c rmes (nshg,nflow) : modified residual 173c uBrg (nshg,nflow) : step 174c BDiag (nshg,nflow,nflow) : block-diagonal preconditioner 175c iBC (nshg) : BC codes 176c BC (nshg,ndofBC) : BC constraint parameters 177c engBC (nshg) : energy for BC on density or pressure 178c 179c 180c Zdenek Johan, Winter 1989. 181c Zdenek Johan, Winter 1991. (Fortran 90) 182c---------------------------------------------------------------------- 183c 184 include "common.h" 185c 186 dimension y(nshg,ndof), ypre(nshg,ndof), 187 & x(numnp,nsd), 188 & rmes(nshg,nflow), 189 & uBrg(nshg,nflow), BDiag(nshg,nflow,nflow), 190 & iBC(nshg), BC(nshg,ndofBC), 191 & engBC(nshg), ilwork(nlwork), 192 & iper(nshg) 193c 194 dimension ytmp(nshg,ndof), rtmp(nshg,nflow) 195c 196c.... compute the accuracy (cancellation error) -> epsA 197c 198 rtmp = zero 199c 200 ytmp = ypre 201c 202c call tnanq(ytmp,ndof,"ytmp ") 203 iabres = 1 204c 205 call itrRes (ytmp, y, 206 & x, a(mshp), 207 & a(mshgl), a(mwght), iBC, 208 & BC, engBC, a(mshpb), 209 & a(mshglb), a(mwghtb), rtmp, 210 & iper, ilwork) 211c 212 iabres = 0 213c 214 rtmp = rtmp**2 215 call sumgat (rtmp, nflow, summed, ilwork) 216 epsA = (epsM**2) * sqrt(summed) 217c 218c.... compute the norm of the second derivative (truncation error) 219c 220c.... set interval 221c 222 epsSD = sqrt(epsM) 223c 224c.... compute the first residual 225c 226 rtmp = zero 227c 228 ytmp = ypre + epsSD * uBrg 229c 230 call itrRes (ytmp, y, 231 & x, a(mshp), 232 & a(mshgl), a(mwght), iBC, 233 & BC, engBC, a(mshpb), 234 & a(mshglb), a(mwghtb), rtmp, 235 & iper, ilwork) 236c 237c.... compute the second residual and add it to the first one 238c 239 ytmp = ypre - epsSD * uBrg 240c 241 call itrRes (ytmp, y, 242 & x, a(mshp), 243 & a(mshgl), a(mwght), iBC, 244 & BC, engBC, a(mshpb), 245 & a(mshglb), a(mwghtb), rtmp, 246 & iper, ilwork) 247c 248c.... compute the second derivative and its norm 249c 250 rtmp = (( rtmp - two * rmes ) / epsM)**2 251c 252 call sumgat (rtmp, nflow, summed, ilwork) 253 SDnrm = sqrt(summed) 254c 255c.... compute the 'optimum' interval 256c 257 eGMRES = two * sqrt( epsA / SDnrm ) 258c 259c.... flop count 260c 261 flops = flops + 10*nflow*nshg+3*nshg 262c 263c.... end 264c 265 return 266 end 267