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.