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