1*55a74a43SLisandro Dalcintry: 2*55a74a43SLisandro Dalcin execfile 3*55a74a43SLisandro Dalcinexcept NameError: 4*55a74a43SLisandro Dalcin def execfile(file, globals=globals(), locals=locals()): 5*55a74a43SLisandro Dalcin fh = open(file, "r") 6*55a74a43SLisandro Dalcin try: exec(fh.read()+"\n", globals, locals) 7*55a74a43SLisandro Dalcin finally: fh.close() 8*55a74a43SLisandro Dalcin 9*55a74a43SLisandro Dalcinimport petsc4py, sys 10*55a74a43SLisandro Dalcinpetsc4py.init(sys.argv) 11*55a74a43SLisandro Dalcin 12*55a74a43SLisandro Dalcinfrom petsc4py import PETSc 13*55a74a43SLisandro Dalcin 14*55a74a43SLisandro Dalcinexecfile('petsc-mat.py') 15*55a74a43SLisandro Dalcinexecfile('petsc-cg.py') 16*55a74a43SLisandro Dalcin 17*55a74a43SLisandro Dalcinx, b = A.createVecs() 18*55a74a43SLisandro Dalcin 19*55a74a43SLisandro Dalcinksp = PETSc.KSP().create() 20*55a74a43SLisandro Dalcinksp.setType('cg') 21*55a74a43SLisandro Dalcinksp.getPC().setType('none') 22*55a74a43SLisandro Dalcinksp.setOperators(A) 23*55a74a43SLisandro Dalcinksp.setFromOptions() 24*55a74a43SLisandro Dalcin 25*55a74a43SLisandro Dalcinksp.max_it = 100 26*55a74a43SLisandro Dalcinksp.rtol = 1e-5 27*55a74a43SLisandro Dalcinksp.atol = 0 28*55a74a43SLisandro Dalcinx.set(0) 29*55a74a43SLisandro Dalcinb.set(1) 30*55a74a43SLisandro Dalcinksp.solve(b,x) 31*55a74a43SLisandro Dalcinprint("iterations: %d residual norm: %g" % (ksp.its, ksp.norm)) 32*55a74a43SLisandro Dalcin 33*55a74a43SLisandro Dalcinx.set(0) 34*55a74a43SLisandro Dalcinb.set(1) 35*55a74a43SLisandro Dalcinits, norm = cg(A,b,x,100,1e-5) 36*55a74a43SLisandro Dalcinprint("iterations: %d residual norm: %g" % (its, norm)) 37*55a74a43SLisandro Dalcin 38*55a74a43SLisandro DalcinOptDB = PETSc.Options() 39*55a74a43SLisandro Dalcin 40*55a74a43SLisandro Dalcinif OptDB.getBool('plot', True): 41*55a74a43SLisandro Dalcin da = PETSc.DMDA().create([m,n]) 42*55a74a43SLisandro Dalcin u = da.createGlobalVec() 43*55a74a43SLisandro Dalcin x.copy(u) 44*55a74a43SLisandro Dalcin draw = PETSc.Viewer.DRAW() 45*55a74a43SLisandro Dalcin OptDB['draw_pause'] = 1 46*55a74a43SLisandro Dalcin draw(u) 47*55a74a43SLisandro Dalcin 48*55a74a43SLisandro Dalcinif OptDB.getBool('plot_mpl', False): 49*55a74a43SLisandro Dalcin try: 50*55a74a43SLisandro Dalcin from matplotlib import pylab 51*55a74a43SLisandro Dalcin except ImportError: 52*55a74a43SLisandro Dalcin print("matplotlib not available") 53*55a74a43SLisandro Dalcin else: 54*55a74a43SLisandro Dalcin from numpy import mgrid 55*55a74a43SLisandro Dalcin X, Y = mgrid[0:1:1j*m,0:1:1j*n] 56*55a74a43SLisandro Dalcin Z = x[...].reshape(m,n) 57*55a74a43SLisandro Dalcin pylab.figure() 58*55a74a43SLisandro Dalcin pylab.contourf(X,Y,Z) 59*55a74a43SLisandro Dalcin pylab.plot(X.ravel(),Y.ravel(),'.k') 60*55a74a43SLisandro Dalcin pylab.axis('equal') 61*55a74a43SLisandro Dalcin pylab.colorbar() 62*55a74a43SLisandro Dalcin pylab.show() 63