C ALGORITHM 565
C
C PDETWO/PSTEM/GEARB: SOLUTION OF SYSTEMS OF TWO DIMENSIONAL
C NONLINEAR PARTIAL DIFFERENTIAL EQUATIONS
C
C BY D.K. MELGAARD AND R.F. SINCOVEC
C
C ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE 7,1 (MARCH 1981)
C
C CPDE
C MAIN PROGRAM TO ILLUSTRATE THE USAGE OF THE PDETWO PACKAGE CPDE
C USING EXAMPLES 1, 2, AND 3. CPDE
C CPDE
REAL H,S,X,ERR1,DX,Y,TOUT,REALN1,DY,U1,HUSED, CPDE
* EXP,R,T0,EPS,ERMAX,ABS,WORK,DL,DLI,DEV, CPDE
* U2,U3,ERR2,ERR3,REALN2,REALN3 CPDE
INTEGER IX,NPDE,NSTEP,NX,NFE,NODE,MF,NY, CPDE
* NJE,NQUSED,IY,INDEX,I,IWORK,KODE,IK CPDE
COMMON /GEAR3/ HUSED,NQUSED,NSTEP,NFE,NJE CPDE
COMMON /PROB/ DL,DLI,KODE CPDE
DIMENSION U1(10,10),ERR1(10,10),REALN1(10,10) CPDE
DIMENSION U2(5,31),ERR2(5,31),REALN2(5,31) CPDE
DIMENSION U3(2,5,5),ERR3(2,5,5),REALN3(2,5,5) CPDE
DIMENSION WORK(4186),IWORK(155),X(10),Y(31) CPDE
EQUIVALENCE (U1(1,1),U3(1,1,1)),(ERR1(1,1),ERR3(1,1,1)), CPDE
* (REALN1(1,1),REALN3(1,1,1)) CPDE
EQUIVALENCE (U2(1,1),U3(1,1,1)),(ERR2(1,1),ERR3(1,1,1)), CPDE
* (REALN2(1,1),REALN3(1,1,1)) CPDE
C CPDE
DO 1000 IK=1,3 CPDE
KODE = IK CPDE
IF(KODE.EQ.1) GO TO 100 CPDE
IF(KODE.EQ.2) GO TO 300 CPDE
IF(KODE.EQ.3) GO TO 500 CPDE
100 CONTINUE CPDE
C***********************************************************************
C
C EXAMPLE 1. ELLIPTIC PDE
C
C***********************************************************************
WRITE (6,110)
110 FORMAT(1H1,3X,36HOUTPUT FROM EXAMPLE 1. ELLIPTIC PDE//)
C
C DEFINE THE PROBLEM PARAMETERS.
C
NX=10
NY=10
NPDE=1
NODE=NPDE*NX*NY
MF=22
INDEX=1
T0=0.0
H=0.1E-06
EPS=0.1E-06
DX = 1.0/(FLOAT(NX)-1.0)
DY = 1.0/(FLOAT(NY)-1.0)
DO 120 IX=1,NX
Y(IX)=FLOAT(IX)*DY-DY
120 X(IX)=FLOAT(IX)*DX-DX
IWORK(1) = NPDE
IWORK(2) = NX
IWORK(3) = NY
IWORK(4) = 5
IWORK(5) = 4186
IWORK(6) = 100
C
C CALCULATE THE ACTUAL SOLUTION
C
DO 130 IY = 1,NY
DO 130 IX = 1,NX
R = X(IX)
S = Y(IY)
130 REALN1(IX,IY)=3.0*EXP(R)*EXP(S)*(R-R*R)*(S-S*S)
WRITE (6,140) ((REALN1(IX,IY),IX=1,NX),IY=1,NY)
140 FORMAT (/23H STEADY STATE SOLUTION /10(/3H ,10E11.4))
C
C DEFINE THE INITIAL CONDITION
C
DO 160 IY = 1,NY
DO 160 IX = 1,NX
160 U1(IX,IY)=0.0
WRITE (6,170) NODE,T0,H,EPS,MF,U1
170 FORMAT (//6H NODE ,I3,4H T0 ,F9.2,3H H ,E8.1,5H EPS ,E8.1,4H MF ,
* I2 //16H INITIAL VALUES / 10(/3H ,10E11.2))
C
C SET UP THE LOOP FOR CALLING THE INTEGRATOR AT DIFFERENT TOUT VALUES
C
TOUT=.001
DO 260 I=1,5
WRITE (6,200) TOUT
200 FORMAT (// 5H TOUT,E15.6)
C
C CALL THE INTEGRATOR
C
CALL DRIVEP (NODE,T0,H,U1,TOUT,EPS,MF,INDEX,WORK,IWORK,X,Y)
C
C CHECK ERROR RETURN
C
IF (INDEX .NE. 0) GO TO 1000
C
C OUTPUT THE RESULTS
C
WRITE (6,220) HUSED,NQUSED,NSTEP,NFE,NJE
220 FORMAT (/7H HUSED ,E10.4,7H ORDER ,I3,7H NSTEP ,I6,5H NFE ,I5,
* 5H NJE ,I5)
WRITE (6,230) ((U1(IX,IY),IX=1,NX),IY=1,NY)
230 FORMAT ( /9H U VALUES //(3H ,10E11.4))
C
C DETERMINE THE ERROR
C
ERMAX = 0.0
DO 240 IY = 1,NY
DO 240 IX = 1,NX
ERR1(IX,IY) = U1(IX,IY) - REALN1(IX,IY)
IF (ABS(ERMAX) .LT. ABS(ERR1(IX,IY))) ERMAX=ERR1(IX,IY)
240 CONTINUE
WRITE (6,250) ERMAX
250 FORMAT (/51H MAXIMUM DIFFERENCE FROM THE STEADY STATE SOLUTION //
* (3H ,E11.4))
TOUT=TOUT*10.
260 CONTINUE
C
C************************** END OF EXAMPLE 1. **************************
C
GO TO 1000
300 CONTINUE
C***********************************************************************
C
C EXAMPLE 2. BURGERS EQUATION
C
C***********************************************************************
WRITE (6,310)
310 FORMAT(1H1,3X,40HOUTPUT FROM EXAMPLE 2. BURGERS EQUATION//)
C
C DEFINE THE PROBLEM PARAMETERS.
C
DL=.010
DLI=1.0/(2.0*DL)
NX=5
NY=31
NPDE=1
NODE=NX*NY*NPDE
MF=10
INDEX=1
T0=0.0
H=0.1E-06
EPS=0.1E-03
DY = 1.0/(FLOAT(NY)-1.0)
DX = DY
DO 320 I=1,NX
320 X(I)=FLOAT(I)*DX-DX
DO 330 I=1,NY
330 Y(I)=FLOAT(I)*DY-DY
IWORK(1) = NPDE
IWORK(2) = NX
IWORK(3) = NY
IWORK(4) = 12
IWORK(5) = 2687
IWORK(6) = 155
C
C DEFINE THE INITIAL CONDITION
C
DO 360 IY = 1,NY
DO 360 IX = 1,NX
DEV=EXP((-X(IX)-Y(IY)+T0)*DLI)
360 U2(IX,IY)=DEV/(1.0+DEV)
WRITE (6,370) NODE,T0,H,EPS,MF,((U2(IX,IY),IX=1,NX),IY=1,NY)
370 FORMAT (6H NODE ,I3,4H T0 ,F9.2,3H H ,E8.1,5H EPS ,E8.1,4H MF ,I2
* //16H INITIAL VALUES / 31(/3H , 5E11.4))
C
C SET UP THE LOOP FOR CALLING THE INTEGRATOR AT DIFFERENT TOUT VALUES
C
TOUT=0.0
DO 460 I=1,4
TOUT=TOUT+.10
WRITE (6,400) TOUT
400 FORMAT (// 5H TOUT,E15.6)
C
C CALCULATE THE ACTUAL SOLUTION
C
DO 410 IY = 1,NY
DO 410 IX = 1,NX
DEV=EXP((-X(IX)-Y(IY)+TOUT)*DLI)
410 REALN2(IX,IY)=DEV/(1.0+DEV)
C
C CALL THE INTEGRATOR
C
CALL DRIVEP (NODE,T0,H,U2,TOUT,EPS,MF,INDEX,WORK,IWORK,X,Y)
C
C CHECK ERROR RETURN
C
IF (INDEX .NE. 0) GO TO 1000
C
C OUTPUT THE RESULTS
C
WRITE (6,420) HUSED,NQUSED,NSTEP,NFE,NJE
420 FORMAT (/7H HUSED ,E10.4,7H ORDER ,I3,7H NSTEP ,I6,5H NFE ,I5,
* 5H NJE ,I5)
WRITE (6,430) ((U2(IX,IY),IX=1,NX),IY=1,NY)
430 FORMAT ( /9H U VALUES // (3H ,5E11.4))
C
C DETERMINE THE ERROR
C
ERMAX=0.0
DO 440 IY = 1,NY
DO 440 IX = 1,NX
ERR2(IX,IY)=U2(IX,IY)-REALN2(IX,IY)
IF (ABS(ERMAX) .LT. ABS(ERR2(IX,IY))) ERMAX=ERR2(IX,IY)
440 CONTINUE
WRITE (6,450) ERMAX
450 FORMAT (/13H MAX ERROR U ,E11.4)
WRITE (6,455) ((ERR2(IX,IY),IX=1,NX),IY=1,NY)
455 FORMAT (/6H ERROR // (3H ,5E11.4))
460 CONTINUE
C
C************************** END OF EXAMPLE 2. **************************
C
GO TO 1000
500 CONTINUE
C***********************************************************************
C
C EXAMPLE 3. COUPLED SYSTEM OF PDE*S
C
C***********************************************************************
WRITE (6,510)
510 FORMAT(1H1,3X,47HOUTPUT FROM EXAMPLE 3. COUPLED SYSTEM OF PDE*S//
*)
C
C DEFINE THE PROBLEM PARAMETERS. NPDE,NX, AND NY PRIMARILY DETERMINE
C THE DIMENSIONS FOR THE ARRAYS IN PDETWO AND THE MODIFIED GEARB.
C
NX=5
NY=5
NPDE=2
NODE=NX*NY*NPDE
MF=22
INDEX=1
T0=0.0
H=0.1E-06
EPS=0.1E-03
DX=1.0/(FLOAT(NX)-1.0)
DO 520 IX=1,NX
X(IX)=FLOAT(IX)*DX-DX
520 Y(IX)=X(IX)
IWORK(1) = NPDE
IWORK(2) = NX
IWORK(3) = NY
IWORK(4) = 5
IWORK(5) = 2308
IWORK(6) = 50
C
C DEFINE THE INITIAL CONDITIONS
C
DO 560 IY=1,NY
DO 560 IX=1,NX
U3(1,IX,IY)=X(IX)+Y(IY)
560 U3(2,IX,IY)=2.0*X(IX)+3.0*Y(IY)
WRITE (6,570) NODE,T0,H,EPS,MF,((U3(1,IX,IY),IX=1,NX),IY=1,NY),
* ((U3(2,IX,IY),IX=1,NX),IY=1,NY)
570 FORMAT (6H NODE ,I3,4H T0 ,F9.2,3H H ,E8.1,5H EPS ,E8.1,4H MF ,I2
* //19H INITIAL U1 VALUES / 5(/3H ,5E11.4)
* //19H INITIAL U2 VALUES / 5(/3H ,5E11.4))
C
C SET UP THE LOOP FOR CALLING THE INTEGRATOR AT DIFFERENT TOUT VALUES
C
TOUT=.25
DO 660 I=1,4
WRITE (6,600) TOUT
600 FORMAT (//5H TOUT ,E15.6)
C
C CALCULATE THE ACTUAL SOLUTION
C
DO 610 IY=1,NY
DO 610 IX=1,NX
REALN3(2,IX,IY)=TOUT+2.0*X(IX)+3.0*Y(IY)
610 REALN3(1,IX,IY)=TOUT +X(IX)+Y(IY)
C
C CALL THE INTEGRATOR
C
CALL DRIVEP (NODE,T0,H,U3,TOUT,EPS,MF,INDEX,WORK,IWORK,X,Y)
C
C CHECK ERROR RETURN
C
IF (INDEX .NE. 0) GO TO 1000
C
C OUTPUT THE RESULTS
C
WRITE (6,620) HUSED,NQUSED,NSTEP,NFE,NJE
620 FORMAT (/7H HUSED ,E10.4,7H ORDER ,I3,7H NSTEP ,I6,5H NFE ,I5,
* 5H NJE ,I5)
WRITE (6,630) ((U3(1,IX,IY),IX=1,NX),IY=1,NY)
630 FORMAT (/11H U1 VALUES //(3H ,5E11.4))
WRITE (6,635) ((U3(2,IX,IY),IX=1,NX),IY=1,NY)
635 FORMAT (/11H U2 VALUES //(3H ,5E11.4))
C
C DETERMINE THE ERROR
C
DO 640 IY=1,NY
DO 640 IX=1,NX
ERR3(1,IX,IY)=U3(1,IX,IY)-REALN3(1,IX,IY)
ERR3(2,IX,IY)=U3(2,IX,IY)-REALN3(2,IX,IY)
640 CONTINUE
WRITE (6,650) ((ERR3(1,IX,IY),IX=1,NX),IY=1,NY)
650 FORMAT (/9H ERROR U1 //(3H ,5E11.4))
WRITE (6,655) ((ERR3(2,IX,IY),IX=1,NX),IY=1,NY)
655 FORMAT (/9H ERROR U2 //(3H ,5E11.4))
TOUT=TOUT+.25
660 CONTINUE
C
C************************** END OF EXAMPLE 3. **************************
C
1000 CONTINUE
STOP
END
SUBROUTINE BNDRYV (T,X,Y,U,AV,BV,CV,NPDE) BNDRYV
C
C DEFINE THE VERTICAL BOUNDARY CONDITIONS
C
REAL T,U,X,Y,BV,AV,CV
INTEGER NPDE
COMMON /PROB/ DL,DLI,KODE
DIMENSION U(NPDE),AV(NPDE),BV(NPDE),CV(NPDE)
IF(KODE.EQ.1) GO TO 100
IF(KODE.EQ.2) GO TO 200
IF(KODE.EQ.3) GO TO 300
C
C EXAMPLE 1. ELLIPTIC PDE
C
100 CONTINUE
AV(1) = 1.0
BV(1) = 0.0
CV(1) = 0.0
GO TO 400
C
C EXAMPLE 2. BURGERS EQUATION
C
200 CONTINUE
AV(1) = 1.0
BV(1) = 0.0
DEV = EXP((-X-Y+T)*DLI)
CV(1) = DEV/(1.0+DEV)
GO TO 400
C
C EXAMPLE 3. COUPLED SYSTEM OF PDE*S
C
300 CONTINUE
IF (X .NE. 0.0) GO TO 310
AV(1)=1.0
BV(1)=0.0
CV(1)=T+Y
AV(2)=1.0
BV(2)=0.0
CV(2)=T+3.0*Y
GO TO 400
310 CONTINUE
AV(1)=1.0
BV(1)=0.0
CV(1)=T+1.0+Y
AV(2)=0.0
BV(2)=1.0
CV(2)=2.0
400 CONTINUE
RETURN
END
SUBROUTINE BNDRYH (T,X,Y,U,AH,BH,CH,NPDE) BNDRYH
C
C DEFINE THE HORIZONTAL BOUNDARY CONDITIONS
C
REAL T,U,X,Y,BH,AH,CH
INTEGER NPDE
COMMON /PROB/ DL,DLI,KODE
DIMENSION U(NPDE),AH(NPDE),BH(NPDE),CH(NPDE)
IF(KODE.EQ.1) GO TO 100
IF(KODE.EQ.2) GO TO 200
IF(KODE.EQ.3) GO TO 300
C
C EXAMPLE 1. ELLIPTIC PDE
C
100 CONTINUE
AH(1) = 1.0
BH(1) = 0.0
CH(1) = 0.0
GO TO 400
C
C EXAMPLE 2. BURGERS EQUATION
C
200 CONTINUE
AH(1) = 1.0
BH(1) = 0.0
DEV = EXP((-X-Y+T)*DLI)
CH(1) = DEV/(1.0+DEV)
GO TO 400
C
C EXAMPLE 3. COUPLED SYSTEM OF PDE*S
C
300 CONTINUE
IF (Y .NE. 0.0) GO TO 310
AH(1)=0.0
BH(1)=1.0
CH(1)=1.0
AH(2)=1.0
BH(2)=1.0
CH(2)=T+2.0*X+3.0
GO TO 400
310 CONTINUE
AH(1)=0.0
BH(1)=U(2)
CH(1)=T+2.0*X+3.0
AH(2)=1.0
BH(2)=0.0
CH(2)=T+2.0*X+3.0
400 CONTINUE
RETURN
END
SUBROUTINE DIFFH (T,X,Y,U,DH,NPDE) DIFFH
C
C DEFINE THE HORIZONTAL DIFFUSION COEFFICIENTS
C
REAL T,U,X,Y,DH
INTEGER NPDE
COMMON /PROB/ DL,DLI,KODE
DIMENSION U(NPDE),DH(NPDE,NPDE)
IF(KODE.EQ.1) GO TO 100
IF(KODE.EQ.2) GO TO 200
IF(KODE.EQ.3) GO TO 300
C
C EXAMPLE 1. ELLIPTIC PDE
C
100 CONTINUE
DH(1,1) = 1.0
GO TO 400
C
C EXAMPLE 2. BURGERS EQUATION
C
200 CONTINUE
DH(1,1) = 1.0
GO TO 400
C
C EXAMPLE 3. COUPLED SYSTEM OF PDE*S
C
300 CONTINUE
DH(1,1)=0.0
DH(1,2)=0.0
DH(2,1)=U(2)
DH(2,2)=U(1)
400 CONTINUE
RETURN
END
SUBROUTINE DIFFV (T,X,Y,U,DV,NPDE) DIFFV
C
C DEFINE THE VERTICAL DIFFUSION COEFFICIENTS
C
REAL T,U,X,Y,DV
INTEGER NPDE
COMMON /PROB/ DL,DLI,KODE
DIMENSION U(NPDE),DV(NPDE,NPDE)
IF(KODE.EQ.1) GO TO 100
IF(KODE.EQ.2) GO TO 200
IF(KODE.EQ.3) GO TO 300
C
C EXAMPLE 1. ELLIPTIC PDE
C
100 CONTINUE
DV(1,1) = 1.0
GO TO 400
C
C EXAMPLE 2. BURGERS EQUATION
C
200 CONTINUE
DV(1,1) = 1.0
GO TO 400
C
C EXAMPLE 3. COUPLED SYSTEM OF PDE*S
C
300 CONTINUE
DV(1,1)=U(2)
DV(1,2)=U(1)
DV(2,1)=0.0
DV(2,2)=0.0
400 CONTINUE
RETURN
END
SUBROUTINE F(T,X,Y,U,UX,UY,DUXX,DUYY,DUDT,NPDE) F
C
C DEFINE THE PDE
C
REAL T,U,X,Y,UX,UY,DUXX,DUYY,DUDT,EXP
INTEGER NPDE
COMMON /PROB/ DL,DLI,KODE
DIMENSION U(NPDE),UX(NPDE),UY(NPDE),DUXX(NPDE,NPDE),
* DUYY(NPDE,NPDE),DUDT(NPDE)
IF(KODE.EQ.1) GO TO 100
IF(KODE.EQ.2) GO TO 200
IF(KODE.EQ.3) GO TO 300
C
C EXAMPLE 1. ELLIPTIC PDE
C
100 CONTINUE
DUDT(1) = DUXX(1,1) + DUYY(1,1) - 6.0*X*Y*EXP(X)*EXP(Y)*
* (X*Y+X+Y-3.0)
GO TO 400
C
C EXAMPLE 2. BURGERS EQUATION
C
200 CONTINUE
DUDT(1) = -U(1)*(UX(1)+UY(1))+DL*(DUXX(1,1)+DUYY(1,1))
GO TO 400
C
C EXAMPLE 3. COUPLED SYSTEM OF PDE*S
C
300 CONTINUE
DUDT(1)=U(2)*(DUYY(1,1)+DUYY(1,2))-3.0*U(2)*UX(2)+4.0*UY(1)-UY(2)
DUDT(2)=DUXX(2,1)+DUXX(2,2)-UY(2)
400 CONTINUE
RETURN
END
SUBROUTINE PSETM (NPDE,NX,NY,X,Y,U,UMAX,USAVE,DUDTR,DUDT,PW,CON, PSETM
* MITER,IER,NEBAND,WORK,IWORK)
C
C***********************************************************************
C
C PSETM IS INTENDED TO BE USED IN CONJUNCTION WITH THE
C INTERFACE PDETWO FOR THE SOLUTION OF SECOND ORDER TIME DEP-
C ENDENT PARTIAL DIFFERENTIAL EQUATIONS (PDE*S) DEFINED OVER A
C RECTANGULAR REGION. PSETM IS SPECIFICALLY DESIGNED TO REPLACE
C THE ROUTINE PSETB USED IN THE ORDINARY DIFFERENTIAL EQUATION
C INTEGRATOR GEARB (1) IN ORDER TO MINIMIZE THE COMPUTATIONS
C REQUIRED TO GENERATE THE JACOBIAN MATRIX.
C
C PSETM IS CALLED BY STIFFP, A ROUTINE IN THE MODIFIED GEARB,
C WHEN THE INTEGRATION METHOD (MITER=1 OR 2) REQUIRES A JACOBIAN
C MATRIX. WHEN MITER=1, THE JACOBIAN IS DETERMINED BY A USER
C DEFINED ROUTINE, PDB. THIS METHOD SHOULD BE AVOIDED SINCE IT RE-
C QUIRES THE USER TO COMPLETELY UNDERSTAND PDETWO. IF THIS
C METHOD IS USED, PSETM OFFERS NO ADVANTAGE OVER PSETB. IF
C MITER=2, THE JACOBIAN IS APPROXIMATED BY FINITE DIFFERENCES,
C USING THE APPROXIMATION J = (DUDT(U+R) - DUDT(U)) / R,
C WHERE DUDT(U+R) AND DUDT(U) ARE THE VALUES OF THE RIGHT HAND
C SIDE OF THE PDE*S EVALUATED AT U+R AND U RESPECTIVELY, AND R
C IS SOME SMALL NUMBER. PSETM TAKES ADVANTAGE OF THE STRUCTURE
C OF THE JACOBIAN MATRIX REQUIRED BY PDETWO TO REDUCE THE COMPU-
C TATIONS NEEDED TO GENERATE THE MATRIX (REQUIRING ONLY 5 * NPDE
C CALLS TO PDETWO). WITH EACH CALL TO PDETWO, PSETM DETERMINES
C THE ENTRIES IN THE JACOBIAN FOR THE MESH POINTS IN THE FOLLOWING
C PATTERN ..
C
C X O O O O X O O O O
C O O O X O O O O X O
C O X O O O O X O O O
C O O O O X O O O O X
C O O X O O O O X O O
C X O O O O X O O O O
C
C WHERE X REPRESENTS THE POINTS FOR WHICH THE JACOBIAN IS
C APPROXIMATED. MITER=2 IS THE RECOMMENDED METHOD.
C
C TO CREATE A DOUBLE PRECISION VERSION OF PSETM..
C CHANGE THE REAL STATEMENT BELOW AND CHANGE THE SINGLE PRECISION
C FUNCTIONS ABS, SQRT, AND AMAX1.
C
C***********************************************************************
C THE PARAMETERS
C***********************************************************************
C
C NPDE IS THE NUMBER OF PARTIAL DIFFERENTIAL EQUATIONS.
C NX IS THE NUMBER OF MESH POINTS IN THE HORIZONTAL
C DIRECTION.
C NY IS THE NUMBER OF MESH POINTS IN THE VERTICAL
C DIRECTION.
C X ARE THE MESH POINTS IN THE HORIZONTAL DIRECTION.
C Y ARE THE MESH POINTS IN THE VERTICAL DIRECTION.
C U CONTAINS THE CURRENT SOLUTIONS OF THE PARTIAL
C DIFFERENTIAL EQUATIONS.
C UMAX IS THE MAXIMUM U VALUES, WHICH ARE USED TO SCALE
C THE VALUE IN R AND FOR ERROR CONTROL IN GEARB.
C USAVE IS A TEMPORY STORE FOR U VALUES.
C DUDTR IS THE VALUE OF DUDT(U+R) RETURNED FROM PDETWO.
C DUDT IS THE VALUE OF DUDT(U), WHICH IS PASSED TO PSETM
C FROM STIFFP.
C PW STORES THE APPROXIMATION FOR THE JACOBIAN.
C CON IS THE CONSTANT (-H*EL(1)).
C MITER INDICATES THE TYPE OF ITERATION METHOD BEING USED
C BY THE INTEGRATOR.
C MITER = 0 MEANS FUNCTIONAL ITERATION.
C MITER = 1 MEANS THE CHORD METHOD WITH AN ANALYTIC
C JACOBIAN SUPPLIED IN THE USER DEFINED
C ROUTINE PDB. THIS METHOD IS IN GEARB,
C BUT IT SHOULD BE AVOIDED.
C MITER = 2 MEANS THE CHORD METHOD WITH THE JACO-
C BIAN CALCULATED IN PSETM. THIS IS THE
C ONLY METHOD WHERE PSETM IS USEFUL.
C MITER = 3 MEANS THE CHORD METHOD WITH THE JACO-
C BIAN REPLACED BY A DIAGONAL APPROXI-
C MATION BASED ON A DIRECTIONAL DERIV-
C ATIVE. THIS METHOD IS IN GEARB, BUT
C IS NOT GENERALLY USEFUL IN SOLVING PDE*S.
C IER IS AN ERROR INDICATOR USED IN THE ROUTINE DECBR.
C NEBAND = 3 * ML + 1
C WORK PROVIDES WORKING STORAGE FOR DYNAMIC DIMENSIONING
C OF THE ARRAYS IN PDETWO AND THE MODIFIED GEARB.
C IWORK IS THE PIVOT VECTOR USED IN THE LU DECOMPOSITION.
C
C***********************************************************************
C THE VARIABLES
C***********************************************************************
C
C R IS A SMALL INCREMENT USED IN EVALUATING THE
C FUNCTION NEAR THE CURRENT SOLUTION U ( I.E.
C DUDT (U+R)).
C R0 IS A LOWER BOUND ON THE SIZE OF R.
C D IS A SMALL NUMBER USED TO APPROXIMATE THE
C DERIVATIVE.
C NXNPDE = NX * NPDE
C
C T IS THE TIME BEING USED FOR THE INTEGRATION.
C H IS THE CURRENT TIME STEP SIZE BEING USED IN THE
C INTEGRATION.
C UROUND IS THE UNIT ROUNDOFF OF THE MACHINE.
C EPSJ IS SQRT(UROUND).
C ML IS THE WIDTH OF THE LOWER (AND UPPER) HALF OF THE
C BAND OF THE JACOBIAN. ML = (NX + 1) * NPDE - 1.
C IW1 - IW27 ARE THE SUBSCRIPT POINTERS USED TO INDICATE THE
C STORAGE LOCATIONS IN WORK OF ARRAYS IN PDETWO AND
C THE MODIFIED GEARB.
C
C***********************************************************************
C THE ROUTINES CALLED BY PSETM
C***********************************************************************
C
C PDB(NODE,T,U,PW,NODE,ML,ML) IS A USER DEFINED ROUTINE WHICH DEFINES
C THE JACOBIAN IN PW IF MITER=1. SINCE THIS METHOD IS
C NOT RECOMMENDED, A DUMMY ROUTINE IS PROVIDED.
C DECBR(NEBAND,NODE,ML,ML,PW,IWORK,IER) COMPUTES THE LU DECOMPOSITION
C OF THE JACOBIAN FOR THE DIRECT SOLUTION OF THE
C JACOBIAN.
C PDETWO(NPDE,NX,NY,X,Y,T,U,...) IS THE INTERFACE WHICH DESCRETIZES
C THE SPATIAL VARIABLES OF THE TWO DIMENSIONAL SYSTEM OF
C PDE*S, THEREBY EVALUATING THE RIGHT HAND SIDE OF
C THE SYSTEM OF ODE*S.
C
C***********************************************************************
C REFERENCES
C***********************************************************************
C
C 1. A. C. HINDMARSH, GEARB.. SOLUTION OF ORDINARY DIFFERENTIAL
C EQUATIONS HAVING BANDED JACOBIAN, UCID-30059 REV. 2,
C LAWRENCE LIVERMORE LABORATORY, P.O.BOX 808, LIVERMORE,
C CA 94550, JUNE 1977.
C
C***********************************************************************
C
REAL U,UMAX,USAVE,DUDTR,DUDT,PW,CON,R,R0,D,T,H,UROUND
* ,EPSJ,DUMMY,WORK
COMMON /GEAR1/ T,H,DUMMY(3),UROUND,IDUMMY(4)
COMMON /GEAR2/ EPSJ,ML,TDUMY(5)
COMMON /ISTORE/ IW1,IW2,IW3,IW4,IW5,IW6,IW7,IW8,IW9,IW10,IW11,
* IW12,IW13,IW14,IW15,IW16,IW17,IW18,IW19,IW20,IW21,
* IW22,IW23,IW24,IW25,IW26,IW27
DIMENSION UMAX(NPDE,NX,NY),USAVE(NPDE,NX,NY),DUDTR(NPDE,NX,NY),
* DUDT(NPDE,NX,NY),PW(NEBAND,1),U(NPDE,NX,NY),
* WORK(1),IWORK(1),X(1),Y(1)
C
C SET THE NUMBER OF ORDINARY DIFFERENTIAL EQUATIONS
C
NODE = NPDE*NX*NY
C
C IF MITER = 1, CALL PDB TO GENERATE THE JACOBIAN
C
IF (MITER .EQ. 2) GO TO 20
CALL PDB (NODE,T,U,PW,NODE,ML,ML)
DO 10 J=1,NODE
DO 10 I=1,NEBAND
10 PW(I,J)=PW(I,J)*CON
GO TO 150
C
C IF MITER = 2, APPROXIMATE THE JACOBIAN WITH FINITE DIFFERENCES
C
20 D = 0.0
DO 30 K=1,NY
DO 30 J=1,NX
DO 30 I=1,NPDE
D=D+DUDT(I,J,K)**2
30 USAVE(I,J,K)=U(I,J,K)
R0 = ABS(H)*SQRT(D)*1.E03*UROUND
ISIGN=-1
NXNPDE=NX*NPDE
DO 50 I=1,NODE
DO 50 J=1,NEBAND
50 PW(J,I)=0.0
DO 140 M=1,5
C
C APPROXIMATE THE JACOBIAN FOR EACH PDE
C
DO 140 IPDE=1,NPDE
C
C DETERMINE THE MESH POINTS FOR WHICH THE APPROXIMATIONS ARE TO
C BE CALCULATED
C
DO 70 IY=1,NY
IXSTRT=1+MOD(M+2*IY,5)
IF(IXSTRT.GT.NX) GO TO 70
DO 60 IX=IXSTRT,NX,5
R=AMAX1(EPSJ*UMAX(IPDE,IX,IY),R0)
60 U(IPDE,IX,IY)=U(IPDE,IX,IY)+R
70 CONTINUE
C
C DETERMINE DUDT(U+R)
C
CALL PDETWO (NPDE,NX,NY,X,Y,T,U,DUDTR,WORK,WORK(IW1),
* WORK(IW2),WORK(IW3),WORK(IW4),WORK(IW5),WORK(IW6),
* WORK(IW7),WORK(IW8),WORK(IW9),WORK(IW10),WORK(IW11),
* WORK(IW12),WORK(IW13),WORK(IW14),WORK(IW15),WORK(IW16),
* WORK(IW17),WORK(IW18),WORK(IW19),WORK(IW20),WORK(IW21))
DO 140 IYY=1,NY
IXSTRT=1+MOD(M+2*IYY,5)
IF(IXSTRT.GT.NX) GO TO 140
DO 130 IXX=IXSTRT,NX,5
R=AMAX1(EPSJ*UMAX(IPDE,IXX,IYY),R0)
D=CON/R
IY=IYY
IX=IXX
IBEG1=ML+IPDE+1
IBEG2=(IX-1)*NPDE+(IY-1)*NXNPDE
I1=IBEG1
I2=IBEG2
C
C WITH FIVE POINT CENTERED DIFFERENCING, THE ONLY NONZERO
C ENTRIES IN THE JACOBIAN ARE GENERATED FROM THE FIVE MESH
C POINTS USED IN THE DIFFERENCING.
C
DO 120 K=1,5
C
C IF K=1, DETERMINE THE JACOBIAN APPROXIMATIONS FOR THE MESH
C POINT (X(IXX),Y(IYY))
C
GO TO (100,80,80,90,90),K
C
C
C IF K=2 OR K=3, DETERMINE THE JACOBIAN APPROXIMATIONS FOR THE
C MESH POINTS (X(IXX-1),Y(IYY)) AND (X(IXX+1),Y(IYY))
C
80 IX=IXX-ISIGN
IF (IX .EQ. 0 .OR. IX .GT. NX) GO TO 120
I1=IBEG1+NPDE*ISIGN
I2=IBEG2-NPDE*ISIGN
GO TO 100
C
C IF K=4 OR K=5, DETERMINE THE JACOBIAN APPROXIMATIONS FOR THE
C MESH POINTS (X(IXX),Y(IYY-1)) AND (X(IXX),Y(IYY+1))
C
90 IX=IXX
IY=IYY-ISIGN
IF (IY .EQ. 0 .OR. IY .GT. NY) GO TO 120
I1=IBEG1+NXNPDE*ISIGN
I2=IBEG2-NXNPDE*ISIGN
C
C CALCULATE THE FINITE DIFFERENCE APPROXIMATIONS
C (DUDT(U+R) - DUDT(U)) / R
C
100 DO 110 J=1,NPDE
K1=I1-J
K2=I2+J
110 PW(K1,K2)=(DUDTR(J,IX,IY)-DUDT(J,IX,IY))*D
120 ISIGN=-ISIGN
U(IPDE,IXX,IYY)=USAVE(IPDE,IXX,IYY)
130 USAVE(IPDE,IXX,IYY)=0.0
140 CONTINUE
150 DO 160 I=1,NODE
160 PW(ML+1,I)=PW(ML+1,I)+1.0
C
C DO THE LU DECOMPOSITION ON THE JACOBIAN
C
CALL DECBR (NEBAND,NODE,ML,ML,PW,IWORK,IER)
RETURN
END
SUBROUTINE PDETWO (NPDE,NX,NY,X,Y,T,U,DUDT,DVS,DVST,DV,DH, PDETWO
* AH,BH,CH,AV,BV,CV,UX,UY,DXI,DXIR,DXIC,
* XAVG,UDVA,UDHR,UDVB,UDHL,UAVGH,UAVGV)
C
C***********************************************************************
C
C PDETWO IS AN INTERFACE DESIGNED TO SOLVE A SYSTEM OF SECOND
C ORDER TIME DEPENDENT PARTIAL DIFFERENTIAL EQUATIONS (PDE*S)
C DEFINED OVER A RECTANGLE IN CONJUNCTION WITH AN ORDINARY
C DIFFERENTIAL EQUATION (ODE) INTEGRATOR PACKAGE (PRIMARILY A MODI-
C FIED VERSION OF GEARB(1)). IT IS AN EXTENSION OF PDEONE (2,3), AN
C INTERFACE DEVELOPED FOR THE SOLUTION OF ONE-DIMENSIONAL SYSTEMS OF
C PDE*S. THIS EXTENSION CONSISTS OF CHANGING FROM THREE POINT
C DIFFERENCING TO FIVE POINT DIFFERENCING IN ORDER TO APPROXIMATE
C THE SPATIAL DERIVATIVES IN THE TWO DIMENSIONAL PDE SYSTEM. TO
C SOLVE A PDE SYSTEM USING THIS SOFTWARE INTERFACE IN CONJUNCTION
C WITH AN ODE INTEGRATOR, THE USER IS REQUIRED ONLY TO PROVIDE
C THE SPATIAL MESH AND DEFINE THE PROBLEM TO BE EVALUATED. THE
C INTERFACE WILL THEN FORM AND EVALUATE A SEMIDISCRETE APPROXI-
C MATION OF THIS SYSTEM OF PDE*S AND THUS PERMIT ONE TO SOLVE
C THE ORIGINAL PDE SYSTEM AS A SYSTEM OF TIME DEPENDENT ODE*S
C USING AN ODE INTEGRATOR.
C
C TO CREATE A DOUBLE PRECISION VERSION OF PDETWO..
C CHANGE THE REAL STATEMENT BELOW.
C
C***********************************************************************
C PROBLEM DEFINITION
C***********************************************************************
C
C PDETWO IS DESIGNED TO BE USED TO SOLVE A COUPLED SYSTEM OF
C NPDE PDE*S. THE L-TH PDE OF THIS SYSTEM IS OF THE FORM ..
C
C DUDT(L) = F (T, X, Y, U, UX, UY, DUXX, DUYY )
C
C WHERE
C
C DUDT(L) REPRESENTS THE FIRST PARTIAL OF THE L-TH COMPONENT
C OF U WITH RESPECT TO T,
C T IS THE CURRENT TIME,
C X,Y DEFINE THE POSITION IN THE HORIZONTAL AND VERTICAL
C DIRECTIONS RESPECTIVELY,
C
C U = ( U(1), ... , U(NPDE) ),
C UX = ( UX(1), ... , UX(NPDE) ),
C UY = ( UY(1), ... , UY(NPDE) ),
C
C DUXX IS AN NPDE BY NPDE ARRAY SUCH THAT
C D
C DUXX(L,K) = -- (DH(L,K)*UX(K)),
C DX
C THE L-TH ROW OF DUXX CORRESPONDS TO THE L-TH PDE,
C
C DUYY IS AN NPDE BY NPDE ARRAY SUCH THAT
C D
C DUYY(L,K) = -- (DV(L,K)*UY(K)),
C DY
C THE L-TH ROW OF DUYY CORRESPONDS TO THE L-TH PDE,
C
C AND
C
C U(K) IS THE VALUE OF THE SOLUTION FOR THE K-TH PDE AT
C (T, X, Y),
C UX(K) IS THE FIRST PARTIAL DERIVATIVE OF U(K) WITH RESPECT
C TO X,
C UY(K) IS THE FIRST PARTIAL DERIVATIVE OF U(K) WITH RESPECT
C TO Y,
C DH(L,K) IS THE DIFFUSION COEFFICIENT IN THE HORIZONTAL
C DIRECTION FOR THE L-TH PDE ASSOCIATED WITH U(K),
C DV(L,K) IS THE DIFFUSION COEFFICIENT IN THE VERTICAL
C DIRECTION FOR THE L-TH PDE ASSOCIATED WITH U(K).
C
C DH(L,K) AND DV(L,K) MAY BE FUNCTIONS OF T, X, Y AND U.
C
C HORIZONTAL BOUNDARY CONDITIONS
C
C AH(L)*U(L) + BH(L)*UY(L) = CH(L)
C
C VERTICAL BOUNDARY CONDITIONS
C
C AV(L)*U(L) + BV(L)*UX(L) = CV(L)
C
C AH(L), BH(L) AND CH(L) (AV(L), BV(L) AND CV(L)) ARE AT LEAST
C PIECEWISE CONTINUOUS FUNCTIONS OF THEIR RESPECTIVE VARIABLES.
C IF BH(L) .NE. 0 (BV(L) .NE. 0) THEN AH(L), BH(L) AND CH(L)
C (AV(L), BV(L) AND CV(L)) MAY BE FUNCTIONS OF T, X, Y AND U
C BUT OTHERWISE THEY MAY ONLY BE FUNCTIONS OF T, X AND Y. NULL
C BOUNDARY CONDITIONS ARE NOT ALLOWED SINCE PDETWO WAS NOT
C DESIGNED TO SOLVE HYPERBOLIS PDE*S. A ZERO DIVIDE WILL OCCUR
C IF A NULL BOUNDARY CONDITION IS SPECIFIED.
C
C INITIAL CONDITIONS
C
C THE INITIAL SOLUTION IS DEFINED FOR EACH U(K) TO BE A
C KNOWN FUNCTION OF X AND Y FOR THE INITIAL TIME T = T0.
C THE INITIAL CONDITIONS NEED NOT BE CONSISTENT WITH THE
C BOUNDARY CONDITIONS AND MAY CONTAIN DISCONTINUITIES.
C
C***********************************************************************
C USER SUPPLIED ROUTINES
C***********************************************************************
C
C THE USER MUST PROVIDE A MAIN PROGRAM AND FIVE SUBROUTINES
C BNDRYH, BNDRYV, DIFFV, DIFFH AND F. THESE ROUTINES ARE ALL THAT
C IS REQUIRED TO DEFINE COMPLETELY THE PROBLEM GIVEN ABOVE.
C
C 1) THE MAIN PROGRAM DEFINES THE SPATIAL MESH, THE INITIAL
C CONDITIONS AND THE PARAMETERS FOR THE ODE INTEGRATOR, CALLS
C THE INTEGRATOR AND PRINTS OR PLOTS THE RESULTS. THE MAIN
C PROGRAM SHOULD BE CONSTRUCTED AS FOLLOWS ..
C
C DIMENSION U(***,*,**),X(*),Y(**)
C DIMENSION WORK(*****),IWORK(****)
C
C FOR THE DIMENSIONS ABOVE ENTER THE ACTUAL NUMERICAL
C VALUES FOR * = NX, ** = NY, *** = NPDE , **** = NODE, AND
C ***** = NPDE * ( NPDE*(NX*2+3) + 13 + NX ) + NX*4 + 4*NODE
C + LPW + LU
C WHERE
C NODE = NPDE*NX*NY
C LPW = 1 IF MITER = 0
C = (3*(NX+1)*NPDE-2)*NODE IF MITER = 1,2
C = NODE IF MITER = 3
C LU = NODE*(MORDER+1)
C SEE MF AND MORDER BELOW FOR THE DEFINITION OF MITER AND
C MORDER.
C
C DEFINE ..
C
C 1) THE NUMBER OF MESH POINTS IN THE X (HORIZONTAL) DIRECTION,
C NX .GE. 3 AND THE NUMBER OF MESH POINTS IN THE Y (VERTICAL)
C DIRECTION, NY .GE. 3 ( TO CONSERVE STORAGE ORIENT THE
C PROBLEM SO THAT NX .LE. NY ),
C 2) THE MESH POINTS (I.E. X(1) .LT. X(2) .LT. ... .LT. X(NX)
C AND Y(1) .LT. Y(2) .LT. ... .LT. Y(NY)),
C 3) NPDE, THE NUMBER OF PARTIAL DIFFERENTIAL EQUATIONS,
C 4) INITIAL VALUES FOR ALL OF U,
C 5) MORDER, THE MAXIMUM ORDER OF THE METHOD USED IN THE
C MODIFIED GEARB. MORDER MUST BE LESS THAN OR EQUAL TO 12
C IF METH = 1 OR IF METH = 2, IT MUST BE LESS THAN OR EQUAL
C TO 5. SEE MF BELOW FOR THE DEFINITION OF METH.
C 6) THE PARAMETERS FOR THE CALL TO THE INTEGRATOR, AND
C 7) THE FIRST SIX LOCATIONS OF THE ARRAY IWORK AS FOLLOWS..
C IWORK(1) = NPDE
C IWORK(2) = NX
C IWORK(3) = NY
C IWORK(4) = MORDER
C IWORK(5) = NRWK
C IWORK(6) = NRIWK
C WHERE NRWK AND NRIWK ARE THE LENGTHS OF THE ARRAYS WORK
C AND IWORK RESPECTIVELY. THEY SHOULD BE AT LEAST EQUAL
C TO ***** AND **** AS DEFINED IN THE DESCRIPTION OF THE
C DIMENSIONS.
C
C IN DETERMINING THE MESH SPACING, THE USER SHOULD BE AWARE
C THAT THE DIFFERENCE APPROXIMATIONS AT THE INTERIOR POINTS
C ARE SECOND ORDER ONLY FOR UNIFORM MESHES. FOR NONUNIFORM
C MESHES, ONLY FIRST ORDER APPROXIMATIONS SHOULD BE EXPECTED.
C IN UNUSUAL SITUATIONS (SEE REMARK IN SECTION 3 OF THE
C ACCOMPANYING PAPER) INVOLVING COUPLED SYSTEMS OF EQUATIONS,
C CERTAIN APPROXIMATIONS AT THE BOUNDARY ARE ONLY FIRST ORDER.
C IN THESE UNUSUAL SITUATIONS THE ERROR CAN BE MINIMIZED BY
C CHOOSING SMALL MESH SPACINGS NEXT TO THE BOUNDARY. IN ANY
C CASE THE USER IS ADVISED TO CHOOSE A MESH WHICH IS LOCALLY
C NEARLY UNIFORM.
C
C FOR THE MODIFIED GEARB THE PARAMETERS INCLUDE THE DESIRED
C OUTPUT TIME (TOUT), THE DESIRED LOCAL ACCURACY (EPS), THE
C NUMBER OF ODE*S (NODE = NPDE*NX*NY), THE INITIAL TIME (T0),
C THE INITIAL TIME STEP SIZE (H),THE TYPE OF CALL BEING MADE
C (INDEX) AND THE TYPE OF INTEGRATION METHOD DESIRED (MF).
C MF HAS TWO DECIMAL DIGITS, METH AND MITER (MF = 10*METH +
C MITER). METH IS THE BASIC METHOD INDICATOR..
C METH = 1 MEANS THE ADAMS METHODS.
C METH = 2 MEANS THE BACKWARD DIFFERENTIATION FORMULAS
C (BDF), OR STIFF METHODS OF GEAR.
C MITER IS THE ITERATION METHOD INDICATOR..
C MITER = 0 MEANS FUNCTIONAL ITERATION (NO PARTIAL
C DERIVATIVES NEEDED).
C MITER = 1 MEANS THE CHORD METHOD WITH AN ANALYTIC
C JACOBIAN SUPPLIED IN THE USER DEFINED
C ROUTINE PDB. THIS METHOD IS IN GEARB,
C BUT IT SHOULD BE AVOIDED.
C MITER = 2 MEANS THE CHORD METHOD WITH THE JACO-
C BIAN CALCULATED IN PSETM. THIS IS THE
C ONLY METHOD WHERE PSETM IS USEFUL.
C MITER = 3 MEANS THE CHORD METHOD WITH THE JACO-
C BIAN REPLACED BY A DIAGONAL APPROXI-
C MATION BASED ON A DIRECTIONAL DERIV-
C ATIVE. THIS METHOD IS IN GEARB, BUT
C IS NOT GENERALLY USEFUL IN SOLVING PDE*S.
C SEE COMMENTS IN DRIVEP FOR ADDITIONAL INFORMATION ON THESE
C PARAMETERS. NATURALLY THESE PARAMETERS ARE DEPENDENT ON THE
C INTEGRATOR BEING USED. FINALLY THE USER SHOULD CALL THE
C INTEGRATOR.
C
C CALL DRIVEP (NODE,T0,H,U,TOUT,EPS,MF,INDEX,WORK,IWORK,X,Y)
C
C THE INTEGRATOR WILL RETURN THE SOLUTIONS (U) AT T = TOUT
C TO THE MAIN PROGRAM TO BE PRINTED OR PLOTTED. IF A CONTINUA-
C TION TO ANOTHER TOUT IS DESIRED, SIMPLY RESET TOUT AND CALL
C DRIVEP AGAIN.
C
C STOP
C END
C
C 2) THE BOUNDARY SUBROUTINES, BNDRYH AND BNDRYV PROVIDE
C PDETWO WITH THE COEFFICIENTS FOR THE HORIZONTAL AND
C VERTICAL BOUNDARY CONDITIONS RESPECTIVELY. THE CONSTRUCTION
C OF BNDRYV IS ANALOGOUS TO THE CONSTRUCTION OF BNDRYH. THE
C USER SHOULD CONSTRUCT BNDRYH AS FOLLOWS ..
C
C SUBROUTINE BNDRYH (T,X,Y,U,AH,BH,CH,NPDE)
C DIMENSION U(NPDE), AH(NPDE), BH(NPDE), CH(NPDE)
C
C THE INCOMING PARAMETERS X AND Y REPRESENT ANY POINT
C OF X(J) (J = 1,2,...,NX) AND EITHER Y(1) OR Y(NY) RESPECTIVELY.
C DEFINE THE FUNCTIONS AH(K), BH(K) AND CH(K) (K = 1,2,...,NPDE)
C FOR THE LOWER (Y=Y(1)) AND UPPER (Y=Y(NY)) BOUNDARIES. NULL
C BOUNDARY CONDITIONS ARE NOT ALLOWED AND IF SPECIFIED WILL
C RESULT IN A ZERO DIVIDE.
C
C TO INSURE A COMPLETELY ACCURATE DEFINITION OF THE BOUNDARY
C CONDITIONS AT THE CORNER MESH POINTS, THE USER SHOULD
C CONSULT THE APPROPRIATE DIFFERENCE APPROXIMATIONS (SEE
C SECTION 3 OF THE ACCOMPANYING PAPER).
C
C RETURN
C END
C
C 3) DIFFH AND DIFFV DEFINE FOR PDETWO THE HORIZONTAL AND VERTICAL
C DIFFUSION COEFFICIENTS RESPECTIVELY. THE CONSTRUCTION OF
C DIFFV IS ANALOGOUS TO THE CONSTRUCTION OF DIFFH. THE
C USER SHOULD CONSTRUCT DIFFH AS FOLLOWS ..
C
C SUBROUTINE DIFFH (T,X,Y,U,DH,NPDE)
C DIMENSION U(NPDE),DH(NPDE,NPDE)
C
C IN THIS ROUTINE DEFINE THE DH(L,K) COEFFICIENTS (L,K =
C 1,2,...,NPDE). THE INCOMING PARAMTERS X AND Y DENOTE EITHER
C A BOUNDARY POINT OR A MESH MIDPOINT.
C
C RETURN
C END
C
C 4) THE RIGHT HAND SIDE OF THE PDE IS DEFINED FOR PDETWO IN
C THE SUBROUTINE F. THIS ROUTINE IS CONSTRUCTED AS FOLLOWS..
C
C SUBROUTINE F (T,X,Y,U,UX,UY,DUXX,DUYY,DUDT,NPDE)
C DIMENSION U(NPDE), UX(NPDE), UY(NPDE), DUXX(NPDE,NPDE),
C * DUYY(NPDE,NPDE), DUDT(NPDE)
C
C IN THIS ROUTINE, THE INCOMING VALUES X AND Y REPRESENT
C THE MESH POINT BEING EVALUATED AND UX, UY, DUXX AND DUYY
C ARE THE VALUES DENOTED IN THE PROBLEM DEFINITION ABOVE.
C USING THESE VALUES, DEFINE IN DUDT(L) (L = 1,2,...,NPDE)
C THE RIGHT HAND SIDE OF THE PDE*S.
C
C RETURN
C END
C
C***********************************************************************
C THE INPUT PARAMETERS
C***********************************************************************
C
C NPDE IS THE NUMBER OF PARTIAL DIFFERENTIAL EQUATIONS.
C NX IS THE NUMBER OF MESH POINTS IN THE HORIZONTAL
C DIRECTION.
C NY IS THE NUMBER OF MESH POINTS IN THE VERTICAL
C DIRECTION.
C X ARE THE MESH POINTS IN THE HORIZONTAL DIRECTION.
C Y ARE THE MESH POINTS IN THE VERTICAL DIRECTION.
C T IS THE CURRENT TIME.
C U CONTAINS THE CURRENT SOLUTION VALUES FOR ALL THE
C MESH POINTS.
C DVS SAVES THE VERTICAL DIFFUSION COEFFICIENTS FOR
C FUTURE EVALUATIONS.
C DVST SAVES THE VERTICAL DIFFUSION COEFFICIENTS FOR THE
C FUTURE EVALUATIONS OF THE TOP ROW OF MESH POINTS.
C DH RETURNS FROM DIFFH THE HORIZONTAL DIFFUSION COEFFI-
C CIENTS AND CONTAINS DUXX, THE APPROXIMATIONS TO THE
C HORIZONTAL DIVERGENCE TERMS, ON CALLS TO F AND STORES
C THE PREVIOUS VALUE OF DUXX.
C DV RETURNS FROM DIFFV THE VERTICAL DIFFUSION COEFFI-
C CIENTS AND CONTAINS DUYY, THE APPROXIMATIONS TO THE
C VERTICAL DIVERGENCE TERMS, ON CALLS TO F.
C AH, BH, CH CONTAIN THE HORIZONTAL BOUNDARY COEFFICIENTS
C PASSED TO PDETWO FROM BNDRYH.
C AV, BV, CV CONTAIN THE VERTICAL BOUNDARY COEFFICIENTS
C PASSED TO PDETWO FROM BNDRYV.
C UX,UY STORE THE DIFFERENCE APPROXIMATIONS FOR THE FIRST
C PARTIAL OF U WITH RESPECT TO X AND Y RESPECTIVELY.
C DXI,DXIR,
C DXIC STORE INFORMATION ABOUT THE HORIZONTAL MESH SPACING.
C XAVG STORES INFORMATION ABOUT THE AVERAGE BETWEEN TWO MESH
C POINTS ON THE HORIZONTAL AXIS.
C UDVA,UDVB,
C UDHR,UDHL STORE INFORMATION FOR APPROXIMATING DUXX AND DUYY.
C UAVGH,UAVGV CONTAIN U AVERAGES FOR APPROXIMATING THE DIFFUSION
C COEFFICIENTS AT THE MESH MID-POINTS.
C
C***********************************************************************
C THE OUTPUT PARAMETERS
C***********************************************************************
C
C U CONTAINS SOLUTIONS UPDATED TO TIME T FOR THE
C BOUNDARY MESH POINTS DEFINED BY DIRICHLET BOUNDARY
C CONDITIONS.
C DUDT IS THE RIGHT HAND SIDE OF THE SYSTEM OF ODE*S PASSED
C TO THE INTEGRATOR. THESE VALUES ARE CALCULATED
C FROM THE CENTERED DIFFERENCE APPROXIMATIONS OF
C THE SPATIAL VARIABLES.
C
C***********************************************************************
C REFERENCES
C***********************************************************************
C
C 1. A. C. HINDMARSH, GEARB.. SOLUTION OF ORDINARY DIFFERENTIAL
C EQUATIONS HAVING BANDED JACOBIAN, UCID-30059 REV. 2,
C LAWRENCE LIVERMORE LABORATORY, P.O.BOX 808, LIVERMORE,
C CA 94550, JUNE 1977.
C
C 2. R.F. SINCOVEC AND N.K. MADSEN, SOFTWARE FOR NONLINEAR
C PARTIAL DIFFERENTIAL EQUATIONS, ACM TRANSACTIONS ON
C MATHEMATICAL SOFTWARE, SEPT. 1975, PP. 232-260.
C
C 3. R.F. SINCOVEC AND N.K. MADSEN, ALGORITHM 494 PDEONE,
C SOLUTIONS OF SYSTEMS OF PARTIAL DIFFERENTIAL EQUATIONS,
C ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE, SEPT.
C 1975, PP. 262-263.
C
C***********************************************************************
C
REAL T,U,X,Y,UX,UY,DXI,DYI,DXIC,DXIR,DYIA,DYIB,DYIC,
* UDHL,UDHR,DUDT,UDVA,UDVB,XAVG,YAVG,BH,BV,DH,DV,
* UAVGH,UAVGV,AH,AV,DVS,DVST,CH,CV,DX,DUXR,DUXL,DY,DUY
DIMENSION DX(4),DUXR(2),DUXL(2),DY(4),DUY(2,2)
DIMENSION DVS(NPDE,NPDE,NX), DVST(NPDE,NPDE,NX), DV(NPDE,NPDE),
* DH(NPDE,NPDE,2), U(NPDE,NX,NY), DUDT(NPDE,NX,NY), AH(NPDE),
* BH(NPDE), CH(NPDE), AV(NPDE), BV(NPDE), CV(NPDE), UX(NPDE),
* UY(NPDE), DXI(NX), DXIR(NX), DXIC(NX), XAVG(NX), UDVA(NPDE,NX),
* UDHR(NPDE), UDVB(NPDE), UDHL(NPDE), UAVGH(NPDE), UAVGV(NPDE),
* X(NX), Y(NY)
C
C DETERMINE THE MESH SPACING ALONG THE HORIZONTAL AXIS
C
DXI(1)=1./(X(2)-X(1))
DXIR(1)=DXI(1)
DXIC(1)=2.*DXI(1)
XAVG(1)=.5*(X(2)+X(1))
ILIM=NX-1
ILIMY=NY-1
DO 10 I=2,ILIM
XAVG(I)=.5*(X(I+1)+X(I))
DXI(I)=1./(X(I+1)-X(I-1))
DXIR(I)=1./(X(I+1)-X(I))
10 DXIC(I)=2.*DXI(I)
DXI(NX)=DXIR(ILIM)
DXIC(NX)=2.*DXI(NX)
C
C DETERMINE THE DIFFERENCING FOR THE BOUNDARIES
C
DX(1)=X(2)-X(1)
DX(2)=X(3)-X(2)
DX(3)=X(NX-1)-X(NX-2)
DX(4)=X(NX)-X(NX-1)
DUXL(1)=(DX(1)+DX(2))/(DX(1)*DX(2))
DUXL(2)=DX(1)/((DX(1)+DX(2))*DX(2))
DUXR(1)=(DX(4)+DX(3))/(DX(4)*DX(3))
DUXR(2)=DX(4)/((DX(4)+DX(3))*DX(3))
DY(1)=Y(2)-Y(1)
DY(2)=Y(3)-Y(2)
DY(3)=Y(NY-1)-Y(NY-2)
DY(4)=Y(NY)-Y(NY-1)
DUY(2,1)=(DY(1)+DY(2))/(DY(1)*DY(2))
DUY(2,2)=DY(1)/((DY(1)+DY(2))*DY(2))
DUY(1,1)=(DY(4)+DY(3))/(DY(4)*DY(3))
DUY(1,2)=DY(4)/((DY(4)+DY(3))*DY(3))
C
C THE FOLLOWING LOOP DETERMINES THE U APPROXIMATIONS FOR THE BOTTOM
C AND TOP BOUNDARIES. IC HAS THE VALUE OF 1 FOR THE TOP AND 2 FOR THE
C BOTTOM
C
DO 410 IC=1,2
IF (IC .EQ. 1) GO TO 20
M=1
MN=2
MNN=3
SIGN=1.0
GO TO 25
20 M=NY
MN=ILIMY
MNN=MN-1
SIGN=-1.0
25 DYI=1./(Y(MN)-Y(M))*SIGN
DYIA=DYI
DYIC=2.*DYI
YAVG=.5*(Y(M)+Y(MN))
C
C DETERMINE THE BOUNDARY CONDITIONS (LEFT CORNER)
C
CALL BNDRYH (T,X(2),Y(M),U(1,2,M),AH,BH,CH,NPDE)
CALL BNDRYV (T,X(1),Y(MN),U(1,1,MN),AV,BV,CV,NPDE)
ITESTL=0
ITESTH=0
DO 40 K=1,NPDE
IF (BH(K) .NE. 0.0) GO TO 30
U(K,2,M) = CH(K)/AH(K)
ITESTH=ITESTH+1
30 IF (BV(K) .NE. 0.0) GO TO 40
U(K,1,MN) = CV(K)/AV(K)
ITESTL=ITESTL+1
40 CONTINUE
CALL BNDRYH (T,X(1),Y(M),U(1,1,M),AH,BH,CH,NPDE)
CALL BNDRYV (T,X(1),Y(M),U(1,1,M),AV,BV,CV,NPDE)
ITEST=0
DO 70 K=1,NPDE
IF (BV(K) .NE. 0.0 .AND. BH(K) .NE. 0.0) GO TO 70
IF (BH(K) .EQ. 0.0 .AND. BV(K) .NE. 0.0) GO TO 50
IF (BH(K) .NE. 0.0 .AND. BV(K) .EQ. 0.0) GO TO 60
U(K,1,M)=(CH(K)/AH(K)+CV(K)/AV(K))*.5
ITEST=ITEST+1
GO TO 70
50 U(K,1,M)=CH(K)/AH(K)
GO TO 70
60 U(K,1,M)=CV(K)/AV(K)
70 CONTINUE
IF (ITEST .EQ. NPDE) GO TO 130
CALL BNDRYH (T,X(1),Y(M),U(1,1,M),AH,BH,CH,NPDE)
CALL BNDRYV (T,X(1),Y(M),U(1,1,M),AV,BV,CV,NPDE)
C
C EVALUATE THE DIFFUSION COEFFICIENTS (LEFT CORNER)
C
CALL DIFFH (T,X(1),Y(M),U(1,1,M),DH,NPDE)
CALL DIFFV (T,X(1),Y(M),U(1,1,M),DV,NPDE)
C
C EVALUATE DU/DX AND DU/DY (LEFT CORNER)
C
DO 120 K=1,NPDE
IF (BV(K) .NE. 0.0) GO TO 90
UX(K)=(U(K,2,M)-U(K,1,M))*DUXL(1)-(U(K,3,M)-U(K,1,M))*DUXL(2)
GO TO 100
90 UX(K)=(CV(K)-AV(K)*U(K,1,M))/BV(K)
100 IF (BH(K) .NE. 0.0) GO TO 110
UY(K)=((U(K,1,MN)-U(K,1,M))*DUY(IC,1)-(U(K,1,MNN)-U(K,1,M))
* *DUY(IC,2))*SIGN
GO TO 120
110 UY(K)=(CH(K)-AH(K)*U(K,1,M))/BH(K)
120 CONTINUE
C
C CALCULATE THE U AVERAGES (LEFT CORNER)
C
130 DO 140 K=1,NPDE
UAVGH(K)=.5*(U(K,2,M)+U(K,1,M))
UAVGV(K)=.5*(U(K,1,M)+U(K,1,MN))
UDHR(K)=(U(K,2,M)-U(K,1,M))*DXIR(1)
UDVA(K,1)=(U(K,1,MN)-U(K,1,M))*DYIA*SIGN
140 CONTINUE
C
C CALCULATE THE DIFFUSION COEFFICIENTS AT THE MIDPOINTS OF THE
C INTERVALS BETWEEN THE LEFT CORNER AND THE NEIGHBORING POINTS
C
CALL DIFFV (T,X(1),YAVG,UAVGV,DVS(1,1,1),NPDE)
CALL DIFFH (T,XAVG(1),Y(M),UAVGH,DH(1,1,2),NPDE)
IF (ITEST .EQ. NPDE) GO TO 160
C
C CALCULATE DUXX AND DUYY (LEFT CORNER)
C
DO 150 L=1,NPDE
DO 150 K=1,NPDE
DH(K,L,1)=DXIC(1)*(DH(K,L,2)*UDHR(L)-DH(K,L,1)
* *UX(L))
DV(K,L)=DYIC*(DVS(K,L,1)*UDVA(L,1)-DV(K,L)
* *UY(L))*SIGN
150 CONTINUE
C
C EVALUATE THE RIGHT SIDE OF PDE (LEFT CORNER)
C
CALL F (T,X(1),Y(M),U(1,1,M),UX,UY,DH,DV,DUDT(1,1,M),
* NPDE)
C
C SET DUDT = 0 FOR KNOWN BOUNDARY CONDITIONS
C
160 DO 170 K=1,NPDE
IF (BH(K) .EQ. 0.0 .OR. BV(K) .EQ. 0.0) DUDT(K,1,M)=0.0
170 CONTINUE
C
C EVALUATE THE RIGHT CORNER BOUNDARY CONDITIONS
C
CALL BNDRYH (T,X(NX),Y(M),U(1,NX,M),AH,BH,CH,NPDE)
CALL BNDRYV (T,X(NX),Y(M),U(1,NX,M),AV,BV,CV,NPDE)
ITEST = 0
DO 200 K=1,NPDE
IF (BV(K) .NE. 0.0 .AND. BH(K) .NE. 0.0) GO TO 200
IF (BH(K) .EQ. 0.0 .AND. BV(K) .NE. 0.0) GO TO 180
IF (BH(K) .NE. 0.0 .AND. BV(K) .EQ. 0.0) GO TO 190
U(K,NX,M)=(CV(K)/AV(K)+CH(K)/AH(K))*.5
ITEST=ITEST+1
GO TO 200
180 U(K,NX,M)=CH(K)/AH(K)
GO TO 200
190 U(K,NX,M)=CV(K)/AV(K)
200 CONTINUE
IBCK=1
IFWD=2
C
C LOOP ON THE HORIZONTAL BOUNDARY MESH POINTS FROM X(2) TO X(NX-1)
C
DO 300 I=2,ILIM
K=IBCK
IBCK=IFWD
IFWD=K
C
C EVALUATE THE HORIZONTAL BOUNDARY CONDITIONS
C
ITESTC=ITESTH
IF (I .EQ. ILIM) GO TO 220
CALL BNDRYH (T,X(I+1),Y(M),U(1,I+1,M),AH,BH,CH,
* NPDE)
ITESTH=0
DO 210 K=1,NPDE
IF (BH(K) .NE. 0.0) GO TO 210
U(K,I+1,M) = CH(K)/AH(K)
ITESTH=ITESTH+1
210 CONTINUE
220 CONTINUE
CALL BNDRYH (T,X(I),Y(M),U(1,I,M),AH,BH,CH,NPDE)
IF (ITESTC .EQ. NPDE) GO TO 260
C
C CALCULATE DU/DX AND DU/DY (HORIZONTAL BOUNDARY)
C
DO 250 K=1,NPDE
UX(K)=(U(K,I+1,M)-U(K,I-1,M))*DXI(I)
IF (BH(K) .NE. 0.0) GO TO 240
UY(K)=((U(K,I,MN)-U(K,I,M))*DUY(IC,1)-(U(K,I,MNN)-U(K,I,M))
* *DUY(IC,2))*SIGN
GO TO 250
240 UY(K)=(CH(K)-AH(K)*U(K,I,M))/BH(K)
250 CONTINUE
C
C DETERMINE U AVERAGE (HORIZONTAL BOUNDARY)
C
260 DO 270 K=1,NPDE
UAVGV(K)=(U(K,I,M)+U(K,I,MN))*.5
UAVGH(K)=(U(K,I+1,M)+U(K,I,M))*.5
UDHL(K)=UDHR(K)
UDHR(K)=(U(K,I+1,M)-U(K,I,M))*DXIR(I)
UDVA(K,I)=(U(K,I,MN)-U(K,I,M))*DYIA*SIGN
270 CONTINUE
C
C EVALUATE THE DIFFUSION COEFFICIENTS (HORIZONTAL BOUNDARY)
C
CALL DIFFV (T,X(I),YAVG,UAVGV,DVS(1,1,I),NPDE)
CALL DIFFH (T,XAVG(I),Y(M),UAVGH,DH(1,1,IFWD),NPDE)
IF (ITESTC .EQ. NPDE) GO TO 290
CALL DIFFV (T,X(I),Y(M),U(1,I,M),DV,NPDE)
C
C EVALUATE DUXX AND DUYY (HORIZONTAL BOUNDARY)
C
DO 280 L=1,NPDE
DO 280 K=1,NPDE
DH(K,L,IBCK)=DXIC(I)*(DH(K,L,IFWD)*UDHR(L)-
* DH(K,L,IBCK)*UDHL(L))
DV(K,L)=DYIC*(DVS(K,L,I)*UDVA(L,I)-DV(K,L)
* *UY(L))*SIGN
280 CONTINUE
C
C EVALUATE THE RIGHT SIDE OF THE PDE (HORIZONTAL BOUNDARY)
C
CALL F (T,X(I),Y(M),U(1,I,M),UX,UY,DH(1,1,IBCK),
* DV,DUDT(1,I,M),NPDE)
C
C SET DUDT = 0 FOR KNOWN BOUNDARY CONDITIONS
C
290 DO 300 K=1,NPDE
IF (BH(K) .EQ. 0.0) DUDT(K,I,M)=0.0
300 CONTINUE
C
C COMPLETE EVALUATING THE RIGHT CORNER
C
CALL BNDRYV (T,X(NX),Y(MN),U(1,NX,MN),AV,BV,CV,NPDE)
ITESTR=0
DO 305 K=1,NPDE
IF (BV(K) .NE. 0.0) GO TO 305
U(K,NX,MN) = CV(K)/AV(K)
ITESTR=ITESTR+1
305 CONTINUE
CALL BNDRYH (T,X(NX),Y(M),U(1,NX,M),AH,BH,CH,NPDE)
CALL BNDRYV (T,X(NX),Y(M),U(1,NX,M),AV,BV,CV,NPDE)
IF (ITEST .EQ. NPDE) GO TO 360
C
C EVALUATE THE DIFFUSION COEFFICIENTS (RIGHT CORNER)
C
CALL DIFFH (T,X(NX),Y(M),U(1,NX,M),DH(1,1,IBCK),NPDE)
CALL DIFFV (T,X(NX),Y(M),U(1,NX,M),DV,NPDE)
C
C EVALUATE DU/DX AND DU/DY (RIGHT CORNER)
C
DO 350 K=1,NPDE
IF (BV(K) .NE. 0.0) GO TO 320
UX(K)=(U(K,NX,M)-U(K,ILIM,M))*DUXR(1)-(U(K,NX,M)
* -U(K,ILIM-1,M))*DUXR(2)
GO TO 330
320 UX(K)=(CV(K)-AV(K)*U(K,NX,M))/BV(K)
330 IF (BH(K) .NE. 0.0) GO TO 340
UY(K)=((U(K,NX,MN)-U(K,NX,M))*DUY(IC,1)-(U(K,NX,MNN)
* -U(K,NX,M))*DUY(IC,2))*SIGN
GO TO 350
340 UY(K)=(CH(K)-AH(K)*U(K,NX,M))/BH(K)
350 CONTINUE
C
C EVALUATE THE VERTICAL U AVERAGE (RIGHT CORNER)
C
360 DO 370 K=1,NPDE
UAVGV(K)=(U(K,NX,M)+U(K,NX,MN))*.5
UDVA(K,NX)=(U(K,NX,MN)-U(K,NX,M))*DYI*SIGN
370 CONTINUE
C
C EVALUATE THE VERTICAL DIFFUSION COEFFICIENT (ABOVE THE RIGHT CORNER)
C
CALL DIFFV (T,X(NX),YAVG,UAVGV,DVS(1,1,NX),NPDE)
IF (ITEST .EQ. NPDE) GO TO 390
C
C EVALUATE DUXX AND DUYY (RIGHT CORNER)
C
DO 380 L=1,NPDE
DO 380 K=1,NPDE
DH(K,L,IBCK)=DXIC(NX)*(DH(K,L,IBCK)*UX(L)-
* DH(K,L,IFWD)*UDHR(L))
DV(K,L)=DYIC*(DVS(K,L,NX)*UDVA(L,NX)-DV(K,L)
* *UY(L))*SIGN
380 CONTINUE
C
C EVALUATE THE RIGHT SIDE OF THE PDE (RIGHT CORNER)
C
CALL F (T,X(NX),Y(M),U(1,NX,M),UX,UY,DH(1,1,IBCK),DV,
* DUDT(1,NX,M),NPDE)
390 DO 400 K=1,NPDE
IF (BH(K) .EQ. 0.0 .OR. BV(K) .EQ. 0.0) DUDT(K,NX,M)=0.0
400 CONTINUE
IF (IC .NE. 1) GO TO 410
DO 405 IIC=1,NX
DO 405 L=1,NPDE
DO 405 K=1,NPDE
405 DVST(K,L,IIC)=DVS(K,L,IIC)
410 CONTINUE
C
C DETERMINE THE U VALUES FOR THE J-TH ROW
C
IC=2
DO 780 J=2,ILIMY
IF (J .NE. ILIMY) GO TO 500
IC=1
500 DYI=1./(Y(J+1)-Y(J-1))
DYIC=2.*DYI
DYIB=DYIA
DYIA=1./(Y(J+1)-Y(J))
YAVG=(Y(J+1)+Y(J))*.5
C
C DETERMINE THE LEFT BOUNDARY (J-TH ROW)
C
ITESTC=ITESTL
IF (IC .EQ. 1) GO TO 510
CALL BNDRYV (T,X(1),Y(J+1),U(1,1,J+1),AV,BV,CV,NPDE)
ITESTL=0
DO 505 K=1,NPDE
IF (BV(K) .NE. 0.0) GO TO 505
U(K,1,J+1) = CV(K)/AV(K)
ITESTL=ITESTL+1
505 CONTINUE
510 CONTINUE
CALL BNDRYV (T,X(1),Y(J),U(1,1,J),AV,BV,CV,NPDE)
IF (ITESTC .EQ. NPDE) GO TO 560
C
C EVALUATE DU/DX AND DU/DY (LEFT BOUNDARY,J-TH ROW)
C
DO 550 K=1,NPDE
IF (BV(K) .NE. 0.0) GO TO 530
UX(K)=(U(K,2,J)-U(K,1,J))*DUXL(1)-(U(K,3,J)-U(K,1,J))*DUXL(2)
GO TO 540
530 UX(K)=(CV(K)-AV(K)*U(K,1,J))/BV(K)
540 UY(K)=(U(K,1,J+1)-U(K,1,J-1))*DYI
550 CONTINUE
C
C EVALUATE THE HORIZONTAL U AVERAGE (J-TH ROW)
C
560 DO 570 K=1,NPDE
UAVGH(K)=(U(K,2,J)+U(K,1,J))*.5
UDHR(K)=(U(K,2,J)-U(K,1,J))*DXIR(1)
UDVB(K)=UDVA(K,1)
UDVA(K,1)=(U(K,1,J+1)-U(K,1,J))*DYIA
570 CONTINUE
C
C EVALUATE THE DIFFUSION COEFFICIENTS (LEFT BOUNDARY,J-TH ROW)
C
DO 580 L=1,NPDE
DO 580 K=1,NPDE
DV(K,L)=DVS(K,L,1)
580 CONTINUE
IF (IC .EQ. 2) GO TO 582
DO 581 L=1,NPDE
DO 581 K=1,NPDE
581 DVS(K,L,1)=DVST(K,L,1)
582 CONTINUE
CALL DIFFH (T,X(1),Y(J),U(1,1,J),DH,NPDE)
CALL DIFFH (T,XAVG(1),Y(J),UAVGH,DH(1,1,2),NPDE)
IF (IC .EQ. 1) GO TO 590
C
C EVALUATE THE VERTICAL U AVERAGE (J-TH ROW)
C
DO 585 K=1,NPDE
UAVGV(K)=(U(K,1,J+1)+U(K,1,J))*.5
585 CONTINUE
CALL DIFFV (T,X(1),YAVG,UAVGV,DVS(1,1,1),NPDE)
C
C EVALUATE DUXX AND DUYY (LEFT BOUNDARY,J-TH ROW)
C
590 IF (ITESTC .EQ. NPDE) GO TO 610
DO 600 L=1,NPDE
DO 600 K=1,NPDE
DH(K,L,1)=DXIC(1)*(DH(K,L,2)*UDHR(L)-
* DH(K,L,1)*UX(L))
DV(K,L)=DYIC*(DVS(K,L,1)*UDVA(L,1)-DV(K,L)
* *UDVB(L))
600 CONTINUE
C
C EVALUATE THE RIGHT SIDE OF THE PDE (LEFT BOUNDARY,J-TH ROW)
C
CALL F (T,X(1),Y(J),U(1,1,J),UX,UY,DH,DV,DUDT(1,1,J),NPDE)
C
C SET DUDT = 0 FOR KNOWN LEFT BOUNDARY CONDITIONS
C
610 DO 620 K=1,NPDE
IF (BV(K) .EQ. 0.0) DUDT(K,1,J)=0.0
620 CONTINUE
C
C LOOP TO EVALUATE U IN THE CENTER OF THE GRID
C (2.LT.I.LT.NX-1, 2.LT.J.LT.NY-1)
C
IBCK=1
IFWD=2
DO 680 I=2,ILIM
K=IBCK
IBCK=IFWD
IFWD=K
C
C CALCULATE THE HORIZONTAL U AVERAGE, DU/DX AND DU/DY
C (I-TH POINT OF THE J-TH ROW)
C
DO 640 K=1,NPDE
UAVGH(K)=(U(K,I+1,J)+U(K,I,J))*.5
UX(K)=(U(K,I+1,J)-U(K,I-1,J))*DXI(I)
UY(K)=(U(K,I,J+1)-U(K,I,J-1))*DYI
640 CONTINUE
C
C EVALUATE THE DIFFUSION COEFFICIENTS (I-TH POINT OF THE J-TH ROW)
C
DO 650 L=1,NPDE
DO 650 K=1,NPDE
DV(K,L)=DVS(K,L,I)
650 CONTINUE
IF (IC .EQ. 2) GO TO 652
DO 651 L=1,NPDE
DO 651 K=1,NPDE
651 DVS(K,L,I)=DVST(K,L,I)
652 CONTINUE
CALL DIFFH (T,XAVG(I),Y(J),UAVGH,DH(1,1,IFWD),NPDE)
IF (IC .EQ. 1) GO TO 660
DO 655 K=1,NPDE
UAVGV(K)=(U(K,I,J+1)+U(K,I,J))*.5
655 CONTINUE
CALL DIFFV (T,X(I),YAVG,UAVGV,DVS(1,1,I),NPDE)
C
C EVALUATE DUXX AND DUYY (I-TH POINT OF THE J-TH ROW)
C
660 DO 670 L=1,NPDE
UDHL(L)=UDHR(L)
UDHR(L)=(U(L,I+1,J)-U(L,I,J))*DXIR(I)
UDVB(L)=UDVA(L,I)
UDVA(L,I)=(U(L,I,J+1)-U(L,I,J))*DYIA
DO 670 K=1,NPDE
DH(K,L,IBCK)=DXIC(I)*(DH(K,L,IFWD)*UDHR(L)-
* DH(K,L,IBCK)*UDHL(L))
DV(K,L)=DYIC*(DVS(K,L,I)*UDVA(L,I)-DV(K,L)
* *UDVB(L))
670 CONTINUE
C
C EVALUATE THE RIGHT SIDE OF THE PDE (I-TH POINT OF THE J-TH ROW)
C
CALL F (T,X(I),Y(J),U(1,I,J),UX,UY,DH(1,1,IBCK),DV,
* DUDT(1,I,J),NPDE)
680 CONTINUE
C
C EVALUATE THE RIGHT BOUNDARY FOR THE J-TH ROW
C
ITESTC=ITESTR
IF (IC .EQ. 1) GO TO 690
CALL BNDRYV (T,X(NX),Y(J+1),U(1,NX,J+1),AV,BV,CV,NPDE)
ITESTR=0
DO 685 K=1,NPDE
IF (BV(K) .NE. 0.0) GO TO 685
U(K,NX,J+1) = CV(K)/AV(K)
ITESTR=ITESTR+1
685 CONTINUE
690 CONTINUE
CALL BNDRYV (T,X(NX),Y(J),U(1,NX,J),AV,BV,CV,NPDE)
IF (ITESTC .EQ. NPDE) GO TO 730
C
C EVALUATE DU/DX AND DU/DY (RIGHT BOUNDARY,J-TH ROW)
C
DO 710 K=1,NPDE
UY(K)=(U(K,NX,J+1)-U(K,NX,J-1))*DYI
IF (BV(K) .NE. 0.0) GO TO 700
UX(K)=(U(K,NX,J)-U(K,ILIM,J))*DUXR(1)-(U(K,NX,J)
* -U(K,ILIM-1,J))*DUXR(2)
GO TO 710
700 UX(K)=(CV(K)-AV(K)*U(K,NX,J))/BV(K)
710 CONTINUE
C
C EVALUATE THE DIFFUSION COEFFICIENTS (RIGHT BOUNDARY,J-TH ROW)
C
DO 720 L=1,NPDE
DO 720 K=1,NPDE
DV(K,L)=DVS(K,L,NX)
720 CONTINUE
IF (IC .EQ. 2) GO TO 725
DO 721 L=1,NPDE
DO 721 K=1,NPDE
721 DVS(K,L,NX)=DVST(K,L,NX)
725 CONTINUE
CALL DIFFH (T,X(NX),Y(J),U(1,NX,J),DH(1,1,IBCK),NPDE)
730 DO 740 K=1,NPDE
UAVGV(K)=(U(K,NX,J+1)+U(K,NX,J))*.5
UDVB(K)=UDVA(K,NX)
UDVA(K,NX)=(U(K,NX,J+1)-U(K,NX,J))*DYIA
740 CONTINUE
IF (IC .EQ. 1) GO TO 750
CALL DIFFV (T,X(NX),YAVG,UAVGV,DVS(1,1,NX),NPDE)
750 IF (ITESTC .EQ. NPDE) GO TO 770
C
C EVALUATE DUXX AND DUYY (RIGHT BOUNDARY,J-TH ROW)
C
DO 760 L=1,NPDE
DO 760 K=1,NPDE
DH(K,L,IBCK)=DXIC(NX)*(DH(K,L,IBCK)*UX(L)-
* DH(K,L,IFWD)*UDHR(L))
DV(K,L)=DYIC*(DVS(K,L,NX)*UDVA(L,NX)-DV(K,L)
* *UDVB(L))
760 CONTINUE
C
C EVALUATE THE RIGHT SIDE OF THE PDE (RIGHT BOUNDARY,J-TH ROW)
C
CALL F (T,X(NX),Y(J),U(1,NX,J),UX,UY,DH(1,1,IBCK),DV,
* DUDT(1,NX,J),NPDE)
C
C SET DUDT = 0 FOR KNOWN BOUNDARY VALUES
C
770 DO 780 K=1,NPDE
IF (BV(K) .EQ. 0.0) DUDT(K,NX,J)=0.0
780 CONTINUE
RETURN
END
SUBROUTINE STRSET (N,MF,INDEX,IWORK,LOUT,ML,MU) STRSET
C
C***********************************************************************
C
C STRSET IS DESIGNED TO EXTRACT AND CHECK THE INPUT INTEGER
C PARAMETERS FROM THE FIRST SIX LOCATIONS OF THE ARRAY IWORK AND TO
C PROVIDE DYNAMIC DIMENSIONING FOR THE ARRAYS IN THE ROUTINE PDETWO
C AND THE MODIFIED GEARB PACKAGE. STRSET ESTABLISHES THE DIMENSION
C POINTERS, IW1 - IW27 FOR WORK, THE USER DIMENSIONED WORKING STOR-
C AGE ARRAY. THE CORRESPONDENCE BETWEEN THIS WORKING STORAGE AND THE
C ARRAYS IN PDETWO AND THE MODIFIED GEARB IS LISTED BELOW. STRSET
C ALSO SETS THE LOWER AND THE UPPER BANDWIDTHS ML AND MU FOR THE
C JACOBIAN MATRIX.
C
C PARAMETERS ARE SELF-EXPLANATORY IN THE COMMENT STATEMENTS
C BELOW.
C
C***********************************************************************
C CORRESPONDENCE BETWEEN WORK AND THE ARRAYS IN PDETWO AND THE MODI-
C FIED GEARB PACKAGE
C***********************************************************************
C
C CORRESPONDING ARRAY
C LOCATION IN WORK IN PDETWO LENGTH
C ------------------ --------------------- ----------------------
C
C WORK(1) DVS NPDE * NPDE * NX
C WORK(IW1) DVST NPDE * NPDE * NX
C WORK(IW2) DV NPDE * NPDE
C WORK(IW3) DH NPDE * NPDE * 2
C WORK(IW4) AH NPDE
C WORK(IW5) BH NPDE
C WORK(IW6) CH NPDE
C WORK(IW7) AV NPDE
C WORK(IW8) BV NPDE
C WORK(IW9) CV NPDE
C WORK(IW10) UX NPDE
C WORK(IW11) UY NPDE
C WORK(IW12) DXI NX
C WORK(IW13) DXIR NX
C WORK(IW14) DXIC NX
C WORK(IW15) XAVG NX
C WORK(IW16) UDVA NPDE * NX
C WORK(IW17) UDHR NPDE
C WORK(IW18) UDVB NPDE
C WORK(IW19) UDHL NPDE
C WORK(IW20) UAVGH NPDE
C WORK(IW21) UAVGV NPDE
C
C CORRESPONDING ARRAY
C LOCATION IN WORK IN THE MODIFIED GEARB LENGTH
C PACKAGE
C ------------------ --------------------- ----------------------
C
C WORK(IW22) U NODE*(MORDER+1)
C WORK(IW23) UMAX NODE
C WORK(IW24) ERROR NODE
C WORK(IW25) SAVE1 NODE
C WORK(IW26) SAVE2 NODE
C WORK(IW27) PW LPW
C
C WHERE NODE = NPDE*NX*NY, NPDE IS THE NUMBER OF PARTIAL
C DIFFERENTIAL EQUATIONS, NX IS THE NUMBER OF MESH POINTS IN
C THE HORIZONTAL DIRECTION AND NY IS THE NUMBER OF MESH POINTS
C IN THE VERTICAL DIRECTION. SEE COMMENTS IN PDETWO FOR THE
C DEFINITION OF MORDER AND LPW.
C
C***********************************************************************
C
COMMON /IPARAM/ NPDE,NX,NY,MORDER,LWORK,LIWORK
COMMON /ISTORE/ IW1,IW2,IW3,IW4,IW5,IW6,IW7,IW8,IW9,IW10,IW11,
* IW12,IW13,IW14,IW15,IW16,IW17,IW18,IW19,IW20,IW21,
* IW22,IW23,IW24,IW25,IW26,IW27
DIMENSION IWORK(N)
C
C EXTRACT THE INPUT INTEGER PARAMETERS FROM THE ARRAY IWORK
C
NPDE = IWORK(1)
NX = IWORK(2)
NY = IWORK(3)
MORDER = IWORK(4)
NRWK = IWORK(5)
NRIWK = IWORK(6)
C
C CHECK IF NUMBER OF ODE, NX, AND NY ARE CORRECT
C
NODE = NPDE*NX*NY
IF (N .NE. NODE) GO TO 110
IF (NX .LT. 3) GO TO 120
IF (NY .LT. 3) GO TO 130
C
C CHECK IF METH AND MITER OBTAINED FROM THE METHOD FLAG MF ARE WITHIN
C THE LIMIT OF THE ALLOWED VALUES
C
METH = MF/10
IF (METH .LT. 1 .OR. METH .GT. 2) GO TO 140
MITER = MF - 10*METH
IF (MITER .LT. 0 .OR. MITER .GT. 3) GO TO 150
C
C CHECK IF THE MAXIMUM ORDER OF THE METHOD IS WITHIN THE LIMIT
C IF IT EXCEEDS THE UPPER LIMIT, THE MAXIMUM ORDER AVALIABLE IS USED
C
IF (MORDER .LT. 1) GO TO 160
IF (METH .EQ. 1 .AND. MORDER .GT. 12) MORDER = 12
IF (METH .EQ. 2 .AND. MORDER .GT. 5) MORDER = 5
C
C CHECK IF THE LENGTHS OF THE ARRAYS WORK AND IWORK ARE LARGE ENOUGH
C
LU = NODE*(MORDER+1)
LPW = 1
IF(MITER .EQ. 1 .OR. MITER .EQ. 2) LPW = (3*(NX+1)*NPDE - 2)*NODE
IF(MITER .EQ. 3) LPW = NODE
LWORK = NPDE*(NPDE*(2*NX+3) + 13 + NX) + 4*(NX + NODE) + LU + LPW
IF(NRWK .LT. LWORK) GO TO 170
LIWORK = NODE
IF(NRIWK .LT. LIWORK) GO TO 180
GO TO 200
C
C PRINT ERROR MESSAGE
C
110 WRITE (LOUT,115) N,NODE
115 FORMAT(//16H ILLEGAL INPUT..,7H NODE =,I9,2X,28HIS NOT EQUAL TO NP
1DE.NX.NY =,I9//)
INDEX = -4
GO TO 300
C
120 WRITE (LOUT,125) NX
125 FORMAT(//16H ILLEGAL INPUT..,5H NX =,I5,2X,14HIS LESS THAN 3//)
INDEX = -4
GO TO 300
C
130 WRITE (LOUT,135) NY
135 FORMAT(//16H ILLEGAL INPUT..,5H NY =,I5,2X,14HIS LESS THAN 3//)
INDEX = -4
GO TO 300
C
140 WRITE (LOUT,145) METH
145 FORMAT(//16H ILLEGAL INPUT..,7H METH =,I3,2X,14HIS NOT ALLOWED//)
INDEX = -4
GO TO 300
C
150 WRITE (LOUT,155) MITER
155 FORMAT(//16H ILLEGAL INPUT..,8H MITER =,I3,2X,14HIS NOT ALLOWED//)
INDEX = -4
GO TO 300
C
160 WRITE (LOUT,165) MORDER
165 FORMAT(//16H ILLEGAL INPUT..,41H MAXIMUM ORDER OF THE METHODS AVAI
1LABLE =,I3,2X,14HIS NOT ALLOWED//)
INDEX = -4
GO TO 300
C
170 WRITE(LOUT,175) NRWK,LWORK
175 FORMAT(//23H INSUFFICIENT STORAGE..,11H THE LENGTH,I9,2X,60HOF THE
1 WORK ARRAY IS LESS THAN THE REQUIRED MINIMUM STORAGE ,I9//)
INDEX = -6
GO TO 300
C
180 WRITE(LOUT,185) NRIWK,LIWORK
185 FORMAT(//23H INSUFFICIENT STORAGE..,11H THE LENGTH,I9,2X,60HOF THE
1 IWORK ARRAY IS LESS THAN THE REQUIRED MINIMUM STORAGE,I9//)
INDEX = -6
GO TO 300
C
C ASSIGN DIMENSIONS FOR THE ARRAYS
C
200 CONTINUE
NPDE2 = NPDE*NPDE
IW1 = NPDE2*NX + 1
IW2 = IW1 + NPDE2*NX
IW3 = IW2 + NPDE2
IW4 = IW3 + 2*NPDE2
IW5 = IW4 + NPDE
IW6 = IW5 + NPDE
IW7 = IW6 + NPDE
IW8 = IW7 + NPDE
IW9 = IW8 + NPDE
IW10 = IW9 + NPDE
IW11 = IW10 + NPDE
IW12 = IW11 + NPDE
IW13 = IW12 + NX
IW14 = IW13 + NX
IW15 = IW14 + NX
IW16 = IW15 + NX
IW17 = IW16 + NPDE * NX
IW18 = IW17 + NPDE
IW19 = IW18 + NPDE
IW20 = IW19 + NPDE
IW21 = IW20 + NPDE
IW22 = IW21 + NPDE
IW23 = IW22 + NODE*(MORDER+1)
IW24 = IW23 + NODE
IW25 = IW24 + NODE
IW26 = IW25 + NODE
IW27 = IW26 + NODE
C
C SET BANDWIDTH MU AND ML FOR THE JACOBIAN
C
MU = NPDE*(NX+1) - 1
ML = MU
300 RETURN
END
SUBROUTINE DRIVEP (N,T0,H0,U0,TOUT,EPS,MF,INDEX,WORK,IWORK,XM,YM) DRIVEP
C THIS IS A MODIFIED VERSION OF GEARB(1),
C A PACKAGE FOR THE SOLUTION OF THE INITIAL VALUE
C PROBLEM FOR SYSTEMS OF ORDINARY DIFFERENTIAL EQUATIONS,
C DU/DT = F(U,T), U = (U(1),U(2),...,U(N)).
C GEARB IS A VARIANT OF THE GEAR PACKAGE TO BE USED WHEN
C THE JACOBIAN MATRIX DF/DU HAS BANDED OR NEARLY BANDED FORM.
C SUBROUTINE DRIVEP IS A DRIVER ROUTINE FOR THE MODIFIED
C GEARB PACKAGE. THE PACKAGE IS MODIFIED EXCLUSIVELY FOR USE
C WITH PDETWO AND PSETM.
C
C THIS PACKAGE WAS CONSTRUCTED SO AS TO CONFORM TO AS MANY
C ANSI-FORTRAN RULES AS WAS CONVENIENTLY POSSIBLE. THE FORTRAN
C USED VIOLATES ONLY ONE ANSI STANDARD - AN ARRAY NAME APPEARS
C IN A DATA STATEMENT IN THE SUBROUTINE COSET.
C
C
C REFERENCES
C 1. A. C. HINDMARSH, GEARB.. SOLUTION OF ORDINARY DIFFERENTIAL
C EQUATIONS HAVING BANDED JACOBIAN, UCID-30059 REV. 2,
C LAWRENCE LIVERMORE LABORATORY, P.O.BOX 808, LIVERMORE,
C CA 94550, JUNE 1977.
C
C-----------------------------------------------------------------------
C DRIVEP IS TO BE CALLED ONCE FOR EACH OUTPUT VALUE OF T, AND
C IN TURN MAKES REPEATED CALLS TO THE CORE INTEGRATOR, STIFFP.
C
C THE INPUT PARAMETERS ARE..
C N = THE NUMBER OF FIRST-ORDER DIFFERENTIAL EQUATIONS.
C IT SHOULD BE SET TO NODE = NPDE*NX*NY. SEE IWORK
C BELOW FOR THE DEFINITION OF NPDE, NX, NY.
C T0 = THE INITIAL VALUE OF T, THE INDEPENDENT VARIABLE
C (USED ONLY ON FIRST CALL).
C H0 = THE NEXT STEP SIZE IN T (USED FOR INPUT ONLY ON THE
C FIRST CALL).
C U0 = A VECTOR OF LENGTH N CONTAINING THE INITIAL VALUES OF
C U (USED FOR INPUT ONLY ON FIRST CALL).
C TOUT = THE VALUE OF T AT WHICH OUTPUT IS DESIRED NEXT.
C INTEGRATION WILL NORMALLY GO SLIGHTLY BEYOND TOUT
C AND THE PACKAGE WILL INTERPOLATE TO T = TOUT.
C EPS = THE RELATIVE ERROR BOUND (USED ONLY ON THE
C FIRST CALL, UNLESS INDEX = -1). SINGLE STEP ERROR
C ESTIMATES DIVIDED BY UMAX(I) WILL BE KEPT LESS THAN
C EPS IN ROOT-MEAN-SQUARE NORM (I.E. EUCLIDEAN NORM
C DIVIDED BY SQRT(N) ). THE VECTOR UMAX OF WEIGHTS
C IS COMPUTED IN DRIVEP. INITIALLY UMAX(I) IS
C ABS(U(I)), WITH A DEFAULT VALUE OF 1 IF Y(I) = 0
C INITIALLY. THEREAFTER, UMAX(I) IS THE LARGEST VALUE
C OF ABS(U(I)) SEEN SO FAR, OR THE INITIAL UMAX(I) IF
C THAT IS LARGER. TO ALTER EITHER OF THESE, CHANGE THE
C APPROPRIATE STATEMENTS IN THE DO-LOOPS ENDING AT
C STATEMENTS 10 AND 70 BELOW.
C MF = THE METHOD FLAG (USED ONLY ON FIRST CALL, UNLESS
C INDEX = -1). ALLOWED VALUES ARE 10, 11, 12, 13,
C 20, 21, 22, 23. MF HAS TWO DECIMAL DIGITS, METH
C AND MITER (MF = 10*METH + MITER).
C METH IS THE BASIC METHOD INDICATOR..
C METH = 1 MEANS THE ADAMS METHODS.
C METH = 2 MEANS THE BACKWARD DIFFERENTIATION
C FORMULAS (BDF), OR STIFF METHODS OF GEAR.
C MITER IS THE ITERATION METHOD INDICATOR..
C MITER = 0 MEANS FUNCTIONAL ITERATION (NO PARTIAL
C DERIVATIVES NEEDED).
C MITER = 1 MEANS CHORD METHOD WITH ANALYTIC JACOBIAN.
C FOR THIS USER SUPPLIES SUBROUTINE
C PDB (SEE DESCRIPTION BELOW).
C MITER = 2 MEANS CHORD METHOD WITH JACOBIAN CALCULATED
C INTERNALLY BY FINITE DIFFERENCES.
C MITER = 3 MEANS CHORD METHOD WITH JACOBIAN REPLACED
C BY A DIAGONAL APPROXIMATION BASED ON A
C DIRECTIONAL DERIVATIVE.
C INDEX = INTEGER USED ON INPUT TO INDICATE TYPE OF CALL,
C WITH THE FOLLOWING VALUES AND MEANINGS..
C 1 THIS IS THE FIRST CALL FOR THIS PROBLEM.
C 0 THIS IS NOT THE FIRST CALL FOR THIS PROBLEM,
C AND INTEGRATION IS TO CONTINUE.
C -1 THIS IS NOT THE FIRST CALL FOR THE PROBLEM,
C AND THE USER HAS RESET N, EPS, AND/OR MF.
C 2 SAME AS 0 EXCEPT THAT TOUT IS TO BE HIT
C EXACTLY (NO INTERPOLATION IS DONE).
C ASSUMES TOUT .GE. THE CURRENT T.
C 3 SAME AS 0 EXCEPT CONTROL RETURNS TO CALLING
C PROGRAM AFTER ONE STEP. TOUT IS IGNORED.
C SINCE THE NORMAL OUTPUT VALUE OF INDEX IS 0,
C IT NEED NOT BE RESET FOR NORMAL CONTINUATION.
C WORK = FLOATING POINT WORKING ARRAY. SEE STRSET FOR ADDITIONAL
C DETAILS.
C IWORK = INTEGER WORKING ARRAY. THE FIRST SIX LOCATIONS MUST BE
C DEFINED AS FOLLOWS..
C IWORK(1) = NPDE
C IWORK(2) = NX
C IWORK(3) = NY
C IWORK(4) = MORDER
C IWORK(5) = NRWK
C IWORK(6) = NRIWK
C WHERE NPDE IS THE NUMBER OF PARTIAL DIFFERENTIAL EQUA-
C TIONS, NX IS THE NUMBER OF MESH POINTS IN THE HORIZON-
C TAL DIRECTION AND NY IS THE NUMBER OF MESH POINTS IN
C THE VERTICAL DIRECTION, MORDER IS THE MAXIMUM ORDER OF
C METHOD TO BE USED IN THE MODIFIED GEARB, NRWK IS THE
C LENGTH OF THE ARRAY WORK AND NRIWK IS THE LENGTH OF THE
C ARRAY IWORK. MORDER MUST BE LESS THAN OR EQUAL TO 12
C IF METH = 1 OR IF METH = 2, IT MUST BE LESS THAN OR
C EQUAL TO 5.
C XM = ARE THE MESH POINTS IN THE HORIZONTAL DIRECTION.
C YM = ARE THE MESH POINTS IN THE VERTICAL DIRECTION.
C
C AFTER THE INITIAL CALL, IF A NORMAL RETURN OCCURRED AND A NORMAL
C CONTINUATION IS DESIRED, SIMPLY RESET TOUT AND CALL AGAIN.
C ALL OTHER PARAMETERS WILL BE READY FOR THE NEXT CALL.
C A CHANGE OF PARAMETERS WITH INDEX = -1 CAN BE MADE AFTER
C EITHER A SUCCESSFUL OR AN UNSUCCESSFUL RETURN.
C
C THE OUTPUT PARAMETERS ARE..
C H0 = THE STEP SIZE H USED LAST, WHETHER SUCCESSFULLY OR NOT.
C U0 = THE COMPUTED VALUES OF U AT T = TOUT.
C TOUT = THE OUTPUT VALUE OF T. IF INTEGRATION WAS SUCCESSFUL,
C AND THE INPUT VALUE OF INDEX WAS NOT 3, TOUT IS
C UNCHANGED FROM ITS INPUT VALUE. OTHERWISE, TOUT
C IS THE CURRENT VALUE OF T TO WHICH INTEGRATION
C HAS BEEN COMPLETED.
C INDEX = INTEGER USED ON OUTPUT TO INDICATE RESULTS,
C WITH THE FOLLOWING VALUES AND MEANINGS..
C 0 INTEGRATION WAS COMPLETED TO TOUT OR BEYOND.
C -1 THE INTEGRATION WAS HALTED AFTER FAILING TO PASS THE
C ERROR TEST EVEN AFTER REDUCING H BY A FACTOR OF
C 1.E10 FROM ITS INITIAL VALUE.
C -2 AFTER SOME INITIAL SUCCESS, THE INTEGRATION WAS
C HALTED EITHER BY REPEATED ERROR TEST FAILURES OR BY
C A TEST ON EPS. TOO MUCH ACCURACY HAS BEEN REQUESTED.
C -3 THE INTEGRATION WAS HALTED AFTER FAILING TO ACHIEVE
C CORRECTOR CONVERGENCE EVEN AFTER REDUCING H BY A
C FACTOR OF 1.E10 FROM ITS INITIAL VALUE.
C -4 IMMEDIATE HALT BECAUSE OF ILLEGAL VALUES OF INPUT
C PARAMETERS. SEE PRINTED MESSAGE.
C -5 INDEX WAS -1 ON INPUT, BUT THE DESIRED CHANGES OF
C PARAMETERS WERE NOT IMPLEMENTED BECAUSE TOUT
C WAS NOT BEYOND T. INTERPOLATION TO T = TOUT WAS
C PERFORMED AS ON A NORMAL RETURN. TO TRY AGAIN,
C SIMPLY CALL AGAIN WITH INDEX = -1 AND A NEW TOUT.
C -6 INSUFFICIENT STORAGE. THE LENGTH OF THE ARRAY WORK
C OR IWORK IS LESS THAN THE REQUIRED MINIMUM STORAGE.
C
C IN ADDITION TO DRIVEP, THE FOLLOWING ROUTINES ARE PROVIDED IN
C THE PACKAGE..
C INTERP(TOUT,U,N,U0,WORK,XM,YM) INTERPOLATES TO GET THE OUTPUT
C VALUES AT T = TOUT, FROM THE DATA IN THE U
C ARRAY.
C STIFFP(U,UMAX,..) IS THE CORE INTEGRATOR ROUTINE. IT PERFORMS A
C SINGLE STEP AND ASSOCIATED ERROR CONTROL.
C COSET(METH,NQ,EL,TQ,MAXDER) SETS COEFFICIENTS FOR USE IN
C THE CORE INTEGRATOR.
C DECBR(M0,N,ML,MU,B,IP,IER) AND SOLBR(M0,N,ML,MU,B,IP,X)
C ARE USED TO SOLVE BANDED LINEAR SYSTEMS,
C WITH A BAND MATRIX STORED BY ROWS IN B.
C NOTE.. DECBR, AND SOLBR ARE CALLED ONLY IF MITER = 1 OR 2.
C
C THE COMMON BLOCK GEAR3 CAN BE ACCESSED EXTERNALLY BY THE USER
C IF DESIRED. IT CONTAINS THE STEP SIZE LAST USED (SUCCESSFULLY),
C THE ORDER LAST USED (SUCCESSFULLY), THE NUMBER OF STEPS TAKEN
C SO FAR, THE NUMBER OF F EVALUATIONS SO FAR, AND THE NUMBER OF
C JACOBIAN EVALUATIONS SO FAR.
C
C IN THE DATA STATEMENT BELOW, SET..
C LOUT = THE LOGICAL UNIT NUMBER FOR THE OUTPUT OF MESSAGES
C DURING THE INTEGRATION.
C-----------------------------------------------------------------------
INTEGER N,MF,INDEX,ML,MU
INTEGER NC,MFC,KFLAG,JSTART,MLC,MUC,MW,NEBAND,
1 NEB1,NNEB,NQUSED,NSTEP,NFE,NJE
INTEGER LOUT,I,NHCUT,KGO
INTEGER IWORK
REAL T0,H0,U0,TOUT,EPS
REAL T,H,HMIN,HMAX,EPSC,UROUND,EPSJ,HUSED
REAL TOUTP,AYI,D
REAL WORK,XM,YM
COMMON /GEAR1/ T,H,HMIN,HMAX,EPSC,UROUND,NC,MFC,KFLAG,JSTART
COMMON /GEAR2/ EPSJ,MLC,MUC,MW,NEBAND,NEB1,NNEB
COMMON /GEAR3/ HUSED,NQUSED,NSTEP,NFE,NJE
COMMON /IPARAM/ NPDE,NX,NY,MORDER,LWORK,LIWORK
COMMON /ISTORE/ IW1,IW2,IW3,IW4,IW5,IW6,IW7,IW8,IW9,IW10,IW11,
* IW12,IW13,IW14,IW15,IW16,IW17,IW18,IW19,IW20,IW21,
* IW22,IW23,IW24,IW25,IW26,IW27
DIMENSION U0(N),WORK(1),IWORK(1),XM(1),YM(1)
DATA LOUT/6/
C
IF (INDEX .EQ. 0) GO TO 20
IF (INDEX .EQ. 2) GO TO 25
IF (INDEX .EQ. -1) GO TO 30
IF (INDEX .EQ. 3) GO TO 40
IF (INDEX .NE. 1) GO TO 430
IF (EPS .LE. 0.0) GO TO 400
IF (N .LE. 0) GO TO 410
IF ((T0-TOUT)*H0 .GE. 0.0) GO TO 420
C-----------------------------------------------------------------------
C CALL STRSET TO PROVIDE DYNAMIC DIMENSIONING FOR THE ARRAYS IN THE
C ROUTINE PDETWO AND THE MODIFIED GEARB PACKAGE.
C-----------------------------------------------------------------------
CALL STRSET (N,MF,INDEX,IWORK,LOUT,ML,MU)
IF (INDEX .EQ. -4 .OR. INDEX .EQ. -6) RETURN
C-----------------------------------------------------------------------
C IF INITIAL VALUES OF UMAX OTHER THAN THOSE SET BELOW ARE DESIRED,
C THEY SHOULD BE SET HERE. ALL UMAX(I) MUST BE POSITIVE.
C IF VALUES FOR HMIN OR HMAX, THE BOUNDS ON ABS(H), OTHER THAN
C THOSE BELOW ARE DESIRED, THEY SHOULD BE SET BELOW.
C IN THE FOLLOWING STATEMENT, SET
C UROUND = THE UNIT ROUNDOFF OF THE MACHINE, I.E. THE SMALLEST
C POSITIVE UR SUCH THAT 1. + UR .NE. 1. ON THE MACHINE.
C (.711E-14 IS THE UNIT ROUNDOFF FOR CDC 6000/7000 SERIES
C MACHINE).
C-----------------------------------------------------------------------
C
UROUND = 7.11E-14
DO 10 I = 1,N
IM1 = I - 1
IWA = IW22 + IM1
IWB = IW23 + IM1
WORK(IWB) = ABS(U0(I))
IF(WORK(IWB) .EQ. 0.0) WORK(IWB) = 1.0
10 WORK(IWA) = U0(I)
NC = N
T = T0
H = H0
IF ((T+H) .EQ. T) WRITE(LOUT,15)
15 FORMAT(35H WARNING.. T + H = T ON NEXT STEP.)
HMIN = ABS(H0)
HMAX = ABS(T0-TOUT)*10.0
EPSC = EPS
MFC = MF
JSTART = 0
EPSJ = SQRT(UROUND)
MLC = ML
MUC = MU
MW = ML + MU + 1
NEBAND = 2*ML + MU + 1
NEB1 = NEBAND - 1
NNEB = N*NEBAND
NHCUT = 0
GO TO 50
C
C TOUTP IS THE PREVIOUS VALUE OF TOUT FOR USE IN HMAX.-----------------
20 HMAX = ABS(TOUT-TOUTP)*10.0
GO TO 80
C
25 HMAX = ABS(TOUT-TOUTP)*10.0
IF ((T-TOUT)*H .GE. 0.0) GO TO 500
GO TO 85
C
30 IF ((T-TOUT)*H .GE. 0.0) GO TO 440
JSTART = -1
NC = N
EPSC = EPS
MFC = MF
C
40 IF ((T+H) .EQ. T) WRITE(LOUT,15)
C
50 CALL STIFFP (WORK(IW22),WORK(IW23),WORK(IW24),WORK(IW25),
* WORK(IW26),WORK(IW27),N,WORK,IWORK,XM,YM)
C
KGO = 1 - KFLAG
GO TO (60, 100, 200, 300), KGO
C KFLAG = 0, -1, -2, -3
C
60 CONTINUE
C-----------------------------------------------------------------------
C NORMAL RETURN FROM INTEGRATOR.
C
C THE WEIGHTS UMAX(I) ARE UPDATED. IF DIFFERENT VALUES ARE DESIRED,
C THEY SHOULD BE SET HERE. A TEST IS MADE FOR EPS BEING TOO SMALL
C FOR THE MACHINE PRECISION.
C
C ANY OTHER TESTS OR CALCULATIONS THAT ARE REQUIRED AFTER EVERY
C STEP SHOULD BE INSERTED HERE.
C
C IF INDEX = 3, U0 IS SET TO THE CURRENT U VALUES ON RETURN.
C IF INDEX = 2, H IS CONTROLLED TO HIT TOUT (WITHIN ROUNDOFF
C ERROR), AND THEN THE CURRENT U VALUES ARE PUT IN U0 ON RETURN.
C FOR ANY OTHER VALUE OF INDEX, CONTROL RETURNS TO THE INTEGRATOR
C UNLESS TOUT HAS BEEN REACHED. THEN INTERPOLATED VALUES OF U ARE
C COMPUTED AND STORED IN U0 ON RETURN.
C IF INTERPOLATION IS NOT DESIRED, THE CALL TO INTERP SHOULD BE
C REMOVED AND CONTROL TRANSFERRED TO STATEMENT 500 INSTEAD OF 520.
C-----------------------------------------------------------------------
D = 0.0
DO 70 I = 1,N
IM1 = I - 1
IWA = IW22 + IM1
IWB = IW23 + IM1
AYI = ABS(WORK(IWA))
WORK(IWB) = AMAX1(WORK(IWB), AYI)
70 D = D + (AYI/WORK(IWB))**2
D = D*(UROUND/EPS)**2
IF (D .GT. FLOAT(N)) GO TO 250
IF (INDEX .EQ. 3) GO TO 500
IF (INDEX .EQ. 2) GO TO 85
80 IF ((T-TOUT)*H .LT. 0.0) GO TO 40
CALL INTERP (TOUT,WORK(IW22),N,U0,WORK,XM,YM)
GO TO 520
85 IF (((T+H)-TOUT)*H .LE. 0.0) GO TO 40
IF (ABS(T-TOUT) .LE. 100.0*UROUND*HMAX) GO TO 500
IF ((T-TOUT)*H .GE. 0.0) GO TO 500
H = (TOUT - T)*(1.0 - 4.0*UROUND)
JSTART = -1
GO TO 40
C-----------------------------------------------------------------------
C ON AN ERROR RETURN FROM INTEGRATOR, AN IMMEDIATE RETURN OCCURS IF
C KFLAG = -2, AND RECOVERY ATTEMPTS ARE MADE OTHERWISE.
C TO RECOVER, H AND HMIN ARE REDUCED BY A FACTOR OF .1 UP TO 10
C TIMES BEFORE GIVING UP.
C-----------------------------------------------------------------------
100 WRITE (LOUT,105) T
105 FORMAT(//35H KFLAG = -1 FROM INTEGRATOR AT T = ,E16.8/
1 38H ERROR TEST FAILED WITH ABS(H) = HMIN/)
110 IF (NHCUT .EQ. 10) GO TO 150
NHCUT = NHCUT + 1
HMIN = .10*HMIN
H = .10*H
WRITE (LOUT,115) H
115 FORMAT(24H H HAS BEEN REDUCED TO ,E16.8,
1 26H AND STEP WILL BE RETRIED//)
JSTART = -1
GO TO 40
C
150 WRITE (LOUT,155)
155 FORMAT(//44H PROBLEM APPEARS UNSOLVABLE WITH GIVEN INPUT//)
GO TO 500
C
200 WRITE (LOUT,205) T,H
205 FORMAT(//35H KFLAG = -2 FROM INTEGRATOR AT T = ,E16.8,5H H =,
1 E16.8/52H THE REQUESTED ERROR IS SMALLER THAN CAN BE HANDLED//)
GO TO 500
C
250 WRITE (LOUT,255) T
255 FORMAT(//37H INTEGRATION HALTED BY DRIVER AT T = ,E16.8/
1 56H EPS TOO SMALL TO BE ATTAINED FOR THE MACHINE PRECISION/)
KFLAG = -2
GO TO 500
C
300 WRITE (LOUT,305) T
305 FORMAT(//35H KFLAG = -3 FROM INTEGRATOR AT T = ,E16.8/
1 45H CORRECTOR CONVERGENCE COULD NOT BE ACHIEVED/)
GO TO 110
C
400 WRITE (LOUT,405)
405 FORMAT(//28H ILLEGAL INPUT.. EPS .LE. 0.//)
INDEX = -4
RETURN
C
410 WRITE (LOUT,415)
415 FORMAT(//25H ILLEGAL INPUT.. N .LE. 0//)
INDEX = -4
RETURN
C
420 WRITE (LOUT,425)
425 FORMAT(//36H ILLEGAL INPUT.. (T0-TOUT)*H .GE. 0.//)
INDEX = -4
RETURN
C
430 WRITE (LOUT,435) INDEX
435 FORMAT(//24H ILLEGAL INPUT.. INDEX =,I5//)
INDEX = -4
RETURN
C
440 WRITE(LOUT,445) T,TOUT,H
445 FORMAT(//44H INDEX = -1 ON INPUT WITH (T-TOUT)*H .GE. 0./
1 4H T =,E16.8,9H TOUT =,E16.8,6H H =,E16.8/
1 44H INTERPOLATION WAS DONE AS ON NORMAL RETURN./
2 41H DESIRED PARAMETER CHANGES WERE NOT MADE.)
CALL INTERP (TOUT,WORK(IW22),N,U0,WORK,XM,YM)
INDEX = -5
RETURN
C
500 TOUT = T
DO 510 I = 1,N
IM1 = I - 1
IWA = IW22 + IM1
510 U0(I) = WORK(IWA)
520 INDEX = KFLAG
TOUTP = TOUT
H0 = HUSED
IF (KFLAG .NE. 0) H0 = H
RETURN
C----------------------- END OF SUBROUTINE DRIVEP ----------------------
END
SUBROUTINE INTERP (TOUT,U,N,U0,WORK,XM,YM) INTERP
INTEGER N,IDUMMY,JSTART,I,L,J
REAL TOUT,U,U0,T,H,DUMMY,S,S1
DIMENSION U0(N),U(N,1),WORK(1),XM(1),YM(1)
COMMON /GEAR1/ T,H,DUMMY(4),IDUMMY(3),JSTART
COMMON /IPARAM/ NPDE,NX,NY,MORDER,LWORK,LIWORK
COMMON /ISTORE/ IW1,IW2,IW3,IW4,IW5,IW6,IW7,IW8,IW9,IW10,IW11,
* IW12,IW13,IW14,IW15,IW16,IW17,IW18,IW19,IW20,IW21,
* IW22,IW23,IW24,IW25,IW26,IW27
LOGICAL IX,JY
C-----------------------------------------------------------------------
C SUBROUTINE INTERP COMPUTES INTERPOLATED VALUES OF THE DEPENDENT
C VARIABLE U AND STORES THEM IN U0. THE INTERPOLATION IS TO THE
C POINT T = TOUT, AND USES THE NORDSIECK HISTORY ARRAY U, AS FOLLOWS..
C NQ
C U0(I) = SUM U(I,J+1)*S**J ,
C J=0
C WHERE S = -(T-TOUT)/H. FOR THE PACKAGE PDETWO, IT IS MODIFIED TO
C DETERMINE ALSO THE SOLUTIONS FOR CONSTANT, TIME DEPENDENT BOUNDARY
C CONDITIONS.
C-----------------------------------------------------------------------
DO 10 I = 1,N
10 U0(I) = U(I,1)
L = JSTART + 1
S = (TOUT - T)/H
S1 = 1.0
DO 30 J = 2,L
S1 = S1*S
DO 20 I = 1,N
20 U0(I) = U0(I) + S1*U(I,J)
30 CONTINUE
C
C DETERMINE THE SOLUTIONS FOR CONSTANT, TIME DEPENDENT BOUNDARY
C CONDITIONS
C
DO 60 J=1,NY
JNX = (J-1)*NX
JY = J.EQ.1.OR.J.EQ.NY
DO 60 I=1,NX
IJNPDE = (I+JNX-1)*NPDE
IJ1 = 1 + IJNPDE
IX = I.EQ.1.OR.I.EQ.NX
IF (JY) CALL BNDRYH (TOUT,XM(I),YM(J),U0(IJ1),WORK(IW4),
* WORK(IW5),WORK(IW6),NPDE)
IF (IX) CALL BNDRYV (TOUT,XM(I),YM(J),U0(IJ1),WORK(IW7),
* WORK(IW8),WORK(IW9),NPDE)
DO 50 K=1,NPDE
KIJ = K + IJNPDE
KM1 = K - 1
KIW4 = IW4 + KM1
KIW5 = IW5 + KM1
KIW6 = IW6 + KM1
KIW7 = IW7 + KM1
KIW8 = IW8 + KM1
KIW9 = IW9 + KM1
IF (IX.AND.JY) GO TO 40
IF (JY.AND.WORK(KIW5).EQ.0.0) U0(KIJ) = WORK(KIW6)/
* WORK(KIW4)
IF (IX.AND.WORK(KIW8).EQ.0.0) U0(KIJ) = WORK(KIW9)/
* WORK(KIW7)
GO TO 50
40 IF (WORK(KIW5).EQ.0.0.AND.WORK(KIW8).NE.0.0) U0(KIJ) =
* WORK(KIW6)/WORK(KIW4)
IF (WORK(KIW5).NE.0.0.AND.WORK(KIW8).EQ.0.0) U0(KIJ) =
* WORK(KIW9)/WORK(KIW7)
IF (WORK(KIW5).EQ.0.0.AND.WORK(KIW8).EQ.0.0) U0(KIJ) = 0.5*
* ((WORK(KIW6)/WORK(KIW4)) + (WORK(KIW9)/WORK(KIW7)))
50 CONTINUE
60 CONTINUE
RETURN
C----------------------- END OF SUBROUTINE INTERP ----------------------
END
SUBROUTINE STIFFP (U,UMAX,ERROR,SAVE1,SAVE2,PW,N,WORK,IWORK, STIFFP
* XM,YM)
INTEGER N,MF,KFLAG,JSTART,ML,MU,MW,NEBAND,
1 NEB1,NNEB,NQUSED,NSTEP,NFE,NJE
INTEGER I,METH,MITER,NQ,L,IDOUB,MFOLD,NOLD,IRET,MEO,MIO,
1 IWEVAL,MAXDER,LMAX,IREDO,J,NSTEPJ,J1,J2,M,IER,NEWQ
INTEGER IWORK
REAL U,T,H,HMIN,HMAX,EPS,UROUND,UMAX,ERROR,
1 SAVE1,SAVE2,PW,EPSJ,HUSED
REAL EL,OLDL0,TOLD,RMAX,RC,CRATE,EPSOLD,
1 HOLD,FN,EDN,E,EUP,BND,RH,R1,CON,R,HL0,R0,D,PHL0,
1 PR3,D1,ENQ3,ENQ2,PR2,PR1,ENQ1
REAL WORK,XM,YM
REAL TQ
COMMON /GEAR1/ T,H,HMIN,HMAX,EPS,UROUND,NC,MF,KFLAG,JSTART
COMMON /GEAR2/ EPSJ,ML,MU,MW,NEBAND,NEB1,NNEB
COMMON /GEAR3/ HUSED,NQUSED,NSTEP,NFE,NJE
COMMON /IPARAM/ NPDE,NX,NY,MORDER,LWORK,LIWORK
COMMON /ISTORE/ IW1,IW2,IW3,IW4,IW5,IW6,IW7,IW8,IW9,IW10,IW11,
* IW12,IW13,IW14,IW15,IW16,IW17,IW18,IW19,IW20,IW21,
* IW22,IW23,IW24,IW25,IW26,IW27
DIMENSION U(N,1),UMAX(N),ERROR(N),SAVE1(N),SAVE2(N),PW(1),
* WORK(1),IWORK(1),XM(1),YM(1)
C-----------------------------------------------------------------------
C STIFFP PERFORMS ONE STEP OF THE INTEGRATION OF AN INITIAL VALUE
C PROBLEM FOR A SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS.
C STIFFP IS A VERSION FOR BANDED FORM OF THE JACOBIAN MATRIX.
C COMMUNICATION WITH STIFFP IS DONE WITH THE FOLLOWING VARIABLES..
C
C U AN N BY LMAX ARRAY CONTAINING THE DEPENDENT VARIABLES
C AND THEIR SCALED DERIVATIVES. LMAX IS 13 FOR THE ADAMS
C METHODS AND 6 FOR THE GEAR METHODS. LMAX - 1 = MAXDER
C IS THE MAXIMUM ORDER AVAILABLE. SEE SUBROUTINE COSET.
C U(I,J+1) CONTAINS THE J-TH DERIVATIVE OF U(I), SCALED BY
C H**J/FACTORIAL(J) (J = 0,1,...,NQ).
C UMAX AN ARRAY OF N ELEMENTS WITH WHICH THE ESTIMATED LOCAL
C ERRORS IN U ARE COMPARED.
C ERROR AN ARRAY OF N ELEMENTS. ERROR(I)/TQ(2) IS THE ESTIMATED
C ONE-STEP ERROR IN U(I).
C SAVE1, TWO ARRAYS OF WORKING STORAGE,
C SAVE2 EACH OF LENGTH N.
C PW A BLOCK OF LOCATIONS USED FOR PARTIAL DERIVATIVES IF
C MITER IS NOT 0. SEE DESCRIPTION IN DRIVER.
C N THE NUMBER OF FIRST-ORDER DIFFERENTIAL EQUATIONS.
C WORK FLOATING POINT WORKING ARRAY. SEE STRSET FOR ADDITIONAL
C DETAILS.
C IWORK AN INTEGER ARRAY OF LENGTH N USED FOR PIVOT
C INFORMATION IF MITER = 1 OR 2.
C XM ARE THE MESH POINTS IN THE HORIZONTAL DIRECTION.
C YM ARE THE MESH POINTS IN THE VERTICAL DIRECTION.
C T THE INDEPENDENT VARIABLE. T IS UPDATED ON EACH STEP TAKEN.
C H THE STEP SIZE TO BE ATTEMPTED ON THE NEXT STEP.
C H IS ALTERED BY THE ERROR CONTROL ALGORITHM DURING THE
C PROBLEM. H CAN BE EITHER POSITIVE OR NEGATIVE, BUT ITS
C SIGN MUST REMAIN CONSTANT THROUGHOUT THE PROBLEM.
C HMIN, THE MINIMUM AND MAXIMUM ABSOLUTE VALUE OF THE STEP SIZE
C HMAX TO BE USED FOR THE STEP. THESE MAY BE CHANGED AT ANY
C TIME, BUT WILL NOT TAKE EFFECT UNTIL THE NEXT H CHANGE.
C EPS THE RELATIVE ERROR BOUND. SEE DESCRIPTION IN DRIVER.
C UROUND THE UNIT ROUNDOFF OF THE MACHINE.
C MF THE METHOD FLAG. SEE DESCRIPTION IN DRIVER.
C KFLAG A COMPLETION CODE WITH THE FOLLOWING MEANINGS..
C 0 THE STEP WAS SUCCESFUL.
C -1 THE REQUESTED ERROR COULD NOT BE ACHIEVED
C WITH ABS(H) = HMIN.
C -2 THE REQUESTED ERROR IS SMALLER THAN CAN
C BE HANDLED FOR THIS PROBLEM.
C -3 CORRECTOR CONVERGENCE COULD NOT BE
C ACHIEVED FOR ABS(H) = HMIN.
C ON A RETURN WITH KFLAG NEGATIVE, THE VALUES OF T AND
C THE U ARRAY ARE AS OF THE BEGINNING OF THE LAST
C STEP, AND H IS THE LAST STEP SIZE ATTEMPTED.
C JSTART AN INTEGER USED ON INPUT AND OUTPUT.
C ON INPUT, IT HAS THE FOLLOWING VALUES AND MEANINGS..
C 0 PERFORM THE FIRST STEP.
C .GT.0 TAKE A NEW STEP CONTINUING FROM THE LAST.
C .LT.0 TAKE THE NEXT STEP WITH A NEW VALUE OF
C H, EPS, N, AND/OR MF.
C ON EXIT, JSTART IS NQ, THE CURRENT ORDER OF THE METHOD.
C ML,MU THE LOWER AND UPPER HALF BANDWIDTHS, RESPECTIVELY, OF
C THE JACOBIAN. SEE DESCRIPTION IN DRIVER.
C NEBAND THE EXTENDED BANDWIDTH, 2*ML + MU + 1.
C NNEB THE MAXIMUM SIZE OF PW USED, N*NEBAND.
C-----------------------------------------------------------------------
DIMENSION EL(13),TQ(4)
DATA EL(2)/1.0/, OLDL0/1.0/
KFLAG = 0
TOLD = T
IF (JSTART .GT. 0) GO TO 200
IF (JSTART .NE. 0) GO TO 120
C-----------------------------------------------------------------------
C ON THE FIRST CALL, THE ORDER IS SET TO 1 AND THE INITIAL UDOT IS
C CALCULATED. RMAX IS THE MAXIMUM RATIO BY WHICH H CAN BE INCREASED
C IN A SINGLE STEP. IT IS INITIALLY 1.E4 TO COMPENSATE FOR THE SMALL
C INITIAL H, BUT THEN IS NORMALLY EQUAL TO 10. IF A FAILURE
C OCCURS (IN CORRECTOR CONVERGENCE OR ERROR TEST), RMAX IS SET AT 2
C FOR THE NEXT INCREASE.
C-----------------------------------------------------------------------
CALL PDETWO (NPDE,NX,NY,XM,YM, T, U, SAVE1, WORK,WORK(IW1),
* WORK(IW2),WORK(IW3),WORK(IW4),WORK(IW5),WORK(IW6),
* WORK(IW7),WORK(IW8),WORK(IW9),WORK(IW10),WORK(IW11),
* WORK(IW12),WORK(IW13),WORK(IW14),WORK(IW15),WORK(IW16),
* WORK(IW17),WORK(IW18),WORK(IW19),WORK(IW20),WORK(IW21))
DO 110 I = 1,N
110 U(I,2) = H*SAVE1(I)
METH = MF/10
MITER = MF - 10*METH
NQ = 1
L = 2
IDOUB = 3
RMAX = 1.E+4
RC = 0.0
CRATE = 1.0
EPSOLD = EPS
HOLD = H
MFOLD = MF
NOLD = N
NSTEP = 0
NSTEPJ = 0
NFE = 1
NJE = 0
IRET = 3
IF (MITER .NE. 1) GO TO 130
DO 115 I = 1,NNEB
115 PW(I) = 0.0
GO TO 130
C-----------------------------------------------------------------------
C IF THE CALLER HAS CHANGED METH, COSET IS CALLED TO SET
C THE COEFFICIENTS OF THE METHOD. IF THE CALLER HAS CHANGED
C N, EPS, OR METH, THE CONSTANTS E, EDN, EUP, AND BND MUST BE RESET.
C E IS A COMPARISON FOR ERRORS OF THE CURRENT ORDER NQ. EUP IS
C TO TEST FOR INCREASING THE ORDER, EDN FOR DECREASING THE ORDER.
C BND IS USED TO TEST FOR CONVERGENCE OF THE CORRECTOR ITERATES.
C IF THE CALLER HAS CHANGED H, U MUST BE RESCALED.
C IF H OR METH HAS BEEN CHANGED, IDOUB IS RESET TO L + 1 TO PREVENT
C FURTHER CHANGES IN H FOR THAT MANY STEPS.
C-----------------------------------------------------------------------
120 IF (MF .EQ. MFOLD) GO TO 150
MEO = METH
MIO = MITER
METH = MF/10
MITER = MF - 10*METH
MFOLD = MF
IF (MITER .NE. MIO) IWEVAL = MITER
IF (METH .EQ. MEO) GO TO 150
IDOUB = L + 1
IRET = 1
130 CALL COSET (METH, NQ, EL, TQ, MAXDER)
MAXDER = MIN0(MAXDER,MORDER)
LMAX = MAXDER + 1
RC = RC*EL(1)/OLDL0
OLDL0 = EL(1)
140 FN = FLOAT(N)
EDN = FN*(TQ(1)*EPS)**2
E = FN*(TQ(2)*EPS)**2
EUP = FN*(TQ(3)*EPS)**2
BND = FN*(TQ(4)*EPS)**2
GO TO (160, 170, 200), IRET
150 IF ((EPS .EQ. EPSOLD) .AND. (N .EQ. NOLD)) GO TO 160
EPSOLD = EPS
NOLD = N
IRET = 1
GO TO 140
160 IF (H .EQ. HOLD) GO TO 200
RH = H/HOLD
H = HOLD
IREDO = 3
GO TO 175
170 RH = AMAX1(RH,HMIN/ABS(H))
175 RH = AMIN1(RH,HMAX/ABS(H),RMAX)
R1 = 1.0
DO 180 J = 2,L
R1 = R1*RH
DO 180 I = 1,N
180 U(I,J) = U(I,J)*R1
H = H*RH
RC = RC*RH
IDOUB = L + 1
IF (IREDO .EQ. 0) GO TO 690
C-----------------------------------------------------------------------
C THIS SECTION COMPUTES THE PREDICTED VALUES BY EFFECTIVELY
C MULTIPLYING THE U ARRAY BY THE PASCAL TRIANGLE MATRIX.
C RC IS THE RATIO OF NEW TO OLD VALUES OF THE COEFFICIENT H*EL(1).
C WHEN RC DIFFERS FROM 1 BY MORE THAN 30 PERCENT, OR THE CALLER HAS
C CHANGED MITER, IWEVAL IS SET TO MITER TO FORCE THE PARTIALS TO BE
C UPDATED, IF PARTIALS ARE USED. IN ANY CASE, THE PARTIALS
C ARE UPDATED AT LEAST EVERY 20-TH STEP.
C-----------------------------------------------------------------------
200 IF (ABS(RC-1.0) .GT. 0.30) IWEVAL = MITER
IF (NSTEP .GE. NSTEPJ+20) IWEVAL = MITER
T = T + H
DO 210 J1 = 1,NQ
DO 210 J2 = J1,NQ
J = (NQ + J1) - J2
DO 210 I = 1,N
210 U(I,J) = U(I,J) + U(I,J+1)
C-----------------------------------------------------------------------
C UP TO 3 CORRECTOR ITERATIONS ARE TAKEN. A CONVERGENCE TEST IS
C MADE ON THE R.M.S. NORM OF EACH CORRECTION, USING BND, WHICH
C IS DEPENDENT ON EPS. THE SUM OF THE CORRECTIONS IS ACCUMULATED
C IN THE VECTOR ERROR(I). THE U ARRAY IS NOT ALTERED IN THE CORRECTOR
C LOOP. THE UPDATED U VECTOR IS STORED TEMPORARILY IN SAVE1.
C-----------------------------------------------------------------------
220 DO 230 I = 1,N
230 ERROR(I) = 0.0
M = 0
CALL PDETWO (NPDE,NX,NY,XM,YM, T, U, SAVE2, WORK,WORK(IW1),
* WORK(IW2),WORK(IW3),WORK(IW4),WORK(IW5),WORK(IW6),
* WORK(IW7),WORK(IW8),WORK(IW9),WORK(IW10),WORK(IW11),
* WORK(IW12),WORK(IW13),WORK(IW14),WORK(IW15),WORK(IW16),
* WORK(IW17),WORK(IW18),WORK(IW19),WORK(IW20),WORK(IW21))
NFE = NFE + 1
IF (IWEVAL .LE. 0) GO TO 290
C-----------------------------------------------------------------------
C IF INDICATED, THE MATRIX P = I - H*EL(1)*J IS REEVALUATED BEFORE
C STARTING THE CORRECTOR ITERATION. IWEVAL IS SET TO 0 AS AN
C INDICATOR THAT THIS HAS BEEN DONE. IF MITER = 1 OR 2, P IS
C COMPUTED AND PROCESSED IN PSETB. IF MITER = 3, THE MATRIX USED
C IS P = I - H*EL(1)*D, WHERE D IS A DIAGONAL MATRIX.
C-----------------------------------------------------------------------
IWEVAL = 0
RC = 1.0
NJE = NJE + 1
NSTEPJ = NSTEP
GO TO (250, 240, 260), MITER
C240 NFE = NFE + MW
240 NFE = NFE + 5 * NPDE
250 CON = -H*EL(1)
CALL PSETM (NPDE,NX,NY,XM,YM,U,UMAX,ERROR,SAVE1,SAVE2,PW,CON,
* MITER,IER,NEBAND,WORK,IWORK)
IF (IER .NE. 0) GO TO 420
GO TO 350
260 R = EL(1)*.10
DO 270 I = 1,N
270 PW(I) = U(I,1) + R*(H*SAVE2(I) - U(I,2))
CALL PDETWO (NPDE,NX,NY,XM,YM, T, PW, SAVE1, WORK,WORK(IW1),
* WORK(IW2),WORK(IW3),WORK(IW4),WORK(IW5),WORK(IW6),
* WORK(IW7),WORK(IW8),WORK(IW9),WORK(IW10),WORK(IW11),
* WORK(IW12),WORK(IW13),WORK(IW14),WORK(IW15),WORK(IW16),
* WORK(IW17),WORK(IW18),WORK(IW19),WORK(IW20),WORK(IW21))
NFE = NFE + 1
HL0 = H*EL(1)
DO 280 I = 1,N
R0 = H*SAVE2(I) - U(I,2)
PW(I) = 1.0
D = .10*R0 - H*(SAVE1(I) - SAVE2(I))
SAVE1(I) = 0.0
IF (ABS(R0) .LT. UROUND*UMAX(I)) GO TO 280
IF (ABS(D) .EQ. 0.0) GO TO 420
PW(I) = .10*R0/D
SAVE1(I) = PW(I)*R0
280 CONTINUE
GO TO 370
290 IF (MITER .NE. 0) GO TO (350, 350, 310), MITER
C-----------------------------------------------------------------------
C IN THE CASE OF FUNCTIONAL ITERATION, UPDATE U DIRECTLY FROM
C THE RESULT OF THE LAST DIFFUN CALL.
C-----------------------------------------------------------------------
D = 0.0
DO 300 I = 1,N
R = H*SAVE2(I) - U(I,2)
D = D + ( (R-ERROR(I))/UMAX(I) )**2
SAVE1(I) = U(I,1) + EL(1)*R
300 ERROR(I) = R
GO TO 400
C-----------------------------------------------------------------------
C IN THE CASE OF THE CHORD METHOD, COMPUTE THE CORRECTOR ERROR,
C F SUB (M), AND SOLVE THE LINEAR SYSTEM WITH THAT AS RIGHT-HAND
C SIDE AND P AS COEFFICIENT MATRIX, USING THE LU DECOMPOSITION
C IF MITER = 1 OR 2. IF MITER = 3, THE COEFFICIENT H*EL(1)
C IN P IS UPDATED.
C-----------------------------------------------------------------------
310 PHL0 = HL0
HL0 = H*EL(1)
IF (HL0 .EQ. PHL0) GO TO 330
R = HL0/PHL0
DO 320 I = 1,N
D = 1.0 - R*(1.0 - 1.0/PW(I))
IF (ABS(D) .EQ. 0.0) GO TO 440
320 PW(I) = 1.0/D
330 DO 340 I = 1,N
340 SAVE1(I) = PW(I)*(H*SAVE2(I) - (U(I,2) + ERROR(I)))
GO TO 370
350 DO 360 I = 1,N
360 SAVE1(I) = H*SAVE2(I) - (U(I,2) + ERROR(I))
CALL SOLBR (NEBAND,N,ML,MU,PW,IWORK,SAVE1)
370 D = 0.0
DO 380 I = 1,N
ERROR(I) = ERROR(I) + SAVE1(I)
D = D + (SAVE1(I)/UMAX(I))**2
380 SAVE1(I) = U(I,1) + EL(1)*ERROR(I)
C-----------------------------------------------------------------------
C TEST FOR CONVERGENCE. IF M.GT.0, AN ESTIMATE OF THE CONVERGENCE
C RATE CONSTANT IS STORED IN CRATE, AND THIS IS USED IN THE TEST.
C-----------------------------------------------------------------------
400 IF (M .NE. 0) CRATE = AMAX1(.90*CRATE,D/D1)
IF ((D*AMIN1(1.0,2.0*CRATE)) .LE. BND) GO TO 450
D1 = D
M = M + 1
IF (M .EQ. 3) GO TO 410
CALL PDETWO (NPDE,NX,NY,XM,YM, T, SAVE1, SAVE2, WORK,WORK(IW1),
* WORK(IW2),WORK(IW3),WORK(IW4),WORK(IW5),WORK(IW6),
* WORK(IW7),WORK(IW8),WORK(IW9),WORK(IW10),WORK(IW11),
* WORK(IW12),WORK(IW13),WORK(IW14),WORK(IW15),WORK(IW16),
* WORK(IW17),WORK(IW18),WORK(IW19),WORK(IW20),WORK(IW21))
GO TO 290
C-----------------------------------------------------------------------
C THE CORRECTOR ITERATION FAILED TO CONVERGE IN 3 TRIES. IF PARTIALS
C ARE INVOLVED BUT ARE NOT UP TO DATE, THEY ARE REEVALUATED FOR THE
C NEXT TRY. OTHERWISE THE U ARRAY IS RETRACTED TO ITS VALUES
C BEFORE PREDICTION, AND H IS REDUCED, IF POSSIBLE. IF NOT, A
C NO-CONVERGENCE EXIT IS TAKEN.
C-----------------------------------------------------------------------
410 NFE = NFE + 2
IF (IWEVAL .EQ. -1) GO TO 440
420 T = TOLD
RMAX = 2.0
DO 430 J1 = 1,NQ
DO 430 J2 = J1,NQ
J = (NQ + J1) - J2
DO 430 I = 1,N
430 U(I,J) = U(I,J) - U(I,J+1)
IF (ABS(H) .LE. HMIN*1.00001) GO TO 680
RH = .250
IREDO = 1
GO TO 170
440 IWEVAL = MITER
GO TO 220
C-----------------------------------------------------------------------
C THE CORRECTOR HAS CONVERGED. IWEVAL IS SET TO -1 IF PARTIAL
C DERIVATIVES WERE USED, TO SIGNAL THAT THEY MAY NEED UPDATING ON
C SUBSEQUENT STEPS. THE ERROR TEST IS MADE AND CONTROL PASSES TO
C STATEMENT 500 IF IT FAILS.
C-----------------------------------------------------------------------
450 IF (MITER .NE. 0) IWEVAL = -1
NFE = NFE + M
D = 0.0
DO 460 I = 1,N
460 D = D + (ERROR(I)/UMAX(I))**2
IF (D .GT. E) GO TO 500
C-----------------------------------------------------------------------
C AFTER A SUCCESSFUL STEP, UPDATE THE U ARRAY.
C CONSIDER CHANGING H IF IDOUB = 1. OTHERWISE DECREASE IDOUB BY 1.
C IF IDOUB IS THEN 1 AND NQ .LT. MAXDER, THEN ERROR IS SAVED FOR
C USE IN A POSSIBLE ORDER INCREASE ON THE NEXT STEP.
C IF A CHANGE IN H IS CONSIDERED, AN INCREASE OR DECREASE IN ORDER
C BY ONE IS CONSIDERED ALSO. A CHANGE IN H IS MADE ONLY IF IT IS BY A
C FACTOR OF AT LEAST 1.1. IF NOT, IDOUB IS SET TO 10 TO PREVENT
C TESTING FOR THAT MANY STEPS.
C-----------------------------------------------------------------------
KFLAG = 0
IREDO = 0
NSTEP = NSTEP + 1
HUSED = H
NQUSED = NQ
DO 470 J = 1,L
DO 470 I = 1,N
470 U(I,J) = U(I,J) + EL(J)*ERROR(I)
IF (IDOUB .EQ. 1) GO TO 520
IDOUB = IDOUB - 1
IF (IDOUB .GT. 1) GO TO 700
IF (L .EQ. LMAX) GO TO 700
DO 490 I = 1,N
490 U(I,LMAX) = ERROR(I)
GO TO 700
C-----------------------------------------------------------------------
C THE ERROR TEST FAILED. KFLAG KEEPS TRACK OF MULTIPLE FAILURES.
C RESTORE T AND THE U ARRAY TO THEIR PREVIOUS VALUES, AND PREPARE
C TO TRY THE STEP AGAIN. COMPUTE THE OPTIMUM STEP SIZE FOR THIS OR
C ONE LOWER ORDER.
C-----------------------------------------------------------------------
500 KFLAG = KFLAG - 1
T = TOLD
DO 510 J1 = 1,NQ
DO 510 J2 = J1,NQ
J = (NQ + J1) - J2
DO 510 I = 1,N
510 U(I,J) = U(I,J) - U(I,J+1)
RMAX = 2.0
IF (ABS(H) .LE. HMIN*1.000010) GO TO 660
IF (KFLAG .LE. -3) GO TO 640
IREDO = 2
PR3 = 1.E+20
GO TO 540
C-----------------------------------------------------------------------
C REGARDLESS OF THE SUCCESS OR FAILURE OF THE STEP, FACTORS
C PR1, PR2, AND PR3 ARE COMPUTED, BY WHICH H COULD BE DIVIDED
C AT ORDER NQ - 1, ORDER NQ, OR ORDER NQ + 1, RESPECTIVELY.
C IN THE CASE OF FAILURE, PR3 = 1.E20 TO AVOID AN ORDER INCREASE.
C THE SMALLEST OF THESE IS DETERMINED AND THE NEW ORDER CHOSEN
C ACCORDINGLY. IF THE ORDER IS TO BE INCREASED, WE COMPUTE ONE
C ADDITIONAL SCALED DERIVATIVE.
C-----------------------------------------------------------------------
520 PR3 = 1.E+20
IF (L .EQ. LMAX) GO TO 540
D1 = 0.0
DO 530 I = 1,N
530 D1 = D1 + ((ERROR(I) - U(I,LMAX))/UMAX(I))**2
ENQ3 = .50/FLOAT(L+1)
PR3 = ((D1/EUP)**ENQ3)*1.40 + 1.4E-6
540 ENQ2 = .50/FLOAT(L)
PR2 = ((D/E)**ENQ2)*1.20 + 1.2E-6
PR1 = 1.E+20
IF (NQ .EQ. 1) GO TO 560
D = 0.0
DO 550 I = 1,N
550 D = D + (U(I,L)/UMAX(I))**2
ENQ1 = .50/FLOAT(NQ)
PR1 = ((D/EDN)**ENQ1)*1.30 + 1.3E-6
560 IF (PR2 .LE. PR3) GO TO 570
IF (PR3 .LT. PR1) GO TO 590
GO TO 580
570 IF (PR2 .GT. PR1) GO TO 580
NEWQ = NQ
RH = 1.0/PR2
GO TO 620
580 NEWQ = NQ - 1
RH = 1.0/PR1
GO TO 620
590 NEWQ = L
RH = 1.0/PR3
IF (RH .LT. 1.10) GO TO 610
DO 600 I = 1,N
600 U(I,NEWQ+1) = ERROR(I)*EL(L)/FLOAT(L)
GO TO 630
610 IDOUB = 10
GO TO 700
620 IF ((KFLAG .EQ. 0) .AND. (RH .LT. 1.10)) GO TO 610
C-----------------------------------------------------------------------
C IF THERE IS A CHANGE OF ORDER, RESET NQ, L, AND THE COEFFICIENTS.
C IN ANY CASE H IS RESET ACCORDING TO RH AND THE U ARRAY IS RESCALED.
C THEN EXIT FROM 690 IF THE STEP WAS OK, OR REDO THE STEP OTHERWISE.
C-----------------------------------------------------------------------
IF (NEWQ .EQ. NQ) GO TO 170
630 NQ = NEWQ
L = NQ + 1
IRET = 2
GO TO 130
C-----------------------------------------------------------------------
C CONTROL REACHES THIS SECTION IF 3 OR MORE FAILURES HAVE OCCURED.
C IT IS ASSUMED THAT THE DERIVATIVES THAT HAVE ACCUMULATED IN THE
C U ARRAY HAVE ERRORS OF THE WRONG ORDER. HENCE THE FIRST
C DERIVATIVE IS RECOMPUTED, AND THE ORDER IS SET TO 1. THEN
C H IS REDUCED BY A FACTOR OF 10, AND THE STEP IS RETRIED.
C AFTER A TOTAL OF 7 FAILURES, AN EXIT IS TAKEN WITH KFLAG = -2.
C-----------------------------------------------------------------------
640 IF (KFLAG .EQ. -7) GO TO 670
RH = .10
RH = AMAX1(HMIN/ABS(H),RH)
H = H*RH
CALL PDETWO (NPDE,NX,NY,XM,YM, T, U, SAVE1, WORK,WORK(IW1),
* WORK(IW2),WORK(IW3),WORK(IW4),WORK(IW5),WORK(IW6),
* WORK(IW7),WORK(IW8),WORK(IW9),WORK(IW10),WORK(IW11),
* WORK(IW12),WORK(IW13),WORK(IW14),WORK(IW15),WORK(IW16),
* WORK(IW17),WORK(IW18),WORK(IW19),WORK(IW20),WORK(IW21))
NFE = NFE + 1
DO 650 I = 1,N
650 U(I,2) = H*SAVE1(I)
IWEVAL = MITER
IDOUB = 10
IF (NQ .EQ. 1) GO TO 200
NQ = 1
L = 2
IRET = 3
GO TO 130
C-----------------------------------------------------------------------
C ALL RETURNS ARE MADE THROUGH THIS SECTION. H IS SAVED IN HOLD
C TO ALLOW THE CALLER TO CHANGE H ON THE NEXT STEP.
C-----------------------------------------------------------------------
660 KFLAG = -1
GO TO 700
670 KFLAG = -2
GO TO 700
680 KFLAG = -3
GO TO 700
690 RMAX = 10.0
700 HOLD = H
JSTART = NQ
RETURN
C----------------------- END OF SUBROUTINE STIFFP ----------------------
END
SUBROUTINE COSET (METH, NQ, EL, TQ, MAXDER) COSET
C
INTEGER METH,NQ,MAXDER,K
REAL EL
REAL TQ,PERTST
C-----------------------------------------------------------------------
C COSET IS CALLED BY THE INTEGRATOR AND SETS COEFFICIENTS USED THERE.
C THE VECTOR EL, OF LENGTH NQ + 1, DETERMINES THE BASIC METHOD.
C THE VECTOR TQ, OF LENGTH 4, IS INVOLVED IN ADJUSTING THE STEP SIZE
C IN RELATION TO TRUNCATION ERROR. ITS VALUES ARE GIVEN BY THE
C PERTST ARRAY.
C THE VECTORS EL AND TQ DEPEND ON METH AND NQ.
C COSET ALSO SETS MAXDER, THE MAXIMUM ORDER OF THE METHOD AVAILABLE.
C CURRENTLY IT IS 12 FOR THE ADAMS METHODS AND 5 FOR THE BDF METHODS.
C LMAX = MAXDER + 1 IS THE NUMBER OF COLUMNS IN THE U ARRAY.
C THE MAXIMUM ORDER USED MAY BE REDUCED SIMPLY BY DECREASING THE
C NUMBERS IN STATEMENTS 1 AND/OR 2 BELOW.
C
C THE COEFFICIENTS IN PERTST NEED BE GIVEN TO ONLY ABOUT
C ONE PERCENT ACCURACY. THE ORDER IN WHICH THE GROUPS APPEAR BELOW
C IS.. COEFFICIENTS FOR ORDER NQ - 1, COEFFICIENTS FOR ORDER NQ,
C COEFFICIENTS FOR ORDER NQ + 1. WITHIN EACH GROUP ARE THE
C COEFFICIENTS FOR THE ADAMS METHODS, FOLLOWED BY THOSE FOR THE
C BDF METHODS.
C-----------------------------------------------------------------------
DIMENSION PERTST(12,2,3),EL(13),TQ(4)
DATA PERTST(1,1,1),PERTST(2,1,1),PERTST(3,1,1),PERTST(4,1,1)
1 /1., 1., 2., 1./,
2 PERTST(5,1,1),PERTST(6,1,1),PERTST(7,1,1),PERTST(8,1,1)
3 /.3158, .07407, .01391, .002182/,
4 PERTST(9,1,1),PERTST(10,1,1),PERTST(11,1,1),PERTST(12,1,1)
5 /.0002945, .00003492, .000003692, .0000003524/
DATA PERTST(1,2,1),PERTST(2,2,1),PERTST(3,2,1),PERTST(4,2,1)
1 /1., 1., .5, .1667/,
2 PERTST(5,2,1),PERTST(6,2,1),PERTST(7,2,1),PERTST(8,2,1)
3 /.04167, 1., 1., 1./,
4 PERTST(9,2,1),PERTST(10,2,1),PERTST(11,2,1),PERTST(12,2,1)
5 /1., 1., 1., 1./
DATA PERTST(1,1,2),PERTST(2,1,2),PERTST(3,1,2),PERTST(4,1,2)
1 /2., 12., 24., 37.89/,
2 PERTST(5,1,2),PERTST(6,1,2),PERTST(7,1,2),PERTST(8,1,2)
3 /53.33, 70.08, 87.97, 106.9/,
4 PERTST(9,1,2),PERTST(10,1,2),PERTST(11,1,2),PERTST(12,1,2)
5 /126.7, 147.4, 168.8, 191./
DATA PERTST(1,2,2),PERTST(2,2,2),PERTST(3,2,2),PERTST(4,2,2)
1 /2., 4.5, 7.333, 10.42/,
2 PERTST(5,2,2),PERTST(6,2,2),PERTST(7,2,2),PERTST(8,2,2)
3 /13.7, 1., 1., 1./,
4 PERTST(9,2,2),PERTST(10,2,2),PERTST(11,2,2),PERTST(12,2,2)
5 /1., 1., 1., 1./
DATA PERTST(1,1,3),PERTST(2,1,3),PERTST(3,1,3),PERTST(4,1,3)
1 /12., 24., 37.89, 53.33/,
2 PERTST(5,1,3),PERTST(6,1,3),PERTST(7,1,3),PERTST(8,1,3)
3 /70.08, 87.97, 106.9, 126.7/,
4 PERTST(9,1,3),PERTST(10,1,3),PERTST(11,1,3),PERTST(12,1,3)
5 /147.4, 168.8, 191., 1./
DATA PERTST(1,2,3),PERTST(2,2,3),PERTST(3,2,3),PERTST(4,2,3)
1 /3., 6., 9.167, 12.5/,
2 PERTST(5,2,3),PERTST(6,2,3),PERTST(7,2,3),PERTST(8,2,3)
3 /1., 1., 1., 1./,
4 PERTST(9,2,3),PERTST(10,2,3),PERTST(11,2,3),PERTST(12,2,3)
5 /1., 1., 1., 1./
C
GO TO (1,2),METH
1 MAXDER = 12
GO TO (101,102,103,104,105,106,107,108,109,110,111,112),NQ
2 MAXDER = 5
GO TO (201,202,203,204,205),NQ
C-----------------------------------------------------------------------
C THE FOLLOWING COEFFICIENTS SHOULD BE DEFINED TO MACHINE ACCURACY.
C FOR A GIVEN ORDER NQ, THEY CAN BE CALCULATED BY USE OF THE
C GENERATING POLYNOMIAL L(T), WHOSE COEFFICIENTS ARE EL(I)..
C L(T) = EL(1) + EL(2)*T + ... + EL(NQ+1)*T**NQ.
C FOR THE IMPLICIT ADAMS METHODS, L(T) IS GIVEN BY
C DL/DT = (T+1)*(T+2)* ... *(T+NQ-1)/K, L(-1) = 0,
C WHERE K = FACTORIAL(NQ-1).
C FOR THE BDF METHODS,
C L(T) = (T+1)*(T+2)* ... *(T+NQ)/K,
C WHERE K = FACTORIAL(NQ)*(1 + 1/2 + ... + 1/NQ).
C
C THE ORDER IN WHICH THE GROUPS APPEAR BELOW IS..
C IMPLICIT ADAMS METHODS OF ORDERS 1 TO 12,
C BDF METHODS OF ORDERS 1 TO 5.
C-----------------------------------------------------------------------
101 EL(1) = 1.00
GO TO 900
102 EL(1) = 0.50
EL(3) = EL(1)
GO TO 900
103 EL(1) = 4.1666666666666667E-01
EL(3) = 0.750
EL(4) = 1.6666666666666667E-01
GO TO 900
104 EL(1) = 0.3750
EL(3) = 9.1666666666666667E-01
EL(4) = 3.3333333333333333E-01
EL(5) = 4.1666666666666667E-02
GO TO 900
105 EL(1) = 3.4861111111111111E-01
EL(3) = 1.04166666666666670
EL(4) = 4.8611111111111111E-01
EL(5) = 1.0416666666666667E-01
EL(6) = 8.3333333333333333E-03
GO TO 900
106 EL(1) = 3.2986111111111111E-01
EL(3) = 1.1416666666666667E+00
EL(4) = 0.625E+00
EL(5) = 1.7708333333333333E-01
EL(6) = 0.025E+00
EL(7) = 1.3888888888888889E-03
GO TO 900
107 EL(1) = 3.1559193121693122E-01
EL(3) = 1.225E+00
EL(4) = 7.5185185185185185E-01
EL(5) = 2.5520833333333333E-01
EL(6) = 4.8611111111111111E-02
EL(7) = 4.8611111111111111E-03
EL(8) = 1.9841269841269841E-04
GO TO 900
108 EL(1) = 3.0422453703703704E-01
EL(3) = 1.2964285714285714E+00
EL(4) = 8.6851851851851852E-01
EL(5) = 3.3576388888888889E-01
EL(6) = 7.7777777777777778E-02
EL(7) = 1.0648148148148148E-02
EL(8) = 7.9365079365079365E-04
EL(9) = 2.4801587301587302E-05
GO TO 900
109 EL(1) = 2.9486800044091711E-01
EL(3) = 1.3589285714285714E+00
EL(4) = 9.7655423280423280E-01
EL(5) = 0.4171875E+00
EL(6) = 1.1135416666666667E-01
EL(7) = 0.01875E+00
EL(8) = 1.9345238095238095E-03
EL(9) = 1.1160714285714286E-04
EL(10)= 2.7557319223985891E-06
GO TO 900
110 EL(1) = 2.8697544642857143E-01
EL(3) = 1.4144841269841270E+00
EL(4) = 1.0772156084656085E+00
EL(5) = 4.9856701940035273E-01
EL(6) = 0.1484375E+00
EL(7) = 2.9060570987654321E-02
EL(8) = 3.7202380952380952E-03
EL(9) = 2.9968584656084656E-04
EL(10)= 1.3778659611992945E-05
EL(11)= 2.7557319223985891E-07
GO TO 900
111 EL(1) = 2.8018959644393672E-01
EL(3) = 1.4644841269841270E+00
EL(4) = 1.1715145502645503E+00
EL(5) = 5.7935819003527337E-01
EL(6) = 1.8832286155202822E-01
EL(7) = 4.1430362654320988E-02
EL(8) = 6.2111441798941799E-03
EL(9) = 6.2520667989417989E-04
EL(10)= 4.0417401528512640E-05
EL(11)= 1.5156525573192240E-06
EL(12)= 2.5052108385441719E-08
GO TO 900
112 EL(1) = 2.7426554003159906E-01
EL(3) = 1.5099386724386724E+00
EL(4) = 1.2602711640211640E+00
EL(5) = 6.5923418209876543E-01
EL(6) = 2.3045800264550265E-01
EL(7) = 5.5697246105232216E-02
EL(8) = 9.4394841269841270E-03
EL(9) = 1.1192749669312169E-03
EL(10)= 9.0939153439153439E-05
EL(11)= 4.8225308641975309E-06
EL(12)= 1.5031265031265031E-07
EL(13)= 2.0876756987868099E-09
GO TO 900
C
201 EL(1) = 1.0E+00
GO TO 900
202 EL(1) = 6.6666666666666667E-01
EL(3) = 3.3333333333333333E-01
GO TO 900
203 EL(1) = 5.4545454545454545E-01
EL(3) = EL(1)
EL(4) = 9.0909090909090909E-02
GO TO 900
204 EL(1) = 0.48E+00
EL(3) = 0.7E+00
EL(4) = 0.2E+00
EL(5) = 0.02E+00
GO TO 900
205 EL(1) = 4.3795620437956204E-01
EL(3) = 8.2116788321167883E-01
EL(4) = 3.1021897810218978E-01
EL(5) = 5.4744525547445255E-02
EL(6) = 3.6496350364963504E-03
C
900 DO 910 K = 1,3
910 TQ(K) = PERTST(NQ,METH,K)
TQ(4) = .5*TQ(2)/FLOAT(NQ+2)
RETURN
C----------------------- END OF SUBROUTINE COSET -----------------------
END
SUBROUTINE DECBR (MDIM, N, ML, MU, B, IP, IER) DECBR
INTEGER MDIM,N,ML,MU,IP,IER
REAL B
DIMENSION B(MDIM,N), IP(N)
C-----------------------------------------------------------------------
C THIS SUBROUTINE CONSTRUCTS THE LU DECOMPOSITION OF A BAND MATRIX A
C IN THE FORM L*U = P*A , WHERE P IS A PERMUTATION MATRIX, L IS A
C UNIT LOWER TRIANGULAR MATRIX, AND U IS AN UPPER TRIANGULAR MATRIX.
C THE MATRIX A IS ASSUMED HERE TO BE STORED BY ROWS IN THE ARRAY B.
C MDIM = THE FIRST DIMENSION (COLUMN LENGTH) OF THE ARRAY B.
C MDIM MUST BE AT LEAST 2*ML+MU+1.
C N = ORDER OF MATRIX.
C ML,MU = WIDTHS OF THE LOWER AND UPPER PARTS OF THE BAND, NOT
C COUNTING THE MAIN DIAGONAL. TOTAL BANDWIDTH IS ML+MU+1.
C B = MDIM BY N ARRAY CONTAINING THE MATRIX A ON INPUT
C AND ITS FACTORED FORM ON OUTPUT.
C ON INPUT, A IS STORED BY ROWS IN B AS AN (ML+MU+1) BY N
C ARRAY (THE ROWS OF A ARE THE COLUMNS OF B), WITH A COLUMN
C LENGTH OF MDIM FOR B. THUS B(K,I) = B(K+(I-1)*MDIM)
C (1 .LE. K .LE. ML+MU+1) CONTAINS THE PART OF THE I-TH ROW
C OF A WITHIN THE BAND. THAT IS, B(J-I+ML+1,I) CONTAINS
C A(I,J) FOR -ML .LE. (J-I) .LE. MU.
C ON OUTPUT, B CONTAINS THE L AND U FACTORS, WITH U IN ROWS
C 1 TO ML+MU+1, AND -L IN ROWS ML+MU+2 TO 2*ML+MU+1.
C THE DIAGONAL OF L IS NOT STORED, AND THE RECIPROCALS
C OF THE DIAGONAL ELEMENTS OF U ARE STORED.
C IP = ARRAY OF LENGTH N CONTAINING PIVOT INFORMATION ON OUTPUT.
C IER = ERROR INDICATOR..
C = 0 IF NO ERROR,
C = -1 IF MDIM, N, ML, AND/OR MU IS ILLEGAL,
C = K IF THE K-TH PIVOT CHOSEN WAS ZERO (A IS SINGULAR).
C THE INPUT ARGUMENTS ARE MDIM, N, ML, MU, AND B.
C THE OUTPUT ARGUMENTS ARE B, IP, AND IER.
C-----------------------------------------------------------------------
INTEGER I,II,J,K,KK,LL,LR,MX,NP,NR,N1
REAL XM,XX
N1 = N - 1
LL = ML + MU + 1
IER = MIN0(N1,ML,MU,N1-ML,N1-MU,MDIM-LL-ML)
IF (IER .LT. 0) GO TO 110
IER = 0
IF (N .EQ. 1) GO TO 92
IF (ML .EQ. 0) GO TO 32
C INITIALIZE STORAGE. SHIFT AND ZERO OUT ELEMENTS OF B. --------------
DO 30 I = 1,ML
II = MU + I
K = ML + 1 - I
DO 10 J = 1,II
JPK = J + K
10 B(J,I) = B(JPK,I)
K = II + 1
DO 20 J = K,LL
20 B(J,I) = 0.0
30 CONTINUE
32 LR = ML
DO 90 NR = 1,N1
NP = NR + 1
IF (LR .NE. N) LR = LR + 1
C FIND THE PIVOT IN COLUMN NR, WITHIN THE BAND, AT OR BELOW ROW NR. ---
MX = NR
XM = ABS(B(1,NR))
IF (ML .EQ. 0) GO TO 42
DO 40 I = NP,LR
IF (ABS(B(1,I)) .LE. XM) GO TO 40
MX = I
XM = ABS(B(1,I))
40 CONTINUE
42 IP(NR) = MX
IF (MX .EQ. NR) GO TO 60
C INTERCHANGE ROWS NR AND MX. -----------------------------------------
DO 50 I = 1,LL
XX = B(I,NR)
B(I,NR) = B(I,MX)
50 B(I,MX) = XX
60 XM = B(1,NR)
IF (XM .EQ. 0.0) GO TO 100
C CONSTRUCT ELEMENTS OF L AND U. --------------------------------------
B(1,NR) = 1.0/XM
IF (ML .EQ. 0) GO TO 90
XM = -B(1,NR)
KK = MIN0(N-NR,LL-1)
DO 80 I = NP,LR
J = LL + I - NR
XX = B(1,I)*XM
B(J,NR) = XX
DO 70 II = 1,KK
70 B(II,I) = B(II+1,I) + XX*B(II+1,NR)
80 B(LL,I) = 0.0
90 CONTINUE
92 NR = N
IF (B(1,N) .EQ. 0.0) GO TO 100
B(1,N) = 1.0/B(1,N)
RETURN
C
100 IER = NR
RETURN
C
110 IER = -1
RETURN
C----------------------- END OF SUBROUTINE DECBR ---------------------
END
SUBROUTINE SOLBR (MDIM, N, ML, MU, B, IP, Y) SOLBR
INTEGER MDIM,N,ML,MU,IP
REAL B,Y
DIMENSION B(MDIM,N), IP(N), Y(N)
C-----------------------------------------------------------------------
C THIS SUBROUTINE COMPUTES THE SOLUTION OF THE BANDED LINEAR SYSTEM
C A*X = C , GIVEN THE LU DECOMPOSITION OF THE MATRIX A FROM DECBR.
C Y = RIGHT-HAND VECTOR C, OF LENGTH N, ON INPUT,
C = SOLUTION VECTOR X ON OUTPUT.
C ALL OTHER ARGUMENTS ARE AS DESCRIBED IN DECBR.
C ALL THE ARGUMENTS ARE INPUT ARGUMENTS.
C THE OUTPUT ARGUMENT IS Y.
C-----------------------------------------------------------------------
INTEGER I,J,KK,LL,NB,NR,N1
REAL DP,XX
IF (N .EQ. 1) GO TO 60
N1 = N - 1
LL = ML + MU + 1
C IF ML = 0, SKIP THE FIRST BACKSOLUTION PHASE. -----------------------
IF (ML .EQ. 0) GO TO 32
C APPLY PERMUTATION AND L MULTIPLIERS TO Y. ---------------------------
DO 30 NR = 1,N1
IF (IP(NR) .EQ. NR) GO TO 10
J = IP(NR)
XX = Y(NR)
Y(NR) = Y(J)
Y(J) = XX
10 KK = MIN0(N-NR,ML)
DO 20 I = 1,KK
NRPI = NR + I
LLPI = LL + I
20 Y(NRPI) = Y(NRPI) + Y(NR)*B(LLPI,NR)
30 CONTINUE
32 LL = LL - 1
C BACKSOLVE FOR X USING U. --------------------------------------------
Y(N) = Y(N)*B(1,N)
KK = 0
DO 50 NB = 1,N1
NR = N - NB
IF (KK .NE. LL) KK = KK + 1
DP = 0.0
IF (LL .EQ. 0) GO TO 50
DO 40 I = 1,KK
IP1 = I + 1
NRPI = NR + I
40 DP = DP + B(IP1,NR)*Y(NRPI)
50 Y(NR) = (Y(NR) - DP)*B(1,NR)
RETURN
60 Y(1) = Y(1)*B(1,1)
RETURN
C----------------------- END OF SUBROUTINE SOLBR ---------------------
END
SUBROUTINE PDB(NODE,T,U,PW,N0,ML,MU) PDB
DIMENSION U(1),PW(1)
C-----------------------------------------------------------------------
C A USER DEFINED ROUTINE WHICH DEFINES THE JACOBIAN IF MITER = 1.
C SINCE THIS METHOD IS NOT RECOMMENDED, A DUMMY ROUTINE IS PROVIDED.
C THE CALLING PARAMETERS ARE DEFINED IN THE COMMENTS OF SUBROUTINE
C PSETM. NOTE THAT NODE = NPDE*NX*NY, N0 = NODE, AND MU = ML. THE
C DETAILED COMMENTS ON HOW TO WRITE THIS ROUTINE ARE TO BE FOUND IN
C REFERENCE 1 (SEE PSETM).
C-----------------------------------------------------------------------
RETURN
END
.