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

24小时热门版块排行榜    

CyRhmU.jpeg
小木虫论坛-学术科研互动平台 » 计算模拟区 » 程序语言 » Fortran » 【求助】FORTRAN中增加计算区域后出错，
61/1返回列表
查看: 910  |  回复: 5
只看楼主@他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

曾经电工

金虫 (小有名气)

[交流] 【求助】FORTRAN中增加计算区域后出错，已有2人参与

开始我计算点为50×100，程序能够正常编译和运行，结果也还是可以。但是当计算点增加到500×100时，可以运行，而出来的结果完全不正确。看了一下数据，发现是在计算浓度场中，增加计算点后其计算值出现负数。不知怎么回事？另外还有一个问题就是计算过程中本来值为1.000000000000000，但是显示结果却为1.00000000000001，再继续运行下去的话，会出现累计现象，从而影响结果的正确性。参数的精度我都定义为双精度的，这是怎么一个情况？
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
1楼2010-09-14 09:23:25
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

snoopyzhao

至尊木虫 (职业作家)
★ ★
小木虫(金币+0.5):给个红包，谢谢回帖交流
余泽成(金币+1):谢谢参与应助！ 2010-09-14 16:02:06
曾经电工(金币+5):谢谢您的帮助！ 2010-09-14 16:07:21
曾经电工: 回帖置顶 2011-12-06 09:50:19
第一个问题，无源码无真相。

第二个问题，基本上不用担心，所有的浮点计算都有一定的误差，但这点误差基本上可以忽略不计，我想
赞一下(2人) 回复此楼
2楼2010-09-14 14:20:36
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

曾经电工

金虫 (小有名气)
引用回帖:
Originally posted by snoopyzhao at 2010-09-14 14:20:36:
第一个问题，无源码无真相。

第二个问题，基本上不用担心，所有的浮点计算都有一定的误差，但这点误差基本上可以忽略不计，我想

谢谢您的回答，我的浓度计算程序如下，麻烦您帮我看一下是否是程序中的问题呢？我认为这种浮点计算的误差是在积累，最后都影响到结果的正确性了。
module const
implicit none
integer,parameter :: N=100                                
integer,parameter :: L=50                     
integer,parameter :: tmax=900000               
double precision,parameter :: Dis=2.0d-3              
double precision,parameter :: dx=4.0d-6            
double precision,parameter :: Ds=1.68d-9            
double precision,parameter :: Cs0=10.8d0           
double precision,parameter :: Cxm=20.0d0            
double precision,parameter :: Yxs=0.61d0            
double precision,parameter :: Ks=5.204d0              
double precision,parameter :: dts=1.0d0              
end module const


!=====================================
!              主程序
!=====================================
     program main
         use const
     implicit none
         integer :: it,x,z,ci(1:N,1:L)
         double precision :: Sd(1:N,1:L),Cd(1:N,1:L)
         
         call init_Field(Sd,Cd,ci)
         
         do it=1,tmax
         
        call substrate(Sd,Cd,ci)
        
     end do
   
    stop
        end

      subroutine init_Field(Sd,Cd,ci)
      use const
      implicit none
      double precision :: a,Sd(1:N,1:L),Cd(1:N,1:L)  
      integer :: ci(1:N,1:L)
          integer :: x,z
      
          do x=1,N
             do z=1,L
                    ci(x,z)=0
                Sd(x,z)=1.0d0
                Cd(x,z)=0.0d0
         end do
          end do

          call random_seed()
          do x=1,N
         z=1
                 call random_number(a)
         Cd(x,z)=a
         if (Cd(x,z)<=0.5d0) then
                ci(x,z)=1         
         else
                    Cd(x,z)=0.0d0
                 end if
      end do
                  
      return
      end
