[petsc-users] trying to solve a lot of ksp problems quickly

[Matthew Knepley]

On Wed, Aug 6, 2014 at 6:50 PM, Francis Poulin <fpoulin at uwaterloo.ca> 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:           fpoulin at uwaterloo.ca
> Web:            https://uwaterloo.ca/poulin-research-group/
> Telephone:  +1 519 888 4567 x32637
>
>


-- 
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20140806/167ce27a/attachment.html>