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?

  • $\begingroup$ You mention matrix is singular. It hints you may need a regularizer. $\endgroup$ May 21, 2018 at 19:06
  • $\begingroup$ 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. $\endgroup$
    – likask
    May 24, 2018 at 17:08


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.