sparsity pattern + matrix quality

[Guosheng Fu]

Hello,

Two questions:
Q1:
I am wondering how to access info. of the assembled fem matrix properties, e.g., sparsity pattern, non-zero entries, conditioning.

---

Backgroud:
I am using P1 Hdiv-HDG of Cockburn/Sayas (see attached code) to solve the Stokes flow on a unit square with 2128128=32768 elements.
I used stadic condensation so that the global linear system consists of about 0.2M golbal dofs and 0.4M local dofs.
The solver is CG with direct solver preconditioner...
pardiso takes 2 iterations to converge, cost about 1m...
umfpack fails to converge due to shortage in memory...(I killed the process when memory usage hit 7G)

Then, I tried two things:
1, Solve a poisson equation with P5 CG fem, umfpack can solve the linear system of size about 1.6M dof within 20s (memory usage is about 6G)

2, Solve the same Hdiv-HDG stoke problem without static condensation, this time pardiso preconditioner fail to converge!!

I don't know what's going on for my stokes linear system...I guess the matrix quality is really bad

---

Q2:
How about the in house Stokes iterative solver with BDDC preconditioner?
I don't know which preconditioner to use and how to apply the QMRSolver in python.

Best,
Guosheng

[schruste]

Hi,

you can check the dofs with V.ndof (for V a FESpace). This is what you know. If you have a bilinear form a you can get the matrix with a.mat (also clear, I guess). Now you can use a.mat.AsVector() to get the entries of the sparse matrix as a vector. Now checking for the length of this vector gives you the number of nonzero entries. a.mat.AsVector().size should already do the trick.
What do you mean by sparsity pattern? Would you like to have something like spy in matlab?!
For these and many other things you can simply use python tools, e.g. scipy.
With

rows,cols,vals = a.mat.COO()

you get the information of the matrix as row index, column index and entry values of the matrix which you can then feed into scipy

import scipy.sparse as sp
A = sp.csr_matrix((vals,(rows,cols))

Now A has information like A.nnz (number of nonzero entries) and you can plot the sparsity pattern with matplotlib:

import matplotlib.pylab as plt
plt.spy(A)
plt.show()

Getting the condition number is not so easy from scipy (I think!?). If your matrix is small you can do

import numpy as np
np.linalg.cond(A.todense())

for larger systems I don't know the python way to do it.

For NGSolve you can get an estimated condition number of the preconditioned matrix using the "test" flag for a Preconditioner, e.g.

c = Preconditioner(bf=a, type="local", flags = {"test" : True})
c.Update()

will give you the condition number of the jacobi-preconditioned system, approximately.

You can also try to use the (p-version-)BDDC-preditioner from NGSolve. Therefore take

c = Preconditioner(bf=a, type="bddc", flags = {"test" : True})

Note that the BDDC-preconditioner has to be defined before a.Assemble() is called. Further note that the BDDC-preconditioner will use the BDDC-pattern that comes from your (Compound) FESpace. The coarse-grid system is constructed from degrees of freedom that have couplingtype = WIREBASKET. That implies that you should mark your lowest order pressure degrees of freedom not only as COUPLING_TYPE.INTERFACE_DOF but as COUPLING_TYPE.WIREBASKET_DOF so that they are solved for also in the coarse grid problem of the BDDC.

Hope this helps for now.

Best,
Christoph

P.S.: I will have a look at your code later...

[Guosheng Fu]

Christoph,

Yeah, Scipy & Numpy are what I want... :)

But... an assembled fem matrix a.mat is singular due to dirichelt BC (I am thinking of a laplacian + dirichelt bc) not treated explicited in the bilinear form... is ther a way to get the "correct" matrix after applying dirichelt BC? then, I can use numpy to estimate at lease small matrix conditing...

Regarding NGsolve conditioning estimation, it seems to be only working for SPD systems becase it's computing eigenvalues, not singular values....it fails to estimate, say, the mixed laplacian condition number.....

Regrading bddc preconditioner, it looks great for the cg-laplaican solver.
I would like to play with this preconditioner for saddle point systems, eg, the mixed laplacian... but then I need another iterative solver, say QMR as you mentioned or GMRES as I wanted... can you give me a hand on the iterative solver ;)

Lastly, is there a documentation on preconditioners implemented in ngsolver that I can check out? You seems to have a lot!!!

