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.