tvd.f90

使用TVD格式求解流体力学中的激波管问题

module define
  implicit none
  integer,parameter::Nmax=500
  real,parameter::a=0.2,lampta=0.1
  real,parameter::eps=0.00000001
  real,parameter::h=1.0/Nmax,tao=lampta*h/a

  contains

  real function minmod(x,y)
    implicit none
	real::x,y
	minmod=(signal(x)/2.0+signal(y)/2.0)*min(abs(x),abs(y))
	contains
	integer function signal(x)
	  implicit none
	  real::x
	  if(x>=0) then
        signal=1
	  else
	    signal=-1
	  end if
	end function signal
  end function minmod

  real function g(u)
    implicit none
	real::u(-1:1)
    g=minmod(tuta_g(u(-1:0)),tuta_g(u(0:1)))
    contains
	real function tuta_g(u)
	  implicit none
	  real::u(2)
	  tuta_g=(limiter(lampta)-lampta**2)*(u(2)-u(1))/2.0
	end function tuta_g
  end function g

  real function gama(u)
    implicit none
	real::u(-1:2)
	if(u(1)-u(0)==0) then
	  gama=0.0
	else 
	  gama=(g(u(0:2))-g(u(-1:1)))/(u(1)-u(0))
	end if
  end function gama

  subroutine computation_flux(flux,u)
    implicit none
	real::flux(0:Nmax),u(-2:Nmax+2)
    integer::j

	do j=0,Nmax
	  flux(j)=(a*u(j)+h*g(u(j-1:j+1))/tao+a*u(j+1)+h*g(u(j:j+2))/tao&
	          -h*limiter(lampta+gama(u(j-1:j+2)))*(u(j+1)-u(j))/tao)/2.0
    end do
  end subroutine computation_flux
  
  subroutine initial(u,x)
    implicit none 
	real::u(-2:Nmax+2),x(-2:Nmax+2)
	integer::j

	do j=-2,Nmax+2
	  if(x(j)<=0.2) then
	    u(j)=0.0
	  else if(x(j)>0.2.and.x(j)<=0.5) then
	    u(j)=sin((x(j)-0.2)*10*3.1415926)
	  else if(x(j)>0.5.and.x(j)<=0.7) then
	    u(j)=7.5*(x(j)-0.5)
	  else 
	    u(j)=-1.0
	  end if
    end do
  end subroutine
  real function limiter(x)
    implicit none
	real::x
	if(abs(x)>=eps) then
	  limiter=abs(x)
	else
	  limiter=(x**2+eps**2)/(2*eps)
    end if
  end function

end module

program main
  use define
  implicit none
  real::flux(0:Nmax),u(-2:Nmax+2),u_1(-2:Nmax+2),x(-2:Nmax+2),T
  integer::j,n,Step
  T=1.0
  Step=int(T/tao)
  do j=-2,Nmax+2
    x(j)=0+j*h
  end do

  call initial(u,x)
  do n=1,Step
    print*,n
    call computation_flux(flux,u)
	
	do j=1,Nmax
	  u_1(j)=u(j)-tao*(flux(j)-flux(j-1))/h
	end do

	u_1(0)=0.0
	u_1(-1)=u_1(0)
	u_1(-2)=u_1(0)
	u_1(Nmax+1)=u_1(Nmax)
	u_1(Nmax+2)=u_1(Nmax+1)

    call computation_flux(flux,u_1)
    do j=1,Nmax
	  u_1(j)=(u_1(j)-tao*(flux(j)-flux(j-1))/h)
	end do

	u_1(0)=0.0
	u_1(-1)=u_1(0)
	u_1(-2)=u_1(0)
	u_1(Nmax+1)=u_1(Nmax)
	u_1(Nmax+2)=u_1(Nmax+1)
	u=(u+u_1)/2.0
  end do
  open(2006,file="u.txt")
  do j=0,Nmax
    write(2006,*)u(j)
  end do
  open(2007,file="x.txt")
  do j=0,Nmax
    write(2007,*)x(j)
  end do
end program