Dealing with non-physical negative ODE solutions using ODEPACK

Hi and thank you all again. I am solving the reaction-diffusion-advection equation as follows

$\partial_{t} n\left(t,z\right) = -\partial_{z}\left(\Phi\left(t,z\right)\right) + p\left(t,z\right) -n\left(t,z\right)l\left(t,z\right)$

The terms $p\left(t,z\right)$ and $l\left(t,z\right)$ are computed according to big sets of chemical, and photo-chemical reactions. These reactions are characterized by rates which span huge different scales. Hence the system is stiff. I clarify this point because otherwise a simple Euler-Implicit scheme is perfectly valid for any time step if only transport $\Phi\left(t,z\right)$ is present. Nevertheless the solution absolutely fails involving chemistry if I cannot handle the truncation error. Since I need to integrate for great amounts of times (up to millions of years) to achieve the stationary state, hence I need an adaptive time step scheme. I solve the stationary state by pure time evolution because of further research requirements.

I found the best solution was using ODEPACK. Taking the spatial discretization we end up with a system of ordinary differential equations in which each dependent variable is the density of a certain chemical species at a certain height.

https://computation.llnl.gov/casc/odepack/

Here I attach a publication which basically solves the same kind of system as mine (planetary atmosphere), and uses ODEPACK.

In order to use ODEPACK within my fortran90 code I work with a modern fortran ODEPACK interface to call the old fashion fortran77 ODEPACK routines. In particular the one known as DLSODAR.

I faced the problem that sometimes DLSODAR yields negative densities which does not make sense in this physical context.

My solution was simple. The integration step is made through calling another routine within DLSODAR code, named DSTODA. I introduce some code to just detect if any dependent variable becomes negative

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) statement basically looks for each Y(I) possible and when it is less than zero, the routine returns the last Y with all values greater than zero, altogether with the attempted time step. Both stored before DSTODA call in YOLD HOLD respectively.

After that returned unfinished solution, I reduce the last attempted time step by a factor $1/10$ and proceed from the last step. Furthermore, when from DSTODA routine the time derivative and jacobians are called, if some density is negative, those auxiliar routines do not change the previous computed time derivative and jacobian. I found this important because sometimes DSTODA may yield a positive integrated solution but not correct. Hence not changing the time derivative in any case will change the solution's sign. I think this is due to the fact that sometimes the time derivative or jacobian are computed with intermediate Y solutions that not being positive may get a positive one.

--------------------------------------------------------------------------

It works, except for a certain kind of reactions. The physical origin of the equation does not care because as other reactions is characterized by a forward and reverse reaction rates. Anyway the forward mechanism has $k_{f} \approx {10^{-24},10^{-12}} cm^{3} s^{-1}$ and reverse rate $k_{r} \approx {10^{-38},10^{-27}} s^{-1}$. I will give the particular example as an independent simulation.

The reaction is as follows, considering we have D = 500 height layers (index j) and nElems = 3 (index i) chemical species

y(1+(j-1)*nElems) + y(2+(j-1)*nElems) <---> y(3+(j-1)*nElems)

And the code

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')
close(10)

! Importing intial water density profile
n(:,:) = 0.0d0
open(unit=10,file='inProf.dat')
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')
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')
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

And part of the output is as follows (I removed from the source code some prints* for the sake of being clearer)

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

Basically the time step is stuck between 1e-3 s and 1e-1 s. It is not worth to wait more than a day just to see if ODEPACK finally decides to increase dt. I tried changing tolerances too. But the biggest issue is that the solution is wrong too. As an example a t = 100 s I display the simulation output In the red solution (corresponding to Y(3)) those peaks are wrong. The initial condition is taken from a well-mixed hydrostatic atmosphere with a given composition.

As an extension, if we let DLSODAR to run on its own, without any modifications at al, we have the following output, for t = 1e12 s Where the white discontinuity is how gnuplot deals when we have negative values on a logarithmic scale. I should not accept negative solutions at any time value. This output only takes a few seconds to be computed.

The final (hopefully correct) stationary state solution if anyway we let DLSODAR to keep computing no matter negative numbers looks like Any suggestions when ODEPACK yields negative non physical solutions?. Is this a known problem common to occur?. I feel the problem should be so easy to solve, but my mind cannot see it clearer.

Again, thank you so much.

• Does the original equation only admit "physical" solutions? – nicoguaro Jan 7 '18 at 17:45
• Are the negative values comparable in magnitude to the absolute error tolerance? If so, then it's just the solver not being able to distinguish values less than abstol apart ($-\epsilon\approx0$), and I'd say it's pretty normal. – Kirill Jan 7 '18 at 20:33
• Even if the negative values are larger than tolerance, it's worth pointing out that in general, even if the exact ODE guarantees that the solution will remain positive, there is no such guarantee for the (time-)*discretized* ODE. That's just a fact of life. – Wolfgang Bangerth Jan 7 '18 at 22:36
• There may be ways to guarantee positivity of the numerical solution, but they don't involve using ODEPACK, so that's probably not what you're looking for. – David Ketcheson Jan 8 '18 at 4:00
• If you add @ before the username, people will be notified. – nicoguaro Jan 8 '18 at 14:40