Fortran嵌套循环变量无法随空间更新的技术求助

输入参数 #

输入文件proj1pt1data.txt内容：

13.00 -11.00 0.00002 0.0143 400000.00 11

模拟代码 #

program my_project
    implicit none

    ! Variable Declarations
    real :: Q, z0, y0, slope, mann, domain, S_b, S_f, Fr, dx, dt, g, U, C, p, ps, ds, d, dz_dx
    integer :: i, n_steps, n_time_steps, t_index, x_index
    real, dimension(:), allocatable :: X, z_vals, y_vals, U_vals, q_t, dz_vals

    ! Read Input File
    open(11, file='proj1pt1data.txt')
    read(11, *) Q, z0, slope, mann, domain, y0
    close(11)

    ! Constants
    g = 9.81
    dx = 100.0
    dt = 1.0  
    n_steps = int(domain / dx)  
    n_time_steps = 100  

    ! Morphodynamic Constants
    p = 0.4
    ps = 2650
    ds = 0.0001
    d = (ps - 1000) / 1000

    ! Allocate Arrays
    allocate(X(n_steps), z_vals(n_steps), y_vals(n_steps), U_vals(n_steps), q_t(n_steps), dz_vals(n_steps))

    ! Open Output File
    open(12, file='mydata.dat')
    write(12,*) 'Time (s)', 'X (m)', 'z (m)', 'Depth (m)', 'U (m/s)', 'qt (m2/s)'

    ! ***********************
    ! Initialization (Runs Once)
    ! ***********************
    do i = 1, n_steps
        X(i) = (i - 1) * dx
        z_vals(i) = z0 + slope * X(i)  
        y_vals(i) = y0  
        U_vals(i) = Q / y0   
    end do

    ! ***********************
    ! Start Time Loop 
    ! ***********************
    do t_index = 1, n_time_steps  
        print *, "t=", t_index, "Before update: y(1)=", y_vals(1), "y(2)=", y_vals(2)
        ! Start Spatial Loop
        do x_index = 2, n_steps-1

            ! Compute Froude Number
            Fr = U_vals(x_index) / sqrt(g * y_vals(x_index))

            ! Compute energy dissipation (S_f)
            S_f = mann**2 * U_vals(x_index)**2 / y_vals(x_index)**(4.0/3.0)
            
            ! Set known bed slope as S_b 
            S_b = slope

            ! Compute Change in Depth 
            y_vals(x_index) = y_vals(x_index) + dx * ((S_b - S_f) / (1.0 - Fr**2))
       
            ! Apply physical constraints
            if (y_vals(x_index) < 0.0) y_vals(x_index) = 0.0

            ! Compute velocity based on new depth
            U_vals(x_index) = Q / y_vals(x_index)
          
            ! Compute Sediment Transport Rate
            C = 1.0 / mann * sqrt(g / y_vals(x_index))
            q_t(x_index) = 0.05 * (U_vals(x_index)**5) / (C**3 * sqrt(g) * d**2 * ds)
            
            ! Compute Bed Elevation Change (Now Coupled)
            dz_dx = - (q_t(x_index+1) - q_t(x_index)) / ((1.0 - p) * 2 * dx)
           
            z_vals(x_index) = z_vals(x_index) + dz_dx * dt
            
            ! Write to output file
            write(12,*) t_index * dt, X(x_index), z_vals(x_index), y_vals(x_index), U_vals(x_index), q_t(x_index)
      
        end do ! End Spatial Loop

    end do ! End Time Loop  

    print *, "t=", t_index, "Before update: y(1)=", y_vals(1), "y(2)=", y_vals(2)

    ! Close Output File
    close(12)

    ! Deallocate Arrays
    deallocate(X, z_vals, y_vals, U_vals, q_t, dz_vals)

end program my_project

当前问题 #

