[Numpy-discussion] need help with simple conjugate gradient Laplace solver
rob
rob at pythonemproject.com
Sat Mar 30 09:00:03 EST 2002
After I post, I always see the dumb error. I am not including the 6x
term in my finite difference equation. It now converges, but I get
wierd looking V map. Rob.
Here is the fixed 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
N=30 # size of box
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:N,0:N,0]=0 # box at 1V with one side at 0V
u[0:N,0,0:N]=1
u[0,0:N,0:N]=1
u[0:N,0:N,N]=1
u[0:N,N,0:N]=1
u[N,0:N,0:N]=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]-
6*u[1:NX-1,1:NY-1,1:NZ-1])
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=0.0; #set residuals=0
while(1):
q[1:NX-1,1:NY-1,1:NZ-1]=(6*p[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
More information about the NumPy-Discussion
mailing list