solvey.for

2D有限差分法程序,未加边界条件,适合一般情况下的波动方程模拟

!     copyright:wangaifeng 
!     DATE     :2007-07-25
!     Company  :GUCAS
      subroutine Solvey(u,gt,nx,ny,ighost,gama)
	implicit real(a-h,o-z)
	real u(-ighost:nx+ighost,-ighost:ny+ighost,0:3)
      real gl(-ighost:nx+ighost,-ighost:ny+ighost,0:3)
	real gr(-ighost:nx+ighost,-ighost:ny+ighost,0:3)
	real dgl(-ighost:nx+ighost,-ighost:ny+ighost,0:3)
	real dgr(-ighost:nx+ighost,-ighost:ny+ighost,0:3)
	real ggl(-ighost:nx+ighost,-ighost:ny+ighost,0:3)
	real ggr(-ighost:nx+ighost,-ighost:ny+ighost,0:3)
	real gt(-ighost:nx+ighost,-ighost:ny+ighost,0:3)
	real minmod


	eps=1e-6

      do i=-ighost,nx+ighost
	do j=-ighost,ny+ighost
	
      p=(gama-1.0)*(u(i,j,3)-0.5*u(i,j,1)**2/u(i,j,0)
     & -0.5*u(i,j,2)**2/u(i,j,0))
      
	uu=u(i,j,1)/u(i,j,0)
	c=sqrt(gama*p/u(i,j,0))
	vv=u(i,j,2)/u(i,j,0)
	H=c*c/(gama-1)+0.5*(uu**2+vv**2)

	a1=vv
	a2=vv-c
	a3=vv+c

	a1l=0.5*(a1+sqrt(a1*a1+eps*eps))
	a2l=0.5*(a2+sqrt(a2*a2+eps*eps))
	a3l=0.5*(a3+sqrt(a3*a3+eps*eps))

	a1r=0.5*(a1-sqrt(a1*a1+eps*eps))
	a2r=0.5*(a2-sqrt(a2*a2+eps*eps))
	a3r=0.5*(a3-sqrt(a3*a3+eps*eps))

	gl(i,j,0)=u(i,j,0)/2/gama*(2*(gama-1)*a1l+a2l+a3l)
	gl(i,j,2)=u(i,j,0)/2/gama*(2*vv*(gama-1)*a1l
     &	+(vv-c)*a2l+(vv+c)*a3l)
	gl(i,j,1)=u(i,j,0)*uu/2/gama*(2*(gama-1)*a1l
     &	+a2l+a3l)
	gl(i,j,3)=u(i,j,0)/2/gama*(2*((gama-1)*H-c*c)*a1l
     &	+(H-c*vv)*a2l+(H+vv*c)*a3l)


	gr(i,j,0)=u(i,j,0)/2/gama*(2*(gama-1)*a1r+a2r+a3r)
	gr(i,j,2)=u(i,j,0)/2/gama*(2*vv*(gama-1)*a1r
     &	+(vv-c)*a2r+(vv+c)*a3r)
	gr(i,j,1)=u(i,j,0)*uu/2/gama*(2*(gama-1)*a1r
     &	+a2r+a3r)
	gr(i,j,3)=u(i,j,0)/2/gama*(2*((gama-1)*H-c*c)*a1r
     &	+(H-c*vv)*a2r+(H+vv*c)*a3r)



	enddo

      

	do j=-ighost+1,ny+ighost
	 do k=0,3
	  dgl(i,j,k)=gl(i,j,k)-gl(i,j-1,k)
	  dgr(i,j,k)=gr(i,j,k)-gr(i,j-1,k)
	 enddo
	enddo



	do j=-ighost+1,ny+ighost-1
	 do k=0,3
	y=minmod(dgl(i,j,k),dgl(i,j+1,k))
	ggl(i,j,k)=gl(i,j,k)+0.5*y
	 enddo
	enddo


	do j=-ighost+1,ny+ighost-2
	 do k=0,3
	y=minmod(dgr(i,j+1,k),dgr(i,j+2,k))
	ggr(i,j,k)=gr(i,j+1,k)-0.5*y
	 enddo
	enddo

	do j=-ighost+1,ny+ighost-2
	 do k=0,3
	gt(i,j,k)=ggl(i,j,k)+ggr(i,j,k)
	enddo
	enddo

      
      enddo
	end