-
Notifications
You must be signed in to change notification settings - Fork 0
/
acm-toms-690.f
2374 lines (2370 loc) · 92.9 KB
/
acm-toms-690.f
1
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
45
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
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
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
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
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
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
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
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
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
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
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
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
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
C ----- NEW VERSION --------------------
C PDECHEB DISCRETISATION MODULE
C ****************************
C THIS MODULE DISCRETISES MIXED SYSTEMS OF PARTIAL DIFFERENTIAL
C EQUATIONS IN ONE SPACE VARIABLE AND ORDINARY DIFFERENTIAL EQUATIONS.
C THIS IS THE MARK 1 VERSION OF PDECHEB 10TH AUGUST 1987 AS WRITTEN BY
C DR MARTIN BERZINS
C DEPARTMENT OF COMPUTER STUDIES
C THE UNIVERSITY
C LEEDS LS2 9JT
C ALL RIGHTS RETAINED.
C ( DOCUMENTATION MODIFIED ON 25/2/90 )
C THE CLASS OF EQUATIONS THAT CAN BE HANDLED IS GIVEN BY
C
C .
C Q (X,T, U, U , U , U , V , V )
C I - - X -T - XT - -
C
C -M M .
C = X (X R (X,T, U, U , U , U , V, V ))
C I - -X - XT -T - - X
C WHERE
C U = ( U , ... , U ) TRANSPOSE AND I = 1,... , NPDE.
C - 1 NPDE
C
C THE P.D.E. FLUX FUNCTION R(.......) IS ASSUMED TO BE CONTINUOUS
C W.R.T. THE SPACE VARIABLE R BUT THE FUNCTION Q(.... ) IS
C ALLOWED TO BE ONLY PIECEWISE CONTINUOUS PROVIDED THAT THE
C DISCONTINUITIES ARE PLACED AT SPATIAL MESH POINTS.
C THE OTHER VECTORS ARE DEFINED IN THE SAME WAY AS U EXCEPT THAT
C . -
C V = ( V , ... , V ) TRANSPOSE AND V IS SIMILARLY DEFD.
C - 1 NV -
C
C .
C WHERE V AND V ARE THE SOLUTION OF A COUPLED SYSTEM OF
C - - ORDINARY DIFFERENTIAL EQUATIONS OF
C DIMENSION NV.
C IN THE CASE WHEN NV > 0 THIS SYSTEM OF DIFFERENTIAL EQUATIONS IS
C ASSUMED TO HAVE THE FORM.
C .
C FV ( V, V , XI, UI, UI , RI , UI , UI , ) = 0
C -- - - - - - X - - T - XT
C
C .
C WHERE ALL THE VECTORS APART FROM V AND V ARE OF LENGTH NPDE*NXI.
C - -
C AND CONTAIN THE VALUES OF THE P.D.E. VARIABLES U , U , R, U ,U
C X T XT
C AT THE SPATIAL O.D.E. /P.D.E. COUPLING POINTS DEINED BY
C THE VECTOR XI(NXI) .
C
C THE SPATIAL MESH IS BOUNDED BY A AND B A < X < B.
C THE BOUNDARY CONDITIONS HAVE THE FORM
C . .
C B (T) R(X,T, U, U , V, V ) = G (T, U , U , V, V) AT X = A
C A - - -X - - -A - -X - -
C AND
C . .
C B (T) R(X,T, U, U , V, V ) = G (T, U , U , V, V) AT X = B.
C B - - -X - - -B - -X - -
C
C WHERE NOT ALL OF THE FUNCTIONS B AND G ARE SET TO ZERO.
C
C THE INITIAL CONDITIONS ARE GIVEN BY
C U (X, 0) = K (X) AND V(O) = K
C - -1 - -2 0
C THE DISCRETISATION METHOD USED BY THIS MODULE IS BASED ON A C -
C COLLOCATION METHOD AND EVALUATES THE P.D.E. FUNCTIONS INBETWEEN
C THE USER SUPPLIED MESH POINTS . ANY DISCONTINUITIES IN THE P.D.E.
C DEFINING FUNCTION Q MUST THEREFORE BE AT THE USER SUPPLIED
C BREAK POINTS .
C REFERENCES
C ----------
C BERZINS M. AND DEW P.M.
C CHEBYSHEV POLYNOMIAL SOFTWARE FOR ELLIPTIC PARABOLIC P.D.ES
C ACM TRANS. ON MATH. SOFT. 1990 PP XX - YY.
C
C BERZINS M. AND DEW P.M.
C A NOTE ON C0 CHEBYSHEV METHODS FOR PARABOLIC EQUATIONS.
C IMA JOURNAL OF NUMERICAL ANALYSIS (1987) 7, 15-37.
C
C BERZINS M. AND DEW P.M. DEPT OF COMPUTER STUDIES
C THE UNIVERSITY , LEEDS , LS2 9JT ,REPORT NO 180 .
C
C----------------------------------------------------------------------
C
C HOW TO USE THIS MODULE
C **********************
C (1) DECIDE ON THE FORM OF THE SPATIAL DISCRETISATION METHOD TO BE
C USED. THIS MODULE ALLOWS YOU TO DEFINE A SET OF (NEL + 1) SPATIAL
C BREAKPOINTS , XBK(I) I = 1,NEL+1 . THESE BREAKPOINTS IN TURN
C DEFINE NEL SPATIAL ELEMENTS. I.E.
C XBK(I) =< ELEMENT I =< XBK(I+1).
C IN THE CASE WHEN THE FUNCTION Q(.....) IN THE P.D.E. DEFINITION
C HAS DISCONTINIUITIES YOU MUST PLACE A BREAKPOINT AT EACH
C DISCONTINUITY.
C PDECHEB WILL APPROXIMATE THE SOLUTION TO THE P.D.E. IN THE SPACE
C DIMENSION BY USING A PIECEWISE CHEBYSHEV POLYNOMIAL BETWEEN EACH
C PAIR OF BREAKPOINTS. THE DEGREE OF THIS POLYNOMIAL IS SPECIFIED
C BY THE VARIABLE NPOLY WHICH MUST BE GREATER THAN OR
C EQUAL TO 1. WHEN IT IS 1 THE SPATIAL MESH CONSISTS ONLY OF THE
C BREAKPOINTS AND A LINEAR POLYNOMIAL IS USED TO APPROXIMATE THE
C SOLUTION BETWEEN THESE POINTS.
C THE ONLY PRE-SET UPPER LIMIT TO THE DEGREE OF POLYNOMIAL THAT
C CAN BE USED IS A SOME WHAT ARBITRARY LIMIT OF 50.
C
C
C (2) SET NPDE = NUMBER OF P.D.E.S TO BE SOLVED , MUST BE >= 0 .
C SET NEL = NUMBER OF SPATIAL ELEMENTS TO BE USED , MUST BE >= 0.
C DEFINE AN ARRAY OF BREAKPOINTS IN THE DOUBLE PRECISION ARRAY
C XBK(IBK) WHERE IBK = NEL + 1 AND
C XBK(I) < XBK(I+1) , I = 1,NEL
C SET NPOLY TO THE DEGREE OF THE POLYNOMIALIN EACH ELEMENT , > 1 .
C NOTES ON THE CHOICE OF NPOLY AND XBK(IBK)
C -----------------------------------------
C IT SHOULD BE NOTED THAT THE PDECHEB SOFTWARE HAS NO MEANS OF
C ESTIMATING OR CONTROLLING THE SPATIAL DISCRETISATION ERROR.
C THE ERROR INCURRED WILL DEPEND ON THE NUMBER AND POSITION
C OF THE BREAK POINTS AND ON THE DEGREE OF POLYNOMIAL USED.
C THE GENERAL ADVICE IS USE AS FEW BREAK POINTS AS POSSIBLE
C AND AS HIGH A DEGREE OF POLYNOMIAL AS SEEMS SENSIBLE FOR
C THE PROBLEM AT HAND. THE APPROPRIATENESS OF A GIVEN DEGREE
C OF POLYNOMIAL CAN BE JUDGED BY THE FACT THAT THE HIGHER
C DEGREE COEFFICIENTS OF THE POLYNOMIAL EXPANSION SHOULD BE
C SMALL IN COMPARISON WITH THE LOWER POLYNOMIAL DEGREE.
C IT SHOULD ALSO BE NOTED THAT THE LOCAL TIME ERROR TOLERANCES
C ( THE PARAMETERS RTOL AND ATOL - SEE SECTION 5 BELOW)
C PASSED TO THE DASSL CODE SHOULD BE AN ORDER OF MAGNITUDE
C SMALLER THAN THE EXPECTED SPATIAL ERROR. INEVITABLY THESE
C AREAS OF UNCERTAINTY MEAN THAT SOME EXPERIMENTATION WITH
C A GIVEN PROBLEM IS NECESSARY BEFORE AN CONFIDENCE CAN BE
C PLACED IN THE NUMERICAL RESULTS.
C
C
C SET M FOR SPACE CO-ORDINATE TYPE
C = 0 FOR CARTESIAN, = 1 FOR CYLINDRICAL, = 2 FOR SPHERICAL
C SET NV = NUMBER OF O.D.E.S COUPLED TO THE P.D.E.S
C SET NXI = NUMBER OF SPACE POINTS WHERE THE O.D.E.S ARE COUPLED
C SET XI(NXI) TO THE VALUES OF THE COUPLING POINTS
C FOR USE BY THE ROUTINE PDECHB WHICH DEFINES THE O.D.E. SYSTEM
C BEING SOLVED BY THE INTEGRATOR.
C SET NPTS = NEL * NPTL-1 , THIS IS THE TOTAL NUMBER OF SPATIAL
C DISCRETISATION POINTS USED BY THIS MODULE.
C DECLARE A DOUBLE PRECISION ARRAY X, OF DIMENSION NPTS ; X(NPTS)
C AND PUT X(1) = XBK(1) = A AND X(NPTS) = XBK(NEL+1) = B
C WHERE A AND B ARE THE LEFT AND RIGHT EDGES OF THE SPTAIAL MESH.
C THIS ARRAY WILL BE USED TO RETURN TO YOU THE SPATIAL MESH POINTS
C USED BY THIS SPATIAL DISCRETISATION MODULE.
C DECLARE A DOUBLE PRECISION ARRAY OF LENGTH NEQN WHERE
C NEQN = NPDE * NPTS + NV
C THAT IS USED TO HOLD THE SOLUTION VECTOR . NEQN IS THE NUMBER OF
C ORDINARY DIFFERENTIAL EQUATIONS THAT MUST BE PASSED ACROSS TO
C THE DASSL PACKAGE . THE SOLUTION TO THIS SYSTEM OF ORDINARY
C DIFFFERENTIAL EQUATIONS THAT IS GENERATED BY DISCRETISING THE
C CLASS OF P.D.E.S DEFINED ABOVE IS ORDERED IN ,SAY, U(NEQN) AS
C FOLLOWS.
C U(I) , I = (J-1) * NPDE + K , K = 1,...,NPDE , J = 1 ,.. ,NPTS
C CONTAINS THE SOLUTION FOR P.D.E. K AT MESH POINT X(J).
C U(L) , L = NPDE*NPTS + L1 , L1 = 1,..., NV
C CONTAINS THE COUPLED O.D.E. COMPONENT V(L1)
C DEFINE A DOUBLE PRECISION WORKSPACE OF LENGTH IWK WHERE
C IWK = NPTL*(3*NPTL + 2 + 7*NPDE + NEL) +NXI*(5*NPDE+1) + NV + 2
C THIS IS THE WORKSPACE THAT MUST BE PASSED ACROSS TO THE DASSL
C ROUTINE AS THE WORKSPACE FOR THE O.D.E. RESIDUAL DEFINING ROUTINE
C PDECHB.
C
C
C
C (3) PROVIDE A SET OF ROUTINES WHICH DESCRIBE THE PRECISE FORM OF THE
C P.D.E. TO BE SOLVED. FOUR ROUTINES MUST BE PROVIDED AND THE NAMES
C OF THESE ROUTINES ARE FIXED. THESE ROUTINES ARE:
C
C SPDEFN : FORMS THE FUNCTIONS Q AND R IN THE P.D.E. DESCRIBED
C ABOVE. THIS ROUTINE FORMS THE VALUES OF THE FUNCTIONS
C Q AND R OVER SEVERAL MESHPOINTS SIMULTANEOUSLY.
C IN FACT AT THE X(NPTL) POINTS IN ONE ELEMENT AT A TIME.
C THE MESH POINTS USED ARE INTERNALLY
C GENERATED BY THE DISCRETISATION ROUTINE AND ARE
C BETWEEN THE USER DEFINED BREAKPOINTS.
C SBNDR : FORMS THE FUNCTIONS B AND G ASSOCIATED WITH THE
C BOUNDARY CONDITIONS FOR THE P.D.E. ABOVE.
C UVINIT : SUPPLIES THE INITIAL VALUES OF THE P.D.E. PART AND ALSO
C SUPPLIES THE INITIAL VALUES OF THE O.D.E. PART.
C SODEFN : SUPPLIES THE ODE RESIDUAL AS DEFINED BY THE FUNCTION
C FV ABOVE.
C NOTE: THE P.D.E. SOLUTION VALUES AT THE COUPLING POINTS
C PASSED INTO SODEFN ARE DEFINED BY POLYNOMIAL INTERPOLATION
C ON THE VALUES AT THE P.D.E. SPATIAL MESH POINTS.
C
C N.B. EXAMPLES OF THESE ROUTINES FOR THREE PROBLEMS ARE PROVIDED
C IN THE EXAMPLE PROBLEMS SECTION BELOW.
C
C
C (4) CALL THE INITIALISATION ROUTINE INICHB, USING THE FORM
C
C CALL INICHB(NEQN,NPDE,NPTS,X,U,WK,IWK,M,TS,IBAND,ITIME,XBK,
C * IBK,NEL,NPOLY,NV,NXI,XI,IDEV)
C***********************************************************************
C ROUTINE FOR INITIALISATION OF CHEBYSHEV GENERALIZED COLLOCATION METHOD
C
C PARAMETER LIST
C ----------------
C NEQN: EMPTY ON ENTRY, ON EXIT IT CONTAINS THE NUMBER OF
C O.D.E.S GENERATED BY THE DISCRETISED FORM OF THE
C P.D.E. , GIVEN BY NPDE*NEL*(NPTL-1) + NPDE + NV.
C
C NPDE NUMBER OF PARABOLIC P.D.E.S IN ONE SPACE DIMENSION
C
C NPTS NUMBER OF SPATIAL GRID POINTS USED IN M.O.L. SOLUTION.
C NOTE THIS SHOULD BE EQUAL TO (NPTL-1)*NEL + 1
C
C X(NPTS) EMPTY ARRAY ON ENTRY . ON EXIT THIS ARRAY
C CONTAINS THE MESH USED IN SEMI-DISCRETISATION
C
C M =0,1,2 IF CARTESIAN CYLINDRICAL OR SPHERICAL POLARS.
C
C U(NEQN) SOLUTION VECTOR EMPTY ON ENTRY CONTAINS INITIAL
C VALUES ON EXIT. THIS ARRAY IS ORDERED AS FOLLOWS.
C U(1) - U(NPDE*NPTS) P.D.E. SOLUTION COMPONENTS.
C U(NPDE*NPTS+1) - U(NEQN) O.D.E. COMPONENTS THAT ARE
C COUPLED TO THE P.D.E.
C
C WK(IWK) REAL WORKSPACE USED TO PASS FOUR MATRICES AND VARIOUS
C USEFUL VECTORS TO THE O.D.E.FUNCTION CALL ROUTINE
C RESID SEE BELOW FOR A DETAILED DESCRIPTION.
C
C TS STARTING LEVEL OF TIME INTEGRATION.
C
C XBK(IBK) REAL ARRAY OF BREAK POINTS IBK = NEL +1 WHERE
C NEL IS THE NUMBER OF SPATIAL ELEMENTS.
C XBK(1) = XLEFT
C XBK(I) =< XBK(I+1) I = 1,...,NEL .
C XBK(IBK) = XRIGHT.
C NEL THE NUMBER OF SPATIAL ELEMENTS , >= 1
C
C NPOLY THE DEGREE OF THE APPROXIMATING POLYNOMIAL USED
C BETWEEN EACH PAIR OF BREAKPOINTS .
C ITIME THIS MUST BE SET = 1 ON THE CALL OF THIS MODULE
C PRIOR TO THE DASSL PACKAGE BEING CALLED.
C ONCE DASSL HAS RETURNED THIS ROUTINE MAY BE CALLED
C WITH ITIME = 2 TO RECOVER THE SPATIAL MESH USED
C (THIS IS PLACED IN THE ARRAY X(NPTS) ).
C
C NV THE NUMBER OF AUXILARY O.D.E.S THAT ARE COUPLED TO
C THE P.D.E.
C NXI THE NUMBER OF COUPLING POINTS AT WHICH P.D.E. VALUES
C ARE USED TO DEFINE THE O.D.E.S
C
C XI(NXI) A VECTOR SPECIFYING THE POSITION OF THESE POINTS.
C NOTE THAT THESE POINTS MUST BE DISTINCT FROM THE
C BREAK - POINTS .
C
C IDEV NUMBER OF OUTPUT CHANNEL ON WHICH ERROR MESSAGES TO
C DO WITH THE COLLOCATION DISCRETISATION WILL APPEAR.
C
C
C
C (5) SET TS AND TOUT FOR START AND END INTEGRATION TIMES
C SET INFO, WORK ARRAYS AS REQUIRED FOR TIME INTEGRATION
C AND CALL THE DASSLD ROUTINE AS FOLLOWS
C CALL DDASSL (PDECHB, NEQ, T, Y, YDOT, TOUT, INFO, RTOL, ATOL,
C * , IDID, RWORK, LRW, IWORK, LIW, WKRES, IRESWK, DGEJAC)
C
C DASSL CODE SOLVES A SYSTEM OF DIFFERENTIAL/ALGEBRAIC EQUATIONS OF
C THE FORM G(T,Y,YDOT) = 0.
C
C VALUES FOR Y AND YPRIME AT THE INITIAL TIME MUST BE GIVEN AS INPUT
C THESE VALUES MUST BE CONSISTENT, (THAT IS. IF T,Y,YDOT ARE THE GIVEN
C INITIAL VALUES, THEY MUST SATISFY G(T,Y,YDOT) = 0.)
C THE SUBROUTINE SOLVES THE SYSTEM FROM T TO TOUT. IT IS EASY TO
C CONTINUE THE SOLUTION TO GET RESULTS AT ADDITIONAL TOUT. THIS IS THE
C INTERVAL MODE OF OPERATION. INTERMEDIATE RESULTS CAN ALSO BE
C OBTAINED EASILY BY USING THE INTERMEDIATE-OUTPUT CAPABILITY.
C
C ------------AN OVERVIEW OF ARGUMENTS TO DASSL-----------------------
C THE PARAMETERS ARE
C
C PDECHB -- THIS IS A SUBROUTINE PROVIDED BY PDECHEB TO DEFINE THE
C DIFFERENTIAL/ALGEBRAIC SYSTEM
C
C NEQ -- THE NUMBER OF DIFFERENTIAL/ALGEBRAIC EQUATIONS TO BE SOLVED
C
C T -- THIS IS THE CURRENT VALUE OF THE INDEPENDENT VARIABLE.
C
C TOUT -- THIS IS A POINT AT WHICH A SOLUTION IS DESIRED.
C
C INFO(*) -- THE BASIC TASK OF THE CODE IS TO SOLVE THE System
C FROM T TO TOUT AND RETURN AN ANSWER AT TOUT. INFO(*) IS
C INTEGER ARRAY WHICH IS USED TO COMMUNICATE EXACTLY HOW
C YOU WANT THIS TASK TO BE CARRIED OUT.
C
C Y(*) -- THIS ARRAY CONTAINS THE SOLUTION COMPONENTS AT T
C
C YDOT(*) -- THIS ARRAY CONTAINS THE DERIVATIVES OF Y(*) AT T
C
C RTOL,ATOL -- THESE QUANTITIES REPRESENT ABSOLUTE AND Relative error
C TOLERANCES WHICH YOU PROVIDE TO INDICATE HOW ACCURATELY
C YOU WISH THE SOLUTION TO BE COMPUTED. YOu may choose
C THEM TO BE BOTH SCALARS OR ELSE BOTH VECtors.
C
C IDID -- THIS SCALAR QUANTITY IS AN INDICATOR REPORTING WHAT THE CODE
C YOU MUST MONITOR THIS INTEGER VARIABLE TO DECIDE WHAT ACTION
C TO TAKE NEXT.
C
C RWORK(*),LRW -- RWORK(*) IS A REAL WORK ARRAY OF Length lrw which
C PROVIDES THE CODE WITH NEEDED STORAGE SPACE.
C
C IWORK(*),LIW -- IWORK(*) IS AN INTEGER WORK ARRAY OF LENGTH LIW
C WHICH PROVIDES THE CODE WITH NEEDED storage space.
C
C WKRES,IRESWK -- THESE ARE REAL AND INTEGER PARAMETER ARRAYS WHICH
C ARE USED TO COMMUNICATE INFORMATION FROM THE
C INITIALISATION ROUTINE INICHB TO THE SPATIAL
C DISCRETISATION SUBROUTINE PDECHB.
C
C DGEJAC -- THIS IS THE NAME OF A DUMMY SUBROUTINE WHICH IS PROVIDED
C BY THE PDECHEB SOFTWARE. IT MUST BE DECLARED AS EXTERNAL
C IN THE CALLING PROGRAM.
C
C QUANTITIES WHICH ARE USED AS INPUT ITEMS ARE
C NEQ, T Y(*), YDOT(*), TOUT, INFO(*), RTOL, ATOL , RWORK(1),
C RWORK(2), RWORK(3), LRW, IWORK(1), IWORK(2),IWORK(3),AND LIW.
C
C QUANTITIES WHICH MAY BE ALTERED BY THE CODE ARE
C T, Y(*), YDOT(*), INFO(1), RTOL, ATOL, IDID, RWORK(*) AND IWORK(*)
C
C (6) POST PROCESS THE SOLUTION.
C
C THE SOLUTION VECTOR RETURNED BY DASSL CAN BE USED FOR POST-
C PROCESSING IN A NUMBER OF WAYS.
C
C SPATIAL INTERPOLATION
C ---------------------
C THE VECTOR Y(T) RETURNED BY DASSL CONSISTS ONLY OF SOLUTION
C VALUES AT THE MESH POINTS DEFINED BY INICHB. THE FOLLOWING
C INTERPOLATION ROUTINE ENABLES SOLUTION VALUES AT OTHER POINTS
C TO BE OBTAINED.
C
C SUBROUTINE INTERC(XP,UP,NP,U,NEQ,NPDE,IFLAG,ITYPE,WK,IWK)
C********************************************************************
C
C SPACE INTERPOLATION ROUTINE FOR POST-PROCESSING OF SOLUTION
C PRODUCED BY DASSL.
C THIS ROUTINES PROVIDES VALUES OF THE SOLUTION AND POSSIBLY THE
C FIRST DERIV IN SPACE AND THE FLUX ON THE MESH XP(NP).
C
C PARAMETERS
C --------------
C NPDE ON ENTRY MUST CONTAIN NO OF PARABOLIC EQUATIONS
C NPTS ON ENTRY MUST CONTAIN THE NUMBER OF SPATIAL
C MESH POINTS USED IN TIME INTEGRATION.
C NP ON ENTRY MUST CONTAIN THE NUMBER OF SPATIAL
C INTERPOLATION POINTS
C XP(NP) ARRAY WHICH ON ENTRY
C CONTAINS THE SPATIAL INTERPOLATION POINTS
C WE ASSUME THAT
C XP(I) < XP(I+1) , I = 1,...,NP-1
C UP(NPDE,NP,ITYPE) EMPTY ARRAY FOR THE INTERPOLATED VALUES AT
C THE CURRENT TIME LEVEL. THE VALUES OF THIS
C ARRAY ON EXIT DEPEND ON THE PARAMETER ITYPE.
C U(NPDE,NPTS) THE CURRENT SOLUTION VECTOR COMPUTED BY THE ODE
C TIME INTEGRATOR MUST BE SUPPLIED IN THIS VECTOR.
C IFLAG ERROR FLAG = 0 ON SUCCESSFUL RETURN
C = 1 IF EXTRAPOLATION TRIED.
C = 2 IF WORKSPACE NOT INITIAL
C ISED ON ENTRY BY INICHB.
C = 3 ILLEGAL VALUE OF ITYPE.
C ITYPE = 1 ONLY THE SOLUTION IS OUTPUT IN THE ARRAY UP
C UP(J,K,1) HOLDS U(XP(K),T) FOR PDE J
C 2 AS FOR 1 BUT THE FIRST DERIV IS ALSO OUTPUT.
C UP(J,K,2) HOLDS D/DX U(XP(K),T).
C
C WK(IWK) THE WORKSPACE USED BY THE CHEBYSHEV METHOD. THIS
C MUST BE THE WORKSPACE INITIALISED BY INICHB.
C
C
C !*********************************************************!
C ! IN THE CASE WHEN THE EXACT SOLUTION IS NOT KNOWN IT MAY !
C ! STILL BE NECESSARY TO SUPPLY A DUMMY ROUTINE EXACT TO !
C ! SATISFY LOADER REQUIREMENTS (SEE THE NEXT SECTION FOR A !
C ! DESCRIPTION OF EXACT. !
C !*********************************************************!
C
C ESTIMATING THE ERROR WHEN THE EXACT SOLUTION IS KNOWN.
C ------------------------------------------------------
C THIS CAN BE DONE BY THE FOLLOWING CALL
C
C CALL ERROR(U,NPDE,NPTS,X,M,ENORM,GERR,T,RELERR,ABSERR,
C * ITRACE,RWK,IWK)
C
C**********************************************************************
C THE FOLLOWING ROUTINE COMPUTES THE ERROR ENORM IN THE NUMERICAL
C SOLUTION BY USING A COMBINATION OF THE L2 FUNCTION AND VECTOR
C NORMS. GERR IS THE MAXIMUM ERROR AT THE GRID POINTS
C THE EXACT SOLUTION IS ASSUMED TO BE GIVEN BY THE USER PROVIDED
C SUBROUTINE EXACT(T,NPDE, NP, XP, US)
C DOUBLE PRECISION US(NPDE, NP),XP(NP),T
C WHERE US(J,I) ON EXIT CONTAINS THE SOLUTION AT TIME T
C FOR NPDE J AT THE MESH POINT XP(I)
C
C PARAMETER LIST
C --------------
C U(NEQN) SOLUTION VECTOR COMPUTED BY DASSL AT TIME T . ON
C ENTRY THIS ARRAY IS ASSUMED TO BE ORDERED AS FOLLOWS
C U(1) - U(NPDE*NPTS) P.D.E. SOLUTION COMPONENTS.
C U(NPDE*NPTS+1) - U(NEQN) O.D.E. COMPONENTS THAT ARE
C COUPLED TO THE P.D.E.
C
C NPDE NUMBER OF PARABOLIC P.D.E.S IN ONE SPACE DIMENSION
C
C NPTS NUMBER OF SPATIAL GRID POINTS USED IN M.O.L. SOLUTION.
C NOTE THIS SHOULD BE EQUAL TO (NPTL-1)*NEL + 1
C
C X(NPTS) ON ENTRY THIS ARRAY MUST
C CONTAIN THE MESH USED IN SEMI-DISCRETISATION
C
C M =0,1,2 IF CARTESIAN CYLINDRICAL OR SPHERICAL POLARS.
C
C ENORM L2 ERROR NORM ESTIMATED BY USING TRAPEZOIDAL RULE
C WITH 100 EVENLY SPACED POINTS IS OUTPUT IN ENORM
C
C GERR MAXIMUM GRID ERROR OVER THE ARRAY OF SPATIAL GRID
C POINTS X(NPTS) IS OUTPUT IN GERR
C
C T CURRENT TIME LEVEL OF TIME INTEGRATION ( INPUT).
C
C RELERR RELATIVE ERROR TOLERANCE SUPPLIED TO DASSL (RTOL IN
C THE CALL TO THAT ROUTINE) (INPUT)
C
C ABSERR ABSOLUTE ERROR TOLERANCE SUPPLIED TO DASSL (ATOL IN
C THE CALL TO THAT ROUTINE). (INPUT)
C
C ITRACE INTEGER TRACE LEVEL SET TO ZERO FOR NO TRACE SET =1
C FOR TRACE INFORMATION. (INPUT)
C
C RWK(IWK) REAL WORKSPACE INITIALISED BT INICHB AND PASSED TO
C THE D.A.E.FUNCTION CALL ROUTINE RESID
C SEE BELOW FOR A DETAILED DESCRIPTION.(INPUT)
C**********************************************************************
C EXAMPLE PROBLEM ONE
C SOLUTION OF MOVING BOUNDARY PROBLEM BY CO-ORDINATE TRANSFORMATION.
C********************************************************************
C THIS PROBLEM IS THE ONE PHASE STEFAN PROBLEM (HOFFMAN (1977) ) SEE
C FURZELAND R.M. A COMPARATIVE STUDY OF NUMERICAL METHODS FOR MOVING
C BOUNDARY PROBLEMS. J.I.M.A. (1977) ,26, PP 411 - 429.
C THE PROBLEM HAS MELTING DUE TO HEAT INPUT AT THE FIXED
C BOUNDARY . THE P.D.E. IS DEFINED BY THE EQUATIONS
C U = U 0 < Y < S(T) , 0.1 < T < 1
C T YY
C U = - EXP(T) , Y = 0
C Y .
C U = 0 AND S(T) = - U ON THE MOVING BOUNDARY Y = S(T).
C Y
C AND THE INITIAL SOLUTION VALUES AT T = 0.1 ARE GIVEN BY THE ANALYTIC
C SOLUTION
C U = EXP(T-Y) - 1 , S(T) = T.
C THE PROBLEM IS REWRITTEN BY USING THE CO-ORDINATE TRANSFORMATION
C GIVEN BY X(T) = Y / S(T) . THE EQUATIONS THEN READ
C .
C S * S * U - S * S * X * U = U , X IN (0,1).
C T X XX
C WITH THE NEUMANN TYPE BOUNDARY CONDITIONS
C .
C U = - EXP(T) AT X=0 AND U = - S(T) * S(T) AT X = 1
C X X
C AND THE O.D.E. COUPLING POINT EQUATION AT X = 1 WHICH IMPLICITLY
C DEFINES S(T) IS GIVEN BY
C U(1,T) = 0
C THE EXACT SOLUTION IS NOW DEFINED BY
C U(X,T) = EXP((T - X*S(T)) , S(T) = T
C
C WE SHALL NOW DETAIL THE ROUTINES NEEDED TO DESCRIBE THIS PROBLEM.
C PROBLEM DESCRIPTION ROUTINES
C ******************************
C EXACT SOLUTION
C SUBROUTINE EXACT( TIME, NPDE, NPTS, X, U)
C ROUTINE FOR P.D.E. EXACT VALUES (IF KNOWN)
C INTEGER NPDE, NPTS
C DOUBLE PRECISION X(NPTS), U(NPDE,NPTS), TIME, P
C P=DSQRT(2.0D0)*0.5D0
C DO 10 I = 1,NPTS
C10 U(1,I) = DEXP( TIME * (1 - X(I))) - 1.0D0
C RETURN
C END
C SUBROUTINE UVINIT( NPDE, NPTS, X, U, NV, V)
C ROUTINE FOR P.D.E. INITIAL VALUES.
C INTEGER NPDE, NPTS, NV
C DOUBLE PRECISION X(NPTS), U(NPDE,NPTS), TIME, V(NV)
C TIME=0.1D0
C V(1)=0.1D0
C CALL EXACT(TIME,NPDE,NPTS,X,U)
C RETURN
C END
C
C SUBROUTINE SPDEFN(T, X, NPTL, NPDE, U, DUDX, UDOT, UTDX, Q, R,
C 1 NV, V, VDOT, IRES)
C PROBLEM INTERFACE FOR THE MOVING BOUNDARY PROBLEM.
C INTEGER NPTL, NPDE, NV, I, IRES
C DOUBLE PRECISION X(NPTL), U(NPDE,NPTL), DUDX(NPDE,NPTL), T,
C 1 V(1), VDOT(1), F(NPDE,NPTL), Q(NPDE,NPTL) ,R(NPDE,NPTL),
C 2 UDOT(NPDE,NPTL), UTDX(NPDE,NPTL)
C DO 10 I = 1,NPTL
C R(1,I) = DUDX(1,I)
C Q(1,I) = V(1)*V(1)*UDOT(1,I) -X(I)*VD(I)*DUDX(1,I) * V(1)
C10 CONTINUE
C RETURN
C END
C SUBROUTINE SBNDR( T, BETA, GAMMA, U, UX, UDOT, UTDX, NPDE,
C 1 LEFT, NV, V, VDOT, IRES)
C THIS SUBROUTINE PROVIDES THE LEFT AND RIGHT BOUNDARY VALUES
C FOR THE MOVING BOUNDARY PROBLEM IN THE FORM.
C BETA(I) * DU/DX(I) = GAMMA(I)
C WHERE I = 1,NPDE AND GAMMA IS A FUNCTION OF U,X AND T
C
C INTEGER NPDE, NV, IRES
C LOGICAL LEFT
C DOUBLE PRECISION BETA(NPDE), GAMMA(NPDE), U(NPDE), UX(NPDE)
C - ,T, V(1), VDOT(1), UDOT(NPDE), UTDX(NPDE)
C BETA(1) = 1.0D0
C IF(LEFT)THEN
C GAMMA(1) = -V(1)*DEXP(T)
C ELSE
C GAMMA(1) = -V(1)*VD(1)
C END IF
C RETURN
C END
C
C SUBROUTINE SODEFN(T, NV, V, VDOT, NPDE, NXI, XI, UI, UXI, RI,
C 1 UTI, UTXI, VRES, IRES)
C ROUTINE TO PROVIDE RESIDUAL OF COUPLED ODE SYSTEM FOR THE
C MOVING BOUNDARY PROBLEM.
C NOTE HOW IRES CAN BE RESET TO COPE WIH ILLEGAL VALUES OF THE
C MOVING BOUNDARY POSITION V(1).
C INTEGER NPDE, NXI, NV, IRES
C DOUBLE PRECISION T, XI(NXI), UI(NPDE,NXI), UXI(NPDE,NXI),
C 1 RI(NPDE,NXI), UTI(NPDE,NXI), UTXI(NPDE,NXI), VRES(NV),
C 2 V(NV), VDOT(NV), TEM
C VRES(1) = UI(1,1)
C TEM = 1.0D0
C IF(IRES .EQ. -1)TEM = 0.0D0
C IF(V(1) .LT. 0.0D0)IRES = -1
C VRES(1) = TEM * UI(1,1)
C RETURN
C END
C
C EXAMPLE PROGRAM ONE ....................
C
C C0 COLLOCATION PARAMETERS
C PARAMETER ( IBK = 21, NEL = IBK-1 , NPDE = 1, NV = 1,
C 1 NPOLY = 2, NPTS = NEL*NPOLY+1, NXI = 1,
C 2 NEQ = NPTS * NPDE + NV,
C 3 NWKRES= (NPOLY+1) * (5*NXI + 3*NPOLY+NEL+5+7*NPDE) +
C 4 NPDE * 8 + 3 + NV + NXI,
C DDASSL TIME INTEGRATION PARAMETERS
C 5 MAXORD = 5, LRW = 40 + (MAXORD+4) * NEQ + NEQ**2,
C 6 LIW = 20 + NEQ )
C
C INTEGER IWORK(LIW), INFO(15), IBAND, M, ITIME, I, IDID, IRESWK,
C 1 IDEV, ITRACE
C DOUBLE PRECISION XBK(IBK), X(NPTS), Y(NEQ), YDOT(NEQ),
C 1 WKRES(NWKRES), RWORK(LRW), XI(1), T, TOUT, RTOL, ATOL,
C 2 ENORM, GERR, VERROR, CTIME, TOL
C EXTERNAL PDECHB, DGEJAC
C COMMON /SDEV2/ ITRACE, IDEV
C COMMON /PROB1/ TOL
C TOL = 0.1D-5/50.D0
C M = 0
C T = TOL
C IDEV = 4
C ITRACE = 1
C WRITE(IDEV,9)NPOLY, NEL
C9 FORMAT(' TEST PROBLEM 1'/' ***********'/' POLY OF DEGREE =',I4,
C 1 ' NO OF ELEMENTS = ',I4)
C XI(1) = 1.0D0
C DO 10 I = 1,IBK
C10 XBK(I) = (I-1.0D0)/(IBK-1.0D0)
C INITIALISE THE P.D.E. WORKSPACE
C ITIME = 1
C CALL INICHB(NEQ, NPDE, NPTS, X, Y, WKRES, NWKRES, M, T, IBAND,
C 1 ITIME, XBK, IBK, NEL, NPOLY, NV, NXI, XI, IDEV)
C IF(ITIME .EQ. -1)THEN
C WRITE(IDEV, 15)
C15 FORMAT(' INICHB ROUTINE RETURNED ITIME = -1 - RUN HALTED ')
C GOTO 100
C END IF
C SETUP DASSL PARAMETERS
C RTOL = TOL
C ATOL = TOL
C DO 20 I = 1,11
C20 INFO(I) = 0
C BANDED MATRIX OPTION WHEN INFO(6) = 1
C IF(INFO(6) .EQ. 1)THEN
C IWORK(1) = IBAND
C IWORK(2) = IBAND
C END IF
C30 TOUT = T * 10.0D0
C IF(TOUT .GE. 2.D0)TOUT =2.0D0
C CALL DDASSL( PDECHB, NEQ, T, Y, YDOT, TOUT, INFO, RTOL, ATOL,
C 1 IDID, RWORK, LRW, IWORK, LIW, WKRES, IRESWK, DGEJAC)
C IF( IDID .LT. 0 )THEN
C DASSL FAILED TO FINISH INTEGRATION.
C WRITE(IDEV,40)T,IDID
C40 FORMAT(' AT TIME T = ',D11.3,' DASSL RETURNED IDID =',I3)
C GOTO 100
C ELSE
C DASSL INTEGRATED TO T = TOUT
C CALL TO POST PROCESSING HERE E.G. SPACE INTERPOLATION.
C ITRACE = 1
C CALL ERROR( Y, NPDE, NPTS, X, M, ENORM, GERR, T, RTOL, ATOL,
C 1 ITRACE, WKRES, NWKRES)
C ITRACE = 0
C VERROR = Y(NEQ) - T
C WRITE(IDEV,50)Y(NEQ),VERROR
C50 FORMAT(' MOVING BOUNDARY IS AT ',D12.4,' WITH ERROR=',D12.4)
C IF(TOUT .LT. 1.99D0 ) GOTO 30
C END IF
C100 CONTINUE
C WRITE(IDEV,110)IWORK(11),IWORK(12),IWORK(13)
C110 FORMAT(' NSTEPS =',I5,' NRESID =',I5,' JAC = ',I4)
C STOP
C END
C
C
C
C EXAMPLE PROBLEM TWO
C ********************
C THIS PROBLEM IS DEFINED BY
C -2 2 2
C U U = X ( X U U ) + 5 U + 4 X U U , X IN (0,1)
C T X X X
C
C THE LEFT BOUNDARY CONDITION AT X = 0 (LEFT = .TRUE. ) IS GIVEN BY
C U (0,T) = 0.0
C X
C THE RIGHT BOUNDARY CONDITION IS (LEFT = .FALSE.)
C U( 1,T) = EXP ( -T )
C
C THE INITIAL CONDITION IS GIVEN BY THE EXACT SOLUTION ;
C U( X, T ) = EXP ( 1 - X*X - T ) , X IN ( 0,1)
C 2
C**********************************************************************
C
C C0 COLLOCATION PARAMETERS
C PARAMETER ( IBK = 2, NEL = IBK-1 , NPDE = 1, NV = 0,
C 1 NPOLY = 10, NPTS = NEL*NPOLY+1, NXI = 0,
C 2 NEQ = NPTS * NPDE + NV,
C C NWKRES= 2*(NPOLY+1)*(NPOLY+NEL+2) + 2 + NV +
C 3 NWKRES= (NPOLY+1) * (5*NXI + 3*NPOLY+NEL+5+7*NPDE) +
C 4 NPDE * 8 + 3 + NV + NXI,
C C NPDE * (7 * (NPOLY+1+NXI) + 8),
C DDASSL TIME INTEGRATION PARAMETERS
C 5 MAXORD = 5, LRW = 40 + (MAXORD+4) * NEQ + NEQ**2,
C 6 LIW = 20 + NEQ )
C
C INTEGER IWORK(LIW), INFO(15), IBAND, M, ITIME, I, IDID, IRESWK,
C 1 IDEV, ITRACE
C DOUBLE PRECISION XBK(IBK), X(NPTS), Y(NEQ), YDOT(NEQ),
C 1 WKRES(NWKRES), RWORK(LRW), XI(1), T, TOUT, RTOL, ATOL,
C 2 ENORM, GERR, CTIME
C EXTERNAL PDECHB, DGEJAC
C COMMON /SDEV2/ ITRACE, IDEV
C M = 2
C T = 0.0D0
C IDEV = 4
C ITRACE = 1
C WRITE(IDEV,9)NPOLY, NEL
C9 FORMAT(' TEST PROBLEM 1'/' ***********'/' POLY OF DEGREE =',I4,
C 1 ' NO OF ELEMENTS = ',I4)
C DO 10 I = 1,IBK
C10 XBK(I) = (I-1.0D0) / (IBK - 1.0D0)
C INITIALISE THE P.D.E. WORKSPACE
C ITIME = 1
C CALL INICHB(NEQ, NPDE, NPTS, X, Y, WKRES, NWKRES, M, T, IBAND,
C 1 ITIME, XBK, IBK, NEL, NPOLY, NV, NXI, XI, IDEV)
C IF(ITIME .EQ. -1)THEN
C WRITE(IDEV, 15)
C15 FORMAT(' INICHB ROUTINE RETURNED ITIME = -1 - RUN HALTED ')
C GOTO 100
C END IF
C SETUP DASSL PARAMETERS
C RTOL = 1.0D-8
C ATOL = 1.0D-8
C DO 20 I = 1,11
C20 INFO(I) = 0
C INFO(11)= 1
C BANDED MATRIX OPTION WHEN INFO(6) = 1
C IF(INFO(6) .EQ. 1)THEN
C IWORK(1) = IBAND
C IWORK(2) = IBAND
C END IF
C30 TOUT = T + 0.1D0
C CALL DDASSL( PDECHB, NEQ, T, Y, YDOT, TOUT, INFO, RTOL, ATOL,
C 1 IDID, RWORK, LRW, IWORK, LIW, WKRES, IRESWK, DGEJAC)
C IF( IDID .LT. 0 )THEN
C DASSL FAILED TO FINISH INTEGRATION.
C WRITE(IDEV,40)T,IDID
C40 FORMAT(' AT TIME T = ',D11.3,' DASSL RETURNED IDID =',I3)
C GOTO 100
C ELSE
C DASSL INTEGRATED TO T = TOUT
C CALL TO POST PROCESSING HERE E.G. SPACE INTERPOLATION.
C CALL ERROR( Y, NPDE, NPTS, X, M, ENORM, GERR, T, RTOL, ATOL,
C 1 ITRACE, WKRES, NWKRES)
C IF(TOUT .LT. 0.99D0 ) GOTO 30
C END IF
C100 CONTINUE
C WRITE(IDEV,110)IWORK(11),IWORK(12),IWORK(13)
C110 FORMAT(' NSTEPS =',I5,' NRESID =',I5,' JAC = ',I4)
C STOP
C END
C SUBROUTINE UVINIT( NPDE, NPTS, X, U, NV,V)
C ROUTINE FOR P.D.E. INITIAL VALUES.
C INTEGER NPDE, NPTS, NV
C DOUBLE PRECISION X(NPTS), U(NPDE,NPTS), V(1), T
C T = 0.0D0
C V(1) IS A DUMMY VARIABLE AS THERE ARE NO COUPLED O.D.E.S
C CALL EXACT( T, NPDE, NPTS, X, U )
C RETURN
C END
C
C SUBROUTINE SPDEFN( T, X, NPTL, NPDE, U, DUDX, UDOT, UTDX, Q, R,
C 1 NV, V, VDOT, IRES)
C ROUTINE TO DESCRIBE THE BODY OF THE P.D.E.
C THE P.D.E. IS WRITEN AS -M M
C Q(X,T,U, U , U , U ) = X (X R(X,T,U,U , U , U ))
C X T TX X T TX X
C THE FUNCTIONS Q AND R MUST BE DEFINED IN THIS ROUTINE.
C DEFINITIONS FOR THE MODEL PROBLEM ARE GIVEN
C NOTE NV = 0 : THERE IS NO O.D.E PART.
C INTEGER NPDE, NPTL, I, J, NV, IRES
C DOUBLE PRECISION T, X(NPTL), U(NPDE,NPTL), DUDX(NPDE,NPTL),
C 1 UDOT(NPDE,NPTL), Q(NPDE,NPTL), R(NPDE,NPTL), V, VDOT,
C 2 UTDX(NPDE,NPTL)
C DO 10 I = 1,NPTL
C R(1,I) = U(1,I) * DUDX(1,I)
C Q(1,I) = U(1,I) * UDOT(1,I) - 5.0D0 * U(1,I)**2
C 1 - 4.0D0 * U(1,I)*DUDX(1,I)*X(I)
C10 CONTINUE
C RETURN
C END
C
C SUBROUTINE SBNDR( T, BETA, GAMMA, U, UX, UDOT, UTDX, NPDE, LEFT,
C 1 NV, V, VDOT, IRES)
C BOUNDARY CONDITIONS ROUTINE
C INTEGER NPDE, NV, IRES
C DOUBLE PRECISION T, BETA(NPDE), GAMMA(NPDE), U(NPDE), C2,
C 1 UX(NPDE), V, VDOT, UDOT(NPDE), UTDX(NPDE)
C LOGICAL LEFT
C IF(LEFT) THEN
C BETA (1) = 1.0D0
C GAMMA(1) = 0.0D0
C ELSE
C BETA (1) = 0.0D0
C GAMMA(1) = U(1) - DEXP( -T )
C BETA (1) = 1.0D0
C GAMMA(1) = - 2.D0 *U(1)**2
C END IF
C RETURN
C END
C
C DUMMY O.D.E. ROUTINE AS NV IS ZERO
C SUBROUTINE SODEFN
C RETURN
C END
C EXACT SOLUTION
C SUBROUTINE EXACT( TIME, NPDE, NPTS, X, U)
C ROUTINE FOR P.D.E. EXACT VALUES (IF KNOWN)
C INTEGER NPDE, NPTS, I
C DOUBLE PRECISION X(NPTS), U(NPDE,NPTS), TIME
C DO 10 I = 1,NPTS
C10 U(1,I) = DEXP( 1.0D0 - X(I)**2 - TIME)
C RETURN
C END
C
C EXAMPLE PROBLEM 3
C *********************
C THIS PROBLEM IS DEFINED BY
C -1
C U = ( C U ) - C * EXP(-2U) + EXP(-U) , X IN (-1,0)
C T 1 X X 1
C AND
C -1
C U = ( C U ) - C * EXP(-2U) + EXP(-U) , X IN (0,1)
C T 2 X X 2
C WHERE
C C = 0.1 AND C = 1.0
C 1 2
C
C THE LEFT BOUNDARY CONDITION AT X =-1 (LEFT = .TRUE. ) IS GIVEN BY
C U(-1,T) = LOG ( - C + T + P)
C 1
C THE RIGHT BOUNDARY CONDITION IS (LEFT = .FALSE.)
C U( 1,T) + (C + T + P ) U = LOG ( - C + T + P) + 1.0D0
C X
C
C THE INITIAL CONDITION IS GIVEN BY THE EXACT SOLUTION ;
C U( X, T ) = LOG ( C X + T + P ) , X IN ( -1, 0)
C 1
C U( X, T ) = LOG ( C X + T + P ) , X IN ( 0, 1)
C 2
C**********************************************************************
C SUBROUTINE UVINIT( NPDE, NPTS, X, U, NV,V)
C ROUTINE FOR P.D.E. INITIAL VALUES.
C INTEGER NPDE, NPTS, NV
C DOUBLE PRECISION X(NPTS), U(NPDE,NPTS), V(1), T
C T = 0.0D0
C V(1) IS A DUMMY VARIABLE AS THERE ARE NO COUPLED O.D.E.S
C CALL EXACT( T, NPDE, NPTS, X, U )
C RETURN
C END
C
C SUBROUTINE SPDEFN( T, X, NPTL, NPDE, U, DUDX, UDOT, UTDX, Q, R,
C 1 NV, V, VDOT, IRES)
C ROUTINE TO DESCRIBE THE BODY OF THE P.D.E.
C THE P.D.E. IS WRITEN AS -M M
C Q(X,T,U, U , U , U ) = X (X R(X,T,U,U , U , U ))
C X T TX X T TX X
C THE FUNCTIONS Q AND R MUST BE DEFINED IN THIS ROUTINE.
C DEFINITIONS FOR THE MODEL PROBLEM ARE GIVEN
C NOTE NV = 0 : THERE IS NO O.D.E PART.
C INTEGER NPDE, NPTL, I, J, NV, IRES
C DOUBLE PRECISION T, X(NPTL), U(NPDE,NPTL), DUDX(NPDE,NPTL),
C 1 UDOT(NPDE,NPTL), Q(NPDE,NPTL), R(NPDE,NPTL), V, VDOT,
C 2 UTDX(NPDE,NPTL), C
C IF(X(1) .LT. 0.0D0 .AND. X(NPTL) .LE. 0.0D0)THEN
C ELEMENT TO LEFT OF THE INTERFACE AT 0.0
C C = 0.1D0
C ELSE
C C = 1.0D0
C END IF
C DO 10 I = 1,NPTL
C R(1,I) = DUDX(1,I) /C
C Q(1,I) = UDOT(1,I) - DEXP(-U(1,I))- DEXP(-2.0D0*U(1,I))* C
C10 CONTINUE
C RETURN
C END
C
C SUBROUTINE SBNDR( T, BETA, GAMMA, U, UX, UDOT, UTDX, NPDE, LEFT,
C 1 NV, V, VDOT, IRES)
C BOUNDARY CONDITIONS ROUTINE
C INTEGER NPDE, NV, IRES
C DOUBLE PRECISION T, BETA(NPDE), GAMMA(NPDE), U(NPDE), C2,
C 1 UX(NPDE), V, VDOT, UDOT(NPDE), UTDX(NPDE)
C LOGICAL LEFT
C IF(LEFT) THEN
C BETA (1) = 0.0D0
C GAMMA(1) = U(1) - DLOG( -0.1 + T + 1.0D0)
C ELSE
C C2 = 1.0D0
C BETA (1) = C2 * ( C2 + T + 1.0D0)
C GAMMA(1) = U(1) - DLOG( C2 + T + 1.0D0) + 1.0D0
C END IF
C RETURN
C END
C
C DUMMY O.D.E. ROUTINE AS NV IS ZERO
C SUBROUTINE SODEFN
C RETURN
C END
C EXACT SOLUTION
C SUBROUTINE EXACT( TIME, NPDE, NPTS, X, U)
C ROUTINE FOR P.D.E. EXACT VALUES (IF KNOWN)
C INTEGER NPDE, NPTS, I, IDERIV
C DOUBLE PRECISION X(NPTS), U(NPDE,NPTS), TIME, C
C COMMON /PROB3/ IDERIV
C IF(IDERIV .EQ. 0)THEN
C DO 10 I = 1,NPTS
C C = 1.0D0
C IF(X(I) .LT. 0.0D0)C = 0.1D0
C10 U(1,I) = DLOG( TIME + 1.0D0 + C * X(I))
C ELSE
C DO 20 I = 1,NPTS
C C = 1.0D0
C IF(X(I) .LT. 0.0D0)C = 0.1D0
C U(1,I) = C / ( TIME + 1.0D0 + C * X(I))
C IF(X(I) .EQ. 0.0D0) U(1,I) = 0.55D0 / ( TIME + 1.0D0 )
C20 CONTINUE
C END IF
C RETURN
C END
C
C
C OTHER EXAMPLE PROBLEMS
C **********************
C
C EXAMPLE PROGRAMS ARE PROVIDED FOR THE POOL EVAPORATION PROBLEM
C DESCRIBED IN THE PAPER THAT ACCOMPANIES THIS CODE AND FOR THE
C
C FOURTH ORDER P.D.E. PROBLEM WRITTEN AS ELLIPTIC-PARABOLIC SYSTEM.
C
C U = K U + UU - U U
C XXT XXXX XXX X XX
C
C-----------------------------------------------------------------------
C
C SOFTWARE STRUCTURE
C ******************
C
C THE SUBROUTINES IN THIS MODULE CAN BE BROKEN DOWN INTO THREE
C PARTS.
C
C 1) INITIALISATION ROUTINES 2) DEFINITION OF D.A.E.S
C
C ---------- ------------
C ! INICHB ! ! PDECHB !
C ---------- ------------
C ! ! ! !
C ---------- ---------------- ! ----------
C ! CSET ! ! DRES OR CRES ! ! ! CHINTR !
C ---------- ---------------- ! ----------
C ! ! ! !
C ---------- ---------- --------- ----------
C ! UVINIT ! ! SPDEFN ! ! SBNDR ! ! SODEFN !
C ---------- ---------- --------- ----------
C
C THE FOUR BOTTOM LEVEL ROUTINES ARE THE USER DEFINED PROBLEM
C DEFINITION ROUTINES.
C
C 3) POST PROCESSING (SPACE INTERPOLATION ).
C
C ----------
C ! INTERC ! THIS ROUTINE CAN BE CALLED BY THE USER.
C ----------
C !
C ----------
C ! INTRCH ! SYSTEM INTERPOLATION ROUTINE ONLY.
C ----------
C
C 4) ERROR MESSAGE HANDLER
C ALL THE ABOVE ROUTINES MAY CALL A GENERAL ERROR MESSAGE
C HANDLING ROUTINE CALLED
C ----------
C ! SCHERR !
C ----------
C
C 5) OPTIONAL ERROR CALCULATION ROUTINE FOR ANALYTIC SOLUTION
C PROBLEMS IN THIS CASE THE USER MAY CALL AN ERROR CALCULATION
C ROUTINE CALLED ERROR WHICH IN TURN CALLS A USER DEFINED
C ROUTINE TO SUPPLY THE ANALYTIC SOLUTION (MUST BE NAMED EXACT)
C ----------
C ! ERROR !
C ----------
C !
C ----------
C ! EXACT !
C ----------
C----------------------------------------------------------------------
C
C INTERFACES TO THE INDIVIDUAL ROUTINES NOW FOLLOW IN THE FOLLOWING
C ORDER
C (1) INICHB
C (2) CSET
C (3) PDECHB
C (4) CHINTR
C (5) INTERC
C (6) INTRCH
C (7) SCHERR
C (8) ERROR
C THE LAST ROUTINE IS A P.D.E. UTILITY TO EVALUATE THE P.D.E.
C ERROR NORMS AND GRID ERRORS FOR PROBLEMS WITH ANALYTIC SOL.
C
C
C***********************************************************************
SUBROUTINE INICHB(NEQN,NPDE,NPTS,X,U,WK,IWK,M,TS,IBAND,ITIME,XBK,
* IBK,NEL,NPOLY,NV,NXI,XI,IDEV)
C***********************************************************************
C ROUTINE FOR INITIALISATION OF CHEBYSHEV GENERALIZED COLLOCATION METHOD
C
C PARAMETER LIST
C ----------------
C NEQN: EMPTY ON ENTRY, ON EXIT IT CONTAINS THE NUMBER OF
C O.D.E.S GENERATED BY THE DISCRETISED FORM OF THE
C P.D.E. , GIVEN BY NPDE*NEL*(NPTL-1) + NPDE + NV.
C
C NPDE NUMBER OF PARABOLIC P.D.E.S IN ONE SPACE DIMENSION
C
C NPTS NUMBER OF SPATIAL GRID POINTS USED IN M.O.L. SOLUTION.
C NOTE THIS SHOULD BE EQUAL TO (NPTL-1)*NEL + 1
C
C X(NPTS) EMPTY ARRAY ON ENTRY . ON EXIT THIS ARRAY
C CONTAINS THE MESH USED IN SEMI-DISCRETISATION
C
C M =0,1,2 IF CARTESIAN CYLINDRICAL OR SPHERICAL POLARS.
C
C U(NEQN) SOLUTION VECTOR EMPTY ON ENTRY CONTAINS INITIAL
C VALUES ON EXIT. THIS ARRAY IS ORDERED AS FOLLOWS.
C U(1) - U(NPDE*NPTS) P.D.E. SOLUTION COMPONENTS.
C U(NPDE*NPTS+1) - U(NEQN) O.D.E. COMPONENTS THAT ARE
C COUPLED TO THE P.D.E.
C
C WK(IWK) REAL WORKSPACE USED TO PASS FOUR MATRICES AND VARIOUS
C USEFUL VECTORS TO THE O.D.E.FUNCTION CALL ROUTINE
C RESID SEE BELOW FOR A DETAILED DESCRIPTION.
C
C TS STARTING LEVEL OF TIME INTEGRATION.
C
C XBK(IBK) REAL ARRAY OF BREAK POINTS IBK = NEL +1 WHERE
C NEL IS THE NUMBER OF SPATIAL ELEMENTS.
C XBK(1) = XLEFT