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