SUBROUTINE INITAL C... C... H2 - HEAT CONDUCTION IN ONE DIMENSIONAL SOLIDS C... C... DOUBLE PRECISION CODING IS USED IMPLICIT DOUBLE PRECISION (A-H,O-Z) C... C... COMMON AREA WITH THE COORDINATE SYSTEM COMMON/I/ NC C... C... SELECT COORDINATE SYSTEM C... C... NC = 0, CARTESIAN COODINATES C... C... NC = 1, CYLINDRICAL COORDINATES C... C... NC = 2, SPHERICAL COORDINATES C... NC=0 C... C... INITIAL CONDITIONS IF(NC.EQ.0)CALL INIT10 IF(NC.EQ.1)CALL INIT11 IF(NC.EQ.2)CALL INIT12 RETURN END SUBROUTINE DERV COMMON/I/ NC C... C... CALCULATION OF TEMPORAL DERIVATIVES IF(NC.EQ.0)CALL DERV10 IF(NC.EQ.1)CALL DERV11 IF(NC.EQ.2)CALL DERV12 RETURN END SUBROUTINE PRINT(NI,NO) COMMON/I/ NC C... C... WRITING AND PLOTTING OF NUMERICAL SOLUTIONS IF(NC.EQ.0)CALL PRT10(NI,NO) IF(NC.EQ.1)CALL PRT11(NI,NO) IF(NC.EQ.2)CALL PRT12(NI,NO) RETURN END SUBROUTINE INIT10 C... C... INITIAL CONDITIONS FOR CARTESIAN COORDINATES C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NX=21) COMMON/T/ TIME, NSTOP, NORUN + /Y/ T(NX) + /F/ TT(NX) + /S/ TX(NX), TXX(NX) + /P/ HX0K + /I/ NC C... C... SELECT BIOT NUMBER, H*X0/K IF(NORUN.EQ.1)HX0K=0.0D0 IF(NORUN.EQ.2)HX0K=1.0D0/5.0D0 IF(NORUN.EQ.3)HX0K=1.0D0/2.0D0 IF(NORUN.EQ.4)HX0K=1.0D0/1.0D0 IF(NORUN.EQ.5)HX0K=1.0D0/0.5D0 IF(NORUN.EQ.6)HX0K=1.0D0/0.05D0 C... C... INITIAL CONDITION DO 1 I=1,NX T(I)=1.0D0 1 CONTINUE RETURN END SUBROUTINE DERV10 C... C... TEMPORAL DERIVATIVES FOR CARTESIAN COORDINATES C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NX=21) COMMON/T/ TIME, NSTOP, NORUN + /Y/ T(NX) + /F/ TT(NX) + /S/ TX(NX), TXX(NX) + /P/ HX0K + /I/ NC C... C... BOUNDARY CONDITION AT X = 0 NL=2 TX(1)=0.0D0 C... C... BOUNDARY CONDITION AT X = X0 IF(NORUN.NE.6)THEN NU=2 TX(NX)=-HX0K*T(NX) ELSE IF(NORUN.EQ.6)THEN NU=1 T(NX)=0.0D0 END IF C... C... TXX CALL DSS044(0.0D0,1.0D0,NX,T,TX,TXX,NL,NU) C... C... FOURIER'S SECOND LAW DO 1 I=1,NX TT(I)=TXX(I) 1 CONTINUE RETURN END SUBROUTINE PRT10(NI,NO) C... C... WRITING AND PLOTTING SOLUTION FOR CARTESIAN COORDINATES C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NX=21) COMMON/T/ TIME, NSTOP, NORUN + /Y/ T(NX) + /F/ TT(NX) + /S/ TX(NX), TXX(NX) + /P/ HX0K + /I/ NC C... C... MONITOR OUTPUT IF(TIME.LT.0.001D0)WRITE(*,2)HX0K 2 FORMAT(//,' H*X0/K = ',F8.3,' FROM PRT10') WRITE(*,1)TIME,T(1),T(NX) C... C... OPEN A FILE FOR MATLAB PLOTTING IF((NORUN.EQ.1).AND.(TIME.LT.0.001D0))THEN OPEN(8,FILE='h2.out') END IF C... C... HEADING FOR THE NUMERICAL SOLUTION IF(TIME.LT.0.001D0)WRITE(NO,3) 3 FORMAT(8X,'t ',3X,' x = 0 ',3X,'x = 0.2',3X,'x = 0.4', + 3X,'x = 0.6',3X,'x = 0.8',3X,' x = 1 ') C... C... WRITE NUMERICAL SOLUTION FOR 21 GRID POINTS NXI=4 C... C... WRITE SOLUTION ON FILE OUTPUT WRITE(NO,1)TIME,(T(I),I=1,NX,NXI) C... C... WRITE SOLUTION ON FILE H2.OUT WRITE( 8,1)TIME,(T(I),I=1,NX,NXI) 1 FORMAT(F10.2,6F10.5) RETURN END SUBROUTINE INIT11 C... C... INITIAL CONDITIONS FOR CYLINDRICAL COORDINATES C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NR=21) COMMON/T/ TIME, NSTOP, NORUN + /Y/ T(NR) + /F/ TT(NR) + /S/ TR(NR), TRR(NR) + /P/ HR0K + /I/ NC C... C... SELECT BIOT NUMBER, H*R0/K IF(NORUN.EQ.1)HR0K=0.0D0 IF(NORUN.EQ.2)HR0K=1.0D0/5.0D0 IF(NORUN.EQ.3)HR0K=1.0D0/2.0D0 IF(NORUN.EQ.4)HR0K=1.0D0/1.0D0 IF(NORUN.EQ.5)HR0K=1.0D0/0.5D0 IF(NORUN.EQ.6)HR0K=1.0D0/0.05D0 C... C... INITIAL CONDITION DO 1 I=1,NR T(I)=1.0D0 1 CONTINUE RETURN END SUBROUTINE DERV11 C... C... TEMPORAL DERIVATIVES FOR CYLINDRICAL COORDINATES C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NR=21) COMMON/T/ TIME, NSTOP, NORUN + /Y/ T(NR) + /F/ TT(NR) + /S/ TR(NR), TRR(NR) + /P/ HR0K + /I/ NC C... C... DIRICHLET BOUNDARY CONDITION AT R = R0 IF(NORUN.EQ.6)THEN NU=1 T(NR)=0.0D0 END IF C... C... TR CALL DSS004(0.0,1.0D0,NR,T,TR) C... C... BOUNDARY CONDITION AT R = 0 NL=2 TR(1)=0.0D0 C... C... NEUMANN BOUNDARY CONDITION AT R = R0 IF(NORUN.NE.6)THEN NU=2 TR(NR)=-HR0K*T(NR) END IF C... C... TRR CALL DSS044(0.0D0,1.0D0,NR,T,TR,TRR,NL,NU) C... C... FOURIER'S SECOND LAW DO 1 I=1,NR IF(I.EQ.1)THEN TT(I)=2.0D0*TRR(I) ELSE R=DFLOAT(I-1)/DFLOAT(NR-1) TT(I)=TRR(I)+(1.0D0/R)*TR(I) END IF 1 CONTINUE RETURN END SUBROUTINE PRT11(NI,NO) C... C... WRITING AND PLOTTING SOLUTION FOR CYLINDRICAL COORDINATES C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NR=21) COMMON/T/ TIME, NSTOP, NORUN + /Y/ T(NR) + /F/ TT(NR) + /S/ TR(NR), TRR(NR) + /P/ HR0K + /I/ NC C... C... MONITOR OUTPUT IF(TIME.LT.0.001D0)WRITE(*,2)HR0K 2 FORMAT(//,' H*R0/K = ',F8.3,' FROM PRT11') WRITE(*,1)TIME,T(1),T(NR) C... C... OPEN A FILE FOR MATLAB PLOTTING IF((NORUN.EQ.1).AND.(TIME.LT.0.001D0))THEN OPEN(8,FILE='h2.out') END IF C... C... HEADING FOR THE NUMERICAL SOLUTION IF(TIME.LT.0.001D0)WRITE(NO,3) 3 FORMAT(8X,'t ',3X,' r = 0 ',3X,'r = 0.2',3X,'r = 0.4', + 3X,'r = 0.6',3X,'r = 0.8',3X,' r = 1 ') C... C... WRITE NUMERICAL SOLUTION FOR 21 GRID POINTS NRI=4 C... C... WRITE SOLUTION ON FILE OUTPUT WRITE(NO,1)TIME,(T(I),I=1,NR,NRI) C... C... WRITE SOLUTION ON FILE H2.OUT WRITE( 8,1)TIME,(T(I),I=1,NR,NRI) 1 FORMAT(F10.2,6F10.5) RETURN END SUBROUTINE INIT12 C... C... INITIAL CONDITIONS FOR SPHERICAL COORDINATES C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NR=21) COMMON/T/ TIME, NSTOP, NORUN + /Y/ T(NR) + /F/ TT(NR) + /S/ TR(NR), TRR(NR) + /P/ HR0K + /I/ NC C... C... SELECT BIOT NUMBER, H*R0/K IF(NORUN.EQ.1)HR0K=0.0D0 IF(NORUN.EQ.2)HR0K=1.0D0/5.0D0 IF(NORUN.EQ.3)HR0K=1.0D0/2.0D0 IF(NORUN.EQ.4)HR0K=1.0D0/1.0D0 IF(NORUN.EQ.5)HR0K=1.0D0/0.5D0 IF(NORUN.EQ.6)HR0K=1.0D0/0.05D0 C... C... INITIAL CONDITION DO 1 I=1,NR T(I)=1.0D0 1 CONTINUE RETURN END SUBROUTINE DERV12 C... C... TEMPORAL DERIVATIVES FOR SPHERICAL COORDINATES C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NR=21) COMMON/T/ TIME, NSTOP, NORUN + /Y/ T(NR) + /F/ TT(NR) + /S/ TR(NR), TRR(NR) + /P/ HR0K + /I/ NC C... C... DIRICHLET BOUNDARY CONDITION AT R = R0 IF(NORUN.EQ.6)THEN NU=1 T(NR)=0.0D0 END IF C... C... TR CALL DSS004(0.0,1.0D0,NR,T,TR) C... C... BOUNDARY CONDITION AT R = 0 NL=2 TR(1)=0.0D0 C... C... NEUMANN BOUNDARY CONDITION AT R = R0 IF(NORUN.NE.6)THEN NU=2 TR(NR)=-HR0K*T(NR) END IF C... C... TRR CALL DSS044(0.0D0,1.0D0,NR,T,TR,TRR,NL,NU) C... C... FOURIER'S SECOND LAW DO 1 I=1,NR IF(I.EQ.1)THEN TT(I)=3.0D0*TRR(I) ELSE R=DFLOAT(I-1)/DFLOAT(NR-1) TT(I)=TRR(I)+(2.0D0/R)*TR(I) END IF 1 CONTINUE RETURN END SUBROUTINE PRT12(NI,NO) C... C... WRITING AND PLOTTING SOLUTION FOR SPHERICAL COORDINATES C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NR=21) COMMON/T/ TIME, NSTOP, NORUN + /Y/ T(NR) + /F/ TT(NR) + /S/ TR(NR), TRR(NR) + /P/ HR0K + /I/ NC C... C... MONITOR OUTPUT IF(TIME.LT.0.001D0)WRITE(*,2)HR0K 2 FORMAT(//,' H*R0/K = ',F8.3,' FROM PRT12') WRITE(*,1)TIME,T(1),T(NR) C... C... OPEN A FILE FOR MATLAB PLOTTING IF((NORUN.EQ.1).AND.(TIME.LT.0.001D0))THEN OPEN(8,FILE='h2.out') END IF C... C... HEADING FOR THE NUMERICAL SOLUTION IF(TIME.LT.0.001D0)WRITE(NO,3) 3 FORMAT(8X,'t ',3X,' r = 0 ',3X,'r = 0.2',3X,'r = 0.4', + 3X,'r = 0.6',3X,'r = 0.8',3X,' r = 1 ') C... C... WRITE NUMERICAL SOLUTION FOR 21 GRID POINTS NRI=4 C... C... WRITE SOLUTION ON FILE OUTPUT WRITE(NO,1)TIME,(T(I),I=1,NR,NRI) C... C... WRITE SOLUTION ON FILE H2.OUT WRITE( 8,1)TIME,(T(I),I=1,NR,NRI) 1 FORMAT(F10.2,6F10.5) RETURN END