Bulirsch-Stoer algorithm to solve simple chemistry. Spanning long time intervals after stationarity

Question

Hi all and thank you in advance.

I am working on a time-dependent transport-chemistry model to study the composition of planetary atmospheres.

The equations are the following $$\frac{\partial n(z,t)}{\partial t} = -\frac{\partial \Phi (z,t)}{\partial z} + p(z,t) -n(z,t)l(z,t)$$ for each chemical component, characterized by a density profile $n(z,t)$.

Let's focus in the chemistry part, because the problem I am going to explain does not take place for transport alone (modeled as a thermodynamic diffusion flux $\Phi$). $$\frac{\partial n(z,t)}{\partial t} = p(z,t) -n(z,t)l(z,t).$$

Let's suppose we have three molecules interacting following the following reactions

H2O + O --> 2OH

2OH --> H2O + O

The set of equations to be modeled are

$$\partial n_{H2O}/\partial{t} = \partial n_{O}/\partial{t} = -An_{H2O}n_{O} +Bn_{OH}^{2}\\ \partial n_{OH}/\partial{t} = 2An_{H2O}n_{O} -2Bn_{OH}^{2}$$

The algorithm I am testing is a semi-implicit Bulirsch-Stoer method (see Numerical recipes in Fortran, page 735, for a detailed explanation). I basically adapted the code is provided within the book.

This algorithm is suitable for my project because it makes an error estimation so that time step is constantly adapted. I need to integrate up to million of years in about 1 hour, i.e. my simulation up to a certain point has to integrate with $h > 10^5 s$.

Next I am displaying a screenshot from this "toy model" typical realization (This would be a first step before including hundreds of reactions of the same kind)

The system clearly achieves stationary state, but the time step does not increase beyond than $h \approx 10^{-3}$. The idea is that in a Stiff Systems time step increases as reactions contributions as these become stationary.

So far, I think the core of my problem is the piece of code devoted to determine the error associated for each integrated step

errmax=SMALL
do i=1,nv
  errmax=max(errmax,abs(yerr(i)/yscal(i)))
end do
errmax=errmax/eps

Where what I less understand is the "yscal" factor. In the previous reference (Numerical Recipes in Fortran, page 735) it is mentioned

We now mention an important point: "It is absolutely crucial to scale your variables properly when integrating stiff problems with automatic stepsize adjustment. As in our nonstiff routines, you will be asked to supply a vector $y_{scal}$ with which the error is to be scaled. For example, to get constant fractional errors, simply set $y_{scal} = |y|$. You can get constant absolute errors relative to some maximum values by setting $y_{scal}$ equal to those maximum values. In stiff problems, there are often strongly decreasing pieces of the solution which you are not particularly interested in following once they are small. You can control the relative error above some threshold $C$ and the absolute error below the threshold by setting

$y_{scal} = max(C,|y|)$

If you are using appropriate non-dimensional units, then each component of $C$ should be of order unity. If you are not sure what values to take for $C$, simply try setting each component equal to unity. We strongly advocate the previous expression for stiff problems

"Playing" with such $y_{scal}$ definitely the behavior changes so that I believe including the proper expression for it I can get the desired result. When I include "transport" alone, I obtain the desired evolution in $h$ the time step. Something happens with chemistry.

I am open to any kind of suggestions. Even if you can provide me other ideas to implement. I experimented a bit with Rosenbrock methods with the same exact result. I am really stuck with such model.

Answer 1

Numerical Recipes is well-known as a bad resource for implementations of methods for DiffEqs, and you'll note that the mention that extrapolation methods are generally worse can be found in there (and it's pretty easy to verify on most test equations). So this whole thing somewhat starts off in the wrong direction.

Anyways, if you take a look at Hairer II you'll notice there are a lot things that could be going on. First and foremost, if you're using the standard proportional error control that NR seems to recommend then you'll run into something called "the hump behavior" where the error difference can take a bit to decrease after going past a stiff increase. This is why it's generally recommended that you use a Gustafsson accelerated type of controller for stiff ODE solvers. So maybe that has something to do with it.

Another thing that can be going on is just floating point error. Your numbers are far too large and spread out in order of magnitude to perform operations well without getting a lot of floating point error. You might want to rescale the units so that things are closer to one.

Note that Hairer's dopri5 and dop853 actually display this same behavior that if the error is too low then their timesteps end up becoming constant. Usually this isn't a big deal but there are some test equations which can really show this, and I haven't found the time to investigate why it does this (maybe large exponents do something weird in Fortran?). So maybe you could be seeing something similar.

Or maybe that error scaling is the issue. Honestly, I have never seen that in a production ODE solver or in a book/publication about ODE solver methods from a professional in the field. If you want to do the standard algorithm, it's something like:

qtmp = EEst^(1/(order(alg)+1))/gamma
q = max(inv(qmax),min(inv(qmin),qtmp))
dtnew = dt/q

q is the scaling factor for dt. gamma = 0.9 is the safety factor. qmax is the maximum amplification factor usually taken to be between 2-10, qmin is the minimum de-amplification usually 0.2. It's done in the inversion to make it more numerically stable. EEst is the error estimate, which is:

EEst = norm(err./(abstol .+ max.(abs.(uprev),abs.(u)).*reltol))

where err is the the difference between embedded methods (or whatever estimate you use) the norm is the semi-norm Hairer writes about:

norm(u) = sqrt(sum(u.^2)/length(u)

In this case, EEst < 1 only if the tolerances are satisfied, in which case you reject steps where it's > 1. I can guarantee you that basic solvers are using this, and more advanced solvers are using a play on this with a PI-controller or acceleration.

But anyways, this is probably all besides the point because as stated in the comments, there are well-built ODE solvers so unless you are specifically doing something with the intent on improving them... don't recreate them. In fact, here's a well known stiff extrapolation code in Fortran already:

http://www.unige.ch/~hairer/prog/stiff/seulex.f