[問題] 解暫態方程式

[lprsdmn (ziv)]

我要使用差分法解一個波動方程式

假設是二維的

這樣要解的是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
→ D49361128:你可以預設一個超大的U   EX:U(100,100,1000000)
→ D49361128:然後設定dt=1.d0/t_max  tmax隨你設不要超過預先最大的
→ D49361128:然後跑的時候  do t=1,t_max  這樣就解決了  不過要記得
→ D49361128:讀取的時候最後一格為 t_max 不是1000000   不然會讀到0
→ D49361128:多的格子就不要理他了XD  (小心不要超過記憶體)
→ D49361128:不然另一個方法就是把決定U的地方包成另一個副程式  由
→ D49361128:主程式宣告並給她t_max  EX:U(100,100,t_max)
→ D49361128:兩者來說  第一種只要不要超過記憶體就沒啥問題
→ D49361128:第二種 在t_max 太大的時候可能U會開不成
D大你說到我的問題點了,第一種的話 太容易超出過記憶體
第二種時間run太長也不行,照我目前的需求,時間可能至少需要10萬格
兩種方法我都有試,但是都無法順利解決

推 terryys:真的需要全部t的網格都存下來嗎? 計算上最多應該只要有
→ terryys:幾個delta t就夠了吧
→ D49361128:可能他要繪圖之類的吧  有些暫態模擬好像會討論
t大說到我的需求了,我並不想把所有的t網格都存下來
我想改成只看需要的那幾個dt就夠了
ex.設定跑1000步,每100步輸出資料


推 jsb:不是用動態宣告就行了喔. . . 還是我搞錯問題了？
※ 編輯: lprsdmn         來自: 140.117.58.27        (04/27 10:39)
推 terryys:每一步輸出資料不行嗎 假如你需要t=n*dt和t=(n+1)*dt來
→ terryys:算t=(n+2)*dt,那你算完以後就可以輸出t=(n+2)*dt然後用
→ terryys:t=(n+1)*dt,t=(n+2)*dt就可以算t=(n+3)*dt,這樣只需要3個
→ terryys:(100,100)的array
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)