Fortran 90迭代求解偏微分方程：内嵌计算不收敛，子程序调用正常

Hey there, let's break down why your inlined iteration code isn't converging while the subroutine version works. Looking at the code snippet you shared, there's a clear red flag with how you're handling array indices and loop structure—let's dig into it.

First, the Critical Bug in Your Snippet #

Take a look at this part:

DO IT = 2,ITMAX
  DO I = 1,IMAX
    PHIN(IT-1,I,1) = PHIN(IT-1,I,2) - (Y(2) - Y(1))*UINF*PHIY(I)
  END DO
  PHIN(IT,I,1) = PHIN(IT-1,I,1)  ! <-- Here's the problem!
DO...

After the inner I loop finishes, the value of I is IMAX + 1 (Fortran increments the loop variable after the last iteration). That means when you assign PHIN(IT,I,1), you're either:

• Writing to an index beyond your array bounds (if IMAX is the upper limit of the I dimension), which causes undefined behavior, or
• Only updating the very last (out-of-range or unintended) element of PHIN(IT,:,1) instead of every I from 1 to IMAX.

In your working subroutine version, I'd bet this assignment is inside the inner I loop, so it updates every element for each I—that's why the subroutine converges correctly.

Fixes to Get Your Inlined Code Working #

Let's fix the loop structure first:

1. Move the PHIN(IT,I,1) Assignment Inside the I Loop #

This ensures you update every I index for each time step IT:

DO IT = 2,ITMAX
  DO I = 1,IMAX
    ! Update the previous time step's boundary value
    PHIN(IT-1,I,1) = PHIN(IT-1,I,2) - (Y(2) - Y(1))*UINF*PHIY(I)
    ! Now copy the updated value to the current time step
    PHIN(IT,I,1) = PHIN(IT-1,I,1)
  END DO
  ! Rest of your iteration logic goes here...
END DO

2. Double-Check Variable Scope and Update Order #

Another common issue when inlining vs using subroutines is variable persistence:

• In a subroutine, local variables are reinitialized each call (unless saved), while inlined code uses variables from the parent scope. Make sure you're not accidentally overwriting values that should be preserved across iterations.
• Verify that all array updates follow the correct PDE stencil order—for example, if you need to use old values for all neighbors before updating, inlining might have you overwriting a value too early (subroutines often encapsulate this to avoid overwriting).

3. Add Bounds Checking for Debugging #

Fortran can help catch array index errors if you enable bounds checking. Compile with flags like:

• For gfortran: -fbounds-check
• For Intel Fortran: -check bounds

This will throw an error if you're accessing PHIN with I = IMAX + 1, confirming our initial suspicion.

Confirm the Fix #

After adjusting the loop structure, run your code again. The convergence should match the subroutine version because you're now correctly updating every element of PHIN(IT,:,1) instead of a single invalid index.

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