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
