[Numpy-discussion] need help with simple conjugate gradient Laplace solver

[rob]

I'm experimenting with electrostatics now.  I have an iterative Jacobian
Laplace solver working but it can be slow.  It creates a beautiful 3D
Animabob image.  
So I decided to try out a conjugate-gradient solver, which should be an
order of mag better.  It runs but doesn't converge.  One thing I am
wondering, where is the conjugate?  In my FEM code, the solver realy
does use a conjugate, while this one here that I pieced together from
several other programs does not.  Why is it called conjugate gradient
without a conjugate ?  :)   Here is the code:


from math import *
from Numeric import *

#
#***  ENTER DATA
filename= "out"
#
bobfile=open(filename+".bob","w")

print "\n"*30

NX=30    # number of cells

NY=30

NZ=30

resmax=1e-3    # conj-grad tolerance

#allocate arrays

##Ex=zeros((NX+2,NY+2,NZ+2),Float)
##Ey=zeros((NX+2,NY+2,NZ+2),Float)
##Ez=zeros((NX+2,NY+2,NZ+2),Float)

p=zeros((NX+1,NY+1,NZ+1),Float)
q=zeros((NX+1,NY+1,NZ+1),Float)
r=zeros((NX+1,NY+1,NZ+1),Float)
u=zeros((NX+1,NY+1,NZ+1),Float)

u[0:30,0:30,0]=0    # box at 1V with one side at 0V
u[0:30,0,0:30]=1
u[0,0:30,0:30]=1
u[0:30,0:30,30]=1
u[0:30,30,0:30]=1
u[30,0:30,0:30]=1

r[1:NX-1,1:NY-1,1:NZ-1]=(u[1:NX-1,0:NY-2,1:NZ-1]+            #initialize
r matrix
                                          u[1:NX-1,2:NY,1:NZ-1]+
                                          u[0:NX-2,1:NY-1,1:NZ-1]+
                                          u[2:NX,1:NY-1,1:NZ-1]+
                                          u[1:NX-1,1:NY-1,0:NZ-2]+
                                          u[1:NX-1,1:NY-1,2:NZ])

p[...]=r[...]   #initialize p matrix



#
#**** START ITERATIONS
#
N=(NX-2)*(NY-2)*(NZ-2)   # left over from Jacobi solution, not used
KK=0     # iteration counter

res=resmax=0.0;    #set  residuals=0


while(1):     

    q[1:NX-1,1:NY-1,1:NZ-1]=(p[1:NX-1,0:NY-2,1:NZ-1]+       # finite
difference eq
                                              p[1:NX-1,2:NY,1:NZ-1]+
                                              p[0:NX-2,1:NY-1,1:NZ-1]+
                                              p[2:NX,1:NY-1,1:NZ-1]+
                                              p[1:NX-1,1:NY-1,0:NZ-2]+
                                              p[1:NX-1,1:NY-1,2:NZ])
    
    
#    Calculate r dot p and p dot q  

    rdotp = 0.0
    pdotq = 0.0
    
    rdotp = add.reduce(ravel( r[1:NX-1,1:NY-1,1:NZ-1] *
p[1:NX-1,1:NY-1,1:NZ-1]))
    pdotq = add.reduce(ravel( p[1:NX-1,1:NY-1,1:NZ-1] *
q[1:NX-1,1:NY-1,1:NZ-1]))

#    Set alpha value  

    alpha = rdotp/pdotq

#   Update solution and residual  

    u[1:NX-1,1:NY-1,1:NZ-1] += alpha*p[1:NX-1,1:NY-1,1:NZ-1]
    r[1:NX-1,1:NY-1,1:NZ-1] +=  - alpha*q[1:NX-1,1:NY-1,1:NZ-1]
    
#    calculate beta  

    rdotq = 0.0
    
    rdotq = 
add.reduce(ravel(r[1:NX-1,1:NY-1,1:NZ-1]*q[1:NX-1,1:NY-1,1:NZ-1]))
    
    beta = rdotq/pdotq

#   Set the new search direction  

    p[1:NX-1,1:NY-1,1:NZ-1] = r[1:NX-1,1:NY-1,1:NZ-1] -
beta*p[1:NX-1,1:NY-1,1:NZ-1]
    
   
    res = sort(ravel(r[1:NX-1,1:NY-1,1:NZ-1]))[-1]    #find largest
residual
#    resmax = max(resmax,abs(res))


    KK+=1
    #      
    print "Iteration Number %d  Residual %1.2e" %(KK,abs(res))
    if (abs(res)<=resmax): break     # if residual is small enough break
out


print "Number of Iterations ",KK

-- 
-----------------------------
The Numeric Python EM Project

www.pythonemproject.com