Fortran Coder

标题: 数组越界 [打印本页]

---

作者: IDCY    时间: 2023-9-8 15:58
标题: 数组越界
自己编了一个求激波管问题精确解的代码 但是运行不出来 求各位大佬帮我看看哪里有问题
给下修改意见

D(ZQZ{DW1`7B]Z9)2E%PFOV.png (14.8 KB, 下载次数: 131)

下载附件

2023-9-8 15:57 上传

main.f90

2023-9-8 15:55 上传

点击文件名下载附件

2.44 KB, 下载次数: 2

---

作者: IDCY    时间: 2023-9-8 15:58

[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

---

作者: Transpose    时间: 2023-9-8 18:23

[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         ! 空间网格数

---

作者: IDCY    时间: 2023-9-9 15:01

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 Coder (http://bbs.fcode.cn/)

Powered by Discuz! X3.2