aliouying 发表于 2014-5-16 17:02 好的,谢谢热心关注,祝你答辩顺利 |
算法看得不是很明白,最近毕业的事情都快疯了~ 没太多时间帮忙看 |
我的神啊,怎么犯了这么个低级错误啊,我要头撞南墙了,无地自容啊,跪谢大神 |
[Fortran] 纯文本查看 复制代码 module BLSM implicit none integer, parameter:: N=1000 real(kind=8), parameter:: dx=1000.d0 real(kind=8), parameter:: dy=1000.d0 real(kind=8), parameter:: dis=200000.d0 integer(kind=8), parameter:: M=-int(dis/dx) real(kind=8), parameter:: tao=0.00036d0 real(kind=8), parameter:: v_F=5.556d0 real(kind=8), parameter:: h=10.d0 real(kind=8), parameter:: detp=00.d0/(1030.d0*9.8d0) real(kind=8), parameter:: pi=3.1415926535898d0 real(kind=8), parameter:: f=0.0000773336d0 real(kind=8), parameter:: rho=1030.d0 real(kind=8), parameter:: g=9.8d0 real(kind=8), parameter:: ro=100000.d0 real(kind=8), parameter:: alpha=1E-10 real(kind=8), parameter:: lamda=1E-10 real(kind=8), parameter:: EPS=1E-28 endmodule BLSM module YDHS use BLSM implicit none contains FunctionB_u(v,zeta,i,j,cishu) implicitnone real(kind=8)::B_u,P_Dx,P_Dx_ZH,P_Dx_DD,ii,jj,P_Dxx real(kind=8)::v(M:N+M,-N-1:N+1) real(kind=8)::zeta(M:N+M,-N:N),ii_2,jj_2 integer::i,j,cishu ii=i*dx jj=j*dy ii_2=ii*ii jj_2=jj*jj B_u=-dy*f/(2.d0*V_F)*(v(i,j+1)+v(i-1,j+1)+v(i-1,j)+v(i,j))& +2.d0*g*dy/(V_F*dx)*(zeta(i,j)-zeta(i-1,j))& +(4.d0*dy*g*alpha*detp/V_F)*ii& *dexp(-alpha*ii_2-lamda*jj_2)& +(2.d0*dy*tao*ro**4.d0/(rho*h*V_F))/((ro*ro+ii_2+jj_2)**2.d0)& *(ii/2.d0+(dsqrt(3.d0)/2.d0)*jj) return endFunctionB_u FunctionB_v(u,zeta,i,j,cishu) implicitnone real(kind=8)::B_v,P_Dy_ZH,ii,jj,ii_2,jj_2 real(kind=8)::u(M:N+M+1,-N:N) real(kind=8)::zeta(M:N+M,-N:N) integer::i,j,cishu,cishu_S_P_Dy ii=i*dx jj=j*dy ii_2=ii*ii jj_2=jj*jj B_v=f*dy/(2.d0*V_F)*(u(i,j)+u(i,j-1)+u(i+1,j-1)+u(i+1,j))& +(2.d0*g/V_F)*(zeta(i,j)-zeta(i,j-1))& +4.d0*g*lamda*detp*(dy/V_F)*jj*dexp(-alpha*ii_2-lamda*jj_2)& -(2.d0*tao*ro**4.d0*dy/(rho*h*V_F))/((ro*ro+ii_2+jj_2)**2.0)& *(dsqrt(3.d0)/2.d0*ii-jj/2.d0) return endFunctionB_v Function B_zeta(u,v,i,j) implicit none real(kind=8)::B_zeta real(kind=8)::u(M:N+M+1,-N:N) real(kind=8)::v(M:N+M,-N-1:N+1) integer::i,j B_zeta=(2.d0*h*dy/(V_F*dx))*(u(i+1,j)-u(i,j))& +(2.d0*h/V_F)*(v(i,j+1)-v(i,j)) return end Function B_zeta !______________ end module YDHS program main use BLSM use YDHS implicit none real(kind=8)::ii,jj,P_Dx real(kind=8)::u0(M:N+M+1,-N:N),v0(M:N+M,-N-1:N+1) ,zeta0(M:N+M,-N:N) real(kind=8)::u1(M:N+M+1,-N:N),v1(M:N+M,-N-1:N+1) ,zeta1(M:N+M,-N:N) integer::i,j,cishu,k real(kind=8)::dx2 integer::cishu_u,cishu_v,cishu_zeta character*80::filename open (30,file='P_Dy_ZH.dat') open (20,file='P_Dx_ZH.dat') !初值全部设定为0 cishu=0 u0=0.d0 v0=0.d0 zeta0=0.d0 ! 先赋初值再使用 u1=0.d0 v1=0.d0 zeta1=0.d0 !----------------------- u1(M,:) =0.d0 !u的岸边界值 u1(M+1,:) =0.d0 !u的第一次迭代边界值 v1(:,-N) =0.d0 !v的水边界值 v1(:,-N+1)=0.d0 !v第一次迭代的边界值 zeta1(:,:)=0.d0 100 write(*,*) "迭代次数",cishu cishu=cishu+1 cishu_u=cishu+10000 write(filename,FMT='(I6.6,A4)')cishu_u,'.dat' filename=trim(filename) open(100,file=filename,status='unknown') doi=M+1,N+M ii=i*dx if(i==300.or.i==600)then write(*,*) i,'I''m_u' end if u1(i,-N)=u1(i,-N+1) !输运条件导数为0 doj=-N+1,N-1 jj=j*dy u1(i,j+1)=u1(i,j-1)+B_u(v0,zeta0,i,j,cishu) write(100,*)ii,jj, u1(i,j+1) end do end do close (100) !---------------------- cishu_v=cishu+20000 write(filename,FMT='(I6.6,A4)')cishu_v,'.dat' filename=trim(filename) open(101,file=filename,status='unknown') do i=M,N+M ii=i*dx if(i==300.or.i==600)then write(*,*) i,'I''m_v' end if v1(:,-N-1)=v1(:,-N) !提输运条件 do j=-N+1,N jj=j*dy v1(i,j+1)=v1(i,j-1)+B_v(u0,zeta0,i,j,cishu) write(101,*) ii,jj, v1(i,j+1) end do end do close (101) !------------------------------------------------ cishu_zeta=cishu+30000 write(filename,FMT='(I6.6,A4)')cishu_zeta,'.dat' filename=trim(filename) open(1001,file=filename,status='unknown') do i=M,N+M ii=i*dx if(i==300.or.i==600)then write(*,*) i,'Im_zeta' end if doj=-N+1,N-1 jj=j*dy zeta1(i,j+1)=zeta1(i,j-1)+B_zeta(u1,v1,i,j) write(1001,*) ii,jj, zeta1(i,j+1) !030001.dat end do end do close (1001) !________判断前后两次差别[/align] [align=left] do i=M,N+M do j=-N,N dx2=dx2+(u1(i,j)-u0(i,j))**2.d0+(v1(i,j)-v0(i,j))**2.d0 +(zeta1(i,j)-zeta0(i,j))**2.d0 end do end do dx2=dsqrt(dx2) write(*,*) "前后差距",dx2 if(dx2>=eps)then u0=u1 v0=v1 zeta0=zeta1 goto 100 end if stop endprogrammain 这就是只算了最后一个。不是累加。累加应该是 dx2 = dx2 + ... |
vvt 发表于 2014-5-13 22:37 我的dx2就是累加啊?不是算最后一个,dx2是u1,v1,zeta1 三个数组对应的数字与前面一次计算的u0,v0,zeta0的差值啊。还是谢谢你的关注。 |
太长了,实在没精力看。 我觉得应该是楼主的代码写得不好,而不是算法问题。 你确定 dx2 就只算最后一个?而不是累加? |
感谢各位,帮忙看看的大神,这个破问题纠结了2个月了,郁闷死了 |
因为长度受到限制,所以就没有吧计算结果的图发上来,大家可以运行一下,迭代下去会逐步变大的 |
捐赠本站|Archiver|关于我们 About Us|小黑屋|Fcode ( 京ICP备18005632-2号 )
GMT+8, 2024-12-26 05:17