Dear all, I want to comment some extrange stuff I'm experiencing with numpy. Please, let me know if this is expected and known.
I'm trying to solve a model nonlinear PDE, 2D Bratu problem (-Lapacian u - alpha * exp(u), homogeneus bondary conditions), using the simple finite differences with a 5-point stencil. I implemented the finite diference scheme in pure-numpy, and also in a F90 subroutine, next wrapped with f2py. Next, I use PETSc (through petsc4py) to solve the problem with a Newton method, a Krylov solver, and a matrix-free technique for the Jacobian (that is, the matrix is never explicitelly assembled, its action on a vector is approximated again with a 1st. order finite direrence formula). And the, surprise! The pure-numpy implementation accumulates many more inner linear iterations (about 25%) in the complete nonlinear solution loop than the one using the F90 code wrapped with f2py. Additionally, PETSc have in its source distribution a similar example, but implemented in C and using some PETSc utilities for managing structured grids. In short, this code is in C and completelly unrelated to the previously commented code. After running this example, I get almost the same results that the one for my petsc4py + F90 code. All this surprised me. It seems that for some reason numpy is accumulating some roundoff, and this is afecting the acuracy of the aproximated Jacobian, and then the linear solvers need more iteration to converge. Unfortunatelly, I cannot offer a self contained example, as this code depends on having PETSc and petsc4py. Of course, I could write myself the nonlinear loop, and a CG solver, but I am really busy. Can someone comment on this? Is all this expected? Have any of you experienced somethig similar? -- 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 _______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
