Faster way to invert block-diagonal systems from DG discretization

Hello again,

I would like to get a faster inversion of a local block-diagonal system, which I am using for time-stepping.
The bilinear form is a += (grad(u)*grad(v)+lam*v+u*mu)*dx
where u, v belong to DG-Pk space, and lam and mu belong to DG-P0 space.
The matrix is symmetric and block diagonal. What is the fastest way to invert this matrix?
The following direct approach is a bit to slow

inva = a.mat.Inverse(fes.FreeDofs())

Attached is the code.



Here I found a solution, might be 3-4 times faster than sparse-direct solver:
First create a block Jacobi preconditioner (which is actually the exact solver), and then apply one-step smoothing. This seems to be the fastest way that I can think of…

Yes, I matrix-vector multiply with the block-Jacobi should be the fastest choice, it uses dense matrix-vector for the small blocks, and is parallel.

But double-check whether the geom_free for your bilinear-form is really doing what you want.


The geom_free flag is very mysterious to me. I may have to post another question on the subject. As for the current matrix problem, I think it is of no use, and should be turned off.

Well, the matrix-free approach looks quite promising for time domain DG simulation as was demonstrated in the i-tutorials 3.3 Discontinuous Galerkin Discretizations for linear transport — NGS-Py 6.2.2302 documentation
But, I need to restructure the code quite a bit to do geom_free…