# 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{\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.