Fortran Coder

查看: 8014|回复: 2
打印 上一主题 下一主题

[求助] 显示数组溢出,实在找不到问题了,请大家帮忙

[复制链接]

20

帖子

8

主题

0

精华

熟手

F 币
115 元
贡献
71 点
跳转到指定楼层
楼主
发表于 2015-3-1 08:12:33 | 只看该作者 回帖奖励 |倒序浏览 |阅读模式
[Fortran] 纯文本查看 复制代码
program DVM
!
    implicit none   !-->all the variables MUST be declared
!
    integer,parameter :: n = 10
    !Declaration of variable
    real(8) :: U,a,dt
    real(8) :: xs,ys   
    integer :: p,q
    integer :: i
    !p for which vortex it is, q for how long it exist
    real(8),dimension(n-1):: gama
    !the circulation of each vortex
    real(8),dimension(n,n) :: xx,yy,xxi,yyi
    !the location of the vortices and their corresponding images, the first index refers to which vortex it is, and the second one refers to how long it have existed
    
    print*,'input the Velocity of the free stream, the Radius of the cylinder'
    read*,U,a
    print*,'input the Unit Time Step'
    read*,dt
    xs = 0
    ys = a
        
    !uu = U(1-((xx(p,q)**2 - yy(p,q)**2)*a**2)/((xx(p,q)**2 + yy(p,q)**2)**2))
    !vv = -2*U*xx(p,q)*yy(p,q)*a**2/(xx(p,q)**2+yy(p,q)**2)**2
    !formula (14) (15) in the lecture note
    
    do p = 1,n-1
        !set q = 1 to get the new pop-up vortex---------------------
        i = 1
        q = 1
        !define the value of x of every pop-up vortex---------------
        xx(p,q) = 0
        !define the value of y of every pop-up vortex--------------
        if (MOD(p,2)/=0) then
            yy(p,q) = 1.1*a
        else 
            yy(p,q) = -1.1*a
        end if
        !calculate gamma of the pop-up vortex----------------------
        q = p
        call getgama() 
        q = p
        !calculate the new location of the vortices at present-----
        
            do i = 1,p
                call getlocation()
                q = q - 1
            end do

    end do  
    print*,gama

    
    