!=========================================================================================
! 浓度场迭代计算: 在每个坐标方向上分别采用追赶法进行求解三对角方程组，
! 其他方向则按照显示处理的交替方向隐式方法ADI(alternating direction implicit)
!=========================================================================================
          subroutine substrate(Sd,Cd,ci)
          use const
          implicit none
          integer :: x,z,ci(1:N,1:L),Zm,Zb,k(1:L)
          double precision :: Sd(1:N,1:L),Cd(1:N,1:L),Ud(1:N,1:L),Md(1:N,1:L),a1(1:N,1:L),a2(1:N,1:L),a3(1:N,1:L),aa1(1:N,1:L),aa2(1:N,1:L),aa3(1:N,1:L)  ! a1,a2,a3,b1为三对角阵方程系数(x方向)；ps为底物消耗量
      double precision :: ps2(1:N,1:L),f(1:N,1:L),fe(1:N,1:L),y(1:N,1:L),m(1:N,1:L),me(1:N,1:L),r(1:N,1:L),tt,kk,I,b1(1:N,1:L),b2(1:N,1:L),ps1(1:N,1:L)
          double precision :: um,ms
          
          ms=0.562137d0
      um=0.25986d0
          Zm=0
          k(1:L)=0
          do z=1,L
             do x=1,N
                    k(z)=k(z)+ci(x,z)
                 end do
                 if (k(z)==0) then
                        Zm=z-1
                        goto 300
                 else
                    continue
                 end if
           end do
          
300   Zb=Zm+10
      
          do x=1,N
             do z=Zb+1,L
                    Sd(x,z)=1.0d0
                 end do
          end do
          
          do x=1,N
             do z=1,L
                    Ud(x,z)=Sd(x,z)
                 end do
          end do
          
          
          !!! X-direction
          tt=um
          kk=ms
          a1(1:N,1:Zb)=-1.0d0
          a3(1:N,1:Zb)=-1.0d0
          a2(1:N,1:Zb)=2.0d0+(2.0d0*(Dis**2)*(dx**2))/(dts*Ds)
          do z=2,Zb
             do x=2,N-1
                    b1(x,z)=Sd(x,z+1)+((2.0d0*(Dis**2)*(dx**2))/(dts*Ds)-2.0d0)*Sd(x,z)+Sd(x,z-1)-ps1(x,z)
                        ps1(x,z)=((Dis**2)*(dx**2)/Ds)*((tt*Cxm*Cd(x,z)*Sd(x,z))/(Yxs*Cs0*(Ks/Cs0+Sd(x,z)))+Cxm*Cd(x,z)*kk/Cs0)
                 end do
          end do
          
                                                 
    do z=2,Zb
           f(2,z)=a3(2,z)/a2(2,z)
           fe(2,z)=a2(2,z)
            do x=3,N-2
                   fe(x,z)=a2(x,z)-(a1(x,z)*f(x-1,z))
                   f(x,z)=a3(x,z)/fe(x,z)
            end do
                fe(N-1,z)=a2(N-1,z)-(a1(N-1,z)*f(N-2,z))
!----------------------------------------------------------
            y(2,z)=(b1(2,z)-a1(2,z)*Sd(1,z))/a2(2,z)                        
                do x=3,N-2
                   y(x,z)=(b1(x,z)-a1(x,z)*y(x-1,z))/fe(x,z)
                end do
        y(N-1,z)=(b1(N-1,z)-a3(N-1,z)*Sd(N,z)-a1(N-1,z)*y(N-2,z))/fe(N-1,z)
!----------------------------------------------------------
                Ud(N-1,z)=y(N-1,z)
                do x=N-2,2,-1
                   Ud(x,z)=y(x,z)-f(x,z)*Ud(x+1,z)
        end do
