I've been using Petsc for a few years with reasonably good success. My application is 3D EM forward modeling and inversion.
What has been working well is basically an adaptation of what I did in serial mode, by solving the following system of equations: |Mxx Mxy Mxz| |Hx| |bx| |Myx Myy Myz| |Hy| = |by| |Mzx Mzy Mzz| |Hz| |bz| Because this system is very stiff and ill-conditioned, the preconditioner that has been successfully used is the ILU(k) of the diagonal sub-blocks of the coefficient matrix only, not the ILU(k) of the entire matrix. I've tried, with less success, converting to distributed arrays, because in reality my program requires alternating between solving the system above, and solving another system based on the divergences of the magnetic field, and so this requires a lot of message passing between the nodes. I think that using DA's would require passing only the ghost values instead of all the values as I am now doing. So I coded up DA's, and it works, but only okay, and not as well as the case above. The main problem is that it seems to work best for ILU(0), whereas the case above works better and better with the more fill-in specified. I think I've tried every possible option, but I can't get any improvement over just using ILU(0), but the problem is that with ILU(0), it takes too many iterations, and so the total time is more than when I use my original method with, say, ILU(8). The differences between using DA's and the approach above is that the fields are interlaced in DA's, and the boundary values are included as unknowns (with coefficient matrix values set to 1.0), whereas my original way the boundary values are incorporated in the right-hand side. Maybe that hurts the condition of my system. I would really like to use DA's to improve on the efficiency of the code, but I can't figure out how to do that. It was suggested to me last year to use PCFIELDSPLIT, but there are no examples, and I'm not a c programmer so it's hard for me to look at the source code and know what to do. (I'm only just now able to get back to this). Does anyone do any prescaling of their systems? If so, does anyone have examples of how this can be done? Any advice is greatly appreciated. Thanks, Randy -- Randall Mackie GSY-USA, Inc. PMB# 643 2261 Market St., San Francisco, CA 94114-1600 Tel (415) 469-8649 Fax (415) 469-5044 California Registered Geophysicist License No. GP 1034
