Stabilization of the SIESTA MHD Equilibrium Code Using Rapid Cholesky Factorization

The SIESTA MHD equilibrium code solves the discretized nonlinear MHD force $F\equiv J$X$B-\nabla p$ for a 3D plasma which may contain islands and stochastic regions. At each nonlinear evolution step, it solves a set of linearized MHD equations which can be written $r\equiv $\textit{Ax}$-b=$0, where $A$ is the linearized MHD Hessian matrix. When the solution norm \textbar $x$\textbar is small enough, the nonlinear force norm will be close to the linearized force norm \textbar $r$\textbar $\approx $ 0 obtained using preconditioned GMRES. In many cases, this procedure works well and leads to a vanishing nonlinear residual (equilibrium) after several iterations in SIESTA. In some cases, however, \textbar $x$\textbar \textgreater 1 results and the SIESTA code has to be restarted to obtain nonlinear convergence. In order to make SIESTA more robust and avoid such restarts, we have implemented a new rapid QR factorization of the Hessian which allows us to rapidly and accurately solve the least-squares problem $ A^{T}r =$ 0, subject to the condition \textbar $x$\textbar \textless 1. This avoids large contributions to the nonlinear force terms and in general makes the convergence sequence of SIESTA much more stable. The innovative rapid QR method is based on a pairwise row factorization of the tri-diagonal Hessian. It provides a complete Cholesky factorization while preserving the memory allocation of $A$.