Fortran Coder

查看: 27491|回复: 9

[求助] 算出来的结果全是NAN,也没有发现被除数是0,对数是负的情...

[复制链接]

20

帖子

8

主题

0

精华

熟手

F 币
126 元
贡献
73 点
发表于 2021-2-4 17:49:24 | 显示全部楼层 |阅读模式
本帖最后由 xiao. 于 2021-2-4 17:58 编辑

不知道为啥那个变量的值一直是NAN,打断点找原因也找不出来,求助大家!
大噶们帮忙试一试,谢谢啦,再次谢谢大家!

ipl_r 变量在第96行左右!


[Fortran] 纯文本查看 复制代码
!***********************参数设置*************************************************************
      implicit real*8(i-n)   
      real*8 x,a,b,answer                   
      real*8 B_NUM(2966296),AB_CO(2966296)     !波数,吸收系数
      real*8 fk(10000)                         !k的概率分布函数
      real*8 fk_0                              !记录吸收系数的对数小于-22的概率分布函数
      real*8 gk(10000)                         !k的累积概率分布函数
      real*8 fki(10000)                        !kIg的概率分布函数
      real*8 gki(10000)                        !kIg的累积概率分布函数 
      REAL*8 Ig                                 !g空间上的光谱辐射强度Ig=Ik/f(T,k),Ik=∫Iηδ(k-kη)dη,cf(T,k)=1/Ib(∫Ibηδ(k-kη)dη)
      REAL*8 PL_IB                             !黑体总辐射强度
      REAL*8 PL_IBG(2966296)                   !黑体光谱辐射强度
      REAL*8,PARAMETER:: PI=3.141592653589793
      REAL*8,PARAMETER:: C1=3.742*10.0**(-8) !普朗克第一常数,(W·cm^4/m^2)
      REAL*8,PARAMETER:: C2=1.4388           !普朗克第二常数,(cm·K)
      REAL*8,PARAMETER:: E=2.718281828       !自然常数
      REAL*8,PARAMETER:: tg=1800.0           !气体温度,K
      REAL*8 DELTAK                           !吸收系数间隔
      real*8 maxk                             !吸收系数最大值
      real*8 gaussh                           !剖分步长
      real*8 sum_gk                           ! 利用累计k分布直接计算k的积分值
      real*8 sum_lbl                          !利用lbl直接计算k的积分值
      real*8 sum_gauss                        !利用高斯求积计算k的积分值

      real*8 ipl_r(21)                        !第i个辐射源项到第j个辐射热流点的辐射热流,从左到右
      real*8 ipr_l(21)                        !第i个辐射源项到第j个辐射热流点的辐射热流,从右到左
      real*8 delta_q(20)                      !辐射源项
      real*8 delta_x                          !控制体的长度,delta_x/划分区间数,cm
      real*8 distance                         !总长度,cm
      real*8 temp_sum_l(21),temp_sum_r(21),temp(21) 
      character(10) line_distance
!**************************************************************
!********外部声明函数******************************************
      external expi3                          !3阶指数积分函数
