作者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)