dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Tue, Mar 14, 2006 at 07:21:26PM +0100, Garth N. Wells wrote:
> On Tue, 2006-03-14 at 11:39 -0600, Anders Logg wrote:
> > On Tue, Mar 14, 2006 at 06:32:27PM +0100, Garth N. Wells wrote:
> > 
> > > > Cleanup KrylovSolver:
> > > > 
> > > >  - Remove setFoo() functions in public interface
> > > >  - Cleanup and restructure implementation
> > > > 
> > > > Please check that everything looks ok. I think it the code is easier
> > > > to follow now, but I may be wrong.
> > > > 
> > > 
> > > Looks a lot nicer, and now solver parameters can be set through objects
> > > which use use KrylovSolver (like PDE and NewtonSolver). To allow full
> > > functionality, I think that LU should become a special case of a
> > > KrylovSolver so that an LU solver can be chosen when using  PDE,
> > > NewtonSolver or any other class which uses a solver. 
> > 
> > I think that's already possible:
> > 
> >     pde.set("solver", "direct");
> > 
> > That will choose the LU solver, at least in LinearPDE.
> > 
> 
> I was thinking more about NewtonSolver, for which the solver is a
> private member object. I shpuld change this to be a pointer to the
> solver, and create the appropriate solver type in the NewtonSolver
> constructor.

ok.

> > The LU solver is a special case of Krylov in PETSc. We could make it a
> > special case of the KrylovSolver, but there's enough code in LU.cpp
> > that I'm not so keen on throwing into KrylovSolver.cpp.
> > 
> 
> I had a quick look at the code in LU.cpp. How much of it is being used
> at the moment? It looks like it could be due for a clean-up. (That said,
> some dense matrix functionality would be useful which I see floating
> around in there.)

I haven't looked at the code for a while, but it looks pretty ok to
me. There is some debugging that could be removed, but feel free to
clean it up if you find something that can be improved.

Most of the stuff in there, in particular copyToDense() is for solving
a system with a VirtualMatrix (only defined by the product) with a
direct LU factorization, which is a little special. The copyToDense()
function extracts a dense representation by applying the matrix-vector
product to unit vectors. This is only used by the ODE solvers as far
as I know. It's also good to have for testing.

> > > I see that the private member functions setSolver() and
> > > setPreconditioner() are still there. They could be removed and the
> > > necessary code placed in KrylovSolver::init(..). Alternatively, they
> > > could be kept and called from within KrylovSolver::init(..) rather than
> > > KrylovSolver::solve(..).
> > > 
> > > Garth
> > 
> > I thought about doing this (calling them from the init() function) but
> > I didn't want to mess with the logic of the init() function (return
> > when not needing to initialize). We could move the stuff in init() to
> > say initKSP() and then call all three from init()?
> > 
> 
> If it's only possible to set the solver type in the constructor, it
> would be cleaner to set the solver type in init(), or in a function
> called from init(), rather than setting the solver every time solve() is
> called. At first glance, this gives the impression that the solver type
> can be changed. 

I agree, but I think it must be done each time since if the system
changes size between solves, then the KSP pointer must be recreated
and thus the solver and preconditioner reset.

When I think of it, this would make it ideal to just throw in the
setSolver() and setPreconditioner() in the init() function. I'll take
a look...

> The other option is to enter the solver type into the parameter system,
> and let setSolver() pull it out of the database and set it every time
> solve() is called. 

Perhaps. It's difficult to say what should go in the parameter system
and what should be an option. The main principle could be that
essential stuff should be arguments and non-essential stuff should go
in the parameter system. My feeling is that the type of solver and
preconditioner is essential, people would expect to be able to specify
it to the constructor, whereas the number of iterations etc is some
parameter.

In the context of solving a PDE, the type of solver is a parameter
("direct" or "iterative") since one could expect the solver type to be
chosen automatically. On the other hand, when working with the class
KrylovSolver, the type of iterative method and preconditioner are
essential and should perhaps not be parameters.

/Anders