# Heuristic help with preconditioning large system ODEs

I'm looking for some general insight on preconditioning. In particular, relevant references/resources/comments would be greatly appreciated. Note, I have been through some of the literature, but am having trouble finding a reference for a problem similar to the one I am solving.

The system I'm solving takes the form

$Q_{ij} \ddot{y}_j + S_{ijk} \dot{y}_j\dot{y}_k +V_i =0$

where Einstein summation is assumed, and the arrays $Q_{ij}, S_{ijk},V_i$ are all low order polynomials in the dependent variables $\{ y_1, ... ,y_N\}$. Here, $N$ is the 'resolution' of the model, and $(i,j,k)\in(1,..,N)$.

Now, because $Q_{ij}$ is non trivial to invert, I've been solving this using an initial differential algebraic equation solver, namely the IDA solver that is a part of the SUNDIALS suite, via the FORTRAN interface. This solver uses a "variable-order, variable-coefficient" backwards differentiation formula. I also change variables, to make the system 1st order, i.e. $\dot{y}_k = v_k$.

The solver, in general, works as follows. Our governing equation can be written as

$F(t,\vec{y},\dot{\vec{y}}) = 0; \quad \vec{y}(t_o) = \vec{y}_o; \quad \dot{\vec{y}}=\dot{\vec{y}}_o$.

where $\vec{y} = (y_1,...,y_N)$. Now, the backwards differentiation formula, of qth order, implies

$\sum_{i=0}^q \alpha(n,i) \vec{y}(n-i) =h(n) \dot{\vec{y}}(n)$,

where $h(n)=t(n)-t(n-1)$ is the step size, and $\alpha(n,i)$ are functions of q and the history of the time steps. Substituting this back into the governing equation, we have equations of the form

$G(\vec{y}(n)) = F\left(t(n), \vec{y}(n),h(n)^{-1}\sum_{i=0}^q \alpha(n,i) \vec{y}(n-i)\right) = 0$.

The solver works by finding the solution to this equation via some form of Newton iteration. That is, it looks at

$J[\vec{y}(n)^{m+1}-\vec{y}(n)^m] = -G(\vec{y}(n)^m)$,

where the superscript $(m,m+1)$ refers to the $m$th, and $m+1$th approximation to the variable $\vec{y}$ and the Jacobian $J$ is given by

$J= \frac{\partial G}{\partial y} = \frac{\partial F}{\partial y}+ \alpha \frac{\partial F}{\partial \dot{y}}; \quad \alpha =\frac{\alpha(n,0)}{h(n)}$

Now, as the system approaches the phenomena I wish to study, the equations take much longer to integrate. I have optimized the 'low hanging fruit', e.g. the linear algebra is done using OpenBLAS, OpenMP is implemented etc...

I suspect that I need to use an iterative linear solver for these larger problems, and this calls for a preconditioner. My understanding is that a preconditioner $P$ is a matrix such $P^{-1}J$ gives the same solutions as would be the case for $P=I$, but is potentially easier to solve.

Now, I know that there exists a permanent progressive solution to my governing equation, such that $y_n = a_n e^{inct}$ where $a_n,c$ are numerically computed constants. The solutions I'm interested in modeling are going to be perturbations from these. So atleast for short times, the calculated solutions will be small deviations from the permanent progressive solutions.

So, I am interested in using information from these permanent progressive solutions to see if I can improve the efficiency of solving the more general equations.

In a sentence, I am unclear of how to do this, and have struggled to find relevant resources where specific examples are carried out in detail.

Naively, I would like to output the Jacobian of my permanent progressive solution so that I could get a better idea of its structure, and from here I could pick a preconditioner. I cannot seem to figure out how to do this in SUNDIALS.

Any information would be greatly appreciated,

Nick

Assuming you have the necessary Jacobian routine, I'd probably use a package that will compute an incomplete LU factorization as a "black-box" preconditioner to get started. In SUNDIALS, you will have to compute the factorization yourself; for that, I'd use something like SuperLU (you could also use MUMPS, or other routines). You can then supply the factored Jacobian matrix as a "user-supplied preconditioner". The IDA manual should have some information on how to do that in IDA.

Frankly, I find it easier to do these types of things in PETSc. It does have a better mechanism for run-time control and composition of

• time-steppers (including DAE solvers)
• nonlinear equation solvers
• nonlinear preconditioners
• Krylov subspace methods
• linear preconditioners (includes LU decomposition, multigrid, ILU, etc.)

and that run-time capability makes it much easier to experiment with "black-box"-type preconditioners.

PETSc has not yet implemented variable-step, variable-order BDF methods (for ODEs or DAEs), nor has it implemented sensitivity analysis (forward or adjoint), so if you need those features, you are still better off using SUNDIALS. You can call SUNDIALS from PETSc, but I believe that only CVODE is implemented, and it calls CVODE using the GMRES linear solver (I don't believe there is any flexibility in using different linear solvers when using SUNDIALS from PETSc).

If you know something about the physics of your system (for instance, removing coupling terms is a decent approximation of the physics), you might be able to construct a "physics-based" preconditioner that yields good performance. Good examples of physics-based preconditioners can be found in "Jacobian-Free Newton–Krylov methods: a survey of approaches and applications" by Dana A. Knoll and David E. Keyes, Journal of Computational Physics, Volume 193, Issue 2, pages 357-397.

• Thanks for the thoughtful reply, it's appreciated. I'll spend some time going through your suggestions, but I figure I'll throw out a naive question. 1) Do you know how one outputs the Jacobian, and $G$, in SUNDIALS? I suppose this a question for their mailing list. – Nick P Jun 11 '14 at 22:14
• Yes, although I'd have to look it up. You can use the difference quotient Jacobian routines (if you grep through the source or look through the manual, I believe they have DQJac in the name). It is more preferable to input the Jacobian information yourself, if you have that information; there is also information on that approach in the manual. I'm not sure if there is public access to $G$, and I don't think you'd need it. You should just be able to supply a preconditioner yourself through SUNDIALS interfaces, so you don't need to construct the Newton iterations yourself. – Geoff Oxberry Jun 11 '14 at 23:36
• I've found significant improvement with using a physically motivated subset of my jacobian as a preconditioner. However, my system is seriously limited by costly linear algebra routines (i.e. a large number of matrix-matrix operations), so I want to minimize the amount of calls to the residual function. I could input the entire jacobian, but it's even more costly than the residual function execution, as it depends on higher order arrays. Do you have any suggestions on how to minimize the # of calls? Or is this where picking a better preconditioner would come in handy? – Nick P Jul 3 '14 at 6:34
• If you want to minimize residual function calls, you probably want to avoid using difference quotients to calculate the Jacobian matrix (or any approximations). Good, cheap preconditioners will also help you. Linear system solves in the Newton step will involve two function evaluations for each Jacobian-vector product, so fewer linear system iterations = fewer residual evaluations. Failing that, you might try other integrators (such as in PETSc), or nonlinear preconditioners (also see PETSc). – Geoff Oxberry Jul 3 '14 at 7:15
• thanks for the tips. Can I ask you what a reasonable (I know these are the types of vague questions people hate, so I do apologize - but never the less I'll throw this out there) number of iterations I should expect between time steps for a generic non-stiff problem with a suitable time step? – Nick P Jul 6 '14 at 22:14