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.
