opkdemo7.f 16 KB
 rtoy committed Jan 12, 2006 1 `````` program opkdemo7 `````` rtoy committed Aug 14, 2004 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 ``````c----------------------------------------------------------------------- c Demonstration program for the DLSODI package. c This is the version of 14 June 2001. c c This version is in double precision. c C this program solves a semi-discretized form of the Burgers equation, c c u = -(u*u/2) + eta * u c t x xx c c for a = -1 .le. x .le. 1 = b, t .ge. 0. c Here eta = 0.05. c Boundary conditions: u(-1,t) = u(1,t) = 0. c Initial profile: square wave c u(0,x) = 0 for 1/2 .lt. abs(x) .le. 1 c u(0,x) = 1/2 for abs(x) = 1/2 c u(0,x) = 1 for 0 .le. abs(x) .lt. 1/2 c c An ODE system is generated by a simplified Galerkin treatment c of the spatial variable x. c c Reference: c R. C. Y. Chin, G. W. Hedstrom, and K. E. Karlsson, c A Simplified Galerkin Method for Hyperbolic Equations, c Math. Comp., vol. 33, no. 146 (April 1979), pp. 647-658. c c The problem is run with the DLSODI package with a 10-point mesh c and a 100-point mesh. In each case, it is run with two tolerances c and for various appropriate values of the method flag mf. c Output is on unit lout, set to 6 in a data statement below. c----------------------------------------------------------------------- external res, addabd, addafl, jacbd, jacfl integer i, io, istate, itol, iwork, j, 1 lout, liw, lrw, meth, miter, mf, ml, mu, 2 n, nout, npts, nerr, 3 nptsm1, n14, n34, n14m1, n14p1, n34m1, n34p1 integer nm1 double precision a, b, eta, delta, 1 zero, fourth, half, one, hun, 2 t, tout, tlast, tinit, errfac, 3 atol, rtol, rwork, y, ydoti, elkup double precision eodsq, r4d `````` rtoy committed Aug 17, 2004 45 `````` dimension n(1) `````` rtoy committed Aug 14, 2004 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 `````` dimension y(99), ydoti(99), tout(4), atol(2), rtol(2) dimension rwork(2002), iwork(125) c Pass problem parameters in the Common block test1. common /test1/ r4d, eodsq, nm1 c c Set problem parameters and run parameters data eta/0.05d0/, a/-1.0d0/, b/1.0d0/ data zero/0.0d0/, fourth/0.25d0/, half/.5d0/, one/1.0d0/, 1 hun/100.0d0/ data tinit/0.0d0/, tlast/0.4d0/ data tout/.10d0,.20d0,.30d0,.40d0/ data ml/1/, mu/1/, lout/6/ data nout/4/, lrw/2002/, liw/125/ data itol/1/, rtol/1.0d-3, 1.0d-6/, atol/1.0d-3, 1.0d-6/ c iwork(1) = ml iwork(2) = mu nerr = 0 c c Loop over two values of npts. do 300 npts = 10, 100, 90 c c Compute the mesh width delta and other parameters. delta = (b - a)/npts r4d = fourth/delta eodsq = eta/delta**2 nptsm1 = npts - 1 n14 = npts/4 n34 = 3 * n14 n14m1 = n14 - 1 n14p1 = n14m1 + 2 n34m1 = n34 - 1 n34p1 = n34m1 + 2 `````` rtoy committed Aug 17, 2004 79 80 `````` n(1) = nptsm1 nm1 = n(1) - 1 `````` rtoy committed Aug 14, 2004 81 82 83 84 85 86 87 88 89 90 91 92 93 94 ``````c c Set the initial profile (for output purposes only). c do 10 i = 1,n14m1 10 y(i) = zero y(n14) = half do 20 i = n14p1,n34m1 20 y(i) = one y(n34) = half do 30 i = n34p1,nptsm1 30 y(i) = zero c if (npts .gt. 10) write (lout,1010) write (lout,1000) `````` rtoy committed Aug 17, 2004 95 96 `````` write (lout,1100) eta,a,b,tinit,tlast,ml,mu,n(1) write (lout,1200) zero, (y(i), i=1,n(1)), zero `````` rtoy committed Aug 14, 2004 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 ``````c c The j loop is over error tolerances. c do 200 j = 1,2 c c Loop over method flag loop (for demonstration). c do 100 meth = 1,2 do 100 miter = 1,5 if (miter .eq. 3) go to 100 if (miter .le. 2 .and. npts .gt. 10) go to 100 if (miter .eq. 5 .and. npts .lt. 100) go to 100 mf = 10*meth + miter c c Set the initial profile. c do 40 i = 1,n14m1 40 y(i) = zero y(n14) = half do 50 i = n14p1,n34m1 50 y(i) = one y(n34) = half do 60 i = n34p1,nptsm1 60 y(i) = zero c t = tinit istate = 0 c write (lout,1500) rtol(j), atol(j), mf, npts c c Output loop for each case c do 80 io = 1,nout c c call DLSODI if (miter .le. 2) call dlsodi (res, addafl, jacfl, n, y, 1 ydoti, t, tout(io), itol, rtol(j), atol(j), 2 1, istate, 0, rwork, lrw, iwork, liw, mf) if (miter .ge. 4) call dlsodi (res, addabd, jacbd, n, y, 1 ydoti, t, tout(io), itol, rtol(j), atol(j), 2 1, istate, 0, rwork, lrw, iwork, liw, mf) `````` rtoy committed Aug 17, 2004 138 `````` write (lout,2000) t, rwork(11), iwork(14),(y(i), i=1,n(1)) `````` rtoy committed Aug 14, 2004 139 140 141 142 143 144 145 146 147 148 149 150 151 152 ``````c c If istate is not 2 on return, print message and loop. if (istate .ne. 2) then write (lout,4000) mf, t, istate nerr = nerr + 1 go to 100 endif c 80 continue c write (lout,3000) mf, iwork(11), iwork(12), iwork(13), 1 iwork(17), iwork(18) c c Estimate final error and print result. `````` rtoy committed Aug 17, 2004 153 154 `````` itemp = n(1) errfac = elkup( itemp, y, rwork(21), itol, rtol(j), atol(j) ) `````` rtoy committed Aug 14, 2004 155 156 157 158 159 160 161 162 163 164 165 `````` if (errfac .gt. hun) then write (lout,5001) errfac nerr = nerr + 1 else write (lout,5000) errfac endif 100 continue 200 continue 300 continue c write (lout,6000) nerr `````` rtoy committed Jan 12, 2006 166 ``````c stop `````` rtoy committed Aug 14, 2004 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 ``````c 1000 format(20x,' Demonstration Problem for DLSODI') 1010 format(///80('*')///) 1100 format(/10x,' Simplified Galerkin Solution of Burgers Equation'// 1 13x,'Diffusion coefficient is eta =',d10.2/ 2 13x,'Uniform mesh on interval',d12.3,' to ',d12.3/ 3 13x,'Zero boundary conditions'/ 4 13x,'Time limits: t0 = ',d12.5,' tlast = ',d12.5/ 5 13x,'Half-bandwidths ml = ',i2,' mu = ',i2/ 6 13x,'System size neq = ',i3/) c 1200 format('Initial profile:'/17(6d12.4/)) c 1500 format(///80('-')///'Run with rtol =',d12.2,' atol =',d12.2, 1 ' mf =',i3,' npts =',i4,':'//) c 2000 format('Output for time t = ',d12.5,' current h =', 1 d12.5,' current order =',i2,':'/17(6d12.4/)) c 3000 format(//'Final statistics for mf = ',i2,':'/ 1 i4,' steps,',i5,' res,',i4,' Jacobians,', 2 ' rwork size =',i6,', iwork size =',i6) c 4000 format(///80('*')//20x,'Final time reached for mf = ',i2, 1 ' was t = ',d12.5/25x,'at which istate = ',i2////80('*')) 5000 format(' Final output is correct to within ',d8.1, 1 ' times local error tolerance') 5001 format(' Final output is wrong by ',d8.1, 1 ' times local error tolerance') 6000 format(//80('*')// 1 'Run completed. Number of errors encountered =',i3) c c end of main program for the DLSODI demonstration problem. end subroutine gfun (n, t, y, g) c This subroutine computes the right-hand side function g(y,t). c It uses r4d = 1/(4*delta), eodsq = eta/delta**2, and nm1 = n - 1 c from the Common block test1. c integer i, n, nm1 double precision t, y, g, r4d, eodsq, two dimension g(n), y(n) common /test1/ r4d, eodsq, nm1 data two/2.0d0/ c g(1) = -r4d*y(2)**2 + eodsq*(y(2) - two*y(1)) c do 20 i = 2,nm1 g(i) = r4d*(y(i-1)**2 - y(i+1)**2) 1 + eodsq*(y(i+1) - two*y(i) + y(i-1)) 20 continue c g(n) = r4d*y(nm1)**2 + eodsq*(y(nm1) - two*y(n)) c return c end of subroutine gfun for the DLSODI demonstration problem. end `````` rtoy committed Aug 17, 2004 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 `````` subroutine res (n, t, y, v, r, ires) c This subroutine computes the residual vector c r = g(t,y) - A(t,y)*v . c It uses nm1 = n - 1 from Common. c If ires = -1, only g(t,y) is returned in r, since A(t,y) does c not depend on y. c integer i, ires, n, nm1 double precision t, y, v, r, r4d, eodsq, one, four, six, 1 fact1, fact4 dimension y(*), v(*), r(*) dimension n(*) integer itemp common /test1/ r4d, eodsq, nm1 data one /1.0d0/, four /4.0d0/, six /6.0d0/ c itemp = n(1) call gfun (itemp, t, y, r) if (ires .eq. -1) return c fact1 = one/six fact4 = four/six r(1) = r(1) - (fact4*v(1) + fact1*v(2)) do 10 i = 2, nm1 10 r(i) = r(i) - (fact1*v(i-1) + fact4*v(i) + fact1*v(i+1)) r(n(1)) = r(n(1)) - (fact1*v(nm1) + fact4*v(n(1))) return c end of subroutine res for the DLSODI demonstration problem. end `````` rtoy committed Aug 14, 2004 256 257 258 259 260 261 262 `````` subroutine addabd (n, t, y, ml, mu, pa, m0) c This subroutine computes the matrix A in band form, adds it to pa, c and returns the sum in pa. The matrix A is tridiagonal, of order n, c with nonzero elements (reading across) of 1/6, 4/6, 1/6. c integer i, n, m0, ml, mu, mup1, mup2 double precision t, y, pa, fact1, fact4, one, four, six `````` rtoy committed Aug 17, 2004 263 264 `````` dimension y(*), pa(m0,*) dimension n(*) `````` rtoy committed Aug 14, 2004 265 266 267 268 269 270 271 272 273 `````` data one/1.0d0/, four/4.0d0/, six/6.0d0/ c c Set the pointers. mup1 = mu + 1 mup2 = mu + 2 c Compute the elements of A. fact1 = one/six fact4 = four/six c Add the matrix A to the matrix pa (banded). `````` rtoy committed Aug 17, 2004 274 `````` do 10 i = 1,n(1) `````` rtoy committed Aug 14, 2004 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 `````` pa(mu,i) = pa(mu,i) + fact1 pa(mup1,i) = pa(mup1,i) + fact4 pa(mup2,i) = pa(mup2,i) + fact1 10 continue return c end of subroutine addabd for the DLSODI demonstration problem. end subroutine addafl (n, t, y, ml, mu, pa, m0) c This subroutine computes the matrix A in full form, adds it to c pa, and returns the sum in pa. c It uses nm1 = n - 1 from Common. c The matrix A is tridiagonal, of order n, with nonzero elements c (reading across) of 1/6, 4/6, 1/6. c integer i, n, m0, ml, mu, nm1 double precision t, y, pa, r4d, eodsq, one, four, six, 1 fact1, fact4 `````` rtoy committed Aug 17, 2004 293 294 `````` dimension y(*), pa(m0,8) dimension n(1) `````` rtoy committed Aug 14, 2004 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 `````` common /test1/ r4d, eodsq, nm1 data one/1.0d0/, four/4.0d0/, six/6.0d0/ c c Compute the elements of A. fact1 = one/six fact4 = four/six c c Add the matrix A to the matrix pa (full). c do 110 i = 2, nm1 pa(i,i+1) = pa(i,i+1) + fact1 pa(i,i) = pa(i,i) + fact4 pa(i,i-1) = pa(i,i-1) + fact1 110 continue pa(1,2) = pa(1,2) + fact1 pa(1,1) = pa(1,1) + fact4 `````` rtoy committed Aug 17, 2004 311 312 `````` pa(n(1),n(1)) = pa(n(1),n(1)) + fact4 pa(n(1),nm1) = pa(n(1),nm1) + fact1 `````` rtoy committed Aug 14, 2004 313 314 315 316 317 318 319 320 321 322 323 324 325 326 `````` return c end of subroutine addafl for the DLSODI demonstration problem. end subroutine jacbd (n, t, y, s, ml, mu, pa, m0) c This subroutine computes the Jacobian dg/dy = d(g-a*s)/dy c and stores elements c i j c dg /dy in pa(i-j+mu+1,j) in band matrix format. c It uses r4d = 1/(4*delta), eodsq = eta/delta**2, and nm1 = n - 1 c from the Common block test1. c integer i, n, m0, ml, mu, mup1, mup2, nm1 double precision t, y, s, pa, diag, r4d, eodsq, two, r2d `````` rtoy committed Aug 17, 2004 327 328 `````` dimension y(*), s(*), pa(m0,*) dimension n(*) `````` rtoy committed Aug 14, 2004 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 `````` common /test1/ r4d, eodsq, nm1 data two/2.0d0/ c mup1 = mu + 1 mup2 = mu + 2 diag = -two*eodsq r2d = two*r4d c 1 1 c Compute and store dg /dy pa(mup1,1) = diag c c 1 2 c Compute and store dg /dy pa(mu,2) = -r2d*y(2) + eodsq c do 20 i = 2,nm1 c c i i-1 c Compute and store dg /dy pa(mup2,i-1) = r2d*y(i-1) + eodsq c c i i c Compute and store dg /dy pa(mup1,i) = diag c c i i+1 c Compute and store dg /dy pa(mu,i+1) = -r2d*y(i+1) + eodsq 20 continue c c n n-1 c Compute and store dg /dy pa(mup2,nm1) = r2d*y(nm1) + eodsq c c n n c Compute and store dg /dy `````` rtoy committed Aug 17, 2004 365 `````` pa(mup1,n(1)) = diag `````` rtoy committed Aug 14, 2004 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 ``````c return c end of subroutine jacbd for the DLSODI demonstration problem. end subroutine jacfl (n, t, y, s, ml, mu, pa, m0) c This subroutine computes the Jacobian dg/dy = d(g-a*s)/dy c and stores elements c i j c dg /dy in pa(i,j) in full matrix format. c It uses r4d = 1/(4*delta), eodsq = eta/delta**2, and nm1 = n - 1 c from the Common block test1. c integer i, n, m0, ml, mu, nm1 double precision t, y, s, pa, diag, r4d, eodsq, two, r2d `````` rtoy committed Aug 17, 2004 381 382 `````` dimension y(*), s(8), pa(m0,*) dimension n(*) `````` rtoy committed Aug 14, 2004 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 `````` common /test1/ r4d, eodsq, nm1 data two/2.0d0/ c diag = -two*eodsq r2d = two*r4d c c 1 1 c Compute and store dg /dy pa(1,1) = diag c c 1 2 c Compute and store dg /dy pa(1,2) = -r2d*y(2) + eodsq c do 120 i = 2,nm1 c c i i-1 c Compute and store dg /dy pa(i,i-1) = r2d*y(i-1) + eodsq c c i i c Compute and store dg /dy pa(i,i) = diag c c i i+1 c Compute and store dg /dy pa(i,i+1) = -r2d*y(i+1) + eodsq 120 continue c c n n-1 c Compute and store dg /dy `````` rtoy committed Aug 17, 2004 414 `````` pa(n(1),nm1) = r2d*y(nm1) + eodsq `````` rtoy committed Aug 14, 2004 415 416 417 ``````c c n n c Compute and store dg /dy `````` rtoy committed Aug 17, 2004 418 `````` pa(n(1),n(1)) = diag `````` rtoy committed Aug 14, 2004 419 420 421 422 423 424 425 426 427 428 429 430 431 ``````c return c end of subroutine jacfl for the DLSODI demonstration problem. end double precision function elkup (n, y, ewt, itol, rtol, atol) c This routine looks up approximately correct values of y at t = 0.4, c ytrue = y9 or y99 depending on whether n = 9 or 99. These were c obtained by running DLSODI with very tight tolerances. c The returned value is c elkup = norm of ( y - ytrue ) / ( rtol*abs(ytrue) + atol ). c integer n, itol, i `````` rtoy committed May 04, 2006 432 `````` double precision y, ewt, rtol, atol, y9, y99 `````` rtoy committed Aug 17, 2004 433 434 ``````c double precision y99a, y99b, y99c, c 1 y99d, y99e, y99f, y99g, dvnorm `````` rtoy committed Aug 14, 2004 435 `````` dimension y(n), ewt(n), y9(9), y99(99) `````` rtoy committed Aug 17, 2004 436 437 438 439 440 441 `````` dimension rtol(*), atol(*) c dimension y99a(16), y99b(16), y99c(16), y99d(16), y99e(16), c 1 y99f(16), y99g(3) c equivalence (y99a(1),y99(1)), (y99b(1),y99(17)), c 1 (y99c(1),y99(33)), (y99d(1),y99(49)), (y99e(1),y99(65)), c 1 (y99f(1),y99(81)), (y99g(1),y99(97)) `````` rtoy committed Aug 14, 2004 442 443 444 445 `````` data y9 / 1 1.07001457d-01, 2.77432492d-01, 5.02444616d-01, 7.21037157d-01, 1 9.01670441d-01, 8.88832048d-01, 4.96572850d-01, 9.46924362d-02, 1-6.90855199d-03 / `````` rtoy committed Aug 17, 2004 446 `````` data y99 / `````` rtoy committed Aug 14, 2004 447 448 449 `````` 1 2.05114384d-03, 4.19527452d-03, 6.52533872d-03, 9.13412751d-03, 1 1.21140191d-02, 1.55565301d-02, 1.95516488d-02, 2.41869487d-02, 1 2.95465081d-02, 3.57096839d-02, 4.27498067d-02, 5.07328729d-02, `````` rtoy committed Aug 17, 2004 450 `````` 1 5.97163151d-02, 6.97479236d-02, 8.08649804d-02, 9.30936515d-02, `````` rtoy committed Aug 14, 2004 451 452 453 `````` 1 1.06448659d-01, 1.20933239d-01, 1.36539367d-01, 1.53248227d-01, 1 1.71030869d-01, 1.89849031d-01, 2.09656044d-01, 2.30397804d-01, 1 2.52013749d-01, 2.74437805d-01, 2.97599285d-01, 3.21423708d-01, `````` rtoy committed Aug 17, 2004 454 `````` 1 3.45833531d-01, 3.70748792d-01, 3.96087655d-01, 4.21766871d-01, `````` rtoy committed Aug 14, 2004 455 456 457 `````` 1 4.47702161d-01, 4.73808532d-01, 5.00000546d-01, 5.26192549d-01, 1 5.52298887d-01, 5.78234121d-01, 6.03913258d-01, 6.29252015d-01, 1 6.54167141d-01, 6.78576790d-01, 7.02400987d-01, 7.25562165d-01, `````` rtoy committed Aug 17, 2004 458 `````` 1 7.47985803d-01, 7.69601151d-01, 7.90342031d-01, 8.10147715d-01, `````` rtoy committed Aug 14, 2004 459 460 461 `````` 1 8.28963844d-01, 8.46743353d-01, 8.63447369d-01, 8.79046021d-01, 1 8.93519106d-01, 9.06856541d-01, 9.19058529d-01, 9.30135374d-01, 1 9.40106872d-01, 9.49001208d-01, 9.56853318d-01, 9.63702661d-01, `````` rtoy committed Aug 17, 2004 462 `````` 1 9.69590361d-01, 9.74555682d-01, 9.78631814d-01, 9.81840924d-01, `````` rtoy committed Aug 14, 2004 463 464 465 `````` 1 9.84188430d-01, 9.85656465d-01, 9.86196496d-01, 9.85721098d-01, 1 9.84094964d-01, 9.81125395d-01, 9.76552747d-01, 9.70041743d-01, 1 9.61175143d-01, 9.49452051d-01, 9.34294085d-01, 9.15063568d-01, `````` rtoy committed Aug 17, 2004 466 `````` 1 8.91098383d-01, 8.61767660d-01, 8.26550038d-01, 7.85131249d-01, `````` rtoy committed Aug 14, 2004 467 468 469 `````` 1 7.37510044d-01, 6.84092540d-01, 6.25748369d-01, 5.63802368d-01, 1 4.99946558d-01, 4.36077986d-01, 3.74091566d-01, 3.15672765d-01, 1 2.62134958d-01, 2.14330497d-01, 1.72640946d-01, 1.37031155d-01, `````` rtoy committed Aug 17, 2004 470 471 `````` 1 1.07140815d-01, 8.23867920d-02, 6.20562432d-02, 4.53794321d-02, 1 3.15789227d-02, 1.98968820d-02, 9.60472135d-03 / `````` rtoy committed Aug 14, 2004 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 ``````c if (n .eq. 99) go to 99 c c Compute local error tolerance using correct y (n = 9). c call dewset( n, itol, rtol, atol, y9, ewt ) c c Invert ewt and replace y by the error, y - ytrue. c do 20 i = 1, 9 ewt(i) = 1.0d0/ewt(i) 20 y(i) = y(i) - y9(i) go to 200 c c Compute local error tolerance using correct y (n = 99). c 99 call dewset( n, itol, rtol, atol, y99, ewt ) c c Invert ewt and replace y by the error, y - ytrue. c do 120 i = 1, 99 ewt(i) = 1.0d0/ewt(i) 120 y(i) = y(i) - y99(i) c c Find weighted norm of the error and return. c 200 elkup = dvnorm (n, y, ewt) return c end of function elkup for the DLSODI demonstration program. end``````