Solving nonlinear PDE with finite difference based on Newton-Krylov

I am now working on solving MHD equations with finite difference method, which include nonlinear equations:
$$
frac{partialrho}{partial t}+nablacdotleft[left(rho_0+rhoright){v}right]-nablacdotleft(dnablarhoright)=0
frac{partial{v}}{partial t}+nabla{{v}}cdot{{v}}-frac{1}{rho+rho_0}left[{j}_0times{B}+left(nablatimes{B}right)timesleft({B}+{B}_0right)right]-nablacdotleft(nunabla{v}right)=0
frac{partial{B}}{partial t}-nablatimesleft({v}times{B}right)-etanabla^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.