1*55a74a43SLisandro Dalcinimport sys, petsc4py 2*55a74a43SLisandro Dalcinpetsc4py.init(sys.argv) 3*55a74a43SLisandro Dalcin 4*55a74a43SLisandro Dalcinfrom petsc4py import PETSc 5*55a74a43SLisandro Dalcinimport Bratu3D as Bratu3D 6*55a74a43SLisandro Dalcin 7*55a74a43SLisandro DalcinOptDB = PETSc.Options() 8*55a74a43SLisandro Dalcin 9*55a74a43SLisandro DalcinN = OptDB.getInt('N', 16) 10*55a74a43SLisandro Dalcinlambda_ = OptDB.getReal('lambda', 6.0) 11*55a74a43SLisandro Dalcindo_plot = OptDB.getBool('plot', False) 12*55a74a43SLisandro Dalcin 13*55a74a43SLisandro Dalcinda = PETSc.DMDA().create([N, N, N], stencil_width=1) 14*55a74a43SLisandro Dalcin#app = App(da, lambda_) 15*55a74a43SLisandro Dalcin 16*55a74a43SLisandro Dalcinsnes = PETSc.SNES().create() 17*55a74a43SLisandro DalcinF = da.createGlobalVec() 18*55a74a43SLisandro Dalcinsnes.setFunction(Bratu3D.formFunction, F, 19*55a74a43SLisandro Dalcin args=(da, lambda_)) 20*55a74a43SLisandro DalcinJ = da.createMat() 21*55a74a43SLisandro Dalcinsnes.setJacobian(Bratu3D.formJacobian, J, 22*55a74a43SLisandro Dalcin args=(da, lambda_)) 23*55a74a43SLisandro Dalcin 24*55a74a43SLisandro Dalcinsnes.setFromOptions() 25*55a74a43SLisandro Dalcin 26*55a74a43SLisandro DalcinX = da.createGlobalVec() 27*55a74a43SLisandro DalcinBratu3D.formInitGuess(X, da, lambda_) 28*55a74a43SLisandro Dalcinsnes.solve(None, X) 29*55a74a43SLisandro Dalcin 30*55a74a43SLisandro DalcinU = da.createNaturalVec() 31*55a74a43SLisandro Dalcinda.globalToNatural(X, U) 32*55a74a43SLisandro Dalcin 33*55a74a43SLisandro Dalcin 34*55a74a43SLisandro Dalcindef plot(da, U): 35*55a74a43SLisandro Dalcin scatter, U0 = PETSc.Scatter.toZero(U) 36*55a74a43SLisandro Dalcin scatter.scatter(U, U0, False, PETSc.Scatter.Mode.FORWARD) 37*55a74a43SLisandro Dalcin rank = PETSc.COMM_WORLD.getRank() 38*55a74a43SLisandro Dalcin if rank == 0: 39*55a74a43SLisandro Dalcin solution = U0[...].reshape(da.sizes, order='f').copy() 40*55a74a43SLisandro Dalcin try: 41*55a74a43SLisandro Dalcin from matplotlib import pyplot 42*55a74a43SLisandro Dalcin pyplot.contourf(solution[:, :, N//2]) 43*55a74a43SLisandro Dalcin pyplot.axis('equal') 44*55a74a43SLisandro Dalcin pyplot.show() 45*55a74a43SLisandro Dalcin except: 46*55a74a43SLisandro Dalcin pass 47*55a74a43SLisandro Dalcin PETSc.COMM_WORLD.barrier() 48*55a74a43SLisandro Dalcin scatter.destroy() 49*55a74a43SLisandro Dalcin U0.destroy() 50*55a74a43SLisandro Dalcin 51*55a74a43SLisandro Dalcin 52*55a74a43SLisandro Dalcinif do_plot: plot(da, U) 53*55a74a43SLisandro Dalcin 54*55a74a43SLisandro Dalcin 55*55a74a43SLisandro DalcinU.destroy() 56*55a74a43SLisandro DalcinX.destroy() 57*55a74a43SLisandro DalcinF.destroy() 58*55a74a43SLisandro DalcinJ.destroy() 59*55a74a43SLisandro Dalcinda.destroy() 60*55a74a43SLisandro Dalcinsnes.destroy() 61