一维黎曼问题求助

ligenFortran · 发表于 2020-6-12 22:36:22

求助：我的计算结果，压力算的不对，我检查代码好多次，都不知道哪错了。这是我根据网上的代码略微修改，基本没改，求大神帮助。

!    MAIN (Explicit method)（我的代码）!    MacCormack_1D_Rimenn problem
!===============================================================

   PROGRAM MacCormack_1D_Rimenn
   IMPLICIT NONE
   INTEGER, PARAMETER :: NMAX = 1001
   REAL(KIND=8), DIMENSION(NMAX) :: X
   REAL(KIND=8), DIMENSION(NMAX,0:2) :: U, UB, UBB, E, EB, EV
   REAL(KIND=8) :: rou, uu, p, Te, a, M, h
   REAL(KIND=8) :: DT, t, TT=0.4                                                       !声明时间离散变量
   INTEGER :: n
   REAL(KIND=8) :: rou1, u1, p1, rou2, u2, p2
   REAL(KIND=8) :: DX
   INTEGER :: NX, IX                                                                      !声明空间离散变量
   REAL(KIND=8) :: sigma, PI = 4.0d0*DATAN(1.0d0), GAMA = 1.4, R = 278.0
   REAL(KIND=8) :: maxvel = 1e-10, vel, eta=0.25, theta
   INTEGER :: I, K
   WRITE (*, *) "Enter NX"
   !READ  (*, *) NX
! DT = adjustable
   NX = 1000                ! 网格数
   IX = NX+1
   DX = 2.0d0/NX             ! 计算域（0，2）
   DO I = 1, IX
x(i) = dble((I-1)*DX)
   end do                   ! 网格
!-------------------------------------------------------
! Shock tube initialtion condition
!-------------------------------------------------------
   rou1=1.0
   u1=0.0
   p1=1.0
   rou2=0.125
   u2=0.0
   p2=0.1

   DO I = 1, (IX-1)/2
U(i,0)=rou1
U(i,1)=rou1*u1
U(i,2)=p1/(GAMA-1)+0.5*rou1*u1*u1
   END DO
   DO I = (IX-1)/2+1, IX
U(i,0)=rou2
U(i,1)=rou2*u2
U(i,2)=p2/(GAMA-1)+0.5*rou2*u2*u2
   END DO

   t = 0
! March
! time step(capture signal---stability)
10    DO I = 2, IX-1
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. maxvel) then
maxvel = vel
end if
   END DO
   DT = 0.8*DX/maxvel
   sigma = DT/DX
!---------------------------------------------------------
! Computation
   t = t + DT
   write(*,*)'T=',T,'dt=',dt
   DO I = 2,IX-1
DO K = 0,2
         !开关函数
theta=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)
         !人工粘性修正
EV(i,k) = U(i,k)+0.5*eta*theta*(U(i+1,k)-2*U(i,k)+U(i-1,k))
END DO
   END DO
   DO K = 0,2
DO I = 2,IX-1
U(i,k) = EV(i,k)
END DO
   END DO

! Predictor
! UB_{i} = U^n_i - sigma * (E^n_(i+1) - E^n_i)
   do i=1,IX
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
   end do

   DO I = 1, IX-1
DO K = 0, 2
UB(I,K) = U(I,K) - sigma*( E(I+1,K) - E(I,K) )
END DO
   END DO
! Corrector
   do i=1,IX-1
uu=UB(i,1)/UB(i,0)
      p=(GAMA-1)*UB(i,2)-0.5*UB(i,1)*UB(i,1)/UB(i,0)
      EB(i,0)=UB(i,1)
      EB(i,1)=UB(i,0)*uu*uu+p
      EB(i,2)=(UB(i,2)+p)*uu
   end do

   DO I = 2, IX-1
DO k = 0,2
U(i,k)=0.5*(U(i,k)+UB(i,k))-0.5*sigma*(EB(i,k)-EB(i-1,k))
END DO
   END DO

! Boundary condition
   DO K = 0, 2
U(1,k) = U(2,k)
U(IX,k) = U(IX-1,k)
   END DO
!-------------------------------------------------------
   if (  t .lt. 0.4  ) goto 10
!-------------------------------------------------------
! Plot
   open(1,file='MacCormack_1D_Rimenn.txt',status='unknown')
   DO i=1,IX
   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)
   END DO
   close(1)

81 format(D20.10,D20.10,D20.10,D20.10,D20.10)


   END PROGRAM MacCormack_1D_Rimenn

ligenFortran · 发表于 2020-6-12 22:40:07

这是网上的代码

! 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

li913 · 发表于 2020-6-13 10:41:18

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

ligenFortran · 发表于 2020-6-15 08:06:42

