作者YmemY (*＊米)
看板Fortran
标题[请益] 回圈跳不出来??
时间Tue Dec 18 02:55:26 2012
请教大家，我在做一个一维扩散的模拟，
在我一步一步debug之後，知道问题在回圈的部分，好像会无法结束??
可以成功编译成执行档，但一旦开始执行，就无法结束了，也无法输出结果..
以下是程式码，麻烦版有帮我看看哪边出错了??
program diffusion !1-D diffusion
implicit none
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!宣告区
real L,D,J,C,T,F,smallt
!L=length,D=diffusivity,J=diffusion flux,C=initial concentration
!F=Fourier's number, smallt=time interval
real, parameter :: defaultL = 10
real, parameter :: defaultD = 1.03
real, parameter :: defaultJ = 1E-5
real, parameter :: defaultC = 1E-4
real array1(-1:10001),array2(-1:10001)
!the length is divided into 10000 pieces, 10001 nodes, 2 virtual nodes
integer*8 p, pmax, x
!p=pace, x=location
character mode
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!输入设定参数
write(*,*) "press ""1"" to select the default values, or""2"" to input
them"
read(*,*) mode
if(mode=="1") then
L = defaultL
D = defaultD
J = defaultJ
C = defaultC
else
write(*,*) "please input the length:"
read(*,*)L
write(*,*) "please input the diffusivity:"
read(*,*)D
write(*,*) "please input the diffusion flux:"
read(*,*)J
write(*,*) "please input the initial concentration:"
read(*,*)C
end if
write(*,*) "please input the time:"
read(*,*)T
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!计算用到的参数
smallt = 0.25*L*L/D/10**8
pmax = T/smallt
F = D*smallt/(L*L)*10**8
array1 = C
array2 = C
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!这里就是出错的开始!!!
do p=1,pmax
do x=0,10000
array1(-1) = J*L/10000/D + array1(0)
array1(10001) = array1(10000) + array1(0) - array1(-1)
array2(x)=F*(array1(x+1)+array1(x-1))+(1-2*F)*array1(x)
end do
array1 = array2
end do
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!将阵列输出到文字档
open(unit=10,file="diffusion1.txt")
do x=-1,10001
write(10,*) array1(x)
end do
STOP
END
--
※ 发信站: 批踢踢实业坊(ptt.cc)
◆ From: 140.112.211.150
1F：→ Cypresslin:应该不是无法结束，是回圈太大，pmax还要乘上10000 12/18 09:52
2F：→ kkarthur:smallt = 0.25*L*L/D/10**8 检查你的smallt的值 12/18 12:25
3F：→ kkarthur:如果没有意外应该是0,导致pmax极大,回圈跑不完 12/18 12:26
4F：→ kkarthur:改成smallt = 0.25*L*L/D/10.0**8.0 试试看 12/18 12:26
5F：→ kkarthur:不仅是这个地方,例如array1(-1) = J*L/10000/D 12/18 12:30
6F：→ kkarthur:这都会'让你的结果与预期的不同 12/18 12:30
谢谢~
目前我把回圈数设定成100的话都可以执行成功，所以应该是回圈太大的问题，
0.25*L*L/D/10.0**8.0有改过但还是跑不完~
这问题满棘手的..因为有些参数就是得这麽小阿......QQ
※ 编辑: YmemY 来自: 140.112.4.209 (12/18 13:46)
7F：→ mk650:t=10,pamx=41200000,x=412000000000,粉大的回圈,如果参数没 12/18 13:56
8F：→ mk650:错,也只能丢着给它跑啦 12/18 13:56
9F：→ Cypresslin:你可以再把回圈里跟x无关的拿出来，可能可以快一咪咪 12/18 13:57
10F：→ mk650:呵三楼的精简化本人常干,不过它这个程式很简单,只有两行多余 12/18 14:00
11F：→ mk650:差别不大啦.重点在t要多大,t=0.0001够吗? 12/18 14:00
t最多可能要到100...囧
不过我想或许我的D数字错了，他可能是10^-10这麽小，所以smallt可以不用这麽小..
虽然试了一次还是行不通 但方向应该没错 再搞搞看...谢谢
12F：→ mk650:t=1,十分钟跑完,I7 2600套装很便宜电脑,更好的电脑得会更快 12/18 14:17
13F：→ mk650:这样很快了,我的程式都要跑一周以上,有三周还没跑完的 12/18 14:18
挖赛!!!!!
原来真的要这麽久XD
※ 编辑: YmemY 来自: 140.112.211.150 (12/18 22:53)
14F：推 george1224:建议把你所计算的数值的order提高一点，电脑在进行计算 12/19 22:32
15F：→ george1224:的时候，其精度会有差别... 12/19 22:33
感谢各位的提示，经由将参数调整地平易近人一点後，模拟出正确的结果了！
※ 编辑: YmemY 来自: 140.112.211.150 (12/28 02:08)