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

[Francis Poulin]

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:           fpoulin at uwaterloo.ca
Web:            https://uwaterloo.ca/poulin-research-group/
Telephone:  +1 519 888 4567 x32637

________________________________
From: Matthew Knepley [knepley at gmail.com]
Sent: Wednesday, August 06, 2014 8:30 PM
To: Francis Poulin
Cc: petsc-users at mcs.anl.gov
Subject: Re: [petsc-users] trying to solve a lot of ksp problems quickly

On Wed, Aug 6, 2014 at 6:50 PM, Francis Poulin <fpoulin at uwaterloo.ca<mailto: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<mailto:fpoulin at uwaterloo.ca>
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20140807/679e1355/attachment-0001.html>