SUBROUTINE INITAL IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NX=201) COMMON/T/ T, NFIN, NORUN 1 /Y/ U(NX) 2 /F/ UT(NX) 3 /S/ UX(NX), UXX(NX), X(NX) 4 /C/ VIS 5 /A/ A1(NX), A2(NX), A3(NX) C... C... PROBLEM PARAMETERS VIS=0.003D+00 C... C... SPATIAL GRID AND THE INITIAL CONDITION DO 1 I=1,NX X(I)=DFLOAT(I-1)/DFLOAT(NX-1) U(I)=PHI(X(I),0.D+00) 1 CONTINUE C... C... INITIAL DERIVATIVES CALL DERV RETURN END SUBROUTINE DERV IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NX=201) COMMON/T/ T, NFIN, NORUN 1 /Y/ U(NX) 2 /F/ UT(NX) 3 /S/ UX(NX), UXX(NX), X(NX) 4 /C/ VIS 5 /A/ A1(NX), A2(NX), A3(NX) C... C... BOUNDARY CONDITIONS U(1) =PHI(X(1) ,T) U(NX)=PHI(X(NX),T) UT(1) =0.D+00 UT(NX)=0.D+00 C... C... UX CALL DSS038(NX,X,U ,UX ,A1,A2,A3) C... C... UXX CALL DSS038(NX,X,UX,UXX,A1,A2,A3) C... C... BURGERS EQUATION DO 2 I=2,NX-1 UT(I)=VIS*UXX(I)-U(I)*UX(I) 2 CONTINUE RETURN END DOUBLE PRECISION FUNCTION PHI(X,T) C... C... FUNCTION PHI(X,T) COMPUTES THE EXACT SOLUTION FOR COMPARISON C... WITH THE NUMERICAL SOLUTION. IT IS ALSO USED TO DEFINE THE C... INITIAL AND BOUNDARY CONDITIONS FOR THE NUMERICAL SOLUTION C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/C/ VIS C... C... ANALYTICAL SOLUTION TO BURGERS EQUATION A=(0.05D+00/VIS)*(X-0.5D+00+4.95D+00*T) B=(0.25D+00/VIS)*(X-0.5D+00+0.75D+00*T) C=( 0.5D+00/VIS)*(X-0.375D+00) IF(T.LT.0.9D+00)THEN EA=DEXP(-A) EB=DEXP(-B) EC=DEXP(-C) PHI=(0.1D+00*EA+0.5D+00*EB+EC)/(EA+EB+EC) ELSE ECMA=DEXP(C-A) ECMB=DEXP(C-B) PHI=(0.1D+00*ECMA+0.5D+00*ECMB+1.0D+00)/(ECMA+ECMB+1.0D+00) END IF RETURN END SUBROUTINE PRINT(NI,NO) IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NX=201) COMMON/T/ T, NFIN, NORUN 1 /Y/ U(NX) 2 /F/ UT(NX) 3 /S/ UX(NX), UXX(NX), X(NX) 4 /C/ VIS 5 /A/ A1(NX), A2(NX), A3(NX) 6 /P/ UN(6,NX), IP C... C... INITIALIZE A COUNTER FOR THE PLOTTING DATA IP/0/ C... C... MONITOR THE SOLUTION IP=IP+1 WRITE(*,*)IP,T C... C... STORE THE NUMERICAL SOLUTION FOR SUBSEQUENT PLOTTING VIA SUBROU- C... TINE PLOT1 DO 10 I=1,NX UN(IP,I)=U(I) 10 CONTINUE C... C... PRINT THE NUMERICAL AND ANALYTICAL SOLUTIONS, THE DIFFERENCE C.. BETWEEN THE SOLUTIONS AND THE TERMS IN BURGERS EQUATION WRITE(NO,2)T 2 FORMAT(1H ,//,5H T = ,F5.2,//, 1 9X,1HX,5X,5HU NUM,4X,6HU ANAL,5X,5HERROR,3X,7HVIS*UXX,5X,5H-U*UX, 2 8X,2HUT) DO 3 I=1,NX UANAL=PHI(X(I),T) ERROR=U(I)-UANAL DIFF=VIS*UXX(I) CONV=-U(I)*UX(I) WRITE(NO,4)X(I),U(I),UANAL,ERROR,DIFF,CONV,UT(I) 4 FORMAT(F10.3,6F10.5) 3 CONTINUE C... C... PLOT THE SOLUTION CALL PLOT1 C... C... REINITIALIZE THE COUNTER FOR PLOTTING IF(IP.EQ.6)IP=0 RETURN END SUBROUTINE PLOT1 IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NX=201) COMMON/T/ T, NFIN, NORUN 1 /Y/ U(NX) 2 /F/ UT(NX) 3 /S/ UX(NX), UXX(NX), X(NX) 4 /C/ VIS 5 /A/ A1(NX), A2(NX), A3(NX) 6 /P/ UN(6,NX), IP C... C... WRITE A FILE FOR TOP DRAWER PLOTTING IF(IP.EQ.1)THEN OPEN(2,FILE='T.TOP',STATUS='NEW') WRITE(2,12) 12 FORMAT(' SET LIMITS X FROM 0 TO 1 Y FROM 0 TO 1.25',/, 1 ' SET FONT DUPLEX') END IF C... C... WRITE THE ANALYTICAL SOLUTION DO 16 I=1,NX UANAL=PHI(X(I),T) WRITE(2,11)X(I),UANAL 11 FORMAT(F10.3,F10.5) 16 CONTINUE WRITE(2,13) 13 FORMAT(' JOIN 1') C... C... WRITE THE NUMERICAL SOLUTION AND PLOT LABELS IF(IP.EQ.6)THEN WRITE(2,1) 1 FORMAT(' SET SYMBOL 0') WRITE(2,11)(X(I),UN(1,I),I=1,NX) WRITE(2,2) 2 FORMAT(' SET SYMBOL 2') WRITE(2,11)(X(I),UN(2,I),I=1,NX) WRITE(2,3) 3 FORMAT(' SET SYMBOL 4') WRITE(2,11)(X(I),UN(3,I),I=1,NX) WRITE(2,4) 4 FORMAT(' SET SYMBOL 6') WRITE(2,11)(X(I),UN(4,I),I=1,NX) WRITE(2,5) 5 FORMAT(' SET SYMBOL 8') WRITE(2,11)(X(I),UN(5,I),I=1,NX) WRITE(2,6) 6 FORMAT(' SET SYMBOL 1') WRITE(2,11)(X(I),UN(6,I),I=1,NX) WRITE(2,14)NORUN 14 FORMAT( +' SET FONT DUPLEX' +,/,' TITLE 5.0 9.5', + ' "Figure ',I2,'a: Burgers Equation"' +,/,' TITLE LEFT " u(x,t)"' +,/,' TITLE BOTTOM "x"' +,/,' TITLE 6.25 0.75 "u(x,t) vs x"' +,/,' TITLE 4.25 0.30 "t = 0, 0.2,..., 1, NX = 201', + ' Uniform Grid"') C... C... NEXT PLOT WRITE(2,15) 15 FORMAT('NEW FRAME') END IF RETURN END