Hello,

I am looking for a library which has a direct or iterative solver for singular linear systems. This kind of systems results when the Pressure Poisson Equation with Neumann BC are discretized. I have been using an approach of setting a Dirichlet condition at an interior node and after that solve another Poisson equation in order to remove the errors around that node at which the Dirichlet condition was imposed.

Although this scheme works very weel, I would like to solve the Singular PPE without this approach. I do not know if a Least Square scheme is suitable and often used for this kind of equation.

Please, I need some advice and suggestion regards how to solve the PPE.

This is still on the first page.

http://www.cfd-online.com/Forum/main.cgi?read=61695

I have already seen and answered that post, but I would like to know if a Least Square procedure can resolve the singular problem present in the Pressure Poisson Equation.

Thanks.

Sorry, I missed that.

Least squares could be used, but it's wildly overkill since you know the null space.

A standard procedure that works well is to apply an iterative method as described here:

http://www.cfd-online.com/Forum/main.cgi?read=61700

Some packages (like PETSc) provide support for solving singular systems where you know the null space using the procedure I described generalized to arbitrary solvers and preconditioners. With a direct solver, I don't think there is any disadvantage to fixing the pressure at a point, the main disadvantage with an iterative solver is that it perturbs the spectrum and breaks symmetry.

For more complicated singular systems there are more elaborate schemes such as in this paper:

<p style="font-family: monospace"> @article{reichel2005bfg, title={{Breakdown-free GMRES for Singular Systems}}, author={Reichel, L. and Ye, Q.}, journal={SIAM JOURNAL ON MATRIX ANALYSIS AND APPLICATIONS}, volume={26}, number={4}, pages={1001}, year={2005}, publisher={SIAM SOCIETY FOR INDUSTRIAL AND APPLIED} }

Hi Jed,

In the PETSc package, there is a object called MatNullSpace which must be created with a set of orthonormal vectors which forms the null space of the operator A in Ax = b.

However, I can not figure out how to find the orthonormal base vectors which form the null space of the operator L from Lp = DV*/dt (the resulting linear system of the PPE laplacian(p) = div(V*)/dt).

Is it straightforward?

Many thanks.

HI Jed,

Correct me if I am wrong. For Laplacian operators with full Neumann BC, I have only a constant vector as basis of the nullspace? Which means that in the PETSc setup, I have the parameter "has_cnst" equal true, the parameter "n" equal zero and the set of basis vector "vecs" as a NULL pointer, right? See below.

--- MatNullSpaceCreate --- Creates a data structure used to project vectors out of null spaces.

--- Synopsis ---

MatNullSpaceCreate(MPI_Comm comm,PetscTruth has_cnst,PetscInt n,const Vec vecs[],MatNullSpace *SP)

--- Input Parameters ---

comm - the MPI communicator associated with the object

has_cnst - PETSC_TRUE if the null space contains the constant vector; otherwise PETSC_FALSE

n - number of vectors (excluding constant vector) in null space

vecs - the vectors that span the null space (excluding the constant vector); these vectors must be orthonormal. These vectors are NOT copied, so do not change them after this call. You should free the array that you pass in.

I have the parameter "has_cnst" equal true, the parameter "n" equal zero and the set of basis vector "vecs" as a NULL pointer, right?

Yes, this is correct.

Dear Jed,

I am really glad for your help. Now, I must setup the PETSc and start to work.

I have installed the most recent PETSc package in Ubuntu via the Synaptic Package, but I am having some problems when compiling some examples. I will try to find out what is going on.

Anyway, many thanks.

Allan

I actually recommend building from source, even if you are on a platform that has a packaged PETSc. Use the options <p style="font-family: monospace">--enable-shared --download-ml --download-hypre --download-umfpack --download-mumps
to give you a reasonable setup (two parallel algebraic multigrid preconditioners and two direct solvers, Mumps is parallel). Now run any PETSc program with <p style="font-family: monospace">-pc_type hypre -pc_hypre_type boomeramg
or <p style="font-family: monospace">-pc_type ml -mat_no_inode
for AMG preconditioning. For production runs or benchmarks, build an optimized version by compiling a different PETSC_ARCH with <p style="font-family: monospace">--with-debugging=0
I hope this helps.

Hi Jed,

I'm used to implement my C/C++ codes in the following way:

- include headers files - specify the path for the header and library files - compile and link the code

I have spent a lot of time trying to do it with the PETSc library but unfortunatelly it seems that I can't do that. It seems that I always have to use the makefile schemes which I would like to avoid.

So, is there any way that I can implement my code in an IDE (e.g. Eclipse), include the necessary header files in the code, specify the path for the header and library files in the IDE and finally be happy?

The tricky bit is that PETSc can optionally link against a huge number (> 40) of external libraries which may be in different locations. It's not okay to ask a user to deal with all these directories by hand. I have a FindPETSc.cmake module which I've been using for a while, but just learned it's not as portable as I had thought (I'll post a link here when that's fixed). CMake can use the native build system of whatever platform you are on and integrates with many IDE.

If you run make getlinklibs and make getincludedirs you'll get working flags. Often these can be compressed a bit (and they're better in PETSc-dev, 2.4 is due out in a month or so). If you have a bare-bones install, you may be able to just set -I${PETSC_DIR}/bmake/${PETSC_ARCH} -I${PETSC_DIR}/include and -L${PETSC_DIR}/lib/${PETSC_ARCH} -lpetscts -lpetscsnes -lpetscksp -lpetscdm -lpetscmat -lpetscvec -lpetsc, use the mpi compilers, and make sure Lapack/BLAS gets linked in somewhere. The portable thing is to dump the whole gory output of `make getlinklibs' on your linker (or IDE).

Hi Jed,

In previous posts, you have shown how to remove the null space of the singular coefficient matrix A from the singular linear system Ax = b.

However, the post was closed (I'm wondering why!) and would like to know if you have some note about it. Or if you could write the procedure again.

Thanks in advance.

Allan

Sorry this went so long, I didn't see your reply. I can still find my old post in the archives, but here it is again.

Suppose the operator is <code>A</code> and the orthonormal colums of <code>Z</code> form the null space of <code>A</code>. Let <code>Q</code> be the projection into the orthogonal complement of <code>Z</code>:


Q x = x - (x'Z)Z

Now suppose you are solving the right preconditioned system

P A x = P b
with GMRES, hence we will be expanding the solution <code>x</code> in the space

{Pb, (PA)Pb, (PA)^2Pb, ...}
and we would like <code>x</code> to have no component in the null space <code>(x = Qx)</code>.

The natural way to achieve this is to replace every preconditioner application

y = P b
with

y = Q P b
hence the Krylov space becomes

{QPb, (QPA) QPb, (QPA)^2 QPb, ...}
While this iteration can break down (there are a few papers on robust GMRES for singular systems) it is reasonably reliable. Note that the right hand side should be consistent (usually just solve <code>Ax=Qb</code>) if you expect the converged solution to actually solve <code>Ax=b</code>. For the pressure Poisson, <code>Z</code> is a single constant column vector.

[olatayo]

Pls I need help. How †? derive pressure poisson equation from normalized momentum equation in cylindrical coordinate. I have spent so much time ? this ad my carreer is a stake. I m limited because I had an accident. Pls help M?????e????? .olatayo005@yah.com