作者lprsdmn (ziv)
看板Fortran
标题[问题] 解暂态方程式
时间Tue Apr 26 19:19:12 2011
我要使用差分法解一个波动方程式
假设是二维的
这样要解的是u
所以宣告u(100,100,100)
前两个分别是X方向和Y方向的网格数,後面那个是t
也就是u(x,y,t)
假设我给dt是0.01
总时间长T给定1
这样就有了100格,一切都很顺利
我想改的问题来了~
今天我想写一个副程式sub来解上面的波动方程式
但是时间方面我想改成以下这样
do it=1, 1000
call sub
end do
it每跑一步就是跑了一个dt
我想让他跑1000步,也就是不会被限制宣告那边就要给他设定t的网格数
所以就不用设定时间总长
请问我该如何做呢?
--
※ 发信站: 批踢踢实业坊(ptt.cc)
◆ From: 140.117.58.27
1F：→ D49361128:你可以预设一个超大的U EX:U(100,100,1000000) 04/26 21:43
2F：→ D49361128:然後设定dt=1.d0/t_max tmax随你设不要超过预先最大的 04/26 21:44
3F：→ D49361128:然後跑的时候 do t=1,t_max 这样就解决了 不过要记得 04/26 21:45
4F：→ D49361128:读取的时候最後一格为 t_max 不是1000000 不然会读到0 04/26 21:46
5F：→ D49361128:多的格子就不要理他了XD (小心不要超过记忆体) 04/26 21:46
6F：→ D49361128:不然另一个方法就是把决定U的地方包成另一个副程式 由 04/26 21:49
7F：→ D49361128:主程式宣告并给她t_max EX:U(100,100,t_max) 04/26 21:49
8F：→ D49361128:两者来说 第一种只要不要超过记忆体就没啥问题 04/26 21:51
9F：→ D49361128:第二种 在t_max 太大的时候可能U会开不成 04/26 21:54
D大你说到我的问题点了,第一种的话 太容易超出过记忆体
第二种时间run太长也不行,照我目前的需求,时间可能至少需要10万格
两种方法我都有试,但是都无法顺利解决
10F：推 terryys:真的需要全部t的网格都存下来吗? 计算上最多应该只要有 04/26 22:16
11F：→ terryys:几个delta t就够了吧 04/26 22:16
12F：→ D49361128:可能他要绘图之类的吧 有些暂态模拟好像会讨论 04/26 23:18
t大说到我的需求了,我并不想把所有的t网格都存下来
我想改成只看需要的那几个dt就够了
ex.设定跑1000步,每100步输出资料
13F：推 jsb:不是用动态宣告就行了喔. . . 还是我搞错问题了？ 04/27 09:44
※ 编辑: lprsdmn 来自: 140.117.58.27 (04/27 10:39)
14F：推 terryys:每一步输出资料不行吗 假如你需要t=n*dt和t=(n+1)*dt来 04/27 11:04
15F：→ terryys:算t=(n+2)*dt,那你算完以後就可以输出t=(n+2)*dt然後用 04/27 11:05
16F：→ terryys:t=(n+1)*dt,t=(n+2)*dt就可以算t=(n+3)*dt,这样只需要3个 04/27 11:06
17F：→ terryys:(100,100)的array 04/27 11:06
t大可以帮我看看我程式该如何修改吗
你真的了解
program wave
implicit none
real(8) u(0:101,0:101,0:10001)
real(8) lamda,vel,c,dx,dy,dt,x,y,z,tfinal
integer i,j,k,nxstep,nystep,ntime,L
real, parameter :: pi=3.14159
c=2.9d8
L=1
tfinal=2.1d-9
dx=0.01*L
dy=0.01*L
dt=2.1d-13
nxstep=(L/dx)
nystep=(L/dy)
ntime=int(tfinal/dt)+1
!b.c
do k=0,ntime
do j=0,nystep
u(0,j,k)=0.
u(nxstep,j,k)=0.
end do
end do
do k=0,ntime
do i=0,nxstep
u(i,0,k)=0
u(i,nystep,k)=0
end do
end do
!i.c
do i=1,nxstep-1
do j=1,nystep-1
x=dx*i
y=dy*j
u(i,j,0)=sin(pi*x)*sin(pi*y)
u(i,j,1)=sin(pi*x)*sin(pi*y)
end do
end do
!==============================
do i=1,nxstep-1
do j=1,nystep-1
do k=1,ntime
u(i,j,k+1)=2*u(i,j,k)-u(i,j,k-1)+c*c*((dt/dx)**2*(u(i+1,j,k)-2*u(i,j,k)&
+u(i-1,j,k))+(dt/dy)**2*(u(i,j+1,k)-2*u(i,j,k)+u(i,j-1,k)))+100*dt*dt
end do
end do
end do
!=========================================================
open (7, file='ya1.plt')
open (8, file='ya2.plt')
open (9, file='ya3.plt')
do i=0,nxstep
write(7,*) i , u(i,50,5)
write(8,*) i , u(i,50,20)
write(9,*) i , u(i,50,100)
end do
stop
end
※ 编辑: lprsdmn 来自: 140.117.58.27 (04/27 14:54)