Solving Newton-Raphson step with ill-conditioned sparse matrix

Question

I am trying to build a complex simulator for the transport of mass & heat in porous media. I am currently following coarsely the algorithm laid out in an older software and have got the simulator to work, but the iterative solvers that I use run into trouble at large timesteps and prevent the convergence of the system to steady state. I was wondering if somebody could help me to improve my method.

The simulator builds the Jacobian matrix by numeric differentiation using double precision numbers. The desired final matrix size is in the order of 500,000x500,000 with each row having about 15-20 entries. The structure of the Jacobian is as follows:

• The diagonal entries stem from the thermodynamic properties of the elements in the system.
• The off-diagonal elements describe the transport of mass/heat between the elements. They scale with the time step. Due to the formulation currently used they are non-symmetric but only differ by a volume ratio between the elements; if need be they could be made symmetric.

The system J * (-dx) = r is solved using iterative solvers; in particular BICGSTAB with ILU0 works ok but I have tried others as well.

The problem arises when the timesteps get larger. The diagonal elements stay the same, but the off-diagonal elements scale with the timestep and the solvers used appear to run into problems once the timestep reaches ~1e+8 seconds.

What is the best way to proceed? I was thinking about one or more of the following:

• Use a solver with higher precision. But will this work, given that my derivative is done numerically to only ~8 significant digits?
• Is decomposing the matrix an option? If so, what decomposition would work?
• Are other solvers an option?

Any help is highly appreciated.

Cheers

Peter

Answer 1

For a general list of strategies, see:

• Why is Newton's method not converging?

• Why is my iterative linear solver not converging?

I ask a bunch of questions throughout the answer; please respond to them in your question as edits, because some of the answers to those questions are information that will help other users answer your question.

What is the best way to proceed? I was thinking about one or more of the following:

Use a solver with higher precision. But will this work, given that my derivative is done numerically to only ~8 significant digits?

I believe the main consequence of using higher precision in the solver will be that you preserve more significant digits within the solver. If the finite differencing is done outside the solver, this will (obviously) be unaffected; if the derivative is calculated in higher precision, it gives you more options in terms of selecting the numerical differencing parameter. Have you tried varying that differencing parameter and examining its effect on the conditioning of your system?

Is using an analytical Jacobian matrix (or a reasonably good analytical approximation thereof) an option? Something that can be problematic is when the numerical differencing parameter is of the same order as relevant numerical parameters (or scales) in your problem. This issue was one I encountered in debugging a reaction-diffusion simulation: a matrix-free approximation of the Jacobian was not sufficiently accurate and caused ill-conditioning issues. Replacing the matrix-free approximation of the Jacobian matrix solved my issues. (See also: When is Newton-Krylov not an appropriate solver?)

Is decomposing the matrix an option? If so, what decomposition would work?

No, it's not. Using something like an sparse LU decomposition will likely require too much memory (you already have ~10 million nonzeroes, and fill-in will only increase that number).

Are other solvers an option?

Probably. What have you tried? Something that would also help for analysis would be an estimate of the conditioning of your system. GMRES (without restarts) provides one such estimate, and it would be helpful to know if the linear systems you're solving are ill-conditioned.

Specifically, it would be helpful to know what:

• time-steppers
• preconditioners
• linear solvers
• libraries

you are using. (Mentioning ILU(0) and BiCGStab is a good start.)

It would also be helpful to see equations, and understand more about your problem. Is it parabolic? Are there likely to be the same sorts of errors at multiple scales? Multigrid methods might be helpful for preconditioning, and scale better than ILU(0).

1e8 is a really large time step, without any other context. Why such a large number? Is there a scale issue? Can the problem be rescaled or nondimensionalized?

Answer 2

Is decomposing the matrix an option? If so, what decomposition would work?

Yes, I think that is a very reasonable option. 500000 is not a huge number of equations and your system is quite sparse. For your current, unsymmetric case you would need to perform an $LU$ decomposition. If you can symmetrize it, as you suggest, you can save almost 1/2 the memory by performing an $LL^T$ or $LDL^T$ decomposition.

I suggest trying the MUMPS sparse solver:

http://mumps.enseeiht.fr/

It has options for all three types of factorizations. It is designed to run in parallel on distributed memory computers but also runs quite efficiently on single-CPU or multicore machines. It has an out-of-core option that you can use to reduce memory usage.