question on DA's and performance

[Randall Mackie]

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