改动了所有的变量声明，并且取消了所有的子程序

ligenFortran · 发表于 2020-6-15 08:08:54

li913 发表于 2020-6-13 10:41
你需要写清楚，改动了哪里？

您好，我改动了所有变量声明，并取消了子程序，可以帮我看看吗。很的感谢

necrohan · 发表于 2020-6-15 12:38:08

我的建议是先从原始程序开始，确认原始程序没问题，把你改的地方一点一点加进去，比如先改初始化部分，编译执行看有没有问题，然后改其它的。
原始程序可能有些不严谨的地方，你可以修改得更严谨一些，但是不建议改动下标范围、变量名称，使用别人的程序最好适应别人的思路，小程序好改，大程序就会很难办。
最好先改不容易错的地方，多做备份。

li913 · 发表于 2020-6-15 20:37:47

亲，我测试结果，两个代码的结果是一样的，都有台阶。

ligenFortran · 发表于 2020-6-15 21:16:11

带子程序的我的代码，和我的一整个主程序都算的不对

!    MAIN (Explicit method)
!    MacCormack_1D_Rimenn problem
!===============================================================

   PROGRAM MacCormack_1D_Rimenn
   IMPLICIT NONE
   INTEGER, PARAMETER :: NMAX = 1001
   REAL(KIND=8), DIMENSION(NMAX) :: X
   REAL(KIND=8), DIMENSION(NMAX,0:2) :: U, UB, UBB, E, EB
   REAL(KIND=8) :: uu, p
   REAL(KIND=8) :: DT, t, TT=0.4                                                       !声明时间离散变量
   INTEGER :: n
   REAL(KIND=8) :: DX
   INTEGER :: NX, IX                                                                      !声明空间离散变量
   REAL(KIND=8) :: sigma, PI = 4.0d0*DATAN(1.0d0), GAMA = 1.4, R = 278.0
   REAL(KIND=8) :: maxvel = 1e-10, vel
   INTEGER :: I, K
   WRITE (*, *) "Enter NX"
   READ  (*, *) NX
! DT = adjustable
! NX = 1000
   IX = NX+1
   DX = 2.0d0/NX             ! 计算域（0，2）
   DO I = 1, IX
x(i) = dble((I-1)*DX)
   end do                   ! 网格

! Initial Condition
   call Init(IX, U)
   t = 0

! March
! time step(capture signal---stability)
10    DO I = 2, IX-1
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. maxvel) then
maxvel = vel
end if
   END DO
   DT = 0.8*DX/maxvel
   sigma = DT/DX
! Computation
   t = t + DT
   call MacCormack_1D_Solver(U,sigma,IX)
   if (  t .lt. 0.4  ) goto 10
! Plot
   call Results(IX, U, DX)

   WRITE(*, *) "Numerical Solution is in MACCORMACK-1D-Rimenn.txt"
   WRITE(*, *) "Calculations are successfully completed. "
   WRITE(*, *) "Hit any key to close DOS window!"

   STOP
   END PROGRAM MacCormack_1D_Rimenn

!-------------------------------------------------------
! Shock tube initialtion condition
!-------------------------------------------------------
   subroutine Init(IIX, U)
   IMPLICIT NONE
   INTEGER, PARAMETER :: NMAX = 1001
   INTEGER, INTENT(IN) :: IIX
   REAL(KIND=8), INTENT(OUT), DIMENSION(NMAX,0:2) :: U
   REAL(KIND=8), DIMENSION(NMAX) :: X
   REAL(KIND=8) :: GAMA=1.4, rou1, u1, p1, rou2, u2, p2
   REAL(KIND=8) :: dx
   INTEGER :: i

   rou1=1.0
   u1=0.0
   p1=1.0
   rou2=0.125
   u2=0.0
   p2=0.1

   DO I = 1, IIX
if ( i .le. (IIX-1)/2 ) then
U(i,0)=rou1
U(i,1)=rou1*u1
U(i,2)=p1/(GAMA-1)+0.5*rou1*u1*u1
else
U(i,0)=rou2
U(i,1)=rou2*u2
U(i,2)=0.5*rou2*u2*u2 + p2/(GAMA-1)
end if
   END DO
   END subroutine Init

!-------------------------------------------------------
! 根据U计算E
! Input：U, 当前U矢量
! Export: E，计算得到的E矢量( flux term )
!-------------------------------------------------------
   Subroutine U2E(U,E,ist,ie)
   IMPLICIT NONE
   INTEGER, PARAMETER :: NMAX = 1001
   REAL(KIND=8), DIMENSION(NMAX) :: X
   INTEGER, INTENT(IN) :: ist,ie
   REAL(KIND=8), INTENT(IN), DIMENSION(NMAX,0:2) :: U
   REAL(KIND=8), INTENT(OUT), DIMENSION(NMAX,0:2) :: E
   REAL(KIND=8) :: uu, p
   REAL(KIND=8) :: gama=1.4
   INTEGER :: i

   do i=ist,ie
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
   end do
   end Subroutine U2E

