Step-51 shows instability when implemented with IPHDG and DGQ spaces

[Greg Wang]

Hi there,

I have recently implemented a hybridizable discontinuous Galerkin method for the convection-diffusion equation adapted from step-51. Instead of the LHDG method used therein, I used IPHDG.

Some major differences between the two include that IPHDG doesn't have a vector-valued auxiliary variable q as in step-51 and that an interior penalty term <alpha * kappa / hK (u_h - uhat_h),(v_h - vhat_h)> is present in IPHDG where kappa is the diffusion parameter, hK is the mesh size and alpha is the interior penalty parameter, usually set to be 4.0*p^2 (p = fe.degree).

When solving an example 2D problem on a unit square domain and a uniformly refined mesh, and when using FE_DGP space paired with FE_FaceP space, things worked out fine across the board with kappa = 1, 1e-2, 1e-6.

When FE_DGQ and FE_FaceQ are used, however, things went off the rails where I had to tune up the IP parameter alpha to at least 4.0*p^2/h_K to retain the correct rates of convergence. In fact, the bigger the diffusion parameter kappa, the earlier the instability shows which seems to confirm that the issue lies in the interior penalty term.

All of this happened when I didn't change the code at all except for switching *P spaces to *Q spaces.

I'm wondering whether there is anyone here who has experience with IPHDG and/or using FE_DGQ vs FE_DGP in general. Any suggestions and insights would be highly appreciated!

Best regards,

Greg

[Greg Wang]

Hi again,

Just an update that the issue has been located at this chunk of code in step-51 (issue for step-51 when adapted to IPHDG, not step-51 itself):

for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face) for (unsigned int i=0; i<fe_local.dofs_per_cell; ++i) { if (fe_local.has_support_on_face(i,face)) fe_local_support_on_face[face].push_back(i); } for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face) for (unsigned int i=0; i<fe.dofs_per_cell; ++i) { if (fe.has_support_on_face(i,face)) fe_support_on_face[face].push_back(i); }

This bit serves to save the program from assembling a number of zeros terms on each cell particularly when we use LHDG with Q-polynomials. However, for IPHDG, this will assemble the consistency term and the symmetry term < u_h - uhat_h, nabla*n v_h > + < v_h - vhat_h, nabla*n u_h > incorrectly since nabla*n u_h does not necessarily vanish on a face when u_h vanishes on that face. Using P-polynomials dodges this bullet since they have support on every face.

Best,

Greg