应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 程序语言 » 【求助】sqrt domain error（fortran）【解决】
51/1返回列表
查看: 1006  |  回复: 6
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看
当前主题已经存档。
当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

joshzrr

至尊木虫 (著名写手)

[交流] 【求助】sqrt domain error（fortran）【解决】

程序简介： 一维喷管 无粘 maccormack格式
程序运行总是出现标题里的错误 但是build也没问题 找不到错误了 改到网格为31  喷管长度为3 就可以运行了 不知道为什么
下面是程序：
!nozzle
PARAMETER N=101
DIMENSION ROU(N),V(N),T(N),P(N),X(N),A(N),DELTT(N),H(N)
DIMENSION U1(N),U2(N),U3(N),F1(N),F2(N),F3(N)
DIMENSION DU1(N),DU2(N),DU3(N),DF1(N),DF2(N),DF3(N)
DIMENSION RU1(N),RU2(N),RU3(N),RF1(N),RF2(N),RF3(N),RROU(N),RT(N)
DIMENSION DRU1(N),DRU2(N),DRU3(N),DRF1(N),DRF2(N),DRF3(N)
DIMENSION FDU1(N),FDU2(N),FDU3(N),FU1(N),FU2(N),FU3(N)
REAL L
C=0.5
R=1.4
L=10.0
DELTX=L/(N-1)
X(1)=0.0
OPEN(1,FILE='ZCE.DAT')
DO 10 I=1,N-1
X(I+1)=X(I)+DELTX
10 CONTINUE
!---------------NOZZLE SHAPE AND INITIAL CONDITION----------------------
!A(I)=(1.0-0.661514*exp(-(alog(2.))*((X(I)-1.5)/0.6)**2))/0.338486
DO 20 I=1,N
!A(I)=1.+1.1*(X(I)-1.5)**2
A(I)=(1.0-0.661514*exp(-(alog(2.))*((X(I)-2.5)/0.6)**2))/0.338486
20 CONTINUE
DO 22 I=(N+1)/2,N
A(I)=(0.536572-0.198086*exp(-(alog(2.))*((X(I)-2.5)/0.6)**2))/0.338486
22 CONTINUE
ROU(1)=1
T(1)=1
DO 50 I=1,N
T(I)=1-0.06*X(I)
50 CONTINUE
DO 40 I=2,N-1
ROU(I)=1-0.023*X(I)
V(I)=(0.05+0.11*X(I))
T(I)=1-0.006*X(I)
40 CONTINUE
V(N)=0.4*SQRT(T(N))
ROU(N)=1/R
!H(I)=(1/A(I))*SQRT((2/(R+1))*(1+0.5*(R-1)*(H(I)**2))**((R+1)/(R-1)))
!ROU(I)=(1+0.5*(R-1)*(H(I)**2))**(1/(1-R))
!T(I)=(1+0.5*(R-1)*(H(I)**2))**(-1)
!40 CONTINUE
!DO 50 I=1,N
!V(I)=H(I)*SQRT(T(I))
!50 CONTINUE
!DO 50 I=1,N (1/A(I))*SQRT((2/(R+1))*(1+0.5*(R-1)*(H(I)**2))**((R+1)/(R-1)))
!T(I)=1/(1+0.5*(R-1)*H(I)**2)
!50 CONTINUE
!50 CONTINUE
!DO 60 I=1,N
!V(I)=0.59/(ROU(I)*A(I))
!60 CONTINUE
DO 70 I=1,N
U1(I)=ROU(I)*A(I)
U2(I)=ROU(I)*A(I)*V(I)
U3(I)=ROU(I)*A(I)*(T(I)/(R-1.)+R/2.*V(I)**2)
F1(I)=U2(I)
F2(I)=U2(I)**2/U1(I)+(R-1.)/R*(U3(I)-R/2.*U2(I)**2/U1(I))
F3(I)=R*U2(I)*U3(I)/U1(I)-R*(R-1.)/2.*U2(I)**3/(U1(I)**2)
70 CONTINUE

WRITE(1,*)"*****NOZZLE SHAPE AND THE INITIAL CONDITIONS****"
WRITE(1,21)
21        FORMAT(/13X,'X',5X,'A',7X,'ROU',5X,'V',7X,'T',5X,'H'/)
DO 31 I=1,N
31 WRITE(1,41) X(I),A(I),ROU(I),V(I),T(I),H(I)
41        FORMAT(8X,F6.2,2X,F6.3,2X,F6.3,2X,F6.3,2X,F6.3,2X,F6.3)
TIME=0.
        !___________________________________________________________
