もう一歩のところで進まない

 一次元浅水流方程式(波動方程式)の数値計算がうまくいかない。実行可能ファイルは作成できるが、計算結果がおかしい。途中の変数を出力してみるのだが、プログラムの改善点を見つけるには至らない。
 以下がプログラムだ。

module globals

!implicit none

integer,save:: n,itrmax,pintv,i,fopen

real(8),save:: dt,xl,fr,gr,h0,dx,dh

real(8),allocatable,save:: x(:),u(:),v(:),w(:),h(:),p(:),q(:),pn(:),qn(:)

!public n,dt,xl,dx,fr,gr,itrmax,pintv,h0,dh,x,u,v,w,h,p,q,pn,qn

end module globals




module subprogs

use globals

implicit none
contains

subroutine set_init(n,dt,xl,dx,fr,gr,itrmax,pintv,h0,dh,&

x,u,v,w,h,p,q,pn,qn)

!implicit none

real(8),allocatable:: x(:),u(:),v(:),w(:),h(:),p(:),q(:),pn(:),qn(:)

integer n,itrmax,pintv

real(8) dt,xl,fr,gr,h0,dx,dh

open(10,file='data.d')

write(*,*)'data読み込み'

read(10,*) n,itrmax,pintv

read(10,*) dt,xl,fr,gr,h0

close(10)

write(*,*)'data読み込み終了'

dx=xl/dble(n-1)

dh=0.1d0*h0

allocate (x(n),u(n),v(n),w(n),h(n),p(n),q(n),pn(n),qn(n))

x(:)=(/(dx*dble(i-1),i=1,n)/)

u(:)=sqrt(gr*h0)*fr

h(:)=h0+dh*exp(-(x(:)-0.5d0*xl)**2/1.0d2)

call uh2pqvw(u,h,n)

end subroutine set_init


subroutine uh2pqvw(u,h,n)

integer n,i

real(8) u(n),v(n),w(n),h(n),p(n),q(n)

!allocate(u(n),v(n),w(n),h(n),pn(n),qn(n))

!real(8) gr

do i=1,n

v(i) = u(i) + sqrt(gr * h(i))

w(i) = u(i) - sqrt(gr * h(i))

p(i) = sqrt(gr * h(i)) + 1.0 / 2.0 * u(i)

q(i) = sqrt(gr * h(i)) - 1.0 / 2.0 * u(i)

enddo

write(*,*) '初期条件'

end subroutine uh2pqvw




subroutine print_uh(x,h,u)

!implicit none

real(8),allocatable::x(:),h(:),u(:)

!allocate (x(n),h(n),u(n))

!real(8) gr

!integer fopen

!fopen = 0

if(fopen==1)then

open(20,file='uh.d')

else if (fopen==-1) then


close(20)

return

endif

do i=1,n

write(20,'(10e16.8)')x(i),h(i),u(i),u(i)/sqrt(gr*h(i))

enddo

write(20,*)' '

end subroutine print_uh


subroutine chk_cno(dx,dt,v,w)

!implicit none

real(8) cno1,cno2,cno,dt,dx

real(8) v(:),w(:)

cno1=maxval(abs(v(:))*dt/dx)

cno2=maxval(abs(w(:))*dt/dx)

cno=max(cno1,cno2)

if(cno>=1)then

write(*,*) 'stop, cno >= 1, cno = ',cno

stop

endif

end subroutine chk_cno


subroutine cm1d(pn,p,i,v,dt,dx)

!implicit none

integer i,n

real(8) ,allocatable::v(:),pn(:),p(:)

!real(8) v(n),p(n),q(n),w(n),pn(n),qn(n)

real(8) dt,dx,cno

cno=v(i)*dt/dx

write(*,*)'クーラン数',cno

if(cno>=0.0d0)then

pn(i)=p(i) - cno * (p(i ) - p(i-1))

else

pn(i)=p(i) - cno * (p(i+1) - p(i ))

endif

end subroutine cm1d


subroutine bc_thru(pn,qn,p,q,v,w,n,dt,dx)

!implicit none

integer i,n

real(8),allocatable::p(:),q(:),pn(:),qn(:),v(:),w(:)

real(8) cno1,cno2,dt,dx

do i=1,n,n-1

pn(i)=p(i)

qn(i)=q(i)

cno1=v(i)*dt/dx

cno2=w(i)*dt/dx

if(i==1)then

if (cno1 < 0.0d0) pn(i) = p(i) - cno1 * (p(i+1) - p(i))

if (cno2 < 0.0d0) qn(i) = q(i) - cno2 * (q(i+1) - q(i))

endif

enddo

end subroutine bc_thru


subroutine pq2uhvw(pn,qn,gr,u,h,v,w,n)

!implicit none

!real(8) ,allocatable:: u(:),p(:),q(:),v(:),w(:),h(:),pn(:),qn(:)

real(8) gr,cno

integer i,n

real(8) x(n),u(n),v(n),w(n),h(n),p(n),q(n),pn(n),qn(n)

do i=1,n

u(i) = p(i) - q(i)

cno = 1.0 / 2.0 * (p(i) + q(i))

h(i) = cno**2 / gr

v(i) = u(i) + cno

w(i) = u(i) - cno

enddo

end subroutine pq2uhvw

end module subprogs



program main

use globals

use subprogs

implicit none

!real(8),allocatable :: x(:),u(:),v(:),w(:),h(:)

!real(8),allocatable :: p(:),q(:),pn(:),qn(:)

!real(8) dt,xl,dx,fr,gr,h0,dh

integer itr

call set_init(n,dt,xl,dx,fr,gr,itrmax,pintv,h0,dh,&

x,u,v,w,h,p,q,pn,qn)

write(*,*)'set_init ok'

fopen=1

call print_uh(x,h,u)

write(*,*)'初期値'

do itr=1 ,itrmax

call chk_cno(dx,dt,v,w)

write(*,*)'check クーラン数'

do i=2,n-1

call cm1d(pn,p,i,v,dt,dx)

call cm1d(qn,q,i,w,dt,dx)

enddo

write(*,*)'cm1d ok'

call bc_thru(pn,qn,p,q,v,w,n,dt,dx)

write(*,*) 'bc_thru ok'

call pq2uhvw(pn,qn,gr,u,h,v,w,n)

write(*,*)'pq2uhvw ok'

p(:)=pn(:)

q(:)=qn(:)

if(mod(itr,pintv)==0)then

fopen=0

call print_uh(x,h,u)

write(*,*)'経過出力'

endif

enddo

fopen=-1

call print_uh(x,h,u)

write(*,*)'計算終了'

deallocate(x,u,v,w,h,p,q,pn,qn)

end program main






!data.d


100 200 20

1.0d-1 1.0d2 0.5d0 9.8d0 1.0d-1

ブログ気持玉

クリックして気持ちを伝えよう!

ログインしてクリックすれば、自分のブログへのリンクが付きます。

→ログインへ

なるほど(納得、参考になった、ヘー)
驚いた
面白い
ナイス
ガッツ(がんばれ!)
かわいい

気持玉数 : 0

この記事へのコメント

この記事へのトラックバック