Fortran 中的变量识别问题

Posted 2023-03-27

【中文标题】Fortran 中的变量识别问题

【英文标题】：Problem with variable recognition in Fortran

【发布时间】：2019-07-29 04:57:47

【问题描述】：

我是 Fortran 的初学者，正在编写代码，通过使用微分方程求解器在不同时间步求解微分方程来求解动力学机制。 这是一个下载整个项目的 zip 文件的链接：

http://www.filedropper.com/fortranmicropyrolyzersetup

微分求解器的输入变量在代码中定义如下：

     ! Declaration of variables
     implicit none
     EXTERNAL                   ::  FEXSB_AUTO, JEX_SB
     integer                    ::  neq,Mf,lrw,liw,iwork,itol,itask,istate,iopt     !solver parameters                
     integer                    ::  j,jjk,m                                         !Counters
     double precision           ::  ATOL,RTOL, RWORK                                !solver parameters
     double precision           ::  T, TOUT                                         !starting time (s), timestep time (s)
     double precision           ::  Y, w                                            !molar fraction of biomass (-), mass fraction of biomass (-)
     double precision           ::  y_gas, w_gas, Conc                              !molar fraction in gas phase (-), concentration in gas phase (mol/m³), mass fraction in gas phase (-)
     double precision           ::  n(speciescount)                                 !number of moles (used temporarily to calculate initial molar fracions)  
     character*5                ::  simulationNumber


  !Setting the solver parameters
     neq = SpeciesCount         !number of equations
     ITOL = 1                   !RTOL and ATOL are integers
     RTOL = 1.0D-8               !Relative tolerance
     ATOL = 1.0D-15              !Absolute tolerance
     ITASK = 1                  
     ISTATE = 1
     LRW = 22 +  9*SpeciesCount + 2*SpeciesCount**2                  !Array sizing (see VODE)
     LIW = 30+SpeciesCount                                         !Array sizint (see VODE)
     MF = 22                    !Use BDF with finite difference Jacobian
     IOPT = 1                   !Optional input specified
     Iwork(6) = 7000            !Increase maximum iteration steps from 500 to 2000 otherwise the solver does not converge

然后在每个时间步的 do 循环中调用微分求解器，如下所示：

  ! Solve reactor equations (see FEXSB) and advance time, until stop criterium is met
     TimeStep = 0   ! dimensionless time step
     !DO while(wm(1).lt.0.9999) -> previously used stop criterium
     DO while((y_gas(1).lt.0.99999).OR.(TimeStep.lt.10))        !stop criterium: molar fraction Helium = 0.9999 AND do at least 100 timesteps
         TimeStep = TimeStep + 1
         write(*,*)"Doing iteration for time step",TimeStep
         call DVODE(FEXSB_AUTO,NEQ,Y,T,TOUT,ITOL,RTOL,ATOL,ITASK,ISTATE,IOPT,RWORK,LRW,IWORK,LIW,JEX_SB,MF)     !solve reactor equations to Y
        !CALL DVODE(FEXSB_AUTO,NEQ,Y,T,TOUT,ITOL,RTOL,ATOL,ITASK,ISTATE,IOPT,RWORK,LRW,IWORK,LIW,JEX_SB,MF,RPAR,IPAR)

        ! calculate w, y_gas, w_gas and Conc from Y
        do j = 1, SpeciesCount

            if (Y(j) .lt. 1.0D-10) then      ! round to zero below 10e-10 to avoid negative numbers and numerical problems with small numbers
                Y(j) = 0.0D0
            endif
            if (MolarFlowRate(j) .lt. 1.0D-10) then
                MolarFlowRate(j) = 0.0D0
            endif
            w(j) = Y(j)*n0_tot*amms(j) / (mass_sample)
            y_gas(j) = MolarFlowRate(j) / TotalMolarFlowRate
            w_gas(j) = MolarFlowRate(j)*amms(j) / (1000*MassFlowRate) ! factor 1000 to put molar mass in kg/mol instead of g/mol
            Conc(j) = MolarFlowRate(j) / VolumetricFlowRate

        end do

微分方程求解器没有成功完成行的do循环：

DO while((y_gas(1).lt.0.99999).OR.(TimeStep.lt.10))

当在该行中调用 ODE Solver 时，变量似乎在第一个时间步被识别：