!--------------DELTT----------------
121 DO 80 I=1,N
DELTT(I)=C*DELTX/(SQRT(T(I))+V(I))
FU1(I)=U1(I)
FU2(I)=U2(I)
FU3(I)=U3(I)
F1(I)=U2(I)
F2(I)=U2(I)**2/U1(I)+(R-1.)/R*(U3(I)-R/2.*U2(I)**2/U1(I))
F3(I)=R*U2(I)*U3(I)/U1(I)-R*(R-1.)/2.*U2(I)**3/(U1(I)**2)
80 CONTINUE
DO 90 I=1,N-1
IF(DELTT(I).LT.DELTT(I+1))  THEN
DT=DELTT(I)
DELTT(I)=DELTT(I+1)
DELTT(I+1)=DT
ENDIF
90 CONTINUE
DT=DELTT(N)
!------------------PREDICTOR STEP------------------
DO 100 I=1,N-1
DU1(I)=-(F1(I+1)-F1(I))/DELTX
!DU2(I)=-(F2(I+1)-F2(I))/DELTX+ROU(I)*T(I)*(A(I+1)-A(I))/DELTX/R
DU2(I)=-(F2(I+1)-F2(I))/DELTX+(R-1)/R*(FU3(I)-R/2.*FU2(I)**2/FU1(I))*(ALOG(A(I+1))-ALOG(A(I)))/DELTX
DU3(I)=-(F3(I+1)-F3(I))/DELTX
RU1(I)=FU1(I)+DU1(I)*DT
RU2(I)=FU2(I)+DU2(I)*DT
RU3(I)=FU3(I)+DU3(I)*DT
RF1(I)=RU2(I)
RF2(I)=RU2(I)**2/RU1(I)+(R-1.)/R*(RU3(I)-R/2*RU2(I)**2/RU1(I))
RF3(I)=R*RU2(I)*RU3(I)/RU1(I)-R*(R-1.)/2.*RU2(I)**3/(RU1(I)**2)
RROU(I)=RU1(I)/A(I)
RT(I)=(R-1)*(RU3(I)/RU1(I)-R/2.*(RU2(I)/RU1(I))**2)
100 CONTINUE
!----------------CORRECTOR STEP------------------
DO 110 I=2,N
DRU1(I)=-(RF1(I)-RF1(I-1))/DELTX
!DRU2(I)=-(RF2(I)-RF2(I-1))/DELTX+1./R*RROU(I)*RT(I)*(A(I)-A(I-1))/DELTX
DRU2(I)=-(RF2(I)-RF2(I-1))/DELTX+(R-1)/R*(RU3(I)-R/2*RU2(I)**2/RU1(I))*(ALOG(A(I))-ALOG(A(I-1)))/DELTX
DRU3(I)=-(RF3(I)-RF3(I-1))/DELTX
110 CONTINUE
DO 120 I=2,N-1
FDU1(I)=0.5*(DU1(I)+DRU1(I))
FDU2(I)=0.5*(DU2(I)+DRU2(I))
FDU3(I)=0.5*(DU3(I)+DRU3(I))
U1(I)=FU1(I)+FDU1(I)*DT
U2(I)=FU2(I)+FDU2(I)*DT
U3(I)=FU3(I)+FDU3(I)*DT
120 CONTINUE
U2(1)=2.*U2(2)-U2(3)
U3(1)=U1(1)*(T(1)/(R-1)+R/2.*(U2(1)/U1(1))**2)
T(N)=2*T(N-1)-T(N-2)
!U1(N)=2*U1(N-1)-U1(N-2)
!U2(N)=2*U2(N-1)-U2(N-2)
!U3(N)=2*U3(N-1)-U3(N-2)
DO 130 I=1,N
ROU(I)=U1(I)/A(I)
V(I)=U2(I)/U1(I)
T(I)=(R-1)*(U3(I)/U1(I)-R/2*V(I)**2)
P(I)=ROU(I)*T(I)
H(I)=V(I)/SQRT(T(I))
130 CONTINUE

TIME=TIME+1
       
     !WRITE(1,*)"**********INMEDIA VALUE***************"
     ! WRITE(1,99)TIME
    !99    FORMAT(5X,"*******TIMES OF CALCULATION=",F9.1)
        !WRITE(1,85)
   !85        FORMAT(10X,'K',6X,'ROU',6X,'V',5X,'T',7X,'P',7X,'U1',7X,'U2',7X,'U3')
     !DO 200 I=1,N
        !WRITE(1,86)I,ROU(I),V(I),T(I),P(I),U1(I),U2(I),U3(I)
   !200        CONTINUE
   !86 FORMAT(7X,I4,2X,F6.3,2X,F6.3,2X,F6.3,2X,F6.3,2X,F6.3,2X,F6.3,2X,F6.3,2X,F6.3)
     
