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

Question

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.

Answer 1

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.