velocity field breaks hybrid DG

Apologies in advance if this question is off-topic, but it seemed like the kind of thing someone here might be able to answer.

I’m trying to solve a time-dependent problem that involves both convection and a non-linear source term. I had previously implemented a simpler version of the problem (omitting the nonlinear term) using a standard DG formulation. I’d like to add the non-linear term, but unfortunately, linearization for DG doesn’t seem to work. So, as suggested here, I’ve been working on switching over to HDG.

The problem I’m having is that my velocity field seems to make the matrix singular. I’m not sure if this is because of a bug in my code, or because the method itself isn’t suitable. I’m currently leaning towards the latter, because everything behaves as expected when I use a constant velocity (except, curiously, when that constant is zero. That also makes the matrix singular). For context: my velocity field is zero in some parts of the domain, it has constant positive divergence in others, and it is discontinuous on the boundary separating those two regions.

I’m not terribly familiar with hybrid methods, so I was hoping someone could give me some rough indication of how robust they’re expected to be. Should we expect HDG to fail with such a velocity field, or is it more likely that there’s a bug in my code?

Any guidance would be greatly appreciated.

Best regards,


Hi Julien,

how are you doing the linearisation? As far as I am aware BilinearForm.AssembleLinearization does not work with DG-like terms (jumps ect.). As a result, I always do the linearised form by hand in a separate BFI. Maybe thats worth a try for you?

Best wishes

Hi Henry,

Yes, that’s definitely worth trying. If I can’t get HDG to work I’ll try that next.



for pure convection problem, if your (normal) velocity is zero in certain part of facets, the standard uowinding hdg numerical flux lead to zero entries in the martix for the coeersponding facet dofs. as a workaround, you may use a modified flux with a small positive stab term, like
uhatup = IfPos(v.n, u, uhat)+stab*(u-uhat)
stab shall be activated only when v.n is nearly zero

Hi Guosheng,

Thank you for the suggestion. Unfortunately, it doesn’t seem to have resolved the issue. The matrix inversion still fails when I use a zero velocity field or the more complicated velocity field I described earlier.

I’ve attached a simplified version of my code (minus the calculation of the velocity field) in case anyone wants to take a look and tell me if there’s something obviously wrong with it.

best regards,