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 _______________________________________________ Libmesh-users mailing list Libmesh-users@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/libmesh-users
