Hello Matt, Thanks for the help.
I had some some tests in Matlab using the two methods and because I am trying to compute the growth rates of a system that are very small, it seemed to be faster to find the inverse but I will try both and see what works best. Yes, my inexperience let me forget about what kind of matrix I needed to define. I hope I've learned my lesson for next time. 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 ________________________________ From: Matthew Knepley [[email protected]] Sent: Wednesday, August 06, 2014 8:30 PM To: Francis Poulin Cc: [email protected] Subject: Re: [petsc-users] trying to solve a lot of ksp problems quickly On Wed, Aug 6, 2014 at 6:50 PM, Francis Poulin <[email protected]<mailto:[email protected]>> wrote: 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. 1) For this to make sense, you have to be doing many many more than N timesteps, where N is the number of unknowns. Why not just solve this equation at each timestep? 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]) 2) Storing the inverse in a sparse matrix does not make sense. Its likely dense. Also, you have not preallocated so it will be slow. Use a dense matrix and both problems will go away. Thanks, Matt 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]<mailto:[email protected]> Web: https://uwaterloo.ca/poulin-research-group/ Telephone: +1 519 888 4567 x32637<tel:%2B1%20519%20888%204567%20x32637> -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener
