Dear Charles, As I said, I have no time to code the pure Python+numpy nonlinear and linear loops, and the matrix-free stuff to mimic the PETSc implementation. However, I post the F90 code and the numpy code, and a small script for testing with random input. When I have some spare time, I'll try to do the complete application in pure python.
Regards, On 3/1/08, Charles R Harris <[EMAIL PROTECTED]> wrote: > > > > 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 > > -- 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
bratu2dlib.f90
Description: Binary data
import numpy
def bratu2d(alpha, x, f):
from numpy import exp
nx, ny = x.shape
hx = 1.0/(nx-1)
hy = 1.0/(ny-1)
#
u = x[1:-1, 1:-1] # center
uN = x[1:-1, :-2] # north
uS = x[1:-1, 2: ] # south
uE = x[ :-2, 1:-1] # east
uW = x[2:, 1:-1] # west
#
f[:,:] = x
f[1:-1, 1:-1] = \
(2*u - uE - uW) * (hy/hx) \
+ (2*u - uN - uS) * (hx/hy) \
- alpha * exp(u) * (hx*hy)
import numpy as N
import bratu2dlib # f90 version
import bratu2dnpy # numpy version
for i in range(10):
m, n = 32, 32
alpha = 6.8
X = N.random.rand(m*n).reshape(m,n,order='f')
F1 = N.zeros([m,n], order='f')
F2 = N.zeros([m,n], order='f')
bratu2dlib.bratu2d(alpha, X, F1)
bratu2dnpy.bratu2d(alpha, X, F2)
D = F1 - F2
print D.max()
_______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
