[petsc-users] Rigid body nullspace for Stokes operator

[Jed Brown]

There is MatSetNearNullSpace, which is used to ensure an approximation property in coarse levels of multigrid. That should always be set when dealing with problems like this, regardless of boundary conditions. Separately, there is MatSetNullSpace, which you should only use if you are solving a singular system with that null space.

Note that you need at least three non-colinear points to constrain the null space so if you, for example, have a Dirichlet condition at only one point or along one straight line, there will still be a null space.

Amneet Bhalla <mail2amneet at gmail.com> writes:

> I think Mark mentioned this earlier, but I want to make sure that the rigid
> body null vectors should be specified only when Neumann boundary conditions
> are used on all boundaries of the domain, correct? Alternatively, if a
> Dirichlet boundary condition is used (on any part of the domain boundary)
> then there is no null space, i.e., the operator is a full rank matrix?
>
> If the above is true, then I think I do not need to specify the rigid body
> null modes because I am using Dirichlet boundary conditions for the
> velocity solver.
>
> On Wed, Oct 30, 2024 at 12:28 PM Jed Brown <jed at jedbrown.org> wrote:
>
>> Yes to 6 null vectors in 3D, 3 null vectors in 2D. The center of mass does
>> not need to be identified because you can algebraically orthogonalize
>> (lines 411-420 here).
>>
>>
>> https://urldefense.us/v3/__https://gitlab.com/petsc/petsc/-/blob/main/src/dm/impls/plex/plexfem.c?ref_type=heads*L377-425__;Iw!!G_uCfscf7eWS!cUe-WWTFiWh16JyZXhgGBBJ9kXbQpDaKFqY17p826i1tQk6dtu_ccEOtB9z_w_rXWK8hYv5SqNXhXSmZGcc$ 
>>
>> See also this implementation with raw coordinates. GAMG orthogonalizes
>> within each aggregate (in a later phase of the algorithm) so global
>> orthogonalization is not necessary.
>>
>>
>> https://urldefense.us/v3/__https://gitlab.com/petsc/petsc/-/blob/main/src/ksp/pc/impls/gamg/agg.c?ref_type=heads*L387__;Iw!!G_uCfscf7eWS!cUe-WWTFiWh16JyZXhgGBBJ9kXbQpDaKFqY17p826i1tQk6dtu_ccEOtB9z_w_rXWK8hYv5SqNXhcaBBq7s$ 
>>
>> Amneet Bhalla <mail2amneet at gmail.com> writes:
>>
>> > I think the nullspace for the velocity operator is of the form
>> >
>> > vnull = U + ω × r
>> > in which U is a rigid body velocity and ω is the rigid body rotational
>> > velocity, and r is the radius vector from the center of mass. I believe I
>> > need to construct 6 nullspace vectors in 3D and 3 nullspace vectors in
>> 2D.
>> > Sounds correct? Also does the center of mass coordinates matter when
>> > defining r?
>> >
>> > On Wed, Oct 30, 2024 at 7:53 AM Amneet Bhalla <mail2amneet at gmail.com>
>> wrote:
>> >
>> >> @Mark: Is there some document/paper that I can follow to check the
>> algebra
>> >> of these zero eigenvectors/null space modes?
>> >>
>> >> @Jed : We use a projection method preconditioner to solve the coupled
>> >> velocity pressure system as described here (
>> >> https://urldefense.us/v3/__https://www.sciencedirect.com/science/article/pii/S0021999123004205__;!!G_uCfscf7eWS!cUe-WWTFiWh16JyZXhgGBBJ9kXbQpDaKFqY17p826i1tQk6dtu_ccEOtB9z_w_rXWK8hYv5SqNXhOpVmAEw$ ).
>> It
>> >> is an approximation of the Schur complement. As a part of projection
>> >> preconditioner, we need to solve just the momentum equation separately
>> >> without considering the pressure part.
>> >>
>> >> On Tue, Oct 29, 2024 at 8:03 PM Jed Brown <jed at jedbrown.org> wrote:
>> >>
>> >>> And to be clear, we recommend using fieldsplit Schur to separate the
>> >>> pressure and velocity part (there are many examples). Applying AMG
>> directly
>> >>> to the saddle point problem will not be a good solver because the
>> >>> heuristics assume positivity and do not preserve inf-sup stability
>> (nor do
>> >>> standard smoothers).
>> >>>
>> >>> Mark Adams <mfadams at lbl.gov> writes:
>> >>>
>> >>> > This is linear elasticity and there are 6 "null" vectors (if you
>> removed
>> >>> > Dirichlet boundary conditions these are eigenvectors with zero
>> >>> eigenvalue):
>> >>> > 3 translations, x, y, z, and three rotatiions xx, yy ,zz.
>> >>> > x = (1,0,0,1,0,0,1,0 ...)
>> >>> > and xx is something like (0, z_1, -y_1, 0, z_2, -y_2, ...) where z_1
>> is
>> >>> the
>> >>> > z coordinate of the first vertex, etc.
>> >>> >
>> >>> > Mark
>> >>> >
>> >>> > On Tue, Oct 29, 2024 at 3:47 PM Amneet Bhalla <mail2amneet at gmail.com
>> >
>> >>> wrote:
>> >>> >
>> >>> >> Hi Mark,
>> >>> >>
>> >>> >> Thanks! I am not sure how to construct null space and zero energy
>> modes
>> >>> >> manually for this operator. Is there some theory or documentation I
>> can
>> >>> >> follow to figure out what the null space and zero energy modes look
>> >>> like
>> >>> >> for this operator? Once I know what these are in symbolic form, I
>> >>> think I
>> >>> >> should be able to construct them manually.
>> >>> >>
>> >>> >> Best,
>> >>> >> --Amneet
>> >>> >>
>> >>> >> On Tue, Oct 29, 2024 at 7:35 AM Mark Adams <mfadams at lbl.gov> wrote:
>> >>> >>
>> >>> >>> Oh my mistake. You are using staggered grids. So you don't have a
>> >>> block
>> >>> >>> size that hypre would use for the "nodal" methods.
>> >>> >>> I'm not sure what you are doing exactly, but try hypre and you
>> could
>> >>> >>> create the null space, zero energy modes, manually, attach to the
>> >>> matrix
>> >>> >>> and try GAMG.
>> >>> >>> You can run with '-info :pc' and grep on GAMG to see if GAMG is
>> >>> picking
>> >>> >>> the null space up (send this output if you can't figure it out).
>> >>> >>>
>> >>> >>> Thanks,
>> >>> >>> Mark
>> >>> >>>
>> >>> >>> On Tue, Oct 29, 2024 at 9:28 AM Mark Adams <mfadams at lbl.gov>
>> wrote:
>> >>> >>>
>> >>> >>>> This coordinate interface is just a shortcut for vertex based
>> >>> >>>> discretizations with 3 dof per vertex, etc. (maybe works in 2D).
>> >>> >>>> You will need to construct the null space vectors manually and
>> >>> attach it
>> >>> >>>> to the matrix. Used by GAMG.
>> >>> >>>>
>> >>> >>>> Note, for hypre you want to use the "nodal" options and it does
>> not
>> >>> use
>> >>> >>>> these null space vectors. That is probably the way you want to go.
>> >>> >>>> eg: -pc_hypre_boomeramg_nodal_coarsen
>> >>> >>>>
>> >>> >>>> I would run with hypre boomerang and -help and grep on nodal to
>> see
>> >>> all
>> >>> >>>> the "nodal" options and use them.
>> >>> >>>>
>> >>> >>>> Thanks,
>> >>> >>>> Mark
>> >>> >>>>
>> >>> >>>>
>> >>> >>>> On Mon, Oct 28, 2024 at 8:06 PM Amneet Bhalla <
>> mail2amneet at gmail.com
>> >>> >
>> >>> >>>> wrote:
>> >>> >>>>
>> >>> >>>>> Hi Folks,
>> >>> >>>>>
>> >>> >>>>> I am trying to solve the momentum equation in a projection
>> >>> >>>>> preconditioner using GAMG or Hypre solver. The equation looks
>> like
>> >>> for
>> >>> >>>>> velocity variable *v* looks like:
>> >>> >>>>>
>> >>> >>>>>
>> >>> >>>>>  [image: Screenshot 2024-10-28 at 4.15.17 PM.png]
>> >>> >>>>>
>> >>> >>>>> Here, μ is spatially varying dynamic viscosity and λ is spatially
>> >>> >>>>> varying bulk viscosity. I understand that I need to specify rigid
>> >>> body
>> >>> >>>>> nullspace modes to the multigrid solver in order to accelerate
>> its
>> >>> >>>>> convergence. Looking into this routine
>> >>> MatNullSpaceCreateRigidBody() (
>> >>> >>>>>
>> >>>
>> https://urldefense.us/v3/__https://petsc.org/release/manualpages/Mat/MatNullSpaceCreateRigidBody/__;!!G_uCfscf7eWS!bVR6duCoDqPhZrWS-sm1c5qxsFPjZMhdT86AqLpPzWgVy5qoRhd4_Jue2LJOIS6LRrtV2cHGrqger1Yvb-Y5f-0$
>> >>> >>>>> <
>> >>>
>> https://urldefense.us/v3/__https://petsc.org/release/manualpages/Mat/MatNullSpaceCreateRigidBody/__;!!G_uCfscf7eWS!eKqgIJjCdMzIU76f7X65AmGxrU_-lC7W02BMWafJ77DNf_IuQk6O1X3qU1x9Ez8NJ20vZEL-mF6T1yNmDnwv0eWa2w$
>> >>> >),
>> >>> >>>>> I see that I need to provide the coordinates of each node. I am
>> >>> using
>> >>> >>>>> staggered grid discretization. Do I need to provide coordinates
>> of
>> >>> >>>>> staggered grid locations?
>> >>> >>>>>
>> >>> >>>>> Thanks,
>> >>> >>>>> --
>> >>> >>>>> --Amneet
>> >>> >>>>>
>> >>> >>>>>
>> >>> >>>>>
>> >>> >>>>>
>> >>> >>
>> >>> >> --
>> >>> >> --Amneet
>> >>> >>
>> >>> >>
>> >>> >>
>> >>> >>
>> >>>
>> >>
>> >>
>> >> --
>> >> --Amneet
>> >>
>> >>
>> >>
>> >>
>> >
>> > --
>> > --Amneet
>>
>
>
> -- 
> --Amneet