嗨，再次感谢大家。我正在求解反应-扩散-平流方程如下
和项是根据大量化学和光化学反应计算得出的。这些反应的特征在于跨越巨大不同尺度的速率。因此系统是僵硬的。我澄清了这一点，因为如果仅存在传输那么简单的欧拉隐式方案对于任何时间步都是完全有效的。然而，如果我无法处理截断错误，那么涉及化学的解决方案绝对会失败。由于我需要积分大量时间（长达数百万年）才能达到静止状态，因此我需要一个自适应时间步长方案。由于进一步的研究需要，我通过纯时间演化来解决稳态。
我发现最好的解决方案是使用 ODEPACK。采用空间离散化，我们最终得到一个常微分方程系统，其中每个因变量是某个化学物质在某个高度的密度。
https://computation.llnl.gov/casc/odepack/
在这里我附上一个出版物，它基本上解决了与我的系统（行星大气）相同的系统，并使用了ODEPACK。
为了在我的 fortran90 代码中使用 ODEPACK，我使用现代的 fortran ODEPACK 接口来调用老式的 fortran77 ODEPACK 例程。特别是被称为 DLSOAR 的那个。
我面临的问题是，有时 DLSODAR 会产生负密度，这在这种物理环境中是没有意义的。
我的解决方案很简单。集成步骤是通过调用 DLODAR 代码中名为 DSTODA 的另一个例程来完成的。我介绍了一些代码来检测是否有任何因变量变为负数
C-----------------------------------------------------------------------
C CALL DSTODA(NEQ,Y,YH,NYH,YH,EWT,SAVF,ACOR,WM,IWM,F,JAC,DPRJA,DSOLSY)
C-----------------------------------------------------------------------
DO 295 I = 1,N
295 YOLD(I) = Y(I)
TOLD = TN
HOLD = H
SIGNVAL = .TRUE.
CALL DSTODA (NEQ, Y, RWORK(LYH), NYH, RWORK(LYH), RWORK(LEWT),
1 RWORK(LSAVF), RWORK(LACOR), RWORK(LWM), IWORK(LIWM),
2 F, JAC, DPRJA, DSOLSY)
CALL POSITIVENESS (SIGNVAL,NEQ,Y) ! (30 - 11 - 2017)
IF(SIGNVAL.EQV..FALSE.) KFLAG = - 3
KGO = 1 - KFLAG
GO TO (300, 530, 540, 400), KGO ! IF KFLAG = -3, solution is negative and DLSODAR returns
! from the last non negative solution to repeat if desired
! with a smaller time step H. (14 - 12 - 2017)
C-----------------------------------------------------------------------
CALL POSITIVENESS(SIGNVAL,NEQ,Y)语句基本上查找每个Y(I)可能的值，当它小于零时，例程返回Y所有值大于零的最后一个，以及尝试的时间步长。两者分别在DSTODA调用之前存储YOLD HOLD。
在返回未完成的解决方案之后，我将最后尝试的时间步减少并从最后一步开始。此外，当从例程调用时间导数和雅可比时，如果某个密度为负，则这些辅助例程不会更改先前计算的时间导数和雅可比。我发现这很重要，因为有时 DSTODA 可能会产生积极的集成解决方案但不正确。因此，在任何情况下不改变时间导数都会改变解的符号。我认为这是因为有时时间导数或雅可比是用中间解决方案计算的，不是正数可能会得到正数。DSTODAY
-------------------------------------------------- ----------------------
它有效，除了某种反应。方程的物理起源并不关心，因为其他反应的特征是正向和反向反应速率。无论如何，正向机制有和反向速率。我将给出特定示例作为独立模拟。
反应如下，考虑到我们有 D = 500 个高度层（索引 j）和 nElems = 3（索引 i）化学物质
y(1+(j-1)*nElems) + y(2+(j-1)*nElems) <---> y(3+(j-1)*nElems)
和代码
use dlsodar_f95
implicit none
integer i,j,jj,u, D, nElems
integer lrw, liw, MU, ML
double precision p(500),kf(500),kr(500)
double precision n(3,500),sol_f(3*500),dummy1(7),dummy2(2)
double precision sol(3*500)
double precision t, xend, dt, dti, tf
double precision Rtol, Atol(3*500)
character(len=1024) forward, reverse
D = 500
nElems = 3
! Importing reaction constants
open(unit=10,file='reaction.dat')
read(10,*) p(:)
read(10,*) forward, kf(:)
read(10,*) reverse, kr(:)
close(10)
! Importing intial water density profile
n(:,:) = 0.0d0
open(unit=10,file='inProf.dat')
read(10,*)
do i = 1,D
read(10,*) dummy1(:), n(3,i), n(1,i), dummy2(:), n(2,i)
end do
close(10)
! Solving using DLSODAR
t = 0.0d0
dt = 0.0d0
Rtol = 1e-1
Atol(:) = 1e-2
MU = 3
ML = 3
lrw = 22+(10 + (2*(MU)+ ML))*nElems*D
liw = 20+nElems*D
tf = 1e1
dti = 1d-15 ! dti is the trial initial time step
u = 0
do i = 1,D
do j=1,nElems
u = u + 1
sol(u) = n(j,i)
end do
end do
do while (t.lt.tf)
xend = tf
call ODE_DLSODAR(t,xend,sol,derivs,jac,
* DLSODAR_dummy_Constraint,0,
* xend,sol_f,dti,dt,nElems,nElems,liw,lrw,
* relative_tolerance=RTOL,
* absolute_tolerances=ATOL,
* max_steps=int(1e9), ! If max_steps is too big, maybe there is an allocation error
* store_steps=.false.,verbose_level=1)
!
! If solution is negative this loop will repeat from last non negative
! iteration but with a smaller time step
!
t = xend
if (t.lt.tf) then
dti = dt/10.
u = 0
do i = 1,D
do j=1,nElems
u = u + 1
sol(u) = sol_f(u)
end do
end do
end if
end do
u = 0
do i = 1,D
do j=1,nElems
u = u + 1
n(j,i) = sol_f(u)
end do
end do
open(unit=10,file='finProf.dat')
do i = 1,D
write(10,*) p(i), n(1,i), n(2,i), n(3,i),
* kf(i),kr(i)
end do
close(10)
contains
c --------------------------------------------------------------
subroutine jac(NEQ, X, Y, ML, MU, PD, NROWPD)
implicit none
! Arguments
integer, intent(in) :: NEQ,ML,MU,NROWPD
double precision, intent(in) :: X
double precision, dimension(NEQ), intent(in) :: Y
double precision, dimension(NROWPD,NEQ), intent(out) :: PD
! Local variables
integer i, j,u
double precision kf(500),kr(500)
character(len=1024) forward, reverse
do u = 1, NEQ
if (Y(u).NE.Y(u)) then
print*, "yields NAN"
stop
end if
if (Y(u).lt.0.0d0) then
return
end if
end do
PD(:,:) = 0.0d0
! Importing reaction constants
open(unit=10,file='reaction.dat')
read(10,*)
read(10,*) forward, kf(:)
read(10,*) reverse, kr(:)
close(10)
print*, "dentro jac:",x
do i=1,500
PD(2,3+(i-1)*3) = kr(i)
PD(3,2+(i-1)*3) = -kf(i)*Y((i-1)*3+1)
PD(3,3+(i-1)*3) = kr(i)
PD(4,1+(i-1)*3) = -kf(i)*Y((i-1)*3+2)
PD(4,2+(i-1)*3) = -kf(i)*Y((i-1)*3+1)
PD(4,3+(i-1)*3) = -kr(i)
PD(5,1+(i-1)*3) = -kf(i)*Y((i-1)*3+2)
PD(5,2+(i-1)*3) = kf(i)*Y((i-1)*3+1)
PD(6,1+(i-1)*3) = kf(i)*Y((i-1)*3+2)
end do
return
end subroutine jac
c --------------------------------------------------------------
subroutine derivs(NEQ,X,Y,DYDX)
implicit none
integer, intent(in) :: NEQ! number of variables
double precision, intent(in) :: X! independent variable, time
double precision, dimension(NEQ), intent(in) :: Y ! dependent variables. Unknowns algorithm is solving
double precision, dimension(NEQ), intent(out) :: DYDX ! time derivative
integer u ! running through nElems*D positions
integer i, j, counter
double precision dummy(2) ! dummy variable
double precision kf(500),kr(500)
character(len=1024) forward, reverse
do u = 1, NEQ
if (Y(u).NE.Y(u)) then
print*, "yields NAN"
stop
end if
if (Y(u).lt.0.0d0) then
return
end if
end do
DYDX(:) = 0.0D0
! Importing reaction constants
open(unit=10,file='reaction.dat')
read(10,*)
read(10,*) forward, kf(:)
read(10,*) reverse, kr(:)
close(10)
print*, "dentro dydx:",x
do i=1,500
DYDX((i-1)*3+1) = -kf(i)*Y((i-1)*3+1)*Y((i-1)*3+2)
* + kr(i)*Y((i-1)*3+3)
DYDX((i-1)*3+2) = -kf(i)*Y((i-1)*3+1)*Y((i-1)*3+2)
* + kr(i)*Y((i-1)*3+3)
DYDX((i-1)*3+3) = kf(i)*Y((i-1)*3+1)*Y((i-1)*3+2)
* - kr(i)*Y((i-1)*3+3)
end do
return
end subroutine derivs
c --------------------------------------------------------------
end
部分输出如下（为了更清楚，我从源代码中删除了一些打印*）
dentro dydx: 4.1576048498046223
dentro dydx: 4.1577169032470946
dentro dydx: 4.1577169032470946
dentro dydx: 4.1578289566895670
dentro dydx: 4.1578289566895670
dentro dydx: 4.2698823991657884
dentro dydx: 4.2698823991657884
dentro dydx: 4.1858423173086221
dentro dydx: 4.1858423173086221
Y(I) : -1.4252481539273311E-013
Where I : 25
Ha salido negativo
el xend alcanzado es : 4.1578289566895670
pasos temporales, dti : 1.1205344247195033E-004 , dt : 0.11205344247622159
dentro dydx: 4.1578289566895670
dentro dydx: 4.1690343009371889
dentro dydx: 4.1690343009371889
Y(I) : -8.0043582938772017E-007
Where I : 4
Ha salido negativo
el xend alcanzado es : 4.1578289566895670
pasos temporales, dti : 1.1205344247622160E-002 , dt : 1.1205344247622160E-002
dentro dydx: 4.1578289566895670
dentro dydx: 4.1589494911143294
dentro dydx: 4.1589494911143294
Y(I) : -1.0391305382066468E-010
Where I : 4
Ha salido negativo
el xend alcanzado es : 4.1578289566895670
pasos temporales, dti : 1.1205344247622159E-003 , dt : 1.1205344247622159E-003
dentro dydx: 4.1578289566895670
dentro dydx: 4.1579410101320429
dentro dydx: 4.1579410101320429
dentro dydx: 4.1580530635745188
dentro dydx: 4.1580530635745188
dentro dydx: 4.2701065061024490
dentro dydx: 4.2701065061024490
dentro dydx: 4.1860664242065013
dentro dydx: 4.1860664242065013
Y(I) : -1.0170965354084183E-013
Where I : 25
Ha salido negativo
el xend alcanzado es : 4.1580530635745188
pasos temporales, dti : 1.1205344247622159E-004 , dt : 0.11205344252793006
基本上时间步长在 1e-3 s 和 1e-1 s 之间。只看 ODEPACK 是否最终决定增加 dt 是不值得等待超过一天的。我也尝试改变公差。但最大的问题是解决方案也是错误的。以 = 100 s 为例，我显示模拟输出
在红色解决方案（对应于 Y(3)）中，这些峰是错误的。初始条件取自具有给定成分的充分混合的流体静力学气氛。
作为扩展，如果我们让 DLSODAR 自行运行，无需任何修改，我们将得到以下输出，对于 t = 1e12 s
当我们在对数刻度上有负值时，白色不连续性是 gnuplot 的处理方式。我不应该接受任何时间价值的负面解决方案。计算此输出仅需几秒钟。
最终的（希望是正确的）静止状态解决方案，如果无论如何我们让 DLSODAR 继续计算，无论负数看起来像
当 ODEPACK 产生负面的非物理解决方案时有什么建议吗？这是一个常见的已知问题吗？. 我觉得这个问题应该很容易解决，但我的头脑无法更清楚地看到它。
再次，谢谢你。