Hi Lisandro, Thanks for your recommendation. Btw, does the poisson eqn arising from fractional step gives a matrix which is SPD? Because my grid's are non-uniform in both x,y directions. Shouldn't that result in a non-symmetric matrix? But I think it's still PD, positive definite. Correct me if I'm wrong.
Thanks Lisandro Dalcin wrote: > Well, after taking into accout Barry's comments, you have have the > following choices. > > * You can use a direct method based on LU factorization using > '-ksp_type preonly -pc_type lu' . This way, PETSc will compute the LU > factors the fist time they are needed; after that, every call to > KSPSolve will reuse those factors. This will work only in sequential > with a default PETSc build, but you could also build PETSc with MUMPS, > and it will let you do the parallel factorization. For MUMPS to > actually work in your matrix, I believe you have to add the following > line: > > MatConvert(A, MATAIJMUMPS, MAT_REUSE_MATRIX, &A); > > after assembling (ie. MatAssembleBegin/End calls) your Poisson matrix. > > > * You can use CG with '-ksp_type cg' (I assume your matrix is SPD, as > it is in a standard fractional step method), and a preconditioner. And > then, I believe the best choice for your application will bee > BoomerAMG. It has a rather high setup cost, but solves are fast. Or > your could use ML, it has less setup costs, but the solvers are a bit > slower. So if you make many timesteps, I would say that BoomerAMG will > pay. > > Finally, if you use the last option, perhaps you can try Paul Fischer > tricks. I tried to add this to KSP's some time ago, but I stoped for > many reasons (the main one, lack of time). You can take a look at > this: > > http://citeseer.ist.psu.edu/492082.html > > A similar (equivalent?) approach is this other one (perhaps a bit > easier to implement, depending on your taste) > doi.wiley.com/10.1002/cnm.743 > > > On 2/5/08, Ben Tay <zonexo at gmail.com> wrote: > >> Hi Lisandro, >> >> I'm using the fractional step mtd to solve the NS eqns as well. I've >> tried the direct mtd and also boomerAMG in solving the poisson eqn. >> Experience shows that for smaller matrix, direct mtd is slightly faster >> but if the matrix increases in size, boomerAMG is faster. Btw, if I'm >> not wrong, the default solver will be GMRES. I've also tried using the >> "Struct" interface solely under Hypre. It's even faster for big matrix, >> although the improvement doesn't seem to be a lot. I need to do more >> tests to confirm though. >> >> I'm now doing 2D simulation with 1400x2000 grids. It's takes quite a >> while to solve the eqns. I'm wondering if it'll be faster if I get the >> inverse and then do matrix multiplication. Or just calling KSPSolve is >> actually doing something similar and there'll not be any speed >> difference. Hope someone can enlighten... >> >> Thanks! >> >> Lisandro Dalcin wrote: >> >>> Ben, some time ago I was doing some testing with PETSc for solving >>> incompressible NS eqs with fractional step method. I've found that in >>> our software and hardware setup, the best way to solve the pressure >>> problem was by using HYPRE BoomerAMG. This preconditioner usually have >>> some heavy setup, but if your Poison matrix does not change, then the >>> sucessive solves at each time step are really fast. >>> >>> If you still want to use a direct method, you should use the >>> combination '-ksp_type preonly -pc_type lu' (by default, this will >>> only work on sequential mode, unless you build PETSc with an external >>> package like MUMPS). This way, PETSc computes the LU factorization >>> only once, and at each time step, the call to KSPSolve end-up only >>> doing the triangular solvers. >>> >>> The nice thing about PETSc is that, if you next realize the >>> factorization take a long time (as it usually take in big problems), >>> you can switch BoomerAMG by only passing in the command line >>> '-ksp_type cg -pc_type hypre -pc_hypre_type boomeramg'. And that's >>> all, you do not need to change your code. And more, depending on your >>> problem you can choose the direct solvers or algebraic multigrid as >>> you want, by simply pass the appropriate combination options in the >>> command line (or a options file, using the -options_file option). >>> >>> Please, if you ever try HYPRE BoomerAMG preconditioners, I would like >>> to know about your experience. >>> >>> Regards, >>> >>> On 2/5/08, Ben Tay <zonexo at gmail.com> wrote: >>> >>> >>>> Hi everyone, >>>> >>>> I was reading about the topic abt inversing a sparse matrix. I have to >>>> solve a poisson eqn for my CFD code. Usually, I form a system of linear >>>> eqns and solve Ax=b. The "A" is always the same and only the "b" changes >>>> every timestep. Does it mean that if I'm able to get the inverse matrix >>>> A^(-1), in order to get x at every timestep, I only need to do a simple >>>> matrix multiplication ie x=A^(-1)*b ? >>>> >>>> Hi Timothy, if the above is true, can you email me your Fortran code >>>> template? I'm also programming in fortran 90. Thank you very much >>>> >>>> Regards. >>>> >>>> Timothy Stitt wrote: >>>> >>>> >>>>> Yes Yujie, I was able to put together a parallel code to invert a >>>>> large sparse matrix with the help of the PETSc developers. If you need >>>>> any help or maybe a Fortran code template just let me know. >>>>> >>>>> Best, >>>>> >>>>> Tim. >>>>> >>>>> Waad Subber wrote: >>>>> >>>>> >>>>>> Hi >>>>>> There was a discussion between Tim Stitt and petsc developers about >>>>>> matrix inversion, and it was really helpful. That was in last Nov. >>>>>> You can check the emails archive >>>>>> >>>>>> http://www-unix.mcs.anl.gov/web-mail-archive/lists/petsc-users/2007/11/threads.html >>>>>> >>>>>> >>>>>> Waad >>>>>> >>>>>> */Yujie <recrusader at gmail.com>/* wrote: >>>>>> >>>>>> what is the difference between sequantial and parallel AIJ matrix? >>>>>> Assuming there is a matrix A, if >>>>>> I partitaion this matrix into A1, A2, Ai... An. >>>>>> A is a parallel AIJ matrix at the whole view, Ai >>>>>> is a sequential AIJ matrix? I want to operate Ai at each node. >>>>>> In addition, whether is it possible to get general inverse using >>>>>> MatMatSolve() if the matrix is not square? Thanks a lot. >>>>>> >>>>>> Regards, >>>>>> Yujie >>>>>> >>>>>> >>>>>> On 2/4/08, *Barry Smith* <bsmith at mcs.anl.gov >>>>>> <mailto:bsmith at mcs.anl.gov>> wrote: >>>>>> >>>>>> >>>>>> For sequential AIJ matrices you can fill the B matrix >>>>>> with the >>>>>> identity and then use >>>>>> MatMatSolve(). >>>>>> >>>>>> Note since the inverse of a sparse matrix is dense the B >>>>>> matrix is >>>>>> a SeqDense matrix. >>>>>> >>>>>> Barry >>>>>> >>>>>> On Feb 4, 2008, at 12:37 AM, Yujie wrote: >>>>>> >>>>>> > Hi, >>>>>> > Now, I want to inverse a sparse matrix. I have browsed the >>>>>> manual, >>>>>> > however, I can't find some information. could you give me >>>>>> some advice? >>>>>> > >>>>>> > thanks a lot. >>>>>> > >>>>>> > Regards, >>>>>> > Yujie >>>>>> > >>>>>> >>>>>> >>>>>> >>>>>> ------------------------------------------------------------------------ >>>>>> Looking for last minute shopping deals? Find them fast with Yahoo! >>>>>> Search. >>>>>> <http://us.rd.yahoo.com/evt=51734/*http://tools.search.yahoo.com/newsearch/category.php?category=shopping> >>>>>> >>>>>> >>>>> >>> >>> >> > > >