opkdemo7.f 16 KB
Newer Older
1
      program opkdemo7
rtoy's avatar
rtoy committed
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
45
      dimension n(1)
rtoy's avatar
rtoy committed
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
79
80
      n(1) = nptsm1
      nm1 = n(1) - 1
rtoy's avatar
rtoy committed
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)
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's avatar
rtoy committed
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)
138
          write (lout,2000) t, rwork(11), iwork(14),(y(i), i=1,n(1))
rtoy's avatar
rtoy committed
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.
153
154
        itemp = n(1)
        errfac = elkup( itemp, y, rwork(21), itol, rtol(j), atol(j) )
rtoy's avatar
rtoy committed
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
166
c      stop
rtoy's avatar
rtoy committed
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

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's avatar
rtoy committed
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
263
264
      dimension y(*), pa(m0,*)
      dimension n(*)
rtoy's avatar
rtoy committed
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).
274
      do 10 i = 1,n(1)
rtoy's avatar
rtoy committed
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
293
294
      dimension y(*), pa(m0,8)
      dimension n(1)
rtoy's avatar
rtoy committed
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
311
312
      pa(n(1),n(1)) = pa(n(1),n(1)) + fact4
      pa(n(1),nm1) = pa(n(1),nm1) + fact1
rtoy's avatar
rtoy committed
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
327
328
      dimension y(*), s(*), pa(m0,*)
      dimension n(*)
rtoy's avatar
rtoy committed
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
365
      pa(mup1,n(1)) = diag
rtoy's avatar
rtoy committed
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
381
382
      dimension y(*), s(8), pa(m0,*)
      dimension n(*)
rtoy's avatar
rtoy committed
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
414
      pa(n(1),nm1) = r2d*y(nm1) + eodsq
rtoy's avatar
rtoy committed
415
416
417
c
c                     n   n
c Compute and store dg /dy
418
      pa(n(1),n(1)) = diag
rtoy's avatar
rtoy committed
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's avatar
rtoy committed
432
      double precision y, ewt, rtol, atol, y9, y99
433
434
c      double precision y99a, y99b, y99c,
c     1   y99d, y99e, y99f, y99g, dvnorm
rtoy's avatar
rtoy committed
435
      dimension y(n), ewt(n), y9(9), y99(99)
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's avatar
rtoy committed
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 /
446
      data y99 /
rtoy's avatar
rtoy committed
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,
450
     1 5.97163151d-02, 6.97479236d-02, 8.08649804d-02, 9.30936515d-02,
rtoy's avatar
rtoy committed
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,
454
     1 3.45833531d-01, 3.70748792d-01, 3.96087655d-01, 4.21766871d-01,
rtoy's avatar
rtoy committed
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,
458
     1 7.47985803d-01, 7.69601151d-01, 7.90342031d-01, 8.10147715d-01,
rtoy's avatar
rtoy committed
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,
462
     1 9.69590361d-01, 9.74555682d-01, 9.78631814d-01, 9.81840924d-01,
rtoy's avatar
rtoy committed
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,
466
     1 8.91098383d-01, 8.61767660d-01, 8.26550038d-01, 7.85131249d-01,
rtoy's avatar
rtoy committed
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,
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's avatar
rtoy committed
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