1
$\begingroup$

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')
  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

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

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

enter image description here

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

enter image description here

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.

$\endgroup$
  • 2
    $\begingroup$ Does the original equation only admit "physical" solutions? $\endgroup$ – nicoguaro Jan 7 '18 at 17:45
  • 1
    $\begingroup$ 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. $\endgroup$ – Kirill Jan 7 '18 at 20:33
  • 5
    $\begingroup$ 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. $\endgroup$ – Wolfgang Bangerth Jan 7 '18 at 22:36
  • $\begingroup$ 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. $\endgroup$ – David Ketcheson Jan 8 '18 at 4:00
  • $\begingroup$ If you add @ before the username, people will be notified. $\endgroup$ – nicoguaro Jan 8 '18 at 14:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.