On Sat, Mar 1, 2008 at 12:43 PM, Lisandro Dalcin <[EMAIL PROTECTED]> wrote:
> 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? > > Could you attach the pure numpy solution along with a test case (alpha=?). Chuck
_______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
