petsc4py python contexts
Ethan Coon
So I thought some might enjoy this nifty bit of hooks into PETSc's
matrix-free world, using petsc4py.
We all know PETSc is good for solving linear systems with an assortment
of solvers and preconditioners. We also know that python rocks for
high-level control code. petsc4py allows us to combine those things.
But what isn't so obvious is that petsc4py also enables us to use
PETSc's matrix-free concepts to solve systems of matrix-like equations,
without ever actually writing down full matrices.
Say, for instance, in an example that showed up in my work, you wanted
to apply boundary conditions on a Galerkin-FEM-discretized PDE via
Lagrange multipliers. You end up with a block system of matrix
equations which is also an indefinite saddle-point problem. Saddle
point problems are difficult to solve, because they're very poorly
preconditioned by traditional preconditioners. However, the discretized
interior PDE operator is well-preconditioned by traditional methods, and
the Lagrange Multiplier equations are well-preconditioned by the Schur
Complement.
Say our system is given by:
[A B^T ] u = v
[B 0 ] lmbda = mu
Then, by block elimination, we get the Schur complement, S =
-(B*(A^-1)*B^T). However, we'd never want to actually form S -- B, A,
and B^T are sparse, but A^-1 (and therefore S) are likely dense. So
we'd prefer to keep S as a "matrix-free" object... while we never
actually calculate S, we can determine the action of S on a vector by
using the various components.
The attached code does this through a bit of petsc4py hooks into PETSc
known as PythonContexts. Basically these tie into the PETSc matrix-free
world to enable using python code to apply the action of the matrix-free
Schur complement (and full, block system, which we also never fully form
as a matrix but keep in the block pieces using index sets).
I'll put this up on the groups page later for perusal... enjoy,
Ethan
--
-------------------------------------------
Ethan Coon
DOE CSGF - Graduate Student
Dept. Applied Physics & Applied Mathematics
Columbia University
212-854-0415
http://www.ldeo.columbia.edu/~ecoon/
-------------------------------------------
Lisandro Dalcin
> Hi all,
>
> But what isn't so obvious is that petsc4py also enables us to use
> PETSc's matrix-free concepts to solve systems of matrix-like equations,
> without ever actually writing down full matrices.
>
It is obvious, but just for the author :-) ... You know, there are so
much code there in mpi4py/petsc4py/slepc4py/tao4py, and so little
documentation...
>
> The attached code does this through a bit of petsc4py hooks into PETSc
> known as PythonContexts. Basically these tie into the PETSc matrix-free
> world to enable using python code to apply the action of the matrix-free
> Schur complement (and full, block system, which we also never fully form
> as a matrix but keep in the block pieces using index sets).
>
> I'll put this up on the groups page later for perusal... enjoy,
>
Mmm... perhaps something like this should also go to the demos/
directory for the next petsc4py release...
--
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
Ethan Coon
>
> Mmm... perhaps something like this should also go to the demos/
> directory for the next petsc4py release...
>
I think Braxton said there was something that used PythonContexts in
there (not sure where), but he couldn't understand much of it. You're
welcome to steal my code, which was just my way of experimenting with it
this summer. I can try to add some more comments if anyone has
questions/parts that are unclear.