Thank you Jed & John! That's extremely helpful. John - thank you so much for the insight and pointer to reciprocal() and pointwise_mult()! As for adding additional NumericVector to the system - do you happen to know any example or code snippet that did this? What I am doing now is creating a new ExplicitSystem and using its solution vector as the storage, which handles the initialization quite nicely but maybe is an overkill.
Best, Shawn On Mon, May 20, 2019 at 7:11 AM Jed Brown <j...@jedbrown.org> wrote: > John Peterson <jwpeter...@gmail.com> writes: > > > On Sun, May 19, 2019 at 9:15 PM Yuxiang Wang <yw...@virginia.edu> wrote: > > > >> Dear all, > >> > >> Sorry for the spam and being new to the field. > >> > >> I am currently trying to implement an elastodynamics problem with > explicit > >> method (central difference method, to be specific). I am planning to use > >> lumped mass matrix (and Rayleigh damping when needed), so the system > matrix > >> will be simply a diagonal matrix. Solving is therefore trivial - I just > >> need to do per-element division of the rhs by the diagonal entries to > get > >> the solution vector. > >> > >> For this problem, I have the option of just treating it as an implicit > >> system - fill the system.matrix with only diagonal components and call > the > >> PETSc LU solver. This is very easy thanks to a lot of the libmesh > >> infrastructure available. If I do so, should I be concerned about a > >> significant slowdown? Or would PETSc be smart enough to realize that > this > >> is already diagonal matrix and be efficient in solving it? > >> > > > > Yes, I'm pretty sure this will be significantly slower than doing a > > "manual" inversion by taking the reciprocal of the matrix diagonal. I > don't > > know of anything in PETSc's LU solver that will detect this particular > > special case. For an explicit code where you want every timestep to be as > > fast as possible, it will likely be prohibitively slow. > > If you have a sparse matrix with only the diagonal nonzero, then LU will > create no fill and it'll actually be pretty fast (likely hard to measure > compared to the cost of evaluating the explicit part), but -pc_type > jacobi would also be an exact solver and is more precisely what you > need. Note that the default bjacobi/ilu is also an exact solver in this > circumstance, as are many other preconditioners. > > > My other choice would be to create a NumericVector myself to store the > >> diagonal system matrix entries, and perform the per-element division. > This > >> would take more work and will not be using the already well-tested > libmesh > >> infrastructure, so I am trying to see whether I can get away with doing > >> this without compromising on the performance. > > > > > > This would probably work best. You can add additional vectors to Systems > > and assemble into them similarly to the way you assemble the right-hand > > side vector. Also note that NumericVectors have the reciprocal() API, > which > > will allow you to quickly compute the inverse, as well as > pointwise_mult() > > which should allow you to quickly apply it. > > > > -- > > John > > > > _______________________________________________ > > Libmesh-users mailing list > > Libmesh-users@lists.sourceforge.net > > https://lists.sourceforge.net/lists/listinfo/libmesh-users > -- Yuxiang "Shawn" Wang, PhD yw...@virginia.edu +1 (434) 284-0836 _______________________________________________ Libmesh-users mailing list Libmesh-users@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/libmesh-users