!-----------------------------------------------------------
                Ud(N,z)=(4.0d0*Ud(2,z)+4.0d0*Ud(N-1,z)-Ud(N-2,z)-Ud(3,z))/6.0d0
                Ud(1,z)=Ud(N,z)
        end do
   
        do x=1,N
           Ud(x,1)=(4.0d0*Ud(x,2)-Ud(x,3))/3.0d0
       Ud(1,1)=Ud(N,1)
    end do
       
        do x=1,N
           do z=1,L
              Md(x,z)=Ud(x,z)
           end do
        end do

        

      !!! Z-direction
      aa1(1:N,1:Zb)=-1.0d0
          aa3(1:N,1:Zb)=-1.0d0
          aa2(1:N,1:Zb)=2.0d0+(2.0d0*(Dis**2)*(dx**2))/(dts*Ds)
          do x=2,N-1
             do z=2,Zb
                    b2(x,z)=Md(x+1,z)+((2.0d0*(Dis**2)*(dx**2))/(dts*Ds)-2.0d0)*Md(x,z)+Md(x-1,z)-ps2(x,z)
                        ps2(x,z)=((Dis**2)*(dx**2)/Ds)*((tt*Cxm*Cd(x,z)*Md(x,z))/(Yxs*Cs0*(Ks/Cs0+Md(x,z)))+Cxm*Cd(x,z)*kk/Cs0)
         end do
          end do
          
             do x=2,N-1      
                m(x,2)=aa3(x,2)/aa2(x,2)
                me(x,2)=aa2(x,2)
            do z=3,Zb-1   
                       me(x,z)=aa2(x,z)-(aa1(x,z)*m(x,z-1))
                       m(x,z)=aa3(x,z)/me(x,z)
                end do
                        me(x,Zb)=aa2(x,Zb)-(aa1(x,Zb)*m(x,Zb-1))
!------------------------------------------------------------------------------------
                    r(x,2)=(b2(x,2)-aa1(x,2)*Md(x,1))/aa2(x,2)
                        do z=3,Zb-1
                           r(x,z)=(b2(x,z)-aa1(x,z)*r(x,z-1))/me(x,z)
                        end do
            r(x,Zb)=(b2(x,Zb)-aa3(x,Zb)*Md(x,Zb+1)-aa1(x,Zb)*r(x,Zb-1))/me(x,Zb)
!-----------------------------------------------------------------------------------
                        Sd(x,Zb)=r(x,Zb)
                        do z=Zb-1,2,-1
                       sd(x,z)=r(x,z)-m(x,z)*Sd(x,z+1)
            end do
!-----------------------------------------------------------------------------------
                        Sd(x,1)=(4.0d0*Sd(x,2)-Sd(x,3))/3.0d0
             end do
                 do z=1,Zb
                    Sd(N,z)=(4.0d0*Sd(N-1,z)+4.0d0*Sd(2,z)-Sd(3,z)-Sd(N-2,z))/6.0d0
                        Sd(1,z)=Sd(N,z)
                 end do
         
         return
             end

[ Last edited by 曾经电工 on 2010-9-14 at 16:13 ]
赞一下 回复此楼
3楼2010-09-14 16:04:57
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

snoopyzhao

至尊木虫 (职业作家)
你这个程序的输出在哪里呢？
赞一下 回复此楼
4楼2010-09-14 19:54:31
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

曾经电工

金虫 (小有名气)
输出的话，我可以在主程序里面加入的：
program main
         use const
     implicit none
         integer :: it,x,z,ci(1:N,1:L)
         double precision :: Sd(1:N,1:L),Cd(1:N,1:L)
         
         call init_Field(Sd,Cd,ci)
         
         do it=1,tmax
         
        call substrate(Sd,Cd,ci)
        do x=1,N
                   do z=1,L
                      write(*,*) x,z,Sd(x,z)
                   end do
                end do
     end do
   
    stop
        end
赞一下 回复此楼
5楼2010-09-15 08:35:07
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zyj8119

木虫 (著名写手)
★
曾经电工(金币+1):谢谢参与
曾经电工(金币+5):谢谢您的帮助！ 2010-11-10 08:21:13
引用回帖:
Originally posted by 曾经电工 at 2010-09-14 16:04:57:

谢谢您的回答，我的浓度计算程序如下，麻烦您帮我看一下是否是程序中的问题呢？我认为这种浮点计算的误差是在积累，最后都影响到结果的正确性了。
module const
implicit none
integer,parameter :: N=100   ...

我觉得你最好在程序里面打开一个文件，然后把所有的数据都存放在文件中，我用CVF6.6运行你的程序什么东西都没有。
赞一下(1人) 回复此楼
好好学习，天天向上。
6楼2010-11-09 19:53:29
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 曾经电工 的主题更新
61/1返回列表
普通表情 龙 兔 虎 猫 高级回复(可上传附件)
信息提示
关闭
请填处理意见
关闭