Fortran Coder

查看: 177|回复: 11

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

[复制链接]

17

帖子

5

主题

0

精华

入门

F 币
98 元
贡献
61 点
发表于 2020-6-12 22:36:22 | 显示全部楼层 |阅读模式
2.png 1.png 微信截图_20200612221729.png
求助:我的计算结果,压力算的不对,我检查代码好多次,都不知道哪错了。这是我根据网上的代码略微修改,基本没改,求大神帮助。


!     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




2.png
1.png
回复

使用道具 举报

17

帖子

5

主题

0

精华

入门

F 币
98 元
贡献
61 点
 楼主| 发表于 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       

531

帖子

1

主题

0

精华

大宗师

F 币
2639 元
贡献
1684 点
发表于 2020-6-13 10:41:18 | 显示全部楼层
你需要写清楚,改动了哪里?

17

帖子

5

主题

0

精华

入门

F 币
98 元
贡献
61 点
 楼主| 发表于 2020-6-15 08:06:42 来自移动端 | 显示全部楼层
改动了所有的变量声明,并且取消了所有的子程序

17

帖子

5

主题

0

精华

入门

F 币
98 元
贡献
61 点
 楼主| 发表于 2020-6-15 08:08:54 | 显示全部楼层
li913 发表于 2020-6-13 10:41
你需要写清楚,改动了哪里?

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

128

帖子

2

主题

0

精华

大师

F 币
769 元
贡献
400 点

规矩勋章

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

531

帖子

1

主题

0

精华

大宗师

F 币
2639 元
贡献
1684 点
发表于 2020-6-15 20:37:47 | 显示全部楼层
亲,我测试结果,两个代码的结果是一样的,都有台阶。

17

帖子

5

主题

0

精华

入门

F 币
98 元
贡献
61 点
 楼主| 发表于 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

17

帖子

5

主题

0

精华

入门

F 币
98 元
贡献
61 点
 楼主| 发表于 2020-6-15 21:20:10 | 显示全部楼层
necrohan 发表于 2020-6-15 12:38
我的建议是先从原始程序开始,确认原始程序没问题,把你改的地方一点一点加进去,比如先改初始化部分,编译 ...

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

17

帖子

5

主题

0

精华

入门

F 币
98 元
贡献
61 点
 楼主| 发表于 2020-6-16 10:19:07 | 显示全部楼层
li913 发表于 2020-6-15 20:37
亲,我测试结果,两个代码的结果是一样的,都有台阶。

不一样的,仔细看第三个变量的值,我的代码多了一个台阶
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

QQ|捐赠本站|Archiver|关于我们 About Us|群聊|Fcode

GMT+8, 2020-7-8 19:44

Powered by Discuz! X3.2

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表