[petsc-users] Solving singular systems with petsc
Nidish
nb25 at rice.edu
Mon Aug 17 13:42:44 CDT 2020
Thank you for the email, Zak,
I have not looked into continuation methods for this problem yet, and
I'd love to hear your thoughts on it!
The problem I have is two-fold:
1. Obtaining the eigenvectors of the system
K.v = lam*M.v
where the K and M are stiffness and mass matrices coming from a
finite element model without any Dirichlet boundary conditions
(i.e., having 6 Rigid body modes) and with a few RBE3 constraints
(introducing degrees of freedom with "zero mass" in the mass
matrix). So the system has, in addition to the 6 rigid body modes, a
few "spurious" null vectors coming from these RBE3 constraints
(which are enforced using Lagrange multipliers).
2. Conducting a linear solve of the system:
K.x = b
where K is from above.
What I'm trying to do with both of these is to conduct a
Hurty/Craig-Bampton Component Mode Synthesis, which is like a Schur
Condensation with a few "fixed interface modal" DoFs added to the
reduced system.
So both can be solved by obtaining the vectors that span the null space
of the system.
I've created two separate threads for the two problems because I felt
both are slightly different questions. But I understand now that the
whole thing boils down to obtaining the nullspaces.
Thank you,
Nidish
On 8/17/20 1:33 PM, zakaryah wrote:
> Hi Nidish,
>
> I may not fully understand your problem, but it sounds like you could
> benefit from continuation methods. Have you looked into this? If it's
> helpful, I have some experience with this and I can discuss with you
> by email.
>
> Cheers, Zak
>
> On Mon, Aug 17, 2020 at 2:31 PM Nidish <nb25 at rice.edu
> <mailto:nb25 at rice.edu>> wrote:
>
> Thankfully for this step, the matrix is not dense. But thank you.
>
> Nidish
>
> On 8/17/20 8:52 AM, Jose E. Roman wrote:
> >
> >> El 17 ago 2020, a las 14:27, Barry Smith <bsmith at petsc.dev
> <mailto:bsmith at petsc.dev>> escribió:
> >>
> >>
> >> Nidish,
> >>
> >> Your matrix is dense, correct? MUMPS is for sparse matrices.
> >>
> >> Then I guess you could use Scalapack
> http://netlib.org/scalapack/slug/node48.html#SECTION04323200000000000000
> to do the SVD. The work is order N^3 and parallel efficiency may
> not be great but it might help you solve your problem.
> >>
> >> I don't know if SLEPc has an interface to Scalapack for
> SVD or not.
> > Yes, SLEPc (master) has interfaces for ScaLAPACK and Elemental
> for both SVD and (symmetric) eigenvalues.
> >
> > Jose
> >
> >> Barry
> >>
> >>
> >>
> >>
> >>
> >>
> >>
> >>> On Aug 17, 2020, at 2:51 AM, Jose E. Roman <jroman at dsic.upv.es
> <mailto:jroman at dsic.upv.es>> wrote:
> >>>
> >>> You can use SLEPc's SVD to compute the nullspace, but it has
> pitfalls: make sure you use an absolute convergence test (not
> relative); for the particular case of zero singular vectors,
> accuracy may not be very good and convergence may be slow (with
> the corresponding high computational cost).
> >>>
> >>> MUMPS has functionality to get a basis of the nullspace, once
> you have computed the factorization. But I don't know if this is
> easily accessible via PETSc.
> >>>
> >>> Jose
> >>>
> >>>
> >>>
> >>>> El 17 ago 2020, a las 3:10, Nidish <nb25 at rice.edu
> <mailto:nb25 at rice.edu>> escribió:
> >>>>
> >>>> Oh damn. Alright, I'll keep trying out the different options.
> >>>>
> >>>> Thank you,
> >>>> Nidish
> >>>>
> >>>> On 8/16/20 8:05 PM, Barry Smith wrote:
> >>>>> SVD is enormously expensive, needs to be done on a full
> dense matrix so completely impractical. You need the best tuned
> iterative method, Jose is the by far the most knowledgeable about
> that.
> >>>>>
> >>>>> Barry
> >>>>>
> >>>>>
> >>>>>> On Aug 16, 2020, at 7:46 PM, Nidish <nb25 at rice.edu
> <mailto:nb25 at rice.edu>> wrote:
> >>>>>>
> >>>>>> Thank you for the suggestions.
> >>>>>>
> >>>>>> I'm getting a zero pivot error for the LU in slepc while
> calculating the rest of the modes.
> >>>>>>
> >>>>>> Would conducting an SVD for just the stiffness matrix and
> then using the singular vectors as bases for the nullspace work? I
> haven't tried this out just yet, but I'm wondering if you could
> provide me insights into whether this will.
> >>>>>>
> >>>>>> Thanks,
> >>>>>> Nidish
> >>>>>>
> >>>>>> On 8/16/20 2:50 PM, Barry Smith wrote:
> >>>>>>> If you know part of your null space explicitly (for
> example the rigid body modes) I would recommend you always use
> that information explicitly since it is extremely expensive
> numerically to obtain. Thus rather than numerically computing the
> entire null space compute the part orthogonal to the part you
> already know. Presumably SLEPc has tools to help do this, naively
> I would just orthogonalized against the know subspace during the
> computational process but there are probably better ways.
> >>>>>>>
> >>>>>>> Barry
> >>>>>>>
> >>>>>>>
> >>>>>>>
> >>>>>>>
> >>>>>>>> On Aug 16, 2020, at 11:26 AM, Nidish <nb25 at rice.edu
> <mailto:nb25 at rice.edu>> wrote:
> >>>>>>>>
> >>>>>>>> Well some of the zero eigenvectors are rigid body modes,
> but there are some more which are introduced by
> lagrange-multiplier based constraint enforcement, which are non
> trivial.
> >>>>>>>>
> >>>>>>>> My final application is for a nonlinear simulation, so I
> don't mind the extra computational effort initially. Could you
> have me the suggested solver configurations to get this type of
> eigenvectors in slepc?
> >>>>>>>>
> >>>>>>>> Nidish
> >>>>>>>> On Aug 16, 2020, at 00:17, Jed Brown <jed at jedbrown.org
> <mailto:jed at jedbrown.org>> wrote:
> >>>>>>>> It's possible to use this or a similar algorithm in
> SLEPc, but keep in mind that it's more expensive to compute these
> eigenvectors than to solve a linear system. Do you have a
> sequence of systems with the same null space?
> >>>>>>>>
> >>>>>>>> You referred to the null space as "rigid body modes".
> Why can't those be written down? Note that PETSc has convenience
> routines for computing rigid body modes from coordinates.
> >>>>>>>>
> >>>>>>>> Nidish <
> >>>>>>>> nb25 at rice.edu <mailto:nb25 at rice.edu>
> >>>>>>>>> writes:
> >>>>>>>>
> >>>>>>>> I just use the standard eigs function
> (https://www.mathworks.com/help/matlab/ref/eigs.html
> >>>>>>>> ) as a black box. I think it uses a lanczos type method
> under the hood.
> >>>>>>>>
> >>>>>>>> Nidish
> >>>>>>>>
> >>>>>>>> On Aug 15, 2020, 21:42, at 21:42, Barry Smith <
> >>>>>>>> bsmith at petsc.dev <mailto:bsmith at petsc.dev>
> >>>>>>>>> wrote:
> >>>>>>>>
> >>>>>>>> Exactly what algorithm are you using in Matlab to get the
> 10 smallest
> >>>>>>>> eigenvalues and their corresponding eigenvectors?
> >>>>>>>>
> >>>>>>>> Barry
> >>>>>>>>
> >>>>>>>>
> >>>>>>>>
> >>>>>>>> On Aug 15, 2020, at 8:53 PM, Nidish <nb25 at rice.edu
> <mailto:nb25 at rice.edu>
> >>>>>>>>> wrote:
> >>>>>>>> The section on solving singular systems in the manual
> starts with
> >>>>>>>>
> >>>>>>>> assuming that the singular eigenvectors are already known.
> >>>>>>>>
> >>>>>>>>
> >>>>>>>> I have a large system where finding the singular
> eigenvectors is not
> >>>>>>>>
> >>>>>>>> trivially written down. How would you recommend I proceed
> with making
> >>>>>>>> initial estimates? In MATLAB (with MUCH smaller
> matrices), I conduct an
> >>>>>>>> eigensolve for the first 10 smallest eigenvalues and take the
> >>>>>>>> eigenvectors corresponding to the zero eigenvalues from
> this. This
> >>>>>>>> approach doesn't work here since I'm unable to use SLEPc
> for solving
> >>>>>>>>
> >>>>>>>>
> >>>>>>>> K.v = lam*M.v
> >>>>>>>>
> >>>>>>>> for cases where K is positive semi-definite (contains a
> few "rigid
> >>>>>>>>
> >>>>>>>> body modes") and M is strictly positive definite.
> >>>>>>>>
> >>>>>>>>
> >>>>>>>> I'd appreciate any assistance you may provide with this.
> >>>>>>>>
> >>>>>>>> Thank you,
> >>>>>>>> Nidish
> >>>>>>>>
> >>>>>> --
> >>>>>> Nidish
> >>>> --
> >>>> Nidish
> --
> Nidish
>
--
Nidish
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20200817/e84f11dd/attachment-0001.html>
More information about the petsc-users
mailing list