contains 
    subroutine image(xxx,yyy,xxxi,yyyi)
        implicit none
        real(8) :: xxx,yyy,xxxi,yyyi
        xxxi = a**2*(xxx/(xxx**2 + yyy**2))
        yyyi = a**2*(yyy/(xxx**2 + yyy**2))
   
        return
    end subroutine image




    subroutine getgama()
        implicit none
        real(8), dimension(p) :: simp
        !to simplify the calculation, use simp to represent the iterative part of formula (18)
        real(8), parameter :: Pi = 3.1415927      
        if(p==1) then  
            simp(1) = -2*Pi*U*(1-((xs**2 - ys**2)*(a**2))/((xs**2 + ys**2)**2))
            call image(xx(1,1),yy(1,1),xxi(1,1),yyi(1,1))
            gama(1) = simp(1)/(-(ys-yy(1,1))/((xs-xx(1,1))**2+(ys-yy(1,1))**2)+(ys-yyi(1,1))/((xs-xxi(1,1))**2+(ys-yyi(1,1))**2))
        else
            if (MOD(p,2)/=0) then
                do while (i<p)
                    call image(xx(i,q),yy(i,q),xxi(i,q),yyi(i,q))
                    simp(i+1) = simp(i)+gama(i)*(ys-yy(i,q))/((xs-xx(i,q))**2+(ys-yy(i,q))**2)-gama(i)*(ys-yyi(i,q))/((xs-xxi(i,q))**2+(ys-yyi(i,q))**2)
                    q = q-1
                    i = i+1
                end do
                gama(p) = simp(p)/(-(ys-yy(p,1))/((xs-xx(p,1))**2+(ys-yy(p,1))**2)+(ys-yyi(p,1))/((xs-xxi(p,1))**2+(ys-yyi(p,1))**2))
            else 
                do while (i<p)
                    call image(xx(i,q),yy(i,q),xxi(i,q),yyi(i,q))
                    simp(i+1) = simp(i)+gama(i)*(-ys-yy(i,q))/((xs-xx(i,q))**2+(-ys-yy(i,q))**2)-gama(i)*(-ys-yyi(i,q))/((xs-xxi(i,q))**2+(-ys-yyi(i,q))**2)
                    q = q-1
                    i = i+1
                end do  
                gama(p) = simp(p)/(-(-ys-yy(p,1))/((xs-xx(p,1))**2+(-ys-yy(p,1))**2)+(-ys-yyi(p,1))/((xs-xxi(p,1))**2+(-ys-yyi(p,1))**2))
            end if
        end if
           
        return 
    end subroutine getgama

    subroutine getlocation()
        implicit none
        !p tells how many vortices do we have now
        !we should move (p-1) vortices
        !q tells which stage is each vortex at
        real, parameter :: Pi = 3.1415927
        real(8) :: Ustream,Vstream
        real(8) :: simpxm = 0,simpxp = 0,simpym = 0,simpyp = 0
        integer :: j,m
        m = q
        
            do j = 1,p
                if(j/=i) then
                    q = p + 1 - j
                    simpxm = simpxm + dt*( Gama(j) / (2*Pi) ) * (yy(i,q) - yy(j,q)) / ((xx(i,q)-xx(j,q))**2 + (yy(i,q)-yy(j,q))**2)
				end if
			end do
            
            do j = 1,p
                q = p + 1 - j
				simpxp = simpxp + dt*( Gama(j) / (2*Pi) ) * (yy(i,q) - yyi(j,q)) / ((xx(i,q)-xxi(j,q))**2 +(yy(i,q)-yyi(j,q))**2)
			end do
            
            do j = 1,p
                if(j/=i) then
                    q = p + 1 - j
                    simpym = simpym + dt*( Gama(j) / (2*Pi) ) * (xx(j,q) - xx(i,q)) / ((xx(i,q)-xx(j,q))**2 + (yy(i,q)-yy(j,q))**2)
				end if
			end do
            
            do j = 1,p
                    q = p + 1 - j
				    simpyp = simpyp + dt*( Gama(j) / (2*Pi) ) * (xxi(j,q) - xx(i,q)) / ((xx(i,q)-xxi(j,q))**2 +(yy(i,q)-yyi(j,q))**2)
			end do
            
            Ustream = U * (1-(xx(i,m)**2-yy(i,m)**2) * a**2 / (xx(i,m)**2+yy(i,m)**2)**2)      
        
		    Vstream = -2 * U * xx(i,m) * yy(i,m) * a**2 / (xx(i,m)**2+yy(i,m)**2)**2 
            
		    xx(i,m+1) = xx(i,m) + Ustream * dt - simpxm + simpxp

		    yy(i,m+1) = yy(i,m) + Vstream * dt - simpym + simpyp
            

        
        
        return
    end subroutine getlocation
end program DVM



这是留的大作业,研究好久了,感觉不难,但做起来很繁琐,麻烦大家了,另外请问有没有调试程序的好方法,一个一个找太累了!!!
    

分享到:  微信微信
收藏收藏 点赞点赞 点踩点踩

954

帖子

0

主题

0

精华

大师

F 币
184 元
贡献
75 点

规矩勋章元老勋章新人勋章水王勋章热心勋章

QQ
沙发
发表于 2015-3-1 11:51:40 | 只看该作者
因为不知道你的算法,所以只是猜想。

主程序里 43 到 49 行,q 在循环中变小。每次循环 q = q - 1

但是 getlocation 函数里,又改变了 q,比如 q = p + 1 - j,这样很危险,不是吗?

46

帖子

8

主题

0

精华

熟手

F 币
211 元
贡献
131 点
板凳
发表于 2015-3-25 22:00:02 | 只看该作者
我觉得问题也许是在subroutine getgama中。一般我的经验是按步骤测试,首先去掉sub,运行主程序,然后一步一步加入call。加入getgamma后,print*, gama 没有正确显示数值(出现NaN)。

纯属个人意见,仅供参考。
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

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

GMT+8, 2024-4-27 17:59

Powered by Tencent X3.4

© 2013-2024 Tencent

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