DO    140 I=1,N
      IF(ABS(FU1(I)-U1(I)).GT.1.E-5)     GOTO 121
      IF(ABS(FU2(I)-U2(I)).GT.1.E-5)     GOTO 121
      IF(ABS(FU3(I)-U3(I)).GT.1.E-5)     GOTO 121
140          CONTINUE

!         WRITE(*,*)"**********FINAL VALUE***************"
!         WRITE(1,*)"**********FINAL VALUE***************"
!     WRITE(1,82)
!          WRITE(*,82)
!82    FORMAT(10X,'I',6X,'X',7X,'A',6X,'ROU',6X,'V',6X,'T',7X,'P')
!    DO 300 I=1,N
  !   WRITE(*,91)I,X(I),A(I),ROU(I),V(I),T(I),P(I)
!          WRITE(1,91)I,X(I),A(I),ROU(I),V(I),T(I),P(I)
!91    FORMAT(7X,I4,2X,F6.2,2X,F6.3,2X,F6.3,2X,F6.3,2X,F6.3,2X,F6.3,2X,F6.3)
!300          CONTINUE
         ! WRITE(*,92)TIME
!92    FORMAT(5X,"*******TIMES OF CALCULATION=",F9.1)
WRITE(1,*)'************ROU*************'
WRITE(1,123)( ROU(I) ,I=1,N)
WRITE(1,*)'**********T*****************'
WRITE(1,123)( T(I) ,I=1,N)
WRITE(1,*)'*************P**************'

WRITE(1,123)( P(I) ,I=1,N)
WRITE(1,*)'***********H****************'

WRITE(1,123)( H(I) ,I=1,N)

123 FORMAT(1X,F9.6)
      END       

这个板块不如材料的火啊  

[ Last edited by woshilsh on 2008-12-11 at 14:12 ]
回复此楼

» 猜你喜欢

» 本主题相关商家推荐: (我也要在这里推广)
1楼 2008-12-09 10:30:37
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

joshzrr

至尊木虫 (著名写手)
收敛了 哈哈 改了c的值
赞一下 回复此楼
7楼2008-12-11 10:21:30
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 7 个回答

grant.tgb

木虫 (小有名气)

开方遇到零或小于零的数了。

呵呵，这是可以肯定的事。

查查你程序里sqrt()中的变量吧。
赞一下 回复此楼
2楼2008-12-10 09:39:23
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

grant.tgb

木虫 (小有名气)
★ ★
coldwind042(金币+2,VIP+0):谢谢提供帮助！
DO 130 I=1,N
ROU(I)=U1(I)/A(I)
V(I)=U2(I)/U1(I)
T(I)=(R-1)*(U3(I)/U1(I)-R/2*V(I)**2)
P(I)=ROU(I)*T(I)
!！！！！！！！！！！！！！ 问题出在这个T(I)=(R-1)*(U3(I)/U1(I)-R/2*V(I)**2)处，当算到time=20，在I=51处，T(i)<0了，那么肯定就是U3(I)/U1(I) < R/2*V(I)**2了！你自己查看是U3还是U1还是V的问题吧。反正就这3个变量。
H(I)=V(I)/SQRT(T(I))
130 CONTINUE
赞一下(10人) 回复此楼
3楼2008-12-10 09:47:34
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

joshzrr

至尊木虫 (著名写手)
引用回帖:
Originally posted by grant.tgb at 2008-12-10 09:47:
DO 130 I=1,N
ROU(I)=U1(I)/A(I)
V(I)=U2(I)/U1(I)
T(I)=(R-1)*(U3(I)/U1(I)-R/2*V(I)**2)
P(I)=ROU(I)*T(I)
!！！！！！！！！！！！！！ 问题出在这个T(I)=(R-1)*(U3(I)/U1(I)-R/2*V(I)**2)处，当算到time ...

首先谢谢你的帮助
你是怎么知道的错误出在这里的啊？

但是我往回查了下没感觉错啊 还有为什么改成31个节点 就行了呢？

[ Last edited by joshzrr on 2008-12-10 at 18:32 ]
赞一下 回复此楼
4楼2008-12-10 18:31:26
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
查看全部 7 个回答
普通表情 龙 兔 虎 猫 高级回复 (可上传附件)
信息提示
关闭
请填处理意见
关闭