Fortran Coder»论坛 › 研讨区 › Fortran算法研究 › 一维黎曼问题求助

一维黎曼问题求助

查看数: 12345 | 评论数: 11 | 收藏 0

关灯 | 提示：支持键盘翻页<-左右->

帖子模式

ligenFortran

发布时间: 2020-6-12 22:36

正文摘要:

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

ligenFortran 发表于 2020-6-21 20:51:16

已解决，谢谢大家啦。真的是很粗心，漏掉了一个括号。通过输出中间结果发现问题的。

ligenFortran 发表于 2020-6-16 10:20:09

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

一维黎曼问题求助

正文摘要:

回复