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). 281
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? 935