!-------------------------------------------------------
! 1D_MacCormack差分格式求解器
! Input: U，上一时刻U矢量
! Export: U，计算得到得当前时刻U矢量( conserved variable )
!-------------------------------------------------------
   Subroutine MacCormack_1D_Solver(U,sigma,IX)
   IMPLICIT NONE
   INTEGER, PARAMETER :: NMAX = 1001
   REAL(KIND=8), INTENT(IN) :: sigma
   INTEGER, INTENT(IN) :: IX
   REAL(KIND=8), INTENT(INOUT), DIMENSION(NMAX,0:2) :: U
   REAL(KIND=8), DIMENSION(NMAX) :: X
   REAL(KIND=8), DIMENSION(NMAX,0:2) :: UB, UBB, E, EB,EV
   INTEGER :: I, K
   REAL(KIND=8) :: eta=0.25, theta

   DO I = 2,IX-1
DO K = 0,2
         !开关函数
theta=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)
         !人工粘性修正
EV(i,k) = U(i,k)+0.5*eta*theta*(U(i+1,k)-2*U(i,k)+U(i-1,k))
END DO
   END DO
   DO K = 0,2
DO I = 2,IX-1
U(i,k) = EV(i,k)
END DO
   END DO

! Predictor
! UB_{i} = U^n_i - sigma * (E^n_(i+1) - E^n_i)
   call U2E(U,E,1,IX)
   DO I = 1, IX-1
DO K = 0, 2
UB(I,K) = U(I,K) - sigma*( E(I+1,K) - E(I,K) )
END DO
   END DO
! Corrector
   call U2E(UB,EB,1,IX-1)
   DO I = 2, IX-1
DO k = 0,2
!UBB(I,K) = UB(I,K) - sigma * (EB(I,K) - EB(I-1,K))
U(i,k)=0.5*(U(i,k)+UB(i,k))-0.5*sigma*(EB(i,k)-EB(i-1,k))
END DO
   END DO
! Updating
   !DO k = 0,2
!DO I = 2, IX-1
! U(I,k) =  0.5d0*( U(I,k) + UBB(I,k) )
!END DO
   !END DO
! Boundary condition
   DO K = 0, 2
U(1,k) = U(2,k)
U(IX,k) = U(IX-1,k)
   END DO
   END Subroutine MacCormack_1D_Solver

!-------------------------------------------------------
! Plot Data
! Input: U，上一时刻U矢量
! Export: ，计算得到得当前时刻原始变量( conserved variable )
!-------------------------------------------------------
   SUBROUTINE Results(IX, U, DX)
   IMPLICIT NONE
   INTEGER, PARAMETER :: NMAX = 1001
   INTEGER, INTENT(IN) :: IX
   REAL(KIND=8), INTENT(IN) :: DX
   REAL(KIND=8), INTENT(IN), DIMENSION(NMAX,0:2) :: U
   REAL(KIND=8), DIMENSION(NMAX) :: X
   REAL(KIND=8) :: rou, uu, p, Te,a,M,h
   INTEGER :: I
   REAL(KIND=8) :: GAMA = 1.4, R = 278.0

   open(1,file='MacCormack_1D_Rimenn.txt',status='unknown')
   DO i=1,IX
rou=U(i,0)
      uu=U(i,1)/rou
      p=(GAMA-1)*(U(i,2)-0.5*U(i,0)*uu*uu)
Te=p/(rou*R)
a=sqrt( gama*R*Te )
M=uu/a
h=gama*R*Te/(gama-1)
      write(1,81) (I-1)*DX,rou,uu,p,Te,a,M,h
      81 FORMAT(8(F10.5, 1x))
   End DO
   close(1)
   END SUBROUTINE Results

ligenFortran · 发表于 2020-6-15 21:20:10

necrohan 发表于 2020-6-15 12:38
我的建议是先从原始程序开始，确认原始程序没问题，把你改的地方一点一点加进去，比如先改初始化部分，编译 ...

谢谢，我就先按您的方法找找问题在哪，加油

ligenFortran · 发表于 2020-6-16 10:19:07

li913 发表于 2020-6-15 20:37
亲，我测试结果，两个代码的结果是一样的，都有台阶。

不一样的，仔细看第三个变量的值，我的代码多了一个台阶

[通用算法] 一维黎曼问题求助

规矩勋章