Hello, I am trying to figure out a faster way of doing in petsc4py and thought I would ask to see how this should be done in petsc. If I can figure that out then I think I can figure out how to translate it but I am having issues with the basic algorithm I should be using.
I am solving a linear system of ODEs of the form B du/dt = A u using Crank-Nicholson. So the update equation can be written as (B - dt/2*A)*u^(n+1) = (B + dt/2*A)*u^n or, if you don't mind using the backslash operator, u^(n+1) = (B - dt/2*A)\(B + dt/2*A)*u^n The integration is faster of course after I build the matrix C = (B - dt/2*A)\(B + dt/2*A) however building this is very slow and I am not sure how to build it efficiently. At the moment the majority of my code is just doing that. In petsc4py I do the following: (you can find the full code at https://github.com/francispoulin/qg_1L_ts/blob/master/qg_1L_ts.py) AeT = PETSc.Mat().createAIJ([Ny-1,Ny-1]) AeT.setUp(); AeT.assemble() ksp = PETSc.KSP().create(PETSc.COMM_WORLD) ksp.setOperators(B-dto2*A) ksp.setTolerances(1e-16) pc = ksp.getPC() pc.setType('none') ksp.setFromOptions() btemp = B+dto2*A btemp.assemble() bcol = PETSc.Vec().createMPI(Ny-1) xcol = PETSc.Vec().createMPI(Ny-1) bcol.setUp(); xcol.setUp() sc,ec = bcol.getOwnershipRange() for i in range(0,Ny-1): bcol[sc:ec] = btemp[sc:ec,i] bcol.assemble(); xcol.assemble() ksp.solve(bcol,xcol) AeT[sc:ec,i] = xcol[sc:ec] AeT.assemble() We are building the matrix one row at a time and, as accurate as that maybe, it is very slow. Is there perhaps a fast way of doing this kind of operation in petsc? Sorry for the bother but any advice would be greatly appreciated. Cheers, Francis ------------------ Francis Poulin Associate Professor Associate Chair, Undergraduate Studies Department of Applied Mathematics University of Waterloo email: [email protected] Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637
