solvex.for

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

!     copyright:wangaifeng 
!     DATE     :2007-07-25
!     Company  :GUCAS
      subroutine Solvex(u,ft,nx,ny,ighost,gama)
	implicit real(a-h,o-z)
	real u(-ighost:nx+ighost,-ighost:ny+ighost,0:3)
      real fl(-ighost:nx+ighost,-ighost:ny+ighost,0:3)
	real fr(-ighost:nx+ighost,-ighost:ny+ighost,0:3)
	real dfl(-ighost:nx+ighost,-ighost:ny+ighost,0:3)
	real dfr(-ighost:nx+ighost,-ighost:ny+ighost,0:3)
	real ffl(-ighost:nx+ighost,-ighost:ny+ighost,0:3)
	real ffr(-ighost:nx+ighost,-ighost:ny+ighost,0:3)
	real ft(-ighost:nx+ighost,-ighost:ny+ighost,0:3)
	real minmod


	eps=1e-6


      do j=-ighost,ny+ighost
	
	do i=-ighost,nx+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=uu
	a2=uu-c
	a3=uu+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))

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

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


	enddo

	do i=-ighost+1,nx+ighost
	 do k=0,3
	dfl(i,j,k)=fl(i,j,k)-fl(i-1,j,k)
	dfr(i,j,k)=fr(i,j,k)-fr(i-1,j,k)
	 enddo
	enddo

	do i=-ighost+1,nx+ighost-1
	 do k=0,3
	y=minmod(dfl(i,j,k),dfl(i+1,j,k))
	ffl(i,j,k)=fl(i,j,k)+0.5*y
	 enddo
	enddo

	do i=-ighost+1,nx+ighost-2
	 do k=0,3
	y=minmod(dfr(i+1,j,k),dfr(i+2,j,k))
	ffr(i,j,k)=fr(i+1,j,k)-0.5*y
	 enddo
	enddo

	do i=-ighost+1,nx+ighost-2
	 do k=0,3
	ft(i,j,k)=ffl(i,j,k)+ffr(i,j,k)
	enddo
	enddo

      
      enddo
	end