Apologies for the tongue-in-cheek title.

According to this discussion,

native QNDF julia solvers now consistently outperform CVODE_BDF on the test suite. However, my experience solving reasonably stiff sets of ~100’s of ODEs that arise from method of lines problems is that CVODE_BDF is typically faster than native Julia solvers, and often more stable. To give some examples, here’s a simple method of lines solution of the single diffusion-reaction equation:

\frac{\partial c}{\partial t} = D \frac{\partial^2 c}{\partial x^2} - kc

with fixed Dirichlet boundary conditions at both ends.

```
using DifferentialEquations
using LinearAlgebra
using Sundials
using BenchmarkTools
# Code Method of Lines Solver for Film Theory Model
# This code is modified from code published by Moore et al. in the SI of 10.1039/D2EE03237F.
#Inputs:
# - n, number of species, integer
# - R, function for reaction rate. Function input is local concentrations (mol/m3). Function output is rate of change of local concentrations (mol/m3.s). Input and ouput are both vectors of length n.
# - D, vector of diffusivities (m2/s). Vector of length n
# - c0, Initial concentrations of all species, mol/m3.
# - csurf, surface concentration of species 1 (dissolving gas, most often CO2), mol/m3
# - cbulk, concentrations of all species in bulk. Must satisfy R(cbulk) = 0, otherwise an error will be returned.
# - N, number of cells to divide domain into (non-uniform mesh is used)
# - L, length of the domain. Typically about 100 μm.
# - tfinal, the endpoint of the time integration. This should be infinite, but in practice a very long time (e.g. 1e5 seconds) is sufficient.
# - OM, number of orders of magnitude of a over which to grow cells
#Outputs:
# - Solution object to diffusion-reaction problem at all times.
function AbsorptionFilmTheorySolver(n, R, D::Array{Float64}, c0, csurf, cbulk; N = 50, L = 300e-6,
tfinal = 100.0, OM = 2.0, solver = CVODE_BDF())
#Check that the bulk conditions are at equilibrium
if norm(R(cbulk)) > 1e-6
error("Bulk conditions are not at equilibrium.")
end
#Preliminary Calculations
#Boundary Points of Cells
b = 10 .^(LinRange(0.0, OM, N+1))
b = b .- b[1]
b = b .* (L/b[end])
#Cell Center Points & Cell Widths
x = (b[2:N+1] .+ b[1:N]) ./ 2
Δb = b[2:N+1] .- b[1:N]
Δx = x[2:N] .- x[1:N-1]
#Set Initial Conditions (average values within cells)
c0_array = zeros(N,n)
for i = 1:N
for j = 1:n
c0_array[i, j] = c0[j]
end
end
#ODE Function
function odefun!(dcdt, c, (JLeft, JRight, n, N, csurf, cbulk, D, R, Δb, Δx), t)
#Calculate dcdt for cell adjacent to gas-liquid interface
#Start with CO2 at gas-liquid interface
JLeft = -D[1]*(c[1,1] - csurf)/(Δb[1]/2)
JRight = -D[1]*(c[2, 1] - c[1, 1])/Δx[1]
dcdt[1, 1] = ((JLeft) - (JRight))/Δb[1] + R(@view c[1, :])[1]
#Then do the same for all other species at the gas-liquid interface
if n > 1
for j = 2:n
JLeft = 0.0
JRight = -D[j]*(c[2, j] - c[1, j])/Δx[1]
dcdt[1, j] = ((JLeft) - (JRight))/Δb[1] + R(@view c[1, :])[j]
end
end
#Calculate dcdt in cell adjacent to end of domain
for j = 1:n
JLeft = -D[j]*(c[N, j] - c[N-1, j])/Δx[end]
JRight = -D[1]*(cbulk[j] - c[N, j])/(Δb[N]/2)
dcdt[N, j] = ((JLeft) - (JRight))/Δb[end] + R(@view c[N, :])[j]
end
#Calculate rate of change of concentration in interior cells
for i = 2:N-1
for j = 1:n
JLeft = -D[j]*(c[i, j] - c[i-1, j])/Δx[i-1]
JRight = -D[j]*(c[i+1, j] - c[i, j])/Δx[i]
dcdt[i, j] = ((JLeft) - (JRight))/Δb[i] + R(@view c[i, :])[j]
end
end
end
#Solve ODE
tspan = (0.0, tfinal)
prob = ODEProblem(odefun!, c0_array, tspan, (0.0, 0.0, n, N, csurf, cbulk, D, R, Δb, Δx))
sol = solve(prob, solver)
return (x, sol)
end
n, D, c0, csurf, cbulk, L, N, tfinal = (1, [1e-9], [0.0], 10.0, [0.0], 100e-6, 50, 3.0)
R(c) = [-c[1]]
@btime x, sol = AbsorptionFilmTheorySolver(n, R, D, c0, csurf, cbulk; N = N, L = L,
tfinal = 100.0, OM = 2.0, solver = CVODE_BDF());
@btime x, sol = AbsorptionFilmTheorySolver(n, R, D, c0, csurf, cbulk; N = N, L = L,
tfinal = 100.0, OM = 2.0, solver = Rosenbrock23());
@btime x, sol = AbsorptionFilmTheorySolver(n, R, D, c0, csurf, cbulk; N = N, L = L,
tfinal = 100.0, OM = 2.0, solver = QNDF());
@btime x, sol = AbsorptionFilmTheorySolver(n, R, D, c0, csurf, cbulk; N = N, L = L,
tfinal = 100.0, OM = 2.0, solver = QBDF());
```

