Fortran Coder

查看: 1142|回复: 3

[数值问题] 数组越界

[复制链接]

12

帖子

4

主题

0

精华

入门

F 币
52 元
贡献
26 点
发表于 2023-9-8 15:58:04 | 显示全部楼层 |阅读模式
自己编了一个求激波管问题精确解的代码 但是运行不出来 求各位大佬帮我看看哪里有问题
给下修改意见
D(ZQZ{DW1`7B]Z9)2E%PFOV.png

main.f90

2.44 KB, 下载次数: 2

12

帖子

4

主题

0

精华

入门

F 币
52 元
贡献
26 点
 楼主| 发表于 2023-9-8 15:58:42 | 显示全部楼层
[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

160

帖子

2

主题

1

精华

大师

Vim

F 币
965 元
贡献
470 点

规矩勋章

发表于 2023-9-8 18:23:32 | 显示全部楼层
[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         ! 空间网格数

12

帖子

4

主题

0

精华

入门

F 币
52 元
贡献
26 点
 楼主| 发表于 2023-9-9 15:01:40 | 显示全部楼层
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 ...

十分感谢大佬
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

捐赠本站|Archiver|关于我们 About Us|小黑屋|Fcode ( 京ICP备18005632-2号 )

GMT+8, 2024-6-24 04:45

Powered by Tencent X3.4

© 2013-2024 Tencent

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