Yep - vector_mult, reciprocal, and pointwise_mult are how we do explicit updating in MOOSE:
https://github.com/idaholab/moose/blob/next/framework/src/timeintegrators/ActuallyExplicitEuler.C#L182 Derek On Mon, May 20, 2019 at 2:38 PM Yuxiang Wang <yw...@virginia.edu> wrote: > 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 > _______________________________________________ Libmesh-users mailing list Libmesh-users@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/libmesh-users
