| 查看: 1029 | 回复: 5 | |||
曾经电工金虫 (小有名气)
|
[交流]
【求助】FORTRAN中增加计算区域后出错， 已有2人参与
|
| 开始我计算点为50×100，程序能够正常编译和运行，结果也还是可以。但是当计算点增加到500×100时，可以运行，而出来的结果完全不正确。看了一下数据，发现是在计算浓度场中，增加计算点后其计算值出现负数。不知怎么回事？另外还有一个问题就是计算过程中本来值为1.000000000000000，但是显示结果却为1.00000000000001，再继续运行下去的话，会出现累计现象，从而影响结果的正确性。参数的精度我都定义为双精度的，这是怎么一个情况？ |
回复此楼» 猜你喜欢
277求调剂 数一104分
已经有7人回复
-
288环境专硕,求调材料方向
已经有16人回复
-
277工科求调剂
已经有8人回复
-
266求调剂
已经有19人回复
-
环境专硕调剂
已经有3人回复
-
081200-11408-276学硕求调剂
已经有5人回复
-
材料调剂
已经有13人回复
-
工科 22408 267求推荐
已经有4人回复
-
305求调剂
已经有3人回复
-
270分求调剂
已经有5人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- 急求fortran运行错误原因，在线等
已经有7人回复
- 请教 fortran 运行错误的原因
已经有13人回复
- fortran转换CHGCAR出错
已经有4人回复
- FORTRAN新手 求助主程序循环问题
已经有10人回复
- 我编的Simpson积分法fortran程序给不出结果，大侠们看看哪里出了问题？
已经有4人回复
- 请教一个fortran小程序编译出错的问题，谢谢
已经有9人回复
- 【求助】Fortran语言赋值问题？
已经有3人回复
- Fortran的格式化输入输出问题
已经有14人回复
- 有关fortran的一次而问题，希望大家能帮帮忙，谢谢
已经有4人回复
- 【求助】用fortran解矩阵问题【已解决】
已经有5人回复
snoopyzhao
至尊木虫 (职业作家)
- 程序强帖: 16
- 应助: 157 (高中生)
- 贵宾: 0.02
- 金币: 18844.7
- 红花: 29
- 帖子: 3803
- 在线: 1422.4小时
- 虫号: 183750
- 注册: 2006-02-13
- 专业: 污染生态化学
2楼2010-09-14 14:20:36
曾经电工
金虫 (小有名气)
- 应助: 0 (幼儿园)
- 金币: 2090.4
- 散金: 68
- 帖子: 119
- 在线: 224.2小时
- 虫号: 582687
- 注册: 2008-07-25
- 性别: GG
- 专业: 燃烧学
|
谢谢您的回答，我的浓度计算程序如下，麻烦您帮我看一下是否是程序中的问题呢？我认为这种浮点计算的误差是在积累，最后都影响到结果的正确性了。 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
snoopyzhao
至尊木虫 (职业作家)
- 程序强帖: 16
- 应助: 157 (高中生)
- 贵宾: 0.02
- 金币: 18844.7
- 红花: 29
- 帖子: 3803
- 在线: 1422.4小时
- 虫号: 183750
- 注册: 2006-02-13
- 专业: 污染生态化学
4楼2010-09-14 19:54:31
曾经电工
金虫 (小有名气)
- 应助: 0 (幼儿园)
- 金币: 2090.4
- 散金: 68
- 帖子: 119
- 在线: 224.2小时
- 虫号: 582687
- 注册: 2008-07-25
- 性别: GG
- 专业: 燃烧学
5楼2010-09-15 08:35:07
6楼2010-11-09 19:53:29
