C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C --- DRIVER FOR RADAU5 AT DISCHARGE PRESSURE CONTROL PROBLEM C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C LINK dr_DPcontrol radau5m vecteurw decsol dc_decsol OR C LINK dr_DPcontrol radau5m vecteurw lapack lapackc dc_lapack IMPLICIT DOUBLE PRECISION (A-H,O-Z) C --- PARAMETERS FOR RADAU5 (FULL JACOBIAN) PARAMETER (ND=4,LWORK=4*ND*ND+23*ND+20,LIWORK=3*ND+20) DIMENSION Y(ND),WORK(LWORK),IWORK(LIWORK) CHARACTER*20 NOMFICY,NOMFICZ EXTERNAL FEX,JEX,SOLOUT,SOLOUTZ,MAS C --- DEFINITION OF THE RESULT FILES WRITE(6,*)'NAME OF THE RESULT FILE FOR THE DIFFERENTIAL COMPONENTS?' READ(*,*) NOMFICY WRITE(6,*)'THE RESULT FILE FOR Y IS : ',NOMFICY NFICY = 10 OPEN(NFICY,FORM='FORMATTED',STATUS='UNKNOWN',FILE=NOMFICY) WRITE(6,*)'NAME OF THE RESULT FILE FOR THE ALGEBRAIC COMPONENTS?' READ(*,*) NOMFICZ WRITE(6,*)'THE RESULT FILE FOR Z IS : ',NOMFICZ NFICZ = 11 OPEN(NFICZ,FORM='FORMATTED',STATUS='UNKNOWN',FILE=NOMFICZ) C --- PARAMETER IN THE DIFFERENTIAL EQUATION IPAR=0 C --- DIMENSION OF THE SYSTEM N=4 C --- COMPUTE THE JACOBIAN ANALYTICALY IJAC=1 C --- JACOBIAN IS A FULL MATRIX MLJAC=N C --- DIFFERENTIAL EQUATION IS IN IMPLICIT FORM IMAS=1 MLMAS=0 MUMAS=0 C --- OUTPUT ROUTINE IS USED DURING INTEGRATION IOUT=1 C --- INITIAL VALUES X=0.0D0 Y(1)= 0.25D0 Y(2)= 6.8566812196182226D0 Y(3)=734.047759838158626176D0 Y(4)= 9.998240239254900579D0 C --- ENDPOINT OF INTEGRATION XEND=40.0D0 C --- REQUIRED TOLERANCE WRITE(*,*)'RTOL?' READ(*,*)RTOL ATOL=1.0D0*RTOL ITOL=0 C --- INITIAL STEP SIZE H=1.0D-7 C --- SET DEFAULT VALUES DO I=1,20 IWORK(I)=0 WORK(I)=0.D0 END DO IWORK(5)=3 IWORK(6)=1 WORK(9)=2.0D0 C --- CALL OF THE SUBROUTINE RADAU5 CALL RADAU5M(N,FEX,X,Y,XEND,H, & RTOL,ATOL,ITOL, & JEX,IJAC,MLJAC,MUJAC, & MAS,IMAS,MLMAS,MUMAS, & SOLOUT,SOLOUTZ,IOUT,NFICY,NFICZ, & WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID) C --- CLOSE THE RESULT FILES CLOSE(NFICY) CLOSE(NFICZ) C --- PRINT FINAL SOLUTION WRITE (6,99) X,Y(1),Y(2) WRITE (6,97) X,Y(3) WRITE (6,98) X,Y(4) 99 FORMAT(1X,'X =',F5.2,' Y =',2E22.14) 97 FORMAT(1X,'X =',F5.2,' Y =',E22.14) 98 FORMAT(1X,'X =',F5.2,' Z =',E22.14) C --- PRINT STATISTICS WRITE (6,90) RTOL 90 FORMAT(' RTOL=',D8.2) WRITE (6,91) (IWORK(J),J=14,20) 91 FORMAT(' FCN=',I5,' JAC=',I4,' STEP=',I4,' ACCPT=',I4, & ' REJCT=',I3,' DEC=',I4,' SOL=',I5) STOP END C C SUBROUTINE SOLOUT (NFICY,NFICZ,NR,XOLD,X,Y,CONT,LRC,N, & RPAR,IPAR,IRTRN) C --- PRINTS SOLUTION AT EQUIDISTANT OUTPUT-POINTS BY USING "CONTR5" IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION Y(N),CONT(LRC),RPAR(*) COMMON /INTERN/XOUTY,XOUTZ IF (NR.EQ.1) THEN WRITE (NFICY,99) X,Y(1),Y(2),Y(3),NR-1 WRITE (NFICZ,98) X,Y(4),NR-1 XOUTY=0.5D0 ELSE 10 CONTINUE IF (X.GE.XOUTY) THEN C --- CONTINUOUS OUTPUT FOR RADAU5M WRITE (NFICY,99) XOUTY,CONTR5(1,XOUTY,CONT,LRC), & CONTR5(2,XOUTY,CONT,LRC), & CONTR5(3,XOUTY,CONT,LRC),NR-1 WRITE (NFICZ,98) XOUTY,CONTR5(4,XOUTY,CONT,LRC),NR-1 XOUTY=XOUTY+0.5D0 GOTO 10 END IF END IF 99 FORMAT(1X,F5.2,3E22.14,I4) 98 FORMAT(1X,F5.2,E22.14,I4) RETURN END C C SUBROUTINE SOLOUTZ (NFICY,NFICZ,NR,XOLD,X,Y,N,VALG,LRV, & RPAR,IPAR,IRTRN,ICAS,R1,R2,HTOT,H2,H3,LAST) C --- PRINTS SOLUTION AT EQUIDISTANT OUTPUT-POINTS BY USING "VECTWY" AND "VECTWZ" IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION Y(N),RPAR(*) DIMENSION VWT(9),WY(6),VALG(LRV) COMMON/INTERN/XOUTY,XOUTZ EXTERNAL VECTWZ,VECTWY IF (NR.EQ.1) THEN WRITE(NFICY,99) X,Y(1),Y(2),Y(3),NR-1 WRITE(NFICZ,98) X,Y(4),NR-1 XOUTZ=0.5D0 XOUTY=0.5D0 ELSE IF (LAST) THEN XPAR=X ELSE XPAR=XOLD ENDIF 51 CONTINUE IF (XPAR.GE.XOUTZ) THEN T=(XOUTZ+HTOT-X)/HTOT CALL VECTWZ(ICAS,VWT,R1,R2,T) W1=VWT(1) W2=VWT(2) W3=VWT(3) W4=VWT(4) W5=VWT(5) W6=VWT(6) W7=VWT(7) W8=VWT(8) W9=VWT(9) YN4=W1*VALG(4)+W2*VALG(16)+W3*VALG(28)+ + W4*VALG(8)+W5*VALG(20)+W6*VALG(32)+ + W7*VALG(12)+W8*VALG(24)+W9*VALG(36) WRITE(NFICZ,98) XOUTZ,YN4,NR-1 XOUTZ=XOUTZ+0.5D0 GOTO 51 ENDIF 52 CONTINUE IF (X.GE.XOUTY) THEN IF ((XOLD-H2).GE.XOUTY) THEN HTOTY=HTOT-H3 T=(XOUTY-XOLD+HTOTY)/HTOTY R=(HTOT-H3-H2)/HTOTY CALL VECTWY(WY,R,T) WY1=WY(1) WY2=WY(2) WY3=WY(3) WY4=WY(4) WY5=WY(5) WY6=WY(6) YN1=WY1*VALG(1)+WY2*VALG(13)+WY3*VALG(25) YN1=YN1+WY4*VALG(5)+WY5*VALG(17)+WY6*VALG(29) YN2=WY1*VALG(2)+WY2*VALG(14)+WY3*VALG(26) YN2=YN2+WY4*VALG(6)+WY5*VALG(18)+WY6*VALG(30) YN3=WY1*VALG(3)+WY2*VALG(15)+WY3*VALG(27) YN3=YN3+WY4*VALG(7)+WY5*VALG(19)+WY6*VALG(31) ELSE HTOTY=H2+H3 T=(XOUTY-X+HTOTY)/HTOTY R=H2/HTOTY CALL VECTWY(WY,R,T) WY1=WY(1) WY2=WY(2) WY3=WY(3) WY4=WY(4) WY5=WY(5) WY6=WY(6) YN1=WY1*VALG(5)+WY2*VALG(17)+WY3*VALG(29) YN1=YN1+WY4*VALG(9)+WY5*VALG(21)+WY6*VALG(33) YN2=WY1*VALG(6)+WY2*VALG(18)+WY3*VALG(30) YN2=YN2+WY4*VALG(10)+WY5*VALG(22)+WY6*VALG(34) YN3=WY1*VALG(7)+WY2*VALG(19)+WY3*VALG(31) YN3=YN3+WY4*VALG(11)+WY5*VALG(23)+WY6*VALG(35) ENDIF WRITE(NFICY,99) XOUTY,YN1,YN2,YN3,NR-1 XOUTY=XOUTY+0.5D0 GOTO 52 ENDIF ENDIF RETURN 99 FORMAT(1X,F5.2,3E22.14,I4) 98 FORMAT(1X,F5.2,E22.14,I4) END C C SUBROUTINE FEX(N,X,Y,DY,RPAR,IPAR) C --- RIGHT-HAND SIDE OF DISCHARGE PRESSURE CONTROL PROBLEM IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION MU DIMENSION Y(N),DY(N),RPAR(*) MU=15.0D0+5.0D0*DTANH(X-10.0D0) AUX=Y(3)*5.0D-02*(3.35D0-0.075D0*Y(4)+1.0D-03*Y(4)*Y(4)) DY(1)=5.0D-02*(Y(2)-Y(1))-AUX/3.0D+02 DY(2)=-1.0D0/75.0D0*(AUX-99.1D0) DY(3)=MU-Y(4) DY(4)=Y(3)*Y(3)*2.5D-03-2458.1764+MU*MU/(1.44D0*Y(1)*Y(1)) RETURN END C C SUBROUTINE JEX(N,X,Y,DFY,LDFY,RPAR,IPAR) C --- JACOBIAN OF DISCHARGE PRESSURE CONTROL PROBLEM IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION MU DIMENSION Y(N),DFY(LDFY,N),RPAR(*) DO I=1,LDFY DO J=1,N DFY(I,J)=0.0D0 END DO END DO AUX=DTANH(X-10.0D0) MU=15.0D0+5.0D0*AUX DMU=5.0D0*(1-AUX*AUX) AUX1=-1.0D0/1500.0D0*(3.35D0-0.075D0*Y(4)+1.0D-3*Y(4)*Y(4)) AUX2=-Y(3)/1500.0D0*(2.0D-3*Y(4)-0.075D0) AUX3=MU/(0.72D0*Y(1)*Y(1)) DFY(1,1)=-5.0D-2 DFY(1,2)=5.0D-2 DFY(1,3)=0.25D0*AUX1 DFY(1,4)=0.25D0*AUX2 DFY(2,3)=AUX1 DFY(2,4)=AUX2 DFY(3,4)=-1.0D0 DFY(4,1)=-AUX3*MU/Y(1) DFY(4,3)=5.0D-03*Y(3) RETURN END C C SUBROUTINE MAS(N, AM, LMAS,RPAR,IPAR) C --- MASS MATRIX IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION AM(LMAS,N),RPAR(*) AM(1,1) = 1.0D0 AM(1,2) = 1.0D0 AM(1,3) = 1.0D0 AM(1,4) = 0.0D0 RETURN END