• 程序无编译/运行报错，但水深(y)、流速(U)、输沙率(q_t)仅随时间更新，沿空间x轴所有位置值完全相同；床面高程(z)可时空变化，但其他变量空间分布无差异。
• 若跳过单独的初始化循环，输出文件无内容。
• 期望输出：100个时间步、4000个空间点，所有变量（z、y、U、q_t）均呈现时空变化——沿x轴有分布差异，且每个时间步各x位置变量向系统平衡更新。

控制方程 #

稳态一维水深方程 #

$$\frac{dy}{dx} = \frac{S_b - S_f}{1 - Fr^2}$$
欧拉离散格式：

y^(n+1) = y^n + dx * ((S_b - S_f)/(1 - Fr^2))_n

Engulend-Hansen输沙率方程 #

$$q^T = \frac{0.05 u^5}{C3 \sqrt{g} \Delta^2 d_s}$$
其中 $C = \frac{1}{n} \sqrt{\frac{g}{y}}$（n为曼宁系数）

Exner床面演变方程 #

$$(1-p)\left(\frac{dz_b}{dt}\right)i = -\frac{q_i^T(t) - q{i,UPWIND}^T(t)}{\Delta x}$$

问题根源与修正方案 #

1. 水深更新的空间依赖缺失 #

当前代码中，水深更新仅使用当前x位置的旧值，未利用上游（或下游）的空间梯度信息——而稳态水深方程本质是沿程的渐变流方程，需要从边界开始逐段递推更新，而非每个位置独立计算。

修正方法：将空间循环改为从上游到下游的递推式更新，固定上游边界条件后逐段计算：

! 在时间循环内替换原空间更新逻辑：
do x_index = 2, n_steps
    ! 用上一位置的已更新值计算当前位置的水力参数
    Fr = U_vals(x_index-1) / sqrt(g * y_vals(x_index-1))
    S_f = mann**2 * U_vals(x_index-1)**2 / y_vals(x_index-1)**(4.0/3.0)
    ! 递推计算当前位置水深
    y_vals(x_index) = y_vals(x_index-1) + dx * ((S_b - S_f) / (1.0 - Fr**2))
    if (y_vals(x_index) < 0.0) y_vals(x_index) = 0.0
    ! 更新流速
    U_vals(x_index) = Q / y_vals(x_index)
end do

2. 输沙率与床面更新的迎风离散错误 #

当前Exner方程使用中心差分，未考虑输沙率的传输方向，会导致数值不稳定或空间分布无变化。需改用迎风差分（根据流速方向选择上游值）。

修正Exner方程的离散逻辑：

! 先计算所有位置的输沙率（水深和流速更新完成后）
do x_index = 1, n_steps
    C = 1.0 / mann * sqrt(g / y_vals(x_index))
    q_t(x_index) = 0.05 * (U_vals(x_index)**5) / (C**3 * sqrt(g) * d**2 * ds)
end do

! 迎风法更新床面高程（假设流速向下游，取上游输沙率）
do x_index = 2, n_steps
    dz_dt = - (q_t(x_index) - q_t(x_index-1)) / ((1.0 - p) * dx)
    z_vals(x_index) = z_vals(x_index) + dz_dt * dt
end do

3. 输出逻辑优化 #

当前代码在空间循环内逐位置写入，未包含边界数据且输出顺序混乱。建议在每个时间步结束后，统一遍历所有空间位置写入：

! 替换原空间循环内的write语句，改为在时间循环末尾执行：
do x_index = 1, n_steps
    write(12,*) t_index * dt, X(x_index), z_vals(x_index), y_vals(x_index), U_vals(x_index), q_t(x_index)
end do

4. 初始化循环的必要性 #

初始化循环是必须的——它为所有空间位置的变量赋予合法初始值（如沿程床面高程、均匀水深），若跳过会导致数组未初始化，引发未定义行为，输出空内容是正常现象，无需修改初始化步骤。

内容的提问来源于stack exchange，提问作者user30026192