dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Thu, Mar 16, 2006 at 01:48:13PM +0100, Garth N. Wells wrote:
> On Thu, 2006-03-16 at 06:22 -0600, Anders Logg wrote: 
> > On Thu, Mar 16, 2006 at 06:09:23AM -0600, Anders Logg wrote:
> > 
> > > > This appears to be a PETSc problem when destroying a solver before it
> > > > has been used to solve a linear system.
> > > > 
> > > > The problem can be avoided by not calling init(0,0) (as well as
> > > > readParameters()) in the KrylovSolver constructor. As far as I can see,
> > > > there is no need to call init(0,0) in the constructor. It creates a
> > > > PETSc solver, which is then deleted by KrylovSolver::init() the first
> > > > time KrylovSolver::solve() is called. Can the call be removed from the
> > > > constructor?
> > > > 
> > > > Garth
> > > 
> > > Aha. We had the following in the init() function before. I couldn't
> > > understand why it was there in the first place so I removed it. Seems
> > > we need it.
> > > 
> > >   // Don't reinitialize on first solve
> > >   if ( ksp != 0 && this->M == 0 && this->N == 0 )
> > >   {
> > >     this->M = M;
> > >     this->N = N;
> > >     return;
> > >   }
> > > 
> > > Does it work if you put it back?
> > > 
> > > I think we need to do init() in the constructor so that it is possible
> > > to call KrylovSolver::solver() and do stuff with the KSP pointer if
> > > someone should want to do that.
> > > 
> > > /Anders
> > 
> > When I think of it, maybe we can remove the KrylovSolver::solver()
> > function? The DOLFIN KrylovSolver is anyway doing pretty much what it
> > wants to do by itself and would perhaps not obey all options specified
> > by a user, especially since ksp is sometimes destroyed and then
> > created again in a solve.
> > 
> 
> Yes, this is an issue. Also, if a user asks for the pointer to the PETSc
> solver and makes some changes (different solver or preconditioner for
> example), this may be rolled-back by KrylovSolver::init() which applies
> a different solver or preconditioner. At the moment, to avoid this
> problem the user needs to create a KrylovSolver with the
> KrylovSolver::default_solver and KrylovSolver::default_pc
> 
> On the flip side, the PETSc solvers have a lot of options which are
> needed for particular problems, so I think that we still need to access
> this somehow. Could we provide two solver interfaces - one standard with
> the regular features, and an advanced interface where the user is left
> more on their own to do things correctly?
> 
> > If we remove it, we don't need to do init() in the constructor and the
> > problem will go away. Would that work?
> > 
> > We should probably add some comments that (1) PETSc gets confused if
> > the system changes size between solves (which is why we need to do
> > init() between solves) and that (2) PETSc gets confused if KSP is
> > created, destroyed and created again before a solve (which we avoid).
> 
> On the last point, I've only seen this problem in combination with the
> Hypre preconditioners. I'm hoping that this will be fixed in PETSc. 
> 
> Garth

My suggestion is that we remove KrylovSolver::solver(). If someone
wants to do something very special with PETSc, they might as well do
it through PETSc anyway. For Matrix and Vector, we have the mat() and
vec() functions which make it possible to access the PETSc data for
use with PETSc functions. I think this is all we need to do. If
someone wants to do something very special, just ask for the PETSc
data and call PETSc.

It would also make life much easier for us.

Of course, we should make sure to add the most common options as
parameters. We don't need to cover the whole range of PETSc options,
but it's very easy to just throw in a new parameter in
DefaultParameters.h and just read the value in
KrylovSolver::readParameters().

Any opinions?

/Anders