一维黎曼问题求助

li913

你需要写清楚，改动了哪里？

ligenFortran

这是网上的代码


! MacCormack1D.for
!------------------------------------------------------------------------------------------------------------
!利用差分格式求解一维激波管问题（语言版本）(网上的代码)        
!--------------------------------------------------------------------------------------------------------------*/
      program MacCormack1D
      implicit double precision (a-h,o-z)
      parameter (M=1000)
      common /G_def/ GAMA,PI,J,JJ,dL,TT,Sf
       
      dimension U(0:M+1,0:2),Uf(0:M+1,0:2)
      dimension Ef(0:M+1,0:2)
!气体常数
      GAMA=1.4
      PI=3.1415926
!网格数
      J=M
!计算区域  
      dL=2.0
!总时间
      TT=0.45
!时间步长因子
      Sf=0.8
       
      call Init(U,dx)
      
      T=0

1     dt=CFL(U,dx)
      T=T+dt
      write(*,*)'T=',T,'dt=',dt
      call MacCormack_1D_Solver(U,Uf,Ef,dx,dt)
        call bound(U,dx)
      if(T.lt.TT)goto 1
                  
      call Output(U,dx)
      end
!-------------------------------------------------------
!计算时间步长
!入口: U， 当前物理量，dx, 网格宽度；
!返回: 时间步长。
!-------------------------------------------------------
      double precision function CFL(U,dx)
      implicit double precision (a-h,o-z)
      common /G_def/ GAMA,PI,J,JJ,dL,TT,Sf
      dimension U(0:J+1,0:2)

      dmaxvel=1e-10
      do 10 i=1,J
       uu=U(i,1)/U(i,0)
       p=(GAMA-1)*(U(i,2)-0.5*U(i,0)*uu*uu)  
       vel=dsqrt(GAMA*p/U(i,0))+dabs(uu)
       if(vel.gt.dmaxvel)dmaxvel=vel
10    continue
      CFL=Sf*dx/dmaxvel
      end
   
!-------------------------------------------------------
!初始化
!入口: 无；
!出口: U, 已经给定的初始值，dx，网格宽度。
!-------------------------------------------------------
      subroutine Init(U,dx)
      implicit double precision (a-h,o-z)
      common /G_def/ GAMA,PI,J,JJ,dL,TT,Sf
      dimension U(0:J+1,0:2)

!初始条件
      rou1=1.0
      u1=0
      v1=0
      p1=1.0
      rou2=0.125
      u2=0
      v2=0
      p2=0.1

      dx=dL/J

      do 20 i=0,J/2
       U(i,0)=rou1
       U(i,1)=rou1*u1
       U(i,2)=p1/(GAMA-1)+0.5*rou1*u1*u1
20    continue
      do 21 i=J/2+1,J+1
       U(i,0)=rou2
       U(i,1)=rou2*u2
       U(i,2)=p2/(GAMA-1)+0.5*rou2*u2*u2
21    continue
      end

!-------------------------------------------------------
!边界条件
!入口: dx，网格宽度；
!出口: U， 已经给定边界。
!-------------------------------------------------------
      subroutine bound(U,dx)
      implicit double precision (a-h,o-z)
      common /G_def/ GAMA,PI,J,JJ,dL,TT,Sf
      dimension U(0:J+1,0:2)

!左边界
      do 30 k=0,2
         U(0,k)=U(1,k)
30    continue
!右边界
      do 31 k=0,2
       U(J+1,k)=U(J,k)
31    continue
      end

!-------------------------------------------------------
!根据U计算E
!入口: U，当前U矢量；
!出口: E，计算得到的E矢量，
!     U、E定义见Euler方程组。
!-------------------------------------------------------
      subroutine U2E(U,E,is,in)
      implicit double precision (a-h,o-z)
      common /G_def/ GAMA,PI,J,JJ,dL,TT,Sf
      dimension U(0:J+1,0:2),E(0:J+1,0:2)
      
      do 40 i=is,in
       uu=U(i,1)/U(i,0)
       p=(GAMA-1)*(U(i,2)-0.5*U(i,1)*U(i,1)/U(i,0))
       E(i,0)=U(i,1)
       E(i,1)=U(i,0)*uu*uu+p
       E(i,2)=(U(i,2)+p)*uu
40    continue
      end

!-------------------------------------------------------
!一维差分格式求解器
!入口: U， 上一时刻U矢量，
!     Uf、Ef，临时变量，
!     dx，网格宽度，dt,，时间步长；
!出口: U， 计算得到得当前时刻U矢量。
!-------------------------------------------------------
      subroutine MacCormack_1D_Solver(U,Uf,Ef,dx,dt)
      implicit double precision (a-h,o-z)
      common /G_def/ GAMA,PI,J,JJ,dL,TT,Sf
      dimension U(0:J+1,0:2),Uf(0:J+1,0:2)
      dimension Ef(0:J+1,0:2)
      
      r=dt/dx
      dnu=0.25

      do 60 i=1,J
       do 60 k=0,2
!开关函数
        q=dabs(dabs(U(i+1,0)-U(i,0))-dabs(U(i,0)-U(i-1,0)))/(dabs(U(i+1,0)-U(i,0))+dabs(U(i,0)-U(i-1,0))+1e-10)
!人工黏性项      
        Ef(i,k)=U(i,k)+0.5*dnu*q*(U(i+1,k)-2*U(i,k)+U(i-1,k))
60    continue
      do 61 k=0,2
       do 61 i=1,J
        U(i,k)=Ef(i,k)
61    continue
      
      call U2E(U,Ef,0,J+1)
      
      do 63 i=0,J
       do 63 k=0,2
!U(n+1/2)(i+1/2)
        Uf(i,k)=U(i,k)-r*(Ef(i+1,k)-Ef(i,k))
63    continue                         

!E(n+1/2)(i+1/2)     
      call U2E(Uf,Ef,0,J)
     
      do 64 i=1,J
       do 64 k=0,2
!U(n+1)(i)
        U(i,k)=0.5*(U(i,k)+Uf(i,k))-0.5*r*(Ef(i,k)-Ef(i-1,k))
64    continue
      end

!-------------------------------------------------------
!输出结果, 用数据格式画图
!入口: U， 当前时刻U矢量，
!     dx，网格宽度；
!出口: 无。
!-------------------------------------------------------
      subroutine Output(U,dx)
      implicit double precision (a-h,o-z)
      common /G_def/ GAMA,PI,J,JJ,dL,TT,Sf
      dimension U(0:J+1,0:2)

      open(1,file='result.txt',status='unknown')
             
      do 80 i=0,J+1
       rou=U(i,0)
       uu=U(i,1)/rou
       p=(GAMA-1)*(U(i,2)-0.5*U(i,0)*uu*uu)
       write(1,81)i*dx,rou,uu,p,U(i,2)
80    continue
      close(1)
      
81    format(D20.10,D20.10,D20.10,D20.10,D20.10)
      end