burgersof2d.f90

用ENO格式求解2-DBurgers方程精确捕捉激波

module define 
  implicit none
  integer,parameter::NX=80,NY=80,L=4
  real,parameter::eps=0.001
  real,parameter::hx=2.0/NX
  real,parameter::pi=3.1415926
  real::lampta,tao

  contains
    
	subroutine initial(u,x)!初值函数赋值
	  implicit none
	  real::u(-3:NX+3),x(-3:NX+3)
      integer::i,j

	  do i=-3,NX+3
	    u(i)=2+sin(pi*x(i))
	  end do
	end subroutine initial

	subroutine eno_numerical_divergence(divergence,ux,x,hx)!ENO格式
	  implicit none
	  real::divergence,ux(-3:3),hx,fluxr,fluxl,flux(-3:3),x(-3:3)
	  flux(:)=ux(:)**2/2.0  !流通量
	  call enoflux(fluxr,fluxl,flux,ux,x,hx)

	  divergence=fluxr

	  contains 

	  subroutine enoflux(fluxr,fluxl,flux,u,x,h)
	    use IMSL
		implicit none
		real::fluxr,fluxl,u(-3:3),flux(-3:3),flux1(-3:3),flux2(-3:3),h,d(-3:3),x(-3:3),&
		              A(3,3),b(3,1),c(3,1)
		integer::i,j,kmax,kmin

		flux1(:)=(flux(:)+1.0*u(:))/2.0
        flux2(:)=(flux(:)-1.0*u(:))/2.0

		kmax=0
		kmin=0
		d=flux1
		do i=2,3
		  do j=2,i-3,-1
		    d(j)=(d(j)-d(j-1))/(i-1)*h
		  end do
		  if(abs(d(kmax+1))>abs(d(kmax))) then
		  kmax=kmax
		  kmin=kmin-1
		  else
		  kmax=kmax+1
		  kmin=kmin
		  end if
		end do
        b(:,1)=flux1(kmin:kmax)
		do i=1,3
		  do j=1,3
		    a(i,j)=(kmin+i-1+0.5)**j/j-(kmin+i-1-0.5)**j/j		    
		  end do		  
		end do
		call lin_sol_gen(A,b,c)
		fluxr=c(1,1)+c(2,1)*0.5+c(3,1)*0.5**2

		kmax=1
		kmin=1
		d=flux2
		do i=2,3
		  do j=3,i-2,-1
		    d(j)=(d(j)-d(j-1))/(i-1)*h
		  end do
		  if(abs(d(kmax+1))>abs(d(kmax))) then
		  kmax=kmax
		  kmin=kmin-1
		  else
		  kmax=kmax+1
		  kmin=kmin
		  end if
		end do
		b(:,1)=flux2(kmin:kmax)
		do i=1,3
		  do j=1,3
		    a(i,j)=(kmin+i-2+0.5)**j/j-(kmin+i-2-0.5)**j/j
		  end do
		end do

		call lin_sol_gen(A,b,c)
		fluxr=c(1,1)+c(2,1)*(-0.5)+c(3,1)*(-0.5)**2+fluxr
        
      end subroutine enoflux

    end subroutine eno_numerical_divergence

end module define

program main
  use define
  real::u(-3:NX+3),divergence(0:NX),x(-3:NX+3)
  real::T
  integer::i,j,n,Step
  
  T=0.833
  lampta=0.1
  tao=lampta*hx
  Step=int(T/tao)
  do i=-3,NX+3
    x(i)=-1+i*hx
  end do
          
  call initial(u,x)

  do n=1,Step
    print*,n
	  
	  do i=0,NX
		  call eno_numerical_divergence(divergence(i),u(i-3:i+3),x(i-3:i+3),hx)		  
	  end do	 
	  	  
	  do i=1,NX
		  u(i)=u(i)-tao*(divergence(i)-divergence(i-1))/hx
	  end do
      u(0)=u(NX)
	  u(-1)=u(NX-1)
	  u(-2)=u(NX-2)
	  u(-3)=u(NX-3)
	  u(NX+1)=u(1)
	  u(NX+2)=u(2)
	  u(NX+3)=u(3)
  end do
  open(2006,file="u.txt")
  do i=0,NX
	  write(2006,*)u(i)
  end do
  close(2006)
  open(2007,file="x.txt")
  do i=0,NX
	  write(2007,*)x(i)
  end do
  close(2007)

end program