I am now working on solving MHD equations with finite difference method, which include nonlinear equations: $$ \frac{\partial\rho}{\partial t}+\nabla\cdot\left[\left(\rho_0+\rho\right){v}\right]-\nabla\cdot\left(d\nabla\rho\right)=0\\ \frac{\partial{v}}{\partial t}+\nabla{{v}}\cdot{{v}}-\frac{1}{\rho+\rho_0}\left[{j}_0\times{B}+\left(\nabla\times{B}\right)\times\left({B}+{B}_0\right)\right]-\nabla\cdot\left(\nu\nabla{v}\right)=0\\ \frac{\partial{B}}{\partial t}-\nabla\times\left({v}\times{B}\right)-\eta\nabla^2{{B}}=0 $$ where $\rho$,$v=(v_x,v_y,v_z)$,$B=(B_x,B_y,B_z)$ are the variables need to be solved, $\rho_0$,$d$,$\nu$,$\eta$ are constant scalar fields and $B_0$ is constant vector field. While the equations are definitely nonlinear, I suppose to solve them with Newton method. Finite difference method is used to discretize the equation(1-order on time and 2-order central difference on space). The Jacobian matrix is calculated as follow: $$ DF=[\frac{\partial F_{i,j,k}}{\partial x}], x=(\rho^{n+1}_{0,0,0},\rho^{n+1}_{0,0,1},\cdots,\rho^{n+1}_{0,1,0},\cdots,\rho^{n+1}_{1,0,0},\cdots,{B_x}_{0,0,0}^{n+1},\cdots) $$ where $F$ is the left-hand side of the equation, and I use implicit scheme.

The computation model built for the problem is quit large (num. of nodes > 2,000,000), to solve the huge linear problem in acceptable time, I try to solve it with PetSc library on a parallel platform, GMRES method in KSP is selected as the linear solver.

However, the essential computation time consumption is the linear solving progress. I suspect that the Jacobian matrix is ill-conditioned and caused terrible efficiency.

The structure of the matrix is "multi-diagonal", so is there a way to reduce the solving time? I once try some simple preconditioners like Jacobi which are included by PetSc but get no speed-up. ILU is not support in parallel PetSc program.

Sincerely thanks for your reading.

  • 4
    $\begingroup$ Have you considered to do the linearization for the newton iteration before discretization? This can help spot ill-conditioning or even prevent it, e.g. see my answer here $\endgroup$
    – Bort
    Commented Jun 6, 2019 at 11:15
  • $\begingroup$ You can also try using a Fast Poisson Solver as a preconditioner, but in 3D, this will still be slow $\endgroup$
    – whpowell96
    Commented Mar 3, 2020 at 2:40

1 Answer 1


I'm not familiar with the term "multi diagonal", but presumably you mean either that your matrix is banded, or that it can be reordered into a matrix that is banded. In that case you might try using the SPIKE package:


It's specifically designed for solving linear systems of equations with banded matrices. It can be used to solve banded matrix problems directly, or as a preconditioner.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.