# Solving stiff ODEs: Dealing with Jacobian terms which take too long to compute with finite differences

I have a system of PDEs describing atmospheric chemistry and transport. I use finite-differences to make my system of PDEs into a system of ~10,000 ODEs. I then integrate the ODEs forward in time with the CVODE BDF method (from Sundials). For my current discretization of my PDEs, the jacobian is approximately banded (band width = ~200), so I have been using CVODE's direct banded solver.

However, I'd like to couple additional ODEs describing the evolution atmospheric temperature. Computing the right-hand-side to these new ODEs will be slow because it requires computing how infrared light passes through the atmosphere (which is slow). My worry is that It will become impractical to compute the Jacobian terms related to atmospheric temperature. For example, suppose I computed a jacobian term with finite differences:

$$\frac{d(dT/dt)}{df_\mathrm{CO_2}}=\frac{dF_T}{df_\mathrm{CO_2}}\approx \frac{F(f_\mathrm{CO_2} + \Delta f_\mathrm{CO_2}) - F(f_\mathrm{CO_2})}{\Delta f_\mathrm{CO_2}}$$

Here $$f_\mathrm{CO_2}$$ is the concentration of CO$$_2$$, and $$T$$ is temperature. The problem is that computing just this one jacobian term will require computing $$F(f_\mathrm{CO_2} + \Delta f_\mathrm{CO_2})$$, which is incredibly slow (because infrared radiative transfer calculations are slow). Computing all jacobian terms for tens of species would probably take minutes, which isn't acceptable.

One potential solution is to use a split ODE solver such as ARKODE. Split ODE solvers integrate "fast" solution components with implicit methods and "slow" solution components with explicit methods. I could integrate the atmospheric chemistry with implicit methods, and the temperature change with explicit methods. In this case, I would not need to compute the jacobian terms for temperature (because it is explicitly integrated).

My first question: Will split ODE solvers likely work for my problem?

An alternative potential solution is to use CVODE BDF method (fully implicit), and just omit the jacobian terms which would take forever to compute, and hope that the jacobian is a good enough approximation for the non-linear solves to converge. My second question: Will this work?

Any feedback is appreciated.

• I have no experience with the problem you are trying to compute, but it seems that the infrared light transmission is computed at once across the whole domain (via explicit integration of the transmission equations) at every step, right ? If you had a PDE instead for the transmission, you could add it to extend your previous PDE and thus keep a banded structure with maybe a reasonable Jacobian cost at the end ? Still, this may also be quite impractical because of the code changes it involves... Otherwise, maybe you can use an implicit solver with a Jacobian free Newton Krylov solver ! Sep 4 at 19:27
• @Laurent90 I will be able to preserve the banded jacobian when I couple the PDE describing climate. The challenge is going to be to compute some of the banded jacobian terms. This will take too long. Sep 4 at 23:50
• Are you making use of sparsity pattern? Tools like SparsityDetection.jl and SparseDiffTools.jl can be used to get color vectors of large sparse Jacobians and accelereate the internals of solvers using this information. This can give orders of magnitude performance improvements on systems of ~10,000 ODEs. See this tutorial for details. Sep 5 at 20:01
• @ChrisRackauckas Ya I'm using a banded matrix solver. But ya i've seen these Julia tools. Unfortunately my RHS is in Fortran. But maybe someday I'll move t julia and use them. Sep 6 at 17:34

My first question: Will split ODE solvers likely work for my problem?

From your description, this sounds like a textbook use-case for a split ODE solver. Neither an implicit method nor an explicit method is practical across the entire ODE. IMEX methods sound like a reasonable choice as long as the atmospheric temperature dynamics that you want to treat explicitly is actually nonstiff.

If the difference in timescales or evaluation costs between chemistry+transport and the temperature is extremely large, it might be preferable to use the MRIStep module in ARKODE. This is an implementation of multirate infinitesimal GARK methods which use different timesteps for the different partitions.

Atmospheric chemistry is often very stiff. Even though it would be treated implicitly, it has the potential to cause order reduction and reduced efficiency. Both IMEX and multirate infinitesimal methods are susceptible to this. There is not enough information to say if this would actually be an issue for your problem, though.

My second question: Will [inexact Jacobians with BDF] work?

For sufficiently small timesteps, it will converge, but the rate can be slow.

Newton's method, which is used to solve the nonlinear system to compute the next step, is fairly resilient to inexact Jacobians. In the extreme case of using the zero matrix as the Jacobian, Newton's method degenerates into a fixed point iteration. What you suggest is a quasi-Newton method in which some variables are resolved by Newton's method and some are resolved by fixed point iterations. Whether or not this is practical will depend on the stiffness, timestep, and evaluation cost of the temperature dynamics.

BDF is well-suited for these inexact nonlinear solves. I'm not familiar with the internals of CVODE but typically BDF implementations start the nonlinear solver with an initial guess produced by high order extrapolation. This is already pretty accurate. It doesn't take many iterations or precise Jacobians to reach the desired accuracy; they mostly help with the stability.

Conclusion

Probably the easiest think to do is try CVODE with the partial Jacobian because your code is already set up for this. ARKODE seems like a reasonable second choice, but does require extra work to partition the RHS of the ODE. If both are unsatisfactory, you may want to consider Rosenbrock-W methods, which only solve linear systems of equations and accept completely arbitrary Jacobians. This won't require partitioning, but is not part of SUNDIALS as far as I know. There are more exotic techniques such as waveform relaxation and cosimulation that could be used as well.

• Thanks for the in depth reply. I will try just using CVODE BDF first. However, I just tried using the implicit methods in ARKODE to integrate the chemistry + transport (no climate). And ARKDODE failed to complete any integration. Maybe the ODEs are too stiff for ARKODE? Not sure. Sep 4 at 23:52
• That's a bit surprising to hear an implicit RK method failed like this. I would guess the integrator configuration needs adjustment, but it's hard to say without more information. Sep 5 at 2:29
• Not surprising because it's ARKODE. The ARKODE defaults aren't able to handle very stiff equations: you need to tweak a few things manually to normally get that integrator stable. You even see that in the example problems that come with ARKODE itself, see ROBER for example. The pure Julia KenCarp4 and such are more tuned for handling the highly stiff problems seen in chemical dynamics, as shown in the SciMLBenchmarks, so that might be worth a try. Sep 5 at 19:57
• @ChrisRackauckas. I've read some of your benchmarks on ARKODE. Wish ARKODE would just work. I'll try KenCarp4, and see how it goes. I can do this by making a julia wrapper to my Fortran. Although, in the long run, since I'm working in Fortran, it would be nice if I could link to a compiled version of DifferentialEquations.jl, and call KenCarp4 (or other solvers) from the fortran. I'm pretty sure this is impossible though? Sep 6 at 17:40
• @ChrisRackauckas I tried KenCarp4. Seems to work quite a bit better than ARKODE. With ARKODE the integration only got to 10^5 seconds before getting stuck. KenCarp4 gets to 10^11 seconds, and continues but with pretty small time steps. Both pale in comparison to QNDF or CVODE BDF, which easily carry out the integration. Sep 6 at 19:46