# Getting started with FEM: Ill-conditioned matrix when evaluating flux terms in conservation law?

I have a system of conservation laws of the form

$$\frac{\partial \mathbf{q}}{\partial t} + \nabla \cdot \mathbf{F}\!\left(\mathbf{q}\right) = 0$$

I want to use finite elements to solve this system. As an initial choice, I am using Legendre polynomials $\phi$ for my function basis. Let $j$ be a polynomial index. If I substitute the decomposition of my unknown function $\mathbf{q}$ and the flux function $\mathbf{F}$,

$$\mathbf{q}\!\left(t, x\right) = \sum\limits_j \hat{\mathbf{q}}_j\!\left(t\right) \phi_j\!\left(x\right)$$

$$\mathbf{F}\!\left(\mathbf{q}\right) = \sum\limits_j \hat{\mathbf{F}}\!\left(t\right) \phi_j\!\left(x\right)$$

and finally write out the weak form (given some test volume $\Omega$ with boundary $\partial \Omega$)

$$\sum\limits_j \left(\frac{\partial \hat{\mathbf{q}}_j}{\partial t} \int\limits_\Omega \phi_i \phi_j + \hat{\mathbf{h}}_j \oint\limits_{\partial \Omega} \phi_i \phi_j - \hat{\mathbf{F}}_j \cdot \int\limits_\Omega \phi_j \nabla \phi_i\right) = 0$$

Here, I replaced $\mathbf{F} \cdot \hat{\mathbf{n}} = \mathbf{h}$.

My problem occurs when I actually try to evaluate the surface second and third terms inside the parentheses. Because my problem is 1-dimensional, the second term amounts to just evaluating $\phi_i \phi_j$ at 1 and -1 and taking the difference ($\hat{\mathbf{n}} = \pm \hat{\mathbf{x}}$ is the outward-facing normal). In analogy with the stiffness matrix, we could write this as a matrix $\mathbb{B} = \{b_{ij}\}$. Any given element $b_{ij}$ will be either 0 or 2 because the Legendre polynomials take the value of 1 or -1 at the edges of the domain. This matrix is sparse but singular. The third term exhibits a similar problem if we try to write it as a matrix $\mathbb{C} = \{c_{ij}\}$.

The only way I know to solve for the $\hat{\mathbf{F}}_j$ and $\hat{\mathbf{h}}_j$ coefficients is by solving a matrix equation, but all my attempts to solve them (using PETSc) converge slowly and may not even be correct. I would like to use a preconditioner, but I get errors any time I try (PCSETUP_FAILED due to SUBPC_ERROR), even though I have specified the null spaces of these matrices.

Have I approached this problem incorrectly altogether? Or is there a common way of overcoming this challenge, in PETSc or otherwise?

• You mention matrix is singular. It hints you may need a regularizer. – mathreadler May 21 '18 at 19:06
• You have mix formulation and have two problems here. 1) missing physical equation 2) q has to be in L2 space and F in h-div space to have stability of discrete system. – likask May 24 '18 at 17:08