Transpose 发表于 2023-9-8 18:23
[mw_shl_code=fortran,true]  function flux(rho, u, p, i) result(f)
    real, intent(in) :: rho(0:), u ...

十分感谢大佬

[Fortran] 纯文本查看 复制代码

  function flux(rho, u, p, i) result(f)
    real, intent(in) :: rho(0:), u(0:), p(0:)
    integer, intent(in) :: i
    real :: f
   
    f = rho(i) * u(i)
   
  end function flux

这种传递数组的方式如果没有写下界，默认是1，你这里需要改为默认下界是0

这个最好写成integer

[Fortran] 纯文本查看 复制代码

integer, parameter :: nx = 20         ! 空间网格数

[Fortran] 纯文本查看 复制代码

program shock_tube
  implicit none
  
  ! 常数定义
  real, parameter :: gamma = 1.4       ! 气体比热比
  real, parameter :: L = 10.0           ! 激波管长度
  real, parameter :: nx = 20         ! 空间网格数
  real, parameter :: dx = L / nx       ! 空间步长
  real, parameter :: t_final = 1.6     ! 求解的时间
  real, parameter :: dt = 0.08       ! 时间步长  
  
  ! 变量定义
  integer :: i, n, tsteps
  real :: x, t
  real, dimension(0:nx) :: rho, u, p
  open(100,file='resulr.dat')
  ! 初始条件
  do i = 0, nx
    x = i * dx
   
    if (x <= 0.5 * L) then
      rho(i) =0.445d0
          u(i) = 0.698d0
      p(i) = 3.528d0
    else
      rho(i) = 0.5d0
      u(i) = 0.0d0
      p(i) = 0.571d0
    end if
  end do
  
  ! 求解时间步进
  t = 0.0
  tsteps = t_final / dt
  
  do n = 1, tsteps
    ! Riemann问题求解
    do i = 0, nx-1
      x = (i + 0.5) * dx
      
      if (x <= 0.5 * L) then
       rho(i) =0.445d0
          u(i) = 0.698d0
      p(i) = 3.528d0
    else
      rho(i) = 0.5d0
      u(i) = 0.0d0
      p(i) = 0.571d0
      end if
    end do
   
    ! 更新各个网格的值
    do i = 1, nx-2
      rho(i) = rho(i) - dt / dx * (flux(rho, u, p, i) - flux(rho, u, p, i-1))
      u(i) = u(i) - dt / dx * (flux(rho*u, u*u + p, u*(u(i)+p(i)/(rho(i)*(gamma-1))), i) - &
                               flux(rho*u, u*u + p, u*(u(i-1)+p(i-1)/(rho(i-1)*(gamma-1))), i-1))
      p(i) = p(i) - dt / dx * (flux((p/(gamma-1)+0.5*rho*u*u), (u*(u(i)+p(i)/(rho(i)*(gamma-1)))), &
                               ((rho*u*u+p)*(u(i)+p(i)/(rho(i)*(gamma-1))))/(gamma-1), i) - &
                               flux((p/(gamma-1)+0.5*rho*u*u), (u*(u(i-1)+p(i-1)/(rho(i-1)*(gamma-1)))), &
                               ((rho*u*u+p)*(u(i-1)+p(i-1)/(rho(i-1)*(gamma-1))))/(gamma-1), i-1))
    end do
   
    ! 更新左边界值
    i = 0
    rho(i) = rho(i+1)
    u(i) = -u(i+1)
    p(i) = p(i+1)
   
    ! 更新右边界值
    i = nx-1
    rho(i) = rho(i-1)
    u(i) = -u(i-1)
    p(i) = p(i-1)
   
    ! 更新时间
    t = t + dt
  end do
  
write(100,*) 'variables=', '"x","rho","u","p"'

! 输出结果
  do i = 0, nx
    x = i * dx
        write(100,"(4(e15.7))") x, rho(i), u(i), p(i)
  end do
close(100)

contains
  
  ! 定义通量函数
  function flux(rho, u, p, i) result(f)
    real, intent(in) :: rho(:), u(:), p(:)
    integer, intent(in) :: i
    real :: f
   
    f = rho(i) * u(i)
   
  end function flux

end program shock_tube