1c----------------------------------------------------------------------- 2c Demonstration program for the DLSODIS package. 3c This is the version of 14 June 2001. 4c 5c This version is in double precision. 6c 7c This program solves a semi-discretized form of the Burgers equation, 8c 9c u = -(u*u/2) + eta * u 10c t x xx 11c 12c for -1 .le. x .le. 1, t .ge. 0. 13c Here eta = 0.05. 14c Boundary conditions: u(-1,t) = u(1,t) and du/dx(-1,t) = du/dx(1,t). 15c Initial profile: square wave 16c u(0,x) = 0 for 1/2 .lt. abs(x) .le. 1 17c u(0,x) = 1/2 for abs(x) = 1/2 18c u(0,x) = 1 for 0 .le. abs(x) .lt. 1/2 19c 20c An ODE system is generated by a simplified Galerkin treatment 21c of the spatial variable x. 22c 23c Reference: 24c R. C. Y. Chin, G. W. Hedstrom, and K. E. Karlsson, 25c A Simplified Galerkin Method for Hyperbolic Equations, 26c Math. Comp., vol. 33, no. 146 (April 1979), pp. 647-658. 27c 28c The problem is run with the DLSODIS package with a 12-node mesh, 29c for various appropriate values of the method flag mf. 30c Output is on unit lout, set to 6 in a data statement below. 31c 32c Problem specific data: 33c npts = number of unknowns (npts = 0 mod 4) 34c nnz = number of non-zeros in Jacobian before fill in 35c nnza = number of non-zeros in Jacobian after fill in 36c lrwk = length of real work array (taking into account fill in) 37c liwk = length of integer work array 38c ipia = pointer to ia in iw (ia(j) = iw(ipia+j-1) 39c ipja = pointer to ja in iw (ja(j) = iw(ipja+j-1) 40c ipic = pointer to ic in iw array (ic(j) = iw(ipic+j-1)) 41c ipjc = pointer to jc in iw array (jc(j) = iw(ipjc+j-1)) 42c----------------------------------------------------------------------- 43 parameter (npts = 12, nnz = 3*npts, nnza = 5*npts) 44 parameter (lrwk = 20+3*nnza+28*npts) 45 parameter (liwk = 32+2*nnza+2*npts) 46 parameter (ipia = 31, ipja = 31+npts+1) 47 parameter (ll = 32+npts+nnz, ipic = ll,ipjc = ll+npts+1) 48c 49 external res, addasp, jacsp 50 integer i, io, istate, itol, iw, j, lout, lrw, liw, lyh, 51 1 meth, miter, mf, moss, n, nout, nerr, n14, n34 52 double precision eta, delta, r4d, eodsq, a, b, t, tout, tlast, 53 1 tinit, errfac 54 double precision atol, rtol, rw, y, ydoti, elkup 55 double precision zero, fourth, half, one, hun 56 dimension y(npts), ydoti(npts), tout(4), atol(2), rtol(2) 57 dimension rw(lrwk), iw(liwk) 58c Pass problem parameters in the Common block test1. 59 common /test1/ r4d, eodsq 60c 61c Set problem parameters and run parameters 62 data eta/0.05d0/, a/-1.0d0/, b/1.0d0/ 63 data zero/0.0d0/, fourth/0.25d0/, half/0.5d0/, one/1.0d0/, 64 1 hun/100.0d0/ 65 data tinit/0.0d0/, tlast/0.4d0/ 66 data tout/0.10d0, 0.20d0, 0.30d0, 0.40d0/ 67 data lout/6/, nout/4/ 68 data itol/1/, rtol/1.0d-3, 1.0d-6/, atol/1.0d-3, 1.0d-6/ 69c 70 nerr = 0 71 lrw = lrwk 72 liw = liwk 73c 74c Compute the mesh width delta and other parameters. 75 delta = (b - a)/npts 76 r4d = fourth/delta 77 eodsq = eta/delta**2 78 n14 = npts/4 + 1 79 n34 = 3 * (npts/4) + 1 80 n = npts 81c 82c Set the initial profile (for output purposes only). 83 do 10 i = 1,n 84 10 y(i) = zero 85 y(n14) = half 86 do 20 i = n14+1,n34-1 87 20 y(i) = one 88 y(n34) = half 89c 90 write(lout,1000) 91 write(lout,1100) eta,a,b,tinit,tlast,n 92 write(lout,1200) (y(i), i=1,n) 93c 94c Set the initial sparse data structures for coefficient matrix A 95c and the Jacobian matrix C 96 call struct(iw(ipia), iw(ipja), iw(ipic), iw(ipjc), n) 97c 98c The j loop is over error tolerances. 99 do 200 j = 1,2 100c 101c This method flag loop is for demonstration only. 102 do 200 moss = 0,4 103 do 100 meth = 1,2 104 do 100 miter = 1,2 105 35 mf = 100*moss + 10*meth + miter 106c 107c Set the initial profile. 108 do 40 i = 1,n 109 40 y(i) = zero 110 y(n14) = half 111 do 50 i = n14+1,n34-1 112 50 y(i) = one 113 y(n34) = half 114c 115 t = tinit 116 istate = 0 117c 118 write(lout,1500) itol, rtol(j), atol(j), mf 119c 120c output loop for each case 121 do 80 io = 1,nout 122c 123c Call DLSODIS and print results. 124 call dlsodis (res, addasp, jacsp, n, y, ydoti, t, 125 1 tout(io), itol, rtol(j), atol(j), 1, istate, 0, 126 2 rw, lrw, iw, liw, mf) 127 write(lout,2000) t, rw(11), iw(14), (y(i), i=1,n) 128c 129c If istate is not 2 on return, print message and go to next case. 130 if (istate .ne. 2) then 131 write(lout,4000) mf, t, istate 132 nerr = nerr + 1 133 go to 100 134 endif 135 80 continue 136 write(lout,3000) mf, iw(11), iw(12), iw(13), 137 1 iw(17), iw(18), iw(20), iw(21) 138c 139c Estimate final error and print result. 140 lyh = iw(22) 141 errfac = elkup(n, y, rw(lyh), itol, rtol(j), atol(j), lout) 142 if (errfac .lt. hun) then 143 write(lout,5000) errfac 144 else 145 write(lout,5001) errfac 146 nerr = nerr + 1 147 endif 148 100 continue 149 200 continue 150c 151 write(lout,6000) nerr 152 stop 153c 154 1000 format(20x,' Demonstration Program for DLSODIS' ) 155 1100 format(//10x,'-- Simplified Galerkin solution of ', 156 1 'Burgers equation --'/// 157 1 13x,'Diffusion coefficient is eta =',d10.2/ 158 1 13x,'Uniform mesh on interval',d12.3,' to ',d12.3/ 159 2 13x,'Periodic boundary conditions'/ 160 2 13x,'Initial data are as follows:'//20x,'t0 = ',d12.5/ 161 2 20x,'tlast = ',d12.5/20x,'n = ',i3//) 162c 163 1200 format(/'Initial profile:',/20(6d12.4/)) 164c 165 1500 format(///85('*')///'Run with itol =',i2,' rtol =',d12.2, 166 1 ' atol =',d12.2,' mf = ',i3//) 167c 168 2000 format(' Output for time t =',d12.5,' current h =',d12.5, 169 1 ' current order =',i2/20(6d12.4/)) 170c 171 3000 format(/'Final statistics for mf = ',i3,': ', 172 1 i5,' steps,',i6,' res,',i6,' Jacobians,'/ 173 2 20x,' rw size =',i6,', iw size =',i6/ 174 3 20x,i4,' extra res for each jac,',i4,' decomps') 175c 176 4000 format(/'Final time reached for mf = ',i3, 177 1 ' was t = ',d12.5/25x,'at which istate = ',i2//) 178 5000 format('Final output is correct to within ',d9.2, 179 1 ' times local error tolerance.'/) 180 5001 format('Final output is wrong by ',d9.2, 181 1 ' times local error tolerance.'/) 182 6000 format(///85('*')// 183 1 'Run completed: number of errors encountered =',i3) 184c 185c end of main program for the DLSODIS demonstration program 186 end 187 188 subroutine struct(ia, ja, ic, jc, n) 189c This subroutine computes the initial sparse data structure of 190c the mass (ia,ja) and Jacobian (ic,jc) matrices. 191c 192 integer ia(*), ja(*), ic(*), jc(*), n, jj, k, l, m 193c 194 write(6,1200) 195 k = 0 196 do 33 l = 1,n 197 ia(l) = (l-1)*3+1 198 ic(l) = (l-1)*3+1 199 do 32 m = l,l+2 200 k = k + 1 201 ja(k) = m - 1 202 jc(k) = m - 1 203 32 continue 204 33 continue 205 ia(n+1) = 3*n + 1 206 ic(n+1) = 3*n+1 207 ja(1) = n 208 jc(1) = n 209 ja(k) = 1 210 jc(k) = 1 211c 212 write(6,1300) (ia(jj),jj=1,n+1) 213 write(6,1350) (ja(jj),jj=1,k) 214 write(6,1400) (ic(jj),jj=1,n+1) 215 write(6,1450) (jc(jj),jj=1,k) 216 return 2171200 format('Initial sparse data structures'/) 2181300 format(' ia ',15i4/10(5x,15i4/)) 2191350 format(' ja ',15i4/10(5x,15i4/)) 2201400 format(' ic ',15i4/10(5x,15i4/)) 2211450 format(' jc ',15i4/10(5x,15i4/)) 222 end 223 224 subroutine res (n, t, y, v, r, ires) 225c This subroutine computes the residual vector 226c r = g(t,y) - A(t,y)*v . 227c If ires = -1, only g(t,y) is returned in r, since A(t,y) does 228c not depend on y. 229c No changes need to be made to this routine if n is changed. 230c 231 integer n, ires, i 232 double precision t, y(n), v(n), r(n), r4d, eodsq, one, four, six, 233 1 fact1, fact4 234 common /test1/ r4d, eodsq 235 data one /1.0d0/, four /4.0d0/, six /6.0d0/ 236c 237 call gfun (n, t, y, r) 238 if (ires .eq. -1) return 239c 240 fact1 = one/six 241 fact4 = four/six 242c 243 r(1) = r(1) - (fact4*v(1) + fact1*(v(2) + v(n))) 244 do 10 i = 2,n-1 245 r(i) = r(i) - (fact4*v(i) + fact1*(v(i-1) + v(i+1))) 246 10 continue 247 r(n) = r(n) - (fact4*v(n) + fact1*(v(1) + v(n-1))) 248 return 249c end of subroutine res for the DLSODIS demonstration program 250 end 251 252 subroutine gfun (n, t, y, g) 253c This subroutine computes the right-hand side function g(y,t). 254c It uses r4d = 1/(4*delta) and eodsq = eta/delta**2 255c from the Common block test1. 256c 257 integer n, i 258 double precision t, y(n), g(n), r4d, eodsq, two 259 common /test1/ r4d, eodsq 260 data two/2.0d0/ 261c 262 g(1) = r4d*(y(n)**2 - y(2)**2) + eodsq*(y(2) - two*y(1) + y(n)) 263c 264 do 20 i = 2,n-1 265 g(i) = r4d*(y(i-1)**2 - y(i+1)**2) 266 1 + eodsq*(y(i+1) - two*y(i) + y(i-1)) 267 20 continue 268c 269 g(n) = r4d*(y(n-1)**2 - y(1)**2) + eodsq*(y(1)-two*y(n)+y(n-1)) 270c 271 return 272c end of subroutine gfun for the DLSODIS demonstration program 273 end 274 275 subroutine addasp (n, t, y, j, ip, jp, pa) 276c This subroutine computes the sparse matrix A by columns, adds it to 277c pa, and returns the sum in pa. 278c The matrix A is periodic tridiagonal, of order n, with nonzero elements 279c (reading across) of 1/6, 4/6, 1/6, with 1/6 in the lower left and 280c upper right corners. 281c 282 integer n, j, ip(*), jp(*), jm1, jp1 283 double precision t, y(n), pa(n), fact1, fact4, one, four, six 284 data one/1.0d0/, four/4.0d0/, six/6.0d0/ 285c 286c Compute the elements of A. 287 fact1 = one/six 288 fact4 = four/six 289 jm1 = j - 1 290 jp1 = j + 1 291 if (j .eq. n) jp1 = 1 292 if (j .eq. 1) jm1 = n 293c 294c Add the matrix A to the matrix pa (sparse). 295 pa(j) = pa(j) + fact4 296 pa(jp1) = pa(jp1) + fact1 297 pa(jm1) = pa(jm1) + fact1 298 return 299c end of subroutine addasp for the DLSODIS demonstration program 300 end 301 302 subroutine jacsp (n, t, y, s, j, ip, jp, pdj) 303c This subroutine computes the Jacobian dg/dy = d(g-A*s)/dy by 304c columns in sparse matrix format. Only nonzeros are loaded. 305c It uses r4d = 1/(4*delta) and eodsq = eta/delta**2 from the Common 306c block test1. 307c 308 integer n, j, ip(*), jp(*), jm1, jp1 309 double precision t, y(n), s(n), pdj(n), r4d, eodsq, two, diag, r2d 310 common /test1/ r4d, eodsq 311 data two/2.0d0/ 312c 313 diag = -two*eodsq 314 r2d = two*r4d 315 jm1 = j - 1 316 jp1 = j + 1 317 if (j .eq. 1) jm1 = n 318 if (j .eq. n) jp1 = 1 319c 320 pdj(jm1) = -r2d*y(j) + eodsq 321 pdj(j) = diag 322 pdj(jp1) = r2d*y(j) + eodsq 323 return 324c end of subroutine jacsp for the DLSODIS demonstration program 325 end 326 327 double precision function elkup (n, y, ewt, itol, rtol,atol, lout) 328c This routine looks up approximately correct values of y at t = 0.4, 329c ytrue = y12 or y120 depending on whether n = 12 or 120. 330c These were obtained by running with very tight tolerances. 331c The returned value is 332c elkup = norm of [ (y - ytrue) / (rtol*abs(ytrue) + atol) ]. 333c 334 integer n, itol, lout, i 335 double precision y(n), ewt(n), rtol, atol, y12(12), y120(120), 336 1 y120a(16), y120b(16), y120c(16), y120d(16), y120e(16), 337 2 y120f(16), y120g(16), y120h(8), dvnorm 338 equivalence (y120a(1),y120(1)), (y120b(1),y120(17)), 339 1 (y120c(1),y120(33)), (y120d(1),y120(49)), 340 1 (y120e(1),y120(65)), 341 1 (y120f(1),y120(81)), (y120g(1),y120(97)), 342 1 (y120h(1),y120(113)) 343 data y12/ 344 1 1.60581860d-02, 3.23063251d-02, 1.21903380d-01, 2.70943828d-01, 345 1 4.60951522d-01, 6.57571216d-01, 8.25154453d-01, 9.35644796d-01, 346 1 9.90167557d-01, 9.22421221d-01, 5.85764902d-01, 1.81112615d-01/ 347 data y120a / 348 1 1.89009068d-02, 1.63261891d-02, 1.47080563d-02, 1.39263623d-02, 349 1 1.38901341d-02, 1.45336989d-02, 1.58129308d-02, 1.77017162d-02, 350 1 2.01886844d-02, 2.32742221d-02, 2.69677715d-02, 3.12854037d-02, 351 1 3.62476563d-02, 4.18776225d-02, 4.81992825d-02, 5.52360652d-02/ 352 data y120b / 353 1 6.30096338d-02, 7.15388849d-02, 8.08391507d-02, 9.09215944d-02, 354 1 1.01792784d-01, 1.13454431d-01, 1.25903273d-01, 1.39131085d-01, 355 1 1.53124799d-01, 1.67866712d-01, 1.83334757d-01, 1.99502830d-01, 356 1 2.16341144d-01, 2.33816600d-01, 2.51893167d-01, 2.70532241d-01/ 357 data y120c / 358 1 2.89693007d-01, 3.09332757d-01, 3.29407198d-01, 3.49870723d-01, 359 1 3.70676646d-01, 3.91777421d-01, 4.13124817d-01, 4.34670077d-01, 360 1 4.56364053d-01, 4.78157319d-01, 5.00000270d-01, 5.21843218d-01, 361 1 5.43636473d-01, 5.65330432d-01, 5.86875670d-01, 6.08223037d-01/ 362 data y120d / 363 1 6.29323777d-01, 6.50129662d-01, 6.70593142d-01, 6.90667536d-01, 364 1 7.10307235d-01, 7.29467947d-01, 7.48106966d-01, 7.66183477d-01, 365 1 7.83658878d-01, 8.00497138d-01, 8.16665158d-01, 8.32133153d-01, 366 1 8.46875019d-01, 8.60868691d-01, 8.74096465d-01, 8.86545273d-01/ 367 data y120e / 368 1 8.98206892d-01, 9.09078060d-01, 9.19160487d-01, 9.28460742d-01, 369 1 9.36989986d-01, 9.44763554d-01, 9.51800339d-01, 9.58122004d-01, 370 1 9.63751979d-01, 9.68714242d-01, 9.73031887d-01, 9.76725449d-01, 371 1 9.79811001d-01, 9.82297985d-01, 9.84186787d-01, 9.85466039d-01/ 372 data y120f / 373 1 9.86109629d-01, 9.86073433d-01, 9.85291781d-01, 9.83673704d-01, 374 1 9.81099057d-01, 9.77414704d-01, 9.72431015d-01, 9.65919133d-01, 375 1 9.57609585d-01, 9.47193093d-01, 9.34324619d-01, 9.18631922d-01, 376 1 8.99729965d-01, 8.77242371d-01, 8.50830623d-01, 8.20230644d-01/ 377 data y120g / 378 1 7.85294781d-01, 7.46035145d-01, 7.02662039d-01, 6.55609682d-01, 379 1 6.05541326d-01, 5.53327950d-01, 4.99999118d-01, 4.46670394d-01, 380 1 3.94457322d-01, 3.44389410d-01, 2.97337561d-01, 2.53964948d-01, 381 1 2.14705729d-01, 1.79770169d-01, 1.49170367d-01, 1.22758681d-01/ 382 data y120h / 383 1 1.00271052d-01, 8.13689920d-02, 6.56761515d-02, 5.28075160d-02, 384 1 4.23908624d-02, 3.40811650d-02, 2.75691506d-02, 2.25853507d-02/ 385c 386 if ((n-12)*(n-120) .ne. 0) go to 300 387 if (n .eq. 120) go to 100 388c 389c Compute local error tolerance using correct y (n = 12). 390 call dewset(n, itol, rtol, atol, y12, ewt) 391c 392c Invert ewt and replace y by the error, y - ytrue. 393 do 20 i = 1, 12 394 ewt(i) = 1.0d0/ewt(i) 395 20 y(i) = y(i) - y12(i) 396 go to 200 397c 398c Compute local error tolerance using correct y (n = 120). 399 100 call dewset( n, itol, rtol, atol, y120, ewt ) 400c 401c Invert ewt and replace y by the error, y - ytrue. 402 do 120 i = 1, 120 403 ewt(i) = 1.0d0/ewt(i) 404 120 y(i) = y(i) - y120(i) 405c 406c Find weighted norm of the error. 407 200 elkup = dvnorm (n, y, ewt) 408 return 409c 410c error return 411 300 write(lout,400) n 412 elkup = 1.0d3 413 400 format(/5x,'Illegal use of elkup for n =',i4) 414 return 415c end of function elkup for the DLSODIS demonstration program 416 end 417