adding SNESSetLinearSolve()

[Lisandro Dalcin]

On 11/1/07, Barry Smith <bsmith at mcs.anl.gov> wrote:
> > On Nov 1, 2007 9:24 AM, Lisandro Dalcin <dalcinl at gmail.com> > >> For example, I would love to have SNESSetPresolve/SNESSetPostSolve
> >> and
> >> SNESSetPreStep/SNESSetPostStep, and perhaps a
> >> SNESSetPreLinearSolve/SNESSetPostLinearSolve. Of course, this make
> >> the
> >> API grow with features that are not frecuently needed.
>
>    I'd like to understand why they are ever needed and if they can be
> done
> with current constructs?
>

Barry, sorry for the delay, but I was really busy. I'll try to clarify my needs.

* I'm using PETSc exclusively from Python. In such a dynamic language,
implementing complicated application setups is far easier than in
C/C++.

* Of course, Python is no efficient for low-level stuff. Do I actually
make use of a PETSc-based C++ application (PETSc-FEM) to assemble
vectors and matrices, manage fixations/boundary conditions and linear
constraints between dofs.

* Some application setups will need PETSc solvers (SNES/TS) to provide
a way to set user callback routines (or hooks). That's the reason I
talked about adding SNESSetPresolve/SNESSetPostSolve and
SNESSetPreStep/SNESSetPostStep. This way, PETSc can be adapted to the
application needs if required, and not the other way.

* Taking again SNES as an example: Why it provides SNESSetUpdate()???
It is not used in any other place in PETSc sources. I guess at some
point someone needed a user callback routine to be called just before
computing the Jacobian, so it was added. I'm currently using it to
compute/update auxiliar linear operators to be used in defining a
problem-based preconditioner for NS equations. Of course, I could do
this in computeJacobian() rutine, but then I would have to special
case the use of this preconditioner in order to compute the auxiliar
operator; in short my application has to adapt to PETSc.

* Taking TS as an exaple. Why it provides SetPreStep/SetPostStep???
They are also not even used in a TS example. I guess again that at
some point someone needed it. And then, I'm asking to have more of
that, and no less.

In the end, perhaps the better way for me to go is to implement SNES
and TS directly in Python and do what I want, it is not much hard and
there should not be performace issues at this level. But I would
really prefer to make SNES/TS more configurable (via setting general
pourpose user callback routines) in order to not duplicate code and
efforts.

So the real quention you have to answer is the following: Should I
continue fighting for making PETSc solver more 'user-configurable', or
should I direct my efforts in implemening all those features directly
in Python, where implementing truly polimorfic code is far, far easier
than in C?

In the near future, application written in C/C++/Fortran could take
advantage of petsc4py to implement custom 'shell' matrices, linear
solvers, preconditioners, nonlinear solvers and time steppers.
Something like
SNESCreate(comm,&snes)
SNESSetType(snes, "python")
SNESPythonCreate(snes, "package.module.MyNewtonSolver")

Woud'nt that be a really nice feature?

I'm not sure if all you consider Python (plus wrapped compiled code)
as a feasible way to drive large-scale simulations, but I really
believe it is. Furthermore, from the last three years I worked hard to
make it possible: now Python users can use mpi4py, petsc4py, slepc4py
and tao4py and manage applications with negligible overhead for medium
to large data transfers/ problems.


-- 
Lisandro Dalcín
---------------
Centro Internacional de Métodos Computacionales en Ingeniería (CIMEC)
Instituto de Desarrollo Tecnológico para la Industria Química (INTEC)
Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET)
PTLC - Güemes 3450, (3000) Santa Fe, Argentina
Tel/Fax: +54-(0)342-451.1594