Here, CVODE_BDF outperfoms all alternatives, though QNDF and QBDF are reasonably close,

```
1.744 ms (25720 allocations: 1.74 MiB)
21.744 ms (443587 allocations: 33.78 MiB)
2.655 ms (35104 allocations: 2.47 MiB)
2.903 ms (38038 allocations: 2.66 MiB)
```

If we consider a more complex system (5 simultaneous coupled PDEs with much greater disparity in reaction time-scales, presumably resulting in substantially more stiffness), we run into a different problem: CVODE_BDF is fast and accurate, Rosenbrock23() is slow and accurate, but QNDF and QBDF fail to converge.

```
#Function to solve complete chemistry.
#Physical data references given in SI of manuscript.
#Species 1 - 5 are CO2, OH-, HCO3-, CO32- and H+, respectively.
function FilmTheoryFullChemistry(; cCO2surf = 34, pHbulk = 10, cHCO30 = 1000,
N = 50, OM = 2, L = 100e-6, tfinal = 1e5, solver = CVODE_BDF())
#Physical Parameters
n = 5
D = [1.91e-9, 5.29e-9, 1.185e-9, 0.92e-9, 9.311e-9] #m2/s
k1f = 8.42 #m3/mol.s
k1r = 1.97e-4 #1/s
k3f = 1.254e6 #1/s
k3r = 6e6 #m3/mol.s
k5f = 2.3e7 #L/mol.s
k5r = 2.3e-1 #mol/m3.s
k2f = 3.71e-2 #1/s
k2r = k2f*(k1r/k1f)*(k5f/k5r) #m3/mol.s - Calculated Exactly to Maintain Consistency
k4f = 1.242e9 #m3/mol.s
k4r = k4f*(k3r/k3f)*(k5r/k5f) #1/s - Calculated Exactly to Maintain Consistency
#Reaction Function
function R(c)
return [-k1f*c[1]*c[2] + k1r*c[3] - k2f*c[1] + k2r*c[3]*c[5],
-k1f*c[1]*c[2] + k1r*c[3] + k3f*c[4] - k3r*c[2]*c[3] - k5f*c[2]*c[5] + k5r,
k1f*c[1]*c[2] - k1r*c[3] + k2f*c[1] - k2r*c[3]*c[5] + k3f*c[4] - k3r*c[2]*c[3] + k4f*c[4]*c[5] - k4r*c[3],
- k3f*c[4] + k3r*c[2]*c[3] - k4f*c[4]*c[5] + k4r*c[3],
k2f*c[1] - k2r*c[3]*c[5] - k4f*c[4]*c[5] + k4r*c[3] - k5f*c[2]*c[5] + k5r]
end
#Calculate Bulk Concentrations
cOHbulk = 10.0^(pHbulk - 14)*1000
cHbulk = k5r/(k5f*cOHbulk)
cHCO3bulk = 2*cHCO30*k3f/k3r/(cOHbulk + 2k3f/k3r)
cCO32bulk = 2*cHCO30*cOHbulk/(cOHbulk + 2k3f/k3r)
cCO2bulk = k1r*cHCO3bulk/(k1f*cOHbulk)
cbulk = [cCO2bulk, cOHbulk, cHCO3bulk, cCO32bulk, cHbulk]
c0 = cbulk
#Solve System of Equations
x, sol = AbsorptionFilmTheorySolver(n, R, D, c0, cCO2surf, cbulk;
N = N, L = L, tfinal = tfinal, OM = OM, solver = solver)
ceq = sol(tfinal)
return (x, ceq, sol)
end
@btime x, ceq, sol = FilmTheoryFullChemistry(;solver = CVODE_BDF());
@btime x, ceq, sol = FilmTheoryFullChemistry(;solver = Rosenbrock23());
x, ceq, sol = FilmTheoryFullChemistry(;solver = QNDF());
x, ceq, sol = FilmTheoryFullChemistry(;solver = QBDF());
```

With output:

```
86.481 ms (316204 allocations: 29.61 MiB)
639.581 ms (6328167 allocations: 867.86 MiB)
┌ Warning: dt(1.4551915228366852e-11) <= dtmin(1.4551915228366852e-11) at t=5.708870444591281e-7. Aborting. There is either an error in your model specification or the true solution is unstable.
┌ Warning: dt(1.4551915228366852e-11) <= dtmin(1.4551915228366852e-11) at t=5.226666511312554e-7. Aborting. There is either an error in your model specification or the true solution is unstable.
```

I’m curious if there are any optimizations or tweaks to the Julia solvers I’m missing that could result in improved performance/stability vs CVODE_BDF? Or perhaps other Julia solvers are better suited to these systems (I have played around with ~20 alternatives)? Also, any insight into why QNDF is failing while CVODE_BDF is succeeding in the more complex simulation would be appreciated.

Thanks!