SUBROUTINE INITAL C... C... SOLUTION OF BURGERS' EQUATION BY SPLINE APPROXIMATION ON AN C... ADAPTIVE GRID C... C... SINGLE PRECISION C... C... SUBROUTINE INITAL C... C... (1) DEFINES THE INITIAL ADAPTIVE GRID AS THE BASIC UNIFORM C... GRID VIA A CALL TO SUBROUTINE GRIDE C... C... (2) DEFINES THE INITIAL CONDITION ON THE ADAPTIVE GRID OF C... (1) BY A CALL TO FUNCTION PHI WHICH IMPLEMENTS THE EXACT C... SOLUTION TO BURGERS EQUATION C... C... COMMON/SYSTM1/ IS ACCESSED SO THE NUMBER OF ODES DEFINED BY THE C... ADAPTIVE GRID, N, CAN BE CHANGED AS THE GRID EVOLVES COMMON/SYSTM1/ T0,TF,TP,H,ERROR,N C... C... THE FOLLOWING COMMON AREAS ARE BRIEFLY EXPLAINED BELOW C... C... /T/ /Y/ /F/ STANDARD DSS/2 COMMON AREAS DIMENSIONED FOR UP C... TO 150 ODES C... C... /B/ THE SOLUTION, CB, AND ITS SECOND DERIVATIVE, C... CBXX, DEFINED ON THE BASIC GRID C... C... /SV/ INDEPENDENT VARIABLE DEFINED ON THE ADAPTIVE C... GRID C... C... /SC/ CUBIC SPLINE COEFFICIENTS (LINEAR, QUADRATIC C... AND CUBIC COEFFICIENTS RESPECTIVELY) C... C... /C/ CONSTANTS PERTAINING TO THE PROBLEM C... C... /AB1/ PARAMETERS WHICH CONTROL THE ADAPTIVE GRID C... COMMON /T/ T 1 /Y/ CA(150) 2 /F/ CT(150) 3 /B/ CB(120), CBXX(120) 4 /SV/ XA(150) 5 /SC/ CX(150), CXX(150), CXXX(150) 6 /C/ XL, XU, V, VIS, 7 SN, NB, IP, JP COMMON/AB1/ DL, XFL, XFR, NF, 1 DXF, XB(120), LXB(120), NOB(120) C... C... SET THE VELOCITY AND VISCOSITY OF BURGERS EQUATION V=1.0E+00 VIS=0.003E+00 C... C... INITIALIZE THE COUNTERS FOR PRINTED AND PLOTTED OUTPUT USED C... IN SUBROUTINE PRINT, THE RUNNING SUM OF THE NUMBER OF ODES C... USED IN THE ADAPTIVE GRID IP=0 JP=0 SN=0.0E+00 C... C... DEFINE THE LENGTH OF THE X AXIS XL=0.0E+00 XU=1.0E+00 C... C... DEFINE THE NUMBER OF BASIC GRID POINTS, NUMBER OF SUBINTERVALS C... BETWEEN THE BASIC GRID POINTS WHEN ADAPTIVE POINTS ARE ADDED NB=41 NF=2 C... C... DEFINE AN INITIAL UNIFORM GRID OF NB POINTS CALL GRIDE(NB,XL,XU,DXB,XB) C... C... THE NUMBER OF ADAPTIVE GRID POINTS, N, INITIALLY EQUALS THE C... NUMBER OF BASIC GRID POINTS, NB N=NB C... C... DEFINE THE ADAPTIVE GRID INTERVAL, DXF, FROM THE BASIC GRID C... INTERVAL, DXB DXF=DXB/FLOAT(NF) C... C... DEFINE THE THRESHOLD FOR THE SECOND DERIVATIVE, ABOVE WHICH C... ADAPTIVE GRID POINTS ARE ADDED DL=20.0E+00 C... C... DEFINE THE MAXIMUM DISTANCES TO THE LEFT AND RIGHT OF A BASIC C... GRID POINT, XFL AND XFR, BEYOND WHICH ADAPTIVE GRID POINTS WILL C... NOT BE ADDED SFL=1.0E+00 SFR=1.5E+00 XFL=SFL*V*TP XFR=SFR*V*TP C... C... DEFINE THE INITIAL BASIC AND ADAPTIVE GRIDS, EACH WITH NB POINTS DO 1 I=1,NB C... C... INITIALIZE ARRAY LXB WHICH SPECIFIES THE EXISTENCE OF ADAPTIVE C... GRID POINTS LXB(I)=0 C... C... STORE THE SUBSCRIPTS OF THE BASIC GRID IN ARRAY NOB NOB(I)=I C... C... DEFINE THE INITIAL CONDITIONS OF THE BASIC AND ADAPTIVE GRIDS C... FROM THE EXACT SOLUTION TO BURGERS EQUATION CB(I)=PHI(0.0E+00,XB(I)) CA(I)=CB(I) XA(I)=XB(I) 1 CONTINUE C... C... INITIALIZE THE DERIVATIVE CALCULATIONS IN SUBROUTINE DERV CALL DERV RETURN END SUBROUTINE DERV C... C... SUBROUTINE DERV C... C... (1) COMPUTES THE FIRST AND SECOND DERIVATIVES, UX AND UXX, C... VIA A CALL TO SUBROUTINE NCSPLE C... C... (2) USES DERIVATIVES UX AND UXX TO ASSEMBLE BURGERS EQUATION C... IN THE USUAL METHOD OF LINES FORMAT C... COMMON/SYSTM1/ T0,TF,TP,H,ERROR,N COMMON /T/ T 1 /Y/ CA(150) 2 /F/ CT(150) 3 /B/ CB(120), CBXX(120) 4 /SV/ XA(150) 5 /SC/ CX(150), CXX(150), CXXX(150) 6 /C/ XL, XU, V, VIS, 7 SN, NB, IP, JP COMMON/AB1/ DL, XFL, XFR, NF, 1 DXF, XB(120), LXB(120), NOB(120) C... C... APPLY THE BOUNDARY CONDITIONS FOR BURGERS EQUATION AT X = 0 C... AND X = 1 CA(1)=PHI(T,0.0E+00) CA(N)=PHI(T,1.0E+00) C... C... COMPUTE THE FIRST AND SECOND SPATIAL DERIVATIVES IN BURGERS C... EQUATION, UX AND UXX, BY THE NORMAL CUBIC SPLINE CALL NCSPLE(N,XA,CA,CX,CXX,CXXX) C... C... ASSEMBLE BURGERS EQUATION OVER GRID POINTS 2 TO N-1. COMMON C... /SC/ ACTUALLY CONTAINS THE COEFFICIENTS OF THE CUBIC SPLINE C... C(X) = A0 + A1*(X - XI) + A2*(X - XI)**2 + A3*(X - XI)***3, C... I.E., COMMON/SC/A1(150),A2(150),A3(150), AFTER NCSPLE IS C... CALLED. THUS CX(XI) = A1, CXX(XI) = 2*A2, CXXX(XI) = 6*A3, C... AND THE SECOND ARRAY IN /SC/ MUST BE MULTIPLIED BY 2 TO OBTAIN C... THE SECOND DERIVATIVE, CXX, IN BURGERS EQUATION. ALSO, FROM C... THE EQUATION FOR THE SPLINE, C(XI) = A0 = CA(XI) WHERE ARRAY C... CA IS IN COMMON/Y/. SUBROUTINE NCSPLE THEREFORE NUMERICALLY C... DIFFERENTIATES CA TO CX, CXX AND CXXX AS WELL AS EVALUATING C... THE COEFFICIENTS OF THE CUBIC SPLINE (THESE TWO OPERATIONS C... ARE EQUIVALENT) NM1=N-1 DO 1 I=2,NM1 CT(I)=2.0E+00*VIS*CXX(I)-CA(I)*CX(I) 1 CONTINUE C... C... ZERO THE TIME DERIVATIVES AT I = 1, N TO INSURE THE BOUNDARY C... CONDITIONS ARE OBSERVED CT(1)=0.0E+00 CT(N)=0.0E+00 RETURN END SUBROUTINE PRINT(NI,NO) C... C... SUBROUTINE PRINT C... C... (1) STORES THE NUMERICAL AND EXACT SOLUTIONS TO BURGERS C... EQUATION DURING SUCCESSIVE CALLS FOR PLOTTING C... C... (2) REDEFINES THE ADAPTIVE GRID IN ACCORDANCE WITH CHANGES C... IN THE NUMERICAL SOLUTION VIA A CALL TO SUBROUTINE ANUGB1 C... C... (3) PLOTS THE NUMBER OF POINTS IN THE ADAPTIVE GRID VS T AT C... THE END OF THE RUN AND PRINTS THE AVERAGE NUMBER OF GRID C... POINTS DURING THE RUN C... COMMON/SYSTM1/ T0,TF,TP,H,ERROR,N COMMON /T/ T 1 /Y/ CA(150) 2 /F/ CT(150) 3 /B/ CB(120), CBXX(120) 4 /SV/ XA(150) 5 /SC/ CX(150), CXX(150), CXXX(150) 6 /C/ XL, XU, V, VIS, 7 SN, NB, IP, JP COMMON/AB1/ DL, XFL, XFR, NF, 1 DXF, XB(120), LXB(120), NOB(120) C... C... DIMENSION THE ARRAYS FOR THE PLOTTED SOLUTION DIMENSION GR(150), CAA(150), PN(51), TN(51), 1 CBXXP(2,51) C... C... INCREMENT THE INDEX FOR THE PLOTTED SOLUTION IP=IP+1 C... C... COMPUTE THE RUNNING SUM OF THE NUMBER OF ODES SN=SN+FLOAT(N) C... C... STORE THE NUMBER OF ODES FOR SUBSEQUENT PLOTTING PN(IP)=FLOAT(N) TN(IP)=T C... C... CALL SUBROUTINE NCSPLE TO UPDATE THE SPLINE COEFFICIENTS USED C... IN THE NEW ADAPTIVE GRID CALL NCSPLE(N,XA,CA,CX,CXX,CXXX) C... C... STORE THE DEPENDENT VARIABLE AND ITS SECOND DERIVATIVE ON THE C... ADAPTIVE GRID, CA(I) AND CXX(I), AT THE BASIC GRID POINTS WITH C... SUBSCRIPT K IN PREPARATION FOR UPDATING THE ADAPTIVE GRID. THE C... FACTOR OF 2 MULTIPLYING CXX(K) IS EXPLAINED IN SUBROUTINE DERV DO 1 I=1,NB K=NOB(I) CB(I)=CA(K) CBXX(I)=2.0E+00*CXX(K) 1 CONTINUE C... C... STORE THE MINIMUM AND MAXIMUM VALUES OF THE SECOND DERIVATIVE, C... CBXX, FOR SUBSEQUENT PLOTTING CBMAX=ABS(CBXX(1)) CBMIN=ABS(CBXX(1)) DO 10 I=2,NB CBMAX=AMAX1(ABS(CBXX(I)),CBMAX) CBMIN=AMIN1(ABS(CBXX(I)),CBMIN) 10 CONTINUE CBXXP(1,IP)=CBMAX CBXXP(2,IP)=CBMIN C... C... STORE THE EXACT SOLUTION FOR SUBSEQUENT PLOTTING IF((IP-6)*(IP-16)*(IP-26)*(IP-36)*(IP-46).NE.0)GO TO 3 JP=JP+1 DO 2 I=1,N CAA(I)=PHI(T,XA(I)) GR(I)=-0.05E+00*FLOAT(5-JP) 2 CONTINUE C... C... PLOT THE NUMERICAL SOLUTION WITH THE ADAPTIVE GRID SUPERIMPOSED C... AT SUCCESSIVE VALUES OF THE INDEPENDENT VARIABLE T. INTEGER C... JP SPECIFIES THAT ONLY ADDITIONAL CURVES ARE PLOTTED AFTER THE C... FIRST CALL TO QIKPLT IM=-N IF(JP.GT.1)GO TO 4 C... C... ****************************************************************** C... THE FOLLOWING CODING CALLS A CONTINUOUS (CALCOMP) PLOTTER. IT C... IS NOT STANDARD FORTRAN, AND THEREFORE WILL PROBABLY HAVE TO BE C... CHANGED FOR OTHER COMPUTERS AND PLOTTERS C... CALL QIKSET(4.,0.,0.25,7.,-0.3,0.2) C... CALL QIKSAX(3,3) C... CALL QIKFORM(2,2) C... CALL QIKPLT(XA,CAA,N,3H*X*,8H*U(X,T)*,18H*BURGERS EQUATION*,1) C... CALL PLOT(-5.,1.,-3) C... CALL QLINE(XA,CA,IM,4) C... CALL QLINE(XA,GR,IM,-13) GO TO 6 4 CONTINUE C... CALL QLINE(XA,CAA,N,M0) C... CALL QLINE(XA,CA,IM,4) C... CALL QLINE(XA,GR,IM,-13) C... ****************************************************************** C... C... PRINT THE NUMERICAL SOLUTION CORRESPONDING TO THE PLOTTED SOLUTION 6 WRITE(NO,7)T,N,(XA(I),I=1,N) 7 FORMAT(1H ,//,5H T = ,E11.3,3X,5H N = ,I3,/,3H X,/,(10E11.3)) WRITE(NO,8)(CAA(I),I=1,N) 8 FORMAT(1H ,7H U ANAL,/,(10E11.3)) WRITE(NO,9)( CA(I),I=1,N) 9 FORMAT(1H ,7H U NUM ,/,(10E11.3)) C... C... CONTINUE WITH SUBROUTINE PRINT, PARTICULARLY THE REDEFINITION OF C... THE ADAPTIVE GRID 3 CONTINUE C... C... CALL SUBROUTINE ANUGB1 TO REDEFINE THE ADAPTIVE GRID CALL ANUGB1(NB,CB,CBXX,N,XA,CA) C... C... AT THE END OF THE RUN, PLOT THE NUMBER OF ODES VS T IF(IP.LT.51)RETURN CALL TPLOTS(1,IP,TN,PN) C... C... COMPUTE AND PRINT THE AVERAGE NUMBER OF ODES AN=SN/FLOAT(IP) WRITE(NO,5)AN 5 FORMAT(1H ,//, 1 7H N VS T,3X,8HAVG N = ,F6.2) C... C... PLOT THE NUMBER OF ODES, THE MINIMUM AND MAXIMUM VALUES OF THE C... SECOND DERIVATIVE VS T C... ****************************************************************** C... THE FOLLOWING CODING CALLS A CONTINUOUS (CALCOMP) PLOTTER. IT C... IS NOT STANDARD FORTRAN, AND THEREFORE WILL PROBABLY HAVE TO BE C... CHANGED FOR OTHER COMPUTERS AND PLOTTERS C... CALL QIKSAX(3,3) C... CALL QIKSET(4.,0.,0.25,7.,40.,2.5) C... CALL QIKPLT(TN,PN,IP,3H*T*,3H*N*,16H*NUMBER OF ODES*,1) C... CALL QIKSAX(3,3) C... CALL QIKSET(4.,0.,0.25,7.,0.,400.) C... CALL QIKPLT(TN,CBXXP,IP,3H*T*,6H*CBXX*,19H*SECOND DERIVATIVE*,2) C... ****************************************************************** C... RETURN END REAL FUNCTION PHI(T,X) C... C... FUNCTION PHI(T,X) COMPUTES THE EXACT SOLUTION TO BURGERS C... EQUATION. IT IS USED IN THREE PLACES C... C... (1) TO PROVIDE THE EXACT SOLUTION IN SUBROUTINE PRINT FOR C... COMPARISON WITH THE NUMERICAL SOLUTION. C... C... (2) TO PROVIDE THE INITIAL CONDITION IN SUBROUTINE INITAL. C... C... (3) TO PROVIDE THE BOUNDARY CONDITIONS IN SUBROUTINE DERV. C... C... ARGUMENT LIST C... C... T INITIAL-VALUE INDEPENDENT VARIABLE IN BURGERS EQUATION C... (INPUT) C... C... X BOUNDARY-VALUE INDEPENDENT VARIABLE IN BURGERS EQUATION C... (INPUT) C... C... THE VALUE OF THE FUNCTION (PHI) IS THE EXACT SOLUTION TO BURGERS C... EQUATION AT X AND T. C... C... TYPE REAL VARIABLES AS SINGLE PRECISION C... COMMON /C/ XL, XU, V, VIS A=(0.05E+00/VIS)*(X-0.5E+00+4.95E+00*T) B=(0.25E+00/VIS)*(X-0.5E+00+0.75E+00*T) C=( 0.5E+00/VIS)*(X-0.375E+00) EA=EXP(-A) EB=EXP(-B) EC=EXP(-C) C... C... THE FOLLOWING IF WAS ADDED FOR SHORT WORD LENGTH COMPUTERS, E.G., C... 32-BIT COMPUTERS SUCH AS THE VAX, WHICH CANNOT ACCOMMODATE WIDE C... VARIATIONS IN THE EXP FUNCTION. THE INTENTION IS TO EFFECTIVELY C... AVOID A DIVISION BY ZERO. IT IS BASED ON THE OBSERVATION THAT C... WHEN THE DIVISION BY ZERO OCCURS, EC MLT EB MLT EA MLT (WHERE C... MLT DENOTES "MUCH LESS THAN") SO THE EQUATION FOR PHI BECOMES C... C... PHI = 0.1*EA/EA = 0.1 C... IF(ABS(EA+EB+EC).LT.1.0E-20)THEN PHI=0.1E+00 ELSE PHI=(0.1E+00*EA+0.5E+00*EB+EC)/(EA+EB+EC) END IF RETURN END