call DVODE(FEXSB_AUTO,NEQ,Y,T,TOUT,ITOL,RTOL,ATOL,ITASK,ISTATE,IOPT,RWORK,LRW,IWORK,LIW,JEX_SB,MF)

求解器成功完成了第一次迭代。时间步长再增加一次后，调用函数再次起作用，但这一次变量以某种方式无法识别。当我在第一次步骤后停止代码调试错误时，我意识到 DVODE 所需的变量不再具有设置值，不知何故它们在第一次成功迭代后被删除。

什么可能导致这个问题？

*DECK DVODE
      SUBROUTINE DVODE (F, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL, ITASK,
     1            ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JAC, MF,
     2            RPAR, IPAR)
      EXTERNAL F, JAC
      DOUBLE PRECISION Y, T, TOUT, RTOL, ATOL, RWORK, RPAR
      INTEGER NEQ, ITOL, ITASK, ISTATE, IOPT, LRW, IWORK, LIW,
     1        MF, IPAR
      DIMENSION Y(*), RTOL(*), ATOL(*), RWORK(LRW), IWORK(LIW),
     1          RPAR(*), IPAR(*)
!-----------------------------------------------------------------------
c dvode: Variable-coefficient Ordinary Differential Equation solver,
! with fixed-leading-coefficient implementation.
! This version is in double precision.
!
! DVODE solves the initial value problem for stiff or nonstiff
! systems of first order ODEs,
!     dy/dt = f(t,y) ,  or, in component form,
!     dy(i)/dt = f(i) = f(i,t,y(1),y(2),...,y(NEQ)) (i = 1,...,NEQ).
! DVODE is a package based on the EPISODE and EPISODEB packages, and
! on the ODEPACK user interface standard, with minor modifications.
!-----------------------------------------------------------------------
! Authors:
!               Peter N. Brown and Alan !. Hindmarsh
!               Center for Applied Scientific Computing, L-561
!               Lawrence Livermore National Laboratory
!               Livermore, CA 94551
! and
!               George D. Byrne
!               Illinois Institute of Technology
!               Chicago, IL 60616
!-----------------------------------------------------------------------

任何帮助将不胜感激。请注意，我是 Fortran 新手。如果我应该提供任何其他信息来帮助您回答我的问题，请随时告诉我。提前致谢！

【问题讨论】：

问题出在哪里还不完全清楚。您观察到调试器未能在第二次迭代期间向您显示变量的值。但是您没有告诉我们代码无法解决问题或代码没有运行完成。如果“问题”只是你有一个损坏的调试器，那就再买一个。如果这不是问题，请更明确地说明问题。

我编辑了我的问题，希望这能让您更清楚问题。

如果您能展示一个可编译和可测试的示例并向我们展示这些参数在第一次调用后的值，那就太好了。 DVODE 在文档中是否提到它们中的任何一个都应该被销毁？如果不是，则可能是堆栈损坏。

我在我的问题中添加了一个包含整个项目的 zip 文件的链接，这样您就可以更好地了解确切的问题以及如何解决它。感谢您迄今为止的所有帮助！

请不要在这里这样做。制作适合此处的minimal reproducible example，不要将 zip 文件上传到可能在几个月或几年内失效的外部链接。这里的问题和答案应该对未来的访问者有用。此外，要求研究大型项目的人太多了。编写一个在循环中调用DVODE 的简短测试程序就足够了。它不应该超过，比如说，50 行，可能更少。

【参考方案1】：

第一次迭代后参数应该保持不变的假设是错误的。

在 fortran 函数中，参数总是通过引用传递。因此，fortran 函数永远不能保证您的原始参数集在调用后是相同的。更改参数的功能可能是有意的（或糟糕的）设计。只有文档可以帮助您，而您在问题中提供的内容还不够，这很可能也意味着：可能没有足够的文档用于此案例/功能。

与 fortran 中的 C/C++ 相比，没有用于变量的“const”修饰符可以保证您现在必须假设的内容。

对我来说，似乎唯一的解决方案是在每次调用 DVODE 之前重新初始化所有参数。

【讨论】：

当然不是总是。 “Real Programmers 认可的唯一参数传递机制是按值返回的调用，正如在 IBM/370 Fortran G 和 H 编译器中实现的那样。” 还有更多可能性。

'C/C++ in fortran there's no "const"'我们有intent(in)和value。