!**************************************************************
!*****将数据读取到B_NUM和AB_CO中,并将吸收吸收转为以10为底的对数        
      OPEN(11,FILE='T_1000_P_1.0_CO2Mol_0.1_H2OMol_0.2_COMol_0.03_N2mol
     &_0.67_soot.data')
      do I=1,2966296
      READ(11,*) B_NUM(I),AB_CO(I)
      AB_CO(I)=LOG10(AB_CO(I))
      enddo
      CLOSE(11) 
!*******************计算吸收系数最大值********************************      
      maxk=0
      do i=1,2966295
      if(AB_CO(I+1).gt.AB_CO(I).and.AB_CO(I+1).gt.maxk) then
            maxk=AB_CO(I+1)
      endif
      enddo

!*******************计算吸收系数间隔**********************************
      DELTAK=(maxk+22.0)/10000.0
!**********黑体总辐射强度,温度为2000K*********************************
      PL_IB=(5.67*(10.0**-8.0)*tg**4.0)/PI
!**********黑体光谱辐射强度,单位为(W/(m^2·cm·Sr))*******************
      PL_IBG(1:2966296)=0
      DO I=1,2966296
      PL_IBG(I)=(B_NUM(I)**3.0)*C1/(E**(C2*B_NUM(I)/tg)-1.0)/PI
      enddo   
!**********fk*dk,记录△log10(k)内的点数******************************
      fk(1:10000)=0   !数组初始化为0
      fk_0=0           !初始化为0
      DO I=1,2966296
      if(AB_CO(I).LE.-22.0) then
      fk_0=fk_0+0.01*PL_IBG(I)/PL_IB
      endif
        do J=1,10000
        IF(AB_CO(I).GT.(-22.0+(J-1)*DELTAK)
     &  .AND.AB_CO(I).LT.(-22.0+J*DELTAK) ) THEN
        fk(J)=fk(J)+0.01*PL_IBG(I)/PL_IB
        ENDIF
        enddo
      enddo
!**************************************************************************************************      
!***********************1维气体辐射传递计算,不考虑壁面辐射*****************************************
!**************************************************************************************************    
      ipl_r(1:21)=0                             !第i个辐射热流点的辐射热流,从左到右
      ipr_l(1:21)=0                             !第i个辐射热流点的辐射热流,从右到左
      delta_q(1:20)=0                           !辐射源项     
      distance=0.1                              !总长度,cm
      write(line_distance,"(f10.1)") distance
      delta_x=distance/20.0                     !控制体的长度,delta_x/划分区间数,0.5cm
      temp_sum_l(1:21)=0
      temp_sum_r(1:21)=0
      temp(1:21)=0 
!************************************************************************** 
!                           lbl法计算辐射热流项
!**************************************************************************
      do i=2,21                                 !从左到右计算辐射热流点的值     
         do k=1,2966296
         temp_sum_l(1)=10.0**(AB_CO(k))*delta_x*(i-1)
       
          
        
         ipl_r(i)=ipl_r(i)+PL_IBG(k)*
     &   (1-2.0*expi3(temp_sum_l(1)))*0.01
         
           
         enddo
       
       
         
      enddo
      
      do i=20,1,-1                              !从右到左计算辐射热流点的值
         do k=1,2966296
         temp_sum_r(1)=10.0**(AB_CO(k))*(21-i)*delta_x
         ipr_l(i)=ipr_l(i)+PL_IBG(k)*
     &   (1-2.0*expi3(temp_sum_r(1)))*0.01
         enddo
      enddo 

      
      do i=1,21                 !计算净辐射热流
        temp(i)=pi*(ipl_r(i)-ipr_l(i))
      enddo      
      
      do i=1,20                 !计算辐射源项
        delta_q(i)=(temp(i+1)-temp(i))/delta_x/0.01
      enddo
      
      open(31,file='插值节点辐射热源'//trim(adjustl(line_distance))//'
     &T_1800_P_25.0_H2OMol_0.1_CO2Mol_0.1_N2mol_0.8-kdist.txt') 
      do i=1,20
        write(31,*) delta_q(i)
      enddo    
      close(31)
      
      open(32,file='插值节点净辐射热流'//trim(adjustl(line_distance))//'
     &T_1800_P_25.0_H2OMol_0.1_CO2Mol_0.1_N2mol_0.8-kdist.txt') 
      do i=1,21
        write(32,*) temp(i)
      enddo   
      close(32)

      
 
      
      END

!**************************************************
!**********3阶指数积分计算函数expi3(x)*************
!**************************************************
      function expi3(x) 
      real*8 x,aux,y,expi3
      if(x.lt.1.d0) goto 1
      y=1.d0/x
      aux=1.d0-y*(((y+3.377358d0)*y+2.052156d0)*y+2.709479d-1)/
     &((((y*1.072553d0+5.716943d0)*y+6.945239d0)*y+2.593888d0)*y
     &+2.709496d-1)
      expi3=aux*y*dexp(-x)
      goto 3
1     if(x.le.-3.d0) goto 2
      aux=(((((((7.122452d-7*x-1.766345d-6)*x+2.928433d-5)*x
     &-2.335379d-4)*x+1.664156d-3)*x-1.041576d-2)*x+5.555682d-2)
     &*x-2.500001d-1)*x+9.999999d-1
      expi3=-1.d75
      if(x.ne.0.d0) expi3=x*aux-dlog(dabs(x))-5.772157d-1
      goto 3
2     if(x.gt.-9.d0) then
      aux=1.d0-((((5.176245d-2*x+3.061037d0)*x+3.243665d1)*x
     &+2.244234d2)*x+2.486697d2)/((((x+3.995161d0)*x+3.893944d1)
     &*x+2.263618d1)*x+1.807837d2)
      else
      y=9.d0/x
      aux=1.d0-y*(((y+7.659824d-1)*y-7.271015d-1)*y-1.080693d0)
     &/((((y*2.518750d0+1.122927d1)*y+5.921405d0)*y
     &-8.666702d0)*y-9.724216d0)
      endif
      expi3=aux*dexp(-x)/x
3     continue
      if(x.eq.0.d0) then
      expi3=0.5d0
      else
      expi3=0.5d0*((1.d0-x)*dexp(-x)+x**2*expi3)
      endif
      return       
      end function 




       

读取的数据文件长这样

读取的数据文件长这样
acc3177c1e58befbc1068611355ac52.png

155

帖子

2

主题

1

精华

大师

Vim

F 币
945 元
贡献
461 点

规矩勋章

发表于 2021-2-4 19:32:32 | 显示全部楼层
每次都输出一下ipl_r(i)的值,看看变化。

790

帖子

2

主题

0

精华

大宗师

F 币
3765 元
贡献
2255 点
发表于 2021-2-4 21:27:35 | 显示全部楼层
修改项目属性,使得 nan 报错。
1.png

250

帖子

2

主题

0

精华

宗师

F 币
1730 元
贡献
872 点

规矩勋章

发表于 2021-2-5 08:40:59 | 显示全部楼层
我编译执行了,没有出现NAN

20

帖子

8

主题

0

精华

熟手

F 币
126 元
贡献
73 点
 楼主| 发表于 2021-2-5 10:39:21 | 显示全部楼层
本帖最后由 xiao. 于 2021-2-5 10:45 编辑
li913 发表于 2021-2-4 21:27
修改项目属性,使得 nan 报错。

需要怎么解决呢?第一次遇到这种工问题,不太懂您的意思,能说的详细一点嘛,谢谢啦

20

帖子

8

主题

0

精华

熟手

F 币
126 元
贡献
73 点
 楼主| 发表于 2021-2-5 12:02:19 | 显示全部楼层
necrohan 发表于 2021-2-5 08:40
我编译执行了,没有出现NAN

为啥呢,你是用的我的数据吗,我有两组数据,一组NAN,一组就正常,

250

帖子

2

主题

0

精华

宗师

F 币
1730 元
贡献
872 点

规矩勋章

发表于 2021-2-6 10:29:05 | 显示全部楼层
xiao. 发表于 2021-2-5 12:02
为啥呢,你是用的我的数据吗,我有两组数据,一组NAN,一组就正常,

我就是按你给的数据的前4个数

20

帖子

8

主题

0

精华

熟手

F 币
126 元
贡献
73 点
 楼主| 发表于 2021-2-6 15:25:15 | 显示全部楼层
necrohan 发表于 2021-2-6 10:29
我就是按你给的数据的前4个数

前四个数没问题,算到第五万多行就开始nan了

20

帖子

8

主题

0

精华

熟手

F 币
126 元
贡献
73 点
 楼主| 发表于 2021-2-6 15:26:46 | 显示全部楼层
necrohan 发表于 2021-2-6 10:29
我就是按你给的数据的前4个数

十分感谢!方便加一下微信吗,我把原数据文件发你试一试?,纠结了好长时间了

250

帖子

2

主题

0

精华

宗师

F 币
1730 元
贡献
872 点

规矩勋章

发表于 2021-2-7 08:34:02 | 显示全部楼层
xiao. 发表于 2021-2-6 15:25
前四个数没问题,算到第五万多行就开始nan了

那可能是算法有问题,不懂的人没法帮忙,如果你不懂调试,最简单的方法是在计算语句前用write把几个用到的变量都输出来检查,比如
write(*,*)i,k,ipl_r(i),PL_IBG(k),temp_sum_l(1),expi3(temp_sum_l(1))
看看哪个变量先出现的NaN
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

捐赠本站|Archiver|关于我们 About Us|小黑屋|Fcode ( 京ICP备18005632-2号 )

GMT+8, 2024-3-29 15:25

Powered by Tencent X3.4

© 2013-2024 Tencent

快速回复 返回顶部 返回列表