Applied Linear Algebra

When using Finite Elements, one needs to invert large linear systems that are matrix/vector representations of Bilinear or Linear forms. In the case of MHD, it is well known that these systems are very ill-conditioned. Understanding the different pathologies and their nature is critical for devising ad-hoc iterative linear solvers.

Let us consider the 1d Poisson equation, and the stiffness matrix. The use of Generalized Locally Theory (GLT) allows us to understand the spectral distribution of the eigenvalues through the study of what is known as GLT symbol.

The following plot show the eigenvalues of the stiffness matrix for different B-Splines degrees:

We notice the 0 at the 0-mode. But, we also notice a numerical 0, when increasing the Spline degree. This is in fact a specific pathology when using B-Splines as Finite Elements.

The 0 at 0-mode can be solved by a Multigrid method. For the other numerical zeros, one need an ad-hoc treatment. We are currently developing a parallel Multigrid that fixes this pathology in our framework.

If we take a more complicated, but still simple, problem such as a regularized Maxwell, for which a weak formulation is given by:

Find \mathbf{u} \in \mathbf{H}^1([0,1]^2) such that

(1)   \begin{equation*}\alpha(\nabla \times \mathbf{u},\nabla \times \mathbf{v}) +\beta(\nabla \cdot \mathbf{u},\nabla \cdot \mathbf{v}), \quad \forall \mathbf{v}\in  \mathbf{H}^1([0,1]^2),  \end{equation*}

with 0<\alpha,\beta\le1.

where we recall that \mathbf{H}^1([0,1]^2) = H({\rm curl},[0,1]^2)\cap H({\rm div},[0,1]^2), when the geometry has no re-entrant corners.

In the following tables, I give the number of iterations needed to reach convergence, depending on the values of \alpha and \beta.

In order to construct such solver, we used the GLT symbol related to the weak formulation \label{curl-div}, and performed its study. In general, we can compute the symbol of any bilinear form, and now it is even possible to do it automatically, thanks to the Python library GeLaTo. More about it will be covered in the section “Domain Specific Languages”.

Spectral Analysis and GLT theory

Let us first recall the following definition of a spectral distribution:

Let f:D\subset\mathbb R^d\to\mathbb C,\quad 0<\mbox{mesure}(D)<\infty.
f is a spectral of the sequence of matrices {A_n}_n
the eigenvalues of A_n are approximately a uniform sampling of f over D.
In this case, we write \{A_n\}_n \sim_\lambda f


From these sequences, we can construct an algebra called algebra of GLT sequences:

GLT sequences := *-algebra of matrix-sequences given as a closure under linear combination, product, inversion^{\ast} or conjugation of {T_n(f)}_n, {D_n(a)}_n, and {Z_n}_n.

Since GLT form an involutive algebra, it is very easy to implement a formal language that computes automatically the symbol associated to any weak formulation. Formally speaking, the inference rules for the GLT are

In addition to the previous inference rules, IGA provides as with other rules

The last rule gives an example of how to go from the 1d case to 2d in the case of a weighted mass operator.

Anisotropic diffusion operator

Let us take the example of the anisotropic diffusion operator, written in a weak form

(2)   \begin{align*}  a(v, u) = \int_{\Omega} \kappa_{\parallel}(\textbf{b} \cdot \nabla v) (\textbf{b} \cdot \nabla u) + \kappa_{I} \nabla v \cdot \nabla u d\Omega, \quad \forall u,v \in \mathcal{V}_h  \end{align*}

Using the GeLaTo library (or performing the computations by hand) lead to the following GLT symbol associated to a

    \begin{align*} \frac{\kappa_{\parallel} b_{x}^{2} n_{x} \mathfrak{m}_{p_y} \mathfrak{s}_{p_x}}{n_{y}} + 2 \kappa_{\parallel} b_{x} b_{y} \mathfrak{a}_{p_x} \mathfrak{a}_{p_y} + \frac{\kappa_{\parallel} b_{y}^{2} n_{y} \mathfrak{m}_{p_x} \mathfrak{s}_{p_y}}{n_{x}} \\+ \frac{\kappa_{\perp} n_{x} \mathfrak{m}_{p_y} \mathfrak{s}_{p_x}}{n_{y}} + \frac{\kappa_{\perp} n_{y} \mathfrak{m}_{p_x} \mathfrak{s}_{p_y}}{n_{x}} \label{eq:anidiff-symbol} \end{align*}

Comparing the numerical eigenvalues and the evaluation of the GLT symbol leads is shown in the following plot (cubic B-Splines on 32 \times 32 grid)