Best,
Guosheng

[schruste]

Guosheng,

"correct" matrix:
you have the "full" matrix and you have the bitarray freedofs-array from the finite element space. If you want to remove rows and column you can surely do that in scipy and than analyse the resulting matrix. Is that what you want? If yes, you have to look for the scipy-solution to this problem. I am not an expert on this.

NGSolve conditioning estimates: true

I think you can directly use QMRSolver (similar interface as the CGSolver), check help(QMRSolver) for details. You can also write your own iterative solver, see also the cg-example in py_tutorials.

preconditioners in NGSolve: The source code is the documentation I guess ;( (e.g. preconditioner.hpp,...). The preconditioning system is modular in NGSolve, i.e. you can add different preconditioners yourself and register them to NGSolve. Most preconditioners that are available are for spd-systems only. Here is a list of the most established preconditioners in NGSolve:
* jacobi and block jacobi preconditioners ("local" +/- additional "block" flag)
* multigrid preconditioner ("multigrid") (with many options for the smoothers)
* bddc ("bddc")
* direct solver as preconditioner ("direct")

Problem-dependet preconditioners can be added separately and are not in the core code.

If you want block-preconditioners you would need to write them yourself, but this is actually not too hard, but requires a bit of c++-coding.

Best,
Christoph

[Guosheng Fu]

Hi Christoph,

How about static condensation + CGSolver or QMRSolver?
These iterative solvers need arguments for matrix and preconditioner, it seems not supported yet. Correct me if I am wrong

Best,
Guosheng

[schruste]

Hi Guosheng,

static condensation with CGSolver or any other solver works out of the box. If you use static condensation the resulting matrix will not have the non-zero entries of LOCAL_DOFs, i.e. the global matrix is already the Schur complement matrix. Preconditioners then also act on this Schur complement.

There is one thing to note about the condensated matrix in NGSolve:
If you have N dofs globally (before static condensation) the statically condensated matrix will still have the dimension NxN. However, there will be empty rows and columns for the dofs that are condensated out.
For the Schur complement operations you then have other sparse matrices of dimension NxN for the inner solves/harmonic extensions/... .

Best,
Christoph

[schruste]

of course you have to apply the modifications to your linear system when doing static condensation, see also here:
https://gitlab.asc.tuwien.ac.at/jschoeberl/ngsolve-docu/wikis/ngspy-howto-staticconensation

Note that in the example there, the line

u.vec.data = a.mat.Inverse(fes.FreeDofs(True)) * f.vec

can be replaced with your CGSolver or your QMRSolver...

[Guosheng Fu]

Then, I guess my stupid question is how to replace it...I failed
Can you give me a working sample code :)

[schruste]

try the attached example:

[Guosheng Fu]

Christoph,

Oh... This works... But directly changing CGSolver to QMRSolver does not converge.... I need to look into the QMRSolver implementation at some point...but not right now... ;)

Well, on the other hand, the reason I looked into the solver is because pardiso fail to give me the correct solution for the Hdiv-HDG Stokes/Brinkman problem...Philip suggested umfpack, and it did work, but too memory consuming and I went out of memory for condensed size about 0.1 M (how could this happen?!).... at this point, I realized that direct solver is not what I want.....

What I am doing right now, is using pardiso (which is really fast comparing to umfpack) for a (small) pressure regularized Stokes problem (added a pressure mass term (epspq) to the bilinear form) as preconditioner for the CG solver.... it works fine, converge in three iterations always (well, some times fail, then after adding more pressure mass, it would work), and gives me correct order of convergence, at the price of assemble fem matrices twices...

The conclusion is, I am all fine for now for the solver part...
Thank you for the timely help! ;0

Best always,
Guosheng

[schruste]

A follow up on linear systems and interfacing to python-tools:
I summarized the above mentioned simple tools at
https://gitlab.asc.tuwien.ac.at/jschoeberl/ngsolve-docu/wikis/ngspy-numpy
In the last paragraph I also added a short note on how to use scipys iterative solvers. You could try your problem with their iterative solvers, e.g. gmres.
Note that you can easily incorporate preconditioning as well.
Attached is a small example for a DG Poisson test case.

Best,
Christoph

[Guosheng Fu]

Here is the code...